text:
module Neo4j

using Requests
using JSON

export getgraph, version, createnode, getnode, deletenode, setnodeproperty, getnodeproperty,
       getnodeproperties, updatenodeproperties, deletenodeproperties, deletenodeproperty,
       addnodelabel, addnodelabels, updatenodelabels, deletenodelabel, getnodelabels,
       getnodesforlabel, getlabels, getrel, getrels, getneighbors, createrel, deleterel, getrelproperty,
       getrelproperties, updaterelproperties, cypherQuery
export Connection, Result

const DEFAULT_HOST = "localhost"
const DEFAULT_PORT = 7474
const DEFAULT_URI = "/db/data/"

const JSONObject{T <: AbstractString} = Union{Dict{T,Any},Void}  # UTF8String
const JSONArray = Union{Vector,Void}
const JSONData{T <: AbstractString} = Union{JSONObject,JSONArray,T,Number,Void}

const QueryData = Union{Dict{Any,Any},Void}

# ----------
# Connection
# ----------

"""
   Connection()

### Examples
```julia-repl
julia> con = Neo4j.Connection("localhost")
Neo4j.Connection(false, "localhost", 7474, "/db/data/", "http://localhost:7474/db/data/", "", "")
```
"""
struct Connection
   host::AbstractString #UTF8String
   tls::Bool
   port::Int
   path::AbstractString #UTF8String
   url::AbstractString #UTF8String
   user::AbstractString #UTF8String
   password::AbstractString #UTF8String

   Connection{T <: AbstractString}(host::T; port = DEFAULT_PORT, path = DEFAULT_URI, tls = false, user = "", password = "") =
      new(string(host), tls, port, string(path), string("http://$host:$port$path"), string(user), string(password))
   Connection() = Connection(DEFAULT_HOST)
end

function connurl(c::Connection)
  proto = ifelse(c.tls, "https", "http")
  "$(proto)://$(c.host):$(c.port)$(c.path)"
end

function connurl(c::Connection, suffix::T) where {T <: AbstractString}
  url = connurl(c)
  "$(url)$(suffix)"
end

function connheaders(c::Connection)
  headers = Dict(
    "Accept" => "application/json; charset=UTF-8",
    "Host" => "$(c.host):$(c.port)")
  if c.user != "" && c.password != ""
    payload = "$(c.user):$(c.password)" |> base64encode
    headers["Authorization"] = "Basic $(payload)"
  end
  headers
end

# -----
# Graph
# -----

struct Graph
    # TODO extensions
    node::AbstractString                  #UTF8String
    node_index::AbstractString            #UTF8String
    relationship_index::AbstractString    #UTF8String
    extensions_info::AbstractString       #UTF8String
    relationship_types::AbstractString    #UTF8String
    batch::AbstractString                 #UTF8String
    cypher::AbstractString                #UTF8String
    indexes::AbstractString               #UTF8String
    constraints::AbstractString           #UTF8String
    transaction::AbstractString           #UTF8String
    node_labels::AbstractString           #UTF8String
    version::AbstractString               #UTF8String
    connection::Connection
    relationship::AbstractString          #UTF8String # Not in the spec
end

# UTF8String
Graph{T <: AbstractString}(data::Dict{T,Any}, conn::Connection) = Graph(data["node"], data["node_index"], data["relationship_index"],
    data["extensions_info"], data["relationship_types"], data["batch"], data["cypher"], data["indexes"],
    data["constraints"], data["transaction"], data["node_labels"], data["neo4j_version"], conn,
    "$(conn.url)relationship")

function getgraph(conn::Connection)
    resp = get(conn.url; headers=connheaders(conn))
    if resp.status !== 200
        error("Connection to server unsuccessful: $(resp.status)")
    end
    #Note the Requests lib returns UTF8String, so lets just use that as a standard for REST calls [GearsAD]
    # @show typeof(Requests.json(resp))
    Graph(Dict{AbstractString,Any}(Requests.json(resp)), conn) # UTF8String
end

function getgraph()
    getgraph(Connection())
end

# ----
# Node
# ----

struct Node
    # TODO extensions
    paged_traverse::AbstractString #UTF8String
    labels::AbstractString #UTF8String
    outgoing_relationships::AbstractString #UTF8String
    traverse::AbstractString #UTF8String
    all_typed_relationships::AbstractString #UTF8String
    all_relationships::AbstractString #UTF8String
    property::AbstractString #UTF8String
    self::AbstractString #UTF8String
    outgoing_typed_relationships::AbstractString #UTF8String
    properties::AbstractString #UTF8String
    incoming_relationships::AbstractString #UTF8String
    incoming_typed_relationships::AbstractString #UTF8String
    create_relationship::AbstractString #UTF8String
    data::JSONObject
    metadata::Dict{AbstractString, Any} #UTF8String
    id::Int64
    graph::Graph

    #Constructors
    Node() = new()
    Node(data::JSONObject, graph::Graph) = new(data["paged_traverse"], data["labels"],
         data["outgoing_relationships"], data["traverse"], data["all_typed_relationships"],
         data["all_relationships"], data["property"],
         data["self"], data["outgoing_typed_relationships"], data["properties"],
         data["incoming_relationships"], data["incoming_typed_relationships"],
         data["create_relationship"], data["data"], data["metadata"],
         split(data["self"], "/")[end] |> parse, graph)
end

# ----------
# Statements
# ----------

struct Statement
  statement::AbstractString
  parameters::Dict
end

# -------
# Results
# -------

struct Result
  results::Vector
  errors::Vector
end

# ------------
# Transactions
# ------------

include("transaction.jl")

# --------
# Requests
# --------

function request(url::AbstractString, method::Function, exp_code::Int,
                 headers::Dict{T, T}; json::JSONData = nothing, # ASCIIString,ASCIIString
                 query::QueryData = nothing) where {T <: AbstractString}
    if json == nothing && query == nothing
        resp = method(url; headers = headers)
    elseif json == nothing
        resp = method(url; headers = headers, query = query)
    elseif query == nothing
        resp = method(url; headers = headers, json = json)
    else
        # TODO Figure out if this should ever occur and change it to an error if not
        resp = method(url; headers = headers, json = json, query=query)
    end
    if resp.status !== exp_code
        respdata = Requests.json(resp)
        if respdata !== nothing && "message" in keys(respdata)
            error("Neo4j error: $(respdata["message"])")
        else
            error("Neo4j error: $(url) returned $(resp.status)")
        end
    end
    resp
end

# -----------------
# External requests
# -----------------

function createnode(graph::Graph, props::JSONData = nothing)
    resp = request(graph.node, Requests.post, 201, connheaders(graph.connection); json=props)
    jsrsp = Dict{AbstractString,Any}(Requests.json(resp)) # UTF8String
    # @show typeof(jsrsp)
    Node(jsrsp, graph)
end

function getnode(graph::Graph, id::Int)
    url = "$(graph.node)/$id"
    resp = request(url, Requests.get, 200, connheaders(graph.connection))
    Node(Dict{AbstractString,Any}(Requests.json(resp)), graph) # UTF8String
end

function getnode(node::Node)
    getnode(node.graph, node.id)
end

function deletenode(node::Node)
    request(node.self, Requests.delete, 204, connheaders(node.graph.connection))
end

function deletenode(graph::Graph, id::Int)
    node = getnode(graph, id)
    deletenode(node)
end

function setnodeproperty(node::Node, name::T, value::Any) where {T <: AbstractString}
    url = replace(node.property, "{key}", name)
    request(url, Requests.put, 204, connheaders(node.graph.connection); json=value)
end

function setnodeproperty(graph::Graph, id::Int, name::T, value::Any) where {T <: AbstractString}
    node = getnode(graph, id)
    setnodeproperty(node, name, value)
end

function updatenodeproperties(node::Node, props::JSONObject)
    resp = request(node.properties, Requests.put, 204, connheaders(node.graph.connection); json=props)
end

function getnodeproperty(node::Node, name::T) where {T <: AbstractString}
    url = replace(node.property, "{key}", name)
    resp = request(url, Requests.get, 200, connheaders(node.graph.connection))
    Requests.json(resp)
end

function getnodeproperties(node::Node)
    resp = request(node.properties, Requests.get, 200, connheaders(node.graph.connection))
    Requests.json(resp)
end

function getnodeproperties(graph::Graph, id::Int)
    node = getnode(graph, id)
    getnodeproperties(node)
end

function deletenodeproperties(node::Node)
    request(node.properties, Requests.delete, 204, connheaders(node.graph.connection))
end

function deletenodeproperty(node::Node, name::T) where {T <: AbstractString}
    url = replace(node.property, "{key}", name)
    request(url, Requests.delete, 204, connheaders(node.graph.connection))
end

function addnodelabel(node::Node, label::T) where {T <: AbstractString}
    request(node.labels, Requests.post, 204, connheaders(node.graph.connection); json=label)
end

function addnodelabels(node::Node, labels::JSONArray)
    request(node.labels, Requests.post, 204, connheaders(node.graph.connection); json=labels)
end

function updatenodelabels(node::Node, labels::JSONArray)
    request(node.labels, Requests.put, 204, connheaders(node.graph.connection); json=labels)
end

function deletenodelabel(node::Node, label::T) where {T <: AbstractString}
    url = "$(node.labels)/$label"
    request(url, Requests.delete, 204, connheaders(node.graph.connection))
end

function getnodelabels(node::Node)
    resp = request(node.labels, Requests.get, 200, connheaders(node.graph.connection))
    Requests.json(resp)
end

function getnodesforlabel(graph::Graph, label::T, props::JSONObject=nothing) where {T <: AbstractString}
    # TODO Shouldn't this url be available in the api somewhere?
    url = "$(graph.connection.url)label/$label/nodes"
    resp = request(url, Requests.get, 200, connheaders(graph.connection); query=props)
    [Node(Dict{AbstractString,Any}(nodedata), graph) for nodedata = Requests.json(resp)]
end

function getlabels(graph::Graph)
    resp = request(graph.node_labels, Requests.get, 200, connheaders(graph.connection))
    Requests.json(resp)
end

# -------------
# Relationships
# -------------

struct Relationship
    relstart::AbstractString #UTF8String
    property::AbstractString #UTF8String
    self::AbstractString #UTF8String
    properties::AbstractString #UTF8String
    metadata::Dict{AbstractString ,Any} #UTF8String
    reltype::AbstractString #UTF8String
    relend::AbstractString #UTF8String
    data::JSONObject
    id::Int
    graph::Graph
end

Relationship(data::JSONObject, graph::Graph) = Relationship(data["start"], data["property"],
        data["self"], data["properties"], data["metadata"], data["type"], data["end"], data["data"],
        split(data["self"], "/")[end] |> parse, graph)

function getrels(node::Node; incoming::Bool = true, outgoing::Bool = true)
  rels = Vector{Relationship}()
  if(incoming)
    resp = request(node.incoming_relationships, Requests.get, 200, connheaders(node.graph.connection))
    for rel=Requests.json(resp)
      push!(rels, Relationship(Dict{AbstractString,Any}(rel), node.graph)) #UTF8String
    end
  end
  if(outgoing)
    resp = request(node.outgoing_relationships, Requests.get, 200, connheaders(node.graph.connection))
    for rel=Requests.json(resp)
      push!(rels, Relationship(Dict{AbstractString,Any}(rel), node.graph)) #UTF8String
    end
  end
  rels
end

function getneighbors(node::Node; incoming::Bool = true, outgoing::Bool = true)
  neighbors = Vector{Node}()

  # Do incoming
  if(incoming)
    rels = getrels(node, incoming = true, outgoing = false)
    for rel=rels
      resp = request(rel.relstart, Requests.get, 200, connheaders(node.graph.connection))
      push!(neighbors, Node(Dict{AbstractString,Any}(Requests.json(resp)), node.graph))  # UTF8String
    end
  end
  if(outgoing)
    rels = getrels(node, incoming = false, outgoing = true)
    for rel=rels
      resp = request(rel.relend, Requests.get, 200, connheaders(node.graph.connection))
      push!(neighbors, Node(Dict{AbstractString,Any}(Requests.json(resp)), node.graph))  # UTF8String
    end
  end
  neighbors
end

function getrel(graph::Graph, id::Int)
    url = "$(graph.relationship)/$id"
    resp = request(url, Requests.get, 200, connheaders(graph.connection))
    Relationship(Dict{AbstractString,Any}(Requests.json(resp)), graph)  # UTF8String
end

function createrel(from::Node, to::Node, reltype::AbstractString; props::JSONObject=nothing)
    body = Dict{AbstractString, Any}("to" => to.self, "type" => uppercase(reltype)) # UTF8String
    if props !== nothing
        body["data"] = props
    end
    resp = request(from.create_relationship, Requests.post, 201, connheaders(from.graph.connection), json=body)
    Relationship(Dict{AbstractString,Any}(Requests.json(resp)), from.graph) # UTF8String
end

function deleterel(rel::Relationship)
    request(rel.self, Requests.delete, 204, connheaders(rel.graph.connection))
end

function getrelproperty(rel::Relationship, name::AbstractString)
    url = replace(rel.property, "{key}", name)
    resp = request(url, Requests.get, 200, connheaders(rel.graph.connection))
    Requests.json(resp)
end

function getrelproperties(rel::Relationship)
    resp = request(rel.properties, Requests.get, 200, connheaders(rel.graph.connection))
    Requests.json(resp)
end

function updaterelproperties(rel::Relationship, props::JSONObject)
    request(rel.properties, Requests.put, 204, connheaders(rel.graph.connection); json=props)
end

# ------------
# Cypher query
# ------------
include("cypherQuery.jl")

end # module


text:
"""
    r(t, x, T, v, P)

Returns ``r(t,x) = \\operatorname{grad}_x \\log p(t,x; T, v)`` where
``p`` is the transition density of the process ``P``.
"""
function r(t, x, T, v, P)
    H(t, T, P, V(t, T, v, P) - x)
end

tau(s, t, T) = t + (s.-t).*(2-(s-t)/(T-t))
tau(ss::Vector) = tau(ss, ss[1], ss[end])

# fallback
btilde(t, x, Po) = b(t, x, auxiliary(Po))
btilde!(t, x, out, Po) = b!(t, x, out, auxiliary(Po))
_btilde(t, x, Po) = _b(t, x, auxiliary(Po))
_btilde!(t, x, out, Po) = _b!(t, x, out, auxiliary(Po))

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

"""
    BridgeProp(Target::ContinuousTimeProcess, tt, v, a, cs)

Simple bridge proposal derived from a linear process with time dependent drift given by a [`CSpline`](@ref)
and constant diffusion coefficient `a`.

"""
struct BridgeProp{T} <: ContinuousTimeProcess{T}
    Target
    tt::Vector{Float64}
    v::Tuple{T,T}
    cs::CSpline{T}
    a
    Γ
    BridgeProp{T}(Target::ContinuousTimeProcess{T}, tt, v, a, cs) where T =
        new(Target, tt, v, cs, a, inv(a))
end

startpoint(P::BridgeProp) = P.v[1]
endpoint(_, P::BridgeProp) = P.v[2]

BridgeProp(Target::ContinuousTimeProcess{T}, tt, v, a, cs=CSpline(first(tt), last(tt), zero(T))) where {T} = BridgeProp{T}(Target, tt, v, a, cs)

h(t,x, P::BridgeProp) = P.v[2] - x -  integrate(P.cs, t,  last(P.tt))
b(t, x, P::BridgeProp) = b(t, x, P.Target) + a(t, x, P.Target)*r(t, x, P)

function bderiv(t, x, P::BridgeProp)
    @assert(constdiff(P))
    bderiv(t, x, P.Target) - a(t, x, P.Target)*P.Γ/(last(P.tt) - t)
end

σ(t, x, P::BridgeProp) = σ(t, x, P.Target)
a(t, x, P::BridgeProp) = a(t, x, P.Target)
Γ(t, x, P::BridgeProp) = Γ(t, x, P.Target)
constdiff(P::BridgeProp) = constdiff(P.Target)

btilde(t, x, P::BridgeProp) = P.cs(t)
atilde(t, x, P::BridgeProp) = P.a
ptilde(P::BridgeProp) = Ptilde(P.cs, σ(last(P.tt), P.v[2], P.Target))

function r(t, x, P::BridgeProp)
    P.Γ*h(t, x, P)/(last(P.tt) - t)
end
function H(t, x, P::BridgeProp)
    P.Γ/(last(P.tt) - t)
end

function lptilde(P::BridgeProp{T}) where T
    logpdfnormal(P.v[2] - (P.v[1] + integrate(P.cs, first(P.tt), last(P.tt))), (last(P.tt) - first(P.tt))*P.a)
end

"""
    GuidedProp

General bridge proposal process, only assuming that `Pt` defines `H` and `r` in the right way.
"""
struct GuidedProp{T} <: ContinuousTimeProcess{T}
    Target
    t0; v0; t1; v1
    Pt::ContinuousTimeProcess{T}
    GuidedProp{T}(Target::ContinuousTimeProcess{T}, t0, v0, t1, v1, Pt) where T =
        new(Target, t0, v0, t1, v1, Pt)
end
startpoint(P::GuidedProp) = P.v0
endpoint(_, P::GuidedProp) = P.v1

GuidedProp(Target::ContinuousTimeProcess{T}, t0, v0, t1, v1, Pt::ContinuousTimeProcess{T}) where {T} = GuidedProp{T}(Target, t0, v0, t1, v1, Pt)

V(t, P::GuidedProp) = V(t, P.t1, P.v1, P.Pt)
dotV(t, P::GuidedProp) = dotV(t, P.t1, P.v1, P.Pt)
r(t, x, P::GuidedProp) = r(t, x, P.t1, P.v1, P.Pt)
b(t, x, P::GuidedProp) = b(t, x, P.Target) + a(t, x, P.Target)*r(t, x, P.t1, P.v1, P.Pt)
σ(t, x, P::GuidedProp) = σ(t, x, P.Target)
a(t, x, P::GuidedProp) = a(t, x, P.Target)
Γ(t, x, P::GuidedProp) = Γ(t, x, P.Target)
constdiff(P::GuidedProp) = constdiff(P.Target) && constdiff(P.Pt)

btilde(t, x, P::GuidedProp) = b(t,x,P.Pt)
atilde(t, x, P::GuidedProp) = a(t,x,P.Pt)
ptilde(P::GuidedProp) = P.Pt

function lptilde(P::GuidedProp)
     lp(P.t0, P.v0, P.t1, P.v1, P.Pt)
end

#################################################
"""
    LinearNoiseAppr(tt, P, x, a, direction = forward)

Precursor of the linear noise approximation of `P`. For now no attempt is taken
to add in a linearization around the deterministic path. `direction` can be one of
`:forward`, `:backward` or `:nothing`. The latter corresponds to choosing `β == 0`.
"""
struct LinearNoiseAppr{R,S,T} <: ContinuousTimeProcess{T}
    Target::R
    Y::SamplePath{T}
    a::S
    function LinearNoiseAppr(tt_, P::R, x::T, a::S, direction = forward) where {R,S,T}
        tt = collect(tt_)
        N = length(tt)
        Y = SamplePath(tt, zeros(T, N))
        if direction == :forward
            solve!(R3(), b, Y, x, P)
        elseif direction == :backward
            solvebackward!(R3(), b, Y, x, P)
        elseif direction != :nothing
            throw(ArgumentError("Wrong `direction`"))
        end
        new{R,S,T}(P, Y, a)
    end
end

_b((i,t)::IndexedTime, x, P::LinearNoiseAppr) = β_((max(i,2),t), P)
B(t, P::LinearNoiseAppr) = 0I
β((i,t)::IndexedTime, P::LinearNoiseAppr) = (P.Y.yy[i]-P.Y.yy[i-1])/(P.Y.tt[i]-P.Y.tt[i-1])
a(t, x, P::LinearNoiseAppr) = P.a
constdiff(::LinearNoiseAppr) = true

"""
    GuidedBridge

Guided proposal process for diffusion bridge using backward recursion.

    GuidedBridge(tt, P, Pt, v)

Constructor of guided proposal process for diffusion bridge of `P` to `v` on
the time grid `tt` using guiding term derived from linear process `Pt`.

    GuidedBridge(tt, P, Pt, V, H♢)

Guided proposal process for diffusion bridge of `P` to `v` on
the time grid `tt` using guiding term derived from linear process `Pt`.
Initialize using [`Bridge.gpupdate`](@ref)(H♢, V, L, Σ, v)
"""
struct GuidedBridge{T,S,R2,R} <: ContinuousTimeProcess{T}
    Target::R
    Pt::R2
    tt::Vector{Float64}
    H♢::Vector{S}
    V::Vector{T}

    function GuidedBridge(tt_, P::R, Pt::R2, v::T, h♢::S = Bridge.outer(zero(v))) where {T,R,R2,S}
        tt = collect(tt_)
        N = length(tt)
        H♢ = SamplePath(tt, zeros(S, N))
        V = SamplePath(tt, zeros(T, N))
        gpHinv!(H♢, Pt, h♢)
        gpV!(V, Pt, v)
        new{T,S,R2,R}(P, Pt, tt, H♢.yy, V.yy)
    end
    function GuidedBridge(tt_, P::R, Pt::Union{LinearNoiseAppr, LinearAppr}, v::T, h♢::S = Bridge.outer(zero(v))) where {T,R,S}
        tt = collect(tt_)
        N = length(tt)
        H♢ = SamplePath(tt, zeros(S, N))
        V = SamplePath(tt, zeros(T, N))
        solvebackwardi!(Heun(), ((i,t), K, P) ->  B((i,t), P)*K + K*B((i,t), P)' - a((i,t), P), H♢, h♢, Pt)
        solvebackwardi!(Heun(), b, V, v, Pt)
        new{T,S,typeof(Pt),R}(P, Pt, tt, H♢.yy, V.yy)
    end
end

_b((i,t)::IndexedTime, x, P::GuidedBridge) = b(P.tt[i], x, P.Target) + a(P.tt[i], x, P.Target)*(P.H♢[i]\(P.V[i] - x))
r((i,t)::IndexedTime, x, P::GuidedBridge) = P.H♢[i]\(P.V[i] - x)
H((i,t)::IndexedTime, x, P::GuidedBridge) = inv(P.H♢[i])

σ(t, x, P::GuidedBridge) = σ(t, x, P.Target)
a(t, x, P::GuidedBridge) = a(t, x, P.Target)
Γ(t, x, P::GuidedBridge) = Γ(t, x, P.Target)
constdiff(P::GuidedBridge) = constdiff(P.Target) && constdiff(P.Pt)
btilde(it, x, P::GuidedBridge) = b(it, x, P.Pt)
atilde(it, x, P::GuidedBridge) = a(it, x, P.Pt)

@inline _traceB(t, x, P) = tr(Bridge.B(t, P))
traceB(tt, P) = solve(Bridge.R3(), _traceB, tt, 0.0, P)

lptilde(P::GuidedBridge, u) = logpdfnormal(P.V[1] - u, P.H♢[1]) - traceB(P.tt, P.Pt)

# fallback for testing
lptilde2(P::GuidedBridge, u) = logpdfnormal(P.V[end] - gpmu(P.tt, u, P.Pt), gpK(P.tt, Bridge.outer(zero(u)), P.Pt))

# H♢_ = H♢ -  H♢*L'*inv(Σ + L*H♢*L')*L*H♢
# V_ = H♢_ * (L'*inv(Σ)*v  + H♢*V)

"""
    gpupdate(H♢, V, L, Σ, v)
    gpupdate(P, V, L, Σ, v)

Return updated `H♢, V` when observation `v` at time zero with error `Σ` is observed.
"""
function gpupdate(H♢, V, L, Σ, v)
    if all(diag(H♢) .== Inf)
        H♢_ = SMatrix(inv(L' * inv(Σ) * L))
        V_ = (L' * inv(Σ) * L)\(L' * inv(Σ) *  v)
        H♢_, V_
    else
        Z = I - H♢*L'*inv(Σ*I + L*H♢*L')*L
        Z*H♢, Z*H♢*L'*inv(Σ)*v + Z*V
    end
end

function gpupdate(H♢, V::SVector, L, Σ::Float64, v)
    if all(diag(H♢) .== Inf)
        H♢_ = SMatrix(inv(L' * inv(Σ) * L))
        V_ = (L' * inv(Σ) * L)\(L' * inv(Σ) *  v)
        H♢_, V_
    else
        Z = I - H♢*L'*inv(Σ + L*H♢*L')*L
        Z*H♢, SVector(Z*H♢*L'*inv(Σ)*v) + Z*V
    end
end

gpupdate(P::GuidedBridge, L, Σ, v) = gpupdate(P.H♢[1], P.V[1], L, Σ, v)

# Alternatives to try
#   KT =  PhiTS*(KS + GP2.H♢[1])*PhiTS'
#   KT = PhiTS*KS*PhiTS' + KTS
#   logdet( L*KS*L' + Σ - L*KS*inv(KS + GP2.H♢[1])*KS*L' ) + logdet(KS + GP2.H♢[1]) + 2logdet(PhiTS)
#   logdet( KS + GP2.H♢[1] - KS*L'*inv(L*KS*L' + Σ)*L*KS  ) + logdet(L*KS*L' + Σ ) + 2logdet(PhiTS)
function logdetU(GP1, GP2, L, Σ)
    PhiS = fundamental_matrix(GP1.tt, GP1.Pt)
    PhiTS = fundamental_matrix(GP2.tt, GP2.Pt)
    K =  PhiS*GP1.H♢[1]*PhiS' - GP1.H♢[end]
    H = GP2.H♢[1]
    logdet(inv(K) + L'*inv(Σ)*L + inv(H)) + logdet(Σ) + logdet(H) + logdet(K) + 2logdet(PhiTS)
end

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

struct PBridgeProp{T} <: ContinuousTimeProcess{T}
    Target
    t0; v0; tm; vm; t1; v1
    L; Lt; Σ
    cs::CSpline{T}
    a
    Γ
    PBridgeProp{T}(Target::ContinuousTimeProcess{T}, t0, v0, tm, vm, t1, v1,  L, Σ, a, cs) where T =
        new(Target, t0, v0, tm, vm, t1, v1, L, L', Σ, cs, a, inv(a))
end
PBridgeProp(Target::ContinuousTimeProcess{T}, t0, v0, tm, vm, t1, v1, L, Σ, a, cs=CSpline(t0, t1, zero(T))) where {T} = PBridgeProp{T}(Target, t0, v0, tm, vm, t1, v1,  L, Σ, a, cs)

h1(t, x, P::PBridgeProp) = P.vm - x - integrate(P.cs, t,  P.tm)
h2(t, x, P::PBridgeProp) = P.v1 - x - integrate(P.cs, t,  P.t1)

N(t, P::PBridgeProp) = inv(P.L*P.a*P.Lt*(P.tm - t) + (P.t1 - t)/(P.t1 - P.tm)*P.Σ)
Q(t, P::PBridgeProp) = P.Lt*N(t, P)*P.L

function b(t, x, P::PBridgeProp)
    if t >= P.tm
        b(t, x, P.Target) + a(t, x, P.Target)*P.Γ*h2(t, x, P)/(P.t1 -t)
    else
        b(t, x, P.Target) + a(t, x, P.Target)*(Q(t, P)*h1(t, x, P)  +    (P.Γ - Q(t, P)*(P.tm -t))*h2(t, x, P)/(P.t1 -t))
    end
end
function r(t, x, P::PBridgeProp)
    if t >= P.tm
        P.Γ*h2(t, x, P)/(P.t1 -t)
    else
        (Q(t, P)*h1(t, x, P)  +    (P.Γ - Q(t, P)*(P.tm -t))*h2(t, x, P)/(P.t1 -t))
    end
end
function H(t, x, P::PBridgeProp)
    if t >= P.tm
        P.Γ/(P.t1 - t)
    else
        P.Γ/(P.t1 - t) + Q(t, P)*(P.t1-P.tm)/(P.t1 - t)
    end
end

function lptilde(P::PBridgeProp)
	n = N(P.t0, P)*(P.tm-P.t0)
	U = Any[	(P.t1-P.t0)/(P.t1-P.tm)/(P.tm-P.t0)*n 		-n*P.L/(P.t1-P.tm)
			-P.L'*n/(P.t1-P.tm) 				(P.Γ + P.L'*n*P.L*(P.tm-P.t0)/(P.t1-P.tm))/(P.t1-P.t0)]
	ldm = sumlogdiag(cholupper(U[1,1])') +sumlogdiag(cholupper(U[2,2] - (U[2,1]*inv(U[1,1])*U[1,2]))')

	mu = [P.L*h1(P.t0, P.v0, P); h2(P.t0, P.v0, P)]
	-length(mu)/2*log(2pi) + ldm - 0.5*dot(mu,U*mu)
end

btilde(t, x, P::PBridgeProp) = P.cs(t)
atilde(t, x, P::PBridgeProp) = P.a
σ(t, x, P::PBridgeProp) = σ(t, x, P.Target)
a(t, x, P::PBridgeProp) = a(t, x, P.Target)
constdiff(P::PBridgeProp) = constdiff(P.Target)

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

struct FilterProp{T} <: ContinuousTimeProcess{T}
    Target
    t0; v0; t1; v1
    L; Lt; Σ
    cs
    a
    Γ
    FilterProp{T}(Target::ContinuousTimeProcess{T}, t0, v0, t1, v1,  L, Σ, a, cs) where T =
        new(Target, t0, v0, t1, v1, L, L', Σ, cs, a, inv(a))
end
FilterProp(Target::ContinuousTimeProcess{T}, t0, v0, t1, v1,  L, Σ, a, cs=CSpline(t0, t1, zero(T))) where {T} = FilterProp{T}(Target, t0, v0, t1, v1,  L, Σ, a, cs)

h(t,x, P::FilterProp) = P.v1 - x - integrate(P.cs, t, P.t1)

H(t, P::FilterProp) = P.Lt*inv(P.L*P.a*P.Lt*(P.t1 -t) + P.Σ)*P.L

function b(t, x, P::FilterProp)
    b(t, x, P.Target) + a(t, x, P.Target)*r(t, x, P)
end
function r(t, x, P::FilterProp)
       H(t, P)*h(t, x, P)
end

btilde(t, x, P::FilterProp) = P.cs(t)
atilde(t, x, P::FilterProp) = P.a
σ(t, x, P::FilterProp) = σ(t, x, P.Target)
a(t, x, P::FilterProp) = a(t, x, P.Target)
constdiff(P::FilterProp) = constdiff(P.Target)

function bderiv(t, x, P::FilterProp)
    @assert(constdiff(P))
    bderiv(t, x, P.Target) - a(t, x, P.Target)*H(t, P)
end

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

struct DHBridgeProp{T} <: ContinuousTimeProcess{T}
    Target
    t0; v0::T; t1; v1::T
    DHBridgeProp{T}(Target::ContinuousTimeProcess{T}, t0, v0, t1, v1) where T =
        new(Target, t0, v0, t1, v1)

end
DHBridgeProp(Target::ContinuousTimeProcess{T}, t0, v0, t1, v1) where {T} = DHBridgeProp{T}(Target, t0, v0, t1, v1)

b(t, x, P::DHBridgeProp) = (P.v1 - x )/(P.t1 - t)
σ(t, x, P::DHBridgeProp) = σ(t, x, P.Target)
a(t, x, P::DHBridgeProp) = a(t, x, P.Target)
Γ(t, x, P::DHBridgeProp) = Γ(t, x, P.Target)
constdiff(P::DHBridgeProp) = constdiff(P.Target)

function llikelihood(Xcirc::SamplePath{T}, P::DHBridgeProp{T}) where T
    tt = Xcirc.tt
    xx = Xcirc.yy

    som::Float64 = 0.
    N = length(tt)
    for i in 1:N-1 #skip last value, summing over n-1 elements
        s = tt[i]
        sh = tt[i+1]
        x = xx[i]
        xh = xx[i+1]
        m = b(s,x,P.Target)
        G = Γ(s,x,P.Target)
        Gh = Γ(sh,xh,P.Target)

        som += dot(m,G*(xh - x - 0.5*m*(sh-s))) # girsanov(Xcirc, P.Target, Wiener{Float64}())
        if i < N-1
            y = xh - P.v1
            som -= 0.5*dot(y, (Gh-G)*y)/(P.t1 - tt[i+1])
        end
    end
    som
end

function lptilde(P::DHBridgeProp{T}) where T
    dv = P.v1-P.v0
    -length(P.v1)/2*log(2pi*(P.t1-P.t0)) -0.5*logdet(a(P.t1,P.v1,P)) - (0.5/(P.t1-P.t0))*dot(dv, Γ(P.t0,P.v0,P)*dv)
end

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

# using left approximation #::Union{GuidedProp{T},BridgeProp{T},PBridgeProp{T},FilterProp{T},BesselBridge{T}}
function llikelihoodleft(Xcirc::SamplePath{T}, Po) where T
    tt = Xcirc.tt
    xx = Xcirc.yy

    som::Float64 = 0.
    for i in 1:length(tt)-1 #skip last value, summing over n-1 elements
        s = tt[i]
        x = xx[i]
        r = Bridge.r(s, x, Po)
        som += (dot(b(s,x, target(Po)) - btilde(s,x, Po), r)  ) * (tt[i+1]-tt[i])
        if !constdiff(Po)
            som -= 0.5*tr((a(s,x, target(Po)) - atilde(s, x, Po))*(H(s,x,Po) -  r*r')) * (tt[i+1]-tt[i])
        end
    end
    som
end

# modern implementation
# using left approximation
function llikelihood(::LeftRule, Xcirc::SamplePath, Po::GuidedBridge; skip = 0)
    tt = Xcirc.tt
    xx = Xcirc.yy

    som::Float64 = 0.
    for i in 1:length(tt)-1-skip #skip last value, summing over n-1 elements
        s = tt[i]
        x = xx[i]
        r = Bridge.r((i,s), x, Po)
        som += dot( _b((i,s), x, target(Po)) - _b((i,s), x, auxiliary(Po)), r ) * (tt[i+1]-tt[i])
        if !constdiff(Po)
            H = Hi(i, x, Po)
            som -= 0.5*tr( (a((i,s), x, target(Po)) - atilde((i,s), x, Po))*(H) ) * (tt[i+1]-tt[i])
            som += 0.5*( r'*(a((i,s), x, target(Po)) - atilde((i,s), x, Po))*r ) * (tt[i+1]-tt[i])
        end
    end
    som
end

#using trapezoidal rule
function llikelihoodtrapez(Xcirc::SamplePath{T}, Po::Union{GuidedProp{T},BridgeProp{T},PBridgeProp{T},FilterProp{T}}) where T
    tt = Xcirc.tt
    xx = Xcirc.yy

    som::Float64 = 0.
    i = 1
    s = tt[i]
    x = xx[i]
    r = Bridge.r(s, x, Po)
    som += 0.5*(dot(b(s,x, target(Po)) - btilde(s,x, Po), r)  ) * (tt[i+1]-tt[i])
    if !constdiff(Po)
        som -= 0.25*tr( (a(s,x, target(Po)) - atilde(s, x, Po))*(H(s,x,Po) -  r*r') ) * (tt[i+1]-tt[i])
    end

    for i in 2:length(tt)-1 #skip last value, summing over n-1 elements
        s = tt[i]
        x = xx[i]
        r = Bridge.r(s, x, Po)
        som += 0.5*( dot(b(s,x, target(Po)) - btilde(s,x, Po), r) ) * (tt[i+1]-tt[i-1])
        if !constdiff(Po)
            som -= 0.25*tr( (a(s,x, target(Po)) - atilde(s, x, Po))*(H(s,x,Po) -  r*r') ) * (tt[i+1]-tt[i-1])
        end
    end
    som
end
llikelihood(Xcirc::SamplePath{T}, Po::Union{GuidedProp{T},BridgeProp{T},PBridgeProp{T},FilterProp{T}}) where {T} = llikelihoodleft(Xcirc, Po)
#llikelihood{T}(Xcirc::SamplePath{T}, Po::Union{GuidedProp{T},BridgeProp{T},PBridgeProp{T},FilterProp{T}}) = llikelihoodtrapez(Xcirc, Po)

ext: jl


text:
###
### Finders
###
###
### Finding positions meeting predicates in biological sequence types.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

function Base.findnext(val, seq::BioSequence, start::Integer)
    checkbounds(seq, start)
    v = convert(eltype(seq), val)
    for i in start:lastindex(seq)
        if inbounds_getindex(seq, i) == v
            return i
        end
    end
    return nothing
end

function Base.findnext(f::Function, seq::BioSequence, start::Integer)
    checkbounds(seq, start)
    for i in start:lastindex(seq)
        if f(inbounds_getindex(seq, i))
            return i
        end
    end
    return nothing
end

# No ambiguous sites can exist in a nucleic acid sequence using the two-bit alphabet.
Base.findnext(::typeof(isambiguous), seq::BioSequence{<:NucleicAcidAlphabet{2}}, from::Integer) = nothing

function Base.findprev(val, seq::BioSequence, start::Integer)
    checkbounds(seq, start)
    v = convert(eltype(seq), val)
    for i in start:-1:1
        if inbounds_getindex(seq, i) == v
            return i
        end
    end
    return nothing
end

function Base.findprev(f::Function, seq::BioSequence, start::Integer)
    checkbounds(seq, start)
    for i in start:-1:firstindex(seq)
        if f(inbounds_getindex(seq, i))
            return i
        end
    end
    return nothing
end

Base.findfirst(val, seq::BioSequence) = findnext(val, seq, 1)
Base.findlast(val, seq::BioSequence) = findprev(val, seq, lastindex(seq))


text:
@testset "include_static_included2" begin
    @test true
    push!(Hijack.RUN, 3)
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
"""
Plot flight data from DDR2 on a map
"""

module FlightPlot

include("utility.jl")
include("Projection.jl")
include("Plot.jl")
include("Map.jl")
include("Airspace.jl")
include("Flights.jl")
include("Network.jl")

end # module


text:
module NMRCalibrate # https://github.com/AI4DBiological-Systems/NMRCalibrate.jl

#import JLD

using LinearAlgebra
import NLopt
import BSON, Statistics, JSON, Interpolations

import MultistartOptimization

#import NMRData
import NMRDataSetup # https://github.com/AI4DBiological-Systems/NMRDataSetup.jl
import BoundedLeastSquares # https://github.com/nboyd/BoundedLeastSquares.jl

# dependencies to NMRSpectraSimulator. Need to add explicitly since these packages are not on the Julia public registry.
#import GISSMOReader # https://github.com/AI4DBiological-Systems/GISSMOReader.jl

# dependencies to NMRSpecifyRegions. Need to add explicitly since these packages are not on the Julia public registry.
import NMRSpectraSimulator # https://github.com/AI4DBiological-Systems/NMRSpectraSimulator.jl

#
import NMRSpecifyRegions # https://github.com/AI4DBiological-Systems/NMRSpecifyRegions

import MonotoneMaps # https://github.com/RoyCCWang/MonotoneMaps.jl

include("../src/types.jl")
include("../src/utils.jl")

include("../src/IO/synthetic_setup.jl")

#include("../src/DSP/dsp.jl")

#include("../src/cost/LS.jl")
include("../src/cost/LS_complex.jl")
#include("../src/cost/LS_z.jl")

include("../src/cost/beta_kappa.jl")
include("../src/cost/beta_w.jl")
include("../src/cost/shift.jl")

include("../src/cost/updates.jl")
include("../src/cost/multiplier_updates.jl")

include("../src/cost/prep.jl")
include("../src/front_end.jl")

include("../src/cost/nested.jl")

include("../src/warp/warp_helpers.jl")

include("../src/align/align_calibration_frontend.jl")
include("../src/align/align_quantification_frontend.jl")

export getwraparoundDFTfreqs

end


text:
module Vehicle

# 1. Export/Imports
export go!

# 2. Interface documentation
# A vehicle (v) must implement the following functions:
#
# power_on!(v) - turn on the vehicle's engine
# power_off!(v) - turn off the vehicle's engine
# turn!(v, direction) - steer the vehicle to the specified direction
# move!(v, distance) - move the vehicle by the specified distance
# position(v) - returns the (x,y) position of the vehicle

# 3. Generig definitions for the interface
"""
Power on the vehicle so it is ready to go.
"""
function power_on! end
function power_off! end
function turn! end
function move! end
function position end

# 4. Logic

# Returns a travel plan from current position to destination
function travel_path(position, destination)
    return round(π/6, digits=2), 1_000 # just a test
end

# travel Logic
function go!(vehicle, destination)
    power_on!(vehicle)
    direction, distance = travel_path(position(vehicle), destination)
    turn!(vehicle, direction)
    move!(vehicle, distance)
    power_off!(vehicle)
    nothing
end

end # module


text:
@testset "randomw" begin
for n=2:6
	println("n=$n")
	A0 = randomw(n,1)
	@test isa(A0,Vector)
	@test abs(sum(A0)-1) < 1e-8
	@test all(diff(A0).<=0)

	nSamp = 100
	A0 = randomw(n,nSamp)
	@test size(A0,2)==nSamp
	@test all(abs.(sum(A0,1)-1) .< 1e-8)
	@test all(diff(A0,1).<=0)
end
end


text:
using Test, Sandbox

all_executors = Sandbox.all_executors

# Can we run `sudo` without a password?  If not, don't attempt to test the privileged runner
if !success(`sudo -k -n true`)
    all_executors = filter(exe -> !isa(exe, PrivilegedUserNamespacesExecutor), all_executors)
end

for executor in all_executors
    if !executor_available(executor)
        @error("Skipping $(executor) tests, as it does not seem to be available")
        continue
    end

    rootfs_dir = artifact"AlpineRootfs"
    @testset "$(executor) Sandboxing" begin
        @testset "capturing stdout/stderr" begin
            stdout = IOBuffer()
            stderr = IOBuffer()
            config = SandboxConfig(
                Dict("/" => rootfs_dir);
                stdout,
                stderr,
            )
            @test run(executor(), config, `/bin/sh -c "echo stdout; echo stderr >&2"`)
            @test String(take!(stdout)) == "stdout\n";
            @test String(take!(stderr)) == "stderr\n";
        end

        @testset "ignorestatus()" begin
            config = SandboxConfig(Dict("/" => rootfs_dir))
            @test_throws ProcessFailedException run(executor(), config, `/bin/sh -c "false"`)
            @test !run(executor(), config, ignorestatus(`/bin/sh -c "false"`))
        end

        @testset "environment passing" begin
            # Ensure all those pesky "special" variables make it through
            env = Dict(
                "PATH" => "for",
                "LD_LIBRARY_PATH" => "science",
                "DYLD_LIBRARY_PATH" => "you",
                "SHELL" => "monster",
            )
            stdout = IOBuffer()
            config = SandboxConfig(
                Dict("/" => rootfs_dir),
                Dict{String,String}(),
                env;
                stdout,
            )
            user_cmd = `/bin/sh -c "echo \$PATH \$LD_LIBRARY_PATH \$DYLD_LIBRARY_PATH \$SHELL"`
            @test run(executor(), config, user_cmd)
            @test String(take!(stdout)) == "for science you monster\n";

            # Test that setting some environment onto `user_cmd` can override the `config` env:
            user_cmd = setenv(user_cmd, "DYLD_LIBRARY_PATH" => "my", "SHELL" => "friend")
            @test run(executor(), config, user_cmd)
            @test String(take!(stdout)) == "for science my friend\n";
        end

        @testset "reading from maps" begin
            mktempdir() do dir
                open(joinpath(dir, "note.txt"), write=true) do io
                    write(io, "great success")
                end
                stdout = IOBuffer()
                config = SandboxConfig(
                    Dict("/" => rootfs_dir, "/glados" => dir);
                    stdout,
                )
                @test run(executor(), config, `/bin/sh -c "cat /glados/note.txt"`)
                @test String(take!(stdout)) == "great success";
            end
        end

        @testset "writing to workspaces" begin
            mktempdir() do dir
                stdout = IOBuffer()
                config = SandboxConfig(
                    Dict("/" => rootfs_dir),
                    Dict("/glados" => dir);
                )
                @test run(executor(), config, `/bin/sh -c "echo aperture > /glados/science.txt"`)
                @test isfile(joinpath(dir, "science.txt"))
                @test String(read(joinpath(dir, "science.txt"))) == "aperture\n"
            end
        end

        @testset "pipelining" begin
            pipe = PipeBuffer()
            stdout = IOBuffer()
            first_config = SandboxConfig(
                Dict("/" => rootfs_dir),
                stdout = pipe,
            )
            second_config = SandboxConfig(
                Dict("/" => rootfs_dir),
                stdin = pipe,
                stdout = stdout,
            )
            @test run(executor(), first_config, `/bin/sh -c "echo 'ignore me'; echo 'pick this up foo'; echo 'ignore me as well'"`)
            @test run(executor(), second_config, `/bin/sh -c "grep foo"`)
            @test String(take!(stdout)) == "pick this up foo\n";
        end
    end
end


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module SemanticFlowGraphs
using Reexport

include("Doctrine.jl")
include("ontology/Ontology.jl")
include("RawFlowGraphs.jl")
include("Serialization.jl")
include("Visualization.jl")
include("SemanticEnrichment.jl")
include("CLI.jl")

@reexport using .Doctrine
@reexport using .Ontology
@reexport using .RawFlowGraphs
@reexport using .Serialization
@reexport using .Visualization
@reexport using .SemanticEnrichment

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
using MultipleViewGeometry, Base.Test
using BenchmarkTools
using StaticArrays

ℳ = [Point2DH(x,y,1)  for x=-1000:0.5:1000 for y=-1000:0.5:1000]
#ℳ = [Point2DH(x,y,1)  for x=-10:0.5:10 for y=-10:0.5:10]

transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳ)
ℳʹ, 𝐓 = transform(HomogeneousCoordinates(),CanonicalToHartley(),ℳ)

# @time Λ =  [MMatrix{3,3}(diagm([1.0,1.0,0.0])) for i = 1:length(ℳ)]
# @time transform(CovarianceMatrices(), CanonicalToHartley(), Λ , tuple(𝐓))
# @time Z = transform(CovarianceMatrices(), CanonicalToHartley(), Λ , tuple(𝐓))

@time Λ =  [SMatrix{3,3}(diagm([1.0,1.0,0.0])) for i = 1:length(ℳ)]
@time transform(CovarianceMatrices(), CanonicalToHartley(), Λ , 𝐓)
@time Z = transform(CovarianceMatrices(), CanonicalToHartley(), Λ , 𝐓)

@show Λ[1]
@show Z[1]


text: Dict{String,Any}("1" => Dict{String,Any}("value" => "1","type" => "integer"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
using ArgParse,KUparser
include("conll17.jl")
macro tm(_x) :(if o[:fast]; $(esc(_x)); else; info("$(now()) "*$(string(_x))); $(esc(_x)); end) end
macro msg(x); :(if !o[:fast]; info($x); end); end

# Some default data
grcdev = "/mnt/ai/data/nlp/conll17/ud-treebanks-conll2017/UD_Ancient_Greek/grc-ud-dev.conllu"
grcudp = "/mnt/ai/data/nlp/conll17/UDPipe/udpipe-ud-2.0-conll17-170315/models/ancient_greek-ud-2.0-conll17-170315.udpipe"

function main(args="")
    isa(args, AbstractString) && (args=split(args))
    s = ArgParseSettings()
    s.description="Oracle parser takes a conllu file and parses with the best possible moves."
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        # Parsing options
        ("--parse"; default=grcdev; help="file in conllu format to be parsed")
        ("--ptype"; default="ArcEagerR1"; help="Parser type")
        ("--fast"; action=:store_true; help="fewer messages")
    end

    # Process args
    o = parse_args(args, s; as_symbols=true)
    @msg(s.description)
    @msg(string("opts=",[(k,v) for (k,v) in o]...))
    ptype = eval(parse(o[:ptype]))

    # Parse
    if o[:parse] != nothing
        @tm global corpus = readconllu(o[:parse])
        @tm parses = oparse(ptype, corpus)
        @tm writeconllu(corpus, parses, o[:parse])
    end
end

if basename(PROGRAM_FILE)==basename(@__FILE__); main(ARGS); end


text:
@testset "@map macro(ex)" begin
    n = 3
    A = rand(2, 2);
    b = rand(2)
    identitymap = @map x -> I(n)*x
    @test identitymap == MathematicalSystems.IdentityMap(n)
    linearmap = @map x -> A*x
    @test linearmap == MathematicalSystems.LinearMap(A)
    affinemap = @map x -> A*x + b
    @test affinemap == MathematicalSystems.AffineMap(A, b)
end

@testset "@map macro(ex, args)" begin
    n = 3
    identitymap1 = @map(x->x, dim=n)
    identitymap2 = @map(x->x, dim:3)
    @test identitymap1 == identitymap2 == IdentityMap(n)
    X = BallInf(zeros(3), 1.0)
    identitymap3 = @map(x->x, dim:3, x ∈ X)
    @test identitymap3 == ConstrainedIdentityMap(3, X)
end


text:
# ---
# title: 638. Shopping Offers
# id: problem638
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Depth-first Search
# link: <https://leetcode.com/problems/shopping-offers/description/>
# hidden: true
# ---
#
# In LeetCode Store, there are some kinds of items to sell. Each item has a
# price.
#
# However, there are some special offers, and a special offer consists of one or
# more different kinds of items with a sale price.
#
# You are given the each item's price, a set of special offers, and the number
# we need to buy for each item. The job is to output the lowest price you have
# to pay for **exactly** certain items as given, where you could make optimal
# use of the special offers.
#
# Each special offer is represented in the form of an array, the last number
# represents the price you need to pay for this special offer, other numbers
# represents how many specific items you could get if you buy this offer.
#
# You could use any of special offers as many times as you want.
#
# **Example 1:**
#
#
#
#     Input: [2,5], [[3,0,5],[1,2,10]], [3,2]
#     Output: 14
#     Explanation:
#     There are two kinds of items, A and B. Their prices are /$/2 and /$/5 respectively.
#     In special offer 1, you can pay /$/5 for 3A and 0B
#     In special offer 2, you can pay /$/10 for 1A and 2B.
#     You need to buy 3A and 2B, so you may pay /$/10 for 1A and 2B (special offer #2), and /$/4 for 2A.
#
#
# **Example 2:**
#
#
#
#     Input: [2,3,4], [[1,1,0,4],[2,2,1,9]], [1,2,1]
#     Output: 11
#     Explanation:
#     The price of A is /$/2, and /$/3 for B, /$/4 for C.
#     You may pay /$/4 for 1A and 1B, and /$/9 for 2A ,2B and 1C.
#     You need to buy 1A ,2B and 1C, so you may pay /$/4 for 1A and 1B (special offer #1), and /$/3 for 1B, /$/4 for 1C.
#     You cannot add more items, though only /$/9 for 2A ,2B and 1C.
#
#
# **Note:**
#
#   1. There are at most 6 kinds of items, 100 special offers.
#   2. For each item, you need to buy at most 6 of them.
#   3. You are **not** allowed to buy more items than you want, even if that would lower the overall price.
#
#
## @lc code=start
using LeetCode

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


text:
# fbp2_sino_filter.jl

export fbp2_sino_filter

using FFTW

"""
    sino, Hk, hh, nn = fbp2_sino_filter(how, sino;
    ds=1, dsd=0, extra=0, npad=0, decon1=1, window=:none)

Apply ramp-like filters to sinogram(s) for 2D FBP image reconstruction.
Both parallel-beam and fan-beam tomographic geometries are supported.
This approach of sampling the band-limited ramp avoids the aliasing that
would be caused by sampling the ramp directly in the frequency domain.

in
- `how::Symbol`
  `:arc` (3rd generation CT) or `:flat` (for parallel too)
- `sino::AbstractArray{<:Number}` : `[nb (L)]` sinograms

options
- `ds::RealU` : sample spacing (in distance units, e.g., cm) (default: 1)
- `dsd::RealU` : source-to-detector distance, for `:arc` case only. (default: Inf)
- `extra::Int` : # of extra sinogram radial samples to keep (default: 0)
- `npad::Int` : # of padded samples. (default: 0, means next power of 2)
- `decon1::Int` : deconvolve effect of linear interpolator? (default: 1)
- `window::Symbol` : default: `:none`

out
- `sino::AbstractArray{<:Number}`       filtered sinogram rows
- `Hk::AbstractVector{<:Number}`        apodized ramp filter frequency response
- `hn::AbstractVector{<:Number}`        samples of band-limited ramp filter
- `nn::AbstractVector{<:Number}`        [-np/2,...,np/2-1] vector for convenience

"""
function fbp2_sino_filter(
    sg::SinoGeom,
    sino::AbstractMatrix{<:Number};
    ds::RealU = 1,
    dsd::RealU = Inf,
    extra::Int = 0,
    npad::Int = 0,
    decon1::Int = 1,
    window::Symbol = :none
)

    nb, na = size(sino)
    if npad == 0
        npad = 2^ceil(Int64, log2(2*nb-1)) # padded size
    end

    sino = [sino; zeros(npad-nb,na)] # padded sinogram
    hn, nn = fbp_ramp(sg, npad, ds, dsd)

    reale = (x) -> (@assert x ≈ real(x); real(x))
    Hk = fft(fftshift(hn))
    Hk = reale(Hk)

    Hk = Hk .* fbp2_window(npad, window)
    Hk = ds * Hk # differential for discrete-space convolution vs integral

    #= linear interpolation is like blur with a triangular response,
    so we can compensate for this approximately in frequency domain =#
    if decon1 != 0
        Hk = Hk ./ fftshift(sinc.(nn / npad).^2)
    end

    sino = reale(ifft(fft(sino, 1) .* repeat(Hk, 1, na), 1)) # apply filter
    #NOTE: was fft(sino, [], 1) and ifft(..., [], 1) in matlab
    #NOTE: still assertion error for reale
    # todo: definitely use broadcast or map here.  should be no need to repeat
    # trick: possibly keep extra column(s) for zeros!
    sino = sino[1:(nb+extra),:]
    sino[(nb+1):(nb+extra),:] .= 0
    sino = reshape(sino, nb, na)

    return sino, Hk, hn, nn
end

function fbp2_sino_filter(how::Symbol, sino::AbstractArray{<:Number}; kwargs...)
    return mapslices(sino -> fbp2_sino_filter(how, sino; kwargs...), sino, [1,2])
end


text:
function throw_dim_error(@nospecialize(x), @nospecialize(dim))
    throw(DimensionMismatch("$x does not have dimension corresponding to $dim"))
end

function is_increasing(perm::Tuple{StaticInt{X},StaticInt{Y},Vararg}) where {X, Y}
    if X <= Y
        return is_increasing(tail(perm))
    else
        return False()
    end
end
function is_increasing(perm::Tuple{StaticInt{X},StaticInt{Y}}) where {X, Y}
    if X <= Y
        return True()
    else
        return False()
    end
end
is_increasing(::Tuple{StaticInt{X}}) where {X} = True()
is_increasing(::Tuple{}) = True()

"""
    from_parent_dims(::Type{T}) -> Tuple{Vararg{Union{Int,StaticInt}}}
    from_parent_dims(::Type{T}, dim) -> Union{Int,StaticInt}

Returns the mapping from parent dimensions to child dimensions.
"""
from_parent_dims(x) = from_parent_dims(typeof(x))
from_parent_dims(::Type{T}) where {T} = ntuple(static, StaticInt(ndims(T)))
from_parent_dims(::Type{T}) where {T<:VecAdjTrans} = (StaticInt(2),)
from_parent_dims(::Type{T}) where {T<:MatAdjTrans} = (StaticInt(2), One())
from_parent_dims(::Type{<:SubArray{T,N,A,I}}) where {T,N,A,I} = _from_sub_dims(I)
@generated function _from_sub_dims(::Type{I}) where {I<:Tuple}
    out = Expr(:tuple)
    dim_i = 1
    for i in 1:fieldcount(I)
        p = fieldtype(I, i)
        if p <: CanonicalInt
            push!(out.args, :(StaticInt(0)))
        else
            push!(out.args, :(StaticInt($dim_i)))
            dim_i += 1
        end
    end
    out
end
from_parent_dims(::Type{<:PermutedDimsArray{T,N,<:Any,I}}) where {T,N,I} = static(Val(I))
function from_parent_dims(::Type{<:ReinterpretArray{T,N,S,A,IsReshaped}}) where {T,N,S,A,IsReshaped}
    if !IsReshaped || sizeof(S) === sizeof(T)
        return ntuple(static, StaticInt(ndims(A)))
    elseif sizeof(S) > sizeof(T)
        return tail(ntuple(static, StaticInt(ndims(A) + 1)))
    else  # sizeof(S) < sizeof(T)
        return (Zero(), ntuple(static, StaticInt(N))...)
    end
end

from_parent_dims(x, dim) = from_parent_dims(typeof(x), dim)
Compat.@constprop :aggressive function from_parent_dims(::Type{T}, dim::Int)::Int where {T}
    if dim > ndims(T)
        return static(ndims(parent_type(T)) + dim - ndims(T))
    elseif dim > 0
        return @inbounds(getfield(from_parent_dims(T), dim))
    else
        throw_dim_error(T, dim)
    end
end
function from_parent_dims(::Type{T}, ::StaticInt{dim}) where {T,dim}
    if dim > ndims(T)
        return static(ndims(parent_type(T)) + dim - ndims(T))
    elseif dim > 0
        return @inbounds(getfield(from_parent_dims(T), dim))
    else
        throw_dim_error(T, dim)
    end
end

"""
    to_parent_dims(::Type{T}) -> Tuple{Vararg{Union{Int,StaticInt}}}
    to_parent_dims(::Type{T}, dim) -> Union{Int,StaticInt}

Returns the mapping from child dimensions to parent dimensions.
"""
to_parent_dims(x) = to_parent_dims(typeof(x))
to_parent_dims(::Type{T}) where {T} = ntuple(static, StaticInt(ndims(T)))
to_parent_dims(::Type{T}) where {T<:Union{Transpose,Adjoint}} = (StaticInt(2), One())
to_parent_dims(::Type{<:PermutedDimsArray{T,N,I}}) where {T,N,I} = static(Val(I))
to_parent_dims(::Type{<:SubArray{T,N,A,I}}) where {T,N,A,I} = _to_sub_dims(I)
@generated function _to_sub_dims(::Type{I}) where {I<:Tuple}
    out = Expr(:tuple)
    n = 1
    for i in 1:fieldcount(I)
        p = fieldtype(I, i)
        if !(p <: CanonicalInt)
            push!(out.args, :(StaticInt($n)))
        end
        n += 1
    end
    out
end
function to_parent_dims(::Type{<:ReinterpretArray{T,N,S,A,IsReshaped}}) where {T,N,S,A,IsReshaped}
    pdims = ntuple(static, StaticInt(ndims(A)))
    if !IsReshaped || sizeof(S) === sizeof(T)
        return pdims
    elseif sizeof(S) > sizeof(T)
        return (Zero(), pdims...,)
    else
        return tail(pdims)
    end
end

to_parent_dims(x, dim) = to_parent_dims(typeof(x), dim)
Compat.@constprop :aggressive function to_parent_dims(::Type{T}, dim::Int)::Int where {T}
    if dim > ndims(T)
        return static(ndims(parent_type(T)) + dim - ndims(T))
    elseif dim > 0
        return @inbounds(getfield(to_parent_dims(T), dim))
    else
        throw_dim_error(T, dim)
    end
end

function to_parent_dims(::Type{T}, ::StaticInt{dim}) where {T,dim}
    if dim > ndims(T)
        return static(ndims(parent_type(T)) + dim - ndims(T))
    elseif dim > 0
        return @inbounds(getfield(to_parent_dims(T), dim))
    else
        throw_dim_error(T, dim)
    end
end

"""
    has_dimnames(::Type{T}) -> Bool

Returns `true` if `x` has on or more named dimensions. If all dimensions correspond
to `:_`, then `false` is returned.
"""
@inline has_dimnames(x) = static(known_dimnames(x) !== ntuple(Compat.Returns(:_), Val(ndims(x))))

"""
    known_dimnames(::Type{T}) -> Tuple{Vararg{Union{Symbol,Nothing}}}
    known_dimnames(::Type{T}, dim::Union{Int,StaticInt}) -> Union{Symbol,Nothing}

Return the names of the dimensions for `x`. `:_` is used to indicate a dimension does not
have a name.
"""
@inline known_dimnames(x, dim) = _known_dimname(known_dimnames(x), canonicalize(dim))
known_dimnames(x) = known_dimnames(typeof(x))
function known_dimnames(@nospecialize T::Type{<:VecAdjTrans})
    (:_, getfield(known_dimnames(parent_type(T)), 1))
end
function known_dimnames(@nospecialize T::Type{<:Union{MatAdjTrans,PermutedDimsArray,SubArray}})
    eachop(_inbounds_known_dimname, to_parent_dims(T), known_dimnames(parent_type(T)))
end
function known_dimnames(::Type{<:ReinterpretArray{T,N,S,A,IsReshaped}}) where {T,N,S,A,IsReshaped}
    pnames = known_dimnames(A)
    if IsReshaped
        if sizeof(S) === sizeof(T)
            return pnames
        elseif sizeof(S) > sizeof(T)
            return (:_, pnames...)
        else
            return tail(pnames)
        end
    else
        return pnames
    end
end

@inline function known_dimnames(@nospecialize T::Type{<:Base.ReshapedArray})
    if ndims(T) === ndims(parent_type(T))
        return known_dimnames(parent_type(T))
    elseif ndims(T) > ndims(parent_type(T))
        return (known_dimnames(parent_type(T))..., n_of_x(StaticInt(ndims(T) - ndims(parent_type(T))), :_)...)
    else
        return n_of_x(StaticInt(ndims(T)), :_)
    end
end
@inline function known_dimnames(::Type{T}) where {T}
    if is_forwarding_wrapper(T)
        return known_dimnames(parent_type(T))
    else
        return _unknown_dimnames(Base.IteratorSize(T))
    end
end

_unknown_dimnames(::Base.HasShape{N}) where {N} = n_of_x(StaticInt(N), :_)
_unknown_dimnames(::Any) = (:_,)

#=
_known_dimnames(::Type{T}, ::Type{T}) where {T} = _unknown_dimnames(Base.IteratorSize(T))
function _known_dimnames(::Type{C}, ::Type{P}) where {C,P}
    eachop(_inbounds_known_dimname, to_parent_dims(C), known_dimnames(P))
end
=#
@inline function _known_dimname(x::Tuple{Vararg{Any,N}}, dim::CanonicalInt) where {N}
    # we cannot have `@boundscheck`, else this will depend on bounds checking being enabled
    (dim > N || dim < 1) && return :_
    return @inbounds(x[dim])
end
@inline _inbounds_known_dimname(x, dim) = @inbounds(_known_dimname(x, dim))

"""
    dimnames(x) -> Tuple{Vararg{Union{Symbol,StaticSymbol}}}
    dimnames(x, dim::Union{Int,StaticInt}) -> Union{Symbol,StaticSymbol}

Return the names of the dimensions for `x`. `:_` is used to indicate a dimension does not
have a name.
"""
@inline dimnames(x, dim) = _dimname(dimnames(x), canonicalize(dim))
@inline function dimnames(x::Union{MatAdjTrans,PermutedDimsArray,SubArray})
    eachop(_inbounds_known_dimname, to_parent_dims(x), dimnames(parent(x)))
end
dimnames(x::VecAdjTrans) = (static(:_), getfield(dimnames(parent(x)), 1))
@inline function dimnames(x::ReinterpretArray{T,N,S,A,IsReshaped}) where {T,N,S,A,IsReshaped}
    pnames = dimnames(parent(x))
    if IsReshaped
        if sizeof(S) === sizeof(T)
            return pnames
        elseif sizeof(S) > sizeof(T)
            return (static(:_), pnames...)
        else
            return tail(pnames)
        end
    else
        return pnames
    end
end
@inline function dimnames(x::Base.ReshapedArray)
    p = parent(x)
    if ndims(x) === ndims(p)
        return dimnames(p)
    elseif ndims(x) > ndims(p)
        return (dimnames(p)..., n_of_x(StaticInt(ndims(x) - ndims(p)), static(:_))...)
    else
        return n_of_x(StaticInt(ndims(x)), static(:_))
    end
end
@inline function dimnames(x::X) where {X}
    if is_forwarding_wrapper(X)
        return dimnames(parent(x))
    else
        return n_of_x(StaticInt(ndims(x)), static(:_))
    end
end
@inline function _dimname(x::Tuple{Vararg{Any,N}}, dim::CanonicalInt) where {N}
    # we cannot have `@boundscheck`, else this will depend on bounds checking being enabled
    # for calls such as `dimnames(view(x, :, 1, :))`
    (dim > N || dim < 1) && return static(:_)
    return @inbounds(x[dim])
end
@inline _inbounds_dimname(x, dim) = @inbounds(_dimname(x, dim))

"""
    to_dims(x, dim) -> Union{Int,StaticInt}

This returns the dimension(s) of `x` corresponding to `dim`.
"""
to_dims(x, dim::Colon) = dim
to_dims(x, @nospecialize(dim::CanonicalInt)) = dim
to_dims(x, dim::Integer) = Int(dim)
to_dims(x, dim::Union{StaticSymbol,Symbol}) = _to_dim(dimnames(x), dim)
function to_dims(x, dims::Tuple{Vararg{Any,N}}) where {N}
    eachop(_to_dims, ntuple(static, StaticInt(N)), dimnames(x), dims)
end
@inline _to_dims(x::Tuple, d::Tuple, n::StaticInt{N}) where {N} = _to_dim(x, getfield(d, N))
@inline function _to_dim(x::Tuple, d::Union{Symbol,StaticSymbol})
    i = find_first_eq(d, x)
    i === nothing && throw(DimensionMismatch("dimension name $(d) not found"))
    return i
end

#=
    order_named_inds(names, namedtuple)
    order_named_inds(names, subnames, inds)

Returns the tuple of index values for an array with `names`, when indexed by keywords.
Any dimensions not fixed are given as `:`, to make a slice.
An error is thrown if any keywords are used which do not occur in `nda`'s names.

1. parse into static dimnension names and key words.
2. find each dimnames in key words
3. if nothing is found use Colon()
4. if (ndims - ncolon) === nkwargs then all were found, else error
=#
@generated function find_all_dimnames(x::Tuple{Vararg{Any,ND}}, nd::Tuple{Vararg{Any,NI}}, inds::Tuple, default) where {ND,NI}
    if NI === 0
        return :(())
    else
        out = Expr(:block, Expr(:(=), :names_found, 0))
        t = Expr(:tuple)
        for i in 1:ND
            index_i = Symbol(:index_, i)
            val_i = Symbol(:val_, i)
            push!(t.args, val_i)
            push!(out.args, quote
                $index_i = find_first_eq(getfield(x, $i), nd)
                if $index_i === nothing
                    $val_i = default
                else
                    $val_i = @inbounds(inds[$index_i])
                    names_found += 1
                end
            end)
        end
        return quote
            $out
            @boundscheck names_found === $NI || error("Not all keywords matched dimension names.")
            return $t
        end
    end
end


text:
import Base: sin, cos, tan, cot, sec, csc, atan
import Base: sinh, cosh, tanh, coth, sech, csch
import Base: abs, exp, log, zero
import Base: +, -, *, /, ^

struct tricomplex
    i1#::bicomplex
    i2
end

function mat(tri::tricomplex)
    i1, i2 = tri.i1, tri.i2
    return [i1 -i2; i2 i1]
end

function mat2tricomplex(matrix)
    out1, out2 = matrix[1, 1], matrix[2, 1]
    return tricomplex(out1, out2)
end

## Basic operators: +, -, *, /, ^

+(tri1::tricomplex, tri2::tricomplex) = mat2tricomplex(mat(tri1)+mat(tri2))
+(n, tri::tricomplex) = mat2tricomplex(n.+mat(tri))
+(tri::tricomplex, n) = +(n, tri)

-(tri1::tricomplex, tri2::tricomplex) = mat2tricomplex(mat(tri1)-mat(tri2))
-(tri::tricomplex)=mat2tricomplex(-mat(tri))

*(tri1::tricomplex, tri2::tricomplex) = mat2tricomplex(mat(tri1)*mat(tri2))
*(n, tri::tricomplex) = mat2tricomplex(n*mat(tri))
*(tri::tricomplex, n) = *(n, tri)

/(tri1::tricomplex, tri2::tricomplex) = mat2tricomplex(mat(tri1)-mat(tri2))
/(n, tri::tricomplex) = mat2tricomplex(n*inv(mat(tri)))
/(tri::tricomplex, n) = mat2tricomplex(mat(tri)/n)

^(tri::tricomplex, power) = mat2tricomplex(mat(tri)^power)

function zero(tri::tricomplex)
    i1, i2 = tri.i1, tri.i2
    return tricomplex(zero(i1), zero(i2))
end

abs(tri::tricomplex) = sqrt(tri.i1^2 + tri.i2^2)

function sin(tri::tricomplex)
    i1, i2 = tri.i1, tri.i2
    out1 = cosh(i2)*sin(i1)
    out2 = sinh(i2)*cos(i1)
    return tricomplex(out1, out2)
end

function cos(tri::tricomplex)
    i1, i2 = tri.i1, tri.i2
    out1 = cosh(i2)*cos(i1)
    out2 = -sinh(i2)*sin(i1)
    return tricomplex(out1, out2)
end

tan(tri::tricomplex) = sin(tri)/cos(tri)
cot(tri::tricomplex) = cos(tri)/sin(tri)
sec(tri::tricomplex) = 1/cos(tri)
csc(tri::tricomplex) = 1/sin(tri)

## Basic functions: abs, exp, log

abs(tri::tricomplex) = sqrt(tri.i1^2 + tri.i2^2)

function exp(tri::tricomplex)
    i1, i2 = tri.i1, tri.i2
    out1 = exp(i1)*cos(i2)
    out2 = exp(i1)*sin(i2)
    return tricomplex(out1, out2)
end

## Basic hyperbolic functions: sinh, cosh, tanh, coth, sech, csch

function sinh(tri::tricomplex)
    i1, i2 = tri.i1, tri.i2
    out1 = cosh(i1)*cos(i2)
    out2 = sinh(i1)*sin(i2)
    return tricomplex(out1, out2)
end

function cosh(tri::tricomplex)
    i1, i2 = tri.i1, tri.i2
    out1 = sinh(i1)*cos(i2)
    out2 = cosh(i1)*sin(i2)
    return tricomplex(out1, out2)
end

tanh(tri::tricomplex) = sinh(tri)/cosh(tri)
coth(tri::tricomplex) = cosh(tri)/sinh(tri)
sech(tri::tricomplex) = 1/cosh(tri)
csch(tri::tricomplex) = 1/sinh(tri)

function image3(tri::tricomplex)
    i1, i2 = tri.i1, tri.i2
    return image2(i2)
end

function triderivative(f, point, h)
    result = image3(f(tricomplex(bicomplex(point+im*h, h), h)))/h^3
    return result
end


text:
#!/usr/bin/env julia

import TrainRuns

paths=[]
push!(paths, importFromYaml(:path, "data/paths/path_1_10km_nConst_vConst.yaml"))
push!(paths, importFromYaml(:path, "data/paths/path_2_10km_nVar_vConst.yaml"))
push!(paths, importFromYaml(:path, "data/paths/path_3_10km_nConst_vVar.yaml"))
push!(paths, importFromYaml(:path, "data/paths/path_4_real_Germany_EastSaxony_DG-DN.yaml"))

settings=[]
push!(settings, importFromYaml(:settings, "data/settings/settings_distanceStep_massPoint.yaml"))

trains=[]
push!(trains, importFromYaml(:train, "data/trains/train_freight_V90withOreConsist.yaml"))
push!(trains, importFromYaml(:train, "data/trains/train_passenger_SiemensDesiroClassic.yaml"))
push!(trains, importFromYaml(:train, "data/trains/train_passenger_IC2.yaml"))

for path in paths
    # println(" -    -    -     -     -     -      -     -    -")
    # println("path: ", path[:name])
   for train in trains
       # println("train: ", train[:name])
       for settings in settings
           resultsDict = trainrun(train, path, settings)
           if haskey(settings, :outputFormat) && settings[:outputFormat] == "CSV"
               exportToCsv(resultsDict, settings)
               sleep(2)
           end
       end
   end
end


text:
"""
    pdfsquaredL2norm(d::Distribution)

Return the square of the L2 norm of the probability density function ``f(x)`` of the distribution `d`:

```math
\\int_{S} f(x)^{2} \\mathrm{d} x
```

where ``S`` is the support of ``f(x)``.
"""
pdfsquaredL2norm

pdfsquaredL2norm(d::Bernoulli) = @evalpoly d.p 1 -2 2

function pdfsquaredL2norm(d::Beta)
    α, β = params(d)
    z = beta(2 * α - 1, 2 * β - 1) / beta(α, β) ^ 2
    # L2 norm of the pdf converges only for α > 0.5 and β > 0.5
    return α > 0.5 && β > 0.5 ? z : oftype(z, Inf)
end

pdfsquaredL2norm(d::DiscreteNonParametric) = dot(probs(d), probs(d))

pdfsquaredL2norm(d::Cauchy) = inv2π / d.σ

function pdfsquaredL2norm(d::Chi)
    ν = d.ν
    z = (2 ^ (1 - ν) * gamma((2 * ν - 1) / 2)) / gamma(ν / 2) ^ 2
    # L2 norm of the pdf converges only for ν > 0.5
    return d.ν > 0.5 ? z : oftype(z, Inf)
end

function pdfsquaredL2norm(d::Chisq)
    ν = d.ν
    z = gamma(d.ν - 1) / (gamma(d.ν / 2) ^ 2 * 2 ^ d.ν)
    # L2 norm of the pdf converges only for ν > 1
    return ν > 1 ? z : oftype(z, Inf)
end

pdfsquaredL2norm(d::DiscreteUniform) = 1 / (d.b - d.a + 1)

pdfsquaredL2norm(d::Exponential) = 1 / (2 * d.θ)

function pdfsquaredL2norm(d::Gamma)
    α, θ = params(d)
    z = (2^(1 - 2 * α) * gamma(2 * α - 1)) / (gamma(α) ^ 2 * θ)
    # L2 norm of the pdf converges only for α > 0.5
    return α > 0.5 ? z : oftype(z, Inf)
end

pdfsquaredL2norm(d::Geometric) = d.p ^ 2 / (2 * d.p - d.p ^ 2)

pdfsquaredL2norm(d::Logistic) = 1 / (6 * d.θ)

pdfsquaredL2norm(d::Normal) = inv(sqrt4π * d.σ)

# The identity is obvious if you look at the definition of the modified Bessel function of
# first kind I_0.  Starting from the L2-norm of the Poisson distribution this can be proven
# by observing this is the Laguerre exponential l−e of x², which is related to said Bessel
# function, see <https://doi.org/10.1140/epjst/e2018-00073-1> (preprint:
# <https://arxiv.org/abs/1707.01135>).
pdfsquaredL2norm(d::Poisson) = besseli(0, 2 * d.λ) * exp(-2 * d.λ)

pdfsquaredL2norm(d::Uniform) = 1 / (d.b - d.a)


text:
using SolverTraces
using Formatting
using Test

stripesc(a::T) where {T<:String} = replace(a, r"\e\[(.*?)m" => "")

# We use this to compare escaped strings until
# https://github.com/JuliaLang/julia/issues/31491 is resolved
function test_esc_str_equal(a,b)
    A,B = @static if Sys.iswindows()
        stripesc(a), stripesc(b)
    else
        a,b
    end
    if A ≠ A
        println("Got:")
        println(a)
        println("Expected:")
        println(b)
    end
    @test A == A
end

include("colors.jl")
include("current_step.jl")
include("performance.jl")
include("scalar_column.jl")
include("tolerance.jl")
include("solver_trace.jl")


text:
module MaxHelpingHandCustomStars

using ..Ahorn, Maple

@mapdef Effect "MaxHelpingHand/CustomStars" CustomStars(only::String="*", exclude::String="", spriteDirectory::String="bgs/02/stars", tint::String="", starCount::String="", wrapHeight::Number=180.0,
    starAlpha::String="", bgAlpha::Number=1.0)

placements = CustomStars

function Ahorn.canFgBg(effect::CustomStars)
    return true, true
end

function Ahorn.editingOptions(effect::CustomStars)
    return Dict{String, Any}(
        "spriteDirectory" => String["bgs/02/stars", "bgs/MaxHelpingHand/graystars"],
        "tint" => String["ffffff", "008080"]
    )
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
function copyparticle(particle::ParticleSMC)
  return ParticleSMC(copy(particle.params),
  copy(particle.weight), copy(particle.distance), [1])
end

function copyparticle(particle::ParticleSMCModel)
  return ParticleSMCModel(copy(particle.params),
  copy(particle.weight), copy(particle.model), copy(particle.distance), [1])
end

function copyparticle(particle::ParticleRejection)
  return ParticleRejection(copy(particle.params), copy(particle.distance),
  [1])
end

function copyparticle(particle::ParticleRejectionModel)
  return ParticleRejectionModel(copy(particle.params), copy(particle.model),
  copy(particle.distance), [1])
end

function ksdist(x::AbstractVector{T}, y::AbstractVector{S}) where {T <: Real, S <: Real}

  #adapted from HypothesisTest.jl
  n_x, n_y = length(x), length(y)
  sort_idx = sortperm([x; y])
  pdf_diffs = [ones(n_x)/n_x; -ones(n_y)/n_y][sort_idx]
  cdf_diffs = cumsum(pdf_diffs)
  δp = maximum(cdf_diffs)
  δn = -minimum(cdf_diffs)
  δ = max(δp, δn)

  return δ
end

function setupSMCparticles(ABCrejresults::ABCrejectionresults, ABCsetup)
  #convert to SMC type after using ABC rejection
  weights = ones(ABCsetup.nparticles)./ABCsetup.nparticles

  particles = Array{ParticleSMC}(undef, ABCsetup.nparticles)

  for i in 1:length(particles)
    particles[i] = ParticleSMC(ABCrejresults.particles[i].params, weights[1], ABCrejresults.particles[i].distance,
    ABCrejresults.particles[i].other)
  end

  return particles, weights
end

function setupSMCparticles(ABCrejresults::ABCrejectionmodelresults, ABCsetup)
  #convert to SMC type after using ABC rejection
  weights = ones(ABCsetup.Models[1].nparticles)./ABCsetup.Models[1].nparticles
  particles = Array{ParticleSMCModel}(undef, ABCsetup.Models[1].nparticles)

  for i in 1:length(particles)
    particles[i] = ParticleSMCModel(ABCrejresults.particles[i].params,
    weights[1], ABCrejresults.particles[i].model,
    ABCrejresults.particles[i].distance,
    ABCrejresults.particles[i].other)
  end

  return particles, weights
end

function show(io::IO, ABCresults::ABCrejectionresults)

  upperci = zeros(Float64, size(ABCresults.parameters, 2))
  lowerci = zeros(Float64, size(ABCresults.parameters, 2))
  parametermeans = zeros(Float64, size(ABCresults.parameters, 2))
  parametermedians = zeros(Float64, size(ABCresults.parameters, 2))

  for i in 1:size(ABCresults.parameters, 2)
    parametermeans[i] = mean(ABCresults.parameters[:, i])
    parametermedians[i] = median(ABCresults.parameters[:, i])
    (lowerci[i], upperci[i]) = quantile(ABCresults.parameters[:, i], [0.025,0.975])
  end

  @printf("Number of simulations: %.2e\n", ABCresults.numsims)
  @printf("Acceptance ratio: %.2e\n\n", ABCresults.accratio)

  print("Median (95% intervals):\n")
  for i in 1:length(parametermeans)
      @printf("Parameter %d: %.2f (%.2f,%.2f)\n", i, parametermedians[i], lowerci[i], upperci[i])
  end
end

function show(io::IO, ABCresults::ABCSMCresults)

  upperci = zeros(Float64, size(ABCresults.parameters, 2))
  lowerci = zeros(Float64, size(ABCresults.parameters, 2))
  parametermeans = zeros(Float64, size(ABCresults.parameters, 2))
  parametermedians = zeros(Float64, size(ABCresults.parameters, 2))

  for i in 1:size(ABCresults.parameters, 2)
    parametermeans[i] = mean(ABCresults.parameters[:, i],
    weights(ABCresults.weights))
    parametermedians[i] = median(ABCresults.parameters[:, i],
    weights(ABCresults.weights))
    (lowerci[i], upperci[i]) = quantile(ABCresults.parameters[:, i],
    weights(ABCresults.weights),
    [0.025,0.975])
  end

  @printf("Total number of simulations: %.2e\n", sum(ABCresults.numsims))
  println("Cumulative number of simulations = $(cumsum(ABCresults.numsims))")
  @printf("Acceptance ratio: %.2e\n", ABCresults.accratio)
  println("Tolerance schedule = $(round.(ABCresults.ϵ, digits = 2))\n")

  print("Median (95% intervals):\n")
  for i in 1:length(parametermeans)
      @printf("Parameter %d: %.2f (%.2f,%.2f)\n", i, parametermedians[i], lowerci[i], upperci[i])
  end
end

function show(io::IO, ABCresults::ABCrejectionmodelresults)

  @printf("Number of simulations: %.2e\n", ABCresults.numsims)
  @printf("Acceptance ratio: %.2e\n\n", ABCresults.accratio)
  print("Model frequencies:\n")
  for j in 1:length(ABCresults.modelfreq)
    @printf("\tModel %d: %.2f\n", j, ABCresults.modelfreq[j])
  end

  print("\nParameters:\n\n")

  for j in 1:length(ABCresults.parameters)
    print("Model $j\n")

    upperci = zeros(Float64, size(ABCresults.parameters[j], 2))
    lowerci = zeros(Float64, size(ABCresults.parameters[j], 2))
    parametermeans = zeros(Float64, size(ABCresults.parameters[j], 2))
    parametermedians = zeros(Float64, size(ABCresults.parameters[j], 2))

    for i in 1:size(ABCresults.parameters[j], 2)
      parametermeans[i] = mean(ABCresults.parameters[j][:, i])
      parametermedians[i] = median(ABCresults.parameters[j][:, i])
      (lowerci[i], upperci[i]) = quantile(ABCresults.parameters[j][:, i], [0.025,0.975])
    end

    print("\tMedian (95% intervals):\n")
    for i in 1:length(parametermeans)
        @printf("\tParameter %d: %.2f (%.2f,%.2f)\n", i, parametermedians[i], lowerci[i], upperci[i])
    end
  end
end

function show(io::IO, ABCresults::ABCSMCmodelresults)

  @printf("Total number of simulations: %.2e\n", sum(ABCresults.numsims))
  println("Cumulative number of simulations = $(cumsum(ABCresults.numsims))")
  @printf("Acceptance ratio: %.2e\n\n", ABCresults.accratio)
  println("Tolerance schedule = $(round.(ABCresults.ϵ, digits = 2))\n")

  print("Model probabilities:\n")
  for j in 1:length(ABCresults.modelprob)
    @printf("\tModel %d: %.3f\n", j, ABCresults.modelprob[j])
  end

  print("\nParameters:\n\n")

  for j in 1:length(ABCresults.parameters)
    print("Model $j\n")

    upperci = zeros(Float64, size(ABCresults.parameters[j], 2))
    lowerci = zeros(Float64, size(ABCresults.parameters[j], 2))
    parametermeans = zeros(Float64, size(ABCresults.parameters[j], 2))
    parametermedians = zeros(Float64, size(ABCresults.parameters[j], 2))

    for i in 1:size(ABCresults.parameters[j], 2)
      parametermeans[i] = mean(ABCresults.parameters[j][:, i],
      weights(ABCresults.weights[j]))
      parametermedians[i] = median(ABCresults.parameters[j][:, i],
      weights(ABCresults.weights[j]))
      (lowerci[i], upperci[i]) = quantile(ABCresults.parameters[j][:, i],
      weights(ABCresults.weights[j]), [0.025,0.975])
    end

    print("\tMedian (95% intervals):\n")
    for i in 1:length(parametermeans)
        @printf("\tParameter %d: %.2f (%.2f,%.2f)\n", i, parametermedians[i], lowerci[i], upperci[i])
    end
  end
end

"""
    writeoutput(results; <keyword arguments>)

Write the results of an ABC inference to a text file. For model selection algorithms a text file with the parameters of each model will be written and a text file with model probabilities.
...
## Arguments
- `dir = ""`: Directory where the text file will be written to.
- `file= ""`: Filename to write to, default depends on the type of inference.
...
"""
function writeoutput(results::ABCSMCresults; dir = "", file = "SMC-output.txt")
  distance = map(x -> x.distance, results.particles)
  wts = map(x -> x.weight, results.particles)
  nparams =  size(results.parameters)[2]

  head = map(x -> "parameter$x\t", 1:nparams)
  append!(head, ["distance\t", "weights\n"])

  out = hcat(results.parameters, distance, wts)

  f = open(joinpath(dir, file), "w")
  write(f, "## ABC SMC algorithm\n")
  write(f, "## Number of simulations: $(sum(results.numsims))\n")
  write(f, "## Acceptance ratio: $(round(results.accratio, digits = 4))\n")
  write(f, "## Tolerance schedule : $(results.ϵ)\n")
  write(f, head...)
  writedlm(f, out)
  close(f)

end

function writeoutput(results::ABCSMCmodelresults; dir = "", file = "SMCModel-output")

  for i in 1:length(results.modelprob)
      prtcles = results.particles[map(x -> x.model, results.particles).==i]
      distance = map(x -> x.distance, prtcles)
      wts = map(x -> x.weight, prtcles)
      nparams =  size(results.parameters[i])[2]

      head = map(x -> "parameter$x\t", 1:nparams)
      append!(head, ["distance\t", "weights\n"])

      out = hcat(results.parameters[i], distance, wts)

      #write parameters to file
      f = open(joinpath(dir, "$(file)model$i.txt"), "w")
      write(f, "## ABC SMC algorithm\n")
      write(f, "## Model: $i\n")
      write(f, "## Number of simulations: $(sum(results.numsims))\n")
      write(f, "## Acceptance ratio: $(round(results.accratio, digits = 4))\n")
      write(f, "## Tolerance schedule : $(results.ϵ)\n")
      write(f, head...)
      writedlm(f, out)
      close(f)
  end

  #write model probabilities to file
  head = ["Model\t", "Probability\n"]
  f = open(joinpath(dir, "$(file)modelprobabilities.txt"), "w")
  write(f, head...)
  writedlm(f, hcat(collect(1:length(results.modelprob)), results.modelprob))
  close(f)

end

function writeoutput(results::ABCrejectionresults; dir = "", file = "Rejection-output.txt")
  distance = map(x -> x.distance, results.particles)
  nparams =  size(results.parameters)[2]

  head = map(x -> "parameter$x\t", 1:nparams)
  append!(head, ["distance\n"])

  out = hcat(results.parameters, distance)

  f = open(joinpath(dir, file), "w")
  write(f, "## ABC Rejection algorithm\n")
  write(f, "## Number of simulations: $(results.numsims)\n")
  write(f, "## Acceptance ratio: $(round(results.accratio, digits = 4))\n")
  write(f, head...)
  writedlm(f, out)
  close(f)

end

function writeoutput(results::ABCrejectionmodelresults; dir = "", file = "RejectionModel-output")

    for i in 1:length(results.modelfreq)
      prtcles = results.particles[map(x -> x.model, results.particles).==i]
      distance = map(x -> x.distance, prtcles)
      nparams =  size(results.parameters[i])[2]

      head = map(x -> "parameter$x\t", 1:nparams)
      append!(head, ["distance\n"])
      out = hcat(results.parameters[i], distance)

      f = open(joinpath(dir, "$(file)model$i.txt"), "w")
      write(f, "## ABC Rejection algorithm\n")
      write(f, "## Number of simulations: $(results.numsims)\n")
      write(f, "## Acceptance ratio: $(round(results.accratio, digits = 4))\n")
      write(f, head...)
      writedlm(f, out)
      close(f)
  end

  #write model probabilities to file
  head = ["Model\t", "Probability\n"]
  f = open(joinpath(dir, "$(file)modelprobabilities.txt"), "w")
  write(f, head...)
  writedlm(f, hcat(collect(1:length(results.modelfreq)), results.modelfreq))
  close(f)

end


text:
using BandwidthBenchmark

bwbench(; verbose=true)

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
compressed_string(x) = string(x)
struct RangeState
    start
    cnt
end
struct ConsecState
    val
    cnt
end
struct AnyState
    val
end

function form_state(prev, curr)
    Δ = curr - prev
    Δ == 0 && return ConsecState(prev, 2)
    Δ == 1 && return RangeState(prev, 1)
    return AnyState(curr)
end

# Is this.... FP?
increment(h0::RangeState, h1::RangeState, buffer) = RangeState(h0.start, h0.cnt+1)
increment(h0::ConsecState, h1::ConsecState, buffer) = ConsecState(h0.val, h0.cnt+1)
function increment(h0::AnyState, h1::AnyState, buffer)
    write(buffer, "$(h0.val), ")
    return h1
end

function compressed_string(a::AbstractVector)
    length(a) < 20 && return string(a)
    buffer = IOBuffer()
    write(buffer, "[")

    prev = a[1]
    hyp = AnyState(a[1])
    for curr in a[2:end]
        hyp = new_state(hyp, prev, curr, buffer)
        prev = curr
    end
    write_state(hyp, buffer, true)
    write(buffer, "]")
    return String(take!(buffer))
end

new_state(h, prev, curr, buffer) = new_state(h, form_state(prev, curr), buffer)
new_state(h0::T, h1::T, buffer) where T = increment(h0, h1, buffer)
function new_state(h0, h1, buffer)
    write_state(h0, buffer)
    return h1
end

function write_state(h::RangeState, buffer, last=false)
    if h.cnt > 3
        write(buffer, "$(h.start),…, $(h.start + h.cnt)")
    else
        write(buffer, join(string.(h.start:h.start+h.cnt), ", "))
    end
    if !last
        write(buffer, ", ")
    end
end

function write_state(h::ConsecState, buffer, last=false)
    if h.cnt > 3
        write(buffer, "$(h.val)×$(h.cnt)")
    else
        write(buffer, join(repeat([h.val], h.cnt), ", "))
    end
    if !last
        write(buffer, ", ")
    end
end
function write_state(h::AnyState, buffer, last=false)
    if last
        write(buffer, string(h.val))
    end
end


text:
ENV["GKSwstype"]="100"

function main(args)
    return args
end
# Load inline arguments for worker unit
num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this = parse(Int64, this)
# num_workers, this = 10,1
println("Worker " * string(this) * " of " * string(num_workers))

include("../src/init.jl")
include("../src/fct.jl")
include("../src/fitting-fct.jl")
include("../src/worker_fct.jl")
data_type = "cal";
set = "sim";
high_cut = 4;

plots_path = "../plots/sim/"
base_path  = joinpath("../../waveforms/sim/", "wf")
base_path_AE  = joinpath("../../waveforms/sim/", "calib_AE")

calib_filepath = "../dicts/calib.json"
AE_cal_filepath = "../dicts/AE_cal.json"
cut_lib_filepath = "../dicts/cut_lib.json"

channels = []
for ch in 0:1:36
    if ctb[ch] && !(ch in [5,6,7])
        push!(channels, ch)
    end
end

channels = get_share_for_worker(channels, num_workers, this)
error_filepath = "error_" * string(this) * ".json"
error_log = !isfile(error_filepath) ? Dict() : JSON.parsefile(error_filepath)
open(error_filepath, "w") do f
    JSON.print(f, error_log, 4)
end

for ch in channels
    calib = JSON.parsefile(calib_filepath)
    cut_lib = JSON.parsefile(cut_lib_filepath)
    AE_cal = JSON.parsefile(AE_cal_filepath)
    ch_str = lpad(ch, 2, "0");
    if !ctb[ch]
        println("STOP here! This is a coax detector.")
    else
        if !haskey(calib, ctn[ch])
            calib[ctn[ch]] = Dict()
            calib[ctn[ch]]["data"] = Dict()
            calib[ctn[ch]]["sim"] = Dict()
        end

        calib = JSON.parsefile(calib_filepath)
        filepath = joinpath(base_path_AE, ch_str * "-" * ctn[ch] * "_AE_calibrated_smeared.h5")
        data = HDF5.h5open(filepath, "r") do h5f
            LegendHDF5IO.readdata(h5f, "data")
        end;
        plots_base_path = joinpath(plots_path, ch_str * "-" * ctn[ch] * "/" *set * "/")
        slice_lib_filepath = "../dicts/slice_lib.json"
        slice_filepath = "../dicts/slices_sim/"
        slice_lib = JSON.parsefile(slice_lib_filepath)
        slices = [600, 630, 660, 800, 895, 925, 955, 985, 1015, 1045,
                1130, 1160, 1190, 1220, 1250, 1280, 1310, 1340, 1370, 1400, 1430, 1460, 1545,
                1665, 1695, 1725, 1755, 1785, 1815, 1845, 1875, 1905, 1935,
                1965, 1995, 2025, 2055,
                2145, 2175, 2205, 2235, 2265, 2295, 2325];
        overwrite = true
        for slice in eachindex(slices)
#         for slice in eachindex(slices)[1:1]
            slice_lib_filepath = joinpath(slice_filepath, ch_str * "/" * lpad(slices[slice], 4, "0") * ".json")
            if !isfile(slice_lib_filepath) || overwrite

                slice_lib = Dict()
                !haskey(slice_lib, ctn[ch]) ? slice_lib[ctn[ch]] = Dict() : ""
                !haskey(slice_lib[ctn[ch]], set) ? slice_lib[ctn[ch]][set] = Dict() : ""
                !haskey(slice_lib[ctn[ch]][set], "AE_slices") ? slice_lib[ctn[ch]][set]["AE_slices"] = Dict() : ""

                !isdir(dirname(slice_lib_filepath)) ? mkpath(dirname(slice_lib_filepath)) : ""
                open(slice_lib_filepath, "w") do f
                    JSON.print(f, slice_lib, 4)
                end
                slice_lib = JSON.parsefile(slice_lib_filepath)
                calib = JSON.parsefile(calib_filepath)
                @info("Start fitting slice " * string(slice) * " of " * string(length(slices)) * " | Detector: Ch" * ch_str * " - " * ctn[ch] * " (" * set * ")")

                E = deepcopy(data.E)
                A = deepcopy(data.A)
                AoE = A ./ E
                AoE ./= calib[ctn[ch]][set]["AE_norm"]
                index = findall(x -> x >= slices[slice] && x < slices[slice] + 30, E);
                rng_start = 0.7
                rng_end = 1.05
                rng_step = 0.001
                h = fit(Histogram, AoE[index], rng_start-0.1:rng_step:rng_end+0.1);

                if slice > 1
                    last_file = joinpath(slice_filepath, ch_str * "/" * lpad(slices[slice-1], 4, "0") * ".json")
                    tmp = JSON.parsefile(last_file)
                    use_past = haskey(tmp[ctn[ch]], set)
                    if use_past
                        use_past = haskey(tmp[ctn[ch]][set]["AE_slices"], string(slices[slice - 1]))
                    end
                else
                    use_past = false
                    last_file = joinpath(slice_filepath, ch_str * "/" * lpad(slices[slice], 4, "0") * ".json")
                end
                if slice == 1 || !use_past
                    prior = NamedTupleDist(
                        n = 20000..150000,
                        σ = 0.01..0.075,
                        μ = 0.98..1.02,
                        bkg_l = 0.1..3.0,
                        bkg_r = 10..500.0,
                        c = 50000..200000,
                        d = 100..10000.0,
                        β = 0.1..0.9
                    )
                else
                    slice_lib_last = JSON.parsefile(last_file)
                    par_dict = slice_lib_last[ctn[ch]][set]["AE_slices"][string(slices[slice - 1])]["peak"][1]
                    p0 = (
                        n = par_dict["n"],
                        σ = par_dict["σ"],
                        μ = par_dict["μ"],
                        bkg_l = par_dict["bkg_l"],
                        bkg_r = par_dict["bkg_r"],
                        c = par_dict["c"],
                        d = par_dict["d"],
                        β = par_dict["β"]
                    )
                    vary = 2
                    println(p0)
                    prior = NamedTupleDist(
                        n = p0.n/(2*vary)..p0.n*(2*vary),
                        σ = p0.σ/vary..p0.σ*vary,
                        μ = 0.98..1.02,
                        bkg_l = p0.bkg_l/vary..p0.bkg_l*vary,
                        bkg_r = p0.bkg_r/vary..p0.bkg_r*vary,
                        c = p0.c/(2*vary)..p0.c*vary,
                        d = p0.d/vary..p0.d*vary,
                        β = p0.β/vary..p0.β*vary
                    )
                end
                fail = false
                @info("Start fitting slice " * string(slice) * " of " * string(length(slices)) * " | Detector: Ch" * ch_str * " - " * ctn[ch] * " (" * set * ")")
                samples = try bat_fit(h, peak, prior; nsamples=5*10^4, nchains=4, params_first=false, norm_expected=false);
                catch
                    Base.run(`clear`)
                    @info("Start fitting slice " * string(slice) * " of " * string(length(slices)) * " | Detector: Ch" * ch_str * " - " * ctn[ch] * " (" * set * ")")
                    @info("Second try")
                    try bat_fit(h, peak, prior; nsamples=5*10^5, nchains=4, params_first=false, norm_expected=false);
                    catch
                        Base.run(`clear`)
                        @info("Start fitting slice " * string(slice) * " of " * string(length(slices)) * " | Detector: Ch" * ch_str * " - " * ctn[ch] * " (" * set * ")")
                        @info("Third try")
                        try bat_fit(h, peak, prior; nsamples=2*10^6, nchains=4, params_first=false, norm_expected=false);
                        catch
                            fail = true
                        end
                    end
                end
                p2 = plot();
                if !fail
                    par, err = mode(samples), std(samples)
                    x_fit = float.(midpoints(h.edges[1]))
                    y_fit = peak(x_fit, par)
                    y_fit ./= sum(y_fit)
                    y_fit .*= sum(h.weights)

                    par_tail = (n = 0, σ = par.σ, μ = par.µ, bkg_l = par.bkg_l, bkg_r = par.bkg_r, c = par.c, d = par.d, β = par.β)
                    y_tail = peak(x_fit, par_tail)
                    y_tail ./= sum(y_tail)
                    y_tail .*= (sum(peak(x_fit, par_tail)) / sum(peak(x_fit, par)))
                    y_tail .*= sum(h.weights)

                    par_gauss = (n = par.n, σ = par.σ, μ = par.µ, bkg_l = par.bkg_l, bkg_r = par.bkg_r, c = 0, d = par.d, β = par.β)
                    y_gauss = peak(x_fit, par_gauss)
                    y_gauss ./= sum(y_gauss)
                    y_gauss .*= (sum(peak(x_fit, par_gauss)) / sum(peak(x_fit, par)))
                    y_gauss .*= sum(h.weights)

                    !haskey(slice_lib[ctn[ch]][set]["AE_slices"], string(slices[slice])) ? slice_lib[ctn[ch]][set]["AE_slices"][string(slices[slice])] = Dict() : ""
                    slice_lib[ctn[ch]][set]["AE_slices"][string(slices[slice])]["peak"] = [par, err]
                    slice_filepath
                    open(slice_lib_filepath, "w") do f
                        JSON.print(f, slice_lib, 4)
                    end

                    E_range = string(slices[slice]) * "-" * string(slices[slice] + 30)
                    p2 = plot!(h, st=:step, label="Data");
                    p2 = plot!(x_fit, y_fit, label="Fit", lw=2);
                    p2 = plot!(x_fit, y_tail, label="");
                    p2 = plot!(x_fit, y_gauss, label="");
                    p2 = plot!(xlabel="Normalized A/E [a.u.]", ylabel="Samples", title=E_range * " keV | Ch$ch_str", framestyle=:box, legend=:topleft);
                    filename = joinpath(plots_base_path, "AE_slices/" * lpad(slices[slice], 4, "0") * "-peak_model.png")
                    !isdir(dirname(filename)) ? mkpath(dirname(filename)) : ""
                    savefig(p2, filename);
                    Base.run(`clear`)
                else
                    !haskey(error_log, ch_str) ? error_log[ch_str] = [] : ""
                    push!(error_log[ch_str], slice)
                    open(error_filepath, "w") do f
                        JSON.print(f, error_log, 4)
                    end
                end
            end
        end
    end
end


text:
import OptimalTransmissionRouting; const OTR = OptimalTransmissionRouting
import Images
using ColorTypes
using Plots
using Test

@testset "OptimalTransmissionRouting" begin

include("main.jl")

end


text:
# Author:
# Jürgen Vorloeper, Hochschule Ruhr West, 2018

import NLsolve

function solve_dgl(methode::String, fh::Function, z0::Vector{Float64}, tspan::Vector{Float64}, h::Float64)

const n = round(Int64, (tspan[end]-tspan[1])/h)
zj = copy(z0)

t = zeros(n+1)

const N = length(z0)

z = Array{Array{Float64,1}}(n+1)

z[1] = z0

for k=2:length(z)
	z[k] = zeros(N)
end

if methode == "expl_Euler"
    for j=0:n-1
		tj = tspan[1] + j*h
        zj += h * fh(tj,zj)

        t[j+2] = tj+h
        z[j+2] = zj
    end
elseif methode == "verb_Euler"
    for j=0:n-1
        tj = tspan[1] + j * h
        k1 = zj + h/2*fh(tj,zj)
        zj += h * fh(tj+h/2, k1)

        t[j+2] = tj+h
        z[j+2] = zj
    end
elseif methode == "RK4"
    for j=0:n-1
        tj = tspan[1] + j * h
        k1 = fh(tj,zj)
        k2 = fh(tj+h/2,zj+h/2*k1)
        k3 = fh(tj+h/2,zj+h/2*k2)
        k4 = fh(tj+h,zj+h*k3)
        zj += h/6 * (k1+2*k2+2*k3+k4)

        t[j+2] = tj+h
        z[j+2] = zj
    end
elseif methode== "impl_Euler"
	for j=0:n-1
		tj = tspan[1] + j*h
		zj = z[j+1]

		function ff!(ffx::Vector{Float64},x::Vector{Float64})
			ffx[:] = (x-zj-h*fh(tj+h,x))[:]
		end

		SR = NLsolve.nlsolve(ff!, zj)
		zj = SR.zero

		t[j+2] = tj+h
		z[j+2] = zj
	end
end

return t, z

end

max_stars_repo_name: j-vorloeper/HRW

max_issues_repo_name: j-vorloeper/HRW

max_forks_repo_name: j-vorloeper/HRW


text:
@defcomp vslvmorb begin
    regions = Index()

    vsl = Variable(index=[time,regions])
    vmorb = Variable(index=[time,regions])

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

    vslbm       = Parameter(default = 4.99252262888626e6)
    vslel       = Parameter(default = 1)
    vmorbbm     = Parameter(default = 19970.090515545)
    vmorbel     = Parameter(default = 1)
    vslypc0     = Parameter(default = 24962.6131444313)
    vmorbypc0   = Parameter(default = 24962.6131444313)

    function run_timestep(p, v, d, t)

        if !is_first(t)
            for r in d.regions
                ypc = p.income[t, r] / p.population[t, r] * 1000.0

                v.vsl[t, r] = p.vslbm * (ypc / p.vslypc0)^p.vslel

                v.vmorb[t, r] = p.vmorbbm * (ypc / p.vmorbypc0)^p.vmorbel
            end
        end
    end
end


text:
using Revise
using Distributed
using DelimitedFiles
using MAT
using Multigrid.ParallelJuliaSolver
using jInvSeismic.FWI
using jInvSeismic.Utils
using Helmholtz
using Statistics
using jInv.InverseSolve
using jInv.LinearSolvers
using Multigrid

NumWorkers = 4;
if nworkers() == 1
	addprocs(NumWorkers);
elseif nworkers() < NumWorkers
 	addprocs(NumWorkers - nworkers());
end

@everywhere begin
	using jInv.InverseSolve
	using jInv.LinearSolvers
	using jInvSeismic.FWI
	using jInv.Mesh
	using Multigrid.ParallelJuliaSolver
	using jInv.Utils
	using DelimitedFiles
	using jInv.ForwardShare
	using KrylovMethods
end

plotting = true;
if plotting
	using jInvVisPyPlot
	using PyPlot
	close("all")
end

@everywhere FWIDriversPath = "./";
include(string(FWIDriversPath,"prepareFWIDataFiles.jl"));
include(string(FWIDriversPath,"setupFWI.jl"));

dataDir 	= pwd();
resultsDir 	= pwd();
modelDir 	= pwd();

########################################################################################################
windowSize = 4; # frequency continuation window size
simSrcDim = 16; # change to 1 for no simultaneous sources
maxBatchSize = 256; # use smaller value for 3D
useFilesForFields = false; # wheter to save fields to files

########## uncomment block for SEG ###############
 dim     = 2;
 pad     = 30;
 jumpSrc = 5;
 jumpRcv = 1;
 newSize = [600,300];

 (m,Minv,mref,boundsHigh,boundsLow) = readModelAndGenerateMeshMref(modelDir,
 	"examples/SEGmodel2Dsalt.dat",dim,pad,[0.0,13.5,0.0,4.2],newSize,1.752,2.9);

omega = [3.0,3.3,3.6,3.9,4.2,4.5,5.0,5.5,6.5]*2*pi;
offset  = newSize[1];
println("Offset is: ",offset," cells.")
alpha1 = 5e0;
alpha2 = 5e1;
stepReg = 5e1;
EScycles = 2;
cgit = 7;

freqContSweeps = 5;
freqRanges = [(1,4), (1,4), (4,length(omega)), (4,length(omega)),
		(length(omega), length(omega))];
regularizations = ["high", "high", "low", "low", "low"];
GNiters = [50, 50, 15 ,15, 100];

# ###################################################################################################################
dataFilenamePrefix = string(dataDir,"/DATA_",tuple((Minv.n)...));
resultsFilename = string(resultsDir,"/FWI_",tuple((Minv.n)...));
#######################################################################################################################
writedlm(string(resultsFilename,"_mtrue.dat"),convert(Array{Float16},m));
writedlm(string(resultsFilename,"_mref.dat"),convert(Array{Float16},mref));
resultsFilename = string(resultsFilename,".dat");

println("omega*maximum(h): ",omega*maximum(Minv.h)*sqrt(maximum(1.0./(boundsLow.^2))));
ABLpad = pad + 4;
Ainv  = getParallelJuliaSolver(ComplexF64,Int64,numCores=16,backend=3);

workersFWI = workers();
println(string("The workers that we allocate for FWI are:",workersFWI));

figure(1,figsize = (22,10));
plotModel(m,includeMeshInfo=true,M_regular = Minv,cutPad=pad,limits=[1.5,4.5],figTitle="mref",filename="orig.png");

figure(2,figsize = (22,10));
plotModel(mref,includeMeshInfo=true,M_regular = Minv,cutPad=pad,limits=[1.5,4.5],figTitle="mref",filename="mref.png");

prepareFWIDataFiles(m,Minv,mref,boundsHigh,boundsLow,dataFilenamePrefix,omega,ones(ComplexF64,size(omega)),
									pad,ABLpad,jumpSrc,jumpRcv,offset,workersFWI,maxBatchSize,Ainv,useFilesForFields);

(Q,P,pMis,SourcesSubInd,contDiv,Iact,sback,mref,boundsHigh,boundsLow) =
	setupFWI(m,dataFilenamePrefix,plotting,workersFWI,maxBatchSize,Ainv,SSDFun,useFilesForFields, true);

########################################################################################################
# Setting up the inversion for slowness instead of velocity:
########################################################################################################
function dump(mc,Dc,iter,pInv,PMis,resultsFilename)
	if iter==0
		return;
	end
	fullMc = slowSquaredToVelocity(reshape(Iact*pInv.modelfun(mc)[1] + sback,tuple((pInv.MInv.n)...)))[1];
	Temp = splitext(resultsFilename);
	if iter>0
		Temp = string(Temp[1],iter,Temp[2]);
	else
		Temp = resultsFilename;
	end
	if resultsFilename!=""
		writedlm(Temp,convert(Array{Float16},fullMc));
	end
	if plotting
		figure(888,figsize = (22,10));
		clf();
		filename = splitdir(Temp)[2];
		plotModel(fullMc,includeMeshInfo=true,M_regular = Minv,cutPad=pad,limits=[1.5,4.5],filename=filename,figTitle=filename);
	end
end

#####################################################################################################
# Setting up the inversion for velocity:
#####################################################################################################
mref 		= velocityToSlowSquared(mref)[1];
t    		= copy(boundsLow);
boundsLow 	= velocityToSlowSquared(boundsHigh)[1];
boundsHigh 	= velocityToSlowSquared(t)[1]; t = 0;
modfun 		= identityMod;

########################################################################################################
# Set up Inversion #################################################################################
########################################################################################################

flush(Base.stdout)

GN = "projGN"
maxStep=0.05*maximum(boundsHigh);
regparams = [1.0,1.0,1.0,1e-6];
regfunLow(m,mref,M) 	= wdiffusionReg(m,mref,M,Iact=Iact,C=[]);
regfunHigh(m,mref,M) 	= wFourthOrderSmoothing(m,mref,M,Iact=Iact,C=[]);
if dim==2
	HesPrec=getExactSolveRegularizationPreconditioner();
else
	HesPrec = getSSORCGFourthOrderRegularizationPreconditioner(regparams,Minv,Iact,1.0,1e-8,1000);
end

alpha 	= 1e+2;
pcgTol 	= 1e-1;
maxit 	= 1;

pInv = getInverseParam(Minv,modfun,regfunHigh,alpha,mref[:],boundsLow,boundsHigh,
                         maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
						 minUpdate=1e-3, maxIter = maxit,HesPrec=HesPrec);
mc = copy(mref[:]);

N_nodes = prod(Minv.n.+1);
nsrc = size(Q,2);
p = 16;
Z1 = 2e-4*rand(ComplexF64,(N_nodes, p));

function saveCheckpoint(resultsFilename,mc,Z1,Z2,alpha1,alpha2,pInv,cyc)
	file = matopen(string(splitext(resultsFilename)[1],"_Cyc",cyc,"_checkpoint.mat"), "w");
	write(file,"mc",mc);
	write(file,"Z1",Z1);
	write(file,"Z2",Z2);
	write(file,"alpha1",alpha1);
	write(file,"alpha2",alpha2);
	write(file,"alpha",pInv.alpha);
	write(file,"mref",pInv.mref);
	close(file);
	println("****************************************************************************")
	println("*********************** Saving Checkpoint for cycle ",cyc," ********************")
	println("****************************************************************************")
end

function loadCheckpoint(resultsFilename,cyc)
	file = matopen(string(splitext(resultsFilename)[1],"_Cyc",cyc,"_checkpoint.mat"), "r");
	mc = read(file,"mc");
	Z1 = read(file,"Z1");
	Z2 = read(file,"Z2");
	alpha1 = read(file,"alpha1");
	alpha2 = read(file,"alpha2");
	alpha = read(file,"alpha");
	mref = read(file,"mref");
	close(file);
	return mc,Z1,Z2,alpha1,alpha2,alpha,mref
end

if norm(Z1) == 0.0
	# Standard FWI run
	freqContParams = getFreqContParams(mc, 0, Q,size(P,2), pInv, pMis,
			windowSize, resultsFilename,dump,Iact,sback,
			simSrcDim = simSrcDim);
else
	freqContParams = getFreqContParams(mc, 0, Q,size(P,2), pInv, pMis,
			windowSize, resultsFilename,dump,Iact,sback, Z1=Z1, alpha1=alpha1,
			alpha2Orig=alpha2, stepReg=stepReg,
			simSrcDim = simSrcDim, FWImethod="FWI_ES");
end

for i = 1:freqContSweeps
	freqContParams.cycle = i - 1;
	freqContParams.itersNum = GNiters[i];
	freqContParams.startFrom = freqRanges[i][1];
	freqContParams.endAt = freqRanges[i][2];
	if i > EScycles
		freqContParams.FWImethod = "FWI";
	end
	if regularizations[i] == "low"
		freqContParams.pInv.regularizer = regfunLow;
		freqContParams.updateMref = true;
		freqContParams.pInv.pcgMaxIter = 5;
	else
		freqContParams.pInv.regularizer = regfunHigh;
		freqContParams.updateMref = false;
		freqContParams.pInv.pcgMaxIter = 7;
	end
	mc, = freqCont(freqContParams);
	freqContParams.mc = mc;
end