text:
using Juno
import JuliaInterpreter
using CodeTracking
using CodeTools
import Atom: basepath, handle

const _breakpoints = Dict{Int, Any}()
const _conditions = Dict{Any, Any}()

const _compiledMode = Ref{Any}(JuliaInterpreter.finish_and_return!)

"""
    allbreakpoints()

Updates internal `_breakpoints` dictionary and returns all breakpoints in a vector of dictionaries.
"""
function allbreakpoints()
  bps = JuliaInterpreter.breakpoints()
  simple_bps = []
  empty!(_breakpoints)
  id = 1

  for (bp, cond) in _conditions
    bp in bps && continue
    delete!(_conditions, bp)
  end

  for bp in filter!(bp -> bp isa JuliaInterpreter.BreakpointSignature, bps)
    sbp = simple_breakpoint(bp, id = id)
    if sbp === nothing
      @warn "Not good. Weird breakpoint encountered."
    end
    _breakpoints[id] = bp
    push!(simple_bps, sbp)
    id += 1
  end

  simple_bps
end

function simple_breakpoint(bp::JuliaInterpreter.BreakpointSignature; id = id)
  io = IOBuffer()
  print(io, bp.f)
  if bp.sig !== nothing
    print(io, '(', join("::" .* string.(bp.sig.types), ", "), ')')
  end
  description = String(take!(io))
  file, line = "", 0
  if bp.f isa Method
    _, _, file, line = Base.arg_decl_parts(bp.f)
    file = Atom.basepath(string(file))
    description = replace(description, r" in .* at .*$" => "")
  end

  condition = bp.condition === nothing ? nothing : get(_conditions, bp, "?")

  return Dict(
    :file => file,
    :line => line,
    :description => description,
    :isactive => bp.enabled[],
    :condition => condition,
    :id => id
  )
end

handle("toggleBP") do item
  with_error_message() do
    if haskey(item, "id") && item["id"] ≠ nothing
      id = item["id"]
      if haskey(_breakpoints, id)
        JuliaInterpreter.remove(_breakpoints[id])
      else
        error("Inconsistent internal state.")
      end
    else
      error("We have a problem, Houston.")
    end
    allbreakpoints()
  end
end

handle("addConditionById") do item, cond
  with_error_message() do
    if haskey(item, "id") && item["id"] ≠ nothing
      id = item["id"]
      if haskey(_breakpoints, id)
        bp = _breakpoints[id]
        bp isa JuliaInterpreter.BreakpointSignature || error("Wrong breakpoint type encountered.")

        func = bp.f
        sig = bp.sig
        line = bp.line

        expr = Meta.parse(cond)
        if !(expr isa JuliaInterpreter.Condition)
          error("Breakpoint condition must be an expression or a tuple of a module and an expression.")
        end
        JuliaInterpreter.remove(bp)
        bp = JuliaInterpreter.breakpoint(func, sig, line, expr)
        _conditions[bp] = cond
      else
        error("Inconsistent Internal State. Abort abort aboooort!")
      end
    else
      error("Inconsistent Internal State. Abort abort aboooort!")
    end

    allbreakpoints()
  end
end

handle("getBreakpoints") do
  with_error_message() do
    Dict(
      :breakpoints => allbreakpoints(),
      :onException => JuliaInterpreter.break_on_error[],
      :onUncaught => false
    )
  end
end

isCompileMode() = _compiledMode[] === JuliaInterpreter.Compiled()

handle("toggleCompiled") do
  with_error_message() do
    _compiledMode[] = isCompileMode() ? JuliaInterpreter.finish_and_return! : JuliaInterpreter.Compiled()
    return isCompileMode()
  end
end

handle("toggleException") do
  with_error_message() do
    if JuliaInterpreter.break_on_throw[]
      JuliaInterpreter.break_off(:throw)
    else
      JuliaInterpreter.break_on(:throw)
    end
    return JuliaInterpreter.break_on_throw[]
  end
end

handle("toggleUncaught") do
  with_error_message() do
    if JuliaInterpreter.break_on_error[]
      JuliaInterpreter.break_off(:error)
    else
      JuliaInterpreter.break_on(:error)
    end
    return JuliaInterpreter.break_on_error[]
  end
end

handle("clearbps") do
  with_error_message() do
    JuliaInterpreter.remove()
    allbreakpoints()
  end
end

handle("getbps") do
  with_error_message() do
    ret = []
    for (k, bp) in bps
      push!(ret, Dict(:view => render(Juno.Inline(), bp)))
    end
    ret
  end
end

handle("addArgs") do arg
  with_error_message() do
    bp = add_breakpoint_args(arg)
    bp isa Vector && isempty(bp) && error("""
      Invalid spec or no matching methods found. Please specify as `foo` or `foo(Bar, Baz)`, e.g.
      `sin` or `sin(Int)`. Make sure the function and all types can be reached from `Main`.
      """)
    allbreakpoints()
  end
end

handle("toggleAllActiveBP") do state
  with_error_message() do
    state ? JuliaInterpreter.disable() : JuliaInterpreter.enable()
    allbreakpoints()
  end
end

handle("toggleActiveBP") do item
  with_error_message() do
    if haskey(item, "id") && item["id"] ≠ nothing
      id = item["id"]
      if haskey(_breakpoints, id)
        bp = _breakpoints[id]
        bp.enabled[] ? JuliaInterpreter.disable(bp) : JuliaInterpreter.enable(bp)
      else
        error("Inconsistent internal state.")
      end
    else
      error("I'm broken.")
    end
    allbreakpoints()
  end
end

function with_error_message(f)
  ret, err = nothing, false
  try
    ret = f()
  catch err
    err = sprint(showerror, err)
  end
  Dict(
    :response => ret,
    :error => err
  )
end

"""
    add_breakpoint_args(arg)

Takes a string of the form `foo` or `foo(Bar, Baz)` and sets a breakpoint for the appropriate methods.
"""
function add_breakpoint_args(arg)
  m = match(r"(.*?)(\(.*\))?$", arg)
  m === nothing && return
  if m[1] ≠ nothing
    if m[2] ≠ nothing
      fun = CodeTools.getthing(Main, m[1])
      args = Main.eval(Meta.parse("tuple$(m[2])"))
      meth = which(fun, args)
    else
      meth = CodeTools.getthing(Main, arg)
    end
  end
  JuliaInterpreter.breakpoint(meth)
end

function no_chance_of_breaking()
  bps = JuliaInterpreter.breakpoints()
  !JuliaInterpreter.break_on_error[] && (isempty(bps) || all(bp -> !bp[].isactive, bps))
end


text:
module ITensorChainRules

using ITensors.NDTensors
using ITensors: Indices

using ChainRulesCore
using ..ITensors

import ChainRulesCore: rrule

ITensors.dag(z::AbstractZero) = z

broadcast_notangent(a) = broadcast(_ -> NoTangent(), a)

include(joinpath("NDTensors", "tensor.jl"))
include(joinpath("NDTensors", "dense.jl"))
include("indexset.jl")
include("itensor.jl")
include(joinpath("physics", "sitetype.jl"))
include(joinpath("LazyApply", "LazyApply.jl"))
include(joinpath("mps", "abstractmps.jl"))
include(joinpath("mps", "mpo.jl"))
include(joinpath("Ops", "Ops.jl"))
include("zygoterules.jl")

@non_differentiable broadcast_notangent(::Any)
@non_differentiable Index(::Any...)
@non_differentiable delta(::Any...)
@non_differentiable dag(::Index)
@non_differentiable inds(::Any...)
@non_differentiable commoninds(::Any...)
@non_differentiable noncommoninds(::Any...)
@non_differentiable uniqueinds(::Any...)
@non_differentiable SiteType(::Any)
@non_differentiable ITensors._sitetypes(::Any)
@non_differentiable addtags(::TagSet, ::Any)
@non_differentiable ITensors.filter_inds_set_function(::Function, ::Function, ::Any...)
@non_differentiable ITensors.filter_inds_set_function(::Function, ::Any...)

end

ext: jl


text:
##### Beginning of file

import Pkg # stdlib

struct _TomlFile
    filename::String
    function _TomlFile(path::String)::_TomlFile
        path::String = abspath(strip(path))
        if isfile(path)
            result::_TomlFile = new(path)
            return result
        else
            error("File does not exist")
        end
    end
end

function _parse_toml_file(x::_TomlFile)::Dict{String, Any}
    toml_file_filename::String = x.filename
    toml_file_text::String = read(toml_file_filename, String)
    toml_file_parsed::Dict{String, Any} = Pkg.TOML.parse(toml_file_text)
    return toml_file_parsed
end

function _version_string(x::_TomlFile)::String
    toml_file_parsed::Dict{String, Any} = _parse_toml_file(x)
    version_string::String = toml_file_parsed["version"]
    return version_string
end

function _version_string()::String
    Snapshots_toml_file::_TomlFile = _TomlFile(
        package_directory("Project.toml")
        )
    result_versionstring::String = _version_string(Snapshots_toml_file)
    return result_versionstring
end

function _version_string(m::Method)::String
    m_package_directory::String = package_directory(m)
    m_toml_file::_TomlFile = _TomlFile(
        joinpath(m_package_directory, "Project.toml")
        )
    result_versionstring::String = _version_string(m_toml_file)
    return result_versionstring
end

function _version_string(f::Function)::String
    m_package_directory::String = package_directory(f)
    m_toml_file::_TomlFile = _TomlFile(
        joinpath(m_package_directory, "Project.toml")
        )
    result_versionstring::String = _version_string(m_toml_file)
    return result_versionstring
end

function _version_string(f::Function, types::Tuple)::String
    m_package_directory::String = package_directory(f, types)
    m_toml_file::_TomlFile = _TomlFile(
        joinpath(m_package_directory, "Project.toml")
        )
    result_versionstring::String = _version_string(m_toml_file)
    return result_versionstring
end

function _version_string(m::Module)::String
    m_package_directory::String = package_directory(m)
    m_toml_file::_TomlFile = _TomlFile(
        joinpath(m_package_directory, "Project.toml")
        )
    result_versionstring::String = _version_string(m_toml_file)
    return result_versionstring
end

"""
    version()::VersionNumber

Return the version number of Snapshots.
"""
function version()::VersionNumber
    result_versionstring::String = _version_string()
    result_versionnumber::VersionNumber = VersionNumber(result_versionstring)
    return result_versionnumber
end

"""
    version(m::Method)::VersionNumber

If method `m`
is part of a Julia package, returns the version number of that package.

If method `m`
is not part of a Julia package, throws an error.
"""
function version(m::Method)::VersionNumber
    result_versionstring::String = _version_string(m)
    result_versionnumber::VersionNumber = VersionNumber(result_versionstring)
    return result_versionnumber
end

"""
    version(f::Function)::VersionNumber

If function `f`
is part of a Julia package, returns the version number of
that package.

If function `f`
is not part of a Julia package, throws an error.
"""
function version(f::Function)::VersionNumber
    result_versionstring::String = _version_string(f)
    result_versionnumber::VersionNumber = VersionNumber(result_versionstring)
    return result_versionnumber
end

"""
    version(f::Function, types::Tuple)::VersionNumber

If function `f` with type signature `types`
is part of a Julia package, returns the version number of
that package.

If function `f` with type signature `types`
is not part of a Julia package, throws an error.
"""
function version(f::Function, types::Tuple)::VersionNumber
    result_versionstring::String = _version_string(f, types)
    result_versionnumber::VersionNumber = VersionNumber(result_versionstring)
    return result_versionnumber
end

"""
    version(m::Module)::VersionNumber

If module `m` is part of a Julia package, returns the version number of
that package.

If module `m` is not part of a Julia package, throws an error.
"""
function version(m::Module)::VersionNumber
    result_versionstring::String = _version_string(m)
    result_versionnumber::VersionNumber = VersionNumber(result_versionstring)
    return result_versionnumber
end

##### End of file


text:
#
# types.jl -
#
# Definitions of Tcl constants and types.
#

"""

`TclStatus` is used to represent a status like `TCL_OK` or `TCL_ERROR` returned
by Tcl functions.

"""
struct TclStatus
    code::Cint
end

# Status codes returned by Tcl functions.
const TCL_OK       = TclStatus(0)
const TCL_ERROR    = TclStatus(1)
const TCL_RETURN   = TclStatus(2)
const TCL_BREAK    = TclStatus(3)
const TCL_CONTINUE = TclStatus(4)

Base.show(io::IO, ::MIME"text/plain", x::TclStatus) =
    print(io, (x == TCL_OK       ? "TCL_OK" :
               x == TCL_ERROR    ? "TCL_ERROR" :
               x == TCL_RETURN   ? "TCL_RETURN" :
               x == TCL_BREAK    ? "TCL_BREAK" :
               x == TCL_CONTINUE ? "TCL_CONTINUE" :
               "TclStatus("*string(x.code)*")"))

# Flags for settings the result.
const TCL_VOLATILE = Ptr{Cvoid}(1)
const TCL_STATIC   = Ptr{Cvoid}(0)
const TCL_DYNAMIC  = Ptr{Cvoid}(3)

# Flags for Tcl variables.
const TCL_GLOBAL_ONLY    = Cint(1)
const TCL_NAMESPACE_ONLY = Cint(2)
const TCL_APPEND_VALUE   = Cint(4)
const TCL_LIST_ELEMENT   = Cint(8)
const TCL_LEAVE_ERR_MSG  = Cint(0x200)

# Flags for Tcl processing events.  Set TCL_DONT_WAIT to not sleep: process
# only events that are ready at the time of the call.  Set TCL_ALL_EVENTS to
# process all kinds of events: equivalent to OR-ing together all of the below
# flags or specifying none of them.
const TCL_DONT_WAIT     = Cint(1<<1)
const TCL_WINDOW_EVENTS = Cint(1<<2) # Process window system events.
const TCL_FILE_EVENTS   = Cint(1<<3) # Process file events.
const TCL_TIMER_EVENTS  = Cint(1<<4) # Process timer events.
const TCL_IDLE_EVENTS   = Cint(1<<5) # Process idle callbacks.
const TCL_ALL_EVENTS    = ~TCL_DONT_WAIT      # Process all kinds of events.

# The following values control how blocks are combined into photo images when
# the alpha component of a pixel is not 255, a.k.a. the compositing rule.
const TK_PHOTO_COMPOSITE_OVERLAY = Cint(0)
const TK_PHOTO_COMPOSITE_SET     = Cint(1)

# Flags for evaluating scripts/commands.
const TCL_NO_EVAL       = Cint(0x010000)
const TCL_EVAL_GLOBAL   = Cint(0x020000)
const TCL_EVAL_DIRECT   = Cint(0x040000)
const TCL_EVAL_INVOKE   = Cint(0x080000)
const TCL_CANCEL_UNWIND = Cint(0x100000)
const TCL_EVAL_NOERR    = Cint(0x200000)

struct TclError <: Exception
    msg::String
end

Base.showerror(io::IO, ex::TclError) = print(io, "Tcl/Tk error: ", ex.msg)

# Structure to store a pointer to a Tcl interpreter. (Even though the address
# should not be modified, it is mutable because immutable objects cannot be
# finalized.)
const TclInterpPtr = Ptr{Cvoid}
mutable struct TclInterp
    ptr::TclInterpPtr
    TclInterp(ptr::TclInterpPtr) = new(ptr)
end

"""

A `ManagedObject` is a Tcl or Tk object which manages its reference count and
of which a pointer to a Tcl object can be retrieved by the `__objptr` method.

"""
abstract type ManagedObject end

# Structure to store a pointer to a Tcl object. (Even though the address
# should not be modified, it is mutable because immutable objects cannot be
# finalized.)  The constructor will refuse to build a managed Tcl object with
# a NULL address.
const TclObjPtr = Ptr{Cvoid}
mutable struct TclObj{T} <: ManagedObject
    ptr::TclObjPtr
    function TclObj{T}(ptr::TclObjPtr) where {T}
        ptr != C_NULL || __illegal_null_object_pointer()
        obj = new{T}(Tcl_IncrRefCount(ptr))
        return finalizer(__finalize, obj)
    end
end

struct Callback <: ManagedObject
    intptr::TclInterpPtr # weak reference to interpreter
    obj::TclObj{Function} # command name (possibly fully-qualified)
    func::Function
end

# Tcl wide integer is 64-bit integer.
const WideInt = Int64

# Type used in the signature of a Tcl list object (a.k.a. vector in Julia).
const List = Vector

# Token used by Tcl to identify an object command.
const TclCommand = Ptr{Cvoid}

# Floating-point types.
const FloatingPoint = Union{Irrational,Rational,AbstractFloat}

"""

The abstract type `Trait` is inherited by types indicating specific traits.

See also: [`Tcl.AtomicType`](@ref).

"""
abstract type Trait end

# Defaint the `AtomicType` trait and its 2 singleton sub-types `Atomic` and
# `NonAtomic`.
abstract type AtomicType <: Trait end
for T in (:NonAtomic, :Atomic)
    @eval begin
        struct $T <: AtomicType end
        @doc @doc(AtomicType) $T
    end
end

# Client data used by commands and callbacks.
const ClientData = Ptr{Cvoid}

const TclObjList    = TclObj{List}
const TclObjCommand = TclObj{Function}

# `Name` is anything that can be understood as the name of a variable or of a
# command.
const Name = Union{AbstractString,Symbol,TclObj{String}}

# `StringOrSymbol` can be automatically converted into a `Cstring` by `ccall`.
const StringOrSymbol = Union{AbstractString,Symbol}

# A `Byte` is any bits type that is exactly 8 bits.
const Byte = Union{UInt8,Int8}

# Objects of type `Iterables` are considered as iterators, making an object out
# of them yield a Tcl list.
const Iterables = Union{AbstractVector,Tuple,Set,BitSet}

#------------------------------------------------------------------------------
# Tk widgets and other Tk objects.

abstract type TkObject     <: ManagedObject end
abstract type TkWidget     <: TkObject      end
abstract type TkRootWidget <: TkWidget      end

# An image is parameterized by the image type (capitalized).
mutable struct TkImage{T} <: TkObject
    interp::TclInterp
    path::String
end

# We want to have the object type and path both printed in the REPL but want
# only the object path with the `string` method or for string interpolation.
# Note that: "$w" calls `string(w)` while "anything $w" calls `show(io, w)`.

Base.show(io::IO, ::MIME"text/plain", w::T) where {T<:TkObject} =
    print(io, "$T(\"$(string(w))\")")

Base.show(io::IO, w::TkObject) = print(io, getpath(w))

#------------------------------------------------------------------------------
# Colors

abstract type TkColor end

struct TkGray{T} <: TkColor
    gray::T
end

struct TkRGB{T} <: TkColor
    r::T; g::T; b::T
end

struct TkBGR{T} <: TkColor
    b::T; g::T; r::T
end

struct TkRGBA{T} <: TkColor
    r::T; g::T; b::T; a::T
end

struct TkBGRA{T} <: TkColor
    b::T; g::T; r::T; a::T
end

struct TkARGB{T} <: TkColor
    a::T; r::T; g::T; b::T
end

struct TkABGR{T} <: TkColor
    a::T; b::T; g::T; r::T
end

const TkColorsWithAlpha{T} = Union{TkRGBA{T},TkBGRA{T},TkARGB{T},TkABGR{T}}
const TkColors{T} = Union{TkRGB{T},TkBGR{T},TkColorsWithAlpha{T}}

gray(c::TkGray{T}) where T = c.gray
red(c::TkColors{T}) where T = c.r
green(c::TkColors{T}) where T = c.g
blue(c::TkColors{T}) where T = c.b
alpha(c::TkColorsWithAlpha{T}) where T = c.a

Base.show(io::IO, ::MIME"text/plain", c::TkGray{T}) where {T} =
    print(io,"TkGray{",T,"}(",gray(c),")")

Base.show(io::IO, ::MIME"text/plain", c::TkRGB{T}) where {T} =
    print(io,"TkRGB{",T,"}(",red(c),",",green(c),",",blue(c),")")

Base.show(io::IO, ::MIME"text/plain", c::TkBGR{T}) where {T} =
    print(io,"TkRGB{",T,"}(",blue(c),",",green(c),",",red(c),")")

Base.show(io::IO, ::MIME"text/plain", c::TkRGBA{T}) where {T} =
    print(io,"TkRGBA{",T,"}(",red(c),",",green(c),",",blue(c),",",alpha(c),")")

Base.show(io::IO, ::MIME"text/plain", c::TkBGRA{T}) where {T} =
    print(io,"TkRGBA{",T,"}(",blue(c),",",green(c),",",red(c),",",alpha(c),")")

Base.show(io::IO, ::MIME"text/plain", c::TkARGB{T}) where {T} =
    print(io,"TkARGB{",T,"}(",alpha(c),",",red(c),",",green(c),",",blue(c),")")

Base.show(io::IO, ::MIME"text/plain", c::TkABGR{T}) where {T} =
    print(io,"TkARGB{",T,"}(",alpha(c),",",blue(c),",",green(c),",",red(c),")")

# Extend Base.show for passing colors to Tk (the alpha component, if any, is
# ignored).
Base.show(io::IO, c::TkGray{T}) where {T<:Union{UInt8,UInt16}} =
    (s = _hex(gray(c)); print(io, "#", s, s, s))
Base.show(io::IO, c::TkColors{T}) where {T<:Union{UInt8,UInt16}} =
    print(io, "#", _hex(red(c)), _hex(green(c)), _hex(blue(c)))
Base.show(io::IO, c::TkColors{T}) where {T<:Union{UInt32,UInt64}} =
    print(io, "#", _hex16(red(c)), _hex16(green(c)), _hex16(blue(c)))

_hex(x::UInt8)  = string(x; base=16, pad=2)
_hex(x::UInt16) = string(x; base=16, pad=4)
_hex(x::UInt32) = string(x; base=16, pad=8)
_hex(x::UInt64) = string(x; base=16, pad=16)
_hex16(x::UInt32) = _hex((x >> 16)%UInt16)
_hex16(x::UInt64) = _hex((x >> 48)%UInt16)


text:
"The perfect sidekick to your scientific inquiries"
module DrWatson
import Pkg, LibGit2
const PATH_SEPARATOR = joinpath("_", "_")[2]

# Misc functions for kw-macros
convert_to_kw(ex::Expr) = Expr(:kw,ex.args...)
convert_to_kw(ex) = error("invalid keyword argument syntax \"$ex\"")

# Misc helpers
"""
    readenv(var, default::T) where {T}

Try to read the environment variable `var` and parse it as a `::T`.
If that fails, return `default`.
"""
readenv(var, default::T) where {T} = something(tryparse(T, get(ENV, var, "")), Some(default))

# Pure Julia implementation
include("project_setup.jl")
include("naming.jl")
include("saving_tools.jl")

using UnPack
export @pack!, @unpack

# Functionality that saves/loads
using FileIO
using JLD2
export save, load
export wsave, wload

_wsave(filename, data...; kwargs...) = FileIO.save(filename, data...; kwargs...)

"""
    wsave(filename, data...; kwargs...)

Save `data` at `filename` by first creating the appropriate paths.
Default fallback is `FileIO.save`. Extend this for your types
by extending `DrWatson._wsave(filename, data::YourType...; kwargs...)`.
"""
function wsave(filename, data...; kwargs...)
    mkpath(dirname(filename))
    return _wsave(filename, data...; kwargs...)
end

"Currently equivalent with `FileIO.load`."
wload(data...; kwargs...) = FileIO.load(data...; kwargs...)

include("saving_files.jl")
include("dict_list.jl")

# Functionality that requires Dataframes and other heavy dependencies:
using Requires

# Update messages
using Scratch
const display_update = true
const update_version = "2.7.3"
const update_name = "update_v$update_version"

# Get scratch space for this package
function __init__()
    @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        include("result_collection.jl")
    end
    if display_update
        versions_dir = @get_scratch!("versions")

        if !isfile(joinpath(versions_dir, update_name))

        printstyled(stdout,
        """
        \nUpdate message: DrWatson v$update_version

        * New section "Taking project input-output automation to 11" in the documentation.
          It showcases how to eliminate code duplication and streamline your simulation setup
          and run phase using `savename` and `produce_or_load`.
        * By default now `gitpatch` is NOT saved when calling `tag!` and derivative functions.
          This is due to an unknown problem that causes collecting the git patch to
          never halt, potentially not saving a user's output.
        \n
        """; color = :light_magenta)
        touch(joinpath(versions_dir, update_name))
        end
    end
end

end # Module


text:
module TestLSS

using QuantEcon
using Base.Test
using FactCheck

rough_kwargs = Dict(:atol => 1e-7, :rtol => 1e-7)

# set up
A = .95
C = .05
G = 1.
mu_0 = [.75;]
Sigma_0 = fill(0.000001, 1, 1)

ss = LSS(A, C, G, mu_0)
ss1 = LSS(A, C, G, mu_0, Sigma_0)

vals = stationary_distributions(ss, max_iter=1000, tol=1e-9)

facts("Testing lss.jl") do
    context("test stationarity") do
        vals = stationary_distributions(ss, max_iter=1000, tol=1e-9)
        ssmux, ssmuy, sssigx, sssigy = vals

        @fact ssmux --> roughly(ssmuy; rough_kwargs...)
        @fact sssigx --> roughly(sssigy; rough_kwargs...)
        @fact ssmux --> roughly([0.0]; rough_kwargs...)
        @fact sssigx --> roughly(ss.C.^2 ./ (1 - ss.A .^2); rough_kwargs...)
    end

    context("test replicate") do
        xval1, yval1 = replicate(ss, 100, 5000)
        xval2, yval2 = replicate(ss; t=100, num_reps=5000)
        xval3, yval3 = replicate(ss1; t=100, num_reps=5000)

        for (x, y)  in [(xval1, yval1), (xval2, yval2), (xval3, yval3)]
            @fact x --> roughly(y; rough_kwargs...)
            @fact abs(mean(x)) <= 0.05 --> true
        end

    end

    context("test convergence error") do
        @fact_throws stationary_distributions(ss; max_iter=1, tol=eps())
    end

    context("test geometric_sums") do
        β = 0.98
        xs = rand(10)
        for x in xs
            gsum_x, gsum_y = QuantEcon.geometric_sums(ss, β, [x])
            @fact gsum_x --> roughly([x/(1-β *ss.A[1])])
            @fact gsum_y --> roughly([ss.G[1]*x/(1-β *ss.A[1])])
        end
    end

    context("test constructors") do
        # kwarg version
        other_ss = LSS(A, C, G; mu_0=[mu_0;])
        for nm in fieldnames(ss)
            @fact getfield(ss, nm) --> getfield(other_ss, nm)
        end
    end

end  # facts
end  # module


text:
#=
Compare a 2D E/I linear model to the corresponding Hawkes process
=#

using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:dark) ; plotlyjs()
using FFTW

using ProgressMeter
using Random
Random.seed!(0)

using HawkesSimulator; const global H = HawkesSimulator

function onedmat(x::Real)
  return fill(x,(1,1))
end

## rate model first

function iofunction(v::R,_v0::R,_α::R) where R<:Real
  return max(_α*(v-_v0),zero(R))
end
function iofunction_inv(r::R,_v0::R,_α::R) where R<:Real
  return max(zero(R),r/_α+_v0)
end

global const v0 = -70.0 # -70
global const α = 0.03 # 0.3
global const taus = [20.0, 10.0].*1E-3
global const weights = [1.25 -0.65 ; 1.2 -0.5]
global const input = [50.0, 50.0]
global const dt = 0.01E-3
global const Ttot = 1.0
global const taus_hawk = [3.0,3.0]

##

_ = let vs = range(-90,50.,length=150)
  rs = @. iofunction(vs,v0,α)
  plot(vs,rs;linewidth=2,leg=false)
end

##

function euler_dv(v::Vector{R},_weights::Matrix{R},_input::Vector{R},_taus::Vector{R},_v0::R,_α::R) where R
  r = @. iofunction(v,_v0,_α)
  return (.-v .+ (_weights*r) .+ _input) ./ _taus
end

function run_2D_network()
  v_start_e,v_start_i = 0.0,0.0
  times = range(0.0,Ttot;step=dt)
  ntimes = length(times)
  rates = Matrix{Float64}(undef,2,ntimes)
  v = [v_start_e,v_start_i]
  for tt in eachindex(times)
    rates[:,tt] = iofunction.(v,v0,α)
    dv = euler_dv(v,weights,input,taus,v0,α)
    v += dv*dt
  end
  return times,rates
end

##

runtimes,runrates = run_2D_network()

## Now the Hawkes process
weights_equiv = α^2 .* weights
inputs_equiv =  α .* (input .- (α .* weights+I)*fill(v0,2))

pop_e,pop_i = H.PopulationExp(taus_hawk[1]),H.PopulationExp(taus_hawk[2])

ps_e,ps_i = let _inp = inputs_equiv
  H.PopulationState(pop_e,[_inp[1]]),
  H.PopulationState(pop_i,[_inp[2]])
end

inntw_e,inntw_i = let _weights = weights_equiv
   H.InputNetwork(ps_e,[ps_e, ps_i],[onedmat(_weights[1,1]),onedmat(_weights[1,2])]),
   H.InputNetwork(ps_i,[ps_e, ps_i],[onedmat(_weights[2,1]),onedmat(_weights[2,2])])
end

function run_simulation!(netws,nspikes)
  t_now = 0.0
  H.reset!.(netws) # clear spike trains etc
  @showprogress 1.0 "Running Hawkes process..." for k in 1:nspikes
    t_now = H.dynamics_step!(t_now,netws)
    if k%1_000 == 0
      for netw in netws
        H.clear_trains!(netw.postpops)
      end
    end
  end
  for netw in netws
    H.clear_trains!(netw.postpops,-1.0)
  end
  return t_now
end
##
nspikes = 50_000
t_end=run_simulation!([inntw_e,inntw_i],nspikes)

##

rates_hawkes = [ H.numerical_rates(ps_e)[1],H.numerical_rates(ps_i)[1]]

@info "final rates are $(runrates[:,end])"
@info "Hawkes final rates are $rates_hawkes"
##


text:
module NumFields

export embed_field, reduce!, NumField

using StaticArrays
using StructEquality

abstract type NumField <: Number end
function embed_field end
function reduce! end
macro simple_number_field_lazy(name, field_gen_coeffs, generator)
    p = eval(field_gen_coeffs)
    T = eltype(p)
    N = length(p)

    A = zeros(T, N, N)
    for i = 1:N-1
        A[i+1,i] = 1
    end
    A[:, N] = p
    A_gen = SMatrix{N, N, T}(A)
    powers = [(A_gen^i)[:, 1] for i=0:2*N-2]

    return quote
        mutable struct $name <: NumField
            coeffs :: MVector{$N, $(esc(T))}
            denom :: $(esc(T))
            reduced :: Bool
        end

        Base.promote_rule(::Type{$(esc(name))}, ::Type{$(esc(T))}) = $(esc(name))

        function ($(esc(name)))(n :: Integer)
            n_converted = convert($(esc(T)), n)
            return $(esc(name))([i==1 ? n_converted : 0 for i=1:$N], 1, true)
        end

        function ($(esc(name)))(v :: Vector{<:Integer}, denom :: Integer)
            return $(esc(name))(v, denom, false)
        end

        function ($(esc(name)))(v :: Vector{<:Integer})
            return $(esc(name))(v, 1)
        end

        function ($(esc(name)))(n :: Rational{<:Integer})
            x = ($(esc(name)))(numerator(n))
            x.denom = denominator(n)
            return x
        end

        function NumFields.reduce!(x :: $(esc(name)))
            if !x.reduced
                divisor = gcd(x.coeffs..., x.denom)
                x.coeffs = [div(c,divisor) for c in x.coeffs]
                x.denom = div(x.denom, divisor)
                x.reduced = true
            end
            return x
        end

        const powers = $powers
        const $(esc(generator)) = $(esc(name))([i==2 ? 1 : 0 for i=1:$N],1,true)

        function Base.:+(x :: $(esc(name)), y :: $(esc(name)))
            return $(esc(name))(x.coeffs*y.denom + y.coeffs*x.denom, x.denom*y.denom, false)
        end

        function Base.:-(x :: $(esc(name)), y :: $(esc(name)))
            return $(esc(name))(x.coeffs*y.denom - y.coeffs*x.denom, x.denom*y.denom, false)
        end

        function Base.:-(x :: $(esc(name)))
            return $(esc(name))(-x.coeffs, x.denom, x.reduced)
        end

        function Base.:*(x :: $(esc(name)), y :: $(esc(name)))
            coeffs = zeros(MVector{$N, $(esc(T))})

            @inbounds for i = 0:(2*$N-2)
                factor = 0
                for j = max(0,i-$N+1):min($N-1,i)
                    factor += x.coeffs[j+1] * y.coeffs[i-j+1]
                end
                coeffs += powers[i+1] * factor
            end

            return $(esc(name))(coeffs, x.denom*y.denom, false)
        end

        function Base.:*(x :: $(esc(name)), y :: Number)
            $(esc(name))(x.coeffs*numerator(y),x.denom*denominator(y), false)
        end

        function Base.:*(x :: Number, y :: $(esc(name)))
            $(esc(name))(y.coeffs*numerator(x),y.denom*denominator(x), false)
        end

        function Base.:(==)(x :: $(esc(name)), y :: $(esc(name)))
            x1 = reduce!(x)
            y1 = reduce!(y)
            return x1.coeffs == y1.coeffs && x1.denom == y1.denom
        end

        function Base.://(x :: $(esc(name)), y :: Number)
            return $(esc(name))(x.coeffs*denominator(y),x.denom*numerator(y), false)
        end
        function Base.hash(x :: $(esc(name)))
            x1 = reduce!(x)
            return hash((x1.coeffs, x1.denom))
        end

        function Base.:/(x :: $(esc(name)), y :: Integer)
            return $(esc(name))(x.coeffs,x.denom*y, false)
        end

        function NumFields.embed_field(map_coeff, map_gen, x :: $(esc(name)))
            result = zero(typeof(map_gen))
            for i=1:$N
                result += map_coeff(x.coeffs[i])*map_gen^(i-1)
            end
            return result
        end
    end
end

macro simple_number_field_concrete(name, field_gen_coeffs, generator)
    p = eval(field_gen_coeffs)
    T = eltype(p)
    N = length(p)

    A = zeros(T, N, N)
    for i = 1:N-1
        A[i+1,i] = 1
    end
    A[:, N] = p
    A_gen = SMatrix{N, N, T}(A)
    powers = [(A_gen^i)[:, 1] for i=0:2*N-2]

    return quote
        struct $name <: NumField
            coeffs :: SVector{$N, $(esc(T))}
            denom :: $(esc(T))

            function ($(esc(name)))(v :: AbstractVector{<:Integer}, denom :: Integer)
                divisor = gcd(v..., denom)
                return new(div.(v,divisor), div(denom, divisor))
            end
        end

        # Basic conversion
        Base.promote_rule(::Type{$(esc(name))}, ::Type{$(esc(T))}) = $(esc(name))

        function ($(esc(name)))(n :: Integer)
            n_converted = convert($(esc(T)), n)
            return $(esc(name))([i==1 ? n_converted : 0 for i=1:$N], 1)
        end
        function ($(esc(name)))(v :: AbstractVector{<:Integer})
            return $(esc(name))(v, 1)
        end
        function ($(esc(name)))(n :: Rational{<:Integer})
            n_converted = convert($(esc(T)), numerator(n))
            return $(esc(name))([i==1 ? n_converted : 0 for i=1:$N], denominator(n))
        end

        # Addition, substraction and multiplication
        const powers = $powers
        const $(esc(generator)) = $(esc(name))([i==2 ? 1 : 0 for i=1:$N],1)

        function Base.:+(x :: $(esc(name)), y :: $(esc(name)))
            return ($(esc(name))(x.coeffs*y.denom + y.coeffs*x.denom, x.denom*y.denom))
        end
        function Base.:-(x :: $(esc(name)), y :: $(esc(name)))
            return ($(esc(name))(x.coeffs*y.denom - y.coeffs*x.denom, x.denom*y.denom))
        end
        function Base.:-(x :: $(esc(name)))
            return ($(esc(name))(-x.coeffs, x.denom))
        end
        function Base.:*(x :: $(esc(name)), y :: $(esc(name)))
            coeffs = zeros(SVector{$N, $(esc(T))})

            @inbounds for i = 0:(2*$N-2)
                factor = 0
                for j = max(0,i-$N+1):min($N-1,i)
                    factor += x.coeffs[j+1] * y.coeffs[i-j+1]
                end
                coeffs += powers[i+1] * factor
            end

            return ($(esc(name))(coeffs, x.denom*y.denom))
        end
        function Base.:*(x :: $(esc(name)), y :: Number)
            return ($(esc(name))(x.coeffs*numerator(y),x.denom*denominator(y)))
        end
        function Base.:*(x :: Number, y :: $(esc(name)))
            return ($(esc(name))(y.coeffs*numerator(x),y.denom*denominator(x)))
        end

        # Equality and hashes
        function Base.:(==)(x :: $(esc(name)), y :: $(esc(name)))
            return x.coeffs == y.coeffs && x.denom == y.denom
        end
        function Base.isequal(x :: $(esc(name)), y :: $(esc(name)))
            return x == y
        end
        function Base.hash(x :: $(esc(name)))
            return hash((x.coeffs, x.denom))
        end

        # Division
        function Base.://(x :: $(esc(name)), y :: Number)
            return ($(esc(name))(x.coeffs*denominator(y),x.denom*numerator(y)))
        end
        function Base.:/(x :: $(esc(name)), y :: Number)
            return ($(esc(name))(x.coeffs*denominator(y),x.denom*numerator(y)))
        end

        # Universal property
        function NumFields.embed_field(map_coeff, map_gen, x :: $(esc(name)))
            result = zero(typeof(map_gen))
            for i=1:$N
                result += map_coeff(x.coeffs[i])*map_gen^(i-1)
            end
            return result/map_coeff(x.denom)
        end

    end
end

end


text:
using ProtoBuf

tempdir = mktempdir()
example_path = "tensorflow/core/example"
exampledir = "$tempdir/$example_path"
mkpath(exampledir)
tensorflow_repo_url="https://raw.githubusercontent.com/tensorflow/tensorflow"
tensorflow_tag="v1.15.0"
tensorflow_example_url="$tensorflow_repo_url/$tensorflow_tag/$example_path"
download("$tensorflow_example_url/feature.proto", "$exampledir/feature.proto")
download("$tensorflow_example_url/example.proto", "$exampledir/example.proto")
run(ProtoBuf.protoc(`-I=$tempdir --julia_out=../src/generated $exampledir/example.proto`))


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

JLLWrappers.@generate_main_file_header("LCIO")
JLLWrappers.@generate_main_file("LCIO", UUID("35fd387a-780d-5580-bc3c-2dd06d6487c2"))
end  # module LCIO_jll


text:
# ------------------------------------------------------------------------------
# Geometric objects
# -------------------------------------------------------------------------------

"""
    squaredist(A, B)

Compute the square of the euclidean distance between `A` and `B`.

# Examples
```jldoctest
julia> squaredist(Point(0, 0), Point(2,2))
8
```
"""
function squaredist(A, B)
    d = (A.x-B.x)^2+(A.y-B.y)^2
    simplify(Sym(d))
end

"""
    distance(A, B)

Compute the euclidean distance between `A` and `B`, which is the "ordinary" straight-line distance
between two points.

# Examples
```jldoctest
julia> distance(Point(0, 0), Point(3,4))
5
```
"""
function distance(A, B)
    SymPy.sqrt(squaredist(A, B))
end

"""
    midpoint(A, B)

Find the midpoint between `A` and `B`.

# Examples

```jldoctest
julia> midpoint(Point(0, 0), Point(2,4))
Point(1, 2)
```
"""
function midpoint(A, B)
    x1 = (A.x + B.x)/2
    y1 = (A.y + B.y)/2
    Point(x1, y1)
end

"""
    concurrent(edgelist::Vector{Edge})

Find the point where all edges (lines) in `edgelist` intersect if such point exists. Return nothing otherwise.

# Examples

```jldoctest
julia> concurrent([Edge(Point(0,1), Point(1,1)), Edge(Point(0,1), Point(1,0))])
Point(0, 1)
```
"""
function concurrent(edgelist::Vector{Edge})
    enum = length(edgelist)
    if enum == 1
        throw(ArgumentError("At least two edges are needed. $enum is given."))
    end

    @vars x y

    eqs = Sym[]
    for e in edgelist
        # This should be 0
        eq = (x-e.src.x)*(y-e.dst.y) - (x-e.dst.x)*(y-e.src.y)
        push!(eqs, eq)
    end

    sol = solve(eqs, [x,y])

    if length(sol) == 0
        return nothing
    else
        return Point(simplify(sol[x]), simplify(sol[y]))
    end
end

"""
    concurrent(e1::Edge, e2::Edge)

Find the point the edges `e1` and `e2` intersect if such point exists. Return nothing otherwise.

# Examples

```jldoctest
julia> concurrent(Edge(Point(0,1), Point(1,1)), Edge(Point(0,1), Point(1,0)))
Point(0, 1)
```
"""
function concurrent(e1::Edge, e2::Edge)
    return concurrent([e1, e2])
end


text:
"""
    BCDiag

A structure of diagnostic properties of a binary Classifier.
Facilitates summary plots and tables.
"""
struct BCDiag
    n::Int                  # Number of observations
    n1::Int                 # Number of Class 1
    n0::Int                 # Number of Class 0
    baserate::Float64       # Class 1 incidence rate
    ks::Float64             # Maximum separation
    ksarg::Float64          # max sep occurs at this predictor value
    ksdep::Float64          # max sep occurs at this depth

    concordant::Int         # concordant pairs
    tied::Int               # tied pairs
    discordant::Int         # discordant pairs
    auc::Float64            # Area under ROC = (concordant + 0.5tied)/ total pairs
    gini::Float64           # Gini = 2C - 1

    grp::Vector             # group number [0-groups]
    depth::Vector           # Depth [0, 1]
    cdf1::Vector            # cum1 / n1 = True Positive Rate
    cdf0::Vector            # cum0 / n0 = False Positive Rate

    cntg::Vector            # count this group
    cnt1::Vector            # count class 1 this group
    prd1::Vector            # predicted class 1 this group
    rrObs::Vector           # response rate actual
    rrPrd::Vector           # response rate Predicted

    cumg::Vector            # cum count up to this group
    cum1::Vector            # cum count class 1 up to this group
    cpr1::Vector            # cum predicted class 1 up to this group
    crObs::Vector           # cum response rate actual
    crPrd::Vector           # cum response rate predicted
end

function Base.show(io::IO, ::MIME"text/plain", x::BCDiag)
    print(io,
        "Base rate: ", round(x.baserate, digits = 4),"   n: ", x.n, "   n1: ", x.n1, "   n0: ", x.n0, "\n",
        "ks:        ", round(x.ks, digits = 4), "   occurs at value of ", x.ksarg, " depth of ", x.ksdep, "\n",
        "auroc:     ", round(x.auc, digits = 4), "   concordant: ", x.concordant,
        "   tied: ", x.tied, "   discordant: ", x.discordant, "\n",
        "Gini:      ", round(x.gini, digits = 4)
    )
end

function Base.show(io::IO, x::BCDiag)
    print(io, round(x.ks, digits = 4), "   ", round(x.auc, digits = 4) )
end

"""
    bcdiag(target, pred; groups = 100, rev = true, tie = 1e-6)

Perform diagnostics of a binary classifier.
`target` is a 2 level categorical variable, `pred` is probability of class 1.
`groups` is the number of bins to use for plotting/printing.
`rev` = true orders `pred` from high to low.
`tie` is the tolerance of `pred` where values are considered tied.

Returns a BCDiag struct which can be used for plotting or printing:
- `biasplot` is calibration plot of `target` response rate vs. `pred` response rate
- `ksplot` produces ksplot of cumulative distributions
- `rocplot` plots the Receiver Operating Characteristics curve
- `accuracyplot` plots the accuracy curve with adjustable utility
- `liftcurve` is the lift curve
- `cumliftsurve` is the cumulative lift surve
- `liftable` is the lift table as a DataFrame
- `cumliftable` is the cumulative lift table as a DataFrame
"""
function bcdiag(target::BitArray{1}, pred::Vector; groups = 100, rev = true, tie = 1e-6)
    ks = kstest(target, pred; rev = rev)
    roc = auroc(target, pred; tie = tie)

    grpid = ranks(pred, groups = groups, rev = rev)
    frqtbl = freqtable(grpid, target)
    frq = frqtbl.array
    grp = NamedArrays.names(frqtbl, 1)

    n0, n1 = sum(frq, dims = 1)
    n = n0 + n1

    cntg = vec(sum(frq, dims = 2))          # bin count
    cnt1 = @view frq[:, 2]                  # bin count class 1
    cnt0 = @view frq[:, 1]                  # bin count class 0
    cumg = cumsum(cntg)                     # cum bin count
    cum1 = cumsum(cnt1)                     # cum class 1
    cum0 = cumsum(cnt0)                     # cum class 0

    depth = cumg ./ n                       # Depth including this group
    cdf1  = cum1 ./ n1                      # TPR: True Positive Rate
    cdf0  = cum0 ./ n0                      # FPR: False Positive Rate

    rrObs = cnt1 ./ cntg                    # actual response rate
    crObs = cum1 ./ cumg                    # cum actual response rate

    prd1 = Vector{Float64}(undef, length(grp))  # predicted class 1
    for (i, g) in enumerate(grp)
        prd1[i] = sum(pred[grpid .== g])
    end
    cpr1 = cumsum(prd1)                     # cum predicted class 1

    rrPrd = prd1 ./ cntg                    # predicted response rate
    crPrd = cpr1 ./ cumg                    # cum predicted response rate

    BCDiag( ks..., roc..., grp, depth, cdf1, cdf0,
            cntg, cnt1, prd1, rrObs, rrPrd, cumg, cum1, cpr1, crObs, crPrd )
end

function bcdiag(target::Vector, pred::Vector; groups = 100, rev = true, tie = 1e-6)
    uc = sort(unique(target))
    length(uc) == 2 || error(ArgumentError("target should have 2 levels"))

    bcdiag(target .== uc[2], pred; rev = rev)
end

"""
    biasplot(x::BCDiag)

returns a bias calibration plot of `x` - actual response vs. predicted response
"""
function biasplot(x::BCDiag)
    dmin = min(minimum(x.rrPrd), minimum(x.rrObs))
    dmax = max(maximum(x.rrPrd), maximum(x.rrObs))

    plt = plot(size = (500, 500), aspect_ratio = :equal, legend = :bottomright,
        xguidefontsize = 10, yguidefontsize = 10, titlefontsize = 11)

    title!("Model Calibration Graph", titlefontsize = 11)
    xlabel!("Predicted Response Rate", xguidefontsize = 10)
    ylabel!("Observed Response Rate", yguidefontsize = 10)

    plot!(x.rrPrd, x.rrObs, label = nothing, seriescolor = :blue)
    plot!([dmin, dmax], [dmin, dmax], label = nothing, width = 0.5, seriescolor = :black)

    plt
end

"""
    ksplot(x::BCDiag)

returns a KS plot of `x` - CDF1 (True Positive) and CDF0 (False Positive) versus depth
"""
function ksplot(x::BCDiag)
    plt = plot(size = (500, 500), aspect_ratio = :equal, xlims = (0., 1.), ylims = (0., 1.),
        legend = :bottomright, xguidefontsize = 10, yguidefontsize = 10, titlefontsize = 11)

    title!("KS Graph", titlefontsize = 11)
    xlabel!("Depth", xguidefontsize = 10)
    ylabel!("Empirical CDF", yguidefontsize = 10)

    plot!([0.; x.depth], [0.; x.cdf1], label = "Class 1", width = 1.5, seriescolor = :blue)
    plot!([0.; x.depth], [0.; x.cdf0], label = "Class 0", width = 0.5, seriescolor = :black)
    plot!([0., x.baserate, 1.], [0., 1., 1.], label = "Perfect Class 1", width = 0.5, seriescolor = :blue, linestyle = :dash)
    plot!([0., x.baserate, 1.], [0., 0., 1.], label = "Perfect Class 0", width = 0.5, seriescolor = :black, linestyle = :dash)
    ksidx = findlast(x.depth .<= x.ksdep)
    plot!([x.depth[ksidx], x.depth[ksidx]], [x.cdf0[ksidx], x.cdf1[ksidx]], label = nothing, width = 0.8, linestyle = :dash)

    annotate!((0.03, 0.97, text("Base rate = $(round(x.baserate, digits = 4))", :black, :left, 9)))
    annotate!((0.03, 0.93, text("Depth = $(round(x.ksdep, digits = 4))", :black, :left, 9)))
    annotate!((0.03, 0.89, text("KSarg = $(round(x.ksarg, digits = 4))", :black, :left, 9)))
    annotate!((0.03, 0.85, text("KS = $(round(x.ks, digits = 4))", :black, :left, 9)))

    plt
end

"""
    rocplot(x::BCDiag)

returns a ROC plot of `x` - CDF1 (True Positive) vs. CDF0 (False Positive)
"""
function rocplot(x::BCDiag)
    plt = plot(size = (500, 500), aspect_ratio = :equal, xlims = (0., 1.), ylims = (0., 1.),
        legend = :bottomright, xguidefontsize = 10, yguidefontsize = 10, titlefontsize = 11)

    title!("ROC Graph", titlefontsize = 11)
    xlabel!("False Positive Rate", xguidefontsize = 10)
    ylabel!("True Positive Rate", yguidefontsize = 10)

    plot!([0.; x.cdf0], [0.; x.cdf1], label = "Model", seriescolor = :blue, width = 1.5)
    plot!([0., 1.], [0., 1.], label = "Random", seriescolor = :black, width = 0.3, linestyle = :dash)

    annotate!((0.75, 0.22, text("auc = $(round(x.auc, digits = 4))", :black, :left, 9)))
    annotate!((0.75, 0.18, text("Gini = $(round(2. * x.auc - 1., digits = 4))", :black, :left, 9)))

    plt
end

"""
    accuracyplot(x::BCDiag; util=[1, 0, 0, 1])

Using `util` values for [TP, FN, FP, TN], produce accuracy plot and its [max, argmax, argdep].
Default `util` values of [1, 0, 0, 1] gives the standard accuracy value of (TP+TN)/N.
"""
function accuracyplot(x::BCDiag; util = [1, 0, 0, 1])
    p = x.baserate
    q = 1. - p

    # perfect classifier
    up = dot([p, 0., 0., q], util)                  # TP, FN, FP, TN of perfect classifier

    # 0.5 / 0.5 random classifier
    # ur5 = dot([0.5p, 0.5p, 0.5q, 0.5q], util)

    # baserate random classifier
    urb = dot([p*p, p*q, p*q, q*q], util)

    # model utility value [TP, FN, FP, TN]
    uc = [p.*x.cdf1  p.*(1. .- x.cdf1)  q.*x.cdf0  q.*(1. .- x.cdf0)] * util
    ucmax = findmax(uc)

    plt = plot(size = (500, 500), xlims = (0., 1.), legend = :best,
        xguidefontsize = 10, yguidefontsize = 10, titlefontsize = 11)

    title!("Accuracy Plot", titlefontsize = 11)
    xlabel!("Depth", xguidefontsize = 10)
    ylabel!("Utility", yguidefontsize = 10)

    hline!([up], label = "Perfect $(round(up, sigdigits = 3))", linestyle = :dash)

    plot!(x.depth, uc, width = 1.5, seriescolor = :blue,
        label = "Model $(round(ucmax[1], sigdigits = 3))")
    vline!([x.depth[ucmax[2]]], label = "Depth $(x.depth[ucmax[2]])")

    hline!([urb], label = "Base Random $(round(urb, sigdigits = 3))", linestyle = :dash)

    plt
end

"""
    liftcurve(x::BCDiag)

returns a lift curve plot of `x` - actual and predicted versus depth
"""
function liftcurve(x::BCDiag)
    plt = Plots.plot(size = (500, 500), xlims = (0., 1.),
        legend = :topright, xguidefontsize = 10, yguidefontsize = 10, titlefontsize = 11)

    title!("Lift Curve", titlefontsize = 11)
    xlabel!("Depth", xguidefontsize = 10)
    ylabel!("Response Rate", yguidefontsize = 10)

    plot!(x.depth, x.rrObs, label = "Actual", width = 0.5, linestyle = :dash)
    plot!(x.depth, x.rrPrd, label = "Predicted", width = 1.5, seriescolor = :blue)
    hline!([x.baserate], label = "Base rate", linestyle = :dash)

    plt
end

"""
    cumliftcurve(x::BCDiag)

returns a cumulative lift curve plot of `x` - cumulative actual and predicted vs. depth
"""
function cumliftcurve(x::BCDiag)
    plt = Plots.plot(size = (500, 500), xlims = (0., 1.),
        legend = :topright, xguidefontsize = 10, yguidefontsize = 10, titlefontsize = 11)

    title!("Cumulative Lift Curve", titlefontsize = 11)
    xlabel!("Depth", xguidefontsize = 10)
    ylabel!("Response Rate", yguidefontsize = 10)

    plot!(x.depth, x.crObs, label = "Actual", width = 0.5, linestyle = :dash)
    plot!(x.depth, x.crPrd, label = "Predicted", width = 1.5, seriescolor = :blue)
    hline!([x.baserate], label = "Base rate", linestyle = :dash)

    plt
end

"""
    liftable(x::BCDiag)

returns a lift table of `x` as a DataFrame
"""
function liftable(x::BCDiag)
    out = DataFrame(grp = x.grp, depth = x.depth, count = x.cntg,
                    cntObs = x.cnt1, cntPrd = x.prd1,
                    rrObs = x.rrObs, rrPred = x.rrPrd,
                    liftObs = x.rrObs ./ x.baserate, liftPrd = x.rrPrd ./ x.baserate
                    )

    out
end

"""
    cumliftable(x::BCDiag)

returns a cumulative lift table of `x` as a DataFrame
"""
function cumliftable(x::BCDiag)
    out = DataFrame(grp = x.grp, depth = x.depth, count = x.cumg,
                    cumObs = x.cum1, cumPrd = x.cpr1,
                    crObs = x.cum1 ./ x.cumg,
                    crPrd = x.cpr1 ./ x.cumg,
                    )
    out.liftObs = out.crObs ./ x.baserate
    out.liftPrd = out.crPrd ./ x.baserate

    out
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 ApplicationGatewayProbePropertiesFormat <: SwaggerModel
    protocol::Any # spec type: Union{ Nothing, ApplicationGatewayProtocol } # spec name: protocol
    host::Any # spec type: Union{ Nothing, String } # spec name: host
    path::Any # spec type: Union{ Nothing, String } # spec name: path
    interval::Any # spec type: Union{ Nothing, Int32 } # spec name: interval
    timeout::Any # spec type: Union{ Nothing, Int32 } # spec name: timeout
    unhealthyThreshold::Any # spec type: Union{ Nothing, Int32 } # spec name: unhealthyThreshold
    pickHostNameFromBackendHttpSettings::Any # spec type: Union{ Nothing, Bool } # spec name: pickHostNameFromBackendHttpSettings
    minServers::Any # spec type: Union{ Nothing, Int32 } # spec name: minServers
    match::Any # spec type: Union{ Nothing, ApplicationGatewayProbeHealthResponseMatch } # spec name: match
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState
    port::Any # spec type: Union{ Nothing, Int32 } # spec name: port

    function ApplicationGatewayProbePropertiesFormat(;protocol=nothing, host=nothing, path=nothing, interval=nothing, timeout=nothing, unhealthyThreshold=nothing, pickHostNameFromBackendHttpSettings=nothing, minServers=nothing, match=nothing, provisioningState=nothing, port=nothing)
        o = new()
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("interval"), interval)
        setfield!(o, Symbol("interval"), interval)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("timeout"), timeout)
        setfield!(o, Symbol("timeout"), timeout)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("unhealthyThreshold"), unhealthyThreshold)
        setfield!(o, Symbol("unhealthyThreshold"), unhealthyThreshold)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("pickHostNameFromBackendHttpSettings"), pickHostNameFromBackendHttpSettings)
        setfield!(o, Symbol("pickHostNameFromBackendHttpSettings"), pickHostNameFromBackendHttpSettings)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("minServers"), minServers)
        setfield!(o, Symbol("minServers"), minServers)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("match"), match)
        setfield!(o, Symbol("match"), match)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(ApplicationGatewayProbePropertiesFormat, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        o
    end
end # type ApplicationGatewayProbePropertiesFormat

const _property_map_ApplicationGatewayProbePropertiesFormat = Dict{Symbol,Symbol}(Symbol("protocol")=>Symbol("protocol"), Symbol("host")=>Symbol("host"), Symbol("path")=>Symbol("path"), Symbol("interval")=>Symbol("interval"), Symbol("timeout")=>Symbol("timeout"), Symbol("unhealthyThreshold")=>Symbol("unhealthyThreshold"), Symbol("pickHostNameFromBackendHttpSettings")=>Symbol("pickHostNameFromBackendHttpSettings"), Symbol("minServers")=>Symbol("minServers"), Symbol("match")=>Symbol("match"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("port")=>Symbol("port"))
const _property_types_ApplicationGatewayProbePropertiesFormat = Dict{Symbol,String}(Symbol("protocol")=>"ApplicationGatewayProtocol", Symbol("host")=>"String", Symbol("path")=>"String", Symbol("interval")=>"Int32", Symbol("timeout")=>"Int32", Symbol("unhealthyThreshold")=>"Int32", Symbol("pickHostNameFromBackendHttpSettings")=>"Bool", Symbol("minServers")=>"Int32", Symbol("match")=>"ApplicationGatewayProbeHealthResponseMatch", Symbol("provisioningState")=>"ProvisioningState", Symbol("port")=>"Int32")
Base.propertynames(::Type{ ApplicationGatewayProbePropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayProbePropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayProbePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayProbePropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayProbePropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayProbePropertiesFormat[property_name]

function check_required(o::ApplicationGatewayProbePropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewayProbePropertiesFormat }, name::Symbol, val)
    if name === Symbol("port")
        Swagger.validate_param(name, "ApplicationGatewayProbePropertiesFormat", :maximum, val, 65535, false)
        Swagger.validate_param(name, "ApplicationGatewayProbePropertiesFormat", :minimum, val, 1, false)
    end
end


text:
using AdFem

# create a uniform 10x10 mesh on [0,1]^2
mmesh = Mesh(10,10,0.1)

# xy is a Nx2 matrix for Gauss point coordinates,
# where N is the number of Gauss points
xy = gauss_nodes(mmesh)

# fc_init initializes a DNN weights and biases with 3 hidden layers, 20 neurons per layers
# input dimension is 2 (coordinate), and output is 1 (scalar κθ value)
θ = Variable(fc_init([2, 20,20,20,1]))
κ = squeeze(fc(xy, [20,20,20,1], θ))

# construct the coefficient matrix from κ
A = compute_fem_laplace_matrix1(κ, mmesh)

# assume the source function is given by F(x,y) = x + y
src = xy[:,1] + xy[:,2]
f = compute_fem_source_term1(src, mmesh)

# impose homogenous Dirichlet boundary condition
A, f = impose_Dirichlet_boundary_conditions(A, f, bcnode(mmesh), zeros(length(bcnode(mmesh))))

# solves the linear system A*u = f
u = A\f

# Let's generate some random observation data on each FEM node
uobs = rand(zeros(mmesh.nnode))

loss = sum((u-uobs)^2)

# construct the gradient
g = gradients(loss, θ)

sess = Session(); init(sess)

# a one-liner for L-BFGS-B optimization
BFGS!(sess, loss, method = "L-BFGS-B")


text:
"""
Author: Michael Palumbo
Created: December 2019
Contact: mlp95@psu.edu
Based on code by Alex Wise (aw@psu.edu)
Refactors & optimized by Eric Ford
"""

"""  Functor to estimate RV based on fitting quadratic near minimum of CCF.
TODO: Revist the logic here and see if need to perform transformation first.
"""
struct MeasureRvFromCCFQuadratic <: AbstractMeasureRvFromCCF
    frac_of_width_to_fit::Float64
    measure_width_at_frac_depth::Float64
end

"""
Construct functor to estimate RV based on the CCF.
Optional Arguments:
* `frac_of_width_to_fit`: (0.5)
* `measure_width_at_frac_depth`: (0.5)
"""
function MeasureRvFromCCFQuadratic(; frac_of_width_to_fit::Real = default_frac_of_width_to_fit,
                                    measure_width_at_frac_depth::Real = default_measure_width_at_frac_depth )
    @assert 0.25 <= measure_width_at_frac_depth <= 0.75
    @assert 0.1 <= frac_of_width_to_fit <= 2.0
    MeasureRvFromCCFQuadratic(frac_of_width_to_fit,measure_width_at_frac_depth)
end

using Statistics
function (mrv::MeasureRvFromCCFQuadratic)(vels::A1, ccf::A2 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1} }
    # find the min and fit only the part near the minimum of the CCF
    amin, inds = find_idx_at_and_around_minimum(vels, ccf, frac_of_width_to_fit=mrv.frac_of_width_to_fit, measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)
    if isnan(amin)     return ( rv=NaN, σ_rv=NaN )     end

    mean_v = mean(view(vels,inds))
    X = ones(length(inds),3)
    X[:,2] .= view(vels,inds) .-mean_v
    X[:,3] .= (view(vels,inds) .-mean_v).^2
    (c, b, a)  = (X'*X) \ (X'*view(ccf,inds))
    #=
    # do the polyfit
    pfit = Polynomials.fit(vels[inds], ccf[inds], 2)
    @assert length(Polynomials.coeffs(pfit)) >= 3   # just in case fails to fit a quadratic

    # get center from coeffs
    c, b, a = Polynomials.coeffs(pfit)
    @assert a>0
    =#
    v_at_min_of_quadratic = -b/(2*a) + mean_v
    return ( rv=v_at_min_of_quadratic, σ_rv=NaN )
end

function (mrv::MeasureRvFromCCFQuadratic)(vels::A1, ccf::A2, ccf_var::A2 ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}, T3<:Real, A3<:AbstractArray{T3,1} }
    # find the min and fit only the part near the minimum of the CCF
    amin, inds = find_idx_at_and_around_minimum(vels, ccf, frac_of_width_to_fit=mrv.frac_of_width_to_fit, measure_width_at_frac_depth=mrv.measure_width_at_frac_depth)
    if isnan(amin)     return ( rv=NaN, σ_rv=NaN )     end
    mean_v = mean(view(vels,inds))
    X = ones(length(inds),3)
    X[:,2] .= view(vels,inds) .-mean_v
    X[:,3] .= (view(vels,inds) .-mean_v).^2
    covar = PDiagMat(view(ccf_var,inds))
    denom = (X' * (covar \ X) )
    (c, b, a)  = denom \ (X' * (covar \ view(ccf,inds)) )
    covar_beta = inv(denom)
    v_at_min_of_quadratic = -b/(2*a)
    sigma_rv = abs(v_at_min_of_quadratic) * sqrt(covar_beta[2,2]/b^2+covar_beta[3,3]/a^2)
    v_at_min_of_quadratic += mean_v
    #=
    # do the polyfit
    pfit = Polynomials.fit(vels[inds], ccf[inds], 2)
    @assert length(Polynomials.coeffs(pfit)) >= 3   # just in case fails to fit a quadratic

    # get center from coeffs
    c, b, a = Polynomials.coeffs(pfit)
    =#
    return ( rv=v_at_min_of_quadratic, σ_rv=sigma_rv  )
end


text:
@defcomp impactdiarrhoea begin
    regions = Index()

    diadead = Variable(index=[time,regions])
    diasick = Variable(index=[time,regions])

    diamort = Parameter(index=[regions])
    diamortel = Parameter(default = -1.578625)
    diamortnl = Parameter(default = 1.141714)

    diayld = Parameter(index=[regions])
    diayldel = Parameter(default = -0.418406)
    diayldnl = Parameter(default = 0.699241)

    income = Parameter(index=[time,regions])
    population = Parameter(index=[time,regions])
    gdp90 = Parameter(index=[regions])
    pop90 = Parameter(index=[regions])

    temp90 = Parameter(index=[regions])
    bregtmp = Parameter(index=[regions])
    regtmp = Parameter(index=[time,regions])

    function run_timestep(p, v, d, t)

        for r in d.regions
            ypc = 1000.0 * p.income[t, r] / p.population[t, r]
            ypc90 = 1000.0 * p.gdp90[r] / p.pop90[r]

            # 0.49 is the increase in global temperature from pre-industrial to 1990
            absoluteRegionalTempPreIndustrial = p.temp90[r] - 0.49 * p.bregtmp[r]

            if absoluteRegionalTempPreIndustrial > 0.0
                v.diadead[t, r] = p.diamort[r] * p.population[t, r] * (ypc / ypc90)^p.diamortel * (((absoluteRegionalTempPreIndustrial + p.regtmp[t, r]) / absoluteRegionalTempPreIndustrial)^p.diamortnl - 1.0)

                v.diasick[t, r] = p.diayld[r] * p.population[t, r] * (ypc / ypc90)^p.diayldel * (((absoluteRegionalTempPreIndustrial + p.regtmp[t, r]) / absoluteRegionalTempPreIndustrial)^p.diayldnl - 1.0)
            else
                v.diadead[t, r] = 0.0
                v.diasick[t, r] = 0.0
            end
        end
    end
end


text:
# ---
# title: 470. Implement Rand10() Using Rand7()
# id: problem470
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Random, Rejection Sampling
# link: <https://leetcode.com/problems/implement-rand10-using-rand7/description/>
# hidden: true
# ---
#
# Given the **API** `rand7()` that generates a uniform random integer in the
# range `[1, 7]`, write a function `rand10()` that generates a uniform random
# integer in the range `[1, 10]`. You can only call the API `rand7()`, and you
# shouldn't call any other API. Please **do not** use a language's built-in
# random API.
#
# Each test case will have one **internal** argument `n`, the number of times
# that your implemented function `rand10()` will be called while testing. Note
# that this is **not an argument** passed to `rand10()`.
#
# **Follow up:**
#
#   * What is the [expected value](https://en.wikipedia.org/wiki/Expected_value) for the number of calls to `rand7()` function?
#   * Could you minimize the number of calls to `rand7()`?
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1
#     Output: [2]
#
#
# **Example 2:**
#
#
#
#     Input: n = 2
#     Output: [2,8]
#
#
# **Example 3:**
#
#
#
#     Input: n = 3
#     Output: [3,8,10]
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 105`
#
#
## @lc code=start
using LeetCode

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


text:
istrivial(x) = isbits(x)
istrivial(x::Type) = true
istrivial(x::String) = true

function inline_consts!(ir::IR)
  env = Dict()
  us = IRTools.Inner.usecounts(ir)
  isused(x) = get(us, x, 0) > 0
  for v in reverse(keys(ir))
    st = ir[v]
    effect = effectful(exprtype.((ir,), ir[v].expr.args)...)
    if st.type isa Union{Const,Node} && !effect
      map(v -> v isa Variable && (us[v] -= 1), st.expr.args)
      if st.type isa Node || istrivial(st.type.value) || !isused(v)
        delete!(ir, v)
        env[v] = st.type.value
      else
        ir[v] = st.type.value
      end
    elseif effect && st.type isa Node
      ir[v] = stmt(st, type = widen(st.type))
    end
  end
  return IRTools.prewalk!(x -> get(env, x, x), ir)
end

effectful(f, args...) = false

# TODO: handle control flow
# should be a dataflow analysis similar to blockarg pruning
function deadcode!(ir::IR)
  us = IRTools.Inner.usecounts(ir)
  isused(x) = get(us, x, 0) > 0
  for v in reverse(keys(ir))
    if !isused(v)
      if isexpr(ir[v].expr)
        effectful(exprtype.((ir,), ir[v].expr.args)...) && continue
        map(v -> v isa Variable && (us[v] -= 1), ir[v].expr.args)
      elseif ir[v].expr isa Variable
        us[ir[v].expr] -= 1
      end
      delete!(ir, v)
    end
  end
  return ir
end

function partials!(ir::IR)
  slots = IdDict()
  slot(k, T) = get!(() -> Slot(gensym(:partial), T), slots, k)
  for (v, st) in ir
    ex = st.expr
    if iscall(ex, setindex!)
      exprtype(ir, ex.args[2]) isa Partial{<:Dict} || continue
      key = exprtype(ir, ex.args[2]), exprtype(ir, ex.args[4])::Const
      T = key[1].value[key[2].value]
      insert!(ir, v, Expr(:(=), slot(key, T), ex.args[3]))
    elseif iscall(ex, getindex)
      exprtype(ir, ex.args[2]) isa Partial{<:Dict} || continue
      key = exprtype(ir, ex.args[2]), exprtype(ir, ex.args[3])::Const
      T = key[1].value[key[2].value]
      ir[v] = slot(key, T)
    end
  end
  for (v, st) in ir
    st.type isa Partial{<:Dict} && delete!(ir, v)
  end
  return ir
end

function trimblocks!(ir)
  i = 1
  while i <= length(blocks(ir))
    bl = block(ir, i)
    # assume no implicit branches
    # don't bother with blocks that have args (for now)
    if isempty(bl) && length(branches(bl)) == 1 && isempty(arguments(bl))
      for a in predecessors(bl), i in 1:length(branches(a))
        br = branches(bl)[1]
        branches(a)[i].block == bl.id || continue
        branches(a)[i] = Branch(branches(a)[i], block = br.block, args = br.args)
      end
      deleteblock!(ir, i)
    else
      i += 1
    end
  end
  return ir
end

cleanup!(ir) =
  ir |> inline_consts! |> partials! |> ssa! |> prune! |> deadcode! |> renumber


text:
type Activation
    th::Array{Float64,1}
    tau1::Float64
    tau2::Float64
    m1::Float64
    m2::Float64
    k::Array{Float64,1}

    function Activation(old=Dict("a"=>0),restart="no")
        this = new()
        this.m1 = 1.32;
        this.m2 = 27.4;
        if restart == "no"
            this.th = [0.8];
            this.tau1 = 0.269*this.th[1];
            this.tau2 = 0.452*this.th[1];
            t = linspace(0,this.th[1],10000);
            g1 = (t/this.tau1).^this.m1;
            g2 = (t/this.tau2).^this.m2;
            this.k = [maximum((g1./(1+g1)).*(1./(1+g2)))^-1];
        elseif restart == "yes"
            this.th = [(old["th"])[end]];
            this.tau1 = old["tau1"];
            this.tau2 = old["tau2"];
            this.k = [old["k"][end]];
        end
        return this
    end
end

type LeftVentricle
    V0::Float64
    Emin::Float64
    V::Vector{Float64}
    P::Vector{Float64}
    E::Vector{Float64}
    Emax::Vector{Float64}

    function LeftVentricle(old=Dict("a"=>0),restart="no")
        this = new()
        this.V0 = 10*cm3Tom3;
        this.Emin = 0.0283*mmHgToPa/cm3Tom3;
        if restart == "no"
            this.Emax = [3*mmHgToPa/cm3Tom3];
        elseif restart == "yes"
            this.Emax = [old["Emax"][end]];
        end
        this.V = Vector{Float64}[];
        this.P = Vector{Float64}[];
        this.E = Vector{Float64}[];
        return this
    end
end

type LeftAtrium
    V0::Float64
    E::Float64
    R::Float64
    L::Float64
    V::Vector{Float64}
    P::Vector{Float64}
    Q::Vector{Float64}

    function LeftAtrium()
        this = new()
        this.V0 = 10*cm3Tom3;
        this.R = 3.6e-3*mmHgToPa/cm3Tom3;
        this.L = 3e-5*mmHgToPa/cm3Tom3;
        this.E = 0.13*mmHgToPa/cm3Tom3;
        this.V = Vector{Float64}[];
        this.P = Vector{Float64}[];
        this.Q = Vector{Float64}[];
        return this
    end
end

type RightVentricle
    V0::Float64
    Emin::Float64
    L::Float64
    V::Vector{Float64}
    P::Vector{Float64}
    E::Vector{Float64}
    Q::Vector{Float64}
    Emax::Vector{Float64}

    function RightVentricle(old=Dict("a"=>0),restart="no")
        this = new()
        this.V0 = 10*cm3Tom3;
        this.Emin = 0.0283*mmHgToPa/cm3Tom3;
        if restart == "no"
            this.Emax = [0.4*mmHgToPa/cm3Tom3];
        elseif restart == "yes"
            this.Emax = [old["Emax"][end]];
        end
        this.L = 2.16e-4*mmHgToPa/cm3Tom3;
        this.V = Vector{Float64}[];
        this.P = Vector{Float64}[];
        this.E = Vector{Float64}[];
        this.Q = Vector{Float64}[];
        return this
    end
end

type RightAtrium
    V0::Float64
    E::Float64
    R::Float64
    L::Float64
    V::Vector{Float64}
    P::Vector{Float64}
    Q::Vector{Float64}

    function RightAtrium()
        this = new()
        this.V0 = 10*cm3Tom3;
        this.R = 4.85e-3*mmHgToPa/cm3Tom3;
        this.L = 5e-5*mmHgToPa/cm3Tom3;
        this.E = 0.16*mmHgToPa/cm3Tom3;
        this.V = Vector{Float64}[];
        this.P = Vector{Float64}[];
        this.Q = Vector{Float64}[];
        return this
    end
end

type AorticValve
    Kvo::Float64
    Kvc::Float64
    leff::Float64
    Po::Float64
    Pc::Float64
    Aann::Float64
    Ks::Float64
    zeta::Vector{Float64}

    function AorticValve()
        this = new()
        this.Kvo = 0.16;
        this.Kvc = 0.16;
        this.leff = 0.01;
        this.Po = 0;
        this.Pc = 0;
        this.Aann = 3.8e-4;
        this.Ks = 4e-9/cm3Tom3;
        this.zeta = Vector{Float64}[];
        return this
    end
end

type Heart
    activation::CVModule.Activation
    lv::CVModule.LeftVentricle
    la::CVModule.LeftAtrium
    rv::CVModule.RightVentricle
    ra::CVModule.RightAtrium
    av::CVModule.AorticValve

    function Heart(old=Dict("a"=>0),restart="no")
        this = new()
        if restart == "yes"
            act = old["activation"];
            lv = old["lv"];
            rv = old["rv"];
            this.activation = CVModule.Activation(act,restart);
            this.lv = CVModule.LeftVentricle(lv,restart);
            this.rv = CVModule.RightVentricle(rv,restart);
        elseif restart == "no"
            this.activation = CVModule.Activation();
            this.lv = CVModule.LeftVentricle();
            this.rv = CVModule.RightVentricle();
        end
        this.la = CVModule.LeftAtrium();
        this.ra = CVModule.RightAtrium();
        this.av = CVModule.AorticValve();
        return this
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
module TestDataFrame
    using Base.Test
    using DataArrays
    using DataFrames

    #test_group("Operations on DataFrames that have column groupings")

    x = DataFrame(quote
        a = [1,2,3]
        b = [4,5,6]
    end)
    y = DataFrame(quote
        c = [1,2,3]
        d = [4,5,6]
    end)

    # set_group(x, "group1", ["a", "b"])
    # set_group(y, "group2", ["c", "d"])

    z = deepcopy(x)
    # @test is_group(z, "group1")

    z = hcat(x, y)
    # @test is_group(z, "group1")
    # @test is_group(z, "group2")

    v = DataFrame(quote
        a = [5,6,7]
        b = [8,9,10]
    end)
    z = vcat(DataFrame[v, x])
    # @test is_group(z, "group1")

    z = vcat(v, x)
    # @test is_group(z, "group1")

    # Deleting columns removes any mention from groupings
    delete!(x, "a")
    # @test colnames(x) == ["b"]
    # @test get_groups(x)["group1"] == ["b"]

    ## del calls ref, which properly deals with groupings
    y = without(y, "c")
    @test names(y) == ["d"]
    #@test get_groups(y)["group2"] == ["d"]
    z1 = z[[1]]
    @test names(z1) == ["a"]
    #@test get_groups(z1)["group1"] == ["a"]

    z2 = z[:,[1,1,2]]
    @test names(z2) == ["a", "a_1", "b"]
    #@test get_groups(z2)["group1"] == ["a_1", "b"]

    #test_group("DataFrame assignment")
    df1 = DataFrame(quote
        a = 1:5
        b2 = letters[1:5]
        v2 = randn(5)
    end)
    df2 = DataFrame(quote
        a = reverse([1:5])
        b2 = reverse(letters)[1:5]
        v2 = randn(5)
    end)
    df1[1:2,:] = df2[4:5,:]
    # TODO: Decide on ==
    # @test df1[1:2,:] == df2[4:5,:]

    #test_group("Empty DataFrame constructors")
    df = DataFrame(10, 5)
    @assert nrow(df) == 10
    @assert ncol(df) == 5
    @assert typeof(df[:, 1]) == DataVector{Float64}

    df = DataFrame(Int, 10, 3)
    @assert nrow(df) == 10
    @assert ncol(df) == 3
    @assert typeof(df[:, 1]) == DataVector{Int}
    @assert typeof(df[:, 2]) == DataVector{Int}
    @assert typeof(df[:, 3]) == DataVector{Int}

    df = DataFrame({Int, Float64, ASCIIString}, 100)
    @assert nrow(df) == 100
    @assert ncol(df) == 3
    @assert typeof(df[:, 1]) == DataVector{Int}
    @assert typeof(df[:, 2]) == DataVector{Float64}
    @assert typeof(df[:, 3]) == DataVector{ASCIIString}

    df = DataFrame({Int, Float64, ASCIIString}, ["A", "B", "C"], 100)
    @assert nrow(df) == 100
    @assert ncol(df) == 3
    @assert typeof(df[:, 1]) == DataVector{Int}
    @assert typeof(df[:, 2]) == DataVector{Float64}
    @assert typeof(df[:, 3]) == DataVector{ASCIIString}

    df = DataFrame(zeros(10, 5))
    @assert nrow(df) == 10
    @assert ncol(df) == 5
    @assert typeof(df[:, 1]) == DataVector{Float64}

    df = DataFrame(ones(10, 5))
    @assert nrow(df) == 10
    @assert ncol(df) == 5
    @assert typeof(df[:, 1]) == DataVector{Float64}

    df = DataFrame(eye(10, 5))
    @assert nrow(df) == 10
    @assert ncol(df) == 5
    @assert typeof(df[:, 1]) == DataVector{Float64}

    #test_group("Other DataFrame constructors")
    df = DataFrame([{"a"=>1, "b"=>'c'}, {"a"=>3, "b"=>'d'}, {"a"=>5}])
    @assert nrow(df) == 3
    @assert ncol(df) == 2
    @assert typeof(df[:,"a"]) == DataVector{Int}
    @assert typeof(df[:,"b"]) == DataVector{Char}

    df = DataFrame([{"a"=>1, "b"=>'c'}, {"a"=>3, "b"=>'d'}, {"a"=>5}], ["a", "b"])
    @assert nrow(df) == 3
    @assert ncol(df) == 2
    @assert typeof(df[:,"a"]) == DataVector{Int}
    @assert typeof(df[:,"b"]) == DataVector{Char}

    data = {"A" => [1, 2], "C" => ["1", "2"], "B" => [3, 4]}
    df = DataFrame(data)
    # Specify column_names
    df = DataFrame(data, ["C", "A", "B"])

    # This assignment was missing before
    df = DataFrame(quote Column = ["A"] end)
    df[1, "Column"] = "Testing"

    # flipud() tests
    df = DataFrame(quote A = 1:4 end)
    @assert isequal(flipud(df), DataFrame(quote A = [4, 3, 2, 1] end))
    @assert isequal(df, DataFrame(quote A = [1, 2, 3, 4] end))
    flipud!(df)
    @assert isequal(df, DataFrame(quote A = [4, 3, 2, 1] end))

    # zero-row dataframe and subdataframe test
    df = DataFrame(x=[], y=[])
    @assert nrow(df) == 0
    df = DataFrame(x=[1:3], y=[3:5])
    sdf = sub(df, :(x .== 4) )
    @assert nrow(sdf) == 0

    @assert hash(DataFrame([1 2; 3 4])) == hash(DataFrame([1 2; 3 4]))
    @assert hash(DataFrame([1 2; 3 4])) != hash(DataFrame([1 3; 2 4]))
end


text:
@testset "supported.jl" begin
    @testset "supported_objective" begin
        Dualization.supported_objective(lp1_test()) # ScalarAffineFucntion Objective
        Dualization.supported_objective(lp10_test()) # VariableIndex Objective
        Dualization.supported_objective(qp1_test()) # VariableIndex Objective
        obj_typ_svf = MOI.get(lp1_test(), MOI.ObjectiveFunctionType())
        obj_typ_saf = MOI.get(lp10_test(), MOI.ObjectiveFunctionType())
        obj_typ_qp = MOI.get(qp1_test(), MOI.ObjectiveFunctionType())
        @test Dualization.supported_obj(obj_typ_svf)
        @test Dualization.supported_obj(obj_typ_saf)
        @test Dualization.supported_obj(obj_typ_qp)
    end

    @testset "supported_constraints" begin
        # All supported SAFs
        con_types = MOI.get(lp1_test(), MOI.ListOfConstraintTypesPresent())
        Dualization.supported_constraints(con_types)

        # Intervals Set is not supported
        con_types = MOI.get(lp9_test(), MOI.ListOfConstraintTypesPresent())
        @test_throws ErrorException Dualization.supported_constraints(con_types) # Throws an error if constraint cannot be dualized
    end
end


text:
@inbounds @inline dist2(p, q)::Float64 = sqrt((p[1]-q[1])^2+(p[2]-q[2])^2)
# FMTree class
mutable struct FMTree
    s_init::Vec4f
    s_goal::Vec4f
    N #number of samples
    Pset::Vector{Vec4f} # Point set
    cost::Vector{Float64} #cost
    time::Vector{Float64} #optimal time to connect one node to its parent node
    parent::Vector{Int64} #parent node
    bool_unvisit::BitVector #logical value for Vunvisit
    bool_open::BitVector #logical value for Open
    bool_closed::BitVector #logical value for Closed
    world::World # simulation world config
    itr::Int64 # iteration num
    reach_goal_region::Bool
    end_index::Int64
    idx::Int64
    cannot_find_open_point::Bool

    function FMTree(s_init::Vec4f, s_goal::Vec4f, N, world, right_prob::Float64, left_prob::Float64, keep_prob::Float64, lane_id::Int64, frames)
        # constructer: sampling valid point from the configurationspace
        println("initializing fmt ...")
        reach_goal_region = false
        Pset = Vec4f[]
        push!(Pset, s_init) #inply idx_init = 1
        myrn(min, max) = min + (max-min)*rand()
        gaussianrn(mean, var) = rand(Normal(mean, var), 1)[1]
        distance(x0, x1) = sqrt((x1[2] - x0[2])^2 + (x1[1] - x0[1])^2)
        uniform_percent = 1.0
        N = N - 2
        N_uniform = Int(round(Int(round(N * uniform_percent))/5))
        N_prob = N - (N_uniform*5)
        N_prob = Int(floor(N_prob/11))
        N = N + 2

        # n_scenario = length(scenario)-1
        # sample_each_scenario = convert(Int64, floor(N_prob/n_scenario))
        # rest_point = N - 2 - N_uniform - sample_each_scenario*n_scenario
        for n = 1:N_uniform
            while(true)
                x = myrn(world.x_min[1], world.x_max[1])
                y = myrn(s_init[2], world.x_max[2])
                p = Vec4f(x, y,
                        myrn(world.v_min[1], world.v_max[1]),
                        myrn(world.v_min[2], world.v_max[2]))
                if isValid(world, p)
                    push!(Pset, p)
                    i = 1
                    while i < 5
                        p = Vec4f(x, y,
                            myrn(world.v_min[1], world.v_max[1]),
                            myrn(world.v_min[2], world.v_max[2]))
                        if isValid(world, p)
                            push!(Pset, p)
                            i += 1
                        end
                    end
                    break
                end
            end
        end
#=
        points = sample(s_init, s_goal, lane_id, frames, N_prob, 150.0)
        # println(length(points))
        # print(typeof(x))
        for p in points
            # println(p)
            if isValid(world, p)
                push!(Pset, p)
            end
            if length(Pset) == N-1
                println("lalala")
                break
            end
        end
=#
        N_rest = N - length(Pset) - 1
        for n = 1:N_rest
            while(true)
                x = myrn(world.x_min[1], world.x_max[1])
                y = myrn(s_init[2], world.x_max[2])
                p = Vec4f(x, y,
                        myrn(world.v_min[1], world.v_max[1]),
                        myrn(world.v_min[2], world.v_max[2]))
                if isValid(world, p)
                    push!(Pset, p)
                    break
                end
            end
        end

        #=
        for n = 1:N_prob
            if n != 1
                if scenario[n][2] == scenario[n-1][2]
                    continue
                end
            end
            init = [scenario[n][5],scenario[n][6]]
            init_v = scenario[n][12]
            if n != n_scenario
                next_p_after_init = [scenario[n+1][5],scenario[n+1][6]]
                init_angel = atan(next_p_after_init[2] - init[2], next_p_after_init[1] - init[1])
                x = Vec4f(init[1], init[2], init_v*cos(init_angel), init_v*sin(init_angel))
            else
                previous_p_before_init = [scenario[n-1][5],scenario[n-1][6]]
                init_angel = atan(init[2] - previous_p_before_init[2], init[1] - previous_p_before_init[1])
                x = Vec4f(init[1], init[2], init_v*cos(init_angel), init_v*sin(init_angel))
            end
            # println(x)
            for k = 1:sample_each_scenario
                p_x = gaussianrn(x[1], 1)
                p_y = gaussianrn(x[2], 1)
                x_v = gaussianrn(x[3], 1)
                y_v = gaussianrn(x[4], 1)
                p = Vec4f(p_x, p_y, x_v, y_v)
                push!(Pset, p)
            end
        end
        =#
        push!(Pset, s_goal) #inply idx_goal = N [last]
        cost = zeros(N)
        time = zeros(N)
        parent = ones(Int, N)
        bool_unvisit = trues(N)
        bool_unvisit[1] = false
        bool_closed = falses(N)
        bool_open = falses(N)
        bool_open[1] = true
        println("finish initializing")
        new(s_init, s_goal,
            N, Pset, cost, time, parent, bool_unvisit, bool_open, bool_closed, world, 0, false, 0, 0, false)
    end
end

function show(this::FMTree)
    println("drawing...")
    show(this.world)
    N = length(this.Pset)
    mat = zeros(2, N)
    for idx = 1:N
        mat[:, idx] = this.Pset[idx][1:2]
    end
    idxset_open = findall(this.bool_open)
    idxset_closed = findall(this.bool_closed)
    idxset_unvisit = findall(this.bool_unvisit)
    idxset_tree = setdiff(union(idxset_open, idxset_closed), [1])
    scatter(mat[1, 1], mat[2, 1], c=:blue, s=10, zorder = 100)
    scatter(mat[1, end], mat[2, end], c=:blue, s=10, zorder = 101)
    scatter(mat[1, idxset_open], mat[2, idxset_open], c=:gray, s=2)
    scatter(mat[1, idxset_closed], mat[2, idxset_closed], c=:gray, s=2)
    # println(size(mat))
    scatter(mat[1, idxset_unvisit], mat[2, idxset_unvisit], c=:orange, s=5)
    for idx in idxset_tree
        s0 = this.Pset[this.parent[idx]]
        s1 = this.Pset[idx]
        tau = this.time[idx]
        show_trajectory(s0, s1, tau)
    end

    scatter(mat[1, 1], mat[2, 1], c=:blue, s=20, zorder = 100)
    scatter(mat[1, end], mat[2, end], c=:blue, s=20, zorder = 101)
    println("s_init in show funciton:", this.s_init)
    scatter(this.s_init[1], this.s_init[2], c=:red, s=50)
    scatter(this.s_goal[1], this.s_goal[2], c=:red, s=50)

    # xlim(this.world.x_min[1]-5, this.world.x_max[1]+5)
    # ylim(this.world.x_min[2]-5, this.world.x_max[2]+5)
    println("finish drawing")
end

function solve(this::FMTree, frame_id::Int64, with_savefig = false)
    # keep extending the node until the tree reaches the goal
    println("please set with_savefig=false if you want to measure the computation time" )
    println("start solving")
    while(true)
        extend(this)
        if with_savefig
            if ((this.itr<100) & (this.itr % 20 == 1)) || (this.itr % 200==1)
                close()
                show(this)
                savefig("./fig/test_cvae/"*string(frame_id)*"-"*string(this.itr)*".png")
            end
        end
        if this.itr > 10000 || this.cannot_find_open_point
            opt_idx = 1
            open_set = vcat(findall(this.bool_open), findall(this.bool_closed))
            min_dist = 1500
            for i in open_set
                x = this.Pset[i]
                dist = sqrt((this.s_goal[2] - x[2])^2 + (this.s_goal[1] - x[1])^2)
                if dist < min_dist
                    opt_idx = i
                    min_dist = dist
                end
            end
            if opt_idx != 1
                this.reach_goal_region = true
                this.end_index = opt_idx
            else
                return Int64[]
            end
        end
        if !this.bool_unvisit[end]
            this.idx = this.N
            println("finish normally")
            break
        elseif this.reach_goal_region
            this.idx = this.end_index
            println("finish the goal region")
            break
        end
    end
    idx_solution = Int64[this.idx]
    while(true)
        println("idx:", this.idx)
        this.idx = this.parent[this.idx]
        push!(idx_solution, this.idx)
        this.idx == 1 && break
    end
    println("finish solving")
    return idx_solution
end

function extend(this::FMTree)
    # extend node
    this.itr += 1
    r = 200.0
    # if this.itr == 3
    #     exit()
    # end
    # println(this.itr)
    idxset_open = findall(this.bool_open)
    idxset_unvisit = findall(this.bool_unvisit)

    # println(this.cost)
    this.cannot_find_open_point = false
    try
        idx_lowest = idxset_open[findmin(this.cost[idxset_open])[2]]
    catch ArgumentError
        println("cannot find open point")
        this.cannot_find_open_point = true
        return
    end
    # println(idxset_unvisit)
    idx_lowest = idxset_open[findmin(this.cost[idxset_open])[2]]
    idxset_near, _, _ = filter_reachable(this.Pset, idxset_unvisit,
                                      this.Pset[idx_lowest], r, :F)

    for idx_near in idxset_near
        idxset_cand, distset_cand, timeset_cand = filter_reachable(this.Pset, idxset_open,
                                                    this.Pset[idx_near], r, :B)
        isempty(idxset_cand) && return
        cost_new, idx_costmin = findmin(this.cost[idxset_cand] + distset_cand)
        time_new = timeset_cand[idx_costmin] # optimal time for new connection
        idx_parent = idxset_cand[idx_costmin]
        waypoints = gen_trajectory(this.Pset[idx_parent], this.Pset[idx_near], time_new, 50)
        if isValid(this.world, waypoints)
            this.bool_unvisit[idx_near] = false
            this.bool_open[idx_near] = true
            this.cost[idx_near] = cost_new
            this.time[idx_near] = time_new
            this.parent[idx_near] = idx_parent
            x0 = this.Pset[idx_near]
            # println("goal region:", sqrt((this.s_goal[2] - x0[2])^2 + (this.s_goal[1] - x0[1])^2
            #        + (this.s_goal[3] - x0[3])^2 + (this.s_goal[4] - x0[4])^2))
            if sqrt((this.s_goal[2] - x0[2])^2 + (this.s_goal[1] - x0[1])^2 +
                (this.s_goal[3] - x0[3])^2 + (this.s_goal[4] - x0[4])^2) < 5
                this.reach_goal_region = true
                this.end_index = idx_near
                println("goal region:", sqrt((this.s_goal[2] - x0[2])^2 + (this.s_goal[1] - x0[1])^2
                        + (this.s_goal[3] - x0[3])^2 + (this.s_goal[4] - x0[4])^2))
                break
            end
        end
    end
    this.bool_open[idx_lowest] = false
    this.bool_closed[idx_lowest] = true
end


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

@doc raw"""NodeSystemInfo is a set of ids/uuids to uniquely identify the node.

    IoK8sApiCoreV1NodeSystemInfo(;
        architecture=nothing,
        bootID=nothing,
        containerRuntimeVersion=nothing,
        kernelVersion=nothing,
        kubeProxyVersion=nothing,
        kubeletVersion=nothing,
        machineID=nothing,
        operatingSystem=nothing,
        osImage=nothing,
        systemUUID=nothing,
    )

    - architecture::String : The Architecture reported by the node
    - bootID::String : Boot ID reported by the node.
    - containerRuntimeVersion::String : ContainerRuntime Version reported by the node through runtime remote API (e.g. docker://1.5.0).
    - kernelVersion::String : Kernel Version reported by the node from 'uname -r' (e.g. 3.16.0-0.bpo.4-amd64).
    - kubeProxyVersion::String : KubeProxy Version reported by the node.
    - kubeletVersion::String : Kubelet Version reported by the node.
    - machineID::String : MachineID reported by the node. For unique machine identification in the cluster this field is preferred. Learn more from man(5) machine-id: http://man7.org/linux/man-pages/man5/machine-id.5.html
    - operatingSystem::String : The Operating System reported by the node
    - osImage::String : OS Image reported by the node from /etc/os-release (e.g. Debian GNU/Linux 7 (wheezy)).
    - systemUUID::String : SystemUUID reported by the node. For unique machine identification MachineID is preferred. This field is specific to Red Hat hosts https://access.redhat.com/documentation/en-US/Red_Hat_Subscription_Management/1/html/RHSM/getting-system-uuid.html
"""
mutable struct IoK8sApiCoreV1NodeSystemInfo <: SwaggerModel
    architecture::Any # spec type: Union{ Nothing, String } # spec name: architecture
    bootID::Any # spec type: Union{ Nothing, String } # spec name: bootID
    containerRuntimeVersion::Any # spec type: Union{ Nothing, String } # spec name: containerRuntimeVersion
    kernelVersion::Any # spec type: Union{ Nothing, String } # spec name: kernelVersion
    kubeProxyVersion::Any # spec type: Union{ Nothing, String } # spec name: kubeProxyVersion
    kubeletVersion::Any # spec type: Union{ Nothing, String } # spec name: kubeletVersion
    machineID::Any # spec type: Union{ Nothing, String } # spec name: machineID
    operatingSystem::Any # spec type: Union{ Nothing, String } # spec name: operatingSystem
    osImage::Any # spec type: Union{ Nothing, String } # spec name: osImage
    systemUUID::Any # spec type: Union{ Nothing, String } # spec name: systemUUID

    function IoK8sApiCoreV1NodeSystemInfo(;architecture=nothing, bootID=nothing, containerRuntimeVersion=nothing, kernelVersion=nothing, kubeProxyVersion=nothing, kubeletVersion=nothing, machineID=nothing, operatingSystem=nothing, osImage=nothing, systemUUID=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("architecture"), architecture)
        setfield!(o, Symbol("architecture"), architecture)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("bootID"), bootID)
        setfield!(o, Symbol("bootID"), bootID)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("containerRuntimeVersion"), containerRuntimeVersion)
        setfield!(o, Symbol("containerRuntimeVersion"), containerRuntimeVersion)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("kernelVersion"), kernelVersion)
        setfield!(o, Symbol("kernelVersion"), kernelVersion)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("kubeProxyVersion"), kubeProxyVersion)
        setfield!(o, Symbol("kubeProxyVersion"), kubeProxyVersion)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("kubeletVersion"), kubeletVersion)
        setfield!(o, Symbol("kubeletVersion"), kubeletVersion)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("machineID"), machineID)
        setfield!(o, Symbol("machineID"), machineID)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("operatingSystem"), operatingSystem)
        setfield!(o, Symbol("operatingSystem"), operatingSystem)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("osImage"), osImage)
        setfield!(o, Symbol("osImage"), osImage)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("systemUUID"), systemUUID)
        setfield!(o, Symbol("systemUUID"), systemUUID)
        o
    end
end # type IoK8sApiCoreV1NodeSystemInfo

const _property_map_IoK8sApiCoreV1NodeSystemInfo = Dict{Symbol,Symbol}(Symbol("architecture")=>Symbol("architecture"), Symbol("bootID")=>Symbol("bootID"), Symbol("containerRuntimeVersion")=>Symbol("containerRuntimeVersion"), Symbol("kernelVersion")=>Symbol("kernelVersion"), Symbol("kubeProxyVersion")=>Symbol("kubeProxyVersion"), Symbol("kubeletVersion")=>Symbol("kubeletVersion"), Symbol("machineID")=>Symbol("machineID"), Symbol("operatingSystem")=>Symbol("operatingSystem"), Symbol("osImage")=>Symbol("osImage"), Symbol("systemUUID")=>Symbol("systemUUID"))
const _property_types_IoK8sApiCoreV1NodeSystemInfo = Dict{Symbol,String}(Symbol("architecture")=>"String", Symbol("bootID")=>"String", Symbol("containerRuntimeVersion")=>"String", Symbol("kernelVersion")=>"String", Symbol("kubeProxyVersion")=>"String", Symbol("kubeletVersion")=>"String", Symbol("machineID")=>"String", Symbol("operatingSystem")=>"String", Symbol("osImage")=>"String", Symbol("systemUUID")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1NodeSystemInfo }) = collect(keys(_property_map_IoK8sApiCoreV1NodeSystemInfo))
Swagger.property_type(::Type{ IoK8sApiCoreV1NodeSystemInfo }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1NodeSystemInfo[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NodeSystemInfo }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NodeSystemInfo[property_name]

function check_required(o::IoK8sApiCoreV1NodeSystemInfo)
    (getproperty(o, Symbol("architecture")) === nothing) && (return false)
    (getproperty(o, Symbol("bootID")) === nothing) && (return false)
    (getproperty(o, Symbol("containerRuntimeVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("kernelVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("kubeProxyVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("kubeletVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("machineID")) === nothing) && (return false)
    (getproperty(o, Symbol("operatingSystem")) === nothing) && (return false)
    (getproperty(o, Symbol("osImage")) === nothing) && (return false)
    (getproperty(o, Symbol("systemUUID")) === nothing) && (return false)
    true
end

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


text:
module SpecialPolynomials

using LinearAlgebra
#import LinearAlgebra: dot, ⋅
import SpecialFunctions: gamma

using Polynomials
import Polynomials:
    basis, isconstant, constantterm, assert_same_variable, StandardBasisPolynomial, ⟒
export basis

#import Intervals
#import Intervals: Open, Closed, Unbounded, bounds_types
import Polynomials: domain, Interval, Open, Closed, Unbounded, bounds_types
using QuadGK
using Memoize

using HypergeometricFunctions

include("utils.jl")
include("abstract.jl")

include("Orthogonal/abstract.jl")
include("Orthogonal/orthogonal.jl")
include("Orthogonal/cop.jl")
include("Orthogonal/ccop.jl")

include("Orthogonal/Bessel.jl")
include("Orthogonal/Chebyshev.jl")
include("Orthogonal/Hermite.jl")
include("Orthogonal/Laguerre.jl")
include("Orthogonal/Gegenbauer.jl")
include("Orthogonal/Jacobi.jl")
include("Orthogonal/Legendre.jl")
include("Orthogonal/WeightFunction.jl")

include("Orthogonal/Discrete/discrete-orthogonal.jl")
include("Orthogonal/Discrete/cdop.jl")
include("Orthogonal/Discrete/FallingFactorial.jl")
include("Orthogonal/Discrete/Hahn.jl")
include("Orthogonal/Discrete/Meixner.jl")
include("Orthogonal/Discrete/Krawchouk.jl")
include("Orthogonal/Discrete/Charlier.jl")
include("Orthogonal/Discrete/DiscreteChebyshev.jl")

include("Orthogonal/connection.jl")
include("Orthogonal/glaser-liu-rokhlin.jl")
include("Orthogonal/comrade.jl")

include("Interpolating/interpolating.jl")
include("Interpolating/Lagrange.jl")
include("Interpolating/Newton.jl")

include("Bernstein.jl")

using Requires
function __init__()
    @require FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c" include(
        "fasttransforms.jl",
    )
    @require FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" include(
        "fastgaussquadrature.jl",
    )
end

end # module


text:
module PartiallyTypedTables

import Tables
using ArgCheck: @argcheck

_typed(x) = getfield(x, :typed)
_rest(x) = getfield(x, :rest)

struct PartiallyTypedRowIterator{names,types,Typed<:NamedTuple{names,types},Rest}
    typed::Typed
    rest::Rest
end

_zip(tuples::NTuple{N,Any}...) where {N} = ntuple(i -> map(t -> t[i], tuples), N)

# Convince `julia` that all arguments are non-`nothing` when passed through:
@inline allsomethings(xs...) = _allsomethings((), xs...)
@inline _allsomethings(ys) = ys
# @inline _allsomethings(_, ::Nothing, xs...) = nothing
@inline function _allsomethings(ys, x, xs...)
    x === nothing && return nothing
    return _allsomethings((ys..., x), xs...)
end

_value(::Val{x}) where {x} = x

function Base.iterate(itr::PartiallyTypedRowIterator{names}, prev = nothing) where {names}
    if prev === nothing
        ys0 = map(iterate, Tuple(_typed(itr)))
        y = iterate(_rest(itr))
    else
        (typed_states, rest_state) = prev
        ys0 = map(iterate, Tuple(_typed(itr)), typed_states)
        y = iterate(_rest(itr), rest_state)
    end
    ys = allsomethings(ys0...)
    ys === nothing && return nothing
    y === nothing && return nothing
    typed = foldl(_zip(map(Val, names), ys); init = NamedTuple()) do typed, (name, (x, _))
        merge(typed, (; _value(name) => x))
    end
    state = (map(last, ys), y[2])
    return PartiallyTypedRow(typed, y[1]), state
end

eltypes(::Type{Tuple{}}) = Tuple{}
eltypes(::Type{T}) where {H,T<:Tuple{H,Vararg{Any}}} =
    Base.tuple_type_cons(eltype(H), eltypes(Base.tuple_type_tail(T)))

function Base.eltype(
    ::Type{PartiallyTypedRowIterator{names,types,Typed,Rest}},
) where {names,types,Typed,Rest}
    rtypes = eltypes(types)
    r = eltype(Rest)
    if isconcretetype(r)
        return PartiallyTypedRow{names,rtypes,NamedTuple{names,rtypes},r}
    else
        return PartiallyTypedRow{names,rtypes,NamedTuple{names,rtypes},<:r}
    end
end

Base.IteratorEltype(
    ::Type{PartiallyTypedRowIterator{names,types,Typed,Rest}},
) where {names,types,Typed,Rest} = Base.HasEltype()

Base.length(itr::PartiallyTypedRowIterator) = length(_rest(itr))

Base.IteratorSize(
    ::Type{PartiallyTypedRowIterator{names,types,Typed,Rest}},
) where {names,types,Typed,Rest} = Base.HasLength()
# TODO: "merge" `typed`

struct PartiallyTypedRow{names,types,Typed<:NamedTuple{names,types},Rest}
    typed::Typed
    rest::Rest
end

function Base.getproperty(row::PartiallyTypedRow, name::Symbol)
    haskey(_typed(row), name) && return _typed(row)[name]
    return getproperty(_rest(row), name)
end

function partiallytyped(table, names::Symbol...)
    @argcheck Tables.rowaccess(table)
    @argcheck Tables.columnaccess(table)
    columns = Tables.columns(table)
    return PartiallyTypedRowIterator(
        (; (n => getproperty(columns, n) for n in names)...),
        Tables.rows(table),
    )
end

end # module


text:
using DiffEqBase, Catalyst, OrdinaryDiffEq, Sundials, LinearAlgebra, SparseArrays
using ReactionNetworkImporters
using TimerOutputs

# parameters
const networkname = "testbcrbng"
const tf = 10000.

# finite difference
const build_jac = false
const sparse_jac = false
const zeroout_jac = false

# dense Jacobian
# const build_jac = true
# const sparse_jac = false
# const zeroout_jac = true

# sparse Jacobian
# const build_jac = true
# const sparse_jac = true
# const zeroout_jac = true

# BNG simulation data
datadir  = joinpath(@__DIR__,"../data/bcr")
fname    = joinpath(datadir, "bcr.net")

# we'll time the DiffEq solvers
const to = TimerOutput()
reset_timer!(to)

# BioNetGen network
@timeit to "bionetgen" prnbng = loadrxnetwork(BNGNetwork(),string(networkname,"bng"), fname);
rnbng = prnbng.rn; u0 = prnbng.u₀; p = prnbng.p;
@timeit to "baddodes" addodes!(rnbng; build_jac=build_jac, zeroout_jac=zeroout_jac, sparse_jac=sparse_jac, build_symfuncs=false, build_paramjac=false)
@timeit to "bODEProb" boprob = ODEProblem(rnbng, u0, (0.,tf), p)
u = copy(u0);
du = similar(u);
@timeit to "f1" rnbng.f(du,u,p,0.)
@timeit to "f2" rnbng.f(du,u,p,0.)

function solve_dense(boprob, u, p, build_jac)
    if build_jac
        J = zeros(length(u),length(u))
        @timeit to "J1" rnbng.jac(J,u,p,0.)
        @timeit to "J2" rnbng.jac(J,u,p,0.)
    end
    show(to)
    println()

    @timeit to "KenCarp4-LU-1" begin bsol = solve(boprob,KenCarp4(autodiff=false), abstol=1e-8, reltol=1e-8, saveat=1.); end; show(to)
    @timeit to "KenCarp4-LU-2" begin bsol = solve(boprob,KenCarp4(autodiff=false), abstol=1e-8, reltol=1e-8, saveat=1.); end; show(to)

    bsol
end

function solve_sparse(boprob, u, p)
    J = similar(rnbng.odefun.jac_prototype)
    @timeit to "J1" rnbng.jac(J,u,p,0.)
    @timeit to "J2" rnbng.jac(J,u,p,0.)
    show(to)
    println()

    @timeit to "KenCarp4-SLU-1" begin bsol = solve(boprob,KenCarp4(autodiff=false,linsolve=LinSolveFactorize(lu)), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
    @timeit to "KenCarp4-SLU-2" begin bsol = solve(boprob,KenCarp4(autodiff=false,linsolve=LinSolveFactorize(lu)), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)

    bsol
end

bsol = sparse_jac ? solve_sparse(boprob, u, p) : solve_dense(boprob, u, p, build_jac);

println()