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

# Deprecated functions and objects
#
# Please add new deprecations at the bottom of the file.
# A function deprecated in a release will be removed in the next one.
# Please also add a reference to the pull request which introduced the
# deprecation.
#
# For simple cases where a direct replacement is available, use @deprecate:
# the first argument is the signature of the deprecated method, the second one
# is the call which replaces it. Remove the definition of the deprecated method
# and unexport it, as @deprecate takes care of calling the replacement
# and of exporting the function.
#
# For more complex cases, move the body of the deprecated method in this file,
# and call depwarn() directly from inside it. The symbol depwarn() expects is
# the name of the function, which is used to ensure that the deprecation warning
# is only printed the first time for each call place.

macro deprecate(old, new, ex=true)
    meta = Expr(:meta, :noinline)
    if isa(old, Symbol)
        oldname = Expr(:quote, old)
        newname = Expr(:quote, new)
        Expr(:toplevel,
            ex ? Expr(:export, esc(old)) : nothing,
            :(function $(esc(old))(args...)
                  $meta
                  depwarn($"`$old` is deprecated, use `$new` instead.", Core.Typeof($(esc(old))).name.mt.name)
                  $(esc(new))(args...)
              end))
    elseif isa(old, Expr) && (old.head == :call || old.head == :where)
        remove_linenums!(new)
        oldcall = sprint(show_unquoted, old)
        newcall = sprint(show_unquoted, new)
        # if old.head is a :where, step down one level to the :call to avoid code duplication below
        callexpr = old.head == :call ? old : old.args[1]
        if callexpr.head == :call
            if isa(callexpr.args[1], Symbol)
                oldsym = callexpr.args[1]::Symbol
            elseif isa(callexpr.args[1], Expr) && callexpr.args[1].head == :curly
                oldsym = callexpr.args[1].args[1]::Symbol
            else
                error("invalid usage of @deprecate")
            end
        else
            error("invalid usage of @deprecate")
        end
        Expr(:toplevel,
            ex ? Expr(:export, esc(oldsym)) : nothing,
            :($(esc(old)) = begin
                  $meta
                  depwarn($"`$oldcall` is deprecated, use `$newcall` instead.", Core.Typeof($(esc(oldsym))).name.mt.name)
                  $(esc(new))
              end))
    else
        error("invalid usage of @deprecate")
    end
end

function depwarn(msg, funcsym)
    opts = JLOptions()
    if opts.depwarn == 2
        throw(ErrorException(msg))
    end
    deplevel = opts.depwarn == 1 ? CoreLogging.Warn : CoreLogging.BelowMinLevel
    @logmsg(
        deplevel,
        msg,
        _module=begin
            bt = backtrace()
            frame, caller = firstcaller(bt, funcsym)
            # TODO: Is it reasonable to attribute callers without linfo to Core?
            caller.linfo isa Core.MethodInstance ? caller.linfo.def.module : Core
        end,
        _file=String(caller.file),
        _line=caller.line,
        _id=(frame,funcsym),
        _group=:depwarn,
        caller=caller,
        maxlog=funcsym === nothing ? nothing : 1
    )
    nothing
end

firstcaller(bt::Vector, ::Nothing) = Ptr{Cvoid}(0), StackTraces.UNKNOWN
firstcaller(bt::Vector, funcsym::Symbol) = firstcaller(bt, (funcsym,))
function firstcaller(bt::Vector, funcsyms)
    # Identify the calling line
    found = false
    lkup = StackTraces.UNKNOWN
    found_frame = Ptr{Cvoid}(0)
    for frame in bt
        lkups = StackTraces.lookup(frame)
        for outer lkup in lkups
            if lkup == StackTraces.UNKNOWN || lkup.from_c
                continue
            end
            if found
                found_frame = frame
                @goto found
            end
            found = lkup.func in funcsyms
            # look for constructor type name
            if !found && lkup.linfo isa Core.MethodInstance
                li = lkup.linfo
                ft = ccall(:jl_first_argument_datatype, Any, (Any,), li.def.sig)
                if isa(ft,DataType) && ft.name === Type.body.name
                    ft = unwrap_unionall(ft.parameters[1])
                    found = (isa(ft,DataType) && ft.name.name in funcsyms)
                end
            end
        end
    end
    return found_frame, StackTraces.UNKNOWN
    @label found
    return found_frame, lkup
end

deprecate(m::Module, s::Symbol, flag=1) = ccall(:jl_deprecate_binding, Cvoid, (Any, Any, Cint), m, s, flag)

macro deprecate_binding(old, new, export_old=true, dep_message=:nothing, constant=true)
    dep_message === :nothing && (dep_message = ", use $new instead.")
    return Expr(:toplevel,
         export_old ? Expr(:export, esc(old)) : nothing,
         Expr(:const, Expr(:(=), esc(Symbol(string("_dep_message_",old))), esc(dep_message))),
         constant ? Expr(:const, Expr(:(=), esc(old), esc(new))) : Expr(:(=), esc(old), esc(new)),
         Expr(:call, :deprecate, __module__, Expr(:quote, old)))
end

macro deprecate_stdlib(old, mod, export_old=true, newname=old)
    rename = old === newname ? "" : " as `$newname`"
    dep_message = """: it has been moved to the standard library package `$mod`$rename.
                        Add `using $mod` to your imports."""
    new = GlobalRef(Base.root_module(Base, mod), newname)
    return Expr(:toplevel,
         export_old ? Expr(:export, esc(old)) : nothing,
         Expr(:const, Expr(:(=), esc(Symbol(string("_dep_message_",old))), esc(dep_message))),
         Expr(:const, Expr(:(=), esc(old), esc(new))),
         Expr(:call, :deprecate, __module__, Expr(:quote, old)))
end

macro deprecate_moved(old, new, export_old=true)
    eold = esc(old)
    emsg = string(old, " has been moved to the package ", new, ".jl.\n",
        "Run `Pkg.add(\"", new, "\")` to install it, restart Julia,\n",
        "and then run `using ", new, "` to load it.")
    return Expr(:toplevel,
        :($eold(args...; kwargs...) = error($emsg)),
        export_old ? Expr(:export, eold) : nothing,
        Expr(:call, :deprecate, __module__, Expr(:quote, old), 2))
end

# BEGIN 0.7 deprecations

function promote_eltype_op end

# END 0.7 deprecations

# BEGIN 1.0 deprecations

# END 1.0 deprecations

max_stars_repo_name: DragonXDoom/julia

max_issues_repo_name: DragonXDoom/julia

max_forks_repo_name: DragonXDoom/julia


text:
function nrow(m::Matrix)
  return size(m, 1)
end

function ncol(m::Matrix)
  return size(m, 2)
end

struct Data
  y::Vector{Matrix{Float64}}
  I::Int
  J::Int
  N::Vector{Int}
  m::Vector{Matrix{Int8}}  # TODO: make this Vector{Matrix{Bool}}

  function Data(y)
    @assert all(ncol.(y) .== ncol(y[1]))
    I = length(y)
    J = ncol(y[1])
    N = nrow.(y)
    m = [Int8(1) * isnan.(y[i]) for i in 1:I]
    return new(y, I, J, N, m)
  end
end

include("missing_mechanism.jl")
include("generate_data.jl")

ext: jl

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch


text:
using HTTP
using DataFrames

"""
    network(;apikey, bbox)

Returns shape for every street of the network within the bounding box.

# Arguments
apikey::String   : API key to access HERE
bbox::String     : Bounding box defined as `"lat₁,long₁;lat₂,long₂"`

For more information visit: https://developer.here.com/documentation/traffic/dev_guide/topics/guide.html

Example request: https://developer.here.com/documentation/traffic/dev_guide/topics_v6.1/example-flow-bounding-box.html
"""
function network(;apikey, bbox)
    # Request
    baseurl = "https://traffic.ls.hereapi.com"
    path = "traffic/6.2"
    resource = "flow.xml"
    url = "$baseurl/$path/$resource?apiKey=$apikey&bbox=$bbox&responseattributes=shape"
    r = HTTP.request(:GET, url)
    f = parsexml(String(r.body))
    rootNode = root(f)

    # Build network
    df = DataFrame(lid = Int64[], name = String[], length = Float64[], ffs = Float64[], shape = String[])
    for node in iterate(rootNode)
        if nodename(node) == "FI"
            shape = ""
            for childnode in eachelement(node)
                # Fetch name and length
                if childnode == firstelement(node)
                    push!(df[!, :lid], length(df[!, :lid]) + 1)
                    push!(df[!, :name], childnode["DE"])
                    push!(df[!, :length], parse(Float64, childnode["LE"]))
                # Fetch free flow speed
                elseif childnode == lastelement(node)
                    push!(df[!, :ffs], parse(Float64, childnode["FF"]))
                # Fetch link shape (Note: Link shape is of the same format as the bbox)
                else
                    shape *= replace(nodecontent(childnode), " " => ";")
                end
            end
            push!(df[!, :shape], shape)
        end
    end

    return df
end


text:
using JLD2
using BenchmarkTools
using Benchmarks

# Benchmark parameters

Nx = 256
Ny = 256
Nz = Nx

ranks = (1, 2, 4, 8, 16)

# Run and collect benchmarks

print_system_info()

for r in ranks
    @info "Benchmarking distributed incompressible model strong scaling [N=($Nx, $Ny, $Nz), ranks=$r]..."
    julia = Base.julia_cmd()
    run(`mpiexec -np $r $julia --project strong_scaling_incompressible_model_single.jl $Nx $Ny $Nz`)
end

suite = BenchmarkGroup(["size", "ranks"])
for r in ranks
    jldopen("strong_scaling_incompressible_model_$r.jld2", "r") do file
        suite[((Nx, Ny, Nz), r)] = file["trial"]
    end
end

# Summarize benchmarks

df = benchmarks_dataframe(suite)
sort!(df, :ranks)
benchmarks_pretty_table(df, title="Incompressible model strong scaling benchmark")

suite_Δ = speedups_suite(suite, base_case=((Nx, Ny, Nz), 1))
df_Δ = speedups_dataframe(suite_Δ)
sort!(df_Δ, :ranks)
benchmarks_pretty_table(df_Δ, title="Incompressible model strong scaling speedup")


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

JLLWrappers.@generate_main_file_header("Qt5Svg")
JLLWrappers.@generate_main_file("Qt5Svg", UUID("3af4ccab-a251-578e-a514-ea85a0ba79ee"))
end  # module Qt5Svg_jll


text:
import CompScienceMeshes.chart,CompScienceMeshes.subd_chart,CompScienceMeshes.GSubdMesh,CompScienceMeshes.subdMesh

mutable struct subdBasis{T,M,P} <: Space{T}
  geo::M
  fns::Vector{Vector{Shape{T}}}
  pos::Vector{P}
  data::Tuple{Array{subd_chart},Array{Array{Array{Tuple}}}}
end

function subset(s::S,I) where {S<:subdBasis}
    @warn "No parallel assembly supported when using subdivision bases!"
    return s
end

mutable struct subReferenceSpace{T,D} <: RefSpace{T,D} end

refspace(s::subdBasis{T,M,P}) where {T,M,P} = subReferenceSpace{T,12}()
# 12 only for regular case

function (rs::subReferenceSpace)(nbhd)
    # u = parametric(nbhd)
    # x = cartesian(nbhd)
    s = shapefuns(nbhd)
    scurl = get_shape_curl(nbhd)
    # TODO: this needs fixing, stat!
    return s #,scurl
end

function subdsurface(mesh)
    subd_mesh = GSubdMesh(mesh)
    vertices = mesh.vertices
    subd_elements = subd_mesh.elements
    nvertices = length(vertices)
    nelem = length(subd_elements)
    funs=Array{Vector{BEAST.Shape{Float64}}}(undef,nvertices)
    for i = 1 : nvertices funs[i] = [] end
    for ie = 1:nelem
        inodes = subd_elements[ie].RingNodes
        N = length(inodes)
        for ib = 1:N
            nodeid = inodes[ib]
            coeff = 1.0
            ifun = BEAST.Shape{Float64}(ie,ib,coeff)
            push!(funs[nodeid], ifun)
        end
    end
    data = assembly(subd_mesh)
    return subdBasis(mesh,funs,vertices,data)
end

function assemblydata(subdB::subdBasis)
    return subdB.data
end

function assembly(subdG::subdMesh)
    nelem = length(subdG.elements)
    assemblydata = Array{Array{Array{Tuple}}}(undef,nelem)
    ElementCharts = Array{subd_chart}(undef,nelem)
    for e = 1:nelem
        cha = chart(subdG,e)
        ringnodes = cha.RingNodes
        nnodes = length(ringnodes)
        line=Array{Array{Tuple}}(undef,nnodes)
        for inode = 1:nnodes
            Node = ringnodes[inode]
            entries=[]
            entry=tuple(Node,1.0)
            push!(entries,entry)
            line[inode] = entries
        end
        assemblydata[e] = line
        ElementCharts[e] = cha
    end
    return(ElementCharts,assemblydata)
end


text:
"""
    FrameBuffer{T}

A buffer for images, this type supports efficient updating and calculation of statistics along it's third dimension.

Construct using
    fb = FrameBuffer(img::Matrix, n)
    fb = FrameBuffer{T}(w::Int,h::Int,n::Int)

where `n` is the buffer capacity. `fb` supports `push!, median, mean, sum, var, std, reshape, size, length, capacity, Matrix, isready` and iteration.
"""
mutable struct FrameBuffer{T}
    b::Array{T,3}
    c::Int
    full::Bool
    function FrameBuffer(img::Matrix, n)
        fb = new{eltype(img)}(similar(img, size(img)..., n), 0, false)
        push!(fb, img)
        fb
    end
    function FrameBuffer(buffer::Array{T,3}) where T
        fb = new{T}(buffer, 0, true)
    end
    function FrameBuffer{T}(w::Int,h::Int,d::Int) where T
        new{T}(Array{T,3}(undef,w,h,d), 0, false)
    end
end

capacity(fb::FrameBuffer) = size(fb.b,3)

function Base.push!(b::FrameBuffer, img)
    b.c += 1
    if b.c > capacity(b)
        b.full = true
        b.c = 1
    end
    b.b[:,:,b.c] .= img
end

Base.@propagate_inbounds function Base.getindex(b::FrameBuffer, i::Int)
    Base.@boundscheck if !b.full && i > b.c
        throw(BoundsError(b,i))
    end
    @view b.b[:,:,i]
end

Base.@propagate_inbounds function  Base.getindex(b::FrameBuffer, i,j,k)
    Base.@boundscheck if !b.full && i > b.c
        throw(BoundsError(b,i))
    end
    b.b[i,j,k]
end

Base.Matrix(b::FrameBuffer) = reshape(b.b, :, length(b))
Base.length(b::FrameBuffer) = b.full ? size(b.b, 3) : b.c
imgsize(b::FrameBuffer) = (size(b.b,1),size(b.b,2))

Base.isready(b::FrameBuffer) = b.c > 0
function Base.iterate(b::FrameBuffer, state=1)
    b.full && state > length(b) && return nothing
    !b.full && state > b.c && return nothing
    (b[state],state+1)
end

for f in (median, mean, sum, std, var, reshape, size)
    m = parentmodule(f)
    fs = nameof(f)
    @eval function $m.$fs(b::FrameBuffer{T}, args...)::Matrix{T} where T
        if b.full
            # return dropdims($fs(b.b, args..., dims=3), dims=3)
            return map($fs, Slices(b.b, 3))
        else
            return dropdims($fs(@view(b.b[:,:,1:b.c]), args..., dims=3), dims=3)
        end
    end
end

for f in (diff,)
    m = parentmodule(f)
    fs = nameof(f)
    @eval function $m.$fs(b::FrameBuffer{T}, args...) where T
        if b.full
            return FrameBuffer($fs(b.b, args..., dims=3))
            # return map($fs, Slices(b.b, 3))
        else
            return FrameBuffer($fs(@view(b.b[:,:,1:b.c]), args..., dims=3))
        end
    end
end


text:
using JuMP
using MathOptInterface

const MOI = MathOptInterface

struct Solution
    nodes::Vector{Node}
    lsol::Array{Float64,2}
    πsol::Vector{Float64}
end

struct Result
    status::MOI.TerminationStatusCode
    sol::Solution
    value::Float64
end

include("nlp.jl")
include("soc.jl")


text:
function get_conventions(obj::ADMParser)
    return jcall(obj, "getConventions", IERSConventions, ())
end

function get_data_context(obj::ADMParser)
    return jcall(obj, "getDataContext", DataContext, ())
end

function get_launch_number(obj::ADMParser)
    return jcall(obj, "getLaunchNumber", jint, ())
end

function get_launch_piece(obj::ADMParser)
    return jcall(obj, "getLaunchPiece", JString, ())
end

function get_launch_year(obj::ADMParser)
    return jcall(obj, "getLaunchYear", jint, ())
end

function get_mission_reference_date(obj::ADMParser)
    return jcall(obj, "getMissionReferenceDate", AbsoluteDate, ())
end

function get_mu(obj::ADMParser)
    return jcall(obj, "getMu", jdouble, ())
end

function is_simple_eop(obj::ADMParser)
    return jcall(obj, "isSimpleEOP", jboolean, ())
end

function parse(obj::ADMParser, arg0::InputStream)
    return jcall(obj, "parse", ADMFile, (InputStream,), arg0)
end

function parse(obj::ADMParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "parse", ADMFile, (InputStream, JString), arg0, arg1)
end

function parse(obj::ADMParser, arg0::JString)
    return jcall(obj, "parse", ADMFile, (JString,), arg0)
end

function with_conventions(obj::ADMParser, arg0::IERSConventions)
    return jcall(obj, "withConventions", ADMParser, (IERSConventions,), arg0)
end

function with_data_context(obj::ADMParser, arg0::DataContext)
    return jcall(obj, "withDataContext", ADMParser, (DataContext,), arg0)
end

function with_international_designator(obj::ADMParser, arg0::jint, arg1::jint, arg2::JString)
    return jcall(obj, "withInternationalDesignator", ADMParser, (jint, jint, JString), arg0, arg1, arg2)
end

function with_mission_reference_date(obj::ADMParser, arg0::AbsoluteDate)
    return jcall(obj, "withMissionReferenceDate", ADMParser, (AbsoluteDate,), arg0)
end

function with_mu(obj::ADMParser, arg0::jdouble)
    return jcall(obj, "withMu", ADMParser, (jdouble,), arg0)
end

function with_simple_eop(obj::ADMParser, arg0::jboolean)
    return jcall(obj, "withSimpleEOP", ADMParser, (jboolean,), arg0)
end


text:
using Yao
using QuDiffEq
using LinearAlgebra
using BitBasis
using Random
using Test
using YaoBlocks
using OrdinaryDiffEq

#Linear Diff Equation Unitary M
function f(du,u,p,t)
    du[1] = -3*u[1]^2 + u[2]
    du[2] = -u[2]^2 - u[1]*u[2]
end
@testset "QuNLDE_Test" begin
    Random.seed!(4)
    N = 2
    k = 3
    siz = nextpow(2, N + 1)
    x = normalize!(rand(N))

    A = zeros(ComplexF32,2^(siz),2^(siz))
    A[1,1] = ComplexF32(1)
    A[5,3] = ComplexF32(1)
    A[5,6] = ComplexF32(-3)
    A[9,11] = ComplexF32(-1)
    A[9,7] = ComplexF32(-1)
    tspan = (0.0,0.4)
    qprob = QuLDEProblem(A, x, tspan)
    r, N = func_transform(qprob.A, qprob.b, k)
    out = N*vec(state(r))
    r_out = zero(x)
    f(r_out, x,1,1)
    @test isapprox.(r_out, out[2:3]*sqrt(2), atol = 1e-3) |> all

    prob = ODEProblem(f, x, tspan)
    sol = solve(prob, Euler(), dt = 0.1, adaptive = false)
    r_out = transpose(hcat(sol.u...))
    out = solve(qprob, QuNLDE(3), dt = 0.1)
    @test isapprox.(r_out,real(out), atol = 1e-3) |> all
end


text:
""" Functorial data migration for attributed C-sets.
"""
module DataMigrations
export DataMigration, DeltaMigration, SigmaMigration, migrate, migrate!,
  representable, yoneda, colimit_representables

using ...Syntax, ...Present, ...Theories
using ...Theories: SchemaDesc, ob, hom, dom, codom, attr, adom
using ..Categories, ..FinCats, ..Limits, ..Diagrams, ..FinSets, ..CSets
using ...Graphs, ..FreeDiagrams
import ..Categories: ob_map, hom_map
using ..FinCats: make_map, mapvals

# Data types
############

""" Conjunctive query over schema ``C``.

The diagram comprising the query specifies a finite limit.
"""
const ConjQuery{C<:FinCat} = Diagram{op,C}

""" Gluing or agglomerative query over schema ``C``.

The diagram comprising the query specifies a finite colimit. In the important
special case that the diagram has discrete shape, it specifies a finite
coproduct and the query is called "linear" or "disjunctive".
"""
const GlueQuery{C<:FinCat} = Diagram{id,C}

""" "Gluc query" (gluing of conjunctive queries) over schema ``C``.

The diagram of diagrams comprising the query specifies a finite colimit of
finite limits. In the important special case that the outer diagram has discrete
shape, it specifies a finite coproduct of finite limits and the query is called
a "duc query" (disjoint union of conjunctive queries).
"""
const GlucQuery{C<:FinCat} = Diagram{id,<:TypeCat{<:Diagram{op,C}}}

""" Functor defining a pullback or delta data migration.
"""
const DeltaSchemaMigration{D<:FinCat,C<:FinCat} = FinFunctor{D,C}

""" Functor defining a contravariant data migration using conjunctive queries.
"""
const ConjSchemaMigration{D<:FinCat,C<:FinCat} =
  FinDomFunctor{D,<:TypeCat{<:ConjQuery{C}}}

""" Functor defining a contravariant data migration using gluing queries.
"""
const GlueSchemaMigration{D<:FinCat,C<:FinCat} =
  FinDomFunctor{D,<:TypeCat{<:GlueQuery{C}}}

""" Functor defining a contravariant data migration using gluc queries.
"""
const GlucSchemaMigration{D<:FinCat,C<:FinCat} =
  FinDomFunctor{D,<:TypeCat{<:GlucQuery{C}}}

""" Abstract type for a data migration functor.
"""
abstract type MigrationFunctor{Dom<:ACSet,Codom<:ACSet} <:
  Functor{TypeCat{Dom,ACSetTransformation},TypeCat{Codom,ACSetTransformation}} end

ob_map(F::MigrationFunctor{Dom,Codom}, X) where {Dom,Codom} =
  ob_map(F, Codom, X)

(F::MigrationFunctor)(X::ACSet) = ob_map(F, X)
(F::MigrationFunctor)(α::ACSetTransformation) = hom_map(F, α)

""" Data migration functor given contravariantly.

This type encompasses data migration functors from ``C``-sets to ``D``-sets
given contravariantly by a functor out of the schema ``D``. The simplest such
functor is a pullback data migration ([`DeltaMigration`](@ref)), specified by a
functor ``D → C`` between the schemas. Other important cases include conjunctive
and duc data migrations.
"""
struct DataMigration{Dom,Codom,F<:FinDomFunctor} <: MigrationFunctor{Dom,Codom}
  functor::F
end

DataMigration(functor::F, ::Type{Dom}, ::Type{Codom}) where {F,Dom,Codom} =
  DataMigration{Dom,Codom,F}(functor)
DataMigration(functor::F, ::Type{Codom}) where {F,Codom} =
  DataMigration{ACSet,Codom,F}(functor)

ob_map(F::DataMigration, T::Type, X) = migrate(T, X, F.functor)

""" Delta or pullback functorial data migration between acsets.

Given a functor ``F: D → C``, the delta migration functor ``Δ_F`` from
``C``-sets to ``D``-sets is defined contravariantly by ``Δ_F(X) := X ∘ F``.

See (Spivak, 2014, *Category Theory for the Sciences*) for details.
"""
const DeltaMigration{Dom,Codom} = DataMigration{Dom,Codom,<:DeltaSchemaMigration}

DeltaMigration(args...) = DataMigration(args...)::DeltaMigration

const ConjMigration{Dom,Codom} = DataMigration{Dom,Codom,<:ConjSchemaMigration}
const GlueMigration{Dom,Codom} = DataMigration{Dom,Codom,<:GlueSchemaMigration}
const GlucMigration{Dom,Codom} = DataMigration{Dom,Codom,<:GlucSchemaMigration}

# Contravariant migration
#########################

""" Contravariantly migrate data from the acset `Y` to a new acset of type `T`.

The mutating variant of this function is [`migrate!`](@ref).
"""
function migrate(::Type{T}, X::ACSet, F::FinDomFunctor; kw...) where T <: ACSet
  T(migrate(X, F; kw...))
end
function migrate(X::ACSet, F::FinDomFunctor; kw...)
  migrate(FinDomFunctor(X), F; kw...)
end

""" Contravariantly migrate data from the acset `Y` to the acset `X`.

This is the mutating variant of [`migrate!`](@ref). When the functor on schemas
is the identity, this operation is equivalent to [`copy_parts!`](@ref).
"""
function migrate!(X::ACSet, Y::ACSet, F::FinDomFunctor; kw...)
  copy_parts!(X, migrate(Y, F; kw...))
end

# Delta migration
#----------------

migrate(::Type{T}, X::ACSet, F::DeltaSchemaMigration) where T <: ACSet =
  migrate!(T(), X, F)
migrate(::Type{T}, X::ACSet, FOb, FHom) where T <: ACSet =
  migrate!(T(), X, FOb, FHom)

function migrate!(X::StructACSet{S}, Y::ACSet, F::DeltaSchemaMigration) where S
  partsX = Dict(c => add_parts!(X, c, nparts(Y, nameof(ob_map(F,c))))
                for c in ob(S))
  for (f,c,d) in zip(hom(S), dom(S), codom(S))
    set_subpart!(X, partsX[c], f, partsX[d][subpart(Y, hom_map(F,f))])
  end
  for (f,c) in zip(attr(S), adom(S))
    set_subpart!(X, partsX[c], f, subpart(Y, hom_map(F,f)))
  end
  return X
end

function migrate!(X::ACSet, Y::ACSet, FOb, FHom)
  F = FinFunctor(FOb, FHom, FinCat(Presentation(X)), FinCat(Presentation(Y)))
  migrate!(X, Y, F)
end

function (::Type{T})(X::ACSet, F::FinFunctor) where T <: ACSet
  Base.depwarn("Data migration constructor is deprecated, use `migrate` instead", :ACSet)
  migrate(T, X, F)
end

function (::Type{T})(X::ACSet, FOb::AbstractDict, FHom::AbstractDict) where T <: ACSet
  Base.depwarn("Data migration constructor is deprecated, use `migrate` instead", :ACSet)
  migrate(T, X, FOb, FHom)
end

# Conjunctive migration
#----------------------

function migrate(X::FinDomFunctor, F::ConjSchemaMigration;
                 return_limits::Bool=false, tabular::Bool=false)
  tgt_schema = dom(F)
  limits = make_map(ob_generators(tgt_schema)) do c
    Fc = ob_map(F, c)
    # XXX: Disable domain check because acsets don't store schema equations.
    lim = limit(compose(Fc, X, strict=false),
                alg=SpecializeLimit(fallback=ToBipartiteLimit()))
    if tabular
      J = shape(Fc)
      names = (ob_generator_name(J, j) for j in ob_generators(J))
      TabularLimit(lim, names=names)
    else
      lim
    end
  end
  funcs = make_map(hom_generators(tgt_schema)) do f
    Ff, c, d = hom_map(F, f), dom(tgt_schema, f), codom(tgt_schema, f)
    # XXX: Disable domain check for same reason.
    universal(compose(Ff, X, strict=false), limits[c], limits[d])
  end
  Y = FinDomFunctor(mapvals(ob, limits), funcs, tgt_schema)
  return_limits ? (Y, limits) : Y
end

# Gluing migration
#-----------------

function migrate(X::FinDomFunctor, F::GlueSchemaMigration)
  tgt_schema = dom(F)
  colimits = make_map(ob_generators(tgt_schema)) do c
    Fc = ob_map(F, c)
    # XXX: Force composition to tighten the codomain types.
    colimit(force(compose(Fc, X, strict=false)), alg=SpecializeColimit())
  end
  funcs = make_map(hom_generators(tgt_schema)) do f
    Ff, c, d = hom_map(F, f), dom(tgt_schema, f), codom(tgt_schema, f)
    universal(compose(Ff, X, strict=false), colimits[c], colimits[d])
  end
  FinDomFunctor(mapvals(ob, colimits), funcs, tgt_schema)
end

# Gluc migration
#---------------

function migrate(X::FinDomFunctor, F::GlucSchemaMigration)
  tgt_schema = dom(F)
  colimits_of_limits = make_map(ob_generators(tgt_schema)) do c
    Fc = ob_map(F, c)
    Fc_set, limits = migrate(X, diagram(Fc), return_limits=true)
    (colimit(Fc_set), Fc_set, limits)
  end
  funcs = make_map(hom_generators(tgt_schema)) do f
    Ff, c, d = hom_map(F, f), dom(tgt_schema, f), codom(tgt_schema, f)
    Fc_colim, Fc_set, Fc_limits = colimits_of_limits[c]
    Fd_colim, Fd_set, Fd_limits = colimits_of_limits[d]
    component_funcs = map(ob_generators(dom(Fc_set))) do j
      j′, Ffⱼ = ob_map(Ff, j)
      universal(compose(Ffⱼ, X, strict=false), Fc_limits[j], Fd_limits[j′])
    end
    Ff_set = DiagramHom{id}(shape_map(Ff), component_funcs, Fc_set, Fd_set)
    universal(Ff_set, Fc_colim, Fd_colim)
  end
  FinDomFunctor(mapvals(ob∘first, colimits_of_limits), funcs, tgt_schema)
end

# Sigma migration
#################

""" Sigma or left pushforward functorial data migration between acsets.

Given a functor ``F: C → D``, the sigma data migration ``Σ_F`` is a functor from
``C``-sets to ``D``-sets that is left adjoint to the delta migration functor
``Δ_F`` ([`DeltaMigration`](@ref)). Explicitly, the ``D``-set ``Σ_F(X)`` is
given on objects ``d ∈ D`` by the formula ``Σ_F(x)(d) = \\mathrm{colim}_{F ↓ d}
X ∘ π``, where ``π: (F ↓ d) → C`` is the projection.

See (Spivak, 2014, *Category Theory for the Sciences*) for details.
"""
struct SigmaMigration{Dom,Codom,F<:FinFunctor,CC} <: MigrationFunctor{Dom,Codom}
  functor::F
  comma_cats::CC

  function SigmaMigration(functor::F, ::Type{Dom}, ::Type{Codom}) where
      {F<:FinFunctor, Dom, Codom}
    comma_cats = get_comma_cats(functor)
    new{Dom,Codom,F,typeof(comma_cats)}(functor, comma_cats)
  end
end

SigmaMigration(functor::FinFunctor, ::Type{Codom}) where Codom =
  SigmaMigration(functor, ACSet, Codom)

ob_map(ΣF::SigmaMigration, ::Type{T}, X::ACSet) where T<:ACSet =
  ob_map(ΣF, T, FinDomFunctor(X))

function ob_map(ΣF::SigmaMigration, ::Type{T}, X::FinDomFunctor) where T<:ACSet
  comma_cats, comma_hom_map = ΣF.comma_cats
  diagramD = FreeDiagram(presentation(codom(ΣF.functor)))

  # define Y on objects by taking colimits
  Y = T()
  colimX = map(parts(diagramD, :V)) do i
    F∇d = ob(comma_cats, i)
    Xobs = FinSet{Int,Int}[ ob_map(X, c) for (c,_) in ob(F∇d) ]
    Xhoms = [ hom_map(X, hom(F∇d, g)) for g in parts(F∇d, :E) ]
    colimit(FreeDiagram(Xobs, collect(zip(Xhoms, src(F∇d), tgt(F∇d)))))
  end

  for d in parts(diagramD, :V)
    add_parts!(Y, nameof(ob(diagramD, d)), length(apex(colimX[d])))
  end

  # Define Y on morphisms by the universal property
  for g in parts(diagramD, :E)
    if nparts(Y, nameof(dom(hom(diagramD, g)))) == 0
      continue
    end
    src_colim, tgt_colim = colimX[src(diagramD, g)], colimX[tgt(diagramD, g)]
    ϕ = hom(comma_cats, comma_hom_map[g])
    set_subpart!(Y, nameof(hom(diagramD, g)),
      collect(universal(src_colim,
        Multicospan(apex(tgt_colim), legs(tgt_colim)[collect(ϕ[:V])])
      )))
  end
  return Y
end

""" add_hom!(d::FreeDiagram{Ob, Hom}, src_ob::Ob, tgt_ob::Ob, hom::Hom)

Adds a hom to `d` from the vertex with object `src_ob` to the vertex with object `tgt_ob`.
"""
function add_hom!(d::FreeDiagram, src_ob, tgt_ob, hom)
  src_idx = first(incident(d, src_ob, :ob))
  tgt_idx = first(incident(d, tgt_ob, :ob))
  return add_edge!(d, src_idx, tgt_idx, hom = hom)
end

"""   comma_cats(diagramD::FreeDiagram{FreeSchema.Ob, FreeSchema.Hom}, FOb, FHom)

Given a free category ``\\mathcal{D}`` with no cycles and
a functor represented by the pair `(FOb, FHom)`, returns a diagram
``\\mathcal{D} \\to \\mathsf{Cat}`` which on objects takes ``d \\in \\mathcal{D}`` to the
comma category ``F \\downarrow d``.
"""
function get_comma_cats(F::FinFunctor)
  diagramD = FreeDiagram(presentation(codom(F)))
  FObInv = y -> filter(x -> ob_map(F,x) == y, ob_generators(dom(F)))
  FHomInv = g -> filter(f -> hom_map(F,f) == g, hom_generators(dom(F)))
  comma_cats = FreeDiagram{FreeDiagram, ACSetTransformation}()
  add_vertices!(comma_cats,
    nparts(diagramD, :V),
    ob = map(parts(diagramD, :V)) do _
      FreeDiagram{Pair{FreeSchema.Ob, FreeSchema.Hom}, FreeSchema.Hom}()
    end
  )

  comma_hom_map = Dict{Int,Int}()
  for d in topological_sort(diagramD)
    F∇d = ob(comma_cats, d)
    id_d = id(ob(diagramD, d))

    # If FOb[c] = d, add an objects (c, id(d)) to F∇d
    preimage_d = FObInv(ob(diagramD, d))
    id_obs = add_vertices!(F∇d, length(preimage_d), ob = map(c->c=>id_d, preimage_d))

    # if FHom[h] = id(d), add a hom h: (dom(h), id(d)) -> (codom(h), id(d)) to F∇d
    id_homs = map(FHomInv(id_d)) do h
      add_hom!(F∇d, dom(h) => id_d, codom(h) => id_d, h)
    end

    # for g: d' -> d in D (note that F∇d' must have already been constructed)
    #     populate F∇d with obs and homs coming from F∇d′
    for g in incident(diagramD, d, :tgt)
      d′ = src(diagramD, g)
      F∇g = comma_cat_hom!(F∇d, ob(comma_cats, d′), id_d, hom(diagramD, g), FHomInv)
      comma_hom_map[g] = add_edge!(comma_cats, d′, d, hom=F∇g)
    end
  end

  return comma_cats, comma_hom_map
end

function comma_cat_hom!(F∇d, F∇d′, id_d, g, FHomInv)
  # for (c,f) in F∇d' add ob (c, compose(f,g)) to F∇d
  vs = add_vertices!(F∇d, nparts(F∇d′, :V), ob = map(ob(F∇d′)) do (c,f)
    c => compose(f, g)
  end)

  # for h: (c, f) -> (c', f') in diagram_d', add hom h to F∇d
  es = add_edges!(F∇d, vs[src(F∇d′)], vs[tgt(F∇d′)], hom = hom(F∇d′))

  # for every (c,f) in the recently added objects. If FHom[h] == f, then add the hom
  #       h: (c,f) -> (codom(h), idv)
  # relies on D being free
  for (c, f) in ob(F∇d, vs)
    for h in FHomInv(f)
      dom(h) == c && add_hom!(F∇d, c => f, codom(h) => id_d, h)
    end
  end

  # return the inclusion from F∇d′ to F∇d
  return ACSetTransformation((V = collect(vs), E = collect(es)), F∇d′, F∇d)
end

# Yoneda embedding
#-----------------

""" Construct a representable C-set.

Recall that a *representable* C-set is one of form ``C(c,-): C → Set`` for some
object ``c ∈ C``.

This function computes the ``c`` representable as the left pushforward data
migration of the singleton ``{c}``-set along the inclusion functor ``{c} ↪ C``,
which works because left Kan extensions take representables to representables
(Mac Lane 1978, Exercise X.3.2). Besides the intrinsic difficulties with
representables (they can be infinite), this function thus inherits any
limitations of our implementation of left pushforward data migrations.
"""
function representable(::Type{T}, C::Presentation{Schema}, ob::Symbol) where T <: ACSet
  C₀ = Presentation{Symbol}(FreeSchema)
  add_generator!(C₀, C[ob])
  F = FinFunctor(Dict(ob => ob), Dict(), C₀, C)
  ΣF = SigmaMigration(F, T)

  X = FinDomFunctor(Dict(ob => FinSet(1)),
                    Dict{Symbol,FinFunction{Int}}(), FinCat(C₀))
  ob_map(ΣF, X)
end
representable(::Type{T}, ob::Symbol) where T <: StructACSet =
  representable(T, Presentation(T), ob)

""" Yoneda embedding of category C in category of C-sets.

Because Catlab privileges copresheaves (C-sets) over presheaves, this is the
*contravariant* Yoneda embedding, i.e., the embedding C^op → C-Set.

Returns a `FinDomFunctor` with domain `op(C)`.
"""
function yoneda(::Type{T}, C::Presentation{Schema}) where T <: ACSet
  y_ob = Dict(c => representable(T, C, nameof(c)) for c in generators(C, :Ob))
  y_hom = Dict(Iterators.map(generators(C, :Hom)) do f
    c, d = dom(f), codom(f)
    yc, yd = y_ob[c], y_ob[d]
    initial = Dict(nameof(d) => Dict(1 => yc[1,f]))
    f => homomorphism(yd, yc, initial=initial) # Unique homomorphism.
  end)
  FinDomFunctor(y_ob, y_hom, op(FinCat(C)))
end
yoneda(::Type{T}) where T <: StructACSet = yoneda(T, Presentation(T))

""" Interpret conjunctive data migration as a colimit of representables.

Given a conjunctive data migration (a functor `J → Diag{op}(C)`) and the Yoneda
embedding for `C` (a functor `op(C) → C-Set` computed via [`yoneda`](@ref)),
take colimits of representables to construct a `op(J)`-shaped diagram of C-sets.

Since every C-set is a colimit of representables, this is a generic way of
constructing diagrams of C-sets.
"""
function colimit_representables(F::DeltaSchemaMigration, y)
  compose(op(F), y)
end
function colimit_representables(F::ConjSchemaMigration, y)
  C = dom(F)
  colimits = make_map(ob_generators(C)) do c
    Fc = ob_map(F, c)
    colimit(compose(op(Fc), y))
  end
  homs = make_map(hom_generators(C)) do f
    Ff, c, d = hom_map(F, f), dom(C, f), codom(C, f)
    universal(compose(op(Ff), y), colimits[d], colimits[c])
  end
  FinDomFunctor(mapvals(ob, colimits), homs, op(C))
end

# Schema translation
####################

# FIXME: This function does not belong here.

"""   FreeDiagram(pres::Presentation{FreeSchema, Symbol})

Returns a `FreeDiagram` whose objects are the generating objects of `pres` and
whose homs are the generating homs of `pres`.
"""
function FreeDiagrams.FreeDiagram(pres::Presentation{Schema, Symbol}) where Schema
  obs = Array{FreeSchema.Ob}(generators(pres, :Ob))
  homs = Array{FreeSchema.Hom}(generators(pres, :Hom))
  doms = map(h -> generator_index(pres, nameof(dom(h))), homs)
  codoms = map(h -> generator_index(pres, nameof(codom(h))), homs)
  return FreeDiagram(obs, collect(zip(homs, doms, codoms)))
end

end

max_forks_repo_name: epatters/CompCat


text:
"""
    TextWrap

This module provides the function [`wrap`](@ref) which parses an input text and reorganizes its
white space so that it can be printed with a fixed screen width, optionally indenting it. It also
provides the two convenience functions [`print_wrapped`](@ref) and [`println_wrapped`](@ref).
"""
module TextWrap

export
    wrap,
    print_wrapped,
    println_wrapped

ansi_length(s) = length(replace(s, r"\e\[[0-9]+(?:;[0-9]+)*m" => ""))

function apply_expand_tabs(text::AbstractString, i0::Int)
    out_buf = IOBuffer()
    i = i0 % 8
    for c in text
        if c == '\t'
            while i < 8
                print(out_buf, ' ')
                i += 1
            end
            i = 0
        elseif c == '\n' || c == '\r'
            print(out_buf, c)
            i = 0
        else
            print(out_buf, c)
            i = (i + 1) % 8
        end
    end
    return String(take!(out_buf))
end

function check_width(width::Integer)
    width ≤ 0 && throw(ArgumentError("invalid width $width (must be > 0)"))
    return true
end
function check_indent(indent::Integer, width::Integer)
    0 ≤ indent < width ||
        throw(ArgumentError("invalid intent $indent (must be an integer between 0 and width-1, " *
                            "or an AbstractString)"))
    return true
end
function check_indent(indent::AbstractString, width::Integer)
    length(indent) ≥ width && throw(ArgumentError("invalid intent (must be shorter than width-1)"))
    return true
end

mutable struct State
    cln::Int    # current line number
    cll::Int    # current line length
    bol::Bool   # beginning of line?
    lcise::Bool # last chunk is sentence ending
    soh::String # space on hold

    State() = new(1, 0, true, false, "")
end

struct Params
    width::Int
    iind::String    # initial indent
    sind::String    # subsequent indent
    exp_tabs::Bool  # expand tabs
    rep_white::Bool # replace whitespace
    fix_end::Bool   # fix sentence endings
    brk_hyp::Bool   # break on hyphens
    brk_long::Bool  # break long words
    rec_esc::Bool   # recognize escapes

    function Params(width::Integer,
                    initial_indent::Union{Integer,AbstractString},
                    subsequent_indent::Union{Integer,AbstractString},
                    expand_tabs::Bool,
                    replace_whitespace::Bool,
                    fix_sentence_endings::Bool,
                    break_on_hyphens::Bool,
                    break_long_words::Bool,
                    recognize_escapes::Bool
                    )
        check_width(width)
        check_indent(initial_indent, width)
        check_indent(subsequent_indent, width)

        iind::String = initial_indent isa Integer ? " "^initial_indent : initial_indent
        sind::String = subsequent_indent isa Integer ? " "^subsequent_indent : subsequent_indent

        return new(width, iind, sind, expand_tabs, replace_whitespace, fix_sentence_endings,
                   break_on_hyphens, break_long_words, recognize_escapes)
    end
end

# This function just performs breaks-on-hyphens and passes individual chunks to put_chunk!
function put_chunks!(out_str::IOBuffer, chunk::AbstractString,
                     s::State, p::Params)

    hyphen_re = r"# define a class that matches sequences of word characters
                  # and escape codes, arbitrarily mixed. It's then invoked
                  # with the syntax (?&w)
                  (?(DEFINE)
                      (?<w> (?:\w|\e\[[0-9]+m)*)
                  )

                  # breakdown: 1) possible prefix (can consist of number-only words followed
                  #               by a dash, possibly more than one);
                  #            2) main body: requires at least a letter, includes the dash
                  #            3) rest of the word: also requires at least a letter
                  # notes: the ?: avoids group capturing
                  #        the ?> avoids backtracking as soon as a \p{N} or \p{L} is found
                  #        the ?= is a lookahead

                  ^(?:(?>(?&w)\p{N})(?&w)-)*?   # possible prefix
                      (?>(?&w)\p{L})(?&w)-      # main body
                   (?=(?>(?&w)\p{L})(?&w) )     # rest of the word
                 "x

    while p.brk_hyp
        m = match(hyphen_re, chunk)
        m ≡ nothing && break
        c = m.match
        put_chunk!(out_str, c, s, p)
        s.soh = ""
        chunk = chunk[m.offset+lastindex(c):end]
    end

    put_chunk!(out_str, chunk, s, p)
end

# Writes a chunk to out_str, based on the current state, and updates the state.
# Besides the current position, encoded in (cln, cll, bol) = (current_line_number,
# current_line_length, beginning_of_line), which gets updated, it also sets a flag
# to signal that an end-of-sentence was detected.
# The field s.soh (=space_on_hold) is the spacing which should go in front of chunk,
# and it may be discarded when we're between lines.
function put_chunk!(out_str::IOBuffer, chunk::AbstractString,
                    s::State, p::Params)

    # This is written as a new function rather than a function reference
    # to help type inference
    elength(x)::Int = p.rec_esc ? ansi_length(x) : length(x)

    liind = elength(p.iind)
    lsind = elength(p.sind)
    lchunk = elength(chunk)
    lsoh = length(s.soh)

    if s.cll + lsoh > p.width
        s.soh = ""
        lsoh = 0
        s.cll > 0 && print(out_str, "\n")
        s.cln += 1
        s.cll = 0
        s.bol = true
    end

    if s.bol
        ind, lind = s.cln == 1 ? (p.iind, liind) : (p.sind, lsind)
        print(out_str, ind)
        s.cll += lind
        s.soh = ""
        lsoh = 0
    end

    # Is there enough room for the chunk? or is this the
    # beginning of the text and we cannot break words?
    if s.cll + lsoh + lchunk ≤ p.width || (s.cln == 1 && s.bol && !p.brk_long)
        print(out_str, s.soh, chunk)
        s.cll += lchunk + lsoh
        s.bol = false
    # Does the chunk fit into the next line? or are we
    # forced to put it there?
    elseif lchunk ≤ p.width - lsind || !p.brk_long
        print(out_str, s.bol ? "" : "\n", p.sind, chunk)
        s.cll = lsind + lchunk
        s.cln += 1
        s.bol = false
    # Break it until it fits
    else
        while s.cll + lsoh + lchunk > p.width
            if p.width - s.cll - lsoh > 0
                print(out_str, s.soh, chunk[1:nextind(chunk, 0, p.width-s.cll-lsoh)], "\n", p.sind)
                chunk = chunk[nextind(chunk, 0, p.width-s.cll-lsoh+1):end]
                lchunk = elength(chunk)
            else
                print(out_str, "\n", p.sind)
            end
            s.cll = lsind
            s.cln += 1
            s.soh = ""
            lsoh = 0
        end
        print(out_str, chunk)
        s.cll += lchunk
        s.bol = false
    end

    # Detect end-of-sentences
    s.lcise = occursin(r"\w([\.\!\?…]|\.\.\.)[\"\'´„]?\Z", chunk)

    return out_str, s
end

"""
    wrap(string; keywords...)

Parses `string` and returns a new string in which newlines are inserted as appropriate in order
for each line to fit within a specified width.

The behaviour can be controlled via optional keyword arguments:

* `width` (deafult=`70`): the maximum width of the wrapped text, including indentation.
* `initial_indent` (default=`""`): indentation of the first line. This can
   be any string (shorter than `width`), or it can be an integer number (smaller than `width`).
* `subsequent_indent` (default=`""`): indentation of all lines except the first. Works the same as
   `initial_indent`.
* `break_on_hyphens` (default=`true`): this flag determines whether words can be broken on hyphens,
  e.g. whether "high-precision" can be split into "high-" and "precision".
* `break_long_words` (default=`true`): this flag determines what to do when a word is too long to
  fit in any line. If `true`, the word will be broken, otherwise it will go beyond the desired text
  width.
* `replace_whitespace` (default=`true`): if this flag is `true`, all whitespace characters in the
  original text (including newlines) will be replaced by spaces. Otherwise, they'll be preserved,
  except at the beginning or end of a line.
* `expand_tabs` (default=`true`): if this flag is `true`, tabs will be expanded in-place into
  spaces. Otherwise a tab is counted as a single character. The expansion happens before whitespace
  replacement.
* `fix_sentence_endings` (default=`false`): if this flag is `true`, the wrapper will try to
  recognize sentence endings in the middle of a paragraph and put two spaces before the next
  sentence in case only one is present.
* `recognize_escapes` (default=`true`): if `true`, compute all lengths ignoring ANSI escape codes
  (special character sequences used e.g. to modify the text color or other properties; they look
  e.g. like `\"\\e[94m\"`)
"""
function wrap(text::AbstractString;
              width::Int = 70,
              initial_indent::Union{Integer,AbstractString} = "",
              subsequent_indent::Union{Integer,AbstractString} = "",
              expand_tabs::Bool = true,
              replace_whitespace::Bool = true,
              fix_sentence_endings::Bool = false,
              break_long_words::Bool = true,
              break_on_hyphens::Bool = true,
              recognize_escapes::Bool = true)

    p = Params(width, initial_indent, subsequent_indent, expand_tabs, replace_whitespace,
               fix_sentence_endings, break_on_hyphens, break_long_words, recognize_escapes)

    # A regex to match any sequence of spaces except non-breakable spaces (\xA0)
    spaceregex = r"((?!\xA0)\s)+"

    # Whitespace-only case
    occursin(r"^\s*$", text) && return p.iind

    s = State()

    # We iterate over the text, looking for whitespace where to split.
    i = firstindex(text)
    l = lastindex(text)
    out_str = IOBuffer()

    wsrng = findnext(spaceregex, text, i)
    j, k = wsrng ≢ nothing ?
        (first(wsrng), nextind(text, last(wsrng))) :
        (0, -1)

    while 0 < j ≤ l
        if i < k
            if i < j
                # This is non-whitespace. We write it out according
                # to the current cursor position and the leading space.
                chunk = text[i:prevind(text,j)]

                put_chunks!(out_str, chunk, s, p)
            end
            i = k
        end

        # This is whitespace. We mangle it (expand tabs, fix
        # sentence endings, replace it with single spaces) and
        # then we keep it on hold.
        soh = text[j:prevind(text,k)]
        @assert !isempty(soh)
        if expand_tabs && occursin(r"\t", soh)
            soh = apply_expand_tabs(soh, s.cll)
        end
        if p.fix_end && s.lcise && soh == " "
            soh = "  "
        end
        if p.rep_white
            soh = replace(soh, r"\r?\n"=>"")
            soh = isempty(soh) ? " " : " "^length(soh)
        else
            while (fnl = findfirst(==('\n'), soh)) ≢ nothing # using the `==` form to support julia 1.0-1.2
                ind = s.cln == 1 ? p.iind : p.sind
                print(out_str, ind, '\n')
                s.cll = 0
                s.cln += 1
                s.bol = true
                soh = soh[nextind(soh, fnl):end]
            end
        end
        s.soh = soh

        # Continue the search
        k ≤ j && (k = nextind(text,j))
        wsrng = findnext(spaceregex, text, k)
        j, k = wsrng ≢ nothing ?
            (first(wsrng), nextind(text, last(wsrng))) :
            (0, -1)
    end
    if i ≤ ncodeunits(text)
        # Some non-whitespace is left at the end.
        chunk = text[i:end]
        put_chunks!(out_str, chunk, s, p)
    end
    return String(take!(out_str))
end

"""
    print_wrapped([io,] text...; keywords...)

This is just like the standard `print` function (it prints multiple arguments and accepts an
optional `IO` first argument), except that it wraps the result, and accepts keyword arguments to
pass to [`wrap`](@ref).
"""
print_wrapped(io::IO, text::AbstractString; kwargs...) = print(io, wrap(text; kwargs...))
print_wrapped(args...; kwargs...) = print_wrapped(stdout, args...; kwargs...)
print_wrapped(io::IO; kwargs...) = print_wrapped(io, ""; kwargs...)
print_wrapped(io::IO, args...; kwargs...) = print_wrapped(io, string(args...); kwargs...)

"""
    println_wrapped([io,] text...; keywords...)

Like [`print_wrapped`](@ref), but adds a newline at the end.
"""
function println_wrapped(io::IO, args...; kwargs...)
    print_wrapped(io, args...; kwargs...)
    print(io, "\n")
end
println_wrapped(args...; kwargs...) = println_wrapped(stdout, args...; kwargs...)

end # module TextWrap


text:
export product!, product, *, divide, /

"""
    product!(expod, expod1)

    Add the `expod1` degree to `expod` inplace (equivalent to the monomial product)
"""
function product!(expod::Exponent, expod1::Exponent)
    for (var, deg) in expod1
        if !haskey(expod, var)
            expod.expo[var] = Degree(0,0)
        end
        expod.expo[var].explvar += deg.explvar
        expod.expo[var].conjvar += deg.conjvar
        if (expod.expo[var].explvar, expod.expo[var].conjvar) == (0,0)
            delete!(expod, var)
        end
    end
    update_degree!(expod)
    return expod
end

function product!(expod::Exponent, var::Variable)
    if !haskey(expod, var)
        expod.expo[var] = Degree(0,0)
    end
    expod.expo[var].explvar += 1
    if (expod.expo[var].explvar, expod.expo[var].conjvar) == (0,0)
        delete!(expod, var)
    end

    update_degree!(expod)
    return expod
end

function product(exp1::Exponent, exp2::Exponent)
    expod = Exponent()
    product!(expod, exp1)
    product!(expod, exp2)
    return Exponent(expod)
end

## Polynomial
function product(p1::Polynomial, p2::Polynomial)
    p = Polynomial()
    for (expo1, λ1) in p1
        λ1 != 0 || continue
        for (expo2, λ2) in p2
            λ2 != 0 || continue
            expoprod = product(expo1, expo2)
            add_to_dict!(p.poly, expoprod, λ1 * λ2)
        end
    end
    update_degree!(p)
    return p
end

function product(p1::Polynomial, λ::Number)
    if λ == 0
        return Polynomial()
    end

    pdict = SortedDict{Exponent, Number}()
    for (expo, coeff) in p1
        pdict[expo] = coeff*λ
    end
    return Polynomial(pdict)
end

function product(p1::T, p2::U) where T<:AbstractPolynomial where U<:AbstractPolynomial
    return product(convert(Polynomial, p1), convert(Polynomial, p2))
end
# function product(p1::Number, p2::T) where T<:AbstractPolynomial
#     return product(convert(Polynomial, p1), convert(Polynomial, p2))
# end
# function product(p1::T, p2::Number) where T<:AbstractPolynomial
#     return product(convert(Polynomial, p1), convert(Polynomial, p2))
# end
function product(p1::T, p2::Number) where T<:AbstractPolynomial
    return product(convert(Polynomial, p1), p2)
end
function product(p1::Number, p2::T) where T<:AbstractPolynomial
    return product(convert(Polynomial, p2), p1)
end

function *(λ::Number, expo::Exponent)
    return Polynomial(SortedDict{Exponent, Number}(expo=>λ))
end
function *(expo::Exponent, λ::Number)
    return Polynomial(SortedDict{Exponent, Number}(expo=>λ))
end

function *(p1::T, p2::U) where T<:Union{Number, AbstractPolynomial} where U<:Union{Number, AbstractPolynomial}
    debug(LOGGER, "Inefficient implementation\n$(string(stacktrace()))")
    return product(p1, p2)
end

## Division
function divide(p1::Polynomial, p2::Polynomial)
    if length(p2) != 1
        error(LOGGER, "/(::Polynomial, ::Polynomial): Only allowed for monomial divisor ($(length(p2))-monomial polynomial here).")
    end
    expo, λ = collect(p2)[1]
    if expo.degree != Degree(0,0)
        error(LOGGER, "/(::Polynomial, ::Polynomial): Only allowed for constant divisor ($(expo.degree)-degree monomial here).")
    end
    if λ == 0
        error(LOGGER, "/(::Polynomial, ::Polynomial): Only allowed for non null constant divisor.")
    end
    return p1 * (1/λ)
end

function /(p1::T, p2::U) where T<:AbstractPolynomial where U<:AbstractPolynomial
    return divide(convert(Polynomial, p1), convert(Polynomial, p2))
end
function /(p1::Number, p2::T) where T<:AbstractPolynomial
    return divide(convert(Polynomial, p1), convert(Polynomial, p2))
end
function /(p1::T, p2::Number) where T<:AbstractPolynomial
    return divide(convert(Polynomial, p1), convert(Polynomial, p2))
end

## Point
function *(pt1::Point, λ::Number)
    debug(LOGGER, "Inefficient implementation\n$(string(stacktrace()))")
    pt = Point()
    if λ == 0
        return pt
    end
    for (var, val) in pt1
        pt[var] = λ*val
    end
    return pt
end
*(λ::Number, pt1::Point) = pt1*λ


text:
using Interact
using Colors
s=slider(0:0.1:1, label="Slider X:")
display(s)
c=map(x->RGB(x,0.5,0.5),observe(s))
display(c)
print(c)
using Images, ImageView
[c]

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
0.1
.1
1.
1e-1    # scientific notation
1e+10
1e-10
0x01p-1 # hex float

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
测试两个功能
"""

using Test
using MARY_fRG.Drawers
using MARY_fRG.Triangulated
using MARY_fRG.Fermi
using MARY_fRG.Fermi.Surface
using MARY_fRG.Fermi.Patch

#@testset "正方晶系" begin
#    quad = common_square_lattice(0.)
#    ltris, ladjs = split_square(quad.brillouin, 100)
#    pnum = 44
#    lpats = [find_patch_index_squa(tri.center, quad.brillouin, pnum) for tri in ltris]
#    sf, sfidx = const_energy_line_in_patches(ltris, ladjs, lpats, 0., quad.dispersion[1])
#    plt = draw_lines(quad.brillouin.edges)
#    draw_lines!(plt, sf; lcidx=sfidx)
#    sf2, sf2idx = const_energy_line_in_patches(ltris, ladjs, lpats, -0.2, quad.dispersion[1])
#    draw_lines!(plt, sf2; lcidx=sf2idx)
#    @savepng plt "surface1"
#end
#
#
#@testset "三角格子的费米面" begin
#    ts = common_triangle_lattice(0.)
#    ltris, ladjs = split_hexagon(ts.brillouin, 100)
#    pnum = 42
#    lpats = [find_patch_index_hexa(tri.center, ts.brillouin, pnum) for tri in ltris]
#    sf, sfidx = const_energy_line_in_patches(ltris, ladjs, lpats, 0., ts.dispersion[1])
#    plt = draw_lines(ts.brillouin.edges)
#    draw_lines!(plt, sf; lcidx=sfidx)
#    sf2, sf2idx = const_energy_line_in_patches(ltris, ladjs, lpats, 2.0, ts.dispersion[1])
#    draw_lines!(plt, sf2; lcidx=sf2idx)
#    @savepng plt "surface2"
#end

@testset "kagome格子的费米面" begin
    ts = upperband_kagome_lattice(-0.2)
    ltris, ladjs = split_hexagon(ts.brillouin, 100)
    pnum = 42
    lpats = [find_patch_index_hexa(tri.center, ts.brillouin, pnum) for tri in ltris]
    sf, sfidx = const_energy_line_in_patches(ltris, ladjs, lpats, 0., ts.dispersion[1])
    plt = draw_lines(ts.brillouin.edges)
    draw_lines!(plt, sf; lcidx=sfidx)
    sf2, sf2idx = const_energy_line_in_patches(ltris, ladjs, lpats, 2.0, ts.dispersion[1])
    draw_lines!(plt, sf2; lcidx=sf2idx)
    @savepng plt "surface3"
end


text:
const MONTHS_PER_YEAR = Int64(12)
const HOURS_PER_DAY = Int64(24)
const MINUTES_PER_HOUR = Int64(60)
const MINUTES_PER_DAY = MINUTES_PER_HOUR * HOURS_PER_DAY
const SECONDS_PER_MINUTE = Int64(60)
const SECONDS_PER_HOUR = SECONDS_PER_MINUTE * MINUTES_PER_HOUR
const SECONDS_PER_DAY = SECONDS_PER_HOUR * HOURS_PER_DAY
const MILLISECONDS_PER_SECOND = Int64(1_000)
const MICROSECONDS_PER_MILLISECOND = Int64(1_000)
const NANOSECONDS_PER_MICROSECOND  = Int64(1_000)
const MICROSECONDS_PER_SECOND = Int64(1_000_000)
const NANOSECONDS_PER_MILLISECOND  = Int64(1_000_000)
const NANOSECONDS_PER_SECOND = Int64(1_000_000_000)


text:
export LangJulia

# Data structure for the language.
struct LangJulia
    overwrite_input::Any # Overwrite input
    inline::Any
    dot_fusing::Any # Allow dot fusion
    axpby_header::Any
    alloc_function
    only_overwrite
end
default_alloc_function(k)="similar(A,T)"

"""
    LangJulia(overwrite_input=true,inline=true,dot_fusing=true,axpby_header=:auto,alloc_function,only_overwrite=false)

Code generation in julia language, with optional overwriting of input, inlining
the function and optional usage of dot fusion. The `axpby_header` specifies if axpby function calls should be included in the beginning of the file. The parameter `alloc_function` is a function of three parameters `alloc_function(k)` where `k` is the memory slot (default is `alloc_function(k)=similar(A,T)`). The `only_overwrite` specifies if `f` should be created if the overwrite funtion `f!` contains the actual code.
"""
LangJulia() = LangJulia(true, true, true, :auto, default_alloc_function,false)
LangJulia(overwrite_input) = LangJulia(overwrite_input, true, true, :auto, default_alloc_function,false)
function LangJulia(overwrite_input, inline)
    return LangJulia(overwrite_input, inline, true, :auto, default_alloc_function,false)
end
function LangJulia(overwrite_input, inline, dot_fusing)
    return LangJulia(overwrite_input, inline, dot_fusing, :auto, default_alloc_function,false)
end

# Language specific operations.
comment(::LangJulia, s) = "# $s"

slotname(::LangJulia, i) = "memslots$i"

function assign_coeff(lang::LangJulia, v, i)
    if v == 1
        return (
            "value_one",
            comment(lang, "Saving scalar multiplications using ValueOne()."),
        )
    else
        assign_coeff_basic(lang, v, i)
    end
end

function assign_coeff_basic(lang::LangJulia, v, i)
    T = typeof(v)
    if big(T) == T # High precision coeffs should be parsed as such
        if real(T) == T
            return ("coeff$i", "coeff$i=big\"$v\"")
        else
            vr = real(v)
            vi = imag(v)
            return ("coeff$i", "coeff$i=big\"$vr\" + big\"$vi\"*im")
        end
    else
        return ("coeff$i", "coeff$i=$v")
    end
end

function preprocess_codegen(graph, lang::LangJulia)
    if (lang.dot_fusing)
        return MultiLincombCompgraph(graph) # Merge many lincombs for dot fusion
    else
        return graph
    end
end

# Code generation.
function push_code_matfun_axpby_I!(code)
    # Add code for matfun_axpby! using UniformScaling.
    return push_code_verbatim_string!(
        code,
        "
# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end\n",
    )
end

function push_code_matfun_axpby!(code)
    # Add code for generic matfun_axpby! function.
    return push_code_verbatim_string!(
        code,
        "
# Compute X <- a X + b Y.
function matfun_axpby!(X,a,b,Y)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:m
        @inbounds for j=1:n
            if (b isa ValueOne)
                X[i,j]+=Y[i,j]
            else
                X[i,j]+=b*Y[i,j]
            end
        end
    end
end\n",
    )
end

function function_definition(
    lang::LangJulia,
    graph,
    T,
    funname,
    precomputed_nodes,
)
    code = init_code(lang)
    axpby_header = lang.axpby_header
    if (lang.axpby_header == :auto)
        axpby_header = !lang.dot_fusing
    end

    push_code!(code, "using LinearAlgebra", ind_lvl = 0)
    # If graph has linear combinations, add corresponding axpby functions.
    if any(values(graph.operations) .== :lincomb)
        push_code!(code, "\nstruct ValueOne; end\nValueOne()", ind_lvl = 0)
    end
    lincomb_nodes =
        filter(x -> graph.operations[x] == :lincomb, keys(graph.operations))
    lincomb_with_I =
        filter(y -> any(map(x -> x == :I, graph.parents[y])), lincomb_nodes)
    if (!isempty(lincomb_with_I) && axpby_header)
        push_code_matfun_axpby_I!(code)
    end
    if (!isempty(setdiff!(lincomb_nodes, lincomb_with_I)) && axpby_header)
        push_code_matfun_axpby!(code)
    end

    input_variables = join(precomputed_nodes, ", ")
    inline_string = lang.inline ? "@inline " : ""
    # Generate version a bang version of the function.
    if (lang.overwrite_input)
        if (!lang.only_overwrite)
            push_code!(
                code,
                inline_string * "function $funname($input_variables)",
                ind_lvl = 0,
            )
            #
            eltype_precomputed = join(map(x->"eltype($x)", precomputed_nodes),",");
            push_code!(code,"T=promote_type($eltype_precomputed,$T)");

            for n in precomputed_nodes
                push_code!(code,"$(n)_copy=similar($n,T); copyto!($(n)_copy, $n);");
            end
            copy_input = join(map(x->"$(x)_copy", precomputed_nodes),",");
            push_code!(code, "return $(funname)!($copy_input)")
            push_code!(code, "end", ind_lvl = 0)
        end

        push_code!(code, "", ind_lvl = 0)
        push_code!(
            code,
            inline_string * "function $(funname)!($input_variables)",
            ind_lvl = 0,
        )
    else
        push_code!(
            code,
            inline_string * "function $funname($input_variables)",
            ind_lvl = 0,
        )
    end
    push_code!(
        code,
        "T=promote_type(eltype($(input_variables[1])),$T) " * comment(
            lang,
            "Make it work for many 'bigger' types (matrices and scalars)",
        ),
    )
    return code
end

function function_init(lang::LangJulia, T, mem, graph, precomputed_nodes)
    code = init_code(lang)
    max_nodes = size(mem.slots, 1)

    # Allocate memory for memory slots.
    push_comment!(code, "max_memslots=$max_nodes")
    push_code!(code, "n=size(A,1)")
    # Allocate input.
    push_comment!(
        code,
        "The first slots are precomputed nodes $precomputed_nodes",
    )
    jj = size(precomputed_nodes,1);
    if (!lang.overwrite_input)
        jj = 0;
    end
    for i=jj+1:max_nodes
        thisslotname = slotname(lang, i);
        push_code!(code,"$thisslotname = "*lang.alloc_function(i));
    end

    push_comment!(code,"Assign precomputed nodes memslots ");

    precomp_nodes_string = join(repeat("A",jj),",");
    for (i, n) in enumerate(precomputed_nodes)
        Ak_slot_name = get_slot_name(mem, n)
        if (lang.overwrite_input)
            # Just assign a pointer to the slot to allow overwrite
            push_code!(
                code,
                "$Ak_slot_name=$n " * comment(lang, "overwrite $n"),
            )
        else
            # Otherwise make a copy
            push_code!(code, "copy!($Ak_slot_name,$n)")
        end
    end

    # If needed, allocate identity matrix.
    if has_identity_lincomb(graph)
        push_comment!(code, "Graph has linear combination of identities.")
        push_comment!(code, "The matrix I is explicitly allocated.")
        alloc_slot!(mem, length(precomputed_nodes) + 1, :I)
        I_slot_name = get_slot_name(mem, :I)
        push_code!(code, "$I_slot_name=Matrix{Float64}(I,n,n)")
    else
        push_comment!(code, "Uniform scaling is exploited.")
        push_comment!(code, "No matrix I explicitly allocated.")
    end

    # If needed, initialize variable of type ValueOne for axpby with a=1 or b=1.
    if (any(map(x -> any(x .== 1), values(graph.coeffs))))
        push_code!(code, "value_one=ValueOne()")
    end

    return code
end

function init_mem(lang::LangJulia, max_nof_nodes, precomputed_nodes)
    mem = CodeMem(max_nof_nodes, i -> slotname(lang, i))

    # Make sure the precomputed nodes have memslots
    for (i, n) in enumerate(precomputed_nodes)
        alloc_slot!(mem, i, n)
    end

    return mem
end

function function_end(lang::LangJulia, graph, mem)
    code = init_code(lang)
    retval_node = graph.outputs[end]
    retval = get_slot_name(mem, retval_node)

    push_code!(
        code,
        "return $retval " * comment(lang, "Returning $retval_node"),
    )
    push_code!(code, "end", ind_lvl = 0)
    return code
end

# Add a scaled identity to a matrix in julia code.
function execute_julia_I_op(
    code,
    nodemem,
    non_I_parent_mem,
    non_I_parent_coeff,
    I_parent_coeff,
)
    push_comment!(code, "Add lincomb with identity.")

    if (nodemem != non_I_parent_mem)
        push_code!(code, "copy!($(nodemem),$(non_I_parent_mem))")
    else
        push_comment!(code, "No copy necessary. Inline identity multiple add")
    end

    return push_code!(
        code,
        "matfun_axpby!($(nodemem),$non_I_parent_coeff,$I_parent_coeff,I)",
    )
end

function execute_operation!(
    lang::LangJulia,
    T,
    graph::MultiLincombCompgraph,
    node,
    dealloc_list,
    mem,
)
    op = graph.operations[node]
    if (op != :lincomb)
        return execute_operation_basic!(lang, T, graph, node, dealloc_list, mem)
    else        # Multiple additions goes here

        # Don't overwrite the list
        dealloc_list = deepcopy(dealloc_list)
        setdiff!(dealloc_list, keys(mem.special_names))
        setdiff!(dealloc_list, [:I])

        code = init_code(lang)

        fused_sum = (join("x*" .* string.(graph.parents[node]), '+'))
        push_comment!(code, "Computing $node = $fused_sum")

        # Write the coeffs into appropriate vectors
        coeff_names = Vector()
        coeff_list = graph.coeffs[node]
        parent_mems = Vector()
        id_coeffs = Vector() # Vector of identity additions
        for (i, v) in enumerate(coeff_list)
            (coeff_i, coeff_i_code) = assign_coeff_basic(lang, v, i)
            push_code!(code, coeff_i_code)
            n = graph.parents[node][i]
            if (n != :I)
                push!(coeff_names, coeff_i)
                push!(parent_mems, get_slot_name(mem, n))
            else  # The identites are added in a different vector
                push!(id_coeffs, coeff_i)
            end
        end

        # Determine a target mem: nodemem
        nodemem = nothing
        recycle_parent = nothing
        for n in graph.parents[node]
            if (n in dealloc_list) # If it's about to deallocated
                nodemem = get_slot_name(mem, n)
                recycle_parent = n
                break
            end
        end
        if (isnothing(recycle_parent))
            (nodemem_i, nodemem) = get_free_slot(mem)
            alloc_slot!(mem, nodemem_i, node)
        else
            push_comment!(code, "Smart lincomb recycle $recycle_parent")
        end

        # Print the code
        rhs = join((coeff_names .* ".*") .* parent_mems, " .+ ")
        lhs = nodemem
        if !isempty(rhs)
            push_code!(code, lhs * " .= " * rhs)
        else
            push_code!(code, lhs * " .= 0")
        end

        # Adjust the result with inplace additions of identity
        for c in id_coeffs
            if VERSION > v"1.7.0-DEV.1240"
                push_code!(code, "mul!($nodemem,true,I*$c,true,true)")
            else
                push_comment!(
                    code,
                    "This julia version does not support inplace identity add.",
                )
                push_comment!(
                    code,
                    "Update to Julia 1.7.0 or newer for better performance.",
                )
                push_code!(code, "$nodemem+=I*$c")
            end
        end

        if !(isnothing(recycle_parent))
            # Avoid deallocate recycled parent
            setdiff!(dealloc_list, [recycle_parent])

            # Set the memory pointer
            j = get_slot_number(mem, recycle_parent)
            set_slot_number!(mem, j, node)
        end

        # Deallocate
        for n in dealloc_list
            if n != :I # No memory is ever allocated for the identity matrix.
                i = get_slot_number(mem, n)
                push_comment!(code, "Deallocating $n in slot $i")
                free!(mem, i)
            end
        end

        return (code, nodemem)
    end
end

function execute_operation!(lang::LangJulia, T, graph, node, dealloc_list, mem)
    return execute_operation_basic!(lang, T, graph, node, dealloc_list, mem)
end

# The general base case. Separated for dispatch.
function execute_operation_basic!(
    lang::LangJulia,
    T,
    graph,
    node,
    dealloc_list,
    mem,
)
    op = graph.operations[node]
    parent1 = graph.parents[node][1]
    parent2 = graph.parents[node][2]

    # Don't overwrite the list
    dealloc_list = deepcopy(dealloc_list)
    setdiff!(dealloc_list, keys(mem.special_names))
    setdiff!(dealloc_list, [:I])
    if parent1 != :I || has_identity_lincomb(graph)
        parent1mem = get_slot_name(mem, parent1)
    end
    if parent2 != :I || has_identity_lincomb(graph)
        parent2mem = get_slot_name(mem, parent2)
    end

    code = init_code(lang)
    push_comment!(code, "Computing $node with operation: $op")
    if op == :mult
        # Multiplication has no inline
        (nodemem_i, nodemem) = get_free_slot(mem)
        alloc_slot!(mem, nodemem_i, node)

        push_code!(code, "mul!($nodemem,$parent1mem,$parent2mem)")

    elseif op == :ldiv
        # Left division
        (nodemem_i, nodemem) = get_free_slot(mem)
        alloc_slot!(mem, nodemem_i, node)
        if parent2 == :I
            push_code!(code, "$nodemem=inv($parent1mem)")
        else
            push_code!(code, "$nodemem .=$parent1mem\\$parent2mem")
        end

    elseif op == :lincomb
        coeff_vars = Vector{String}(undef, 2)
        (coeff1, coeff1_code) = assign_coeff(lang, graph.coeffs[node][1], 1)
        push_code!(code, coeff1_code)
        (coeff2, coeff2_code) = assign_coeff(lang, graph.coeffs[node][2], 2)
        push_code!(code, coeff2_code)

        # Use economical memory slots
        if ((parent1 in dealloc_list) || (parent2 in dealloc_list))
            # Smart / inplace: parentX can be used to store newly computed value

            # No allocation needed

            recycle_parent = (parent1 in dealloc_list) ? parent1 : parent2

            push_comment!(code, "Smart lincomb recycle $recycle_parent")

            nodemem = get_slot_name(mem, recycle_parent)

            if parent1 == :I || parent2 == :I
                # BLAS does not work with unform scaling use inplace instead
                if (parent1 == :I)
                    non_I_parent_mem = parent2mem
                    non_I_parent_coeff = coeff2
                    I_parent_coeff = coeff1
                elseif (parent2 == :I)
                    non_I_parent_mem = parent1mem
                    non_I_parent_coeff = coeff1
                    I_parent_coeff = coeff2
                end
                execute_julia_I_op(
                    code,
                    nodemem,
                    non_I_parent_mem,
                    non_I_parent_coeff,
                    I_parent_coeff,
                )
            else
                if (recycle_parent == parent1)
                    push_code!(
                        code,
                        "matfun_axpby!($nodemem,$coeff1,$coeff2,$parent2mem)",
                    )
                else
                    push_code!(
                        code,
                        "matfun_axpby!($nodemem,$coeff2,$coeff1,$parent1mem)",
                    )
                end
            end

            # Avoid deallocated
            setdiff!(dealloc_list, [recycle_parent])

            # Set the memory pointer
            j = get_slot_number(mem, recycle_parent)
            set_slot_number!(mem, j, node)
        else
            # Default behavior:  Allocate new slot
            (nodemem_i, nodemem) = get_free_slot(mem)

            alloc_slot!(mem, nodemem_i, node)

            if (parent1 == :I)
                non_I_parent_mem = parent2mem
                non_I_parent_coeff = coeff2
                I_parent_coeff = coeff1
                execute_julia_I_op(
                    code,
                    nodemem,
                    non_I_parent_mem,
                    non_I_parent_coeff,
                    I_parent_coeff,
                )
            elseif (parent2 == :I)
                non_I_parent_mem = parent1mem
                non_I_parent_coeff = coeff1
                I_parent_coeff = coeff2
                execute_julia_I_op(
                    code,
                    nodemem,
                    non_I_parent_mem,
                    non_I_parent_coeff,
                    I_parent_coeff,
                )
            else
                push_code!(code, "copy!($(nodemem),$parent1mem)")
                push_code!(
                    code,
                    "matfun_axpby!($(nodemem),$coeff1,$coeff2,$parent2mem)",
                )
            end
        end

    else
        error("Unknown operation")
    end

    # Deallocate
    for n in dealloc_list
        if n != :I # No memory is ever allocated for the identity matrix.
            i = get_slot_number(mem, n)
            push_comment!(code, "Deallocating $n in slot $i")
            free!(mem, i)
        end
    end

    return (code, nodemem)
end


text:
# general formula for elementwise activation: output = act(input)
# some examples of act: logistic, tanh, relu, etc
immutable Activation{F,T} <: Transformation
    n::Int
    # α::T  # scaling parameter for some activations
    input::SumNode{T,1}
    output::OutputNode{T,1}

    # construct a new Activation, and link the nodes if it's the identity
    function Activation(n::Int) #, α::T = zero(T))
        input = InputNode(T,n)
        output = OutputNode(T,n)
        # if F == :identity
        #     link_nodes!(output, input)
        # end
        new(n, input, output)
    end
end
Activation(s::Symbol, n::Int) = Activation{s,Float64}(n)

Base.show{F,T}(io::IO, act::Activation{F,T}) = print(io, "$F{$(act.n)}")

input_length(act::Activation) = act.n
output_length(act::Activation) = act.n

# params_length(act::Activation) = 0
# params{F,T}(act::Activation{F,T}) = zeros(T,0)
# grad{F,T}(act::Activation{F,T}) = zeros(T,0)

# ----------------------------------------------------------------------------

const _exp_threshold = 20

# # identity: nothing to do, since we linked the input to output
# transform!(act::Activation{:identity}) = act.output.val
# grad!(act::Activation{:identity}) = act.input.∇
identity′{T<:Number}(x::T, y::T) = one(T)

# for the following, compute the derivative f′(x), where y = act(x) is assumed precomputed
# ref: https://en.wikipedia.org/wiki/Activation_function

logit(x::Real) = log(x / (one(x) - x))
logistic(x::Real) = one(x) / (one(x) + exp(-x))
logistic′{T<:Number}(x::T, y::T) = y * (one(T) - y)

# tanh: act(x) = (eˣ .- e⁻ˣ) ./ (eˣ .+ e⁻ˣ)
tanh′{T<:Number}(x::T, y::T) = one(T) - y^2

softsign{T<:Number}(x::T) = x / (one(T) + abs(x))
softsign′{T<:Number}(x::T, y::T) = (one(T) / (one(T) + abs(x)))^2

relu{T<:Number}(x::T) = max(zero(T), x)
relu′{T<:Number}(x::T, y::T) = x >= zero(T) ? one(T) : zero(T)

# softplus{T<:Number}(x::T) = log(one(T) + exp(x))
# softplus′{T<:Number}(x::T, y::T) = logistic(x)

# set to f(x)=x, f'(x)=1 when x is too big due to floating point math
softplus{T<:Number}(x::T) = x > _exp_threshold ? x : log(one(T) + exp(x))
softplus′{T<:Number}(x::T, y::T) = x > _exp_threshold ? one(T) : logistic(x)

sinusoid{T<:Number}(x::T) = sin(x)
sinusoid′{T<:Number}(x::T, y::T) = cos(x)

gaussian{T<:Number}(x::T) = exp(-(x^2))
gaussian′{T<:Number}(x::T, y::T) = -2x*y

threshold{T<:Number}(x::T) = x > zero(T) ? one(T) : zero(T)
threshold′{T<:Number}(x::T, y::T) = zero(T)

sign′{T<:Number}(x::T, y::T) = zero(T)

# ----------------------------------------------------------------------------

# generic implementations... ensure there's a derivative method of the correct name

const activations = [
    :logistic,
    :tanh,
    :softsign,
    :relu,
    :softplus,
    :sinusoid,
    :gaussian,
    :threshold,
    :sign,
]

for act in activations
    s = string(act)
    f′ = Symbol(s*"′")

    @eval begin
        # elementwise map from input to output
        transform!(act::Activation{Symbol($s)}) = map!($act, act.output.val, act.input.val)

        # backprop gradient calc using specialized derivative
        function grad!(act::Activation{Symbol($s)})
            ∇x = input_grad(act)
            ∇y = output_grad(act)
            x = input_value(act)
            y = output_value(act)
            for i=1:act.n
                ∇x[i] = $f′(x[i], y[i]) * ∇y[i]
            end
        end

        # x-only version
        function $f′(x::Number)
            y = $act(x)
            $f′(convert(typeof(y), x), y)
        end

        value_func(act::Activation{Symbol($s)}) = $act
        deriv_func(act::Activation{Symbol($s)}) = $f′

        # export both functions
        export $act, $f′
    end
end

# ----------------------------------------------------------------------------
# softmax

function transform!{T}(act::Activation{:softmax,T})
    val = act.input.val
    out = act.output.val
    maxinput = maximum(val)
    for i=1:act.n
        # subtract the maxinput to prevent overflow... doesn't change final result
        out[i] = exp(val[i] - maxinput)
    end
    s = one(T) / sum(out)
    for i=1:act.n
        out[i] *= s
    end
    out
end

# the calc is done in the CrossEntropyLoss... just pass that gradient back
function grad!{T}(act::Activation{:softmax,T})
    copy!(act.input.∇, act.output.∇)
end

# ----------------------------------------------------------------------------

default_range(act::Activation) = linspace(-5,5)

# user recipe adds a default x range
@recipe act{F}(act::Activation{F}) = act, default_range(act)

# type recipe converts to a function of xi
@recipe act{A<:Activation}(::Type{A}, act::A) = Transformations.value_func(act)

# ----------------------------------------------------------------------------


text:
export ThresholdGraph, RandomThresholdGraph, CreationSequence, ThresholdRepresentation

"""
`ThresholdGraph(w)` creates a threshold graph using weights from the
vector `w`.

`ThresholdGraph(dw)` creates a threshold graph from a dictionary
mapping vertices to weights.
"""
function ThresholdGraph(w::Array{T,1}) where {T<:Real}
    n = length(w)
    f = Dict{Int,T}()
    for k = 1:n
        f[k] = w[k]
    end
    return ThresholdGraph(f)
end

function ThresholdGraph(dw::Dict{S,T}) where {S,T<:Real}
    G = SimpleGraph{S}()
    vtcs = collect(keys(dw))
    n = length(vtcs)

    for v in vtcs
        add!(G, v)
    end

    for i = 1:n-1
        u = vtcs[i]
        wt_u = dw[u]
        for j = i+1:n
            v = vtcs[j]
            wt_v = dw[v]
            if wt_u + wt_v >= 1
                add!(G, u, v)
            end
        end
    end
    cache_save(G, :ThresholdRepresentation, dw)
    cache_save(G, :name, "Threshold graph")
    return G
end

"""
`RandomThresholdGraph(n)` creates a random threshold graph with `n`
vertices.
"""
function RandomThresholdGraph(n::Int)
    w = rand(n)
    return ThresholdGraph(w)
end

# The following code contributed by Tara Abrishami

"""
`CreationSequence(G)` creates a creation sequence for a threshold
graph `G`. This returns a pair `(seq, vtx_list)` where `seq` is the
creation sequence and `vtx_list` specifies the order in which the
vertices are added when creating `G`. If `G` is not a threshold
graph, then an error is raised.
"""
function CreationSequence(G1::SimpleGraph)
    A = Int[]
    T = eltype(G1)
    V = T[]
    G = deepcopy(G1)
    while length(vlist(G)) != 0
        r = false
        for v in vlist(G)
            if deg(G, v) == 0
                delete!(G, v)
                r = true
                unshift!(A, 0)
                unshift!(V, v)
                break
            end
            if deg(G, v) == NV(G) - 1
                delete!(G, v)
                r = true
                unshift!(A, 1)
                unshift!(V, v)
                break
            end
        end
        if !r
            error("This graph is not a threshold graph")
        end
    end
    return A, V
end

"""
`ThresholdRepresentation(G)` returns a threshold representation of
`G`. This returns a dictionary mapping vertices of `G` to `Rational`
weights. An error is raised if `G` is not a threshold graph.
"""
function ThresholdRepresentation(G::SimpleGraph)
    if cache_check(G, :ThresholdRepresentation)
        return cache_recall(G, :ThresholdRepresentation)
    end
    A, V = CreationSequence(G)

    prev::Int = 0
    prevVal::Rational = 1 // 3
    T = eltype(G)
    D = Dict{T,Rational}()
    D[V[1]] = prevVal
    small::Rational = prevVal
    large::Rational = prevVal
    for i = 2:length(A)
        if A[i] == prev
            D[V[i]] = prevVal
        elseif A[i] == 0
            D[V[i]] = (1 - large) / 2
            prevVal = (1 - large) / 2
            if (1 - large) / 2 < small
                small = (1 - large) / 2
            end
            if (1 - large) / 2 > large
                large = (1 - large) / 2
            end
        elseif A[i] == 1
            D[V[i]] = 1 - small
            prevVal = 1 - small
            if 1 - small > large
                large = 1 - small
            end
            if 1 - small < large
                small = 1 - large
            end
        end
        prev = i
    end
    cache_save(G, :ThresholdRepresentation, D)
    return D
end


text:
export DS, run_lt_test
const DS = DirectSearch

#= Following functions sourced from https://www.sfu.ca/~ssurjano/optimization.html =#
#Valley
#soln: -1.0316 @ [±0.0898, ∓0.7126]
camel6(x)=(4-2.1x[1]^2.0 + x[1]^4/3)x[1]^2.0 + x[1]x[2] + (-4 + 4x[2]^2.0)x[2]^2.0
#soln: 0 @ [1,...,1]
rosenbrock(x;a=1,b=100,d=length(x)) = sum([b*(x[i+1] - x[i]^2)^2 + (x[i]-a)^2 for i=1:d-1])

#Local Minima
ackley(x;a=20,b=0.2,c=2pi,d=length(x)) = -a * exp(-b*√(sum(x.^2)/d)) - exp(sum(cos.(c.*x))/d) + a + exp(1)
eggholder(x) = -(x[2]+47)sin(√abs(x[2] + x[1]/2 + 47)) - x[1]sin(√abs(x[1]-x[2]+47))

#Example Benchmarks from papers
#Audet & Dennis 2006 benchmark one, take x₀=[-2.1,1.7]
#soln: 0 @ [0,0]
bm_1(x;c=[30;40],d=[-30;-40]) = (1-exp(-norm(x)^2)) * max(norm(x-c)^2, norm(x-d)^2)

#Audet & Dennis 2006 benchmark three, should take x₀=[0,...,0]
#soln: -√(3)n @ [-√3,...,-√3]
bm_3(x)=sum(x)
bm_3_con(x)=sum(x.^2)-3*length(x)

function run_lt_test(;f=camel6,n=2,
                     constraints::Vector=[],
                     initial=zeros(Float64,n),
                     lim=1000,
                     lb=-5*ones(Float64,n),
                     ub=5*ones(Float64,n),
                     run=true)
    p = DSProblem(n)
    SetObjective(p, f)
    SetInitialPoint(p, initial)
    SetMaxEvals(p, 1)
    for c in constraints
        AddExtremeConstraint(p, c)
    end
    p.iteration_limit = lim
    SetVariableRanges(p, lb, ub)
    if run
        DS.Optimize!(p)
        println("$(p.iteration)\t$(p.x)\t$(p.x_cost)")
    end
    return p
end


text:
using Images

canny_edges = canny(img, sigma = 1.4, upperThreshold = 0.80, lowerThreshold = 0.20)


text:
###############################################################################
#
#    Copyright (C) 2015  VoxForge
#
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################

function mkclscript(monophones0, tree_hed, folder)
  hmmlist=open(tree_hed,"a");

  monophones0_arr=open(readlines, monophones0)
  for i=2:4
    for phoneln=monophones0_arr
      phone=chomp(phoneln)
      write(hmmlist,"TB 350 \"ST_$phone\_$i\_\" {(\"$phone\",\"*-$phone\+*\",\"$phone\+*\",\"*-$phone\").state[$i]}\n")
    end
  end

  write(hmmlist,"\n")
  write(hmmlist,"TR 1\n")
  write(hmmlist,"\n")
  write(hmmlist,"AU \"$(folder)/fulllist\" \n")
  write(hmmlist,"CO \"$(folder)/tiedlist\" \n")
  write(hmmlist,"\n")
  write(hmmlist,"ST \"$(folder)/trees\" \n")

  close(hmmlist)
end

# if called from command line
if length(ARGS) > 0
  if ! isfile(ARGS[1])
    error("can't find monophones0 file: $(ARGS[1])")
  end
  if ! isfile(ARGS[2])
    error("can't find tree.hed file: $(ARGS[2])")
  end
  if length(ARGS) == 2
    mkclscript(ARGS[1], ARGS[2], "." )
  elseif length(ARGS) == 3
    if ! isdir(ARGS[3])
      error("can't find directory: $(ARGS[3])")
    end

    mkclscript(ARGS[1], ARGS[2], ARGS[3] )
  end

  if length(ARGS) > 3
    error("mkclscript: too many arguments for call from command line")
  end

end

max_stars_repo_name: Hyperparticle/lct-master

max_issues_repo_name: Hyperparticle/lct-master

max_forks_repo_name: Hyperparticle/lct-master


text:
module loadAlessioData
#=
	This module implements convenience functions for
	- loading of Alessio's data
	- reading and writing Fortran observation data files
	- plotting and analysis of the data files
	- interpolating Alessio's data onto appropriate grids
=#

# For:	     data input/output             Printing        Plotting        Data Analysis
#      /--------------------------\  /------------------\  /----\  /----------------------------\
using DelimitedFiles, FortranFiles, Printf, ProgressBars, PyPlot, Peaks, Dierckx, Interpolations

export getPossibleIndices, loadFile, spatialMask, temporalMask, mutualMask, sampleData, plotTrace!, plotState!, recordingMovie, readObsFile, stateInterpolants, generateObservationFileSequence!, spatialMaskCentroid, plotSamplingSquare!

#
# Convenience functions to identify the available data files and load the data consistently
#
function getPossibleIndices(basepath)
	# lists possible data file indices
	files = readdir(basepath)
	indices = []
	for file in files
		if cmp(file[1:6], "Erased") == 0
			push!(indices, parse(Int,file[29:31]))
		end
	end
	sort!(indices)
	unique!(indices)

	return indices
end

function loadData(filename) # loadData function code due to Max Comstock
	# Pre-allocate data array
	total_elements = Int(stat(filename).size / sizeof(UInt16))
	obs_data = Array{UInt16, 1}(undef, total_elements)

	# Read the data and convert from big-endian to system endianness
	read!(filename, obs_data)
	obs_data = ntoh.(obs_data)

	# Shape into frames of 128×128 readings
	obs_data = reshape(obs_data, :, 128, 128)

	# rescale data to [0,1] and convert to Float64 (overkill; Float32 would work, Float16 is too small)
	obs_data = obs_data ./ (maximum(obs_data[:])-minimum(obs_data[:]))

	return obs_data
end

function transformData!(obs)
	# transforms data in-place for consistent orientation; array is now [EPI,ENDO][t,x,y]
	obs[1] .= obs[1][:,:,end:-1:1]
	obs[2] .= obs[2][:,end:-1:1,:]
	return nothing
end

function loadFile(basepath, index::Int)
	# loads data from a pair of files
	basefile = "Erased_2011-05-18_Exp000_Rec"*@sprintf("%03i",index)
	# EPI first, then ENDO
	filenames= [basepath*basefile*"_Cam0-Yellow", basepath*basefile*"_Cam1-Blue"]
	dats = []
	for fname in filenames
		dat = loadData(fname)
		push!(dats, dat)
	end

	# re-orient data arrays for consistency
	transformData!(dats)

	return dats
end

#
# Basic analysis functions to generate masks for the experimental data arrays
#

function spatialMask(obs_data)
	# generates BitArray of size [Nx,Ny] which corresponds to valid spatial indices
	N = size(obs_data)
	mask = zeros(Bool, N[2], N[3])
	for n in 1:N[2], m in 1:N[3]
		for t in 1:N[1]
			if mask[n,m] == false && obs_data[t,n,m] > 0
				mask[n,m] = 1
			end
		end
	end

	return mask
end

function temporalMask(obs_data)
	# generates BitArray of size [Nt] which corresponds to valid temporal indices
	N = size(obs_data)
	mask = zeros(Bool, N[1])
	for t in 1:N[1]
		for n in 1:N[2], m in 1:N[3]
			if mask[t] == false && obs_data[t,n,m] > 0
				mask[t] = 1
			end
		end
	end

	return mask
end

function mutualMask(mask1, mask2)
	# pairwise function to reduce the data ranges for the EPI and ENDO recordings to a consistent subset
	if length(size(mask1)) == 1 && length(size(mask2)) == 1 && length(mask1) != length(mask2)
		println("Temporal masks are of different lengths; truncating to shorter")
		tmp = min(length(mask1), length(mask2))
		mask1 = mask1[1:tmp]
		mask2 = mask2[1:tmp]
	end
	@assert size(mask1) == size(mask2)

	return mask1 .* mask2
end

function mutualMask(masks)
	# array input of masks
	n = length(masks)
	mutmask = masks[1]
	for m = 2:n
		mutmask = mutualMask(mutmask, masks[m])
	end

	return mutmask
end

function spatialMaskCentroid(mask; dx=0.06)

	N = size(mask)
	ctr = 0.0;
	for x in 1:N[1], y in 1:N[2]
		ctr += (x + y*im)*mask[x,y];
	end
	return Int.(round.([imag(ctr), real(ctr)]./sum(mask))).*dx # note the swap in real/imag
end

#
# Sampling functions
#

function sampleData(obs, sampleIndex::Array{T,1}, mutualTimeMask, mutualSpaceMask; dt=2.0) where T <: Integer
	# generates (t,V(t)) arrays for EPI and ENDO recordings in obs
	# first check that the spatial indices are valid
	@assert mutualSpaceMask[sampleIndex[1],sampleIndex[2]] == true
	# generate times from the mutualTimeMask
	t = (1:length(mutualTimeMask))[mutualTimeMask]
	# select obs[t, sampleIndex[1], sampleIndex[2]] for EPI and ENDO
	V = [obs[n][t, sampleIndex[1], sampleIndex[2]] for n=1:2]

	return (t.*dt, V)
end

function sampleData(obs, samplePoint::Array{T,1}, mutualTimeMask, mutualSpaceMask; h=0.02) where T <: AbstractFloat
	sampleIndex = Int.(round.(realSamplePoint ./ h))
	return sampleData(obs, sampleIndex, mutualTimeMask, mutualSpaceMask)
end

#
# Analysis functions for the data
#

function analyzeVoltage(t, V; threshold=0.25)
	#=
	This function uses Dierckx.jl for interpolations and root-finding.
	I would like to use Interpolations.jl, but there's no in-built root-
	finding, which would mean using Roots.jl, so might as well just use Dierkx.jl.
	=#
	# Interpolate the voltage (at a specific spatial point) over time
	Vt = Spline1D(t[:], V.-threshold, k=3);
	R = roots(Vt; maxn=Int(5e3));		# time points R: V(R) == threshold
	D = derivative(Vt, R);		# V'(R)

	# storage for APD, DI, APA for this V(t)
	APD = Float64[]
	DI  = Float64[]
	APA = Float64[]

	# for each root, check whether V'(R)≷ 0 to detect beginning or end of AP
	for n in 1:length(R)-1
	       if D[n] > 0 && D[n+1] < 0
		       push!(APD, R[n+1]-R[n])
		       push!(APA, threshold+maximum(Vt(R[n]:R[n+1])))
	       elseif D[n] < 0 && D[n+1] > 0
		       push!(DI, R[n+1]-R[n])
	       end
	end
	return (APD, DI, APA)
end

function peaks(V; minimumProminence=0.9)

	pks, vals = findmaxima(V)
	pks, proms = peakproms(pks, V)

	inds = (1:length(pks))[proms .>= minimumProminence]

	return (inds, pks, proms)
end

function stateInterpolants(states, mutualTimeMask, mutualSpaceMask; dt=2.0, dx=0.06, dy=0.06)

	tlims = (1:length(mutualTimeMask))[mutualTimeMask][[1;end]]
	t0 = dt.*(tlims[1]:tlims[2])
	x0 = dx.*((1:size(mutualSpaceMask,1)).-1)
	y0 = dy.*((1:size(mutualSpaceMask,2)).-1)

	stateSplines = [ CubicSplineInterpolation((t0, x0, y0), state[(1:length(mutualTimeMask))[mutualTimeMask],:,:]) for state in states ]

	return stateSplines
end

function square(;theta=0.0, l=3.0, center=[3.5;3.5])
	# generates the bounding corners (with first and last repeating) of the square for sampling the states
	if theta ≠ 0.0
		@warn "θ ≠ 0 is currently unimplemented, sorry!"
	end
	if abs(theta) >= pi/2
		@warn "θ is mapped to (-π/2,+π/2)"
		theta = mod(theta, pi/2)
	end
	corners = center .+ l .* ([cos(theta) sin(theta); -sin(theta) cos(theta)] * [-1 1 1 -1; -1 -1 1 1]./2);
	corners = hcat(corners[:,:], corners[:,1])

	return corners
end

function coordinateVectors(corners)
	# corners should be generated from the square function
	# corners = square(theta=theta, l=l, center=center)
	x⃗ = corners[:,2].-corners[:,1]
	y⃗ = corners[:,4].-corners[:,1]
	# s.t.  ∀(a,b) ∈ [-0.5, +0.5]: center .+ a.*x⃗ .+ b.*y⃗ ∈ square
	return (x⃗, y⃗)
end

#
# Data plotting functions
#

function plotTrace!(fig, axs, t, V; threshold=0.5) where T <: Integer

	labels = ["EPI", "ENDO"]
	markers = ["+","x"]

	for ax in axs
		ax.cla()
	end

	for n in 1:length(V) # [EPI, ENDO]
		axs[1].plot(t, V[n], "-", label=labels[n], alpha=0.7)
		inds, pks, proms = peaks(V[n]; minimumProminence=0.1)
		axs[1].plot(t[pks[inds]], V[n][pks[inds]], markers[n]*"k")
		axs[1].plot(t, threshold.*ones(size(t)), "--k", label="", alpha=0.7)

		BCL = diff(t[pks[inds]]);
		binwidth=2.0; bins=((round(minimum(BCL)/binwidth)-1)*binwidth):binwidth:((round(maximum(BCL)/binwidth)+1)*binwidth);
		axs[2].hist(BCL, bins=bins, density=true, log=true, label=labels[n], color="C$(n-1)", alpha=0.7)

		APD, DI, APA = analyzeVoltage(t, V[n]; threshold=threshold)
		axs[3].errorbar(APD, APA, xerr=2.0, yerr=0.05, fmt=".", label=labels[n], alpha=0.7)

		cBCL = min(length(APD),length(DI));
		BCL = APD[1:cBCL] .+ DI[1:cBCL]
		axs[2].hist(BCL, bins=bins, density=true, histtype="step", color="C$(n-1)", log=true, label="APD+DI,"*labels[n], alpha=0.7)
	end
	axs[1].set_xlabel("\$ t \$ [ms]")
	axs[1].set_ylabel("\$ V(t) \$ [a.u.]")
	axs[1].legend(loc=0,edgecolor="none")
	axs[2].set_xlabel("BCL [ms]")
	axs[2].set_ylabel("\$P(\$BCL\$)\$")
	axs[2].legend(loc=0,edgecolor="none", ncol=2)
	axs[3].set_xlabel("APD [ms]")
	axs[3].set_ylabel("APA [a.u.]")
	axs[3].legend(loc=0,edgecolor="none")

	return nothing
end

function plotState!(fig, axs, t, obs, spaceMasks, timeMasks, N, spacemask, timemask; drawColorbar=false)

	if drawColorbar  # only draw the colorbar on the first frame
		imm = axs[1].imshow(obs[1][t,:,:], aspect="equal", vmin=0.0, vmax=1.0, cmap="Oranges", interpolation="none", origin="lower")
		clb = fig.colorbar(imm, ax=axs, aspect=50)
	end
	labels = ["EPI", "ENDO"]

	for n=1:2
		axs[n].cla()
		axs[n].imshow(obs[n][t,:,:], alpha=Float64.(spaceMasks[n]), aspect="equal", vmin=0.0, vmax=1.0, cmap="Oranges", interpolation="none", origin="lower")
		axs[n].contour(Float64.(spacemask), "k", levels=[0.5], linewidths=0.5)
		axs[n].axis("off")
		axs[n].set_xticks([])
		axs[n].set_yticks([])
		axs[n].set_title("$(labels[n]), \$ t=\$"*"$(t*2) ms")
	end
	axs[1].contour(Float64.(spaceMasks[2]), "k", levels=[0.5], linestyles="dashed", linewidths=0.5)
	axs[2].contour(Float64.(spaceMasks[1]), "k", levels=[0.5], linestyles="dashed", linewidths=0.5)

	return nothing
end

function plotSamplingSquare!(fig, axs, corners)

	axs[1].plot(corners[1,:]./0.06, corners[2,:]./0.06, "--k", linewidth=1)
	axs[2].plot(corners[1,:]./0.06, corners[2,:]./0.06, "--k", linewidth=1)

	return nothing
end

function recordingMovie(obs, spaceMasks, timeMasks, spacemask, timemask; figname="recording")

	N = size(obs[1])

	rm("./frames/", force=true, recursive=true)
	mkdir("./frames/")

	println("Making movie frames:")
	fig,axs = PyPlot.subplots(1,2,figsize=(8,3),constrained_layout=true)

	for t in ProgressBar((1:length(timemask))[timemask])

		plotState!(fig, axs, t, obs, spaceMasks, timeMasks, N, spacemask, timemask; drawColorbar=(t==(1:length(timemask))[timemask][1]))

		imname = @sprintf("%04i",t)
		PyPlot.savefig("./frames/$(figname)_$(imname).png")
	end

	PyPlot.close("all")

	println("Running frames through FFMPEG:")
	t = (1:length(timemask))[timemask][1]
	ffmpegCommand = `ffmpeg -r 60 -f image2 -start_number $(t) -i ./frames/$(figname)_%04d.png -c:v hevc -crf 25 -pix_fmt yuv420p ./dynamics/$(figname)_dynamics.mp4`;
	println(ffmpegCommand)
	run(ffmpegCommand);

	rm("./frames/", force=true, recursive=true)

	return nothing
end

#
# Observation file functions
#

function readObsFile(obsFile::String; verbose=false)
	#=
	obsFile : string for observation filename
	=#
	f = FortranFile(obsFile)
	#obs = Array{Float32,1}(undef, 6) # not needed, but useful for reference
	Obs = [] # accumulator
	while !eof(f)
		obs = read(f, (Float32,6));
		if verbose; println("Read observation $(obs)"); end
		push!(Obs, obs)
	end
	if eof(f)
		close(f)
	end
	return Obs
end

function fillObservations!(obs, stateSplines, t; theta=0.0, l=3.0, center=[3.5;3.5], min_x = 0, max_x = 200, min_y = 0, max_y = 200, err = 0.025)
	# states is a Array{Array[Float64,3},1}
	#	states = [EPI, ENDO]
	# Obs is a Array{Array{Float32,1},1} where each element of the top array is of form:
	#	[id, x, y, z, val, err] = [2819.0   1.00  89.00    100.00    2.15E-01    5.00E-02]
	# this updates, in-place, the val and err according to the interpolated state at (x,y,z)

	# the data has shape nt x 2 x 128 x 128
	# while the observations are 1 x 2 x 200 x 200
	# so we need to interpolate the data into two continuous (x,y) functions
	# and evaluate the functions at
	# the interpolated relative observation positions {0,1} x [0,1], x [0,1]
	# where the data has (0,0), (0,1), (1,0), and (1,1) at a prescribed and inscribed square

	corners = square(theta=theta, l=l, center=center)

	for n in 1:length(obs)

		# get relative position (within square)
		#x = corners[1,1] + l*cos(theta)*(obs[n][3]-min_x)/(max_x-min_x)
		#y = corners[2,1] + l*sin(theta)*(obs[n][4]-min_y)/(max_y-min_y)
		x = corners[1,1] + (corners[1,2]-corners[1,1])*(obs[n][3]-min_x)/(max_x-min_x)
		y = corners[2,1] + (corners[2,3]-corners[2,1])*(obs[n][4]-min_y)/(max_y-min_y)

		if obs[n][2] == 1.0		# if on bottom => EPI sample
			obs[n][5] = stateSplines[1](t,y,x)
		elseif obs[n][2] == 50.0	# if on top => ENDO sample
			obs[n][5] = stateSplines[2](t,y,x)
		end

		obs[n][6] = err 		# What is the error in the observation here
	end
	return nothing				# obs is modified in-place
end

function writeObsFile(obsFile::String, Obs)
	#=
	obsFile : string for observation filename
	Obs: observation array
	=#
	f = FortranFile(obsFile, "w")
	#obs = Array{Float32,1}(undef, 6) # not needed, but useful for reference
	for n in 1:length(Obs)
		write(f, Obs[n])
		@debug "Wrote observation $(Obs[n])"
	end
	close(f)
	@debug "Wrote observation file $(obsFile)"
	return nothing
end

function generateObservationFileSequence!(obs, stateSplines, obsDir, sequenceName, sampleTimes; theta=0.0, l=3.0, center=[3.5;3.5], min_x = 0, max_x = 200, min_y = 0, max_y = 200, err = 0.025)
	#=
	We loop over the sampleTimes, modifying the template observation array obs, and writing the array to a new obs file
	=#
	@show center
	println("Writing new observation files to $(obsDir):")
	for (n,t) in ProgressBar(enumerate(sampleTimes))
		fillObservations!(obs, stateSplines, t; theta=theta, l=l, center=center, min_x = min_x, max_x = max_x, min_y = min_y, max_y = max_y, err = 0.025)
		obsFile = @sprintf("%s/%s%04d.dat", obsDir, sequenceName, n)
		writeObsFile(obsFile, obs)
	end
	return nothing
end

end

max_stars_repo_name: cmarcotte/puttingTheDataInDA

max_issues_repo_name: cmarcotte/puttingTheDataInDA

max_forks_repo_name: cmarcotte/puttingTheDataInDA