text:
const LONGER_TESTS = false

if haskey(ENV,"GROUP")
    group = ENV["GROUP"]
else
    group = "All"
end

is_APPVEYOR = ( Sys.iswindows() && haskey(ENV,"APPVEYOR") )

macro testset_module(name, expr)
  if name isa String
    mod = gensym()
    testname = name
  elseif name isa Expr && name.head == :(=) && length(name.args) == 2
    mod, testname = name.args[1]
  else
    error("""
          Use `@testset_module` like the following:
          @time @testset_module "Benchmark Tests" begin include("benchmark_tests.jl") end
          @time @testset_module BenchmarkTests = "Benchmark Tests" begin include("benchmark_tests.jl") end
          """)
  end
  quote
    @eval module $mod
      using Test
      @testset $testname $expr
    end
    nothing
  end
end

#Start Test Script

@time begin
if group == "All" || group == "Interface"
  @time @testset_module "Discrete Algorithm Tests" begin include("discrete_algorithm_test.jl") end
  @time @testset_module "Tstops Tests" begin include("ode/ode_tstops_tests.jl") end
  @time @testset_module "Backwards Tests" begin include("ode/ode_backwards_test.jl") end
  @time @testset_module "Initdt Tests" begin include("ode/ode_initdt_tests.jl") end
  @time @testset_module "Mass Matrix Tests" begin include("mass_matrix_tests.jl") end
  @time @testset_module "Differentiation Trait Tests" begin include("differentiation_traits_tests.jl") end
  @time @testset_module "saveat Tests" begin include("ode/ode_saveat_tests.jl") end
  @time @testset_module "save_idxs Tests" begin include("ode/ode_saveidxs_tests.jl") end
  @time @testset_module "Static Array Tests" begin include("static_array_tests.jl") end
  @time @testset_module "Data Array Tests" begin include("data_array_test.jl") end
  @time @testset_module "u_modifed Tests" begin include("umodified_test.jl") end
  @time @testset_module "Composite Algorithm Tests" begin include("composite_algorithm_test.jl") end
  @time @testset_module "Complex Tests" begin include("complex_tests.jl") end
  @time @testset_module "Stiffness Detection Tests" begin include("stiffness_detection_test.jl") end
  @time @testset_module "Export tests" begin include("export_tests.jl") end
  @time @testset_module "Derivative Utilities Tests" begin include("utility_tests.jl") end
end

if group == "All" || group == "Integrators"
  @time @testset_module "Reinit Tests" begin include("reinit_test.jl") end
  @time @testset_module "Events Tests" begin include("ode/ode_event_tests.jl") end
  @time @testset_module "Cache Tests" begin include("ode/ode_cache_tests.jl") end
  @time @testset_module "Iterator Tests" begin include("iterator_tests.jl") end
  @time @testset_module "Integrator Interface Tests" begin include("integrator_interface_tests.jl") end
  @time @testset_module "Add Steps Tests" begin include("ode/ode_add_steps_tests.jl") end
  @time @testset_module "Error Check Tests" begin include("check_error.jl") end
end

if !is_APPVEYOR && ( group == "All" || group == "Regression" )
  @time @testset_module "Linear Tests" begin include("ode/ode_twodimlinear_tests.jl") end
  @time @testset_module "Dense Tests" begin include("ode/ode_dense_tests.jl") end
end

if !is_APPVEYOR && ( group == "All" || group == "AlgConvergence_I" )
  # ~ 250 s
  @time @testset_module "Partitioned Methods Tests" begin include("partitioned_methods_tests.jl") end
  # ~ 400 s
  @time @testset_module "Convergence Tests" begin include("ode/ode_convergence_tests.jl") end
  # ~ 2 s
  @time @testset_module "Adams Variable Coefficients Tests" begin include("ode/adams_tests.jl") end
  # ~ 50 s
  @time @testset_module "Nordsieck Tests" begin include("ode/nordsieck_tests.jl") end
  #@time @testset "Linear Methods Tests" begin include("linear_method_tests.jl") end
  # ~ 170 s
  @time @testset_module "SSPRK Tests" begin include("ode/ode_ssprk_tests.jl") end
  # ~ 25 s
  @time @testset_module "OwrenZen Tests" begin include("owrenzen_tests.jl") end
  @time @testset_module "Runge-Kutta-Chebyshev Tests" begin include("ode/rkc_tests.jl") end
end

if !is_APPVEYOR && ( group == "All" || group == "AlgConvergence_II" )
  # ~ 110 s
  @time @testset_module "Split Methods Tests" begin include("split_methods_tests.jl") end
  # ~ 550 s
  @time @testset_module "Rosenbrock Tests" begin include("ode/ode_rosenbrock_tests.jl") end
  # ~ 40 s
  @time @testset_module "Linear-Nonlinear Methods Tests" begin include("linear_nonlinear_convergence_tests.jl") end
  # ~ 140 s
  @time @testset_module "Linear-Nonlinear Krylov Methods Tests" begin include("linear_nonlinear_krylov_tests.jl") end
end
end # @time


text:
# ---
# title: 736. Parse Lisp Expression
# id: problem736
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String
# link: <https://leetcode.com/problems/parse-lisp-expression/description/>
# hidden: true
# ---
#
# You are given a string `expression` representing a Lisp-like expression to
# return the integer value of.
#
# The syntax for these expressions is given as follows.
#
# * An expression is either an integer, a let-expression, an add-expression, a mult-expression, or an assigned variable. Expressions always evaluate to a single integer.
#
# * (An integer could be positive or negative.)
#
# * A let-expression takes the form `(let v1 e1 v2 e2 ... vn en expr)`, where `let` is always the string `"let"`, then there are 1 or more pairs of alternating variables and expressions, meaning that the first variable `v1` is assigned the value of the expression `e1`, the second variable `v2` is assigned the value of the expression `e2`, and so on **sequentially** ; and then the value of this let-expression is the value of the expression `expr`.
#
# * An add-expression takes the form `(add e1 e2)` where `add` is always the string `"add"`, there are always two expressions `e1, e2`, and this expression evaluates to the addition of the evaluation of `e1` and the evaluation of `e2`.
#
# * A mult-expression takes the form `(mult e1 e2)` where `mult` is always the string `"mult"`, there are always two expressions `e1, e2`, and this expression evaluates to the multiplication of the evaluation of `e1` and the evaluation of `e2`.
#
# * For the purposes of this question, we will use a smaller subset of variable names. A variable starts with a lowercase letter, then zero or more lowercase letters or digits. Additionally for your convenience, the names "add", "let", or "mult" are protected and will never be used as variable names.
#
# * Finally, there is the concept of scope. When an expression of a variable name is evaluated, **within the context of that evaluation** , the innermost scope (in terms of parentheses) is checked first for the value of that variable, and then outer scopes are checked sequentially. It is guaranteed that every expression is legal. Please see the examples for more details on scope.
#
# **Evaluation Examples:**
#
#
#
#     Input: (add 1 2)
#     Output: 3
#
#     Input: (mult 3 (add 2 3))
#     Output: 15
#
#     Input: (let x 2 (mult x 5))
#     Output: 10
#
#     Input: (let x 2 (mult x (let x 3 y 4 (add x y))))
#     Output: 14
#     Explanation: In the expression (add x y), when checking for the value of the variable x,
#     we check from the innermost scope to the outermost in the context of the variable we are trying to evaluate.
#     Since x = 3 is found first, the value of x is 3.
#
#     Input: (let x 3 x 2 x)
#     Output: 2
#     Explanation: Assignment in let statements is processed sequentially.
#
#     Input: (let x 1 y 2 x (add x y) (add x y))
#     Output: 5
#     Explanation: The first (add x y) evaluates as 3, and is assigned to x.
#     The second (add x y) evaluates as 3+2 = 5.
#
#     Input: (let x 2 (add (let x 3 (let x 4 x)) x))
#     Output: 6
#     Explanation: Even though (let x 4 x) has a deeper scope, it is outside the context
#     of the final x in the add-expression.  That final x will equal 2.
#
#     Input: (let a1 3 b2 (add a1 1) b2)
#     **Output** 4
#     Explanation: Variable names can contain digits after the first character.
#
#
#
# **Note:**
#
# * The given string `expression` is well formatted: There are no leading or trailing spaces, there is only a single space separating different components of the string, and no space between adjacent parentheses. The expression is guaranteed to be legal and evaluate to an integer.
# * The length of `expression` is at most 2000. (It is also non-empty, as that would not be a legal expression.)
# * The answer and all intermediate calculations of that answer are guaranteed to fit in a 32-bit integer.
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
"""
    @DBG expression_to_run_only_when_debugging
"""
macro DBG(ex)
    quote
        isdebugging() && $(esc(ex))
        nothing
    end
end

isdebugging() = false
enable_debug() = (@eval isdebugging() = true; nothing)
disable_debug() = (@eval isdebugging() = false; nothing)


text:
"""
A linear combination of Laguerre polynomials.
"""
struct LaguerrePolynomial{T} <: AbstractLaguerre{T}
    # representation of a linear Combination of Laguerre polynomials
    orders::Vector{Int} # order of Laguerre polynomials.
    coeffs::Vector{T} # coefficients of Laguerre polynomials.
    len::Int # length of orders
    function LaguerrePolynomial(orders::Vector{Int}, coeffs::Vector{T}) where T
        length(orders) == length(coeffs) || throw(DimensionMismatch("orders and coeffs are of different length."))
        all(orders .>= 0) || throw(DomainError("negative orders not allowed."))
        # arrange in a ascending order.
        perm = sortperm(orders)
        orders = orders[perm]
        coeffs = coeffs[perm]
        return new{T}(orders, coeffs, length(orders))
    end
end

function LaguerrePolynomial(order::Int, coef::Number)
    # handle a single Laguerre polynomial.
    return LaguerrePolynomial([order], [coef])
end

function LaguerrePolynomial(order::Int)
    # A single Laguerre Polynomial with coefficient 1.
    return LaguerrePolynomial(order, 1)
end

# can add LaguerrePolynomial to a constant, because a constant is just const*L0
function Base.:+(la::LaguerrePolynomial{T}, num::S) where {S <: Number} where T
    R = promote_type(T, S)
    L0 = LaguerrePolynomial(0, num)
    return la + L0
end

# order doesn't matter in addition.
function Base.:+(num::S, la::LaguerrePolynomial{T}) where T where {S <: Number}
    return la + num
end


text:
import GitLab
myauth = GitLab.authenticate(ENV["GITLAB_AUTH"]) # don't hardcode your access tokens!
println("Authentication successful")
options = Dict("private_token" => myauth.token)

# EventListener settings
mysecret = ENV["MY_SECRET"]
myevents = ["Note Hook", "MergeRequest"]
myrepos = [GitLab.repo_by_name("TestProject1"; headers=options)]
myforwards = [HttpCommon.URI("http://myforward1.com"), "http://myforward2.com"] # can be HttpCommon.URIs or URI strings

# Set up Status parameters
pending_params = Dict(
    "state" => "pending",
    "context" => "Benchmarker",
    "description" => "Running benchmarks..."
)

success_params = Dict(
    "state" => "success",
    "context" => "Benchmarker",
    "description" => "Benchmarks complete!"
)

error_params(err) = Dict(
    "state" => "error",
    "context" => "Benchmarker",
    "description" => "Error: $err"
)

# We can use Julia's `do` notation to set up the listener's handler function
listener = GitLab.EventListener(auth = myauth,
                                secret = mysecret,
                                repos = myrepos,
                                events = myevents,
                                forwards = myforwards) do event
    kind, payload, repo = event.kind, event.payload, event.repository

    if kind == "MergeRequest" && payload["action"] == "closed"
        return HttpCommon.Response(200)
    end

    comment = GitLab.Comment(event.payload["object_attributes"])
    if event.payload["object_attributes"]["noteable_type"] == "push"
        sha = event.payload["after"]
    elseif event.payload["object_attributes"]["noteable_type"] == "MergeRequest"
        sha = "$(comment.id)"
    end

    GitLab.create_status(repo, sha; headers = options, params = pending_params)

    try
        # run_and_log_benchmarks isn't actually a defined function, but you get the point
        ## run_and_log_benchmarks(event, "\$(sha)-benchmarks.csv")
        println("Done !")
    catch err
        GitLab.create_status(repo, sha; headers = options, params = error_params(err))
        return HttpCommon.Response(500)
    end

    GitLab.create_status(repo, sha; headers = options, params = success_params)

    return HttpCommon.Response(200)
end

# Start the listener on localhost at port 8000
## GitLab.run(listener, host=IPv4(127,0,0,1), port=8000)
GitLab.run(listener, host=IPv4(0,0,0,0), port=8000)


text:
struct TestComponent <: InfrastructureSystemsComponent
    name::String
    val::Int
    time_series_container::TimeSeriesContainer
    internal::InfrastructureSystemsInternal
end

struct AdditionalTestComponent <: InfrastructureSystemsComponent
    name::String
    val::Int
    time_series_container::TimeSeriesContainer
    internal::InfrastructureSystemsInternal
end

function TestComponent(name, val)
    return TestComponent(name, val, TimeSeriesContainer(), InfrastructureSystemsInternal())
end

function AdditionalTestComponent(name, val)
    return AdditionalTestComponent(
        name,
        val,
        TimeSeriesContainer(),
        InfrastructureSystemsInternal(),
    )
end

get_internal(component::TestComponent) = component.internal
get_internal(component::AdditionalTestComponent) = component.internal
get_val(component::TestComponent) = component.val

function get_time_series_container(component::TestComponent)
    return component.time_series_container
end

function from_json(io::IO, ::Type{TestComponent})
    data = JSON3.read(io, Dict)
    return deserialize(TestComponent, data)
end

function deserialize(::Type{TestComponent}, data::Dict)
    return TestComponent(
        data["name"],
        data["val"],
        deserialize(TimeSeriesContainer, data["time_series"]),
        deserialize(InfrastructureSystemsInternal, data["internal"]),
    )
end

struct TestEvent <: AbstractRecorderEvent
    common::RecorderEventCommon
    val1::String
    val2::Int
    val3::Float64
end

function TestEvent(val1::String, val2::Int, val3::Float64)
    return TestEvent(RecorderEventCommon("TestEvent"), val1, val2, val3)
end

struct TestEvent2 <: AbstractRecorderEvent
    common::RecorderEventCommon
    val::Int
end

function TestEvent2(val::Int)
    return TestEvent2(RecorderEventCommon("TestEvent2"), val)
end

function runtests(args...)
    test_prefix = "test_"
    for arg in args
        if !startswith(arg, test_prefix)
            arg = test_prefix * arg
        end
        push!(ARGS, arg)
    end

    try
        include("test/runtests.jl")
    finally
        empty!(ARGS)
    end
end


text:
using Arrhenius
using PyCall

```
1. config Conda.jl
2. add PyCall
3. build PyCall
```
ct = pyimport("cantera")

gas = ct.Solution("gri30")

gas.X


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type reboiler
	reboiler()=begin
		PP=outers.PP
		NComp=outers.NComp
		[
			:(LI = Levelpercent_Initial),
			:(OutletLiquid.T = Initial_Temperature),
			:(OutletLiquid.z(1:NComp-1) = Initial_Composition(1:NComp-1)/sum(Initial_Composition)),
		],
		[
			"Initial level Percent","Initial Temperature","Initial Composition",
		],
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Initial liquid height in Percent",
				:Default=>0.70
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Initial Temperature of Reboiler"
			)),
			fill(positive(Dict{Symbol,Any}(
				:Brief=>"Initial Liquid Composition",
				:Lower=>1E-6
			)),(NComp)),
			VesselVolume(Dict{Symbol,Any}(
				:Brief=>"Vessel Geometry",
				:Symbol=>" "
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Liquid inlet stream",
				:PosX=>0.17,
				:PosY=>1,
				:Symbol=>"_{in}^{Liquid}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream",
				:PosX=>0.53,
				:PosY=>1,
				:Symbol=>"_{out}^{Liquid}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream",
				:PosX=>0.17,
				:PosY=>0,
				:Symbol=>"_{out}^{Vapour}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Heat supplied",
				:Protected=>true,
				:PosX=>1,
				:PosY=>0.08,
				:Symbol=>"Q_{in}"
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Temperature  Indicator of Reboiler",
				:Protected=>true,
				:PosX=>0.44,
				:PosY=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Level Indicator of Reboiler",
				:Protected=>true,
				:PosX=>0.53,
				:PosY=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Pressure Indicator of Reboiler",
				:Protected=>true,
				:PosX=>0.35,
				:PosY=>0
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar Holdup in the tray",
				:Protected=>true
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar liquid holdup",
				:Protected=>true
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar vapour holdup",
				:Protected=>true
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy Holdup on tray",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar volume",
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Vapour Density",
				:Protected=>true,
				:Symbol=>"\\rho"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P",
				:Protected=>true
			)),
			[
				:(diff(M)= InletLiquid.F*InletLiquid.z - OutletLiquid.F*OutletLiquid.z - OutletVapour.F*OutletVapour.z),
				:(diff(E) = InletLiquid.F*InletLiquid.h - OutletLiquid.F*OutletLiquid.h - OutletVapour.F*OutletVapour.h + InletQ),
				:(M = ML*OutletLiquid.z + MV*OutletVapour.z),
				:(E = ML*OutletLiquid.h + MV*OutletVapour.h - OutletLiquid.P*Geometry.Vtotal),
				:(sum(OutletLiquid.z)=1.0),
				:(sum(OutletLiquid.z)=sum(OutletVapour.z)),
				:(rhoV = PP.VapourDensity(OutletVapour.T, OutletVapour.P, OutletVapour.z)),
				:(vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(vV = PP.VapourVolume(OutletVapour.T, OutletVapour.P, OutletVapour.z)),
				:(PP.LiquidFugacityCoefficient(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)*OutletLiquid.z = PP.VapourFugacityCoefficient(OutletVapour.T, OutletVapour.P, OutletVapour.z)*OutletVapour.z),
				:(OutletLiquid.P = OutletVapour.P),
				:(OutletLiquid.T = OutletVapour.T),
				:(OutletLiquid.P = InletLiquid.P - Pdrop),
				:(Geometry.Vtotal = ML*vL + MV*vV),
				:(ML * vL = Geometry.Vfilled),
				:(TI * "K" = OutletLiquid.T),
				:(PI * "atm" = OutletLiquid.P),
				:(LI*Geometry.Vtotal= Geometry.Vfilled),
			],
			[
				"Component Molar Balance","Energy Balance","Molar Holdup","Energy Holdup","Mol Fraction Normalisation","Mol fraction Constraint","Vapour Density","Liquid Volume","Vapour Volume","Chemical Equilibrium","Mechanical Equilibrium","Thermal Equilibrium","Pressure Drop","Geometry Constraint","Liquid Level","Temperature Indicator","Pressure Indicator","Level indicator",
			],
			[:PP,:NComp,:Levelpercent_Initial,:Initial_Temperature,:Initial_Composition,],
			[:Geometry,:InletLiquid,:OutletLiquid,:OutletVapour,:InletQ,:TI,:LI,:PI,:M,:ML,:MV,:E,:vL,:vV,:rhoV,:Pdrop,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Levelpercent_Initial::positive
	Initial_Temperature::temperature
	Initial_Composition::Array{positive}
	Geometry::VesselVolume
	InletLiquid::stream
	OutletLiquid::liquid_stream
	OutletVapour::vapour_stream
	InletQ::power
	TI::control_signal
	LI::control_signal
	PI::control_signal
	M::Array{mol}
	ML::mol
	MV::mol
	E::energy
	vL::volume_mol
	vV::volume_mol
	rhoV::dens_mass
	Pdrop::press_delta
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export reboiler
function initial(in::reboiler)
	addEquation(1)
	addEquation(2)
	addEquation(3)
end
function setEquationFlow(in::reboiler)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
end
function atributes(in::reboiler,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Reboiler"
	fields[:Brief]="Model of a dynamic reboiler - kettle with control."
	fields[:Info]="== ASSUMPTIONS ==
* perfect mixing of both phases;
* thermodynamics equilibrium;
* no liquid entrainment in the vapour stream.

== SET ==
*Orientation: vessel position - vertical or horizontal;
*Heads (bottom and top heads are identical)
**elliptical: 2:1 elliptical heads (25% of vessel diameter);
**hemispherical: hemispherical heads (50% of vessel diameter);
*Diameter: Vessel diameter;
*Lenght: Side length of the cylinder shell;

== SPECIFY ==
* the InletLiquid stream;
* the outlet flows: OutletVapour.F and OutletLiquid.F;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).

== OPTIONAL ==
* the reboiler model has three control ports
** TI OutletLiquid Temperature Indicator;
** PI OutletLiquid Pressure Indicator;
** LI Level Indicator of Reboiler;

== INITIAL CONDITIONS ==
* Initial_Temperature :  the reboiler temperature (OutletLiquid.T);
* Levelpercent_Initial : the reboiler liquid level in percent (LI);
* Initial_Composition : (NoComps) OutletLiquid compositions.
"
	drive!(fields,_)
	return fields
end
reboiler(_::Dict{Symbol,Any})=begin
	newModel=reboiler()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
"""
    rse = PartitionTreeMatrix_unbalanced_haar(matrix::Matrix{Float64},jmax::Int;method::Symbol = :Totalvar, p::Float64 = 2.)

Perform the quadtree partition of image matrix depending on total-variation or changes on the boundary.

### Input Arguments
* `matrix`: The image matrix.
* `jmax`: The maximum level of partition.
* `method`: ':Totalvar' indicates using total variation, the other indicates using changes on the boundary.
* `p`: The power used in the criterion to balance the sizes of the sub-regions. The larger `p`, the more weight on balancing the sub-regions.

### Output Arguments
* `rse`: Region start and end. Each pair represents the start and end of a region in the matrix, corresponding to one node in the quadtree.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""

function PartitionTreeMatrix_unbalanced_haar(matrix::Matrix{Float64},jmax::Int;method::Symbol = :Totalvar, p::Float64 = 2.)
    rows,cols = size(matrix)
    N = rows*cols
    rse = Vector{Vector{Int}}(jmax)
    rse[1] = [1,N];

    j = 1
    while j<=jmax-1
        rr = 1
        rse[j+1] = zeros(4*length(rse[j]))
        for r = 1:2:(length(rse[j]) - 1)
            rs_rows,rs_cols = ind2sub((rows,cols),rse[j][r])
            re_rows,re_cols = ind2sub((rows,cols),rse[j][r+1])

            n_rows = re_rows-rs_rows+1
            n_cols = re_cols-rs_cols+1

            n1_rows,n1_cols = PartitionTreeMatrix_unbalanced_haar_subblock(matrix[rs_rows:re_rows,rs_cols:re_cols],method = method, p = p)

            if !(n1_rows==0 || n1_rows == n_rows) && !(n1_cols==0 ||n1_cols==n_cols)
                child_rows = [rs_rows, rs_rows+n1_rows-1, rs_rows+n1_rows, re_rows, rs_rows, rs_rows+n1_rows-1, rs_rows+n1_rows, re_rows];
                child_cols = [rs_cols, rs_cols+n1_cols-1, rs_cols, rs_cols+n1_cols-1, rs_cols+n1_cols, re_cols, rs_cols+n1_cols, re_cols];
                for k = 1:8
                    rse[j+1][rr+k-1] = sub2ind((rows,cols),child_rows[k],child_cols[k])
                end
                rr = rr + 8

            elseif (n1_rows == 0 || n1_rows == n_rows) && !(n1_cols == 0 || n1_cols == n_cols)
                child_rows = [rs_rows, re_rows, rs_rows, re_rows];
                child_cols = [rs_cols, rs_cols+n1_cols-1, rs_cols+n1_cols, re_cols];
                for k = 1:4
                    rse[j+1][rr+k-1] = sub2ind((rows,cols),child_rows[k],child_cols[k])
                end
                rr = rr + 4

            elseif (n1_cols == 0 || n1_cols == n_cols) && !(n1_rows == 0 || n1_rows == n_rows)
                child_rows = [rs_rows, rs_rows+n1_rows-1, rs_rows+n1_rows, re_rows];
                child_cols = [rs_cols, re_cols, rs_cols, re_cols];
                for k = 1:4
                    rse[j+1][rr+k-1] = sub2ind((rows,cols),child_rows[k],child_cols[k])
                end
                rr = rr + 4
            else
                rse[j+1][rr:rr+1] = rse[j][r:r+1]
                rr = rr + 2
            end
        end
        j = j+1
        rse[j] = rse[j][rse[j].!=0]
    end
    return rse
end

"""
     row_cut,col_cut = PartitionTreeMatrix_unbalanced_haar_subblock(matrix::Matrix{Float64};method::Symbol = :Totalvar, p::Float64 = 2.)

Partition a matrix block into four sub-blocks.

### Input Arguments
* `matrix`: The image matrix.
* `method`: ':Totalvar' indicates using total variation, the other indicates using changes on the boundary.
* `p`: The power used in the criterion to balance the sizes of the sub-regions. The larger `p`, the more weight on balancing the sub-regions.

### Output Arguments
* `row_cut`: The cut index on the row direction.
* `col_cut`: The cut index on the col direction.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""

function PartitionTreeMatrix_unbalanced_haar_subblock(matrix::Matrix{Float64};method::Symbol = :Totalvar, p::Float64 = 2.)
    m,n = size(matrix)
    Entropy = Inf

    if m>1 && n>1
        for i=1:m-1
            for j=1:n-1
                if method == :Totalvar
                    Entropy_temp = total_var(matrix[1:i,1:j])/(i*j)^p + total_var(matrix[i+1:m,1:j])/((m-i)*j)^p + total_var(matrix[1:i,j+1:n])/(i*(n-j))^p + total_var(matrix[i+1:m,j+1:n])/((m-i)*(n-j))^p;
                else
                    Entropy_temp = -( norm(matrix[i,:]-matrix[i+1,:],1) + norm(matrix[:,j] - matrix[:,j+1],1) )/((i*j)^(-p) + ((m-i)*j)^(-p) + (i*(n-j))^(-p) + ((m-i)*(n-j))^(-p));
                end
                if Entropy_temp < Entropy
                    Entropy = Entropy_temp
                    row_cut = i
                    col_cut = j
                end
            end
        end
    elseif m == 1 && n>1
        row_cut = 0
        for j = 1:n-1
            if method == :Totalvar
                Entropy_temp = total_var(matrix[1,1:j])/j^p + total_var(matrix[1,j+1:n])/(n-j)^p;
            else
                Entropy_temp = -norm(matrix[:,j] - matrix[:,j+1],1)*(j^(-p) + (n-j)^(-p));
            end
            if Entropy_temp < Entropy
                Entropy = Entropy_temp
                col_cut = j
            end
        end
    elseif n == 1 && m>1
        col_cut = 0
        for i = 1:m-1
            if method == :Totalvar
                Entropy_temp = total_var(matrix[1:i,1])/i^p + total_var(matrix[i+1:m,1])/(m-i)^p;
            else
                Entropy_temp = -norm(matrix[i,:]-matrix[i+1,:],1)*(i^(-p) + (m-i)^(-p));
            end
            if Entropy_temp < Entropy
                Entropy = Entropy_temp
                row_cut = i
            end
        end
    else
        row_cut = 0
        col_cut = 0
    end
    return row_cut,col_cut
end

"""
     v = total_var(matrix::Matrix{Float64})

Compute the total variation of matrix

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function total_var(matrix::Matrix{Float64})
    m,n = size(matrix)
    row_val = abs.(matrix[1:m-1,:] - matrix[2:m,:])
    col_val = abs.(matrix[:,1:n-1] - matrix[:,2:n])
    v = norm(vcat(row_val[:],col_val[:]),1)
    return v
end

function PartitionTreeMatrix_unbalanced_haar_binarytree(matrix::Matrix{Float64}, jmax::Int64, p::Float64)
    N = size(matrix,2)
    rs = fill(0, N+1, jmax)
    rs[1,:] .= 1
    rs[2,1] = N+1

    for j = 1:jmax - 1
        rr = 1
        regioncount = sum(!iszero, rs[:,j]) - 1

        for r = 1:regioncount
            rs1 = rs[r,j]
            rs2 = rs[r+1,j]
            n = rs2 - rs1
            if n > 1
                n1 = PartiitonTreeMatrix_unbalanced_haar_col_partition(matrix[:,rs1:rs2-1],p)
                rs[rr+1,j+1] = rs1+n1
                rs[rr+2,j+1] = rs2
                rr = rr+2

            elseif n == 1
                rs[rr+1,j+1] = rs2
                rr = rr+1
            end
        end
    end

    GP = GraphPart(Vector{Int}(1:N), rs)
    return GP
end

function PartiitonTreeMatrix_unbalanced_haar_col_partition(matrix::Matrix{Float64},p::Float64)
    n = size(matrix,2)
    entropy = Inf
    col_cut = 1
    for j = 1:n-1
        entropy_temp = total_var(matrix[:, 1:j])/(j^p) + total_var(matrix[:,j+1:n])/(n-j)^p
        if entropy_temp < entropy
            entropy = entropy_temp
            col_cut = j
        end
    end
    return col_cut
end


text:
using LazyJSON
using Test

const JSON = LazyJSON

@testset "AbstractString" begin

function validmap(s)
    s = String([s...])
    return map(i->JSON.string_index_isvalid(s, i) ? 1 : 0, 1:sizeof(s))
end

s =                  ['\\' '"']
@test validmap(s) == [  1,  0 ]

s =                  [' ' '\\' '"']
@test validmap(s) == [ 1,   1,  0]

s =                  ['\\' '\\' '\\' '\\' '\\' '\\' '\\']
@test validmap(s) == [  1,   0,   1,   0,   1,   0,   1 ]

s =                  ['\\' '\\' 'n' '\\' '\\' '\\' '\\' '\\']
@test validmap(s) == [  1,   0,  1,   1,   0,   1,   0,   1 ]

s =                  ['\\' '\\' '\\' 'n' '\\' '\\' '\\' '\\' '\\']
@test validmap(s) == [  1,   0,   1,  0,   1,   0,   1,   0,   1 ]

s =                  [' ' '\\' '\\' '\\' '\\' '\\' '\\' '\\']
@test validmap(s) == [ 1,   1,   0,   1,   0,   1,   0,   1 ]

s =                  [' ' '\\' '\\' 'n' '\\' '\\' '\\' '\\' '\\']
@test validmap(s) == [ 1,   1,   0,  1,   1,   0,   1,   0,   1 ]

s =                  [' ' '\\' '\\' '\\' 'n' '\\' '\\' '\\' '\\' '\\']
@test validmap(s) == [ 1,   1,   0,   1,  0,   1,   0,   1,   0,   1 ]

s =                  ['\\' '\\' '\\' '/' '\\' 'n' '\\' 'r' '\\' 't']
@test validmap(s) == [ 1,   0,   1,   0,  1,   0,   1,  0,   1,  0 ]

s =                  ['\\' 'b' '\\' 'f' '\\' 'u']
@test validmap(s) == [ 1,   0,   1,  0,   1,  0 ]

s =                  ['\\' 'a' '\\' 'c' '\\' 'd' '\\' 'e' '\\' 'g' '\\' 'h']
@test validmap(s) == [ 1,   1,   1,  1,   1,  1,   1,  1,   1,  1,   1,  1 ]

s =                  [' ' '\\' '\\' '\\' '/' '\\' 'n' '\\' 'r' '\\' 't']
@test validmap(s) == [ 1,   1,   0,   1,   0,  1,   0,   1,  0,   1,  0 ]

s =                  [' ' '\\' 'b' '\\' 'f' '\\' 'u']
@test validmap(s) == [ 1,   1,   0,   1,  0,   1,  0 ]

s =                  [' ' '\\' 'a' '\\' 'c' '\\' 'd' '\\' 'e' '\\' 'g' '\\' 'h']
@test validmap(s) == [ 1,   1,   1,   1,  1,   1,  1,   1,  1,   1,  1,   1,  1 ]

s =                  ['\\' 'u' '1' '2' '3' '4' '!']
@test validmap(s) == [ 1,   0,  0,  0,  0,  0,  1 ]

s =                  ['\\' 'u' '1' '2' '3' '!']
@test validmap(s) == [ 1,   0,  0,  0,  0,  1 ] # FIXME invalid hex ?

s =                  ['\\' 'u' '!' '2' '3' '4']
@test validmap(s) == [ 1,   0,  1,  1,  1,  1 ] # FIXME invalid hex ?

s =                  ['\\' 'u' '1' '2' '3']
@test validmap(s) == [ 1,   0,  0,  0,  0,] # FIXME truncated ?

s =                  [' ' '\\' 'u' 'a' 'b' 'c' 'd' '!']
@test validmap(s) == [ 1,   1,  0,  0,  0,  0,  0,  1 ]

s =                  ['\\' 'u' 'f' 'g' 'h' 'i' '!']
@test validmap(s) == [  1,  0,  0,  1,  1,  1,  1 ] # FIXME invalid hex ?

s =                  ['\\' 'u' 'd' '8' '0' '0' '\\' 'u' 'd' '8' '0' '0']
@test validmap(s) == [  1,  0,  0,  0,  0,  0,   1,  0,  0,  0,  0,  0 ]

s =                  ['\\' 'u' 'd' '8' '0' '0' '\\' 'u' '1' '2' '3' '4']
@test validmap(s) == [  1,  0,  0,  0,  0,  0,   1,  0,  0,  0,  0,  0 ]

s =                  ['\\' 'u' 'd' '8' '0' '1' '\\' 'u' 'd' 'c' '3' '7'] # pair
@test validmap(s) == [  1,  0,  0,  0,  0,  0,   0,  0,  0,  0,  0,  0 ]

s =                  ['\\' 'u' 'd' '8' '0' '0' '\\' '\\' '\\' '\\' 'u' '1' ]
@test validmap(s) == [  1,  0,  0,  0,  0,  0,   1,   0,   1,   0,  1,  1  ]

s =                  ['\\' 'u' 'd' '8' '0' '0' '\\' '\\' '\\' 'u' '1' ]
@test validmap(s) == [  1,  0,  0,  0,  0,  0,   1,   0,   1,  0,  0  ]

s =                  [' ' '\\' '"' ' ' '\u1234' '!' ]
@test validmap(s) == [ 1,   1,  0,  1,  1, 0, 0, 1  ]

j = "\"Hello!\""

js = JSON.String(j, 1)
@test ncodeunits(js) == 6
@test js == "Hello!"

c, i = iterate(js)
@test c == 'H'
@test i == 2
c, i = iterate(js, i)
@test c == 'e'

@test isvalid(js, 1)
@test codeunit(js, 1) == UInt8('H')
@test js[1] == 'H'
@test isvalid(js, 2)
@test codeunit(js, 2) == UInt8('e')
@test js[2] == 'e'

@test [c for c in js] == ['H', 'e', 'l', 'l', 'o', '!']
@test '"' * String([c for c in js]) * '"' == j

js = JSON.String("""{"foo": "Foo\\u1234Bar"}""", 9)
@test convert(SubString, js) == "FooሴBar"
@test ncodeunits(js) == 12
@test codeunit(js, 1) == UInt8('F')
@test js[1] == 'F'
@test isvalid(js, 2)
@test codeunit(js, 2) == UInt8('o')
@test js[2] == 'o'
@test codeunit(js, 3) == UInt8('o')
@test js[3] == 'o'
@test codeunit(js, 4) == UInt16('\\')
@test js[4] == '\u1234'
@test !isvalid(js, 5)

@test String([c for c in js]) == "Foo\u1234Bar"
@test js == "Foo\u1234Bar"
@test String([c for c in js]) == "FooሴBar"
@test js == "FooሴBar"

js = JSON.String("\"\\ud801\\udc37\"", 1)
@test String([c for c in js]) == "𐐷"
@test js == "𐐷"

js = JSON.String("\"\\ud83e\\udd16\"", 1)
@test String([c for c in js]) == "🤖"

js = JSON.String("\"🤖\"", 1)
@test String([c for c in js]) == "🤖"

js = JSON.String("\"$(Char(129302))\"", 1)
@test String([c for c in js]) == "🤖"
@test js == "🤖"

end # testset


text:
using FactCheck

import NKLandscapes
const NK = NKLandscapes

facts("NKLandscapes") do
  include("genotype.jl")
  include("landscape.jl")
  include("metapopulation.jl")
  include("migration.jl")
  include("mutation.jl")
  include("neighbors.jl")
  include("population.jl")
  include("walks.jl")

  include("fastfunc.jl")
  include("fastnets.jl")
end


text:
################################################################################
# Custom combined place-and-route
################################################################################
function lowtemp_pnr(m::Map, iterations; place_kwargs...)
    # Do a more low-level treatment of the place and route routine - acquire
    # a SAStruct and use it directly. If routing fails, do a low-temp placement
    # to perturb the location of processors and try routing again.
    #
    # Repeat for the specified number of iterations.
    pstruct = SAStruct(m)
    # First placement - capture the final state.
    state = SA.place(pstruct; place_kwargs...)
    # Record the placement back to the Map.
    Mapper2.SA.record(m, pstruct)

    # Try routing
    m = route(m)
    # If it passes a routing check, return.
    check_routing(m, true) && return m

    @info """
        Initial attempt at routing failed. Trying again.
        """

    for i in 1:iterations
        # Low temperature perturbation.
        SA.place(pstruct;
                 place_kwargs...,
                 supplied_state = state,
                )
        SA.record(m, pstruct)
        # Try routing again.
        m = route(m);
        check_routing(m, true) && return m
    end
    @error "Mapping Failed!"

    return m
end

function basic_place(arch::TopLevel, taskgraph::Taskgraph)
    m = NewMap(arch, taskgraph)
    place(m)
    MapperCore.save(m, "temp.jls")
end

function basic_route(arch::TopLevel, taskgraph::Taskgraph)
    m = NewMap(arch, taskgraph)
    MapperCore.load(m, "temp.jls")
    route(m)
    return m.mapping
end

################################################################################
# Custom Placement Functions
################################################################################
function shotgun_placement(arch     ::TopLevel,
                           taskgraph::Taskgraph;
                           nplacements::Int = 1,
                           nsamples::Int    = 1,
                           kwargs...)

    @info "Total Placements: $(nplacements)"
    # closure for parallel placement
    p(x) = low_temp_placement(arch, taskgraph, x;
                              nsamples = nsamples,
                              kwargs...)

    # Parallelize placements
    pmap(x -> p(x), 1:nplacements)
end

function low_temp_placement(arch            ::TopLevel,
                            taskgraph       ::Taskgraph,
                            iteration_number::Int;
                            nsamples = 1,
                            place_kwargs...)

    @info "Sublacement iteration $iteration_number"
    # Construct a new Map object
    m = NewMap(arch, taskgraph)
    # Get the placement structure
    pstruct = SAStruct(m)
    # Scope "state" out of the loop to avoid renaming.
    local state
    for i = 1:nsamples
        if i == 1
            state = SA.place(pstruct; place_kwargs...)
        else
            SA.place(pstruct;
                    place_kwargs...,
                    supplied_state = state,
                    warmer = SA.TrueSAWarm()
                )
        end
        Mapper2.SA.record(m, pstruct)

        # serialize to temp/ dir so routing can find it.
        savename    = "$(iteration_number)_$(i)"
        save_path   = joinpath(PKGDIR, "temp", savename)
        Mapper2.MapperCore.save(m, save_path)
    end
end

################################################################################
# Custom Routing Functions
################################################################################
function low_temp_route(arch, taskgraph;
                        nplacements::Int    = 1,
                        nsamples::Int       = 1,
                        route_kwargs...)

    # create closure for parallelism
    r(x) = multisample_routing(arch, taskgraph, x, nsamples)

    results = pmap(x -> r(x), 1:nplacements)
    return results
end

function multisample_routing(arch, taskgraph, iteration_number, nsamples)
    results = Mapping[]

    m = NewMap(arch, taskgraph)
    for i in 1:nsamples
        # load serialized placement mapping
        load_path = joinpath(PKGDIR, "temp", "$(iteration_number)_$(i)")
        Mapper2.MapperCore.load(m, load_path)
        try
            m = route(m)
        catch
        end
        push!(results, m.mapping)
    end
    return results
end


text:
"""
    to_cartesian(theta, phi)

Calculates the cartesian vector (x, y, z) for a given pair of angles (polar/z
angle theta, azimut/xy angle phi).
"""
function to_cartesian(theta, phi)
    Vec3f0(
        sin(theta) * cos(phi),
        sin(theta) * sin(phi),
        cos(theta)
    )
end

"""
    face_mesh(B::SSHBinner, extrude=0)

Generates a tiled/disconnected mesh where each face represents a bin.

`extrude` allows faces to be moved outwards (≈ (1 + extrude)eᵣ). For
`extrude > 0` the faces will sperate, making their shape easier to see.
"""
function face_mesh(B::SSHBinner, extrude=0)
    vertices = Point3f0[]
    faces = GLTriangleFace[]

    push!(vertices, Point3f0(0, 0, 1) + extrude * Point3f0(0,0,1))
    theta2 = B.thetas[2]
    N = B.phi_divisions[2]
    for k in 0:N-1
        phi1 = 2pi * k/N
        phi2 = 2pi * (k+1) / N
        pos_1 = to_cartesian(theta2, phi1) + extrude * Point3f0(0,0,1)
        pos_2 = to_cartesian(theta2, phi2) + extrude * Point3f0(0,0,1)
        push!(vertices, pos_1)
        push!(vertices, pos_2)
        push!(faces, GLTriangleFace(1, 2+2k, 3+2k))
    end

    for i in 2:length(B.phi_divisions)-1
        theta1 = B.thetas[i]
        theta2 = B.thetas[i+1]
        N = B.phi_divisions[i]

        for k in 0:N-1
            phi1 = 2pi * k / N
            phi2 = 2pi * (k+1) / N
            pos_11 = to_cartesian(theta1, phi1)
            pos_12 = to_cartesian(theta1, phi2)
            pos_21 = to_cartesian(theta2, phi1)
            pos_22 = to_cartesian(theta2, phi2)

            # extrude "rects"
            dir = 0.25(pos_11 + pos_12 + pos_21 + pos_22)
            pos_11 += extrude*dir
            pos_12 += extrude*dir
            pos_21 += extrude*dir
            pos_22 += extrude*dir

            l = length(vertices)
            append!(vertices, [pos_11, pos_12, pos_21, pos_22])
            push!(faces, GLTriangleFace(l+1, l+3, l+2))
            push!(faces, GLTriangleFace(l+2, l+3, l+4))
        end
    end

    push!(vertices, Point3f0(0, 0, -1) + extrude * Point3f0(0,0,-1))
    l = length(vertices)
    theta2 = B.thetas[end-1]
    N = B.phi_divisions[end-1]

    for k in 0:N-1
        phi1 = 2pi * k / N
        phi2 = 2pi * (k+1) / N
        pos_1 = to_cartesian(theta2, phi1) + extrude * Point3f0(0,0,-1)
        pos_2 = to_cartesian(theta2, phi2) + extrude * Point3f0(0,0,-1)
        push!(vertices, pos_1)
        push!(vertices, pos_2)
        push!(faces, GLTriangleFace(l, l+2+2k, l+1+2k))
    end

    # GeometryBasics.Mesh(GeometryBasics.meta(vertices; normals=ns), faces)
    normal_mesh(vertices, faces)
end

"""
    vertex_mesh(B::SSHBinner)

Generates a mesh where each vertex represents (the center of) a bin.

This type of mesh makes it easy to visualize the histogram with colors. Mapping
each bin-filling to a color (e.g. using `to_hue(B)`) gives colors
corresponding to the vertices of this mesh (in order). This can be plotted with
`Makie` using `Makie.mesh(vertexl_mesh(B), color = colors)`.

See also: [`bin_positions`](@ref), [`face_mesh`](@ref)
"""
function vertex_mesh(B::SSHBinner)
    faces = GLTriangleFace[]

    N = 0
    for i in 2:length(B.phi_divisions)
        top_idx = 1
        bot_idx = 1
        top = top_step = 1 / B.phi_divisions[i-1]
        bot = bot_step = 1 / B.phi_divisions[i]
        N2 = N + B.phi_divisions[i-1]

        while (bot < 0.99999) || (top < 0.99999)
            if top < bot
                top += top_step
                top_idx += 1
                i1 = N + mod1(top_idx, B.phi_divisions[i-1])
                i2 = N + top_idx - 1
                i3 = N2 + bot_idx
                push!(faces, GLTriangleFace(i1, i2, i3))
            else
                bot += bot_step
                bot_idx += 1
                i1 = N + top_idx
                i2 = N2 + bot_idx - 1
                i3 = N2 + mod1(bot_idx, B.phi_divisions[i])
                push!(faces, GLTriangleFace(i1, i2, i3))
            end
        end

        if i < length(B.phi_divisions)
            bot_idx += 1
            i1 = N + top_idx
            i2 = N2 + bot_idx - 1
            i3 = N2 + mod1(bot_idx, B.phi_divisions[i])
            push!(faces, GLTriangleFace(i1, i2, i3))
        end

        if i-1 > 1
            top_idx += 1
            i1 = N + mod1(top_idx, B.phi_divisions[i-1])
            i2 = N + top_idx - 1
            i3 = N2 + mod1(bot_idx, B.phi_divisions[i])
            top_idx = mod1(top_idx + 1, B.phi_divisions[i-1])
            push!(faces, GLTriangleFace(i1, i2, i3))
        end
        N = N2
    end

    # GeometryBasics.Mesh(GeometryBasics.meta(vertices; normals=ns), faces)
    normal_mesh(bin_positions(B), faces)
end

"""
    bin_positions(B::SSHBinner)

Returns a list of points corresponding to the center positions of each bin.
This list is equivalent to the vertices used in `vertex_mesh(B)`.

See also: [`vertex_mesh`](@ref)
"""
function bin_positions(B::SSHBinner)
    points = Point3f0[Point3f0(0, 0, 1)]

    for i in 2:length(B.thetas)-2
        theta = 0.5(B.thetas[i] + B.thetas[i+1])
        N = B.phi_divisions[i]

        for k in 0:N-1
            phi = 2pi * (k + 0.5) / N
            push!(points, to_cartesian(theta, phi))
        end
    end
    push!(points, Point3f0(0, 0, -1))

    points
end

"""
    line_segments_minimal(B::SSHBinner)

Generates linesegments marking each bin area. This version is minimal in the
sense that doesn't interpolate any points.

See also: [`line_segments`](@ref)
"""
function line_segments_minimal(B::SSHBinner)
    lines = Point3f0[]

    # first ring
    theta = B.thetas[2]
    N = B.phi_divisions[2]
    for k in 0:N-1
        phi1 = 2pi * k / N
        phi2 = 2pi * (k+1) / N
        append!(lines, [to_cartesian(theta, phi1), to_cartesian(theta, phi2)])
    end

    # Any thing after that is 1 vertical and 1 horizontal line
    for i in 3:length(B.thetas)-1
        theta1 = B.thetas[i-1]
        theta2 = B.thetas[i]
        if B.phi_divisions[i-1] == B.phi_divisions[i]
            N = B.phi_divisions[i]
        elseif B.phi_divisions[i-1] < B.phi_divisions[i]
            N = B.phi_divisions[i-1]
        else
            N = B.phi_divisions[i-1]
        end

        for k in 0:N-1
            phi1 = 2pi * k / N
            phi2 = 2pi * (k+1) / N

            # horizontal lines
            append!(lines, [to_cartesian(theta2, phi1), to_cartesian(theta2, phi2)])
            # vertical line
            append!(lines, [to_cartesian(theta2, phi1), to_cartesian(theta1, phi1)])
        end
    end

    lines
end

"""
    line_segments(B::SSHBinner; N_fragments = 32)

Generates line segments marking each bin area. Extra points may be interpolated
to generate smooth arcs.

Specifically, points are added to an arc, such that it would include at least
`N_fragments` points over an angle 0:pi.

See also: [`line_segments_minimal`](@ref)
"""
function line_segments(B::SSHBinner; N_fragments=32)
    lines = Point3f0[]

    # horizontal rings
    for theta in B.thetas[2:end-1]
        z = cos(theta)
        # Number of steps
        R_xy = sqrt(1.0 - z^2)
        # circumference = 2πR, make more points if longer line
        for k in 0:2N_fragments-1
            phi1 = pi * k / N_fragments
            phi2 = pi * (k+1) / N_fragments
            append!(lines, [to_cartesian(theta, phi1), to_cartesian(theta, phi2)])
        end
    end

    # TODO recast this from phi(theta) to theta(phi)
    substeps = div(N_fragments, length(B.thetas))+1
    # vertical lines
    for i in 3:length(B.thetas)-1
        theta1 = B.thetas[i-1]
        theta2 = B.thetas[i]
        if B.phi_divisions[i-1] == B.phi_divisions[i]
            N = B.phi_divisions[i]
        elseif B.phi_divisions[i-1] < B.phi_divisions[i]
            N = B.phi_divisions[i-1]
        else
            N = B.phi_divisions[i-1]
        end

        for k in 0:N-1
            phi1 = 2pi * k / N
            phi2 = 2pi * (k+1) / N
            for c in 0:substeps-1
                t1 = theta1 + c/substeps * (theta2 - theta1)
                t2 = theta1 + (c+1)/substeps * (theta2 - theta1)

                # vertical line
                append!(lines, [
                    to_cartesian(t2, phi1),
                    to_cartesian(t1, phi1)
                ])
            end
        end
    end

    lines
end


text:
include("iPredict.jl/src/iPredict.jl")
using .iPredict
using DataFrames
using Plots
using DelimitedFiles
using Plots.Measures

function make_hist(path, titleArg)
    Arr = readdlm(path)

    if titleArg == "3 Isoforms" && path == "figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_3_isoforms_array.txt"
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        ylabel = "Probability Density",
        normalize=:pdf,
        bins = 100,
        axis = (font(150)), titlefontsize=150, yticks = 0:5:5,
        margin = 35mm,
        left_margin = 50mm,
        right_margin = 50mm)

    elseif titleArg == "Real"
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        normalize=:pdf,
        xlabel = "Mean No. Isoforms\nper Gene per Cell",
        bins = 100,
        axis = (font(150)), titlefontsize=150, yticks = 0:5:5,
        margin = 35mm,
        left_margin = 50mm,
        right_margin = 50mm)

    else
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        normalize=:pdf,
        bins = 100,
        axis = (font(150)), titlefontsize=150, yticks = 0:5:5,
        margin = 35mm,left_margin = 50mm,
        right_margin = 50mm)
    end
    #vline!([mean(Arr)], color = :black, linewidth = 8)
    return p1
end
#for countsMatrix in scnorm_counts
# "Weibull", "Random", "UniformObserved", "CellVariable"
Weibull_1 = make_hist("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
Weibull_2 = make_hist("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
Weibull_3 = make_hist("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
Weibull_4 = make_hist("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
Weibull_real = make_hist("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_real_isoforms_array.txt", "Real")

Random_1 = make_hist("figures/figure5_data/Random_H1_96_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
Random_2 = make_hist("figures/figure5_data/Random_H1_96_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
Random_3 = make_hist("figures/figure5_data/Random_H1_96_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
Random_4 = make_hist("figures/figure5_data/Random_H1_96_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
Random_real = make_hist("figures/figure5_data/Random_H1_96_scnorm_counts.txt_real_isoforms_array.txt", "Real")

UniformObserved_1 = make_hist("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
UniformObserved_2 = make_hist("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
UniformObserved_3 = make_hist("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
UniformObserved_4 = make_hist("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
UniformObserved_real = make_hist("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_real_isoforms_array.txt", "Real")

CellVariable_1 = make_hist("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
CellVariable_2 = make_hist("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
CellVariable_3 = make_hist("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
CellVariable_4 = make_hist("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
CellVariable_real = make_hist("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_real_isoforms_array.txt", "Real")

p3 = plot( Weibull_1, Random_1, UniformObserved_1, CellVariable_1,
       Weibull_2, Random_2,UniformObserved_2, CellVariable_2,
       Weibull_3, Random_3,UniformObserved_3, CellVariable_3,
       Weibull_4, Random_4,UniformObserved_4, CellVariable_4,
       Weibull_real, Random_real,UniformObserved_real, CellVariable_real,
       layout = (5,4), legend = false, size = (750 *8, 750*8), axis = (font(90)), titlefontsize=90,
       margin = 45mm)

savefig(p3, "figures/supplementary_figs/figure2.png")

###############################################################
using JLD2
using DelimitedFiles
using HypothesisTests
using DataFrames

Weibull_1 = readdlm("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_1_isoforms_array.txt")
Weibull_2 = readdlm("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_2_isoforms_array.txt")
Weibull_3 = readdlm("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_3_isoforms_array.txt")
Weibull_4 = readdlm("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_4_isoforms_array.txt")
Weibull_real = readdlm("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_real_isoforms_array.txt")

Random_1 = readdlm("figures/figure5_data/Random_H1_96_scnorm_counts.txt_1_isoforms_array.txt")
Random_2 = readdlm("figures/figure5_data/Random_H1_96_scnorm_counts.txt_2_isoforms_array.txt")
Random_3 = readdlm("figures/figure5_data/Random_H1_96_scnorm_counts.txt_3_isoforms_array.txt")
Random_4 = readdlm("figures/figure5_data/Random_H1_96_scnorm_counts.txt_4_isoforms_array.txt")
Random_real = readdlm("figures/figure5_data/Random_H1_96_scnorm_counts.txt_real_isoforms_array.txt")

UniformObserved_1 = readdlm("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_1_isoforms_array.txt")
UniformObserved_2 = readdlm("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_2_isoforms_array.txt")
UniformObserved_3 = readdlm("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_3_isoforms_array.txt")
UniformObserved_4 = readdlm("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_4_isoforms_array.txt")
UniformObserved_real = readdlm("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_real_isoforms_array.txt")

CellVariable_1 = readdlm("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_1_isoforms_array.txt")
CellVariable_2 = readdlm("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_2_isoforms_array.txt")
CellVariable_3 = readdlm("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_3_isoforms_array.txt")
CellVariable_4 = readdlm("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_4_isoforms_array.txt")
CellVariable_real = readdlm("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_real_isoforms_array.txt")

df = DataFrame(IsoformsSimulated=Int64[], Comparison = String[], pValue = Float64[])

append!(df, DataFrame(IsoformsSimulated = [1], Comparison = ["All"], pValue=[pvalue(KSampleADTest(vec(Weibull_1), vec(Random_1), vec(UniformObserved_1), vec(CellVariable_1)))]))
append!(df, DataFrame(IsoformsSimulated = [2], Comparison = "All", pValue=[pvalue(KSampleADTest(vec(Weibull_2), vec(Random_2), vec(UniformObserved_2), vec(CellVariable_2)))]))
append!(df, DataFrame(IsoformsSimulated = [3], Comparison = "All", pValue=[pvalue(KSampleADTest(vec(Weibull_3), vec(Random_3), vec(UniformObserved_3), vec(CellVariable_3)))]))
append!(df, DataFrame(IsoformsSimulated = [4], Comparison = "All", pValue=[pvalue(KSampleADTest(vec(Weibull_4), vec(Random_4), vec(UniformObserved_4), vec(CellVariable_4)))]))

append!(df, DataFrame(IsoformsSimulated = [1], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest(vec(UniformObserved_1), vec(CellVariable_1)))]))
append!(df, DataFrame(IsoformsSimulated = [2], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest( vec(UniformObserved_2), vec(CellVariable_2)))]))
append!(df, DataFrame(IsoformsSimulated = [3], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest(vec(UniformObserved_3), vec(CellVariable_3)))]))
append!(df, DataFrame(IsoformsSimulated = [4], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest(vec(UniformObserved_4), vec(CellVariable_4)))]))

jldopen("figures/figure4_data/Supfig2_AD_test_results.jld2", "w") do file
    file["df"] = df
end

###########################################
# overlap
#############################################

function make_hist(path, titleArg)

     Arr = readdlm(path)
     p1 = histogram(Arr, title = titleArg, normalize=:pdf, bins = 100,
     ylims = (0,4), color = palette(:default)[2],
     xlims = (0,1),
     axis = (font(150)), titlefontsize=150, yticks = 0:2:4,
     xticks = 0:0.5:1,
     margin = 35mm,left_margin = 50mm,
     right_margin = 50mm)
     return p1

end

Weibull_1 = make_hist("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
Weibull_2 = make_hist("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
Weibull_3 = make_hist("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
Weibull_4 = make_hist("figures/figure5_data/Weibull_H1_96_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

Random_1 = make_hist("figures/figure5_data/Random_H1_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
Random_2 = make_hist("figures/figure5_data/Random_H1_96_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
Random_3 = make_hist("figures/figure5_data/Random_H1_96_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
Random_4 = make_hist("figures/figure5_data/Random_H1_96_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

UniformObserved_1 = make_hist("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
UniformObserved_2 = make_hist("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
UniformObserved_3 = make_hist("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
UniformObserved_4 = make_hist("figures/figure5_data/UniformObserved_H1_96_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

CellVariable_1 = make_hist("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
CellVariable_2 = make_hist("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
CellVariable_3 = make_hist("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
CellVariable_4 = make_hist("figures/figure5_data/CellVariable_H1_96_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

p3 = plot( Weibull_1, Random_1, UniformObserved_1, CellVariable_1,
       Weibull_2, Random_2,UniformObserved_2, CellVariable_2,
       Weibull_3, Random_3,UniformObserved_3, CellVariable_3,
       Weibull_4, Random_4,UniformObserved_4, CellVariable_4,
       layout = (4,4), legend = false, size = (750 *8, 750*8), axis = (font(90)), titlefontsize=90,
       margin = 45mm)

savefig(p3, "figures/supplementary_figs/figure2_overlap.png")

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
#!/mnt/e/Julia-1.6.2
#########################################################################
# Domain_Parser, Version 1.0.0
# 2021, ShanghaiTech University, Wang Lin
# Email: wanglin3@shanghaitech.edu.cn
#########################################################################

using Base: Float64, Float32
using BioStructures

function get_next_line_frag(chain, frag_num, frag_start::Int64, chain_res_end, frag_matrix, pdb_name, model_num, chain_id)
    truncate_domain = collectatoms(chain, res -> resnumber(res) >= frag_start, calphaselector)
    frag_checkpoints::Int64 = 0
    update_matrix = false
    first_res = true
    residue_num_last::Int64 = 0
    frag_end::Int64 = 0
    for residue in truncate_domain
        b_factor_p::Float32 = tempfactor(residue)
        residue_num_p = resnumber(residue)
        chain_res_end == residue_num_p && begin
            frag_end = residue_num_p
            frag_checkpoints > 0 && ( update_matrix = true )
            break
        end
        first_res == true || begin
            if residue_num_last != residue_num_p - 1
                frag_end = residue_num_last
                (frag_checkpoints > 0) && ( update_matrix = true )
                println("Note: Gap $residue_num_last - $residue_num_p in $pdb_name-$model_num-$chain_id")
                break
            end
        end
        residue_num_last = residue_num_p
        if b_factor_p >= b_factor_min && b_factor_p <= b_factor_max
            ( frag_checkpoints == 0 ) && ( frag_start = residue_num_p )
            ( frag_checkpoints <= disorder_len_max ) && ( frag_checkpoints = frag_checkpoints + 1 )
        else
            if frag_checkpoints > 0
                frag_checkpoints = frag_checkpoints - 1
                ( frag_checkpoints == 0 ) && begin
                    frag_end = residue_num_p
                    update_matrix = true
                    break
                end
            end
        end
        first_res = false
    end
    ( update_matrix == true ) && begin
        frag_len = frag_end - frag_start + 1
        frag_len >= domain_min_len && ( push!(frag_matrix,"$frag_num,$frag_start,$frag_end") )
    end
    return frag_matrix, frag_end
end

function get_terminal_resnum(chain::Chain)
    C_end::Int64 = 1
    N_start::Int64 = 9999999999
    for residue in chain
        residue_num_p::Int64 = resnumber(residue)
        ( C_end <= residue_num_p ) && ( C_end = residue_num_p )
        ( N_start >= residue_num_p ) && ( N_start = residue_num_p )
    end
    return N_start::Int64,C_end::Int64
end

function segment_domain(chain, frag_matrix, pdb_name, model_num, chain_id,index_dir,pdb,output_dir)
    res_domain = length(frag_matrix)
    domain_ID = 0
    domain_info = open("$index_dir/$pdb_name-$model_num-$chain_id-domain.info","w")
    while res_domain > 0
        domain = ""
        domain_ID+=1
        got_domain = []
        frag_id = 0
        for frag in frag_matrix
            frag_id+=1
            frag_start = parse(Int64, split(frag, ",")[2])
            frag_end = parse(Int64, split(frag, ",")[3])
            frag_str = collectresidues(chain, res -> frag_start <= resnumber(res) <= frag_end )
            ( frag_id == 1 ) && begin
                domain = "$frag_start-$frag_end"
                push!(got_domain, frag)
            end
            ( frag_id != 1 ) && begin
                for frag_2 in split(domain,",")
                    frag_2_str = collectresidues(chain, res -> parse(Int64,split(frag_2,"-")[1]) <= resnumber(res) <= parse(Int64,split(frag_2,"-")[2]))
                    Adjacent_residues = 0
                    for residue in frag_2_str
                        frag_dist::Float64 = BioStructures.distance(frag_str, residue, allselector)
                        ( frag_dist <= min_dist_interdomain ) && ( Adjacent_residues+=1 )
                    end
                    ( min_Adjacent_residues <= Adjacent_residues ) && begin
                        domain = "$(domain),$frag_start-$frag_end"
                        push!(got_domain, frag)
                        break
                    end
                end
            end
        end
        for frag in got_domain
            filter!(x->x!=frag, frag_matrix)
        end
        res_domain = length(frag_matrix)
        write(domain_info, "$domain_ID:$domain\n")
    end
    close(domain_info)
    for domain_ddd in readlines("$index_dir/$pdb_name-$model_num-$chain_id-domain.info")
        domain_id = split(domain_ddd,":")[1]
        domain_vec = split(domain_ddd,":")[2]
        com_id = 0
        domain_selection = []
        for domain_range in split(domain_vec,",")
            frag_N = parse(Int64,split(domain_range,"-")[1])
            frag_C = parse(Int64,split(domain_range,"-")[2])
            for resi in frag_N:frag_C
                push!(domain_selection, resi)
            end
        end
        domain_str = collectresidues(chain, res -> resnumber(res) in domain_selection , allselector)
        writepdb("$output_dir/$pdb_name-$model_num-$chain_id-$domain_id.pdb", domain_str)
    end
end

function main(pdb::String, index_dir::String, output_dir::String)
    Structure = read(pdb, PDB)
    pdb_name = chop(basename(pdb), tail = 4)
    for model in Structure
        model_num = modelnumber(model)
        for chain in model
            chain_id = chainid(chain)
            chain_res_start::Int64, chain_res_end::Int64= get_terminal_resnum(chain)
            frag_start::Int64 = chain_res_start
            frag_end::Int64 = chain_res_start
            frag_num = 0
            frag_matrix = []
            while frag_end < chain_res_end
                frag_num+=1
                frag_start = frag_end
                #将 frag_num, frag_start, frag_end 存入 frag_matrix
                frag_matrix, frag_end= get_next_line_frag(chain, frag_num, frag_start, chain_res_end, frag_matrix, pdb_name, model_num, chain_id)
            end
            segment_domain(chain, frag_matrix, pdb_name, model_num, chain_id,index_dir,pdb,output_dir)
        end
    end
end

println("Name of Input Dir for protein pdb files:")
const input_dir=readline(stdin)
base_input = basename(input_dir)
index_dir="$(base_input)-INDEX"
mkdir(index_dir)
println("Output Dir for domain Structure:")
const output_dir=readline(stdin)
global b_factor_max = 100 # max b_factor(pLDDT) to keep a residue.
global b_factor_min = 50 # min b_factor(pLDDT) to keep a residue.
global disorder_len_max = 5 # the max disorder tail in a ordered domain.
global domain_min_len = 10   # this option is means of the min length of fargs. The domain is consist of frags incessant in sequence.
global min_dist_interdomain = 3.50  # 3.00 is a very strict cutoff.
global min_Adjacent_residues = 5 # How many residues are close together between frags will be identified as a domain.

Threads.@threads for pdb_file in readdir(input_dir)
    pdb = "$(input_dir)/$(pdb_file)"
    endswith(pdb, ".pdb") && begin
        try
            main(pdb, index_dir, output_dir)
        catch e1
            println("$e1 Error: $pdb")
        end
    end
end


text:
using DelimitedFiles
using HDF5
using Printf
using StatsBase

include("load_data.jl")

"""
set (overwrite) a description as an attribute for a hdf5 object
"""
function myh5desc(filename, name, description::String)
    name = replace(name, r"\/+" => "/")
    filename = replace(filename, r"\/+" => "/")
    fid = h5open(filename, "cw")
    obj = fid[name]
    if haskey(attributes(obj), "description")
        delete_attribute(obj, "description")
    end
    attributes(obj)["description"] = description
    close(obj)
    close(fid)
end

"""
wrapper for h5write to replace dataset if existing (replace currently does not work)
"""
function myh5write(filename, datasetname, data::AbstractArray)
    datasetname = replace(datasetname, r"\/+" => "/")
    filename = replace(filename, r"\/+" => "/")

    h5open(filename, "cw") do fid
        if haskey(fid, datasetname)
            d = fid[datasetname]
            d = data
        else
            fid[datasetname] = data
        end
    end
end

"""
wrapper for h5write to directly deal with a StatsBase distribution
"""
function myh5write(filename, datasetname, dist::AbstractHistogram)
    myh5write(filename, datasetname, hcat(collect(dist.edges[1][1:end-1]), dist.weights))
end

"""
wrapper for h5write to directly deal with contact train object
"""
function myh5write(filename, groupname::AbstractString, ets::encounter_trains{T}) where T
    myh5write(filename, @sprintf("%s/ids", groupname), ids(ets))
    myh5desc(filename, groupname,
        "Encoutner trains, where `ids` are labels for agents. Agents are numbered, and each agent has a train: `train_agentnumber` with the times of contacts in seconds."
    )

    for (i, train) in enumerate(trains(ets))
        myh5write(filename, @sprintf("%s/train_%d", groupname, i), train)
    end
end


text:
module DiffEqBiological

using DataStructures, MacroTools, Parameters, Reexport, SparseArrays, SymEngine
using DiffEqBase, DiffEqJump, Compat
@reexport using DiffEqBase, DiffEqJump
using DynamicPolynomials, LinearAlgebra, RecipesBase, Latexify

import HomotopyContinuation

import Base: (==)

const ExprValues = Union{Expr,Symbol,Float64,Int}

include("expression_utils.jl")
include("reaction_network.jl")
include("maketype.jl")
include("network_properties.jl")
include("massaction_jump_utils.jl")
include("problem.jl")
include("equilibrate_utils.jl")
include("plot_recipes.jl")
include("latexify_recipes.jl")

# reaction network macro
export @reaction_network, @reaction_func, @min_reaction_network, @empty_reaction_network

# functions to query network properties
export species, params, speciesmap, paramsmap, numspecies, numreactions, numparams
export oderhsfun, jacfun, paramjacfun, odefun, noisefun, sdefun, jumps, regularjumps
export odeexprs, jacobianexprs, noiseexprs, jumpexprs, rateexpr, oderatelawexpr, ssaratelawexpr
export substratestoich, productstoich, netstoich, ismassaction, dependants, dependents, substrates, products, substratesymstoich, productsymstoich
export rxtospecies_depgraph, speciestorx_depgraph, rxtorx_depgraph

# functions to extend empty_ and min_ reaction_networks
export addspecies!, addparam!, add_scale_noise_param!, addreaction!

# functions to add mathematical equations to the network
export addodes!, addsdes!, addjumps!, addequi1!, addequi2!

# problems that can be solved from the network
export ODEProblem, SDEProblem, DiscreteProblem, JumpProblem, SteadyStateProblem

# tolls for finding equilibrium solutions and bifurcation diagrams.
export EquilibrateContent
export @add_constraint, @add_constraints, internal___add___constraint!, fix_parameters, @make_hc_template, make_hc_template!, @add_hc_template, add_hc_template!
export steady_states, stability
export bifurcations, bifurcation_grid, bifurcation_grid_2d, bifurcation_grid_diagram, SimpleHCBifurcationSolver, HCBifurcationSolver, HCSteadyStateSolver

end # module


text:
"""
    space(::SiteType"Qudit";
          dim = 2,
          conserve_qns = false,
          conserve_number = false,
          qnname_number = "Number")

Create the Hilbert space for a site of type "Qudit".

Optionally specify the conserved symmetries and their quantum number labels.
"""
function ITensors.space(
  ::SiteType"Qudit";
  dim=2,
  conserve_qns=false,
  conserve_number=conserve_qns,
  qnname_number="Number",
)
  if conserve_number
    return [QN(qnname_number, n - 1) => 1 for n in 1:dim]
  end
  return dim
end

function ITensors.val(::ValName{N}, ::SiteType"Qudit") where {N}
  return parse(Int, String(N)) + 1
end

function ITensors.state(::StateName{N}, ::SiteType"Qudit", s::Index) where {N}
  n = parse(Int, String(N))
  st = zeros(dim(s))
  st[n + 1] = 1.0
  return st
end

function _op(::OpName"Id", ::SiteType"Qudit"; dim=2)
  mat = zeros(dim, dim)
  for k in 1:dim
    mat[k, k] = 1.0
  end
  return mat
end

function _op(::OpName"Adag", ::SiteType"Qudit"; dim=2)
  mat = zeros(dim, dim)
  for k in 1:(dim - 1)
    mat[k + 1, k] = √k
  end
  return mat
end
_op(::OpName"adag", st::SiteType"Qudit"; kwargs...) = _op(OpName"Adag"(), st; kwargs...)
_op(::OpName"a†", st::SiteType"Qudit"; kwargs...) = _op(OpName"Adag"(), st; kwargs...)

function _op(::OpName"A", ::SiteType"Qudit"; dim=2)
  mat = zeros(dim, dim)
  for k in 1:(dim - 1)
    mat[k, k + 1] = √k
  end
  return mat
end
_op(::OpName"a", st::SiteType"Qudit"; kwargs...) = _op(OpName"A"(), st; kwargs...)

function _op(::OpName"N", ::SiteType"Qudit"; dim=2)
  mat = zeros(dim, dim)
  for k in 1:dim
    mat[k, k] = k - 1
  end
  return mat
end
_op(::OpName"n", st::SiteType"Qudit"; kwargs...) = _op(OpName"N"(), st; kwargs...)

function ITensors.op(on::OpName, st::SiteType"Qudit", s::Index)
  return itensor(_op(on, st; dim=dim(s)), s', dag(s))
end


text:
kernel = "../src/pca/logreg.jl"

function logreg_loss(v, X, b, λ)
    rv = 0.0
    for i in 1:length(b)
        rv += log(1 + exp(-b[i]*(v[1]+dot(X[:, i], view(v, 2:length(v))))))
    end
    rv / length(b) + λ/2 * norm(v)^2
end

# input data set with known optimal solution
X = [1.444786643000158 0.49236792885913283 -0.53258473265429 0.05476455630673194 -1.3473893605265843; 0.48932299731783646 2.0708445447107926 1.2414596020757043 0.9131934117095984 -0.15692043560721075; 0.7774625331093794 0.7234405608945721 -0.037446104354257874 -1.1104987697394342 1.354975413199728]
b = [-1, 1, -1, -1, 1]
v_opt = [-0.3423591553493419, -0.41317049965033387, 0.007294166575956451, 0.6763846515861628]
m, n = size(X)
λ = 1 / n
opt = logreg_loss(v_opt, X, b, λ)

# write test problem to disk
inputfile = tempname()
inputdataset = "X"
outputdataset = "V"
labeldataset = "b"
h5open(inputfile, "w") do file
    file[inputdataset] = X
    file[labeldataset] = b
end

# GD
nworkers = 2
niterations = 200
stepsize = 0.1
outputfile = tempname()
mpiexec(cmd -> run(```
    $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
    --inputdataset $inputdataset
    --outputdataset $outputdataset
    --niterations $niterations
    --saveiterates
    --lambda $λ
    ```))
vs = load_logreg_iterates(outputfile, outputdataset)
v = vs[end]
f = logreg_loss(v, X, b, λ)
@test f < opt * (1+1e-2)

# same as previous, but with nslow > 0
nworkers = 2
nslow = 1
niterations = 100
stepsize = 0.1
outputfile = tempname()
mpiexec(cmd -> run(```
    $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
    --inputdataset $inputdataset
    --outputdataset $outputdataset
    --niterations $niterations
    --saveiterates
    --lambda $λ
    --nslow $nslow
    ```))
vs = load_logreg_iterates(outputfile, outputdataset)
v = vs[end]
f = logreg_loss(v, X, b, λ)
@test f < opt * (1+1e-2)

# same as previous, but with slowprob. > 0
nworkers = 2
slowprob = 0.5
niterations = 100
stepsize = 0.1
outputfile = tempname()
mpiexec(cmd -> run(```
    $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
    --inputdataset $inputdataset
    --outputdataset $outputdataset
    --niterations $niterations
    --saveiterates
    --lambda $λ
    --slowprob $slowprob
    ```))
vs = load_logreg_iterates(outputfile, outputdataset)
v = vs[end]
f = logreg_loss(v, X, b, λ)
@test f < opt * (1+1e-2)

# DSAG
vralgo = "tree"
nworkers = 2
nwait = 1
niterations = 100
stepsize = 0.1
nsubpartitions = 2
outputfile = tempname()
mpiexec(cmd -> run(```
    $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
    --inputdataset $inputdataset
    --nwait $nwait
    --variancereduced
    --vralgo $vralgo
    --nsubpartitions $nsubpartitions
    --outputdataset $outputdataset
    --niterations $niterations
    --saveiterates
    --lambda $λ
    ```))
vs = load_logreg_iterates(outputfile, outputdataset)
v = vs[end]
f = logreg_loss(v, X, b, λ)
@test f < opt * (1+1e-2)

# DSAG w. nwaitschedule < 1.0
nworkers = 2
nwait = 2
niterations = 100
stepsize = 0.1
nsubpartitions = 2
nwaitschedule = 0.9
outputfile = tempname()
mpiexec(cmd -> run(```
    $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
    --inputdataset $inputdataset
    --nwait $nwait
    --variancereduced
    --vralgo $vralgo
    --nsubpartitions $nsubpartitions
    --outputdataset $outputdataset
    --niterations $niterations
    --saveiterates
    --lambda $λ
    --nwaitschedule $nwaitschedule
    ```))
vs = load_logreg_iterates(outputfile, outputdataset)
v = vs[end]
f = logreg_loss(v, X, b, λ)
@test f < opt * (1+1e-2)

# DSAG w. sparse input data
inputfile = tempname()
h5open(inputfile, "cw") do file
    H5SparseMatrixCSC(file, inputdataset, sparse(X))
    file[labeldataset] = b
    flush(file)
end

nworkers = 2
nwait = 1
niterations = 100
stepsize = 0.1
nsubpartitions = 2
outputfile = tempname()
mpiexec(cmd -> run(```
    $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
    --inputdataset $inputdataset
    --nwait $nwait
    --variancereduced
    --vralgo $vralgo
    --nsubpartitions $nsubpartitions
    --outputdataset $outputdataset
    --niterations $niterations
    --lambda $λ
    ```))
vs = load_logreg_iterates(outputfile, outputdataset)
v = vs[end]
f = logreg_loss(v, X, b, λ)
@test f < opt * (1+1e-2)


text:
using MaxwellGuidedMode
using Test

@testset "MaxwellGuidedMode.jl" begin
    # Write your tests here.
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiPolicyV1beta1Eviction <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    deleteOptions::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1DeleteOptions } # spec name: deleteOptions
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata

    function IoK8sApiPolicyV1beta1Eviction(;apiVersion=nothing, deleteOptions=nothing, kind=nothing, metadata=nothing)
        o = new()
        validate_property(IoK8sApiPolicyV1beta1Eviction, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiPolicyV1beta1Eviction, Symbol("deleteOptions"), deleteOptions)
        setfield!(o, Symbol("deleteOptions"), deleteOptions)
        validate_property(IoK8sApiPolicyV1beta1Eviction, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiPolicyV1beta1Eviction, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        o
    end
end # type IoK8sApiPolicyV1beta1Eviction

const _property_map_IoK8sApiPolicyV1beta1Eviction = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("deleteOptions")=>Symbol("deleteOptions"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"))
const _property_types_IoK8sApiPolicyV1beta1Eviction = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("deleteOptions")=>"IoK8sApimachineryPkgApisMetaV1DeleteOptions", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta")
Base.propertynames(::Type{ IoK8sApiPolicyV1beta1Eviction }) = collect(keys(_property_map_IoK8sApiPolicyV1beta1Eviction))
Swagger.property_type(::Type{ IoK8sApiPolicyV1beta1Eviction }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiPolicyV1beta1Eviction[name]))}
Swagger.field_name(::Type{ IoK8sApiPolicyV1beta1Eviction }, property_name::Symbol) =  _property_map_IoK8sApiPolicyV1beta1Eviction[property_name]

function check_required(o::IoK8sApiPolicyV1beta1Eviction)
    true
end

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


text:
# ---
# title: 497. Random Point in Non-overlapping Rectangles
# id: problem497
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search, Random
# link: <https://leetcode.com/problems/random-point-in-non-overlapping-rectangles/description/>
# hidden: true
# ---
#
# Given a list of **non-overlapping**  axis-aligned rectangles `rects`, write a
# function `pick` which randomly and uniformily picks an **integer point** in
# the space covered by the rectangles.
#
# Note:
#
#   1. An **integer point**  is a point that has integer coordinates.
#   2. A point on the perimeter of a rectangle is  **included** in the space covered by the rectangles.
#   3. `i`th rectangle = `rects[i]` = `[x1,y1,x2,y2]`, where `[x1, y1]` are the integer coordinates of the bottom-left corner, and `[x2, y2]` are the integer coordinates of the top-right corner.
#   4. length and width of each rectangle does not exceed `2000`.
#   5. `1 <= rects.length <= 100`
#   6. `pick` return a point as an array of integer coordinates `[p_x, p_y]`
#   7. `pick` is called at most `10000` times.
#
# **Example 1:**
#
#
#
#     Input: ["Solution","pick","pick","pick"]
#     [[[[1,1,5,5]]],[],[],[]]
#     Output: [null,[4,1],[4,1],[3,3]]
#
#
# **Example 2:**
#
#
#
#     Input: ["Solution","pick","pick","pick","pick","pick"]
#     [[[[-2,-2,-1,-1],[1,0,3,0]]],[],[],[],[],[]]
#     Output: [null,[-1,-2],[2,0],[-2,-1],[3,0],[-2,-2]]
#
# **Explanation of Input Syntax:**
#
# The input is two lists: the subroutines called and their arguments.
# `Solution`'s constructor has one argument, the array of rectangles `rects`.
# `pick` has no arguments. Arguments are always wrapped with a list, even if
# there aren't any.
#
#
#
#
## @lc code=start
using LeetCode

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


text:
struct DescriptorPoolState
    allocated::Dictionary{DescriptorType,Int}
    limits::Dictionary{DescriptorType,Int}
    max::Int
end

"""
Entity which manages descriptor pools and descriptor sets.
"""
struct DescriptorAllocator
    device::Device
    pools::Dictionary{Created{DescriptorPool,DescriptorPoolCreateInfo},DescriptorPoolState}
end

DescriptorAllocator(device::Device) = DescriptorAllocator(device, Dictionary())

"""
Free descriptor sets from a descriptor allocator.

!!! warning
    All descriptor sets must have been allocated with the same pool, through the same descriptor allocator.
    They may, however, have been allocated in several times.
"""
function free_descriptor_sets!(da::DescriptorAllocator, sets::Vector{Created{DescriptorSet,DescriptorSetAllocateInfo}})
    pool = unique!(getproperty.(info.(sets), :descriptor_pool))
    length(pool) == 1 || error("The descriptor sets were allocated via different descriptor pools.")
    allocations = compute_allocations(da, sets)
    da.pools[pool].allocated .-= allocations
    free_descriptor_sets(da.device, pool, sets)
end

function new_pool!(da::DescriptorAllocator, allocations::Dictionary{DescriptorType,Int}; next=C_NULL, flags=0)
    limits = [DescriptorPoolSize(t, max(allocations[t], 20)) for t in keys(allocations)]
    maxsize = max(sum(allocations), 50)
    create_info = DescriptorPoolCreateInfo(limits, maxsize; next, flags)
    new = Created(DescriptorPool(da.device, create_info), create_info)
    insert!(da, new, DescriptorPoolState(Dictionary(), limits, maxsize))
    new
end

function compute_allocations(da::DescriptorAllocator, layouts::AbstractVector{Created{DescriptorSetLayout, DescriptorSetLayoutCreateInfo}})
    allocations = Dictionary{DescriptorType,Int}()
    for layout in layouts
        for binding in layout.info.bindings
            allocations[binding.descriptor_type] = get(allocations, binding.descriptor_type, 0) + binding.descriptor_count
        end
    end
    allocations
end

function compute_allocations(da::DescriptorAllocator, sets::AbstractVector{Created{DescriptorSet,DescriptorSetAllocateInfo}})
    merge!(+, compute_allocations.(getproperty.(unique!(info.(sets)), :set_layouts))...)
end

function find_pool!(da::DescriptorAllocator, layouts)
    find_pool!(da, compute_allocations(da, layouts))
end

function find_pool!(da::DescriptorAllocator, allocations::Dictionary{DescriptorType,Int})
    for (pool, state) in zip(keys(da.pools), da.pools)
        n = 0
        has_room = sum(allocations) < sum(state.max) && all(keys(state.allocated)) do t
            state.allocated[t] + allocations[t] < state.limits[t]
        end
        if has_room
            return pool
        end
    end

    # no pool could be reused
    new_pool!(da, allocations)
end

function allocate_descriptor_sets!(da::DescriptorAllocator, layouts::Vector{Created{DescriptorSetLayout,DescriptorSetLayoutCreateInfo}}; next=C_NULL, flags=0)
    pool = find_pool!(da, layouts)
    info = DescriptorSetAllocateInfo(pool, layouts; next, flags)
    allocate_descriptor_sets!(da, pool, info)
end

function allocate_descriptor_sets!(da::DescriptorAllocator, pool::CommandPool, info::DescriptorSetAllocateInfo)
    pool in handle.(keys(da.pools)) || error("Pool $pool was not created with the provided allocator $da.")
    allocations = compute_allocations(da, info.set_layouts)
    merge!(pool.allocated, allocations)
    sets = unwrap(allocate_descriptor_sets(da.device, info))
    Created.(sets, info)
end


text:
module SLAM
export SlamManager, add_image!, add_stereo_image!, get_queue_size
export Params, Camera, run!, to_cartesian, reset!
export Visualizer, ReplaySaver
export set_frame_wc!, process_frame_wc!, set_image!, set_position!

using BSON: @save, @load
using OrderedCollections: OrderedSet, OrderedDict
using GLMakie
using Interpolations
using Images
using ImageDraw
using ImageFeatures
using LeastSquaresOptim
using LinearAlgebra
using Manifolds
using Random
using RecoverPose
using Rotations
using StaticArrays
using SparseArrays
using SparseDiffTools

const Point2 = SVector{2}
const Point2i = SVector{2, Int64}
const Point2f = SVector{2, Float64}
const Point3f = SVector{3, Float64}
const Point3f0 = SVector{3, Float32}

const SE3 = SpecialEuclidean(3)

@inline convert(x::Point2f)::Point2i = x .|> round .|> Int64
@inline convert(x::Vector{Point2f}) = Point2i[xi .|> round .|> Int64 for xi in x]

@inline to_homogeneous(p::SVector{3, T}) where T = SVector{4, T}(p..., one(T))
@inline to_homogeneous(p::SVector{4}) = p

"""
Params:
    x::Point2 Point to convert to CartesianIndex in (row, col) format.
"""
@inline to_cartesian(x) = CartesianIndex(convert(x)...)
@inline function to_cartesian(x::Point2, cell_size::Int64)
    x = convert(x) .÷ cell_size .+ 1
    CartesianIndex(x...)
end

function to_4x4(m::SMatrix{3, 3, T, 9}) where T
    SMatrix{4, 4, T}(
        m[1, 1], m[2, 1], m[3, 1], 0,
        m[1, 2], m[2, 2], m[3, 2], 0,
        m[1, 3], m[2, 3], m[3, 3], 0,
        0,       0,       0,       1)
end
function to_4x4(m::SMatrix{3, 4, T, 12}) where T
    SMatrix{4, 4, T}(
        m[1, 1], m[2, 1], m[3, 1], 0,
        m[1, 2], m[2, 2], m[3, 2], 0,
        m[1, 3], m[2, 3], m[3, 3], 0,
        m[1, 4], m[2, 4], m[3, 4], 1)
end
function to_4x4(m, t)
    SMatrix{4, 4, Float64}(
        m[1, 1], m[2, 1], m[3, 1], 0,
        m[1, 2], m[2, 2], m[3, 2], 0,
        m[1, 3], m[2, 3], m[3, 3], 0,
        t[1],    t[2],    t[3],    1)
end

include("optical_flow/utils.jl")
include("optical_flow/pyramid.jl")
include("optical_flow/lucas_kanade.jl")

include("camera.jl")
include("extractor.jl")
include("tracker.jl")
include("params.jl")
include("frame.jl")
include("motion_model.jl")
include("map_point.jl")
include("map_manager.jl")
include("front_end.jl")
include("estimator.jl")
include("mapper.jl")
include("io/visualizer.jl")
include("io/saver.jl")
include("bundle_adjustment.jl")

"""
```julia
SlamManager(
    params::Params, camera::Camera;
    right_camera::Union{Nothing, Camera} = nothing,
)
```

Slam Manager that is the highest level component in the system.
It is responsible for sending new frames to the other components
and for processing their outputs.

**Note**, that upon creating, SlamManager launches Mapper
in the separate thread.

# Arguments

- `params::Params`: Parameters of the system.
- `image_queue::Vector{Matrix{Gray{Float64}}}`: Queue of the images to be
    processed.
- `right_image_queue::Vector{Matrix{Gray{Float64}}}`: In case of stereo mode,
    queue of images for the right camera.
    It should be in sync with `image_queue` which in this case is for the left
    camera.
- `time_queue::Vector{Float64}`: Queue of timestamps for each of the frame.
    The timestamps are used in the motion model to predict next pose for the
    frame, before it is actually computed.
- `current_frame::Frame`: Current frame that is processed. It is shared among
    all other components in the system.
- `frame_id::Int64`: Id of the current frame.
- `front_end::FrontEnd`: Front-End component that is used for tracking.
- `map_manager::MapManager`: Map manager for the managment of
    keyframes & mappoints.
- `mapper::Mapper`: Mapper that is used for triangulation of keypoints.
    It is launched in the constructor as a separate thread.
- `extractor::Extractor`: Used in extraction of keypoints from the frames.
- `exit_required::Bool`: Set it to `true` to stop SlamManager
    once it is launched.
"""
mutable struct SlamManager
    params::Params

    image_queue::Vector{Matrix{Gray{Float64}}}
    right_image_queue::Vector{Matrix{Gray{Float64}}}
    time_queue::Vector{Float64}

    current_frame::Frame
    frame_id::Int64

    front_end::FrontEnd
    map_manager::MapManager
    mapper::Mapper
    extractor::Extractor

    visualizer::Union{Nothing, Visualizer, ReplaySaver}

    exit_required::Bool

    mapper_thread
    image_lock::ReentrantLock
end

function SlamManager(
    params, camera; right_camera = nothing, visualizer = nothing,
)
    params.stereo && right_camera ≡ nothing &&
        error("[SM] Provide `right_camera` when in stereo mode.")

    avoidance_radius = max(5, params.max_distance ÷ 2)
    image_resolution = (camera.height, camera.width)
    grid_resolution = ceil.(Int64, image_resolution ./ params.max_distance)

    image_queue = Matrix{Gray{Float64}}[]
    right_image_queue = Matrix{Gray{Float64}}[]
    time_queue = Float64[]

    frame = Frame(;camera, right_camera, cell_size=params.max_distance)
    extractor = Extractor(
        params.max_nb_keypoints, avoidance_radius,
        grid_resolution, params.max_distance)

    map_manager = MapManager(params, frame, extractor)
    front_end = FrontEnd(params, frame, map_manager)

    mapper = Mapper(params, map_manager, frame)
    mapper_thread = Threads.@spawn run!(mapper)
    @debug "[SM] Launched mapper thread."

    SlamManager(
        params, image_queue, right_image_queue,
        time_queue, frame, frame.id,
        front_end, map_manager, mapper, extractor,
        visualizer,
        false, mapper_thread, ReentrantLock())
end

"""
```julia
run!(sm::SlamManager)
```

Main routine for the SlamManager. It runs until `exit_required` variable
is set to `true`. After that, it will end its work, wait for other threads
and finish.

Once there is a frame in the queue, it will get it and first send it to the
FrontEnd for tracking. If FrontEnd requires a new Keyframe, then it will also
send it to the mapper thread for Keyframe creation and triangulation.
"""
function run!(sm::SlamManager)
    image::Union{Nothing, Matrix{Gray{Float64}}} = nothing
    right_image::Union{Nothing, Matrix{Gray{Float64}}} = nothing
    time = 0.0

    while !sm.exit_required
        if sm.params.stereo
            image, right_image, time = get_stereo_image!(sm)
        else
            image, time = get_image!(sm)
        end
        if image ≡ nothing
            sleep(1e-2)
            continue
        end

        sm.frame_id += 1
        sm.current_frame.id = sm.frame_id
        sm.current_frame.time = time
        @debug "[SM] Frame $(sm.frame_id) @ $time."

        is_kf_required = track!(sm.front_end, image, time, sm.visualizer)
        if sm.params.reset_required
            reset!(sm)
            continue
        end

        is_kf_required || continue
        try
            add_new_kf!(sm.mapper, KeyFrame(
                sm.current_frame.kfid,
                sm.params.stereo ? deepcopy(sm.front_end.current_pyramid) : nothing,
                sm.params.stereo ? right_image : nothing))
        catch e
            showerror(stdout, e)
            display(stacktrace(catch_backtrace()))
        end
        sleep(1e-2)
    end

    sm.mapper.exit_required = true
    wait(sm.mapper_thread)
    @debug "[SM] Exit required."
end

"""
```julia
add_image!(sm::SlamManager, image, time)
```

Add monocular image and its timestamp to the queue.
"""
function add_image!(sm::SlamManager, image, time)
    lock(sm.image_lock) do
        push!(sm.image_queue, image)
        push!(sm.time_queue, time)
        @debug "[SM] Image queue size $(length(sm.image_queue))."
    end
end

"""
```julia
add_stereo_image!(sm::SlamManager, image, right_image, time)
```

Add stereo image and its timestamp to the queue.
"""
function add_stereo_image!(sm::SlamManager, image, right_image, time)
    lock(sm.image_lock) do
        push!(sm.image_queue, image)
        push!(sm.right_image_queue, right_image)
        push!(sm.time_queue, time)
        @debug "[SM] Stereo image queue size $(length(sm.image_queue))."
    end
end

"""
```julia
get_image!(sm::SlamManager)
```

Get monocular image and its timestamp from the queue.

# Returns:

`(image, timestamp)` if there is an image in the queue.
Otherwise `(nothing, nothing)`.
"""
function get_image!(sm::SlamManager)
    lock(sm.image_lock) do
        isempty(sm.image_queue) && return nothing, nothing
        image = popfirst!(sm.image_queue)
        time = popfirst!(sm.time_queue)
        image, time
    end
end

"""
```julia
get_stereo_image!(sm::SlamManager)
```

Get stereo image and its timestamp from the queue.

# Returns:

`(image, right_image, timestamp)` if there is an image in the queue.
Otherwise `(nothing, nothing, nothing)`.
"""
function get_stereo_image!(sm::SlamManager)
    lock(sm.image_lock) do
        (isempty(sm.image_queue) || isempty(sm.right_image_queue)) &&
            return nothing, nothing, nothing

        image = popfirst!(sm.image_queue)
        right_image = popfirst!(sm.right_image_queue)
        time = popfirst!(sm.time_queue)
        image, right_image, time
    end
end

"""
```julia
get_queue_size(sm::SlamManager)
```

Get size of the queue of images to be processed.
"""
function get_queue_size(sm::SlamManager)
    lock(sm.image_lock) do
        length(sm.image_queue)
    end
end

"""
```julia
reset!(sm::SlamManager)
```

Reset slam manager, front-end and map_manager.
"""
function reset!(sm::SlamManager)
    @warn "[Slam Manager] Reset required."
    sm.params |> reset!
    sm.current_frame |> reset!
    sm.front_end |> reset!
    sm.map_manager |> reset!
    @warn "[Slam Manager] Reset applied."
end

function draw_keypoints!(
    image::Matrix{T}, frame::Frame; right::Bool = false,
) where T <: RGB
    radius = 2
    for kp in values(frame.keypoints)
        right && !kp.is_stereo && continue

        pixel = (right && kp.is_stereo) ? kp.right_pixel : kp.pixel
        in_image(frame.camera, pixel) || continue

        color = kp.is_3d ? T(0, 0, 1) : T(0, 1, 0)
        kp.is_retracked && (color = T(1, 0, 0);)
        draw!(image, CirclePointRadius(to_cartesian(pixel), radius), color)
    end
    image
end

end


text:
#=
# 2D Poisson, mixed bc
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("mixedbc");

# Try making an optional log
@useLog("mixedbclog")

# Set up the configuration (order doesn't matter)
@domain(2)                          # dimension
@solver(CG)                         # Use CG solver
@functionSpace(LEGENDRE, 3)         # basis function, order
@nodes(LOBATTO)                     # elemental node arrangement

# Specify the problem
@mesh(QUADMESH, 10, 4)               # build uniform mesh. 2nd arg=# of elements, (optional)3rd arg=# of BIDs

@variable(u)                        # same as @variable(u, SCALAR)

@testSymbol(v)                      # sets the symbol for a test function

@boundary(u, 1, NEUMANN, "-2*pi*cos(pi/4) * sin(2*pi*(y+0.125))")       # x=0
@boundary(u, 2, NEUMANN, "2*pi*cos(2*pi*1.125) * sin(2*pi*(y+0.125))")  # x=1
@boundary(u, 3, DIRICHLET, "sin(pi/4) * sin(2*pi*(x+0.125))")           # y=0
@boundary(u, 4, DIRICHLET, "sin(2*pi*1.125) * sin(2*pi*(x+0.125))")     # y=1

# Write the weak form
@coefficient(f, "-8*pi*pi*sin(2*pi*(x+0.125))*sin(2*pi*(y+0.125))")
@weakForm(u, "-grad(u)*grad(v) - f*v")

solve(u);

# exact solution
# check error
maxerr = 0;
exact(x,y) = sin(2*pi*(x+0.125))*sin(2*pi*(y+0.125));

for i=1:size(Finch.grid_data.allnodes,2)
    x = Finch.grid_data.allnodes[1,i];
    y = Finch.grid_data.allnodes[2,i];
    err = abs(u.values[i] - exact(x,y));
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

# solution is stored in the variable's "values"
#using Plots
#pyplot();
#display(plot(Finch.grid_data.allnodes[:,1], Finch.grid_data.allnodes[:,2], u.values[:], st=:surface, reuse=false))

# check
log_dump_config();
log_dump_prob();

@finalize()

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
# default plot test

init_plot = plot!()
test_plot = plot()

@test isa(init_plot, SimplePlots.SimplePlot)
@test isa(test_plot, SimplePlots.SimplePlot)

test_is_same_object(init_plot, test_plot)

# lines and scatter test

plot_1 = plot(rand(10))
@test plot_1.index == 1

plot_2 = scatter!(3:7, 10.0 .^ (-2:+2))
@test plot_2.index == 2

test_is_same_object(plot_1, plot_2)


text:
# parameters for the case study

# model and data generation files
const modelFile = "resource_20_30.jl"
const dataFile = "demand_model_iid.jl"

# parameters for estimating the true optimal objective value using SAA
const numMCScenarios = 1000

# degree of nonlinearity in demand model
const degree = Float64[1]

# scaling factor for the additive demand errors
const demand_errors_scaling = 20.0

# starting seed
const startingSeed = 1280

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA