text:
export  excelPopReport,
        excelPopEvolutionReport

function excelPopReport( mpSim::MPsim, timeGrid::Vector{T},
    nodes::AbstractString...; filename::AbstractString = "popReport",
    overwrite::Bool = true, timeFactor::Real = 12.0 ) where T <: Real

    filename, overwrite = setupFile( filename, overwrite )
    tStart = now()
    popReport = nodePopReport( mpSim, timeGrid, nodes... )
    reportGenerationTime = ( now() - tStart ).value / 1000.0

    if isempty( popReport )
        return
    end  # if isempty( popReport )

    writeNodePopData( popReport, mpSim, timeFactor, reportGenerationTime,
        filename, overwrite )

end  # excelPopReport( mpSim, timeGrid, nodes, filename, overwrite, timeFactor )

excelPopReport( mpSim::MPsim, timeRes::Real, nodes::AbstractString...;
    filename::AbstractString = "popReport", overwrite::Bool = true,
    timeFactor::Real = 12.0 ) =
    excelPopReport( mpSim, MP.generateTimeGrid( mpSim, timeRes ), nodes...;
    filename = filename, overwrite = overwrite, timeFactor = timeFactor )

function excelPopEvolutionReport( mpSim::MPsim, timeGrid::Vector{T},
    nodes::AbstractString...; filename::AbstractString = "popReport",
    overwrite::Bool = true, timeFactor::Real = 12.0 ) where T <: Real

    filename, overwrite = setupFile( filename, overwrite )
    tStart = now()
    popReport = nodeEvolutionReport( mpSim, timeGrid, nodes... )
    reportGenerationTime = ( now() - tStart ).value / 1000.0

    if isempty( popReport[1] )
        return
    end  # if isempty( popReport[1] )

    writeNodeEvolutionData( popReport, mpSim, timeFactor, reportGenerationTime,
        filename, overwrite )

end  # excelPopEvolutionReport( mpSim, timeGrid, nodes, filename, overwrite,
     #   timeFactor )

excelPopEvolutionReport( mpSim::MPsim, timeRes::Real, nodes::AbstractString...;
    filename::AbstractString = "popReport", overwrite::Bool = true,
    timeFactor::Real = 12.0 ) =
    excelPopEvolutionReport( mpSim, MP.generateTimeGrid( mpSim, timeRes ),
    nodes...; filename = filename, overwrite = overwrite,
    timeFactor = timeFactor )

include( "private/nodepop.jl" )


text:
import DataFrames

"""
"""
function fit!(
        transformer::ImmutablePredictionsSingleLabelInt2StringTransformer,
        varargs...;
        kwargs...
        )
    if length(varargs) == 1
        return varargs[1]
    else
        return varargs
    end
end

"""
"""
function predict(
        transformer::ImmutablePredictionsSingleLabelInt2StringTransformer,
        single_labelpredictions::AbstractVector,
        varargs...;
        kwargs...
        )
    single_labelpredictions = parse.(Int, single_labelpredictions)
    labelint2stringmap = getlabelint2stringmap(
        transformer.levels,
        transformer.index,
        )
    result = Vector{String}(
        undef,
        length(single_labelpredictions),
        )
    for i = 1:length(result)
        result[i] = labelint2stringmap[single_labelpredictions[i]]
    end
    return result
end

"""
"""
function predict(
        transformer::ImmutablePredictionsSingleLabelInt2StringTransformer,
        single_labelpredictions::DataFrames.AbstractDataFrame,
        varargs...;
        kwargs...
        )
    label_names = DataFrames.names(single_labelpredictions)
    result = DataFrames.DataFrame()
    for i = 1:length(label_names)
        result[label_names[i]] = predict(
            transformer,
            single_labelpredictions[label_names[i]];
            kwargs...
            )
    end
    return result
end

"""
"""
function predict_proba(
        transformer::ImmutablePredictionsSingleLabelInt2StringTransformer,
        varargs...;
        kwargs...
        )
    if length(varargs) == 1
        return varargs[1]
    else
        return varargs
    end
end

ext: jl


text:
# This example demonstrated how to calculate one-loop diagram of free electrons using the Monte Carlo module
# Observable is normalized: Γ₄*N_F where N_F is the free electron density of states

using Distributed
using QuantumStatistics, LinearAlgebra, Random, Printf,  BenchmarkTools, InteractiveUtils, Parameters

const Ncpu = 1 # number of workers (CPU)
const totalStep = 1e7 # MC steps of each worker

include("RPA.jl") # dW0 will be only calculated once in the master, then distributed to other workers. Therefore, there is no need to import RPA.jl for all workers.

addprocs(Ncpu)

@everywhere using QuantumStatistics, Parameters, Random, LinearAlgebra
@everywhere include("parameter.jl")
@everywhere include("interaction.jl")

# claim all globals to be constant, otherwise, global variables could impact the efficiency
########################### parameters ##################################
@everywhere const IsF = false # calculate quasiparticle interaction F or not
@everywhere const AngSize = 16

########################## variables for MC integration ##################
@everywhere const KInL = [kF, 0.0, 0.0] # incoming momentum of the left particle
@everywhere const Qd = [0.0, 0.0, 0.0] # transfer momentum is zero in the forward scattering channel

@everywhere struct Para{Q,T}
    extAngle::Vector{Float64}
    dW0::Matrix{Float64}
    qgrid::Q
    τgrid::T
    function Para(AngSize)
        extAngle = collect(LinRange(0.0, π, AngSize)) # external angle grid
        qgrid = Grid.boseK(kF, 6kF, 0.2kF, 256)
        τgrid = Grid.tau(β, EF / 20, 128)

        vqinv = [(q^2 + mass2) / (4π * e0^2) for q in qgrid.grid]
        dW0 = dWRPA(vqinv, qgrid.grid, τgrid.grid, kF, β, spin, me) # dynamic part of the effective interaction
        return new{typeof(qgrid),typeof(τgrid)}(extAngle, dW0, qgrid, τgrid)
    end
end

@everywhere function phase(tInL, tOutL, tInR, tOutR)
    if (IsF)
        return cos(π * ((tInL + tOutL) - (tInR + tOutR)));
    else
        return cos(π * ((tInL - tOutL) + (tInR - tOutR)))
    end
end

@everywhere function integrand(config)
    if config.curr == 1
        return Weight(1.0, 0.0) # return a weight!
    elseif config.curr == 2
        return eval2(config)
    else
        error("Not implemented!")
    end
end

@everywhere function eval2(config)
    T, K, Ang = config.var[1], config.var[2], config.var[3]
    k1, k2 = K[1], K[1] - Qd
    t1, t2 = T[1], T[2] # t1, t2 both have two tau variables
    θ = config.para.extAngle[Ang[1]] # angle of the external momentum on the right
    KInR = [kF * cos(θ), kF * sin(θ), 0.0]

    vld, wld, vle, wle = vertexDynamic(config, Qd, KInL - k1, t1[1], t1[2])
    vrd, wrd, vre, wre = vertexDynamic(config, Qd, KInR - k2, t2[1], t2[2])

    ϵ1, ϵ2 = (dot(k1, k1) - kF^2) / (2me), (dot(k2, k2) - kF^2) / (2me)
    wd, we = 0.0, 0.0
    # possible green's functions on the top
    gt1 = Spectral.kernelFermiT(t2[1] - t1[1], ϵ1, β)

    # gt2 = Spectral.kernelFermiT(t1[1] - t2[1], ϵ2, β)
    # wd += 1.0 / β * 1.0 / β * gt1 * gt2 / (2π)^3 * phase(t1[1], t1[1], t2[1], t2[1])

    # gt3 = Spectral.kernelFermiT(t1[1] - t2[2], ϵ2, β)
    # G = gt1 * gt3 / (2π)^3 * phase(t1[1], t1[1], t2[2], t2[1])
    # wd += G * (vld * wre)

    # wd += spin * (vld + wld) * (vrd + wrd) * gt1 * gt2 / (2π)^3 * phase(t1[1], t1[1], t2[1], t2[1])
    # println(vld, ", ", wld, "; ", vrd, ", ", wld)

    ############## Diagram v x v ######################
    """
      KInL                      KInR
       |                         |
  t1.L ↑     t1.L       t2.L     ↑ t2.L
       |-------------->----------|
       |       |    k1    |      |
       |   ve  |          |  ve  |
       |       |    k2    |      |
       |--------------<----------|
  t1.L ↑    t1.L        t2.L     ↑ t2.L
       |                         |
      KInL                      KInR
"""
    gd1 = Spectral.kernelFermiT(t1[1] - t2[1], ϵ2, β)
    G = gt1 * gd1 / (2π)^3 * phase(t1[1], t1[1], t2[1], t2[1])
    we += G * (vle * vre)
    ##################################################

    ############## Diagram w x v ######################
    """
      KInL                      KInR
       |                         |
  t1.R ↑     t1.L       t2.L     ↑ t2.L
       |-------------->----------|
       |       |    k1    |      |
       |   we  |          |  ve  |
       |       |    k2    |      |
       |--------------<----------|
  t1.L ↑    t1.R        t2.L     ↑ t2.L
       |                         |
      KInL                      KInR
    """
    gd2 = Spectral.kernelFermiT(t1[2] - t2[1], ϵ2, β)
    G = gt1 * gd2 / (2π)^3 * phase(t1[1], t1[2], t2[1], t2[1])
    we += G * (wle * vre)
    ##################################################

    ############## Diagram v x w ######################
    """
      KInL                      KInR
       |                         |
  t1.L ↑     t1.L       t2.L     ↑ t2.L
       |-------------->----------|
       |       |    k1    |      |
       |   ve  |          |  we  |
       |       |    k2    |      |
       |--------------<----------|
  t1.L ↑    t1.L        t2.R     ↑ t2.R
       |                         |
      KInL                      KInR
    """
    gd3 = Spectral.kernelFermiT(t1[1] - t2[2], ϵ2, β)
    G = gt1 * gd3 / (2π)^3 * phase(t1[1], t1[1], t2[2], t2[1])
    we += G * (vle * wre)
    ##################################################

    ############## Diagram w x w ######################
    """
      KInL                      KInR
       |                         |
  t1.R ↑     t1.L       t2.L     ↑ t2.L
       |-------------->----------|
       |       |    k1    |      |
       |   we  |          |  we  |
       |       |    k2    |      |
       |--------------<----------|
  t1.L ↑    t1.R        t2.R     ↑ t2.R
       |                         |
      KInL                      KInR
"""
    gd4 = Spectral.kernelFermiT(t1[2] - t2[2], ϵ2, β)
    G = gt1 * gd4 / (2π)^3 * phase(t1[1], t1[2], t2[2], t2[1])
    we += G * (wle * wre)
    ##################################################

    return Weight(wd, we)
end

@everywhere function measure(config)
    angidx = config.var[3][1]
    factor = 1.0 / config.reweight[config.curr]
    if config.curr == 1
        config.observable[1][:, angidx] .+= factor
    elseif config.curr == 2
        weight = integrand(config)
        config.observable[2][:, angidx] .+= weight / abs(weight) * factor
    else
        error("Not implemented!")
    end
end

@everywhere normalize(config) = config.observable[2] / sum(config.observable[1]) * AngSize * β

function run(totalStep)
    T = MonteCarlo.TauPair(β, β / 2.0)
    K = MonteCarlo.FermiK(3, kF, 0.2 * kF, 10.0 * kF)
    Ext = MonteCarlo.Discrete(1, AngSize) # external variable is specified

    dof = ([1, 0, 1], [2, 1, 1])
    obs = (zeros(Float64, (2, AngSize)), zeros(Float64, (2, AngSize)))

    para = Para(AngSize)

    avg, std = MonteCarlo.sample(totalStep, (T, K, Ext), dof, obs, integrand, measure, normalize; para=para, print=10)

    NF = TwoPoint.LindhardΩnFiniteTemperature(dim, 0.0, 0, kF, β, me, spin)[1]
    println("NF = $NF")

    avg = avg * NF
    std = std * NF

    println(size(avg))
    for (idx, angle) in enumerate(para.extAngle)
        @printf("%10.6f   %10.6f ± %10.6f  %10.6f ± %10.6f\n", angle, avg[1, idx], std[1,idx], avg[2,idx], std[2,idx])
    end
end

# @btime run(1, 10)
run(totalStep)
# @time run(Repeat, totalStep)


text:
using MCTS
using POMDPToolbox
using POMDPs
using Multilane
using POMCPOW
using D3Trees
using ParallelPOMCPOW

@everywhere begin
    using POMDPs
    using ParallelPOMCPOW
    using Multilane
    POMDPs.actions(p::MLPOMDP, fsrb::ParallelPOMCPOW.FilteredSRBelief) = actions(p, fsrb.srb)
end

tii = true
wup = WeightUpdateParams(smoothing=0.0, wrong_lane_factor=0.5)

solver = POMCPOWSolver(tree_queries=10_000,
                               criterion=MaxUCB(10.0),
                               max_depth=20,
                               max_time=Inf,
                               enable_action_pw=false,
                               k_observation=4.0,
                               alpha_observation=1/8,
                               estimate_value=FORollout(SimpleSolver()),
                               check_repeat_obs=false,
                               node_sr_belief_updater=AggressivenessPOWFilter(wup),
                               tree_in_info=tii
                              )

psolver = ParallelPOMCPOWSolver(solver, 1000)

behaviors = standard_uniform(correlation=0.75)
pp = PhysicalParam(4, lane_length=100.0)
dmodel = NoCrashIDMMOBILModel(10, pp,
                              behaviors=behaviors,
                              p_appear=1.0,
                              lane_terminate=true,
                              max_dist=500.0
                             )
rmodel = SuccessReward(lambda=2,
                       target_lane=4,
                       brake_penalty_thresh=4.0
                      )
pomdp = NoCrashPOMDP{typeof(rmodel)}(dmodel, rmodel, 0.95, false)
mdp = NoCrashMDP{typeof(rmodel)}(dmodel, rmodel, 0.95, false)
is = initial_state(pomdp, Base.GLOBAL_RNG)
ips = MLPhysicalState(is)

agg_up = AggressivenessUpdater(pomdp, 500, 0.1, 0.1, wup, MersenneTwister(50000))
ib = initialize_belief(agg_up, ips)

planner = solve(solver, pomdp)
pplanner = solve(psolver, pomdp)

println("Single Thread")
action(planner, ib)
@time action(planner, ib)

println("Parallel")
action(pplanner, ib)
@time action(pplanner, ib)


text:
# Prints all paths from 's' to 't'
function get_all_paths(data, opt, sp)

  function getAllPathsUtil(arc::Int, u::Int, t::Int, visited::BitArray{1}, path::BitArray{1}, cost::Float64)
    # Mark the current node as visited and store in path
    if cost + sp[u] < opt
      visited[u] = true
      if arc > 0
        path[arc] = true
      end
      # If current vertex is same as destination, then store current path
      if u == t
        push!(X,copy(path))
      else
        for new_arc in data.δ⁺[u]
          if !visited[data.to[new_arc]]
            cost = cost + data.hc[new_arc]
            getAllPathsUtil(new_arc, data.to[new_arc], t, visited, path, cost)
            cost = cost - data.hc[new_arc]
          end
        end
      end
      # Remove current vertex from path[] and mark it as unvisited
      if arc > 0
        path[arc] = false
      end
      visited[u] = false
      return X
    end
  end

  X = Vector{BitArray{1}}()
  path = falses(data.n)
  cost = 0.0

  # Mark all the vertices as not visited
  visited = falses(data.nV)

  # Call the recursive helper function to print all paths
  getAllPathsUtil(0, data.s, data.t, visited, path, cost)
  return X
end

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
@testset "FixedTimeZone" begin
    @test isbitstype(FixedTimeZone)

    @test FixedTimeZone("0123") == FixedTimeZone("UTC+01:23", 4980)
    @test FixedTimeZone("+0123") == FixedTimeZone("UTC+01:23", 4980)
    @test FixedTimeZone("-0123") == FixedTimeZone("UTC-01:23", -4980)
    @test FixedTimeZone("01:23") == FixedTimeZone("UTC+01:23", 4980)
    @test FixedTimeZone("+01:23") == FixedTimeZone("UTC+01:23", 4980)
    @test FixedTimeZone("-01:23") == FixedTimeZone("UTC-01:23", -4980)
    @test FixedTimeZone("01:23:45") == FixedTimeZone("UTC+01:23:45", 5025)
    @test FixedTimeZone("+01:23:45") == FixedTimeZone("UTC+01:23:45", 5025)
    @test FixedTimeZone("-01:23:45") == FixedTimeZone("UTC-01:23:45", -5025)
    @test FixedTimeZone("99:99:99") == FixedTimeZone("UTC+99:99:99", 362439)
    @test FixedTimeZone("UTC") == FixedTimeZone("UTC", 0)
    @test FixedTimeZone("UTC+00") == FixedTimeZone("UTC", 0)
    @test FixedTimeZone("UTC+1") == FixedTimeZone("UTC+01:00", 3600)
    @test FixedTimeZone("UTC-1") == FixedTimeZone("UTC-01:00", -3600)
    @test FixedTimeZone("UTC+01") == FixedTimeZone("UTC+01:00", 3600)
    @test FixedTimeZone("UTC-01") == FixedTimeZone("UTC-01:00", -3600)
    @test FixedTimeZone("UTC+0123") == FixedTimeZone("UTC+01:23", 4980)
    @test FixedTimeZone("UTC-0123") == FixedTimeZone("UTC-01:23", -4980)

    @test FixedTimeZone("+01") == FixedTimeZone("UTC+01:00", 3600)
    @test FixedTimeZone("-02") == FixedTimeZone("UTC-02:00", -7200)
    @test FixedTimeZone("+00:30") == FixedTimeZone("UTC+00:30", 1800)
    @test FixedTimeZone("-00:30") == FixedTimeZone("UTC-00:30", -1800)

    @test_throws ArgumentError FixedTimeZone("1")
    @test_throws ArgumentError FixedTimeZone("01")
    @test_throws ArgumentError FixedTimeZone("123")
    @test_throws ArgumentError FixedTimeZone("012345")
    @test_throws ArgumentError FixedTimeZone("0123:45")
    @test_throws ArgumentError FixedTimeZone("01:2345")
    @test_throws ArgumentError FixedTimeZone("01:-23:45")
    @test_throws ArgumentError FixedTimeZone("01:23:-45")
    @test_throws ArgumentError FixedTimeZone("01:23:45:67")
    @test_throws ArgumentError FixedTimeZone("UTC1")
    @test_throws ArgumentError FixedTimeZone("+1")
    @test_throws ArgumentError FixedTimeZone("-2")

    @testset "broadcastable" begin
        # Validate that FixedTimeZone is treated as a scalar during broadcasting
        fixed_tz = FixedTimeZone("UTC")
        @test size(fixed_tz .== fixed_tz) == ()
    end
end


text:
module PiecewiseOrthogonalPolynomials
using ClassicalOrthogonalPolynomials, LinearAlgebra, BlockArrays, BlockBandedMatrices, ContinuumArrays, QuasiArrays, LazyArrays, LazyBandedMatrices, FillArrays

import BlockArrays: BlockSlice, block, blockindex, blockvec
import BlockBandedMatrices: _BandedBlockBandedMatrix
import ClassicalOrthogonalPolynomials: grid, massmatrix
import ContinuumArrays: @simplify, factorize, TransformFactorization, AbstractBasisLayout, MemoryLayout, layout_broadcasted, ExpansionLayout, basis
import LazyArrays: paddeddata
import LazyBandedMatrices: BlockBroadcastMatrix
import Base: axes, getindex, ==, \, OneTo

export PiecewisePolynomial, ContinuousPolynomial, Derivative, Block

abstract type AbstractPiecewisePolynomial{order,T,P<:AbstractVector} <: Basis{T} end

struct PiecewisePolynomialLayout{order} <: AbstractBasisLayout end
MemoryLayout(::Type{<:AbstractPiecewisePolynomial{order}}) where order = PiecewisePolynomialLayout{order}()

struct PiecewisePolynomial{T,Bas,P<:AbstractVector} <: AbstractPiecewisePolynomial{0,T,P}
    basis::Bas
    points::P
end

PiecewisePolynomial{T}(basis::AbstractQuasiMatrix, points::AbstractVector) where {T} =
    PiecewisePolynomial{T,typeof(basis),typeof(points)}(basis, points)
PiecewisePolynomial(basis::AbstractQuasiMatrix{T}, points::AbstractVector) where {T} =
    PiecewisePolynomial{T}(basis, points)

axes(B::PiecewisePolynomial) = (Inclusion(first(B.points) .. last(B.points)), blockedrange(Fill(length(B.points) - 1, ∞)))

function repeatgrid(ax, g, pts)
    ret = Matrix{eltype(g)}(undef, length(g), length(pts) - 1)
    for j in axes(ret, 2)
        ret[:, j] = affine(ax, pts[j] .. pts[j+1])[g]
    end
    ret
end

function grid(V::SubQuasiArray{T,2,<:PiecewisePolynomial,<:Tuple{Inclusion,BlockSlice}}) where {T}
    P = parent(V)
    _, JR = parentindices(V)
    N = Int(last(JR))
    g = grid(P.basis[:, OneTo(N)])
    repeatgrid(axes(P.basis, 1), g, P.points)
end

struct ApplyFactorization{T, FF, FAC<:Factorization{T}} <: Factorization{T}
    f::FF
    F::FAC
end

\(P::ApplyFactorization, f) = P.f(P.F \ f)

function factorize(V::SubQuasiArray{<:Any,2,<:PiecewisePolynomial,<:Tuple{Inclusion,BlockSlice}})
    P = parent(V)
    _,JR = parentindices(V)
    N = Int(last(JR.block))
    F = factorize(view(P.basis,:,OneTo(N)), length(P.points)-1)
    ApplyFactorization(v -> blockvec(permutedims(v)), TransformFactorization(repeatgrid(axes(P.basis, 1), F.grid, P.points), F.plan))
end

struct ContinuousPolynomial{order,T,P<:AbstractVector} <: AbstractPiecewisePolynomial{order,T,P}
    points::P
end

ContinuousPolynomial{o,T}(pts::P) where {o,T,P} = ContinuousPolynomial{o,T,P}(pts)
ContinuousPolynomial{o}(pts) where {o} = ContinuousPolynomial{o,Float64}(pts)
ContinuousPolynomial{o}(P::ContinuousPolynomial) where {o} = ContinuousPolynomial{o,eltype(P)}(P.points)

PiecewisePolynomial(P::ContinuousPolynomial{0,T}) where {T} = PiecewisePolynomial(Legendre{T}(), P.points)

axes(B::ContinuousPolynomial{0}) where {o} = axes(PiecewisePolynomial(B))
axes(B::ContinuousPolynomial{1}) where {o} =
    (Inclusion(first(B.points) .. last(B.points)), blockedrange(Vcat(length(B.points), Fill(length(B.points) - 1, ∞))))

==(::PiecewisePolynomial, ::ContinuousPolynomial{1}) = false
==(::ContinuousPolynomial{1}, ::PiecewisePolynomial) = false
==(A::ContinuousPolynomial{o}, B::ContinuousPolynomial{o}) where o = A.points == B.points

function getindex(P::PiecewisePolynomial{T}, x::Number, Kk::BlockIndex{1}) where {T}
    K, k = block(Kk), blockindex(Kk)
    b = searchsortedlast(P.points, x)
    if b == length(P.points) == k + 1 # last point
        P.basis[affine(P.points[end-1] .. P.points[end], axes(P.basis, 1))[x], Int(K)]
    elseif b == k
        P.basis[affine(P.points[b] .. P.points[b+1], axes(P.basis, 1))[x], Int(K)]
    else
        zero(T)
    end
end
getindex(P::ContinuousPolynomial{0,T}, x::Number, Kk::BlockIndex{1}) where {T} = PiecewisePolynomial(P)[x, Kk]

function getindex(P::ContinuousPolynomial{1,T}, x::Number, Kk::BlockIndex{1}) where {T}
    K, k = block(Kk), blockindex(Kk)
    if K == Block(1)
        LinearSpline(P.points)[x, k]
    else
        b = searchsortedlast(P.points, x)
        if b == k
            Weighted(Jacobi{T}(1, 1))[affine(P.points[b] .. P.points[b+1], ChebyshevInterval{real(T)}())[x], Int(K)-1]
        else
            zero(T)
        end
    end
end

getindex(P::AbstractPiecewisePolynomial, x::Number, k::Int) = P[x, findblockindex(axes(P, 2), k)]

factorize(V::SubQuasiArray{T,2,<:ContinuousPolynomial{0},<:Tuple{Inclusion,BlockSlice}}) where {T} =
    factorize(view(PiecewisePolynomial(parent(V)), parentindices(V)...))
grid(V::SubQuasiArray{T,2,<:ContinuousPolynomial{0},<:Tuple{Inclusion,BlockSlice}}) where {T} =
    grid(view(PiecewisePolynomial(parent(V)), parentindices(V)...))

function grid(V::SubQuasiArray{T,2,<:ContinuousPolynomial{1},<:Tuple{Inclusion,BlockSlice}}) where {T}
    P = parent(V)
    _, JR = parentindices(Q)
    pts = P.points
    grid(view(PiecewisePolynomial(Weighted(Jacobi{T}(1, 1)), pts), :, JR))
end

#######
# Conversion
#######

function \(P::ContinuousPolynomial{0}, C::ContinuousPolynomial{1})
    T = promote_type(eltype(P), eltype(C))
    # diag blocks based on
    # L = Legendre{T}() \ Weighted(Jacobi{T}(1,1))
    @assert P.points == C.points
    N = length(P.points)
    v = mortar(Fill.((convert(T, 2):2:∞) ./ (3:2:∞), N - 1))
    z = Zeros{T}(axes(v))
    H1 = BlockBroadcastArray(hcat, z, v)
    M1 = BlockVcat(Zeros{T}(N, 2), H1)
    M2 = BlockVcat(Ones{T}(N, 2) / 2, Zeros{T}((axes(v, 1), Base.OneTo(2))))
    H3 = BlockBroadcastArray(hcat, z, -v)
    M3 = BlockVcat(Hcat(Ones{T}(N) / 2, -Ones{T}(N) / 2), H3)
    dat = BlockHcat(M1, M2, M3)'
    _BandedBlockBandedMatrix(dat, axes(P, 2), (1, 1), (0, 1))
end

######
# Gram matrix
######

@simplify function *(Ac::QuasiAdjoint{<:Any,<:ContinuousPolynomial{0}}, B::ContinuousPolynomial{0})
    A = Ac'
    T = promote_type(eltype(A), eltype(B))
    r = A.points
    @assert r == B.points
    N = length(r)
    M = massmatrix(Legendre{T}())
    Diagonal(mortar(Fill.((step(r) / 2) .* M.diag, N - 1)))
end

@simplify function *(Ac::QuasiAdjoint{<:Any,<:ContinuousPolynomial}, B::ContinuousPolynomial)
    A = Ac'
    P = ContinuousPolynomial{0}(A)
    Q = ContinuousPolynomial{0}(B)
    (P \ A)' * (P'Q) * (Q \ B)
end

#####
# Derivative
#####

@simplify function *(D::Derivative, C::ContinuousPolynomial{1})
    T = promote_type(eltype(D), eltype(C))

    # Legendre() \ (D*Weighted(Jacobi(1,1)))
    r = C.points
    N = length(r)
    v = mortar(Fill.((-convert(T, 2):-2:-∞), N - 1))
    z = Zeros{T}(axes(v))
    H = BlockBroadcastArray(hcat, z, v)
    M = BlockVcat(Hcat(Ones{T}(N) / 2, -Ones{T}(N) / 2), H)
    P = ContinuousPolynomial{0}(C)
    P * _BandedBlockBandedMatrix(M', (axes(P, 2), axes(C, 2)), (0, 0), (0, 1))
end

### multiplication

function layout_broadcasted(::Tuple{ExpansionLayout{PiecewisePolynomialLayout{0}},PiecewisePolynomialLayout{0}}, ::typeof(*), a, P)
    @assert basis(a) == P
    _,c = a.args
    T = eltype(c)
    m = length(P.points)
    B = PiecewisePolynomial(P).basis

    ops = [(cₖ = [paddeddata(c)[k:m-1:end]; Zeros{T}(∞)];
            aₖ = B  * cₖ;
            unitblocks(B \ (aₖ .* B))) for k = 1:m-1]

    P * BlockBroadcastMatrix{T}(Diagonal, ops...)
end

function layout_broadcasted(::Tuple{ExpansionLayout{PiecewisePolynomialLayout{0}},PiecewisePolynomialLayout{1}}, ::typeof(*), a, C)
    P = ContinuousPolynomial{0}(C)
    (a .* P) * (P \ C)
end

end # module


text:
function mars_sky_opportunity()
   LSpec(MarsSkyOpp[:,1],MarsSkyOpp[:,2])
end

const MarsSkyOpp=
   [359.818401937046	   0.14974358974359;
    365.508474576271	   0.166324786324786;
    370.774818401937	   0.183931623931624;
    375.799031476998	   0.202051282051282;
    380.762711864407	   0.22034188034188;
    385.605326876513	   0.238974358974359;
    390.387409200969	   0.257777777777778;
    395.108958837772	   0.276410256410256;
    399.769975786925	   0.295384615384615;
    404.370460048426	   0.314358974358974;
    408.910411622276	   0.333504273504274;
    413.450363196126	   0.352820512820513;
    417.990314769976	   0.372136752136752;
    422.469733656174	   0.391452991452991;
    426.949152542373	   0.410598290598291;
    431.489104116223	   0.42957264957265;
    436.089588377724	   0.448547008547009;
    440.690072639225	   0.467521367521368;
    445.351089588378	   0.486495726495726;
    450.133171912833	   0.505299145299145;
    454.915254237288	   0.523931623931624;
    459.818401937046	   0.542393162393162;
    464.903147699758	   0.56034188034188;
    470.169491525424	   0.577948717948718;
    475.617433414044	   0.595213675213675;
    481.246973365617	   0.611965811965812;
    487.179176755448	   0.627692307692308;
    493.353510895884	   0.642735042735043;
    499.709443099274	   0.657264957264957;
    506.125907990315	   0.671623931623932;
    512.542372881356	   0.685982905982906;
    518.837772397094	   0.70051282051282;
    524.891041162228	   0.715897435897436;
    530.702179176755	   0.731965811965812;
    536.210653753027	   0.748888888888889;
    541.295399515738	   0.767008547008547;
    546.077481840194	   0.785811965811966;
    550.799031476998	   0.804786324786325;
    555.399515738499	   0.823760683760684;
    559.818401937046	   0.843076923076923;
    564.176755447942	   0.862564102564103;
    568.65617433414	   0.881880341880342;
    573.317191283293	   0.901025641025641;
    578.099273607748	   0.91965811965812;
    583.123486682809	   0.937777777777778;
    588.329297820823	   0.955555555555556;
    593.958837772397	   0.972649572649573;
    600.375302663438	   0.987179487179487;
    607.639225181598	   0.996752136752137;
    615.508474576271	   0.99965811965812;
    623.498789346247	   0.996410256410256;
    631.186440677966	   0.988547008547008;
    638.450363196126	   0.977606837606838;
    645.411622276029	   0.965128205128205;
    652.070217917676	   0.951794871794872;
    658.668280871671	   0.938290598290598;
    665.266343825666	   0.924615384615385;
    671.985472154964	   0.911452991452991;
    678.886198547215	   0.898974358974359;
    685.847457627119	   0.887008547008547;
    692.990314769976	   0.875726495726496;
    700.133171912833	   0.864615384615385;
    707.397094430993	   0.854188034188034;
    714.782082324455	   0.844102564102564;
    722.167070217918	   0.833846153846154;
    729.430992736078	   0.823418803418803;
    736.694915254237	   0.812820512820513;
    743.898305084746	   0.801880341880342;
    750.9200968523	   0.79008547008547;
    757.760290556901	   0.777606837606838;
    764.418886198547	   0.764273504273504;
    770.956416464891	   0.75042735042735;
    777.191283292978	   0.735555555555556;
    783.18401937046	   0.72;
    788.934624697337	   0.703418803418803;
    794.261501210654	   0.685982905982906;
    799.285714285714	   0.667863247863248;
    804.188861985472	   0.649401709401709;
    808.910411622276	   0.630598290598291;
    813.389830508475	   0.611282051282051;
    817.566585956416	   0.591452991452991;
    821.501210653753	   0.571111111111111;
    825.314769975787	   0.550598290598291;
    829.06779661017	   0.53008547008547]


text:
using LinearAlgebraicRepresentation
using Plasm,SparseArrays
Lar = LinearAlgebraicRepresentation

# input generation
V,EV = Lar.randomcuboids(30, .35)
V = Plasm.normalize(V,flag=true)
Plasm.view(Plasm.numbering(.05)((V,[[[k] for k=1:size(V,2)], EV])))

# subdivision of input edges
W = convert(Lar.Points, V')
cop_EV = Lar.coboundary_0(EV::Lar.Cells)
cop_EW = convert(Lar.ChainOp, cop_EV)
V, copEV = Lar.planar_arrangement_1(W::Lar.Points, cop_EW::Lar.ChainOp)

# compute containment graph of components
bicon_comps = Lar.Arrangement.biconnected_components(copEV)

# visualization of component graphs
EW = Lar.cop2lar(copEV)
W = convert(Lar.Points, V')
hpcs = [ Plasm.lar2hpc(W,[EW[e] for e in comp]) for comp in bicon_comps ]
Plasm.view([ Plasm.color(Plasm.colorkey[(k%12)==0 ? 12 : k%12])(hpcs[k]) for k=1:(length(hpcs)) ])

# computation of 2-cells and 2-boundary
V, copEV, copFE = Lar.planar_arrangement_2(V, copEV, bicon_comps)

# visualization of numbered arrangement
Plasm.view( Plasm.numbering1(0.05)((V, copEV, copFE)) )

# final solid visualization
triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
V = convert(Lar.Points, V')
FVs = convert(Array{Lar.Cells}, triangulated_faces)
Plasm.viewcolor(V::Lar.Points, FVs::Array{Lar.Cells})

# polygonal face boundaries
EVs = Lar.FV2EVs(copEV, copFE)
EVs = convert(Array{Array{Array{Int64,1},1},1}, EVs)
Plasm.viewcolor(V::Lar.Points, EVs::Array{Lar.Cells})

# exploded polygons
model = V,EVs
Plasm.view(Plasm.lar_exploded(model)(1.2,1.2,1.2))


text:
# Extract, preprocess and select attributes from a collection of tweet messages (to be filtered)
# run as:
# JULIA_PROJECT=. julia -p32 -L src/extract-messages-from-tweets.jl
# or
# JULIA_PROJECT=. srun --pty -N1 -xgeoint0 julia -p64 -L src/messages-from-raw-tweets.jl

using CSV, DataFrames, JLD2, StatsBase, TextSearch
include("io.jl")

function extract_messages(filename; lang="es")
    @info filename
    L = []

    valid_cc = Set(VALID_CC)
    eachtweet(filename) do tweet
        if language(tweet) == lang
            cc = countrycode(tweet)
            (isnothing(cc) || !(cc in valid_cc)) && return

            lat, long = try
                mean(tweet["place"]["bounding_box"]["coordinates"][1])
            catch e
                showerror(stderr, e)
                return
            end

            t = preprocess(text(tweet))
            occursin("_url", t) && return
            (startswith(t, "RT") || startswith(t, "I'm at") || count(c->c == ' ', t) <= 6) && return

            id = tweet["id"]
            user_id = tweet["user"]["id"]
            screen_name = tweet["user"]["screen_name"]
            timestamp = parse(Int, tweet["timestamp_ms"])
            push!(L, (id, user_id, screen_name, cc, long, lat, timestamp, t))
        end
    end

    L
end

function process_directory(files, outdir, outname)
    println(stderr, "processing $(length(files)) files => $(outdir)")

    L = DataFrame(
        id=Int64[],
        user_id=Int64[],
        screen_name=String[],
        country_code=String[],
        long=Float64[],
        lat=Float64[],
        timestamp_ms=Int64[],
        text=String[],
    )
    for filename in files
        for r in extract_messages(filename)
            push!(L, r)
        end
    end

    for g in groupby(L, :country_code)
        cc = first(g.country_code)
        output = joinpath(outdir, cc, outname)
        open(output, "w") do f
            gz = GzipCompressorStream(f)
            CSV.write(gz, g, delim='\t')
            close(gz)
        end
    end
end

function process_directory_list(dirlist, outdir, pat="GEO*.log.gz")
    println(stderr, "processing $(length(dirlist)) directories")

    for d in dirlist
        for cc in VALID_CC
            mkpath(joinpath(outdir, cc))
        end
    end

    @sync @distributed for d in dirlist
        outname = replace(d, "/" => "_") * ".tsv.gz"
        process_directory(glob(joinpath(d, pat)), outdir, outname)
    end
end

function main(lang)
    dirlist = glob("DatasetSpanish/16/*/*")
    append!(dirlist, glob("DatasetSpanish/17/*/*"))
    append!(dirlist, glob("DatasetSpanish/18/*/*"))
    append!(dirlist, glob("DatasetSpanish/19/*/*"))
    outdir = "data/$lang/messages"
    process_directory_list(dirlist, outdir)

    #=process_directory_list(glob("DatasetSpanish/17/*/*"), lang)
    process_directory_list(glob("DatasetSpanish/18/*/*"), lang)
    process_directory_list(glob("DatasetSpanish/19/*/*"), lang)=#
end

max_stars_repo_name: sadit/RegionalSpanish

max_issues_repo_name: sadit/RegionalSpanish

max_forks_repo_name: sadit/RegionalSpanish


text:
x = [
    0.0 0.3 0.6
    0.0 0.0 0.0
]
rs = 0.5:0.5:1.

dist = PeriodicEuclidean(ones(2))

@test K(rs, x) ≈ [2 / 3, 1.]
@test K(rs, x; dist) ≈ ones(2)

@test D(rs, x, 2) ≈ [1. 1/3; 1. 1.]
@test D(rs, x, 2; dist) ≈ ones(2, 2)

@test T(rs, x) ≈ [0., 1.]
@test T(rs, x; dist) ≈ [1., 1.]

@test H(rs, x, 0.5) ≈ [4 / 9, 7 / 9]
@test H(rs, x, 0.5; dist) ≈ [8 / 9, 1.]


text:
function PolynomialParametricAcceleration(arg0::Vector3D, arg1::AttitudeProvider, arg2::JString, arg3::AbsoluteDate, arg4::jint)
    return PolynomialParametricAcceleration((Vector3D, AttitudeProvider, JString, AbsoluteDate, jint), arg0, arg1, arg2, arg3, arg4)
end

function PolynomialParametricAcceleration(arg0::Vector3D, arg1::jboolean, arg2::JString, arg3::AbsoluteDate, arg4::jint)
    return PolynomialParametricAcceleration((Vector3D, jboolean, JString, AbsoluteDate, jint), arg0, arg1, arg2, arg3, arg4)
end

function depends_on_position_only(obj::PolynomialParametricAcceleration)
    return jcall(obj, "dependsOnPositionOnly", jboolean, ())
end

function get_parameters_drivers(obj::PolynomialParametricAcceleration)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function init(obj::PolynomialParametricAcceleration, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end


text:
#=
  Where is my bag problem.

  From
  "Probabilistic Reasoning Under Uncertainty with Bayesian Networks"
  https://www.bayesia.com/2018-03-02-probabilistic-reasoning-under-uncertainty

  (Also see The BaysiaLabBook "Bayesian Networks BayesiaLab" v1, page 51f.)

  Travel from Singapore -> Tokyo -> Los Angeles.
  1) The departure of Singapore -> Tokyo is delayed
  2) Just made the Tokyo -> LA trip. There's a 50/50 change the luggage is on the plane to LA.
  3) Arriving at LA, waiting for the luggage in 5 minutes, the luggage has not been seen.
  What is the probability that the bag was in the Tokyo -> LA plane?

  Answer: The probability that the bag was with the LA plane is 33% (as in the talk).

  From Judea Pearl: The Book of Why.

  For p(bag_on_plane) = 0.5
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
  bag_on_plane    0.3364    0.4725     0.0047    0.0053   8953.8759    1.0001     1581.1188

  For p(bag_on_plane) = 0.8
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
  bag_on_plane    0.6641    0.4723     0.0047    0.0051   8363.8461    0.9999     1469.6619

  Cf ~/blog/where_is_my_bag.blog
     ~/webppl/where_is_my_bag.wppl

=#
using Turing, StatsPlots, Distributions
include("jl_utils.jl")

@model function where_is_my_bag(maxTime=10,p=0.5,time=5.0,bag_on_carousel=false)

    # We know that the bag was on the plane with 50% probability.
    # Variant: https://www.youtube.com/watch?v=c71pnonOMkI&t=1074s it's 80/20.
    bag_on_plane ~ flip(p)

    # We assume uniform time step
    time ~ DiscreteUniform(1,maxTime)

    # Probability that the bag is on the carousel given a specific time
    # (if it was on the plane, that is).
    # The probability that the bag is on the carousel on a given time
    # is directly proportional on the time:
    #     time / maxTime
    bag_on_carousel ~ bag_on_plane == true ? flip(1.0*time/maxTime) : Dirac(0.0)

    # He waited 5 minutes without seeing the luggage.
    # Note: It's faster if we observe via the parameters
    # true ~ Dirac(time == 5.0)
    # true ~ Dirac(bag_on_carousel == false)

end

maxTime = 10
p = 0.5
# p = 0.8
time = 5.0
bag_on_carousel = false
model = where_is_my_bag(maxTime,p,time,bag_on_carousel)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
load("flatten.jl")
load("kernels.jl")

code = quote
    A[] = B.*C + D
end

symbols = Dict{Symbol,SymEntry}()
flat_code = {}
receive = ex->append!(flat_code, {ex})
context = Context(symbols, receive)

value = flatten(context, code)

print_code(flat_code)
println("value = $value")
print_context(context)

#indvars = gensym(32) # just some number of dims that should be enough
indvars = (:_i, :_j, :_k, :_l) # todo: use gensyms instead

nd = 2

staged = true
fdef = wrap_kernel(context, flat_code, indvars[1:nd], staged)
kernel = eval(fdef)

A = Array(Float, (2,3))
B = [1 2 3
     4 5 6]
C = [1 0 1
     0 1 0]
D = [ 0  0  0
     10 10 10]
kernel(A, B, C, D)

println("A =\n$A")

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
using MonadFunctions
using Test

@testset "MonadFunctions" begin

# common test functions that handle positive inputs only
incr(x)   = x > 0 ? x + 1 : NONE
decr(x)   = x > 0 ? x - 1 : NONE
double(x) = x > 0 ? 2x    : NONE

@testset "Maybe" begin

    @test 3       |> fmap(incr) == 4
    @test Just(3) |> fmap(incr) == Just(4)
    @test NONE    |> fmap(incr) == NONE
    @test nothing |> fmap(incr) == nothing

    @test 0 |> fmap(incr)               == NONE
    @test 1 |> fmap(decr) |> fmap(decr) == NONE

    @test 1       |> fmap(incr) |> fmap(double) == 4
    @test Just(1) |> fmap(incr) |> fmap(double) == Just(4)

    @test 1  |> fmap(incr) |> or_else(-1) == 2
    @test -1 |> fmap(incr) |> or_else(-1) == -1

    @test cata(()->0, x->x+1, 1)       == 2
    @test cata(()->0, x->x+1, Just(1)) == Just(2)
    @test cata(()->0, x->x+1, NONE)    == 0
    @test cata(()->0, x->x+1, nothing) == 0
end

@testset "Either" begin
    @test left(1)  |> is_left  == true
    @test left(1)  |> is_right == false

    @test right(1) |> is_left  == false
    @test right(1) |> is_right == true

    @test right(1) |> fmap(incr) == right(2)
    @test left(1)  |> fmap(incr) == left(2)

    @test left(1)  |> left_value == 1
    @test right(1) |> right_value == 1

    @test_throws MethodError right_value(left(1))
    @test_throws MethodError left_value(right(1))

    @test_throws ErrorException right(left(1))
    @test_throws ErrorException left(right(1))
end

@testset "Examples" begin

    err_connection = ErrorException("unable to connect to database")
    err_cursor     = ErrorException("unable to get cursor")
    err_query      = ErrorException("cannot run query")

    dbconnect1 = url -> :connected
    dbconnect2 = url -> err_connection

    cursor1 = conn -> :cursor
    cursor2 = conn -> err_cursor

    query1(sql) = cur -> :result
    query2(sql) = cur -> err_query

    url = "localhost:12345"
    sql = "select * from customers"

    @test result(url) |>
            fmap(dbconnect1) |>
            fmap(cursor1) |>
            fmap(query1(sql)) == result(:result)

    @test result(url) |>
            fmap(dbconnect1) |>
            fmap(cursor1) |>
            fmap(query2(sql)) == result(err_query)

    @test result(url) |>
            fmap(dbconnect1) |>
            fmap(cursor2) |>
            fmap(query1(sql)) == result(err_cursor)

    @test result(url) |>
            fmap(dbconnect2) |>
            fmap(cursor1) |>
            fmap(query1(sql)) == result(err_connection)
end

end # top-level testset


text:
"Holds the tableau of a fully implicit Runge-Kutta method."
struct TableauFIRK{T} <: AbstractTableauIRK{T}
    @HeaderTableau

    q::CoefficientsRK{T}

    function TableauFIRK{T}(q) where {T}
        if (q.s > 1 && istrilstrict(q.a)) || (q.s==1 && q.a[1,1] == 0)
            @warn "Initializing TableauFIRK with explicit tableau $(q.name).\nYou might want to use TableauERK instead."
        elseif q.s > 1 && istril(q.a)
            @warn "Initializing TableauFIRK with diagonally implicit tableau $(q.name).\nYou might want to use TableauDIRK instead."
        end

        new(q.name, q.o, q.s, q)
    end
end

function TableauFIRK(q::CoefficientsRK{T}) where {T}
    TableauFIRK{T}(q)
end

function TableauFIRK(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    TableauFIRK{T}(CoefficientsRK(name, order, a, b, c))
end

# TODO function readTableauFIRKFromFile(dir::AbstractString, name::AbstractString)

"Parameters for right-hand side function of fully implicit Runge-Kutta methods."
mutable struct ParametersFIRK{DT, TT, ET <: ODE{DT,TT}, D, S} <: Parameters{DT,TT}
    equ::ET
    tab::TableauFIRK{TT}
    Δt::TT

    t::TT
    q::Vector{DT}
end

function ParametersFIRK(equ::ET, tab::TableauFIRK{TT}, Δt::TT) where {DT, TT, ET <: ODE{DT,TT}}
    ParametersFIRK{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt, 0, zeros(DT, equ.d))
end

"""
Fully implicit Runge-Kutta integrator cache.

### Fields

* `q̃`: initial guess of solution
* `ṽ`: initial guess of vector field
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of solution
* `V`: internal stages of vector field
* `Y`: vector field of internal stages
"""
struct IntegratorCacheFIRK{DT,D,S} <: ODEIntegratorCache{DT,D}
    q̃::Vector{DT}
    ṽ::Vector{DT}
    s̃::Vector{DT}

    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    Y::Vector{Vector{DT}}

    function IntegratorCacheFIRK{DT,D,S}() where {DT,D,S}
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        Y = create_internal_stage_vector(DT, D, S)
        new(zeros(DT,D), zeros(DT,D), zeros(DT,D), Q, V, Y)
    end
end

"Fully implicit Runge-Kutta integrator."
struct IntegratorFIRK{DT, TT, PT <: ParametersFIRK{DT,TT},
                              ST <: NonlinearSolver{DT},
                              IT <: InitialGuessODE{DT,TT}, N, D, S} <: IntegratorRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheFIRK{DT,D,S}
end

function IntegratorFIRK(equation::ODE{DT,TT,FT,N}, tableau::TableauFIRK{TT}, Δt::TT) where {DT,TT,FT,N}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create params
    params = ParametersFIRK(equation, tableau, Δt)

    # create solver
    solver = create_nonlinear_solver(DT, D*S, params)

    # create initial guess
    iguess = InitialGuessODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCacheFIRK{DT,D,S}()

    # create integrator
    IntegratorFIRK{DT, TT, typeof(params), typeof(solver), typeof(iguess), N, D, S}(
                params, solver, iguess, cache)
end

function initialize!(int::IntegratorFIRK, sol::AtomicSolutionODE)
    sol.t̅ = sol.t - timestep(int)

    int.params.equ.v(sol.t, sol.q, sol.v)

    initialize!(int.iguess, sol.t, sol.q, sol.v,
                            sol.t̅, sol.q̅, sol.v̅)
end

function update_params!(int::IntegratorFIRK, sol::AtomicSolutionODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
end

function initial_guess!(int::IntegratorFIRK, sol::AtomicSolutionODE)
    local offset::Int

    # compute initial guess for internal stages
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.v, sol.q̅, sol.v̅, int.cache.q̃, int.cache.ṽ, tableau(int).q.c[i])
        for k in eachindex(int.cache.V[i], int.cache.ṽ)
            int.cache.V[i][k] = int.cache.ṽ[k]
        end
    end
    for i in eachstage(int)
        offset = ndims(int)*(i-1)
        for k in eachdim(int)
            int.solver.x[offset+k] = 0
            for j in eachstage(int)
                int.solver.x[offset+k] += tableau(int).q.a[i,j] * int.cache.V[j][k]
            end
        end
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
                                    params::ParametersFIRK{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}

    local tᵢ::TT

    @assert S == length(Q) == length(V) == length(Y)

    # copy x to Y and compute Q = q + Δt Y
    for i in 1:S
        @assert D == length(Q[i]) == length(V[i]) == length(Y[i])
        for k in 1:D
            Y[i][k] = x[D*(i-1)+k]
            Q[i][k] = params.q[k] + params.Δt * Y[i][k]
        end
    end

    # compute V = v(Q)
    for i in 1:S
        tᵢ = params.t + params.Δt * params.tab.q.c[i]
        params.equ.v(tᵢ, Q[i], V[i])
    end
end

"Compute stages of fully implicit Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersFIRK{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}

    cache = IntegratorCacheFIRK{ST,D,S}()

    quote
        compute_stages!(x, $cache.Q, $cache.V, $cache.Y, params)

        local y1::ST
        local y2::ST

        # compute b = - (Y-AV)
        for i in 1:S
            for k in 1:D
                y1 = 0
                y2 = 0
                for j in 1:S
                    y1 += params.tab.q.a[i,j] * $cache.V[j][k]
                    y2 += params.tab.q.â[i,j] * $cache.V[j][k]
                end
                b[D*(i-1)+k] = - $cache.Y[i][k] + (y1 + y2)
            end
        end
    end
end

"Integrate ODE with fully implicit Runge-Kutta integrator."
function integrate_step!(int::IntegratorFIRK{DT,TT}, sol::AtomicSolutionODE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int, sol)

    # compute initial guess
    initial_guess!(int, sol)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector field at internal stages
    compute_stages!(int.solver.x, int.cache.Q, int.cache.V, int.cache.Y, int.params)

    # compute final update
    update_solution!(sol.q, sol.q̃, int.cache.V, int.params.tab.q.b, int.params.tab.q.b̂, int.params.Δt)

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.v)
end


text:
using Permutations
using DelimitedFiles
using Random

(n, m_pop, m_unp) = (100, 35, 65)

m = m_pop + m_unp

students = vcat(hcat((randperm(m_pop) for i in 1:n)...),
                     m_pop .+ hcat((randperm(m_unp) for i in 1:n)...))
schools = ones(Int64, n, m)

for (i, r) in enumerate(schools)
    if rand() > 0.6
        schools[i] += 1
    end
end

# writedlm("examples/sysopt/students.dat", students)
# writedlm("examples/sysopt/schools.dat", schools)

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
function reconstruct(w::Vector{Float64},c::Vector{Float64},U::Vector{Float64},
    j::Int,k::Float64,σ::Vector{Float64})
    P = function(Pout::Vector{Float64},x::Float64,Pi::Vector{Float64})
        if j == 1 # left boundary polynomial
            σ[1] = U[2] .- U[1];
            σ[2] = U[3] .- U[2];
            Pi[1] = U[1];
            Pi[2] = U[2] .+ σ[1]./k.*(x .- 1.5.*k);
            Pi[3] = 1./c[3].*(-c[1].*U[1] .+ (1-c[2]).*U[2] .- 1./24.*(σ[2] - σ[1]) .+
                (σ[2] .+ σ[1] .- 2.*c[2].*σ[1])./(2.*k).*(x .- 1.5.*k) .+ (σ[2] - σ[1])./(2.*k.^2).*
                (x .- 1.5.*k).^2);
        elseif j == length(U) # right boundary polynomial
            σ[1] = U[end-2] .- U[end-1];
            σ[2] = U[end-1] .- U[end];
            Pi[1] = U[end];
            Pi[2] = U[end-1] .- σ[2]./k.*(x - k.*(length(U)-1.5));
            Pi[3] = 1./c[3].*(-c[1].*U[end] .+ (1-c[2]).*U[end-1] .- 1./24.*(σ[1] - σ[2]) .-
                ((σ[1] .+ σ[2]) - 2.*c[2].*σ[2])./(2.*k).*(x .- k.*(length(U)-1.5)) .+
                (σ[1] - σ[2])./(2.*k.^2).*(x .- k.*(length(U)-1.5)).^2);
        else # interior polynomial
            σ[1] = U[j] .- U[j-1];
            σ[2] = U[j+1] .- U[j];
            Pi[1] = U[j] .+ σ[1]./k.*(x .- k.*(j-0.5));
            Pi[2] = 1./c[2].*((1.-c[1].-c[3]).*U[j] .- 1./24.*(σ[2] .- σ[1]) .+
                ((1 - 2.*c[3]).*U[j+1] .+ 2.*(c[3] .- c[1]).*U[j] .+
                (2.*c[1] .- 1).*U[j-1])./(2.*k).*(x .- k.*(j-0.5)) .+
                (σ[2] .- σ[1])./(2.*k.^2).*(x .- k.*(j-0.5)).^2);
            Pi[3] = U[j] .+ σ[2]./k.*(x .- k.*(j-0.5));
        end
        Pout[1] = dot(w,Pi);
    end
    return P
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
function test_model(m::Model)
    @test length(m.cells) <= m.cfg["cells_max"]
    nconns = 0
    for cell in m.cells
        nconns += length(cell.conns)
        @test cell.input >= -1.0
        @test cell.input <= 1.0
        @test cell.output >= -1.0
        @test cell.output <= 1.0
        @test all(cell.state .>= -1.0)
        @test all(cell.state .<= 1.0)
        @test length(cell.state) == m.cfg["n_cell_state"]
        @test all(cell.params .>= -1.0)
        @test all(cell.params .<= 1.0)
        @test length(cell.params) == m.cfg["n_cell_params"]
        for conn in cell.conns
            @test conn.input >= -1.0
            @test conn.input <= 1.0
            @test all(conn.state .>= -1.0)
            @test all(conn.state .<= 1.0)
            @test length(conn.state) == m.cfg["n_conn_state"]
            @test all(conn.params .>= -1.0)
            @test all(conn.params .<= 1.0)
            @test length(conn.params) == m.cfg["n_conn_params"]
            @test conn.source[] == cell
            if typeof(conn.dest[]) == Cell
                @test conn.dest[] in m.cells
            else
                in_conns = false
                for c2 in m.cells
                    if conn.dest[] in c2.conns
                        in_conns = true
                        break
                    end
                end
                @test in_conns
            end
        end
    end
    @test nconns <= m.cfg["conns_max"]
end

@testset "Initialization tests" begin
    cfg = Config()
    c = Controller(cfg)
    nin = rand(5:10)
    nout = rand(5:10)

    @testset "Layered" begin
        m = Model(cfg, c)
        layered_init!(m, nin, nout)
        test_model(m)
    end

    @testset "Layered single cell" begin
        m = Model(cfg, c)
        layered_init!(m, nin, nout; nhidden=1)
        test_model(m)
    end

    @testset "Random cell" begin
        m = Model(cfg, c)
        random_init!(m, nin, nout)
        test_model(m)
    end

    @testset "Random single cell" begin
        m = Model(cfg, c)
        random_init!(m, nin, nout; nhidden=1)
        test_model(m)
    end
end

@testset "Model tests" begin
    cfg = Config()
    c = Controller(cfg)
    m = Model(cfg, c)
    nin = rand(5:10)
    nout = rand(5:10)
    nhidden = rand(5:10)
    nreward = rand(5:10)

    @testset "Initialize" begin
        random_init!(m, nin, nout; nhidden=nhidden, nreward=nreward)
        test_model(m)
        @test length(m.cells) == nin + nout + nhidden + nreward
        @test length(m.inputs) == nin
        @test length(m.outputs) == nout
    end

    @testset "Input" begin
        set_input!(m, rand(nin))
        test_model(m)
    end

    @testset "Step" begin
        for i in 1:5
            step!(m)
            @test m.state[1] == i
            test_model(m)
        end
    end

    @testset "Output" begin
        outputs = get_output(m)
        test_model(m)
        @test all(outputs .>= -1.0)
        @test all(outputs .<= 1.0)
        @test length(outputs) == nout
    end

    @testset "Reward" begin
        reward!(m, rand(nreward))
        test_model(m)
    end

    @testset "Full step" begin
        for i in 1:5
            println("Base controller step $i")
            @time outputs = step!(m, rand(nin))
            reward!(m, rand(nreward))
            test_model(m)
            @test all(outputs .>= -1.0)
            @test all(outputs .<= 1.0)
            @test length(outputs) == nout
        end
    end
end

function controller_test(cont_constructor::Function; cfg::Dict=Config())
    c = cont_constructor(cfg)
    m = Model(cfg, c)
    nin = rand(5:10)
    nout = rand(5:10)
    random_init!(m, nin, nout; nreward=nout)

    @testset "Full step" begin
        for i in 1:5
            println(string(cont_constructor, " step ", i))
            @time outputs = step!(m, rand(nin))
            reward!(m, rand(nout))
            test_model(m)
            @test all(outputs .>= -1.0)
            @test all(outputs .<= 1.0)
            @test length(outputs) == nout
        end
    end
end

@testset "Default controller" begin
    controller_test(x->Controller(x))
end

@testset "Const controller" begin
    controller_test(const_controller)
end

@testset "Random controller" begin
    controller_test(rand_controller)
end

@testset "Static controller" begin
    controller_test(static_controller;
                    cfg = Config(Config()))
end

# @testset "SNN controller" begin
#     controller_test(snn_controller;
#                     cfg = Config(Config(), "cfg/snn.yaml"))
# end

@testset "CGP controller" begin
    controller_test(cgp_controller)
end


text:
@defcomp emissions begin
    regions = Index()

    E = Variable(index=[time]) # Total CO2 emissions (GtC per year)
    EIND = Variable(index=[time, regions]) # Industrial emissions (GtC per year)
    CCA = Variable(index=[time]) # Cumulative indiustrial emissions
    ABATECOST = Variable(index=[time, regions]) # Cost of emissions reductions  (trillions 2005 USD per year)
    MCABATE = Variable(index=[time, regions]) # Marginal cost of abatement (2005$ per ton CO2)
    CPRICE = Variable(index=[time, regions]) # Carbon price (2005$ per ton of CO2)

    sigma = Parameter(index=[time, regions]) # CO2-equivalent-emissions output ratio
    YGROSS = Parameter(index=[time, regions]) # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
    etree = Parameter(index=[time]) # Emissions from deforestation
    cost1 = Parameter(index=[time, regions]) # Adjusted cost for backstop
    expcost2 = Parameter(index=[regions]) # Exponent of control cost function
    partfract = Parameter(index=[time, regions]) # Fraction of emissions in control regime
    pbacktime = Parameter(index=[time, regions]) # Backstop price
    MIU = Parameter(index=[time, regions]) # Emission control rate GHGs

    function run_timestep(p, v, d, t)

        #Define function for EIND
        for r in d.regions
            v.EIND[t,r] = p.sigma[t,r] * p.YGROSS[t,r] * (1-p.MIU[t,r])
        end

        #Define function for E
        v.E[t] = sum(v.EIND[t,:]) + p.etree[t]

        #Define function for CCA
        if is_first(t)
            v.CCA[t] = sum(v.EIND[t,:]) * 10.
        else
            v.CCA[t] =  v.CCA[t-1] + (sum(v.EIND[t,:]) * 10.)
        end

        #Define function for ABATECOST
        for r in d.regions
        if t.t==1
            v.ABATECOST[t,r] = p.YGROSS[t,r] * p.cost1[t,r] * (p.MIU[t,r]^p.expcost2[r]) * (p.partfract[t,r]^(1 - p.expcost2[r]))
        else
            v.ABATECOST[t,r] = p.YGROSS[t,r] * p.cost1[t,r] * (p.MIU[t,r]^p.expcost2[r] + 50*(p.MIU[t,r] - p.MIU[t-1,r])^4) * (p.partfract[t,r]^(1 - p.expcost2[r]))
        end

        #Define function for MCABATE
        for r in d.regions
            v.MCABATE[t,r] = p.pbacktime[t,r] * p.MIU[t,r]^(p.expcost2[r] - 1)
        end

        #Define function for CPRICE
        for r in d.regions
            v.CPRICE[t,r] = p.pbacktime[t,r] * 1000 * p.MIU[t,r]^(p.expcost2[r] - 1)
        end
    end
end
end

max_stars_repo_name: Environment-Research/Utilitarianism

max_issues_repo_name: Environment-Research/Utilitarianism

max_forks_repo_name: Environment-Research/Utilitarianism


text:
using Test

@testset "Create tests for your project" begin

@testset "Testing linear algebra results" begin
   @test 1 == 1
end;

end; # ex1.jl


text:
# Info about computers code is supposed to run on
#=
Assumes that all computers only differ in a single base dir, such as '/Users/lutz'
=#
@enum Computer ComputerLocal = 1 ComputerLongleaf = 2

const localBaseDir = "/Users/lutz"
const longleafBaseDir = "/nas/longleaf"

# Computer we are running on
function current_computer()
    if isdir(localBaseDir)
        return ComputerLocal
    elseif isdir(longleafBaseDir)
        return ComputerLongleaf
    else
        error("Invalid computer")
    end
end

# Are we on the local computer?
function running_local()
    ci = computer_info()
    return ci.runLocal
end

# Container for info about computers that code may run on
mutable struct ComputerInfo
    compId  ::  Computer
    compName :: String
    runLocal  ::  Bool
    # Base directory that everything hangs off
    # Local: user home dir. Remote: prepend .../lhendri/
    baseDir :: String

    # Incomplete initialization is OK
    ComputerInfo(compId) = new(compId)
end
#

# Settings for each computer are collected here
function computer_info(compId  ::  Computer)
    cS = ComputerInfo(compId)
    cS.runLocal = false

    if compId == ComputerLocal
        cS.compName = "Local"
        cS.baseDir = localBaseDir;
        cS.runLocal = true
    elseif compId == ComputerLongleaf
        cS.compName = "Longleaf"
        cS.baseDir = longleafBaseDir;
    else
        error("invalid")
    end

    return cS
end

function computer_info()
    return computer_info(current_computer())
end

function docu_dir(compId :: Computer)
    ci = computer_info(compId)
    return joinpath(ci.baseDir, "Documents")
end

function docu_dir()
    comp = current_computer()
    return docu_dir(comp)
end

## Dirs hanging off Documents (at least on local)

# Repos live here
function julia_dir(compId :: Computer)
	ci = computer_info(compId);
	if ci.runLocal
	    return joinpath(docu_dir(compId), "julia")
	else
	    return joinpath(homedir(), "julia")
	end
end

function julia_dir()
    comp = current_computer()
    return julia_dir(comp)
end

function project_dir(year1)
    joinpath(docu_dir(), "projects", "p$year1");
end

## Dirs hanging off Julia

function shared_dir(compId :: Computer)
    return joinpath(julia_dir(compId), "shared")
end

function shared_dir()
    comp = current_computer()
    return shared_dir(comp)
end

function test_dir()
    return joinpath(julia_dir(), "test_files")
end


text:
using Documenter, ScenTrees

#const ASSETS = readdir(joinpath(@__DIR__, "src", "assets"))

isCI = get(ENV, "CI", nothing) == "true" #Travis populates this env variable by default

makedocs(
	sitename =  "ScenTrees.jl",
	authors = "Kipngeno Kirui",
	clean = true,
	doctest = true,
	format = Documenter.HTML(prettyurls = isCI),
	pages = ["Home" => "index.md",
		"Tutorials" => Any["tutorial/tutorial1.md",
				    "tutorial/tutorial2.md",
				    "tutorial/tutorial3.md",
				    "tutorial/tutorial31.md",
				    "tutorial/tutorial4.md",
				    "tutorial/tutorial41.md",
				    "tutorial/tutorial5.md"
				]
		]
)

if isCI
    deploydocs(repo = "github.com/kirui93/ScenTrees.jl.git")
end


text:
using HorizonSideRobots
originRobot = Robot(animate = false, "24example.sit")
show!(originRobot)
include("../AbstractType.jl")
include("../functions.jl")

function square(robot::Robot, moveSide::HorizonSide)::Integer
    square = 0
    height = 0
    width = 0
    stepsToBack = 0
    reqBorder = anti_side(moveSide)
    # можно было через while сделать: задать итератор = 4, и при каждом проходе по циклу менять значение side на next_side_pr(side)
    for side in [reqBorder, next_side_pr(reqBorder), next_side_pr(next_side_pr(reqBorder)), anti_side(next_side_pr(reqBorder))]
        while isborder(robot, side)
            if (side == next_side_pr(reqBorder))
                width += 1
            elseif (side == next_side_pr(next_side_pr(reqBorder)))
                height += 1
            elseif (side == reqBorder)
                stepsToBack += 1
            end
            move!(robot, next_side(side))
        end
        move!(robot, side)
    end
    for _ in 1:(height - stepsToBack)
        move!(robot, next_side_pr(moveSide))
    end
    return height * width
end

function master!(robot::Robot)::Integer
    maxSquare = 0
    tempReqSquare = 0
    squares = Set()
    while !(isborder(robot, Nord))
        for side in [Ost, West]
            while true
                if isborder(robot, side)
                    if try_move!(robot, side)
                        tempReqSquare = square(robot, side)
                        push!(squares, tempReqSquare)
                        if tempReqSquare > maxSquare
                            maxSquare = tempReqSquare
                            # можно написать очень похожую на square функцию, которая бы ещё и обвела этит прямоугольник
                        end
                    else
                        break
                    end
                else
                    move!(robot, side)
                end
            end
            if !(isborder(robot, Nord))
                move!(robot, Nord)
            end
        end
    end
    println(squares)
    return maxSquare
end

master!(originRobot)

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
__precompile__()

module CompositeFieldVectors

using Requires

@require Unitful begin
    using Unitful
end

import Base: getindex, setindex!, endof, size, vec, length, copy # read, read!, write

export CompositeFieldVector,
       getindex,
       setindex!,
       endof,
       size,
       vec,
       length,
       copy,
       compositefieldparent,
       compositefieldparents,
       compositefieldname,
       compositefieldnames

include("constructors.jl")
include("methods.jl")

end # module


text:
using ApproxFun, SingularIntegralEquations, DualNumbers, Holomorphic

Disk()

points(Circle()

x=Fun()
w=sqrt(1-x^2)
g=stieltjes(w)
g/(-2π*im)
g|(2Circle())
imag(g|+)
g|-
∂(g)

domain(g)

-∂
cauchy(+1

f(.1+.2im)

Γ=ComplexPlane()\Interval()
S=Space(Γ)

f(z)=stieltjes(w,z)

g=Fun(Fun(x->f(x+0im)-f(x-0im),S.space).coefficients/(-2π*im),S)

g(.1+.2im)
f(.1+.2im)

stieltjes(w).coefficients

g.coefficients

stieltjes(1/w,1.0001)

d=Interval(-1.-1im,1+1.im)

fromcanonical(d,0.1)-fromcanonical(d,0.1+0.im)
fromcanonical(d,0.1)

fromcanonical(d,0.-0.im)
fromcanonical(d,0.+0.im)

using DualNumbers

fromcanonical(d,dual(0.,1.im))
fromcanonical(d,dual(0.,-1.im))

Dual(a::Real,b::Complex)=Dual(a+0*b,b)

sqrt(Dual(-1.,1.im))
sqrt(Dual(

sqrt(-1.-0.im)
sqrt(dual(-1.,-1im))
sqrt(dual(-1.+0.im,+1im))
sqrt(dual(-1.-0.im,-1im))

realpart(sqrt(dual(-1.,-im)))
realpart(sqrt(dual(-1.,im)))

f(z)=sqrt((1.+1.im)*z)

realpart(f(dual(-1/(1.+1.im),-1.im)))

f(-1/(1.+1.im)-0.000001im)

sqrt(-1.)
Pkg.checkout("DualNumbers")


text:
# ---
# title: 43. Multiply Strings
# id: problem43
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, String
# link: <https://leetcode.com/problems/multiply-strings/description/>
# hidden: true
# ---
#
# Given two non-negative integers `num1` and `num2` represented as strings,
# return the product of `num1` and `num2`, also represented as a string.
#
# **Note:**  You must not use any built-in BigInteger library or convert the
# inputs to integer directly.
#
#
#
# **Example 1:**
#
#
#
#     Input: num1 = "2", num2 = "3"
#     Output: "6"
#
#
# **Example 2:**
#
#
#
#     Input: num1 = "123", num2 = "456"
#     Output: "56088"
#
#
#
#
# **Constraints:**
#
#   * `1 <= num1.length, num2.length <= 200`
#   * `num1` and `num2` consist of digits only.
#   * Both `num1` and `num2` do not contain any leading zero, except the number `0` itself.
#
#
## @lc code=start
using LeetCode

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


text:
# Test the extensions of Base.convert
@testset "Convert Extensions" begin
    m = Model()
    @variable(m, x[1:2])
    # GenericAffExpr -> GenericAffExpr
    @testset "Aff->Aff" begin
        aff = x[1] - 3 * x[2] - 2
        @test convert(GenericAffExpr{Float64, VariableRef}, aff) == aff
        @test typeof(convert(GenericAffExpr{Float64, AbstractVariableRef},
                           aff)) == GenericAffExpr{Float64, AbstractVariableRef}
        @test convert(GenericAffExpr{Float64, AbstractVariableRef},
                      aff).terms == aff.terms
        @test convert(GenericAffExpr{Float64, AbstractVariableRef},
                      aff).constant == aff.constant
    end
    # GenericQuadExpr -> GenericQuadExpr
    @testset "Quad->Quad" begin
        quad = 2 * x[1] * x[2] + x[1] - 1
        @test convert(GenericQuadExpr{Float64, VariableRef}, quad) == quad
        @test typeof(convert(GenericQuadExpr{Float64, AbstractVariableRef},
                         quad)) == GenericQuadExpr{Float64, AbstractVariableRef}
        @test convert(GenericQuadExpr{Float64, AbstractVariableRef},
                      quad).aff.terms == quad.aff.terms
        @test convert(GenericQuadExpr{Float64, AbstractVariableRef},
                      quad).aff.constant == quad.aff.constant
        @test convert(GenericQuadExpr{Float64, AbstractVariableRef},
                      quad).terms == quad.terms
    end
    # UnorderedPair -> UnorderedPair
    @testset "UoPair->UoPair" begin
        pair = UnorderedPair(x[1], x[2])
        @test convert(UnorderedPair{VariableRef}, pair) == pair
        @test typeof(convert(UnorderedPair{AbstractVariableRef},
                                    pair)) == UnorderedPair{AbstractVariableRef}
        @test convert(UnorderedPair{AbstractVariableRef}, pair).a == pair.a
        @test convert(UnorderedPair{AbstractVariableRef}, pair).b == pair.b
    end
    # Array -> SparseAxisArray
    @testset "Array->Sparse" begin
        @test convert(JuMPC.SparseAxisArray,
                      x) isa JuMPC.SparseAxisArray
        @test length(convert(JuMPC.SparseAxisArray, x)) == length(x)
        inds = CartesianIndices(x)
        @test convert(JuMPC.SparseAxisArray,
                      x)[inds[1][1]] == x[inds[1]]
        @test convert(JuMPC.SparseAxisArray,
                      x)[inds[2][1]] == x[inds[2]]
    end
    # DenseAxisArray -> SparseAxisArray
    @testset "DenseArray->Sparse" begin
        @variable(m, y[2:3])
        @test convert(JuMPC.SparseAxisArray,
                      y) isa JuMPC.SparseAxisArray
        @test length(convert(JuMPC.SparseAxisArray, y)) == length(y)
        inds = collect(keys(y))
        @test convert(JuMPC.SparseAxisArray,
                      y)[inds[1][1]] == y[inds[1]]
        @test convert(JuMPC.SparseAxisArray,
                      y)[inds[2][1]] == y[inds[2]]
    end
end

# Test ParameterBounds extensions
@testset "Parameter Bounds" begin
    # setup
    m = InfiniteModel()
    m.params[1] = InfOptParameter(IntervalSet(0, 10), Number[], false)
    m.params[2] = InfOptParameter(IntervalSet(0, 10), Number[], false)
    m.params[3] = InfOptParameter(IntervalSet(0, 10), Number[], false)
    par = ParameterRef(m, 1)
    pars = [ParameterRef(m, 2), ParameterRef(m, 3)]
    # test length
    @testset "Base.length" begin
        @test length(ParameterBounds()) == 0
        @test length(ParameterBounds(Dict(par => IntervalSet(0,0)))) == 1
    end
    # test :(==)
    @testset "Base.:(==)" begin
        dict = Dict(par => IntervalSet(0,0))
        @test ParameterBounds() == ParameterBounds()
        @test ParameterBounds() != ParameterBounds(dict)
        @test ParameterBounds(dict) == ParameterBounds(dict)
    end
    # test copy
    @testset "Base.copy" begin
        dict = Dict(par => IntervalSet(0,0))
        @test ParameterBounds(dict) == copy(ParameterBounds(dict))
    end
end

# Test extension of keys
@testset "Base.keys" begin
    m = Model()
    @variable(m, x[1:2], container = SparseAxisArray)
    @test keys(x) == keys(x.data)
end

# Test extension of isapprox
@testset "Base.isapprox" begin
    a = JuMPC.SparseAxisArray(Dict((1,)=>-0, (3,)=>1.1))
    b = JuMPC.SparseAxisArray(Dict((1,)=>-0, (3,)=>1.1 + 1e-30))
    @test isapprox(a, b)
end

# Test the possible_convert functions
@testset "_possible_convert" begin
    m = Model()
    @variable(m, x[1:2])
    # GenericAffExpr conversions
    @testset "GenericAffExpr" begin
        aff = x[1] - 3 * x[2] - 2
        @test typeof(InfiniteOpt._possible_convert(AbstractVariableRef,
                           aff)) == GenericAffExpr{Float64, AbstractVariableRef}
        @test InfiniteOpt._possible_convert(AbstractVariableRef,
                                            aff).terms == aff.terms
        @test InfiniteOpt._possible_convert(AbstractVariableRef,
                                            aff).constant == aff.constant
        @test typeof(InfiniteOpt._possible_convert(GeneralVariableRef,
                                   aff)) == GenericAffExpr{Float64, VariableRef}
    end
    # GenericQuadExpr conversions
    @testset "GenericQuadExpr" begin
        quad = 2 * x[1] * x[2] + x[1] - 1
        @test typeof(InfiniteOpt._possible_convert(AbstractVariableRef,
                         quad)) == GenericQuadExpr{Float64, AbstractVariableRef}
        @test InfiniteOpt._possible_convert(AbstractVariableRef,
                                            quad).aff.terms == quad.aff.terms
        @test InfiniteOpt._possible_convert(AbstractVariableRef,
                                         quad).aff.constant == quad.aff.constant
        @test InfiniteOpt._possible_convert(AbstractVariableRef,
                                            quad).terms == quad.terms
        @test typeof(InfiniteOpt._possible_convert(GeneralVariableRef,
                                 quad)) == GenericQuadExpr{Float64, VariableRef}
    end
end

# Test the _make_vector functions
@testset "_make_vector" begin
    x = [2, 4, 1, 8]
    # test with JuMP containers
    @testset "AbstractArray" begin
        y = convert(JuMPC.SparseAxisArray, x)
        @test sort(InfiniteOpt._make_vector(y)) == [1, 2, 4, 8]
        y = JuMPC.DenseAxisArray(x, [1, 2, 3, 4])
        @test InfiniteOpt._make_vector(y) == x
    end
    # test arrays
    @testset "Array" begin
        @test InfiniteOpt._make_vector(x) == x
        @test InfiniteOpt._make_vector(ones(2, 2)) == ones(2, 2)
    end
    # test other stuff
    @testset "Non-Array" begin
        @test InfiniteOpt._make_vector(2) == 2
    end
end


text:
using Plots, HyperbolicPlane, SimpleDrawing

function tritest(X::HPolygon)
    set_thickness(X,3)
    set_color(X,:blue)
    TT = triangulate(X)
    plot()
    draw(TT)
    draw(X)
    for p in X.plist
        a = HPoint(p)
        set_color(a,:red)
        set_radius(a,1.5)
        draw(a)
    end
    finish()
end


text:
using GLMakie
# using CairoMakie
GLMakie.activate!()
using Makie
using JLD
using Formatting

import Optim
using TopOpt
using TopOpt.TrussTopOptProblems.TrussVisualization: visualize
Nonconvex.@load Ipopt
Nonconvex.@load NLopt

include("problem_defs.jl")
include("../utils.jl")
using .DeflateTruss

RESULT_DIR = joinpath(@__DIR__, "results");

function optimize_truss(problem_name, opt_task; verbose=false, write=false, optimizer="mma", distance="l2",
    deflation_iters=5)
    result_data = Dict()
    objs = Float64[]

    E = 1.0 # Young’s modulus
    v = 0.3 # Poisson’s ratio
    size_ratio = 1
    if problem_name == "dense_graph"
        f = [0.0, 100.0]
        problem = CustomPointLoadCantileverTruss((40*size_ratio,10*size_ratio),Tuple(ones(2)),E,v,f)
    elseif problem == "tim"
        node_points, elements, mats, crosssecs, fixities, load_cases = load_truss_json(
            joinpath(@__DIR__, "tim_2d.json")
        );
        loads = load_cases["0"]
        f = "load case 0"
        problem = TrussProblem(
            Val{:Linear}, node_points, elements, loads, fixities, mats, crosssecs
        );
    else
        error("Unsupported problem $problem_name")
    end

    problem_config = "$(problem_name)_$(opt_task)_$(optimizer)_$distance"
    problem_result_dir = joinpath(RESULT_DIR, problem_config)
    if !ispath(problem_result_dir)
        mkdir(problem_result_dir)
    end

    V = 0.3 # volume fraction
    xmin = 0.0001 # minimum density
    penalty_val = 3.0

    penalty = TopOpt.PowerPenalty(penalty_val)
    solver = FEASolver(Direct, problem; xmin=xmin, penalty=penalty)
    solver()
    result_data["config"] = Dict("E" => E, "ν" => v, "f" => f, "penalty" => penalty_val, "xmin" => xmin, "vol_fraction" => V)

    # * comliance minimization objective
    if occursin("vol_constrained", opt_task) && occursin("min_compliance", opt_task)
        comp = TopOpt.Compliance(problem, solver)
        obj = comp
        volfrac = TopOpt.Volume(problem, solver)
        constr = x -> volfrac(x) - V
        # or equivalently:
        # constr = x -> sum(x) / length(x) - V
    else
        error("Undefined task $(opt_task)")
    end

    x0 = fill(V, length(solver.vars))
    nelem = length(x0)
    println("#elements : $nelem")

    Nonconvex.NonconvexCore.show_residuals[] = verbose

    m = Model(obj)
    addvar!(m, zeros(nelem), ones(nelem))
    add_ineq_constraint!(m, constr)

    if optimizer == "mma"
        alg = MMA02()
        options = MMAOptions(maxiter=100,
            tol=Nonconvex.Tolerance(kkt=0.001),
            # dual_options=Optim.Options(iterations = 100)
            )
    elseif optimizer == "nlopt"
        alg = NLoptAlg(:LD_MMA)
        options = NLoptOptions()
    elseif optimizer == "ipopt"
        alg = IpoptAlg()
        options = IpoptOptions(tol = 1e-4, max_iter=1000)
    else
        error("Undefined optimizer $optimizer")
    end
    convcriteria = Nonconvex.KKTCriteria()

    println("Initial solve:")
    st_time = time()
    r1 = Nonconvex.optimize(m, alg, x0; options=options, convcriteria = convcriteria)
    runtime = time() - st_time
    printfmt("Minimum: {:.3f}; Constraint: {:.3f}; runtime {:.3f}\n", r1.minimum, maximum(constr(r1.minimizer)), runtime)

    fig1 = visualize(problem, topology=r1.minimizer)
    hidedecorations!(fig1.current_axis.x)
    push!(objs, r1.minimum)

    if write
        result_data["init_solve"] = Dict("minimizer" => r1.minimizer, "minimum" => r1.minimum, "runtime" => runtime)
        save(joinpath(problem_result_dir, "r1.png"), fig1)
        # save(joinpath(problem_result_dir, "r1.pdf"), fig1, pt_per_unit = 1)
    else
        display(fig1)
        wait_for_key("Enter to continue...")
    end

    if endswith(opt_task, "deflation")
        power = 4.0
        radius = 20.0
        y_upperbound = 1e2
        result_data["deflation_config"] = Dict("power" => power, "radius" => radius, "dist" => distance,
            "y_upperbound" => y_upperbound)

        if distance == "l2"
            dist_fn = (x,y) -> norm(x-y, 2)
        else
            error("Unimplemented distance $distance")
        end

        solutions = [r1.minimizer]
        for iter in 1:deflation_iters
            println("-"^10)
            println("Iter - $iter")
            function deflation_constr(X::Vector)
                d = zero(eltype(X))
                for sol in solutions
                    d += max(dist_fn(X[1:end-1], sol) - radius, 0)^(-power)
                    # d += dist_fn(X[1:end-1], sol)^(-power) + shift
                end
                println("d = $d, y = $(X[end])")
                return d - X[end]
            end

            # * remake problem because dim changes by 1
            df_obj = x -> obj(x[1:end-1])
            df_constr = x -> constr(x[1:end-1])

            m = Model(df_obj)
            addvar!(m, zeros(nelem), ones(nelem))
            # * deflation slack variable y
            addvar!(m, [0], [y_upperbound])
            add_ineq_constraint!(m, df_constr)
            add_ineq_constraint!(m, deflation_constr)

            df_alg = alg
            df_options = options
            df_convcriteria = convcriteria

            st_time = time()
            df_result = optimize(m, df_alg, vcat(x0, 0.0), options = df_options, convcriteria = df_convcriteria)
            runtime = time() - st_time

            printfmt("Minimum: {:.3f}; Constraint: {:.3f}; DF Constraint: {:.3f}; runtime {:.3f}\n",
                df_result.minimum, maximum(df_constr(df_result.minimizer)),
                maximum(deflation_constr(df_result.minimizer)), runtime)
            final_dist = dist_fn(r1.minimizer, df_result.minimizer[1:end-1])
            println("Distance: $(final_dist)")
            push!(objs, df_result.minimum)

            fig2 = visualize(problem, topology=df_result.minimizer[1:end-1])
            hidedecorations!(fig2.current_axis.x)

            if write
                result_data["deflate_$(iter)"] = Dict("minimizer" => df_result.minimizer, "minimum" => df_result.minimum, "runtime" => runtime, "distance_to_xstar" => final_dist)
                save(joinpath(problem_result_dir, "deflate_r$(iter).png"), fig2)
                # save(joinpath(problem_result_dir, "deflate_r$(iter).pdf"), fig2, pt_per_unit = 1)
            else
                display(fig2)
                wait_for_key("Enter to continue...")
            end
            push!(solutions, df_result.minimizer[1:end-1])
        end

        if write
            obj_fig = Figure(resolution=(800,800), font="Arial")
            ax = Axis(obj_fig[1,1])
            ax.xlabel = "deflation iterations"
            ax.ylabel = "objectives"
            lines!(ax, 0:1:length(objs)-1, objs)
            # save(joinpath(problem_result_dir, "_obj_history.pdf"), obj_fig)
            save(joinpath(problem_result_dir, "_obj_history.png"), obj_fig)
        end
    end # end deflation loop

    if write
        save(joinpath(problem_result_dir, "data.jld"), "result_data", result_data)
    end
end

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

using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--problem"
            help = "problem to run"
            arg_type = String
            default = "dense_graph" # tim
        "--task"
            help = "problem formulation"
            arg_type = String
            default = "min_compliance_vol_constrained_deflation" # _deflation, _buckling_constrained
        "--optimizer"
            help = "optimizer choice"
            arg_type = String
            default = "nlopt" # mma, nlopt, ipopt
        "--distance"
            help = "distance measure choice"
            arg_type = String
            default = "l2" # kl, w2
        "--deflate_iters"
            help = "print extra info."
            arg_type = Int
            default = 5
        "--verbose"
            help = "print extra info."
            action = :store_true
        "--write"
            help = "export result."
            action = :store_true
        "--replot"
            help = "parse and replot saved results."
            action = :store_true
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println("Parsed args:")
    println(parsed_args)
    println("="^10)

    optimize_truss(parsed_args["problem"], parsed_args["task"], verbose=parsed_args["verbose"],
        write=parsed_args["write"], optimizer=parsed_args["optimizer"], distance=parsed_args["distance"],
        deflation_iters=parsed_args["deflate_iters"])
end

main()


text:
@testset "test fft_helpers" begin

    @testset "Optional collect" begin
        y = [1,2,3]
        x = fftshift_view(y, (1))
        @test fftshift(y) == FourierTools.optional_collect(x)
    end

    @testset "Test ft and ift wrappers" begin
        Random.seed!(42)
        testft(arr, dims) = @test(ft(arr, dims) ≈ fftshift(fft(ifftshift(arr, dims), dims), dims))
        testift(arr, dims) = @test(ift(arr, dims) ≈ fftshift(ifft(ifftshift(arr, dims), dims), dims))
        testffts(arr, dims) = @test(ffts(arr, dims) ≈ fftshift(fft(arr, dims), dims))
        testiffts(arr, dims) = @test(iffts(arr, dims) ≈ ifft(ifftshift(arr, dims), dims))
        testrft(arr, dims) = @test(rffts(arr, dims) ≈ fftshift(rfft(arr, dims), dims[2:end]))
        testirft(arr, dims, d) = @test(irffts(arr, d, dims) ≈ irfft(ifftshift(arr, dims[2:end]), d, dims))
        for dim = 1:4
            for _ in 1:3
                s = ntuple(_ -> rand(1:13), dim)
                arr = randn(ComplexF32, s)
                dims = 1:dim
                testft(arr, dims)
                testift(arr, dims)
                dims = 1:rand(1:dim)
                testft(arr, dims)
                testift(arr, dims)
                testffts(arr, dims)
                testiffts(arr, dims)

            end
        end
    end

    @testset "Test 2d fft helpers" begin
        arr = randn((6,7,8))
        dims = [1,2]
        d = 6
        @test(ft2d(arr) == fftshift(fft(ifftshift(arr, (1,2)), (1,2)), dims))
        @test(ift2d(arr) == fftshift(ifft(ifftshift(arr, (1,2)), (1,2)), dims))
        @test(ffts2d(arr) == fftshift(fft(arr, (1,2)), (1,2)))
        @test(iffts2d(arr) == ifft(ifftshift(arr, (1,2)), (1,2)))
        @test(rffts2d(arr) == fftshift(rfft(arr, (1,2)), dims[2:2]))
        @test(rft2d(arr) == fftshift(rfft(ifftshift(arr, (1,2)), (1,2)), dims[2:2]))
        @test(fft2d(arr) == fft(arr, dims))
        @test(ifft2d(arr) == ifft(arr, dims))
        @test(rfft2d(arr) == rfft(arr, (1,2)))
        @test(fftshift2d(arr) == fftshift(arr, (1,2)))
        @test(ifftshift2d(arr) == ifftshift(arr, (1,2)))
        @test(fftshift2d_view(arr) == fftshift_view(arr, (1,2)))
        @test(ifftshift2d_view(arr) == ifftshift_view(arr, (1,2)))

        arr = randn(ComplexF32, (4,7,8))
        @test(irffts2d(arr, d) == irfft(ifftshift(arr, dims[2:2]), d, (1,2)))
        @test(irft2d(arr, d) == irft(arr, d, (1,2)))
        @test(irfft2d(arr, d) == irfft(arr, d, (1,2)))
    end

    @testset "Test ft, ift, rft and irft real space centering" begin
        szs = ((10,10),(11,10),(100,101),(101,101))
        for sz in szs
            @test ft(ones(sz)) ≈ prod(sz) .* delta(sz)
            @test ft(delta(sz)) ≈ ones(sz)
            @test rft(ones(sz)) ≈ prod(sz) .* delta(rft_size(sz), offset=CtrRFT)
            @test rft(delta(sz)) ≈ ones(rft_size(sz))
            @test ift(ones(sz)) ≈ delta(sz)
            @test ift(delta(sz)) ≈ ones(sz) ./ prod(sz)
            @test irft(ones(rft_size(sz)),sz[1]) ≈ delta(sz)
            @test irft(delta(rft_size(sz),offset=CtrRFT),sz[1]) ≈ ones(sz) ./ prod(sz)
        end
    end

    @testset "Test in place methods" begin
        x = randn(ComplexF32, (5,3,10))
        dims = (1,2)
        @test fftshift(fft(x, dims), dims) ≈ ffts!(copy(x), dims)
        @test ffts2d!(copy(x)) ≈ ffts!(copy(x), (1,2))
    end

end


text:
module SossBase

import Base.rand
using NamedTupleTools

include("core/statement.jl")
include("core/model.jl")
include("core/jointdistribution.jl")
include("core/canonical.jl")
include("core/dependencies.jl")
include("core/toposort.jl")
include("core/utils.jl")

include("primitives/rand.jl")
include("primitives/logdensity.jl")
include("primitives/likelihood-weighting.jl")

include("transforms/predictive.jl")
include("transforms/markovblanket.jl")

end # module


text:
#=================================
Methods using Taylor Models
=================================#
using TaylorModels

@inline zeroBox(N) = IntervalBox(0..0, N)
@inline symBox(N) = IntervalBox(-1..1, N)

"""
    enclose_TaylorModels(f::Function, dom::Interval_or_IntervalBox;
                         order::Int=10, normalize::Bool=true)

Compute a range enclosure using Taylor Models substitution.

### Input

- `f`         -- function
- `dom`       -- hyperrectangular domain, either a unidimensional  `Interval` or
                 a multidimensional `IntervalBox`
- `order`     -- (optional, default: `10`) order of the taylor model used to compute
                 an enclosure of `f` over `dom`
- `normalize` -- (optional, default: `true`) if `true`, normalize the taylor model
                 on the unit symmetric box around the origin

### Output

An interval representing the range enclosure (minimum and maximum) of `f` over
its domain `dom`.

### Algorithm

The result is obtained by an evaluation of the function over Taylor model variables
after recentering in the given domain and (optionally) normalizing in the
symmetric ``[-1..1]^n`` domain.

We refer to the documentation and source code of `TaylorModels` for details.
"""
function enclose_TaylorModels(f::Function, dom::Interval_or_IntervalBox;
                              order::Int=10, normalize::Bool=true)
    if normalize
        R = _enclose_TaylorModels_norm(f, dom, order)
    else
        R = _enclose_TaylorModels(f, dom, order)
    end
    return R
end

# univariate
function _enclose_TaylorModels(f::Function, dom::Interval, order::Int)
    x0 = Interval(mid(dom))
    x = TaylorModel1(order, x0, dom)
    return evaluate(f(x-x0), dom)
end

# normalized univariate
function _enclose_TaylorModels_norm(f::Function, dom::Interval, order::Int)
    x0 = Interval(mid(dom))
    x = TaylorModel1(order, x0, dom)
    xnorm = normalize_taylor(x.pol, dom - x0, true)
    xnormTM = TaylorModel1(xnorm, 0..0, 0..0, -1..1)
    return evaluate(f(xnormTM), -1..1)
end

# multivariate
function _enclose_TaylorModels(f::Function, dom::IntervalBox{N}, order::Int) where {N}
    x0 = mid(dom)
    set_variables(Float64, "x", order=2order, numvars=N)
    x = [TaylorModelN(i, order, IntervalBox(x0), dom) for i=1:N]
    return evaluate(f(x...), dom - x0)
end

# normalized multivariate
function _enclose_TaylorModels_norm(f::Function, dom::IntervalBox{N}, order::Int) where {N}
    x0 = mid(dom)
    set_variables(Float64, "x", order=2order, numvars=N)

    zBoxN = zeroBox(N)
    sBoxN = symBox(N)
    x = [TaylorModelN(i, order, IntervalBox(x0), dom) for i=1:N]
    xnorm = [normalize_taylor(xi.pol, dom - x0, true) for xi in x]
    xnormTM = [TaylorModelN(xi_norm, 0..0, zBoxN, sBoxN) for xi_norm in xnorm]
    return evaluate(f(xnormTM...), sBoxN)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/BoundingSphere.jl/blob/master/LICENSE

struct ProjectorStack{P <: AbstractVector}
    # matrix that is decomposed into Σ v_i ⊗ v_i* for
    # an orthonormal system v_i
    vs::Vector{P}
end

function Base.push!(p::ProjectorStack, v)
    @assert norm(v) ≈ 1
    push!(p.vs, v)
    p
end
function Base.pop!(p::ProjectorStack)
    pop!(p.vs)
    p
end

function Base.:*(p::ProjectorStack, v::AbstractVector)
    ret = zero(v)
    for vi in p.vs
        ret = ret + vi*dot(vi, v)
    end
    ret
end

"""
    BoundaryDevice

Finds unique spheres determined by prescribed affine independent
boundary points. In the welzl algorithm this problem needs to be solved
in series, where points are pushed and popped from to the boundary.

Subtypes must implement the following interface:

* push_if_stable!(device, pt)::Bool :
* pop!(device): Remove last point from the boundary.
* get_ball(device)::SqBall : Get the last ball from the device.
* length(device)::Int : Get the current count of boundary points
* ismaxlength(device)::Bool: Check if there are dim+1 boundary points in the device
"""
abstract type BoundaryDevice end

Base.isempty(b::BoundaryDevice) = length(b) == 0

"""
    GaertnerBdry

BoundaryDevice that corresponds to M_B in Section 4 of Gaertners paper.

See also: [BoundaryDevice](@ref)
"""
mutable struct GaertnerBdry{P<:AbstractVector,
                            F<:AbstractFloat} <: BoundaryDevice
    centers::Vector{P}
    square_radii::Vector{F}
    # projection onto of affine space spanned by points
    # shifted such that first point becomes origin
    projector::ProjectorStack{P}
    empty_center::P # center of ball spanned by empty boundary
end

function create_boundary_device(pts, alg)
    P = eltype(pts)
    F = eltype(P)
    projector = ProjectorStack(P[])
    centers = P[]
    square_radii = F[]
    empty_center = F(NaN)*first(pts)
    GaertnerBdry(centers, square_radii, projector, empty_center)
end

function Base.length(b::GaertnerBdry)
    @assert length(b.centers) ==
    length(b.square_radii)

    return length(b.centers)
end

function push_if_stable!(b::GaertnerBdry, pt)
    if isempty(b)
        push!(b.square_radii, zero(eltype(pt)))
        push!(b.centers, pt)
        dim = length(pt)
        return true
    end

    q0 = first(b.centers)
    center = b.centers[end]
    C  = center - q0
    r2 = b.square_radii[end]

    Qm = pt - q0
    M = b.projector
    Qm_bar = M*Qm
    residue = Qm - Qm_bar

    e = sqdist(Qm, C) - r2
    z = 2*sqnorm(residue)

    # TODO should we use norm(residue) instead of z here?
    # seems more intuitive, OTOH z is used in the paper
    tol = eps(eltype(pt)) * max(r2, one(r2))
    @assert !isnan(tol)
    isstable = abs(z) > tol
    if isstable

        center_new  = center + (e/z) * residue
        r2new = r2 + (e^2)/(2z)

        push!(b.projector, residue / norm(residue))
        push!(b.centers, center_new)
        push!(b.square_radii, r2new)
    end
    isstable
end

function Base.pop!(b::GaertnerBdry)
    n = length(b)
    pop!(b.centers)
    pop!(b.square_radii)
    if n >= 2
        pop!(b.projector)
    end
    b
end

function get_ball(b::GaertnerBdry)
    if isempty(b)
        c = b.empty_center
        r2 = zero(eltype(c))
    else
        c = b.centers[end]
        r2 = b.square_radii[end]
    end
    SqBall(c,r2)
end

function ismaxlength(b::GaertnerBdry)
    dim = length(b.empty_center)
    length(b) == dim + 1
end


text:
using Images
import VideoIO

# 字符像素
const pixels = raw"$#@&%ZYXWVUTSRQPONMLKJIHGFEDCBA098765432?][}{/\\)(><zyxwvutsrqponmlkjihgfedcba*+1!:`'-,. "

function video2imgs(video_name::String, img_size::Tuple; framerate=8)
    imgs = Array{Gray{Normed{UInt8,8}},2}[]
    f = VideoIO.openvideo(video_name)
    t = 0.0
    step = 1/framerate
    try
        while !eof(f)
            img = imresize(read(f), img_size)  # 读取 image 并 resize
            push!(imgs, Gray.(img))   # 转换成灰度图，然后 push 到 imgs 里边
            seek(f, t += step)  # t 的单位是秒
        end
    catch EOFError  # 读取到末尾
        # do nothing
    end

    imgs
end

function img2ascii(img::Array{Gray{Normed{UInt8,8}},2})
    res = String[]
    for row in 1:height(img)
        line = ""
        for col in 1:width(img)
            percent = img[row, col]
            index = ceil(Int, percent.val * length(pixels)) + 1  # 计算出灰度对应的字符索引

            # 添加字符像素（最后面加一个空格，是因为命令行有行距却没几乎有字符间距，用空格当间距）
            line *= pixels[index > length(pixels) ? length(pixels)-1 : index] * " "
        end
        push!(res, line)
    end
    join(res, "\n")
end

function play_ascii(ascii_pics::Array{String, 1}, framerate=1/8)
    for frame in ascii_pics
        print(frame)  # 输出
        sleep(framerate)
        run(Cmd(`clear`))  # 清屏
    end
end

function play_audio(video_name::String)
    # 异步播放
    @async run(Cmd(`mpv --no-video $video_name`))
end

# 转换成字符动画
videoname = "BadApple.mp4"
imgs = video2imgs(videoname, (48, 64))
ascii_pics = img2ascii.(imgs)

# 播放
play_audio(videoname)
play_ascii(ascii_pics)


text:
fe1ph(v::Vector{Int64})   = hist(v)[2] ./ size(v)[1]

function fe2p(v::Matrix{Int64}) # frequency estimation of one dimensional probability
  @assert (size(v)[2] == 2) "fe2p requires two columns."
  m1 = maximum(v[:,1])
  m2 = maximum(v[:,2])
  r  = counts(v[:,1], v[:,2], (1:m1, 1:m2))
  l  = size(v)[1]
  r  = r ./ l
  # just to get rid of the numerical inaccuracies and make sure its a probability distribution
  s  = sum(r)
  r ./ s
end

# frequency estimation of one dimensional probability
function fe1p(v::Vector{Int64})
  m = maximum(v)
  l  = size(v)[1]
  r  = counts(v, 1:m)
  r = r ./ l
  # just to get rid of the numerical inaccuracies and make sure its a probability distribution
  s = sum(r)
  r ./ s
end

function fe3p(data::Matrix{Int64})
  @assert (size(data)[2] == 3) "fe3p required three columns."
  x = data[:,1]
  y = data[:,2]
  z = data[:,3]

  mx = maximum(x)
  my = maximum(y)
  mz = maximum(z)

  p  = zeros(mx, my, mz)
  for index = 1:length(x)
    p[x[index], y[index], z[index]] = p[x[index], y[index], z[index]] + 1.0
  end
  p ./ float(length(x))
end

# frequency estimation of conditional probability distribution p(x|y,z)
function fe3c2p(data::Matrix{Int64})
  @assert (size(data)[2] == 3) "fe3c2p required three columns."

  pxyz = fe3p(data)
  pyz  = fe2p(data[:,[2:3]])

  for x = 1:size(data)[1]
    for y = 1:size(data)[2]
      for z = 1:size(data)[3]
        pxyz[x,y,z] = pxyz[x,y,z] / pyz[y,z]
      end
    end
  end
  pxyz
end

function fe2c1p(data::Matrix{Int64})
  @assert (size(data)[2] == 3) "fe2c1p required two columns."
  pxy = fe2p(data)
  py  = fe1p(data[:,1])
  @assert (mininum(py) > 0.0) "p(y) must have full support (no zero-valued entries)"
  for x=1:size(pxy)[1]
    for y=1:size(pxy)[2]
      pxy[x,y] = pxy[x,y] / py[y]
    end
  end
  pxy
end


text:
@testset "bishop.jl" begin
    @printf("> bishop\n")

    @test Bb.BISHOP_MASKS[p2u["a1"]] == Bb.BISHOP_MASKS[p2u["h8"]]
    @test Bb.BISHOP_MASKS[p2u["a8"]] == Bb.BISHOP_MASKS[p2u["h1"]]
    @test Bb.diagoSlide(p2u["e3"]) == Bb.BISHOP_MASKS[p2u["e3"]]
    @test Bb.bishopMasksGen()[p2u["d6"]] == reduce(|, [p2u[s] for s = ["b4", "c5", "c7", "e7", "e5", "f4", "g3"]])
end


text:
function RBMs.zerosum!(rbm::ConvRBM)
    RBMs.zerosum_visible!(rbm)
    RBMs.zerosum_hidden!(rbm)
    return rbm
end

RBMs.zerosum_visible!(rbm::ConvRBM{<:AbstractLayer, <:AbstractLayer}) = rbm
RBMs.zerosum_hidden!(rbm::ConvRBM{<:AbstractLayer, <:AbstractLayer})  = rbm

function RBMs.zerosum_visible!(rbm::ConvRBM{<:RBMs.Potts, <:AbstractLayer})
    RBMs.zerosum!(visible(rbm))
    RBMs.zerosum!(weights(rbm); dims = 1)
    return nothing
end

function RBMs.zerosum_hidden!(rbm::ConvRBM{<:AbstractLayer, <:RBMs.Potts})
    RBMs.zerosum!(hidden(rbm))
    RBMs.zerosum!(weights(rbm); dims = ndims(weights(rbm)) - ndims(hidden(rbm)) + 1)
    return nothing
end

function RBMs.zerosum!(rbm::ConvRBM{<:RBMs.Potts, <:RBMs.Potts})
    RBMs.zerosum!(visible(rbm))
    RBMs.zerosum!(hidden(rbm))
    RBMs.zerosum!(weights(rbm); dims = 1)
    RBMs.zerosum!(weights(rbm); dims = ndims(weights(rbm)) - ndims(hidden(rbm)) + 1)
    return nothing
end


text:
module InterpretableValidation
    using Distributions
    using LinearAlgebra
    using SpecialFunctions
    using NLsolve
    using ExprRules
    using ExprOptimization
    using Parameters
    using Random
    using POMDPs
    using POMDPSimulators
    using POMDPPolicies

    # mvrandn
    export mvrandn_μ, mvrandn
    include("mvrandn.jl")

    # Timeseries Distributions
    export TimeseriesDistribution, IID, GaussianProcess, N_pts, squared_exp_kernel
    include("timeseries_distributions.jl")

    # Inverse Logic
    export sample_constraints, flex_not, not_inv,
            parse_implications, parse_implications!, and_expressions,
            all_before, all_before_inv, all_after, all_after_inv, all_between,
            all_between_inv, any_between, any_between_inv,
            and_inv, or_inv, all_inv, any_inv, bitwise_and_inv, bitwise_or_inv
    include("inverse_logic.jl")

    # Constrained Distributions
    export ConstrainedTimeseriesDistribution, isdiscrete, iscontinuous, isfeasible,
           MvTimeseriesDistribution, InfeasibleConstraint,
           greaterthan!, lessthan!, continuous_equality!, discrete_equality!, constrain_timeseries!
    include("constrained_distributions.jl")

    # Optimization and Grammar
    export loss_fn, sample_comparison, create_stl_grammar, create_policy_grammar,
           default_comparison_distribution, default_comparisons, optimize_stl_policy, optimize_timed_stl, set_global_grammar_params,
           discrete_action_mdp, continuous_action_mdp, sample_history, get_policy, policy_loss_fn
    include("optimization.jl")
end