text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_optimized)" --startup-file=no %s 2>&1 | FileCheck %s

# Intentionally bad code
function matmul!(C, A, B)
    for i in 1:size(A, 1)
        for k in 1:size(A, 2)
            for j in 1:size(B, 2)
                C[i, j] += A[i, k] * B[k, j]
            end
        end
    end
   return nothing
end

emit(matmul!, Matrix{Float64}, Matrix{Float64}, Matrix{Float64})

# CHECK:   func nested @"Tuple{typeof(Main.matmul!), Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}}"(%arg0: !jlir<"typeof(Main.matmul!)">, %arg1: !jlir<"Array{Float64, 2}">, %arg2: !jlir<"Array{Float64, 2}">, %arg3: !jlir<"Array{Float64, 2}">) -> !jlir.Nothing attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %1 = "jlir.arraysize"(%arg2, %0) : (!jlir<"Array{Float64, 2}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %2 = "jlir.sle_int"(%0, %1) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %3 = "jlir.constant"() {value = #jlir<"0">} : () -> !jlir.Int64
# CHECK-NEXT:     %4 = "jlir.ifelse"(%2, %1, %3) : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %5 = "jlir.slt_int"(%4, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%5)[^bb3, ^bb2] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb2:  // pred: ^bb1
# CHECK-NEXT:     %6 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %7 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %8 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%6, %7, %8)[^bb4] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb3:  // pred: ^bb1
# CHECK-NEXT:     %9 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%9, %0, %0)[^bb4] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb4(%10: !jlir.Bool, %11: !jlir.Int64, %12: !jlir.Int64):  // 2 preds: ^bb2, ^bb3
# CHECK-NEXT:     %13 = "jlir.not_int"(%10) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%13, %11, %12)[^bb28, ^bb5] {operand_segment_sizes = dense<[1, 0, 2]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb5(%14: !jlir.Int64, %15: !jlir.Int64):  // 2 preds: ^bb4, ^bb27
# CHECK-NEXT:     %16 = "jlir.constant"() {value = #jlir<"2">} : () -> !jlir.Int64
# CHECK-NEXT:     %17 = "jlir.arraysize"(%arg2, %16) : (!jlir<"Array{Float64, 2}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %18 = "jlir.sle_int"(%0, %17) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %19 = "jlir.ifelse"(%18, %17, %3) : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %20 = "jlir.slt_int"(%19, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%20)[^bb7, ^bb6] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb6:  // pred: ^bb5
# CHECK-NEXT:     %21 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %22 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %23 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%21, %22, %23)[^bb8] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb7:  // pred: ^bb5
# CHECK-NEXT:     %24 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%24, %0, %0)[^bb8] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb8(%25: !jlir.Bool, %26: !jlir.Int64, %27: !jlir.Int64):  // 2 preds: ^bb6, ^bb7
# CHECK-NEXT:     %28 = "jlir.not_int"(%25) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%28, %26, %27)[^bb23, ^bb9] {operand_segment_sizes = dense<[1, 0, 2]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb9(%29: !jlir.Int64, %30: !jlir.Int64):  // 2 preds: ^bb8, ^bb22
# CHECK-NEXT:     %31 = "jlir.arraysize"(%arg3, %16) : (!jlir<"Array{Float64, 2}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %32 = "jlir.sle_int"(%0, %31) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %33 = "jlir.ifelse"(%32, %31, %3) : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %34 = "jlir.slt_int"(%33, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%34)[^bb11, ^bb10] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb10:  // pred: ^bb9
# CHECK-NEXT:     %35 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %36 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %37 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%35, %36, %37)[^bb12] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb11:  // pred: ^bb9
# CHECK-NEXT:     %38 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%38, %0, %0)[^bb12] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb12(%39: !jlir.Bool, %40: !jlir.Int64, %41: !jlir.Int64):  // 2 preds: ^bb10, ^bb11
# CHECK-NEXT:     %42 = "jlir.not_int"(%39) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%42, %40, %41)[^bb18, ^bb13] {operand_segment_sizes = dense<[1, 0, 2]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb13(%43: !jlir.Int64, %44: !jlir.Int64):  // 2 preds: ^bb12, ^bb17
# CHECK-NEXT:     %45 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %46 = "jlir.arrayref"(%45, %arg1, %14, %43) : (!jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Int64, !jlir.Int64) -> !jlir.Float64
# CHECK-NEXT:     %47 = "jlir.arrayref"(%45, %arg2, %14, %29) : (!jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Int64, !jlir.Int64) -> !jlir.Float64
# CHECK-NEXT:     %48 = "jlir.arrayref"(%45, %arg3, %29, %43) : (!jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Int64, !jlir.Int64) -> !jlir.Float64
# CHECK-NEXT:     %49 = "jlir.mul_float"(%47, %48) : (!jlir.Float64, !jlir.Float64) -> !jlir.Float64
# CHECK-NEXT:     %50 = "jlir.add_float"(%46, %49) : (!jlir.Float64, !jlir.Float64) -> !jlir.Float64
# CHECK-NEXT:     %51 = "jlir.arrayset"(%45, %arg1, %50, %14, %43) : (!jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Float64, !jlir.Int64, !jlir.Int64) -> !jlir<"Array{Float64, 2}">
# CHECK-NEXT:     %52 = "jlir.==="(%44, %33) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%52)[^bb15, ^bb14] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb14:  // pred: ^bb13
# CHECK-NEXT:     %53 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %54 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%53, %54, %45)[^bb16] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb15:  // pred: ^bb13
# CHECK-NEXT:     %55 = "jlir.add_int"(%44, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %56 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%55, %55, %56)[^bb16] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb16(%57: !jlir.Int64, %58: !jlir.Int64, %59: !jlir.Bool):  // 2 preds: ^bb14, ^bb15
# CHECK-NEXT:     %60 = "jlir.not_int"(%59) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%60)[^bb18, ^bb17] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb17:  // pred: ^bb16
# CHECK-NEXT:     "jlir.goto"(%57, %58)[^bb13] : (!jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb18:  // 2 preds: ^bb12, ^bb16
# CHECK-NEXT:     %61 = "jlir.==="(%30, %19) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%61)[^bb20, ^bb19] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb19:  // pred: ^bb18
# CHECK-NEXT:     %62 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %63 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %64 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%62, %63, %64)[^bb21] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb20:  // pred: ^bb18
# CHECK-NEXT:     %65 = "jlir.add_int"(%30, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %66 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%65, %65, %66)[^bb21] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb21(%67: !jlir.Int64, %68: !jlir.Int64, %69: !jlir.Bool):  // 2 preds: ^bb19, ^bb20
# CHECK-NEXT:     %70 = "jlir.not_int"(%69) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%70)[^bb23, ^bb22] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb22:  // pred: ^bb21
# CHECK-NEXT:     "jlir.goto"(%67, %68)[^bb9] : (!jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb23:  // 2 preds: ^bb8, ^bb21
# CHECK-NEXT:     %71 = "jlir.==="(%15, %4) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%71)[^bb25, ^bb24] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb24:  // pred: ^bb23
# CHECK-NEXT:     %72 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %73 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %74 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%72, %73, %74)[^bb26] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb25:  // pred: ^bb23
# CHECK-NEXT:     %75 = "jlir.add_int"(%15, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %76 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%75, %75, %76)[^bb26] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb26(%77: !jlir.Int64, %78: !jlir.Int64, %79: !jlir.Bool):  // 2 preds: ^bb24, ^bb25
# CHECK-NEXT:     %80 = "jlir.not_int"(%79) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%80)[^bb28, ^bb27] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb27:  // pred: ^bb26
# CHECK-NEXT:     "jlir.goto"(%77, %78)[^bb5] : (!jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb28:  // 2 preds: ^bb4, ^bb26
# CHECK-NEXT:     %81 = "jlir.constant"() {value = #jlir.nothing} : () -> !jlir.Nothing
# CHECK-NEXT:     "jlir.return"(%81) : (!jlir.Nothing) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: femtomc/brutus

max_issues_repo_name: femtomc/brutus

max_forks_repo_name: femtomc/brutus


text:
using CSV
using Plots
pgfplotsx()

function plot_random16(
    results=CSV.File(joinpath(@__DIR__, "..", "results", "random16.csv"));
    kwargs...
)
    dim = getproperty.(results, :dim)
    te = getproperty.(results, :te) ./ 1e9
    ti = getproperty.(results, :ti) ./ 1e9
    tc = getproperty.(results, :tc) ./ 1e9

    plt = plot(
        dim, tc;
        label="cohomology", markershape=:rect,
        legend=:topleft, xlabel="dim", ylabel="t[s]",
        title="",
        kwargs...,
    )
    plot!(plt, dim, te, label="Eirene", markershape=:circle, color=3)
    plot!(plt, dim, ti, label="involuted", markershape=:diamond, color=2)

    plt_log = plot(
        dim, tc;
        label="cohomology", markershape=:rect,
        legend=false, xlabel="dim", ylabel="t[s]",
        title="",
        yscale=:log10,
        kwargs...,
    )
    plot!(plt_log, dim, te, label="Eirene", markershape=:circle, color=3)
    plot!(plt_log, dim, ti, label="involuted", markershape=:diamond, color=2)
    return plot(plt, plt_log, size=(600, 300))
end

function plot_scaling(
    results=CSV.File(joinpath(@__DIR__, "..", "results", "gcycle.csv"));
    left_ylim=nothing,
    right_ylim=nothing,
    kwargs...,
)
    N = getproperty.(results, :N)
    te = getproperty.(results, :te) ./ 1e9
    ti = getproperty.(results, :ti) ./ 1e9
    tc = getproperty.(results, :tc) ./ 1e9
    if !isnothing(left_ylim)
        left_kwargs = (; kwargs..., ylim=left_ylim)
    else
        left_kwargs = kwargs
    end
    if !isnothing(right_ylim)
        right_kwargs = (; kwargs..., ylim=right_ylim)
    else
        right_kwargs = kwargs
    end

    plt = plot(
        N, tc;
        label="cohomology", markershape=:rect,
        legend=:topleft, xlabel="N", ylabel="t[s]",
        title="",
        left_kwargs...,
    )
    plot!(plt, N, te, label="Eirene", markershape=:circle, color=3)
    plot!(plt, N, ti, label="involuted", markershape=:diamond, color=2)

    plt_log = plot(
        N, tc;
        label="cohomology", markershape=:rect,
        legend=false, xlabel="N", ylabel="t[s]",
        title="",
        yscale=:log10,
        right_kwargs...,
    )
    plot!(plt_log, N, te, label="Eirene", markershape=:circle, color=3)
    plot!(plt_log, N, ti, label="involuted", markershape=:diamond, color=2)
    return plot(plt, plt_log, size=(600,300))
end

ext: jl


text:
module DashTable
    function __init__()
        @warn "The `DashTable` package is deprecated. All DashTable functions have been moved to the Dash package. Please switch to `using Dash` instead."
    end
end

max_stars_repo_name: waralex/dash-table

max_issues_repo_name: waralex/dash-table

max_forks_repo_name: waralex/dash-table


text:
module Graphs

using GZip # NOTE: for debugging purposes
using Random
using Rematch
using ProgressMeter

import JSON

# ---------------------------
# functions to extend in submodules

export pair
function pair(item) end

export sequence
function sequence(obj, name; gaps=false)     end
function sequence(obj; gaps=false)           end
function sequence!(s, obj, name; gaps=false) end
function sequence!(s, obj; gaps=false)       end

export reverse_complement, reverse_complement!
function reverse_complement(item)  end
function reverse_complement!(item) end

export marshal, marshal_fasta, marshal_json, marshal_gfa
function marshal_fasta(io::IO, x; opt=nothing) end
function marshal_json(io::IO, x; opt=nothing) end
function marshal_gfa(io::IO, x; opt=nothing) end

function marshal(io::IO, x; fmt=:fasta, opt=nothing)
    @match fmt begin
        :fasta || :fa => return marshal_fasta(io, x; opt)
        :json         => return marshal_json(io, x; opt)
        :gfa          => return marshal_gfa(io, x; opt)
        _ => error("$fmt not a recognized output format")
    end
end

export unmarshal

export serialize
function serialize(io::IO, x) end

# ------------------------------------------------------------------------
# aux types

Maybe{T} = Union{T,Nothing}

# aliases
"""
	SNPMap = Dict{Int,UInt8}

A sparse array of single nucleotide polymorphisms relative to a consensus.
The key is the locus of the mutation; the value is the modified nucleotide.
"""
const SNPMap = Dict{Int,UInt8}
"""
	InsMap = Dict{Tuple{Int,Int},Array{UInt8,1}}

A sparse array of insertion sequences relative to a consensus.
The key is the (locus(after),offset) of the insertion; the value is the sequence.
"""
const InsMap = Dict{Tuple{Int,Int},Array{UInt8,1}}
"""
	DelMap = Dict{Int,Int}

A sparse array of deletion events relative to a consensus.
The key is the locus (inclusive) of the deletion; the value is the length.
"""
const DelMap = Dict{Int,Int}

export Maybe, SNPMap, InsMap, DelMap

Base.show(io::IO, m::SNPMap) = show(io, [ k => Char(v) for (k,v) in m ])
Base.show(io::IO, m::InsMap) = show(io, [ k => String(Base.copy(v)) for (k,v) in m ])

include("interval.jl")
include("counter.jl")
include("node.jl")
include("util.jl")
include("block.jl")
include("path.jl")
include("edge.jl")
# include("pool.jl")
include("junction.jl")
include("cmd.jl")
include("minimap.jl")
include("mash.jl")

using .Utility:
    read_fasta, write_fasta, name, columns, log,
    make_consensus, alignment_alleles
using .Nodes
using .Blocks
using .Paths
using .Edges
using .Junctions
using .Intervals
# using .Pool

import .Shell: mash, mafft, havecommand
import .Minimap: PanContigs

export Graph
export Shell, Blocks, Nodes, Utility

export graphs, detransitive!, purge!, prune!, finalize!
export pancontigs
export checkblocks

# ------------------------------------------------------------------------
# graph data structure

"""
    struct Graph
        block    :: Dict{String, Block}
        sequence :: Dict{String, Path}
    end

Representation of a multiple sequence alignment. Alignments of homologous sequences
are stored as blocks. A genome is stored as a path, i.e. a list of blocks.
"""
struct Graph
    block    :: Dict{String,Block}   # uuid      -> block
    sequence :: Dict{String,Path}    # isolation -> path
    # TODO: add edge/junction data structure?
end

include("align.jl")
using .Align

# export file formats
include("gfa.jl")

# --------------------------------
# constructors

"""
    Graph(name::String, sequence::Array{UInt8}; circular=false)

Creates a singleton graph from `sequence`. `name` is assumed to be a unique identifier.
If `circular` is unspecified, the sequence is assumed to be linear.
"""
function Graph(name::String, sequence::Array{UInt8}; circular=false)
    block = Block(sequence)
    path  = Path(name, Node(block); circular=circular)

    append!(block, path.node[1], SNPMap(), InsMap(), DelMap())

    return Graph(
         Dict([pair(block)]),
         Dict([pair(path)]),
         # TODO: more items...
    )
end

"""
    graphs(io::IO; circular=false)

Parse a fasta file from stream `io` and return an array of singleton graphs.
If circular is unspecified, all genomes are assumed to be linear.
"""
function graphs(io::IO; circular=false, upper=false)
    case = upper ? (c::UInt8) -> UInt8(uppercase(Char(c))) : (c::UInt8) -> c
    return [Graph(record.name, case.(record.seq); circular=circular) for record in read_fasta(io)]
end

# --------------------------------
# operators

const Link  = NamedTuple{(:block,:strand),Tuple{Block, Bool}}
const Chain = Array{Link, 1}

# XXX: break into smaller functions: too long
"""
    detransitive!(G::Graph)

Find and remove all transitive edges within the given graph.
A transitive chain of edges is defined to be unambiguous: all
sequences must enter on one edge and leave on another. Thus,
this will not perform paralog splitting.
"""
function detransitive!(G::Graph)
    numisos = count_isolates(values(G.sequence))

    # collect all junctions that transitively pass isolates through
    transitives = Junction[]
    for (j, depth) in junctions(values(G.sequence))
        if (numisos[j.left.block] == numisos[j.right.block] == depth)
            push!(transitives, j)
        end
    end

    rev(l::Link)  = (block=l.block,strand=!l.strand)
    rev(c::Chain) = [rev(b) for b in reverse(c)]

    # build chains by threading consecutive transitive junctions
    # TODO: audit each line carefully
    chain = Dict{Block, Chain}()
    for j in transitives
        if j.left.block ∈ keys(chain) && j.right.block ∈ keys(chain)
            c₁, c₂ = chain[j.left.block], chain[j.right.block]

            c₁ == c₂ && continue

            merged =
                if left(j) == last(c₁) && right(j) == first(c₂)
                    cat(c₁, c₂, dims=1)
                elseif left(j) == last(c₁) && rev(right(j)) == last(c₂)
                    cat(c₁, rev(c₂), dims=1)
                elseif rev(left(j)) == first(c₁) && right(j) == first(c₂)
                    cat(rev(c₁), c₂, dims=1)
                elseif rev(left(j)) == first(c₁) && rev(right(j)) == last(c₂)
                    cat(c₂, c₁, dims=1)
                else
                    error("case not covered")
                end

            for b in first.(merged)
                chain[b] = merged
            end

        elseif j.left.block ∈ keys(chain)
            c₀ = chain[j.left.block]
            if left(j) == c₀[end]
                push!(c₀, right(j))
            elseif rev(j.left.block) == c₀[1]
                pushfirst!(c₀, rev(right(j)))
            else
                error("chains should be linear")
            end
            chain[j.right.block] = c₀

        elseif j.right.block ∈ keys(chain)
            c₀ = chain[j.right.block]
            if right(j) == c₀[end]
                push!(c₀, rev(left(j)))
            elseif right(j) == c₀[1]
                pushfirst!(c₀, left(j))
            else
                error("chains should be linear")
            end
            chain[j.left.block] = c₀

        else
            chain[j.left.block]  = [left(j), right(j)]
            chain[j.right.block] = chain[j.left.block]
        end
    end

    # merge chains into one block
    for c in Set(values(chain))
        isos = numisos[c[1].block]
        @assert all([numisos[C.block] == isos for C in c[2:end]])

        new = Block((s ? b : reverse_complement(b) for (b,s) ∈ c)...)
        for iso ∈ keys(isos)
            replace!(G.sequence[iso], c, new)
        end

        for b ∈ first.(c)
            pop!(G.block, b.uuid)
        end

        G.block[new.uuid] = new
    end
end

"""
    prune!(G::Graph)

Remove all blocks from graph `G` that are not currently used by any extant sequence.
Internal function used during guide tree alignment.
"""
function prune!(G::Graph)
    used = Set(n.block.uuid for p in values(G.sequence) for n in p.node)
    filter!((blk)->first(blk) ∈ used, G.block)
end

"""
    purge!(G::Graph)

Remove all blocks from paths found in graph `G` that have zero length.
Internal function used during guide tree alignment.
"""
function purge!(G::Graph)
    for p in values(G.sequence)
        index = Int[]
        for (i, n) in enumerate(p.node)
            if length(n) == 0
                push!(index, i)
                pop!(n.block, n)
            end
        end
        deleteat!(p.node, index)
    end
end

"""
    keeponly!(G::Graph, names::String...)

Remove all sequences from graph `G` that are passed as variadic parameters `names`.
This will marginalize a graph, i.e. return the subgraph that contains only
isolates contained in `names`
"""
function keeponly!(G::Graph, names::T...) where T <: AbstractString
    nameset = Set(names)
    isolate = collect(keys(G.sequence))
    for name in isolate
        name ∉ nameset || continue
        path = G.sequence[name]
        for node in path.node
            pop!(node.block, node)
        end
        delete!(G.sequence, name)
    end

    uuids = collect(keys(G.block))
    for uuid in uuids
        block = G.block[uuid]
        if depth(block) == 0
            delete!(G.block, uuid)
        end
        # TODO: reconsensus?
    end
end

function checkblocks(G::Graph)
    used   = Set(n.block.uuid for p in values(G.sequence) for n in p.node)
    stored = Set(keys(G.block))

    if used != stored
        @show setdiff(used,stored)
        @show setdiff(stored,used)
        # @infiltrate
        error("bad blocks")
    end
end

# ------------------------------------------------------------------------
# i/o & (de)serialization

Base.show(io::IO, G::Graph) = Base.show(io, (paths=values(G.sequence), blocks=values(G.block)))

pancontigs(G::Graph) = let
    uuid = collect(Base.keys(G.block))
    PanContigs(
        uuid,
        [String(sequence(G.block[id])) for id in uuid],
    )
end

# TODO: can we generalize to multiple individuals
#       equivalent to "highway" detection
function serialize(io::IO, G::Graph)
    if length(G.sequence) != 1
        error("only singleton graphs implemented")
    end

    name = collect(keys(G.sequence))[1]
    seq  = collect(values(G.block))[1].sequence

    write_fasta(io, name, seq)
end

"""
    marshal_fasta(io::IO, G::Graph; opt=nothing)

Serialize graph `G` as a fasta format output stream `io`.
Importantly, this will only serialize the consensus sequences for each block
and not the full multiple sequence alignment.

`opt` is currently ignored. It is kept for signature uniformity for other marshal functions
"""
function marshal_fasta(io::IO, G::Graph; opt=nothing)
    for b in values(G.block)
        write_fasta(io, b.uuid, b.sequence)
    end
end

# XXX: think of a way to break up function but maintain graph-wide node lookup table
"""
    marshal_json(io::IO, G::Graph; opt=nothing)

Serialize graph `G` as a json format output stream `io`.
This is the main storage/exported format for PanGraph.
Currently it is the only format that can reconstruct an in-memory pangraph.

`opt` is currently ignored. It is kept for signature uniformity for other marshal functions
"""
function marshal_json(io::IO, G::Graph; opt=nothing)
    NodeID    = NamedTuple{(:id,:name,:number,:strand), Tuple{String,String,Int,Bool}}
    nodes     = Dict{Node{Block}, NodeID}()
    positions = Dict{Block, Dict{NodeID, Tuple{Int,Int}}}()

    # path serialization
    function dict(p::Path)
        blocks = Array{NodeID}(undef, length(p.node))
        counts = Dict{Block,Int}()

        for (i,node) ∈ enumerate(p.node)
            if node.block ∉ keys(counts)
                counts[node.block] = 1
            end
            blocks[i] = (
                id     = node.block.uuid,
                name   = p.name,
                number = counts[node.block],
                strand = node.strand,
            )
            nodes[node] = blocks[i]
            counts[node.block] += 1

            if node.block ∉ keys(positions)
                positions[node.block] = Dict{NodeID, Tuple{Int,Int}}()
            end
            positions[node.block][blocks[i]] = (p.position[i], p.position[i == length(p.node) ? 1 : i+1] - 1)
        end

        return (
            name     = p.name,
            offset   = p.offset,
            circular = p.circular,
            position = p.position,
            blocks   = blocks,
        )
    end

    # block serialization
    pack(d::SNPMap) = [(k,Char(v)) for (k,v) ∈ d]
    pack(d::InsMap) = [(k,String(copy(v))) for (k,v) ∈ d]
    pack(d::DelMap) = [(k,v) for (k,v) ∈ d]

    strip(id) = (name=id.name,number=id.number,strand=id.strand)
    function dict(b::Block)
        return (
            id        = b.uuid,
            sequence  = String(sequence(b)),
            gaps      = b.gaps,
            mutate    = [(strip(nodes[key]), pack(val)) for (key,val) ∈ b.mutate],
            insert    = [(strip(nodes[key]), pack(val)) for (key,val) ∈ b.insert],
            delete    = [(strip(nodes[key]), pack(val)) for (key,val) ∈ b.delete],
            positions = [(strip(key), val) for (key,val) ∈ positions[b]]
        )
    end

    # NOTE: paths must come first as it fills the node lookup table
    paths  = [ dict(path)  for path  ∈ values(G.sequence) ]
    blocks = [ dict(block) for block ∈ values(G.block) ]

    JSON.print(io, (
        paths  = paths,
        blocks = blocks,
    ))
end

# NOTE: only recognizes json input right now
"""
    unmarshal(io::IO)

Deserialize the json formatted input stream `io` into a Graph data structure.
Return a `Graph` type.
"""
function unmarshal(io)
    graph = JSON.parse(io)

    unpack = (
        snp = Dict(),
        ins = Dict(),
        del = Dict(),
    )
    blocks = Dict(map(graph["blocks"]) do blk
        # type wrangling
        b = (
            id       = String(blk["id"]),
            sequence = Array{UInt8}(blk["sequence"]),
            gaps     = Dict{Int,Int}(
                parse(Int,k) => v for (k,v) ∈ blk["gaps"]
            ),
            mutate   = Dict(k=>v for (k,v) ∈ blk["mutate"]),
            insert   = Dict(k=>v for (k,v) ∈ blk["insert"]),
            delete   = Dict(k=>v for (k,v) ∈ blk["delete"]),
        )

        unpack.snp[b.id] = b.mutate
        unpack.ins[b.id] = b.insert
        unpack.del[b.id] = b.delete

        b.id => Block(
            b.id,
            b.sequence,
            b.gaps,
            # empty until we build the required node{block} objects
            Dict{Node{Block},SNPMap}(),
            Dict{Node{Block},InsMap}(),
            Dict{Node{Block},DelMap}(),
        )
    end)

    paths = Dict(map(graph["paths"]) do path
        # type wrangling
        p = (
            name     = String(path["name"]),
            offset   = path["offset"],
            circular = path["circular"],
            blocks   = path["blocks"],
            position = "position" ∈ keys(path) ? path["position"] : [], # to ease upgrade from old pangraph files
        )

        nodes = Node{Block}[]
        sizehint!(nodes, length(p.blocks))

        for (i,blk) ∈ enumerate(p.blocks)
            b = (
                id     = String(blk["id"]),
                strand = blk["strand"],
            )
            push!(nodes, Node(blocks[b.id], b.strand))

            # fill in block variant dictionaries
            filter!(p->p.first!="id", blk)

            blocks[b.id].mutate[nodes[i]] = Dict(
                snp[1] => UInt8(snp[2][1]) for snp ∈ unpack.snp[b.id][blk]
            )

            blocks[b.id].insert[nodes[i]] = Dict(
                let
                    @assert length(ins[1]) == 2
                    key = (ins[1][1], ins[1][2])
                    key => Array{UInt8}(ins[2])
                end for ins ∈ unpack.ins[b.id][blk]
            )

            blocks[b.id].delete[nodes[i]] = Dict(
                del[1] => del[2] for del ∈ unpack.del[b.id][blk]
            )
        end

        path = Path(
            p.name,
            nodes,
            p.offset,
            p.circular,
            p.position,
        )
        length(path.position) > 0 || positions!(path)

        p.name => path
    end)

    return Graph(blocks, paths)
end

# ------------------------------------------------------------------------
# operators

"""
    sequence(G::Graph, name::String)

Return the sequence corresponding to genome `name` within graph `G`
"""
function sequence(g::Graph, name::AbstractString)
    name ∈ keys(g.sequence) || error("'$name' not a valid sequence identifier: have $(collect(keys(g.sequence)))")
    path = g.sequence[name]

    return sequence(path)
end

"""
    sequence(G::Graph)

Return all pairs of `name` => `sequence` encoded within graph `G`
"""
sequence(g::Graph) = [ name => join(String(sequence(node.block, node)) for node ∈ path.node) for (name, path) ∈ g.sequence ]

"""
    realign!(G::Graph; accept)

Realign blocks contained within graph `G`.
Usage of this function requires [MAFFT](https://mafft.cbrc.jp/alignment/software/source.html) to be on the system **PATH**
`accept` should be a function that returns true on blocks you wish to realign.
By default, all blocks are realigned.
"""
function realign!(g::Graph; accept=(_)->true, case=false)
    meter = Progress(length(g.block); desc="polishing progress", output=stderr)
    Threads.@threads for blk in collect(values(g.block))
        if !accept(blk)
            next!(meter)
            continue
        end
        io, node = mafft(blk, case)

        seq = collect(read_fasta(io))
        aln = reduce(hcat, map((r)->r.seq, seq))
        ref = make_consensus(aln)

        nodes = map((r)->node[r.name], seq)
        blk.gaps, blk.mutate, blk.delete, blk.insert, blk.sequence = alignment_alleles(ref, aln, nodes)
        next!(meter)
    end
end

"""
    finalize!(G::Graph)

Compute the position of the breakpoints for each homologous alignment across all sequences within Graph `G`.
Intended to be ran after multiple sequence alignment is complete
"""
function finalize!(g)
    for p in values(g.sequence)
        positions!(p)
    end
end

# ------------------------------------------------------------------------
# main point of entry

"""
    test(path)

Align all sequences found in the fasta file at `path` into a pangraph.
Verifies that after the alignment is complete, all sequences are correctly reconstructed
"""
function test(file="data/marco/mycobacterium_tuberculosis/genomes.fa")
    open = endswith(file,".gz") ? GZip.open : Base.open

    log("> running graph test...")
    log("-> building graph...")

    sequences = String[]
    energy = (aln) -> let
        len = aln.length
        len < 100 && return Inf

        cuts(hit) = (hit.start > 100) + ((hit.length-hit.stop) > 100)

        ncuts = cuts(aln.qry)+cuts(aln.ref)
        nmuts = aln.divergence*aln.length

        return -len + 100*ncuts + 20*nmuts
    end

    graph, isolates = open(file, "r") do io
        isolates  = graphs(io; circular=true, upper=true)
        sequences = [first(sequence(iso)) for iso in isolates]

        reference = Dict(first.(sequence.(isolates)))

        println("-->aligning...")
        align(isolates...;reference=reference,energy=energy,minblock=100), isolates
    end

    log("-> verifying graph...")
    for (i, isolate) ∈ enumerate(isolates)
        name, seq₀ = sequences[i]
        seq₁ = sequence(graph, name)
        if !all(seq₀ .== seq₁)
            path = graph.sequence[name]
            x    = findfirst(seq₁, "ACTTGGCTATCCCGCAGGAC")

            println("> true ($(length(seq₀))):          ", seq₀[1:20])
            println("> reconstructed ($(length(seq₁))): ", seq₁[1:20])
            println("> offset:                          $(path.offset)")
            println("> needed offset:                   $(x)")
            log("--> isolate '$name' incorrectly reconstructed")
        else
            log("--> isolate '$name' correctly reconstructed")
        end
    end

    finalize!(graph)
    graph
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using Omega

function Y(rng, n)
  [normal(rng[@uid], 0.0, 1.0) for i = 1:n]
end

# function test()
#   function x_(ω)
#     a = randexp(ω)
#     a2 = randexp(ω, (1, 2, 3))
#     b = randcycle(ω, 6)
#     c = randn(ω)
#     c2 = randn(ω, (1, 2, 3))
#     d = randstring(ω, 'a':'z', 6)
#     e = randsubseq(ω, collect(1:8), 0.3)
#     f = randperm(ω, 10)
#     (a, a2, b, c, c2, d, e, f)
#   end
#   x = ciid(x_)
#   lω = rand(defΩ())
#   x(lω)
# end

# test()

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using StatsBase
function randeig(t, n; xlims=(-2.1*√n, 2.1*√n), bins=100)
    const xgrid = linspace(xlims..., bins)
    counts = @parallel (+) for _ = 1:t
        fit(Histogram, eigvals(Symmetric(randn(n, n))), xgrid).weights
    end
    Histogram(xgrid, counts)
end


text:
function load_expokit()
return quote

using .Expokit: expmv

end end  # quote / load_expokit

eval(load_expokit())


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

mutable struct ExpressRouteCircuitsArpTableListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{ExpressRouteCircuitArpTable} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function ExpressRouteCircuitsArpTableListResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(ExpressRouteCircuitsArpTableListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(ExpressRouteCircuitsArpTableListResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type ExpressRouteCircuitsArpTableListResult

const _property_map_ExpressRouteCircuitsArpTableListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_ExpressRouteCircuitsArpTableListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{ExpressRouteCircuitArpTable}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ ExpressRouteCircuitsArpTableListResult }) = collect(keys(_property_map_ExpressRouteCircuitsArpTableListResult))
Swagger.property_type(::Type{ ExpressRouteCircuitsArpTableListResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ExpressRouteCircuitsArpTableListResult[name]))}
Swagger.field_name(::Type{ ExpressRouteCircuitsArpTableListResult }, property_name::Symbol) =  _property_map_ExpressRouteCircuitsArpTableListResult[property_name]

function check_required(o::ExpressRouteCircuitsArpTableListResult)
    true
end

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


text:
__precompile__()

module StringCases

  # ============
  #  public api
  # ============

  export upcase, downcase

  export spacecase

  export camelcase, pascalcase
  export snakecase, kebabcase

  export underscore, dasherize

  export camelize, decamelize

  export capitalize, classify

  # ===========
  #  core load
  # ===========

  const upcase = uppercase
  const downcase = lowercase

  include("public/index.jl")
  include("private/index.jl")

  const decamelize = _decamelize

end


text:
using Base.Test
# waveletlayer.jl
import WaveletScattering: WaveletLayer
# bank1d.jl
import WaveletScattering: NullBank, Bank1D
# spec1d.jl
import WaveletScattering: Spec1D

bank = Bank1D(Spec1D())
bottoms = [:fourier]
name = "wavelets"
tops = [:wavelets]
waveletlayer = WaveletLayer(bank=bank,
    bottoms=bottoms, name=name, tops=tops)
@test !Mocha.can_do_bp(waveletlayer)
@test !Mocha.has_neuron(waveletlayer)
@test !Mocha.has_param(waveletlayer)

@test_throws AssertionError WaveletLayer(bank=NullBank(),
    bottoms=bottoms, name=name, tops=tops)
@test_throws AssertionError WaveletLayer(bank=bank,
    bottoms=Symbol[], name=name, tops=tops)
@test_throws AssertionError WaveletLayer(bank=bank,
    bottoms=bottoms, name=name, tops=[:top1, :top2])
@test_throws AssertionError WaveletLayer(bank=bank,
    bottoms=[:bottom1, :bottom2], name=name, tops=tops)


text:
struct NonstationaryVarianceWs
    rΣ_s_s::Matrix{Float64}
    rA1::Matrix{Float64}
    rA2::Matrix{Float64}
    rB1::Matrix{Float64}
    rB2::Matrix{Float64}
    rA2S::Matrix{Float64}
    rB2S::Matrix{Float64}
    rΣ_ns_s::Matrix{Float64}
    rΣ_ns_ns::Matrix{Float64}
    state_stationary_variables::Vector{Bool}
    nonstate_stationary_variables::Vector{Bool}
    function NonstationaryVarianceWs(endogenous_nbr::Int64,
                                        exogenous_nbr::Int64,
                                        state_nbr::Int64,
                                        nonstationary_variables::Vector{Bool},
                                        A2::Matrix{Float64})
        nonstate_nbr = endogenous_nbr - state_nbr
        state_stationary_variables = Vector{Bool}(undef, state_nbr)
        nonstate_stationary_variables = Vector{Bool}(undef, nonstate_nbr)
        state_stationary_variables .= .!nonstationary_variables
        state_stationary_nbr = count(state_stationary_variables)
        fill!(nonstate_stationary_variables, true)
        for i = 1:nonstate_nbr
            for j = 1:state_nbr
                if nonstationary_variables[j] && abs(A2[i, j]) > 1e-10
                    nonstate_stationary_variables[i] = false
                    break
                end
            end
        end
        nonstate_stationary_nbr = count(nonstate_stationary_variables)
        rΣ_s_s = Matrix{Float64}(undef, state_stationary_nbr, state_stationary_nbr)
        rA1 = Matrix{Float64}(undef, state_stationary_nbr, state_stationary_nbr)
        rA2 = Matrix{Float64}(undef, nonstate_stationary_nbr, state_stationary_nbr)
        rB1 = Matrix{Float64}(undef, state_stationary_nbr, exogenous_nbr)
        rB2 = Matrix{Float64}(undef, nonstate_stationary_nbr, exogenous_nbr)
        rA2S = Matrix{Float64}(undef, nonstate_stationary_nbr, state_stationary_nbr)
        rB2S = Matrix{Float64}(undef, nonstate_stationary_nbr, exogenous_nbr)
        rΣ_ns_s = Matrix{Float64}(undef, nonstate_stationary_nbr, state_stationary_nbr)
        rΣ_ns_ns = Matrix{Float64}(undef, nonstate_stationary_nbr, nonstate_stationary_nbr)
        new(rΣ_s_s, rA1, rA2, rB1, rB2,
            rA2S, rB2S, rΣ_ns_s, rΣ_ns_ns,
            state_stationary_variables,
            nonstate_stationary_variables)
    end
end

struct VarianceWs
    B1S::Matrix{Float64}
    B1SB1::Matrix{Float64}
    A2S::Matrix{Float64}
    B2S::Matrix{Float64}
    Σ_s_s::Matrix{Float64}
    Σ_ns_s::Matrix{Float64}
    Σ_ns_ns::Matrix{Float64}
    stationary_variables::Vector{Bool}
    nonstationary_ws::Vector{NonstationaryVarianceWs}
    lre_ws::LinearRationalExpectationsWs
    lyapd_ws::LyapdWs
    function VarianceWs(var_nbr::Int64, state_nbr::Int64,
                           shock_nbr::Int64, lre_ws::LinearRationalExpectationsWs)
        nonstate_nbr = var_nbr - state_nbr
        B1S = Matrix{Float64}(undef, state_nbr, shock_nbr)
        B1SB1 = Matrix{Float64}(undef, state_nbr, state_nbr)
        A2S = Matrix{Float64}(undef, nonstate_nbr, state_nbr)
        B2S = Matrix{Float64}(undef, nonstate_nbr, shock_nbr)
        Σ_s_s = Matrix{Float64}(undef, state_nbr, state_nbr)
        Σ_ns_s = Matrix{Float64}(undef, nonstate_nbr, state_nbr)
        Σ_ns_ns = Matrix{Float64}(undef, nonstate_nbr, nonstate_nbr)
        stationary_variables = Vector{Bool}(undef, var_nbr)
        state_stationary_variables = Vector{Bool}(undef, 0)
        nonstate_stationary_variables = Vector{Bool}(undef, 0)
        nonstationary_ws = Vector{NonstationaryVarianceWs}(undef, 0)
        lyapd_ws = LyapdWs(state_nbr)
        new(B1S, B1SB1, A2S, B2S, Σ_s_s, Σ_ns_s, Σ_ns_ns,
            stationary_variables,
            nonstationary_ws,
            lre_ws, lyapd_ws)
    end
end

function make_stationary_variance!(Σy::Matrix{Float64},
                                   Σ_s_s::Matrix{Float64},
                                   Σ_ns_s::Matrix{Float64},
                                   Σ_ns_ns::Matrix{Float64},
                                   backward_indices)
    n = size(Σy, 1)
    nb = length(backward_indices)
    m1 = m2 = 1
    @inbounds for i = 1:n
        if m1 <= nb && i == backward_indices[m1]
            k1 = k2 = 1
            for j = 1:n
                if k1 <= nb && j == backward_indices[k1]
                    Σy[j,i] = Σ_s_s[k1, m1]
                    k1 += 1
                else
                    Σy[j,i] = Σ_ns_s[k2, m1]
                    k2 += 1
                end
            end
            m1 += 1
        else
            k1 = k2 = 1
            for j = 1:n
                if k1 <= nb && j == backward_indices[k1]
                    Σy[j,i] = Σ_ns_s[m2, k1]
                    k1 += 1
                else
                    Σy[j,i] = Σ_ns_ns[k2, m2]
                    k2 += 1
                end
            end
            m2 += 1

        end
    end
end

function make_nonstationary_variance!(Σy::Matrix{Float64},
                                      Σ_s_s::Matrix{Float64},
                                      Σ_ns_s::Matrix{Float64},
                                      Σ_ns_ns::Matrix{Float64},
                                      state_indices,
                                      state_stationary_indices,
                                      nonstate_indices,
                                      nonstate_stationary_indices)
    n = size(Σy, 1)
    fill!(Σy, NaN)
    m1 = 1
    for (m2, i) in enumerate(state_indices)
        if state_stationary_indices[m2]
            k1 = k2 = 1
            for j in state_indices
                if state_stationary_indices[k2]
                    Σy[j,i] = Σ_s_s[k1, m1]
                    k1 += 1
                end
                k2 += 1
            end
            k1 = k2 = 1
            for j in nonstate_indices
                if nonstate_stationary_indices[k2]
                    Σy[j,i] = Σ_ns_s[k1, m1]
                    k1 += 1
                end
                k2 += 1
            end
            m1 += 1
        end
    end
    m1 = 1
    for (m2, i) in enumerate(nonstate_indices)
        if nonstate_stationary_indices[m2]
            k1 = k2 = 1
            for j in state_indices
                if state_stationary_indices[k2]
                    Σy[j,i] = Σ_ns_s[m1, k1]
                    k1 += 1
                end
                k2 += 1
            end
            k1 = k2 = 1
            for j in nonstate_indices
                if nonstate_stationary_indices[k2]
                    Σy[j,i] = Σ_ns_ns[k1, m1]
                    k1 += 1
                end
                k2 += 1
            end
            m1 += 1
        end
    end
end

function compute_variance!(Σy::Matrix{Float64},
                           A1::AbstractVecOrMat{Float64},
                           A2::AbstractVecOrMat{Float64},
                           B1::AbstractVecOrMat{Float64},
                           B2::AbstractVecOrMat{Float64},
                           Σe::AbstractVecOrMat{Float64},
                           ws::VarianceWs)
    n = size(Σy, 1)
    lre_ws = ws.lre_ws
    mul!(ws.B1S, B1, Σe)
    mul!(ws.B1SB1, ws.B1S, transpose(B1))
    # state variables variance
    extended_lyapd!(ws.Σ_s_s, A1, ws.B1SB1, ws.lyapd_ws)
    if is_stationary(ws.lyapd_ws)
        stationary_variance_blocks!(ws.Σ_ns_s, ws.Σ_ns_ns, A1, A2, B1,
                                    B2, ws.A2S, ws.B2S, ws.Σ_s_s, Σe)
        make_stationary_variance!(Σy, ws.Σ_s_s, ws.Σ_ns_s,
                                  ws.Σ_ns_ns, lre_ws.backward_indices)
        fill!(ws.stationary_variables, true)
    else
        state_nbr = lre_ws.backward_nbr
        state_indices = lre_ws.backward_indices
        if length(ws.nonstationary_ws) == 0
            nonstationary_ws = NonstationaryVarianceWs(lre_ws.endogenous_nbr,
                                                       lre_ws.exogenous_nbr,
                                                       lre_ws.backward_nbr,
                                                       ws.lyapd_ws.nonstationary_variables,
                                                       A2)
            push!(ws.nonstationary_ws, nonstationary_ws)
        else
            nonstationary_ws = ws.nonstationary_ws[1]
        end
        state_stationary_variables = nonstationary_ws.state_stationary_variables
        nonstate_stationary_variables = nonstationary_ws.nonstate_stationary_variables
        fill!(ws.stationary_variables, false)
        m1 = m2 = 1
        for i = 1:size(Σy, 1)
            if m1 <= state_nbr && i == state_indices[m1]
                if state_stationary_variables[m1]
                    ws.stationary_variables[i] = true
                end
                m1 += 1
            else
                if nonstate_stationary_variables[m2]
                    ws.stationary_variables[i] = true
                end
                m2 += 1
            end
        end
        rΣ_ns_s = nonstationary_ws.rΣ_ns_s
        rΣ_ns_ns = nonstationary_ws.rΣ_ns_ns
        rA1    = nonstationary_ws.rA1
        rA2    = nonstationary_ws.rA2
        rB1    = nonstationary_ws.rB1
        rB2    = nonstationary_ws.rB2
        rA2S   = nonstationary_ws.rA2S
        rB2S   = nonstationary_ws.rB2S
        rΣ_s_s = nonstationary_ws.rΣ_s_s
        rΣ_s_s .= view(ws.Σ_s_s, state_stationary_variables, state_stationary_variables)
        rA1 .= view(A1, state_stationary_variables, state_stationary_variables)
        rB1 .= view(B1, state_stationary_variables, :)
        if any(nonstate_stationary_variables)
            rA2 .= view(A2, nonstate_stationary_variables, state_stationary_variables)
            rB2 .= view(B2, nonstate_stationary_variables, :)
        end
        stationary_variance_blocks!(rΣ_ns_s,
                                    rΣ_ns_ns,
                                    rA1,
                                    rA2,
                                    rB1,
                                    rB2,
                                    rA2S,
                                    rB2S,
                                    rΣ_s_s,
                                    Σe)
        make_nonstationary_variance!(Σy,
                                     rΣ_s_s,
                                     rΣ_ns_s,
                                     rΣ_ns_ns,
                                     lre_ws.backward_indices,
                                     state_stationary_variables,
                                     lre_ws.non_backward_indices,
                                     nonstate_stationary_variables)
    end
end

function compute_variance!(
    lreresults::LinearRationalExpectationsResults,
    Σe::AbstractVecOrMat{Float64},
    ws::VarianceWs
)
    A1 = lreresults.gs1
    A2 = lreresults.gns1
    B1 = lreresults.hs1
    B2 = lreresults.hns1
    endogenous_variance = compute_variance!(lreresults.endogenous_variance, A1, A2, B1, B2, Σe, ws)
    empty!(lreresults.stationary_variables)
    append!(lreresults.stationary_variables, ws.stationary_variables)
    return endogenous_variance
end

function compute_variance!(
    variance::AbstractMatrix{Float64},
    lreresults::LinearRationalExpectationsResults,
    Σe::AbstractVecOrMat{Float64},
    ws::VarianceWs
)
    A1 = lreresults.gs1
    A2 = lreresults.gns1
    B1 = lreresults.hs1
    B2 = lreresults.hns1
    endogenous_variance = compute_variance!(variance, A1, A2, B1, B2, Σe, ws)
    return endogenous_variance
end

function stationary_variance_blocks!(Σ_ns_s, Σ_ns_ns, A1, A2, B1, B2, A2S, B2S, Σ_s_s, Σe)
    mul!(B2S, B2, Σe)
    mul!(Σ_ns_s, B2S, transpose(B1))
    mul!(Σ_ns_ns, B2S, transpose(B2))
    mul!(A2S, A2, Σ_s_s)
    mul!(Σ_ns_s, A2S, transpose(A1), 1.0, 1.0)
    mul!(Σ_ns_ns, A2S, transpose(A2), 1.0, 1.0)
end

"""
    correlation!(c::AbstractMatrix{T}, v::AbstractMatrix{T}, sd::AbstractVector{T}) where T <: Real
computes the correlation coefficients c[i,j] = v[i,j]/(sd[i]*sd[j])
"""
function correlation!(c::AbstractMatrix{T}, v::AbstractMatrix{T}, sd::AbstractVector{T}) where T <: Real
    c .= v ./ (sd .* transpose(sd))
end

robustsqrt(x) = sqrt(x + eps())
"""
    correlation(v::AbstractMatrix{T}) where T <: Real
returns the correlation coefficients v[i,j]/sqrt(v[i,i]*v[j,j])
"""
function correlation(v::AbstractMatrix{T}) where T
    sd = robustsqrt.(diag(v))
    c = similar(v)
    correlation!(c, v, sd)
end

"""
autocovariance!(av::Vector{<:AbstractMatrix{T}}, a::AbstractMatrix{T}, v::AbstractMatrix{T}, work1::AbstractMatrix{T}, work2::AbstractMatrix{T},order::Int64)

returns a vector of autocovariance matrices E(y_t y_{t-i}') i = 1,...,i for an vector autoregressive process y_t = Ay_{t-1} + Be_t
"""
function autocovariance!(av::Vector{<:AbstractMatrix{T}},
                         lre_results::LinearRationalExpectationsResults,
                         S1a::AbstractMatrix{T},
                         S1b::AbstractMatrix{T},
                         S2::AbstractMatrix{T},
                         backward_indices::Vector{Int64},
                         stationary_variables::Vector{Bool}) where T <: Real
    S1a .= view(lre_results.endogenous_variance, backward_indices, :)
    n = size(av[1], 1)
    nb = length(backward_indices)
    for p in 1:length(av)
        mul!(S1b, lre_results.gs1, S1a)
        mul!(S2, lre_results.gns1, S1a)
        @inbounds for i = 1:n
            k1 = k2 = 1
            for j = 1:n
                if k1 <= nb && j == backward_indices[k1]
                    av[p][j,i] = S1b[k1, i]
                    k1 += 1
                else
                    av[p][j,i] = S2[k2, i]
                    k2 += 1
                end
            end
        end
        tmp = S1a
        S1a = S1b
        S1b = tmp
    end
    return av
end

"""
autocovariance!(av::Vector{<:AbstractVector{T}}, a::AbstractMatrix{T}, v::AbstractMatrix{T}, work1::AbstractMatrix{T}, work2::AbstractMatrix{T},order::Int64)

returns a vector of autocovariance vector with elements E(y_{j,t} y_{j,t-i}') j= 1,...,n and i = 1,...,i for an vector autoregressive process y_t = Ay_{t-1} + Be_t
"""
function autocovariance!(av::Vector{<:AbstractVector{T}},
                         lre_results::LinearRationalExpectationsResults,
                         S1a::AbstractMatrix{T},
                         S1b::AbstractMatrix{T},
                         S2::AbstractMatrix{T},
                         backward_indices::Vector{Int64},
                         stationary_variables::Vector{Bool}) where T <: Real
    backward_stationary_indices = [i for i in backward_indices if stationary_variables[i]]
    S1a .= view(lre_results.endogenous_variance, backward_stationary_indices, :)
    n = length(av[1])
    nb = length(backward_indices)
    for i in 1:length(av)
        mul!(S1b, lre_results.gs1, S1a)
        mul!(S2, lre_results.gns1, S1a)
        k1 = k2 = 1
        for j = 1:n
            if k1 <= nb && j == backward_indices[k1]
                av[i][j] = S1b[k1, j]
                k1 += 1
            else
                av[i][j] = S2[k2, j]
                k2 += 1
            end
        end
        tmp = S1a
        S1a = S1b
        S1b = tmp
    end
    return av
end

function autocorrelation!(ar::Vector{<:AbstractVecOrMat{T}},
                          lre_results::LinearRationalExpectationsResults,
                          S1a::AbstractMatrix{T},
                          S1b::AbstractMatrix{T},
                          S2::AbstractMatrix{T},
                          backward_indices::Vector{Int64},
                          stationary_variables::Vector{Bool}
                          ) where T <: Real
    autocovariance!(ar, lre_results, S1a, S1b, S2, backward_indices, stationary_variables)
    n = length(stationary_variables)
    dv = Vector{Float64}(undef, n)
    for i = 1:n
        dv[i] = lre_results.endogenous_variance[i, i]
    end
    autocorrelation!(ar, dv)
end

"""
    autocorrelation!(ar::Vector{<:AbstractMatrix{T}}. dv::AbstractVector{T}) where T
returns a vector of autocorrelation matrices corresponging to the autocovariance matrices av
"""
function autocorrelation!(ar::Vector{<:AbstractMatrix{T}}, dv::AbstractVector{T}) where T
    n, m = size(ar[1], 1)
    for a in ar
        for i = 1:n
            for j = 1:m
                a[j, i] ./= sqrt(dv[i]*dv[j])
            end
        end
    end
    return ar
end

"""
    autocorrelation!(ar::Vector{<:AbstractVector{T}}, dv::AbstractVector{T}) where T
returns a vector of autocorrelation vectors corresponging to the autocovariance vector av
"""
function autocorrelation!(ar::Vector{<:AbstractVector{T}}, dv::AbstractVector{T}) where T
    for a in ar
        a ./= dv
    end
    return ar
end

"""
    autocorrelation!(ar::Vector{<:AbstractMatrix{T}}, av::Vector{<:AbstractMatrix{T}}, v::AbstractMatrix{T}) where T
returns a vector of autocorrelation matrices corresponging to the autocovariance matrices av
"""
function autocorrelation!(ar::Vector{<:AbstractMatrix{T}}, av::Vector{<:AbstractMatrix{T}}, v::AbstractMatrix{T}) where T
    n = size(v, 1)
    dv = diag(v)
    for (a1, a2) in zip(ar, av)
        for i = 1:n
            for j = 1:m
                a1[i,j] .= a2[i, j] ./ sqrt(dv[i]*dv[j])
            end
        end
    end
    return ar
end

"""
    autocorrelation!(ar::Vector{<:AbstractVector{T}}, av::Vector{<:AbstractMatrix{T}}, v::AbstractMatrix{T}) where T
returns a vector of autocorrelation vectors corresponging to the autocovariance vector av
"""
function autocorrelation!(ar::Vector{<:AbstractVector{T}}, av::Vector{<:AbstractMatrix{T}}, v::AbstractMatrix{T}) where T
    for (a1, a2) in zip(ar, av)
        a1 .= a2 ./ diag(v)
    end
    return ar
end

function variance_decomposition!(
    VD::AbstractMatrix{<:T},
    LRE_results::LinearRationalExpectationsResults,
    Σe::AbstractMatrix{<:T},
    variance::AbstractVector{<:T},
    work1::AbstractMatrix{<:T},
    work2::AbstractMatrix{<:T},
    lre_variance_ws::LinearRationalExpectations.VarianceWs) where T

    nx = size(Σe, 1)
    # force Σe to be positive definite
    cs = cholesky(Σe + 1e-14*I).L
    for i = 1:nx
        fill!(work1, 0.0)
        # use the vector for ith exogenous variable
        vcs = view(cs, :, i)
        work1 .= vcs*transpose(vcs)
        compute_variance!(work2,
                          LRE_results,
                          work1,
                          lre_variance_ws,
                          )
        k = [3, 4, 6]
        A = LRE_results.g1_1
        B = LRE_results.g1_2
        for j = 1:size(VD, 1)
            VD[j, i] = work2[j, j]/variance[j]
        end
    end
    return VD
end

function variance_decomposition(LREresults::LinearRationalExpectationsResults,
                                LREws::LinearRationalExpectationsWs,
                                Σe::Matrix{Float64},
                                endogenous_nbr::Int64,
                                exogenous_nbr::Int64,
                                state_nbr::Int64)
    VD = Matrix{Float64}(undef, endogenous_nbr, exogenous_nbr)
    work1 = Matrix{Float64}(undef, exogenous_nbr, exogenous_nbr)
    work2 = Matrix{Float64}(undef, endogenous_nbr, endogenous_nbr)
    LREvariance_ws = VarianceWs(
        endogenous_nbr,
        state_nbr,
        exogenous_nbr,
        LREws,
    )
    VD = variance_decomposition!(
        VD,
        LREresults,
        Σe,
        diag(LREresults.endogenous_variance),
        work1,
        work2,
        LREvariance_ws
    )
    return VD
end


text:
using Test

using BcdV3
using StaticArrays
using MimuConstants
using LinearAlgebra
using Printf

function get_callback_maxdiff_prev(r_p, T_p, b_p, show_every = 1)

    callback(k::Int, dlogp, log_p, r, T, b) = begin
        if mod(k, show_every) == 0
            r_e = norm(r - r_p, Inf)
            T_e = [norm(T[k] - T_p[k]) for k = 1:length(T)] |> x  -> norm(x, Inf)
            b_e = [norm(b[k] - b_p[k]) for k = 1:length(b)] |> x  -> norm(x, Inf)
            @printf("%15d logp %.3e  dlogp %.5e r_e %.6e T_e %.6e b_e %.6e \n", k, log_p, dlogp,
            r_e, T_e, b_e)

            # Set old
            r_p = copy(r)
            T_p = deepcopy(T)
            b_p = deepcopy(b)
            # @printf(" C_norm %.8e\n", tr(T[4])/3)
        end
    end

    return callback

end
# Ω(x) = BcdV3.Ω(x)
# Noise less simulations
TT = Float64
active_imus = [3,4,5,6, 27,28, 29, 30]
r_all = InitialPosMIMU4444BT()
r_m = r_all[:, active_imus]
r_ref = r_m[:,1]
r_m .-= r_ref
r_all .-= r_ref

g_m = get_normals_icosahedron()
N_orientions = size(g_m, 2)

sig_a = get_rms_noise_acc()
sig_g = get_rms_noise_gyro()
Na = size(r_m, 2)
Ng = 8

N = 40 # Dyn
N_per_orientation = Int(N/5)
Ns = ones(Int, N_orientions).*N_per_orientation
show_every = 100
N_bcd  = 400_000
tol_dlogp = 1.0e-9
pool = [1]
# pool = workers()
Nt = Na + Ng
g_mag = 9.81

g = [SVector{3,TT}(g_m[:,k]) for k = 1:N_orientions] .|> x -> x/norm(x)*g_mag

b = [(k <= Na ?  0.2 : deg2rad(1.5))*(@SVector randn(3)) for k = 1:Nt]
T = [I + 1.0e-3*@SMatrix randn(3,3) for k = 1:Nt]
T[1] = I + 1.0e-3*SMatrix{3,3, Float64}(triu(randn(3,3)))

r = SMatrix{3,Na, TT}(r_m)

w = [deg2rad(1000.0) * @SVector randn(3) for _ = 1:N]
w_dot = [deg2rad(40e3) * @SVector randn(3) for _ = 1:N]
s = [20.0 * @SVector randn(3) for _ = 1:N]

y_s_true = vcat(
    [T[k]*(-g[n]) + b[k] for k = 1:Na, n = 1:N_orientions],
    [b[k] for k = range(Na+1,length=Ng), n = 1:N_orientions],
)
y_d_true = vcat(
    [T[k]*(Ω(w[n])^2*r[:,k] + Ω(w_dot[n])*r[:,k] + s[n]) + b[k] for k = 1:Na, n = 1:N],
    [T[k]*w[n]+b[k] for k = range(Na+1,length=Ng), n = 1:N]
)

y_s = [y_s_true[k,n] + (k <= Na ? sig_a : sig_g)*randn(3)/sqrt(N_per_orientation) for k = 1:Nt, n = 1:N_orientions]
y_d = [y_d_true[k,n] + ( k <= Na ? sig_a : sig_g)*randn(3) for k = 1:Nt, n = 1:N]

Q_inv = [SMatrix{3,3,Float64}(I)/( k <= Na ? sig_a^2 :  sig_g^2) for k = 1:Nt]

b0 = [b_n + 0.0e-2*randn(3) for  b_n in b]
T0 = [T[k] + 0.0e-3*(k == 1 ? triu(randn(3,3)) : randn(3,3))  for k = 1:Nt]
r0_m = copy(r_m)
r0_m[:,2:end] += 0.0e-3*randn(3,Na-1)
r0 = SMatrix{3,Na}(r0_m)

println("BCD v2")
r_hat, T_hat, b_hat = @time bcd(
    r0, T0, b0, y_s, y_d, Q_inv, Ns, Na , Ng, g_mag, tol_dlogp, N_bcd,
    1.0e-12, 100, 1.0e-13, 200, pool;
    # callback = callback
    callback = get_callback_maxdiff_prev(r, T, b, show_every) #callback
)
println("True")
display(r |> Array)
println("Hat")
display(r_hat |> Array)
println("r0")
display(r0 |> Array)
println("diff hat true")
display(r_hat - r |> Array .|> abs)
println("diff hat 0")
display(r_hat - r0 |> Array .|> abs)
println("sse r_hat")
display(sum((r_hat - r).^2))


text:
module Abstract

using IRTools, MacroTools
using IRTools.All
using IRTools: block

abstract(Ts...) = nothing
partial(Ts...) = abstract(Ts...)

include("utils.jl")
include("cleanup.jl")
include("flow.jl")
include("trace.jl")
include("lib/base.jl")
include("lib/array.jl")
include("lib/struct.jl")
include("lib/random.jl")

end


text:
using FinEtools
using FinEtools.TetRemeshingModule
using NIfTI

File = "Head_256x256x126-256x256x252"
file = niread(File * ".nii")
bs = voxel_size(file.header) .* [size(file,1), size(file,2), size(file,3)-119]
vb = VoxelBoxVolume(file.raw[:, :, 120:end, 1], bs)
vtkexport(File * ".vtk", vb)

fens, fes = T4voximg(vb.data, vec([voxeldims(vb)...]), eltype(vb.data)[255])
t = deepcopy(fes.conn);
v = deepcopy(fens.xyz);
tmid = ones(Int, size(t,1));

desired_ts = 4 * voxel_size(file.header)[1];
bfes = meshboundary(fes);
f = connectednodes(bfes);
bv = zeros(Bool, size(v,1));
# bv[f] = true;

println("Mesh size: initial = $(size(t,1))")
t0 = time()

t, v, tmid = TetRemeshingModule.coarsen(t, v, tmid; bv = bv, desired_ts = desired_ts);

println("Mesh size: final = $(size(t,1)) [$(time() - t0) sec]")

fens.xyz = deepcopy(v)
fes.conn = deepcopy(t)
setlabel!(fes, tmid)

fens = meshsmoothing(fens, fes; npass = 5)
println("count(fes) = $(count(fes))")

File = "Head_mesh.vtk"
vtkexportmesh(File, fens, fes)
@async run(`"paraview.exe" $File`)


text:
# Inside make.jl
push!(LOAD_PATH,"../src/")
using EmpiricalMeasures
using Documenter
makedocs(
         sitename="EmpiricalMeasures.jl",
         modules=[EmpiricalMeasures],
         pages=[
                "Home" => "index.md"
               ])
deploydocs(;
    repo="github.com/davibarreira/EmpiricalMeasures.jl",
)


text:
const shortnames = Dict(
    (".tif", ".tiff") => "GTiff",
    (".img",) => "HFA",
    (".xyz",) => "XYZ",
)

function find_shortname(fn::AbstractString)
    _, ext = splitext(fn)
    for (k, v) in shortnames
        if ext in k
            return v
        end
    end
    error("Cannot determine GDAL Driver for $fn")
end

const GMF = Dict(
    :GMF_ALL_VALID => 0x01,
    :GMF_PER_DATASET => 0x02,
    :GMF_ALPHA => 0x04,
    :GMF_NODATA => 0x08,
    )

"""Takes bitwise OR-ed set of status flags and returns flags."""
function mask_flags(flags::Int32)
    (f.first for f in GMF if (flags & f.second) == f.second)
end
mask_flags(band::ArchGDAL.AbstractRasterBand) = mask_flags(ArchGDAL.GDAL.gdalgetmaskflags(Ptr{Nothing}(band.ptr)))

"""Retrieves nodata value from RasterBand."""
function get_nodata(band::Ptr{Nothing})
    succes_value = Int32(0)
    nodata = ArchGDAL.GDAL.gdalgetrasternodatavalue(band, Ref(succes_value))
    if succes_value == 0
        return nodata
    else
        @warn "Unsuccessful in getting nodata."
        return nothing
    end
end
get_nodata(band::ArchGDAL.AbstractRasterBand) = get_nodata(Ptr{Nothing}(band.ptr))

const gdt_conversion = Dict{DataType,DataType}(
    Bool => UInt8,
    Int8 => UInt8,
    UInt64 => UInt32,
    Int64 => Int32
)

"""Converts type of Array for one that exists in GDAL."""
function cast_to_gdal(A::AbstractArray{<:Real,3})
    type = eltype(A)
    if type in keys(gdt_conversion)
        newtype = gdt_conversion[type]
        @warn "Casting $type to $newtype to fit in GDAL."
        return newtype, convert(Array{newtype}, A)
    else
        error("Can't cast $(eltype(A)) to GDAL.")
    end
end


text:
@enum Version begin

  PEGGED_TO_BENCHMARK      = 102
  MODELS_SUPPORT           = 103
  SEC_DEF_OPT_PARAMS_REQ   = 104
  EXT_OPERATOR             = 105
  SOFT_DOLLAR_TIER         = 106
  REQ_FAMILY_CODES         = 107
  REQ_MATCHING_SYMBOLS     = 108
  PAST_LIMIT               = 109
  MD_SIZE_MULTIPLIER       = 110
  CASH_QTY                 = 111
  REQ_MKT_DEPTH_EXCHANGES  = 112
  TICK_NEWS                = 113
  REQ_SMART_COMPONENTS     = 114
  REQ_NEWS_PROVIDERS       = 115
  REQ_NEWS_ARTICLE         = 116
  REQ_HISTORICAL_NEWS      = 117
  REQ_HEAD_TIMESTAMP       = 118
  REQ_HISTOGRAM            = 119
  SERVICE_DATA_TYPE        = 120
  AGG_GROUP                = 121
  UNDERLYING_INFO          = 122
  CANCEL_HEADTIMESTAMP     = 123
  SYNT_REALTIME_BARS       = 124
  CFD_REROUTE              = 125
  MARKET_RULES             = 126
  PNL                      = 127
  NEWS_QUERY_ORIGINS       = 128
  UNREALIZED_PNL           = 129
  HISTORICAL_TICKS         = 130
  MARKET_CAP_PRICE         = 131
  PRE_OPEN_BID_ASK         = 132
  REAL_EXPIRATION_DATE     = 134
  REALIZED_PNL             = 135
  LAST_LIQUIDITY           = 136
  TICK_BY_TICK             = 137
  DECISION_MAKER           = 138
  MIFID_EXECUTION          = 139
  TICK_BY_TICK_IGNORE_SIZE = 140
  AUTO_PRICE_FOR_HEDGE     = 141
  WHAT_IF_EXT_FIELDS       = 142
  SCANNER_GENERIC_OPTS     = 143
  API_BIND_ORDER           = 144
  ORDER_CONTAINER          = 145
  SMART_DEPTH              = 146
  REMOVE_NULL_ALL_CASTING  = 147
  D_PEG_ORDERS             = 148
  MKT_DEPTH_PRIM_EXCHANGE  = 149
  COMPLETED_ORDERS         = 150
  PRICE_MGMT_ALGO          = 151
  STOCK_TYPE               = 152

end


text:
module Discrete
using Printf
using LinearAlgebra
using ..Spectral
using ..Interp
include("kernel.jl")

function ωQR(kernel, rtol, print::Bool = true)
    # print && println(τ.grid[end], ", ", τ.panel[end])
    # print && println(ω.grid[end], ", ", ω.panel[end])
    ################# find the rank ##############################
    """
    For a given index k, decompose R=[R11, R12; 0, R22] where R11 is a k×k matrix.
    If R11 is well-conditioned, then
    σᵢ(R11) ≤ σᵢ(kernel) for 1≤i≤k, and
    σⱼ(kernel) ≤ σⱼ₋ₖ(R22) for k+1≤j≤N
    See Page 487 of the book: Golub, G.H. and Van Loan, C.F., 2013. Matrix computations. 4th. Johns Hopkins.
    Thus, the effective rank is defined as the minimal k that satisfy rtol≤ σ₁(R22)/σ₁(kernel)
    """
    Nτ, Nω = size(kernel)

    u, σ, v = svd(kernel)
    rank, err = 1, 0.0
    for (si, s) in enumerate(σ)
        # println(si, " => ", s / σ[1])
        if s / σ[1] < rtol
            rank = si - 1
            err = s[1] / σ[1]
            break
        end
    end
    print && println("Kernel ϵ-rank = ", rank, ", rtol ≈ ", err)

    Q, R, p = qr(kernel, Val(true)) # julia qr has a strange design, Val(true) will do a pivot QR
    # size(R) == (Nτ, Nω) if Nω>Nτ
    # or size(R) == (Nω, Nω) if Nω<Nτ

    for idx = rank:min(Nτ, Nω)
        if Nω > Nτ
            R22 = R[idx:Nτ, idx:Nω]
        else
            R22 = R[idx:Nω, idx:Nω]
        end
        u2, s2, v2 = svd(R22)
        # println(idx, " => ", s2[1] / σ[1])
        if s2[1] / σ[1] < rtol
            rank = idx
            err = s2[1] / σ[1]
            break
        end
    end
    print && println("DLR rank      = ", rank, ", rtol ≈ ", err)

    # @assert err ≈ 4.58983288255442e-13

    return p[1:rank]
end

function τnQR(kernel, rank, print::Bool = true)
    ###########  dlr grid for τ  ###################
    print && println("Calculating τ and ωn grid ...")
    @assert rank == size(kernel)[2] #the ω dimension of the τkernel should be the effective rank
    τnqr = qr(transpose(kernel), Val(true)) # julia qr has a strange, Val(true) will do a pivot QR
    return τnqr.p[1:rank]
end

function buildωn(dlrGrid, print::Bool = true)
    ###########  dlr grid for ωn  ###################
    print && println("Calculating ωn grid ...")
    symmetry = dlrGrid.symmetry
    rank = dlrGrid.size

    if symmetry == :ph || symmetry == :pha
        Nωn = Int(ceil(Λ)) * 2 # expect Nω ~ para.Λ/2π, drop 2π on the safe side
        ωnkernel = zeros(Float64, (rank, Nωn + 1))
        ωnGrid = [w for w = 0:Nωn]
        # fermionic Matsubara frequency ωn=(2n+1)π for type==:acorr
        # bosonic Matsubara frequency ωn=2nπ for type==:corr
    else
        Nωn = Int(ceil(Λ)) * 2 # expect Nω ~ para.Λ/2π, drop 2π on the safe side
        ωnkernel = zeros(Complex{Float64}, (rank, 2Nωn + 1))
        ωnGrid = [w for w = -Nωn:Nωn] # fermionic Matsubara frequency ωn=(2n+1)π
    end

    for (ni, n) in enumerate(ωnGrid)
        for r = 1:rank
            ωnkernel[r, ni] = Spectral.kernelΩ(type, n, ωGridDLR[r])
        end
    end
    nqr = qr(ωnkernel, Val(true)) # julia qr has a strange, Val(true) will do a pivot QR
    nGridDLR = sort(ωnGrid[nqr.p[1:rank]])
    return nGridDLR, nqr.p[1:rank]
end

"""
function build(dlrGrid, print::Bool = true)
    Construct discrete Lehmann representation

#Arguments:
- `dlrGrid`: struct that contains the information to construct the DLR grid. The following entries are required:
   Λ: the dimensionless scale β*Euv, rtol: the required relative accuracy, isFermi: fermionic or bosonic, symmetry: particle-hole symmetry/antisymmetry or none
- `print`: print the internal information or not
"""
function build(dlrGrid, print::Bool = true)
    degree = 24 # number of Chebyshev nodes in each panel
    Λ = dlrGrid.Λ
    rtol = dlrGrid.rtol
    τ = τChebyGrid(dlrGrid, degree, print)
    ω = ωChebyGrid(dlrGrid, degree, print)
    print && println("τ grid size = $(τ.ngrid)")
    print && println("ω grid size = $(ω.ngrid)")

    kernel = preciseKernelT(dlrGrid, τ, ω, print)
    testInterpolation(dlrGrid, τ, ω, kernel, print)

    ωIndex = ωQR(kernel, rtol, print)
    rank = length(ωIndex)
    ωGrid = sort(ω.grid[ωIndex])

    τIndex = τnQR(kernel[:, ωIndex], rank, print)
    τGrid = sort(τ.grid[τIndex])

    nFineGrid, nFermiKernel, nBoseKernel = preciseKernelΩn(dlrGrid, ωGrid, print)

    nFermiIndex = τnQR(nFermiKernel, rank, print)
    nFermiGrid = sort(nFineGrid[nFermiIndex])

    nBoseIndex = τnQR(nBoseKernel, rank, print)
    nBoseGrid = sort(nFineGrid[nBoseIndex])

    ########### output  ############################
    print && @printf("%5s  %32s  %32s  %11s  %11s\n", "index", "real freq", "tau", "fermi ωn", "bose ωn")
    for r = 1:rank
        print && @printf("%5i  %32.17g  %32.17g  %16i %16i\n", r, ωGrid[r], τGrid[r], nFermiGrid[r], nBoseGrid[r])
    end

    # dlr = Dict([(:ω, ωGrid), (:τ, τGrid), (:ωn, nFermiGrid)])
    return ωGrid, τGrid, nFermiGrid, nBoseGrid
end
end


text:
using Xpress, Base.Test, MathOptInterface, MathOptInterface.Test, MathOptInterfaceXpress
using CPLEX, MathOptInterfaceCPLEX
using Gurobi, MathOptInterfaceGurobi
using GLPK, MathOptInterfaceGLPK

const MOI = MathOptInterface

const MOIU = MathOptInterface.Utilities

solvers = [XpressOptimizer(),
          CPLEXOptimizer(),
          GLPKOptimizerLP(),
        #   GLPKOptimizerMIP(),
          GurobiOptimizer()]

for solver in solvers
    m = Model()
    x = @variable(m)
    y = @variable(m)
    @objective(m, Min, -x)

    c = @constraint(m, x + y <= 1)
    c1 = @constraint(m, y >= 0)
    c2 = @constraint(m, x >= 0)

    MOIU.resetoptimizer!(m, solver)
    MOIU.attachoptimizer!(m)

    JuMP.optimize(m)

    @test JuMP.hasresultvalues(m)

    @test JuMP.terminationstatus(m) == MOI.Success
    @test JuMP.primalstatus(m) == MOI.FEASIBLE_POINT

    JuMP.resultvalue(x)
    JuMP.resultvalue(y)
    JuMP.resultvalue(x + y)
    JuMP.objectivevalue(m)

    JuMP.dualstatus(m)
    JuMP.resultdual(c)
end


text:
"""
	gmtconvert(cmd0::String="", arg1=[], kwargs...)

Convert, Paste, and/or Extract columns from data tables

Full option list at [`gmtconvert`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html)

Parameters
----------

- **A** : **h_cat** : -- Str or [] --

    The records from the input files should be pasted horizontally, not appended vertically [Default].
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#A)
- **C** : **n_records** : -- Str --  Flags = [+lmin][+umax][+i]

    Only output segments whose number of records matches your given criteria:
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#c)
- **D** : **dump** : -- Str or [] --   Flags = [template[+oorig]]

    For multiple segment data, dump each segment to a separate output file.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#d)
- **E** : **first_last** : -- Str or [] --   Flags = [f|l|m|Mstride]

    Only extract the first and last record for each segment of interest.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#e)
- **F** : **conn_method** : -- Str or [] --   Flags = [c|n|r|v][refpoint]

    Alter the way points are connected (by specifying a scheme) and data are grouped (by specifying a method).
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#f)
- **I** : **invert** : -- Str or [] --      Flags = [tsr]

    Invert the order of items, i.e., output the items in reverse order, starting with the last
    and ending up with the first item.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#i)
- **L** : **list_only** : -- Bool or [] --

    Only output a listing of all segment header records and no data records.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#l)
- **N** : **sort** : -- Str or Number --      Flags = [- |+]col

    Numerically sort each segment based on values in column col.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#n)
- **Q** : **select_num** : -- Str --      Flags =  [~]selection

    Only write segments whose number is included in ``selection`` and skip all others.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#q)
- **S** : **select_hdr** : -- Str --      Flags =  [~]”search string” or [~]/regexp/[i]

    Only output those segments whose header record contains the specified text string.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#s)
- **T** : **suppress ** : -- Str or [] --    Flags = [h|d]

    Suppress the writing of certain records on output. Append h to suppress segment headers
    [Default] or d to suppress duplicate data records. Use T=:hd to suppress both types of records.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#t)
- **W** : **gmtconvert** : -- Str or [] --      Flags =  [+n]

    Attempt to gmtconvert each word in the trialing text to a number and append such values
    to the numerical output columns.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#w)
- **Z** : **range** : -- Str or [] --      Flags =  [first][:last]

    Limit output to the specified record range. If first is not set it defaults to record 0
    (very first record) and if last is not set then it defaults to the very last record.
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/gmtconvert.html#z)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_s)
- $(GMT.opt_swap_xy)
"""
function gmtconvert(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("gmtconvert", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_b(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_g(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_o(cmd, d)
	cmd, = parse_s(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'A', d, [:A :h_cat])
	cmd = add_opt(cmd, 'C', d, [:C :n_records])
	cmd = add_opt(cmd, 'D', d, [:D :dump])
	cmd = add_opt(cmd, 'E', d, [:E :first_last])
	cmd = add_opt(cmd, 'F', d, [:F :conn_method])
	cmd = add_opt(cmd, 'I', d, [:I :invert])
	cmd = add_opt(cmd, 'L', d, [:L :list_only])
	cmd = add_opt(cmd, 'N', d, [:N :sort])
	cmd = add_opt(cmd, 'Q', d, [:Q :select_num])
	cmd = add_opt(cmd, 'S', d, [:S :select_hdr])
	cmd = add_opt(cmd, 'T', d, [:T :suppress])
	cmd = add_opt(cmd, 'W', d, [:W :gmtconvert])
	cmd = add_opt(cmd, 'Z', d, [:Z :range])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "gmtconvert", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtconvert(arg1=[], cmd0::String=""; kw...) = gmtconvert(cmd0, arg1; kw...)


text:
__precompile__()

module UnitfulAstro

import Unitful
using Unitful: @unit, @logscale, @logunit

macro import_from_unitful(args...)
    expr = Expr(:block)
    for arg in args
        use_SI_prefixes, sym = should_we_use_SI_prefixes(arg)
        if use_SI_prefixes
            for prefix in Unitful.si_prefixes
                sym′ = Symbol(prefix, sym)
                push!(expr.args, :(import Unitful.$sym′))
            end
        else
            push!(expr.args, :(import Unitful.$sym))
        end
    end
    expr
end

function should_we_use_SI_prefixes(arg::Expr)
    if arg.head == :(call)
        if arg.args[1] == :(~)
            return true, arg.args[2]
        else
            error("incorrect first argument")
        end
    else
        error("incorrect expression head")
    end
end
should_we_use_SI_prefixes(arg::Symbol) = false, arg

@import_from_unitful ~m ~s ~A ~K ~cd ~yr
@import_from_unitful ~L ~Hz ~N ~Pa ~J ~W ~C ~V ~Ω ~S ~F ~H ~T ~Wb ~lm ~lx ~Bq ~Gy ~Sv ~kat ~eV
@import_from_unitful sr rad ° °C °F Ra minute hr d wk ~bar atm Torr
@import_from_unitful q c0 c μ0 µ0 ε0 ϵ0 Z0 G gn ge h ħ Φ0 me mn mp μB µB Na R k σ
@import_from_unitful inch ft yd mi ac lb oz dr gr lbf

# CGS units have been moved to default Unitful.jl now
@import_from_unitful dyn erg

import UnitfulAngles: arcminute, arcsecond

@unit AU         "AU"       AstronomicalUnit          149_597_870_700.0*m       false  # cf IAU 2012
@unit ly         "ly"       LightYear                 1*c*yr                    false
@unit pc         "pc"       Parsec                    1*AU/arcsecond            true   # cf IAU 2015
# Note that IAU 2015 defines the small angle approximation as exact in the definition of the parsec.
@unit Jy         "Jy"       Jansky                    1e-23erg*s^-1*cm^-2*Hz^-1 true
# Note that Jy uses a Float64 conversion factor because 10^23 overflows Int64.

# SOLAR CONVERSION CONSTANTS (IAU 2015)
@unit Rsun       "R⊙"       SolarRadius               6.957e8*m                 false
@unit Ssun       "S⊙"       SolarIrradiance           1361.0*W*m^-2             false
@unit Lsun       "L⊙"       SolarLuminosity           3.828e26*W                false
# TODO solar effective temperature
@unit GMsun      "GM⊙"      GSolarMass                1.327_124_4e20*m^3*s^-2   false
@unit Msun       "M⊙"       SolarMass                 1*GMsun/G                 false

# PLANETARY CONVERSION CONSTANTS (IAU 2015)
@unit Rearth_e   "R⊕ₑ"      EarthEquatorialRadius     6.3781e6*m                false
@unit Rearth_p   "R⊕ₚ"      EarthPolarRadius          6.3568e6*m                false
@unit Rjup_e     "Rjupₑ"    JupiterEquatorialRadius   7.1492e7*m                false
@unit Rjup_p     "Rjupₚ"    JupiterPolarRadius        6.6854e7*m                false
@unit GMearth    "GM⊕"      GEarthMass                3.986_004e14*m^3*s^-2     false
@unit GMjup      "GMjup"    GJupiterMass              1.266_865_3e17*m^3*s^-2   false
@unit Mearth     "M⊕"       EarthMass                 1*GMearth/G               false
@unit Mjup       "Mjup"     JupiterMass               1*GMjup/G                 false
@unit Rearth     "R⊕"       EarthRadius               1*Rearth_e                false
@unit Rjup       "Rjup"     JupiterRadius             1*Rjup_e                  false
# Note that IAU 2015 states that when the radius is not specified as polar or equatorial, the
# equatorial radius should be used.

# "ten to the Fifty-One Ergs" (approximate energy of a supernova explosition)
@unit foe        "foe"      FiftyOneErgs              1e51*erg                  false

# Solar flux unit
# https://en.wikipedia.org/wiki/Solar_flux_unit
@unit SFU        "SFU"      SolarFluxUnit             10*kJy                    false

# Total electron content unit (used for ionospheric physics and low-frequency radio astronomy)
@unit TECU       "TECU"     TotalElectronContentUnit  1e16*m^-2                 false

@logscale mag "mag" Magnitude 10 -2.5 false

#https://en.wikipedia.org/wiki/AB_magnitude
@logunit    AB_mag  "AB mag"        Magnitude   3631Jy

# Bolometric magnitudes (cf IAU 2015)
# https://en.wikipedia.org/wiki/Absolute_magnitude#Bolometric_magnitude
@logunit    bol_Mag "bol Mag"       Magnitude   3.0128e28*W
@logunit    bol_mag "bol mag"       Magnitude   2.518_021_002e-8*W*m^-2

#https://en.wikipedia.org/wiki/Apparent_magnitude#Standard_reference_values
#https://www.cfa.harvard.edu/~dfabricant/huchra/ay145/mags.html
@logunit    U_mag   "Johnson U mag" Magnitude   1810Jy
@logunit    B_mag   "Johnson B mag" Magnitude   4260Jy
@logunit    V_mag   "Johnson V mag" Magnitude   3640Jy
@logunit    R_mag   "Johnson R mag" Magnitude   3080Jy
@logunit    I_mag   "Johnson I mag" Magnitude   2550Jy
@logunit    J_mag   "Johnson J mag" Magnitude   1600Jy
@logunit    H_mag   "Johnson H mag" Magnitude   1080Jy
@logunit    K_mag   "Johnson K mag" Magnitude    670Jy
@logunit    g_mag   "Gunn g mag"    Magnitude   3730Jy
@logunit    r_mag   "Gunn r mag"    Magnitude   4490Jy
@logunit    i_mag   "Gunn i mag"    Magnitude   4760Jy
@logunit    z_mag   "Gunn z mag"    Magnitude   4810Jy

struct MagnitudeError <: Exception
    x
end
Base.showerror(io::IO, e::MagnitudeError) = print(io, "MagnitudeError: ", e.x)

# https://en.wikipedia.org/wiki/Color_index
# for color indices (subtracting magnitudes of different passbands)
import Base.-
function -(x::Unitful.Level{L,S,T},y::Unitful.Level{L,S,T}) where {L<:Magnitude,S,T<:Unitful.Quantity}
    return Unitful.Level{L,S}(-(x.val, y.val))
end
import Base.-
function -(x::Unitful.Level{L,S,T},y::Unitful.Level{L,R,T}) where {L<:Magnitude,S,R,T<:Unitful.Quantity}
    return Unitful.ustrip(x) - Unitful.ustrip(y)
end
import Base.+
function +(x::Unitful.Level{L,S,T},y::Unitful.Level{L,S,T}) where {L<:Magnitude,S,T<:Unitful.Quantity}
    return Unitful.Level{L,S}(+(x.val, y.val))
end
import Base.+
function +(x::Unitful.Level{L,S,T},y::Unitful.Level{L,R,T}) where {L<:Magnitude,S,R,T<:Unitful.Quantity}
    throw(MagnitudeError("an invalid operation was attempted with magnitudes: $x, $y"))
end
import Base.*
function *(x::Unitful.Level{L,S,T},y::Unitful.Level{L,R,T}) where {L<:Magnitude,S,R,T<:Unitful.Quantity}
    throw(MagnitudeError("an invalid operation was attempted with magnitudes: $x, $y"))
end
import Base./
function /(x::Unitful.Level{L,S,T},y::Unitful.Level{L,R,T}) where {L<:Magnitude,S,R,T<:Unitful.Quantity}
    throw(MagnitudeError("an invalid operation was attempted with magnitudes: $x, $y"))
end
import Base.//
function //(x::Unitful.Level{L,S,T},y::Unitful.Level{L,R,T}) where {L<:Magnitude,S,R,T<:Unitful.Quantity}
    throw(MagnitudeError("an invalid operation was attempted with magnitudes: $x, $y"))
end

const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(UnitfulAstro)
end

end


text:
include("../../src/ZZDiffusionBridge.jl")
using LinearAlgebra
"""
    OUSDE <: AbstractModel

dX_t = ν(μ - X_t)dt + dB_t
ν := intensity
μ := mean reversion
M::Array{Float64, 2} := matrix whose element i,j is equal to ∫_0^T ϕ_i ϕ_j dt
V::Vector{Float64} :=  vector whose element i is equal to ∫_0^T ϕ_i dt
bound1::Vector{Float64} :=  vector whose element i is equal to ∫_0^T upbar{ϕ}_i ϕ_j dt
bound2::Vector{Float64} :=  vector whose element i is equal to ∫_0^T downbar{ϕ}_i ϕ_j dt
"""
struct OUSDE <: AbstractModel
    μ::Float64
    ν::Float64
    M::Array{Float64, 2}
    V::Vector{Float64}
    bound1::Vector{Float64}
    bound2::Vector{Float64}
    function OUSDE(μ, ν, L, T)
        new(μ, ν, generate_matrix(L, T), generate_vector(L, T), generate_bound1(L,T), generate_bound2(L,T))
    end
end

dependence_strucute(::OUSDE) = PartialIndependence()
sampling_scheme(::OUSDE) = Regular()

#Poisson rates
"""
    wait_gengaus(a,b,u)
obtaining waiting time for Inhomogeneous Poisson Process
with rate of the form λ(t) = (a + b*t)^+, `a`,`b` ∈ R, `u` random variable
"""
function wait_gengaus(a,b,u)
    if b > 0
        if a < 0
            τ = sqrt(-log(u)*2.0/b) - a/b
        else #a[i]>0
            τ = sqrt((a/b)^2 - log(u)*2.0/b) - a/b
        end
    elseif  b == 0
        if a > 0
            τ = -log(u)/a
        else #a[i] <= 0
            τ = Inf
        end
    else #b[i] < 0
        if a <= 0
            τ = Inf
        elseif -log(u) <= -a^2/b + a^2/(2*b)
            τ = - sqrt((a/b)^2 - log(u)*2.0/b) - a/b
        else
            τ = Inf
        end
    end
end

"""
    λbar(n, S::System, X::OUSDE , u, v)

Poisson time (real bound) for the coefficient `n`
of model `OUSDE` starting at `u` and ending at `v`

λbar = λ No upper bound and no accept reject step
"""
function λbar(n, S::System, X::OUSDE , u::Float64, v::Float64, t::Float64)
    a = S.θ[n]*(X.ν*X.ν*(dot(X.M[n,:], S.ξ) + X.bound1[n]*v + X.bound2[n]*u - X.μ*X.V[n]) + S.ξ[n])
    b = S.θ[n]*((dot(X.M[:,n], S.θ))*X.ν*X.ν + S.θ[n])
    return wait_gengaus(a, b, rand())
end

####
function runall(SHORT = false)
    Random.seed!(1)
    T = 10.0
    clock = 1000.0
    L = 6
    ν = 1.0
    μ = -5
    u = -1.0
    v = 2.0
    X = OUSDE(μ, ν, L, T)
    XX = zz_sampler(X, T, L, u, v, clock)
    if SHORT == false
        burning = 10.0    #burning
        f = clock - 1.0; n = 100
        db = (f-burning)/n
        b =  burning:db:f
        p = plotmixing(XX, b, T, L, u, v)
        hline!(p, [-5.0], color = :blue)
        xaxis!(p, "t")
        yaxis!(p, "X_t")
        display(p)
        #plot the mean of the process
    end
    return XX
end

error("STOP HERE")

runall()
savefig("../../output/ou_1_m5.pdf")

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge