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

struct SnapshotsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a snapshot.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: snapshotName::String (required)
Param: api_version::String (required)
Param: snapshot::Snapshot (required)
Return: Snapshot
"""
function _swaggerinternal_snapshotsCreateOrUpdate(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String, snapshot; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", Snapshot, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/snapshots/{snapshotName}", ["azure_auth"], snapshot)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "snapshotName", snapshotName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function snapshotsCreateOrUpdate(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String, snapshot; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsCreateOrUpdate(_api, subscriptionId, resourceGroupName, snapshotName, api_version, snapshot; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function snapshotsCreateOrUpdate(_api::SnapshotsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String, snapshot; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsCreateOrUpdate(_api, subscriptionId, resourceGroupName, snapshotName, api_version, snapshot; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes a snapshot.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: snapshotName::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function _swaggerinternal_snapshotsDelete(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/snapshots/{snapshotName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "snapshotName", snapshotName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function snapshotsDelete(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsDelete(_api, subscriptionId, resourceGroupName, snapshotName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function snapshotsDelete(_api::SnapshotsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsDelete(_api, subscriptionId, resourceGroupName, snapshotName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets information about a snapshot.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: snapshotName::String (required)
Param: api_version::String (required)
Return: Snapshot
"""
function _swaggerinternal_snapshotsGet(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Snapshot, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/snapshots/{snapshotName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "snapshotName", snapshotName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function snapshotsGet(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsGet(_api, subscriptionId, resourceGroupName, snapshotName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function snapshotsGet(_api::SnapshotsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsGet(_api, subscriptionId, resourceGroupName, snapshotName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Grants access to a snapshot.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: snapshotName::String (required)
Param: api_version::String (required)
Param: grantAccessData::GrantAccessData (required)
Return: AccessUri
"""
function _swaggerinternal_snapshotsGrantAccess(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String, grantAccessData; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", AccessUri, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/snapshots/{snapshotName}/beginGetAccess", ["azure_auth"], grantAccessData)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "snapshotName", snapshotName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function snapshotsGrantAccess(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String, grantAccessData; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsGrantAccess(_api, subscriptionId, resourceGroupName, snapshotName, api_version, grantAccessData; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function snapshotsGrantAccess(_api::SnapshotsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String, grantAccessData; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsGrantAccess(_api, subscriptionId, resourceGroupName, snapshotName, api_version, grantAccessData; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists snapshots under a subscription.
Param: subscriptionId::String (required)
Param: api_version::String (required)
Return: SnapshotList
"""
function _swaggerinternal_snapshotsList(_api::SnapshotsApi, subscriptionId::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", SnapshotList, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/snapshots", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function snapshotsList(_api::SnapshotsApi, subscriptionId::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsList(_api, subscriptionId, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function snapshotsList(_api::SnapshotsApi, response_stream::Channel, subscriptionId::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsList(_api, subscriptionId, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists snapshots under a resource group.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Return: SnapshotList
"""
function _swaggerinternal_snapshotsListByResourceGroup(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", SnapshotList, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/snapshots", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function snapshotsListByResourceGroup(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsListByResourceGroup(_api, subscriptionId, resourceGroupName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function snapshotsListByResourceGroup(_api::SnapshotsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsListByResourceGroup(_api, subscriptionId, resourceGroupName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Revokes access to a snapshot.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: snapshotName::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function _swaggerinternal_snapshotsRevokeAccess(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/snapshots/{snapshotName}/endGetAccess", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "snapshotName", snapshotName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function snapshotsRevokeAccess(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsRevokeAccess(_api, subscriptionId, resourceGroupName, snapshotName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function snapshotsRevokeAccess(_api::SnapshotsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsRevokeAccess(_api, subscriptionId, resourceGroupName, snapshotName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Updates (patches) a snapshot.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: snapshotName::String (required)
Param: api_version::String (required)
Param: snapshot::SnapshotUpdate (required)
Return: Snapshot
"""
function _swaggerinternal_snapshotsUpdate(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String, snapshot; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", Snapshot, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/snapshots/{snapshotName}", ["azure_auth"], snapshot)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "snapshotName", snapshotName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function snapshotsUpdate(_api::SnapshotsApi, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String, snapshot; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsUpdate(_api, subscriptionId, resourceGroupName, snapshotName, api_version, snapshot; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function snapshotsUpdate(_api::SnapshotsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, snapshotName::String, api_version::String, snapshot; _mediaType=nothing)
    _ctx = _swaggerinternal_snapshotsUpdate(_api, subscriptionId, resourceGroupName, snapshotName, api_version, snapshot; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export snapshotsCreateOrUpdate, snapshotsDelete, snapshotsGet, snapshotsGrantAccess, snapshotsList, snapshotsListByResourceGroup, snapshotsRevokeAccess, snapshotsUpdate


text:
import StatsBase: predict

import DataFrames: AbstractDataFrame, DataFrame

export predict

predict(jlt::JLBoostTreeModel, df) = predict(trees(jlt), df)

predict(jlt::JLBoostTreeModel, df::AbstractDataFrame) = predict(trees(jlt), df)

function predict(jlt::AbstractJLBoostTree, df)
	# TODO a more efficient algorithm. Currently there are too many assignbools being
	# stores the results
	res = Vector{Float64}(undef, nrow(df))
	res .= 0.0

	# stores the assignment array
	assignbool = trues(nrow(df))

    predict!(jlt, df, res, assignbool)
end

function predict(jlts, df)
	mapreduce(x->predict(x, df), +, jlts)
end

function predict!(jlt::JLBoostTree, df, res, assignbool)
	if length(jlt.children) == 2
	    new_assignbool = assignbool .& (getproperty(Tables.columns(df), jlt.splitfeature) .<= jlt.split)
	    predict!(jlt.children[1], df, res, new_assignbool)

	    new_assignbool .= assignbool .& (getproperty(Tables.columns(df), jlt.splitfeature) .> jlt.split)
	    predict!(jlt.children[2], df, res, new_assignbool)
	else length(jlt.children) == 0
	    res[assignbool] .+= jlt.weight
	end
	res
end

function predict!(jlt::WeightedJLBoostTree, df, res, assignbool)
	jlt.eta .* predict!(jlt.tree, df, res, assignbool)
end

ext: jl


text:
using Mimi

include("/Users/lisarennels/.julia/dev/Mimi/wip/create_composite.jl")

import Mimi: components, ComponentPath, ComponentInstance, LeafComponentInstance,
    CompositeComponentInstance, find_comp, dim_names

## Step 1. COMPONENT MENU LIST: Get the Component Structure of the Model
# Create a lightweight structure holding the ComponentPaths of all components,
# which we can then use to look things up on-click.

# Note that this traversal could in essence pull out everything we need at once,
# but this file will break it out into the steps we need to (1) create the hierarchy
# of components and their props (2) on click find all variable and parameter names and
# (3) on that click get the dimensions and values used for VegaLite

# primary function to return a dictionary with keys being the component names and
# values being the ComponentPaths
function _get_all_paths(m::Model)
    all_paths = Dict{Symbol, ComponentPath}()
    for comp in components(m) # iterate over top level ComponentInstances
        _add_paths(m, comp, all_paths)
    end
    return all_paths
end

# a helper function to perform a preorder traversal of a given top-level component
# in model m and add that path, and all sub-component paths, to the paths array
function _add_paths(m::Model, comp::Union{CompositeComponentInstance, LeafComponentInstance}, paths::Dict{Symbol, ComponentPath})
    if isa(comp, CompositeComponentInstance)
        paths[comp.comp_name] = comp.comp_path
        for subcomp in values(comp.comps_dict)
            _add_paths(m, subcomp, paths)
        end
    else # LeafComponentInstance
        paths[comp.comp_name] = comp.comp_path
    end
    return paths
end

## Step 2. ON CLICK OF A COMPONENT DISPLAY VARIABLE AND PARAMETER NAMES
# For any given component, using the information above, which is just the
# component name and path, grab the names of bariables and parameters
# bubbled into that component

paths = _get_all_paths(m)
comp_name = :Comp1
comp_path = paths[comp_name]
comp_def = find_comp(m, comp_path) # returns a ComponentDefinition
vars_names = variable_names(comp_def)
par_names = parameter_names(comp_def)

## Step 3. For each Variable and Parameter get the data and dim names

datum_name = vars_names[1]
dims = dim_names(comp_def, datum_name) # only needs the ComponentDefinition
values = m[comp_path, datum_name] # uses the comp_path to access the ComponentInstance


text:
# solve the proximal problem
# matout = argmin||matout-matin||_F^2, s.t. matout(:,i)<=1
function normcol_lessequal(matin::Matrix{Float64})
    broadcast(/, matin, max(1.0, sqrt(sum(abs2(matin), 1))))
end


text:
struct Output
    seq::String
    seqins::String
    start::Int
    finish::Int
    en::Float64
    sat::Bool
end

mutable struct Alg
    verbose::Bool #verbosity of output
    maxiter::Int  # maximal number of sweeps
    epsconv::Float64 # tolerance for the convergence of T > 0 alg
    mindec::Int # convergence for T = 0 alg, assignment must be consecutively repeated for at least mindec times
    damp::Float64  # damping factor: Pnew = damp*Pold + (1-damp)*Pnew
    μext::Float64 # gap external penalty    exp[-μext], nᵢ ∈ {N+1, 0}, xᵢ = 0
    λo::Vector{Float64} # penalty for opening exp[-λo]I(Δn > 0)
    λe::Vector{Float64} # penalty for extending exp[-λe(Δn - 1)I(Δn >0)]
    μint::Float64 #gap internal penalty exp[-μint], N ≥ n ≥ 1, xᵢ = 0
    nprint::Int # print info every nprint iterations
end

struct Seq
    header::String
    strseq::String
    intseq::Vector{Int}
    ctype::Symbol
    function Seq(h::String,s::String,v::Vector{Int},ctype::Symbol)
        if length(s) != length(v)
            error("length s != length v")
        end
        for i in eachindex(s) #check alignment between strseq and inteseq
            if letter2num(s[i],ctype) != v[i]
                error("intseq not representing strseq at index $i")
            end
        end
        return new(h::String,s::String,v::Vector{Int},ctype::Symbol)
    end

end

function Seq(h::String,s::String,ctype::Symbol)
    v = fill(-1,length(s))
    s = Vector{Char}(s)
    if ctype == :nbase
        replace!(s, 'T' => 'U')
    end
    for i in eachindex(s)
        v[i] = letter2num(s[i],ctype)
    end
    s = String(s)
    Seq(h,s,v,ctype)
end

function Base.show(io::IO, x::Seq)
    println(io,x.header)
    #print_with_color(:cyan,io,x.strseq)
    println(io,x.strseq)
end

struct Jh
    J::Array{Float64,4}  # couplings q × q × L × L (L length of Potts Model)
    h::Array{Float64,2}  # field     q × L
    function Jh(J::Array{Float64,4}, h::Array{Float64,2})
        qJ1,qJ2,NJ1,NJ2 = size(J)
        qh,Nh = size(h)
        qJ1 != qJ2 && error("not square color submatrix")
        NJ1 != NJ2 && error("not square sites submatrix")
        NJ1 != Nh  && error("inconsistent number of sites between fields and couplings")
        qJ1 != qh  && error("inconsistent number of colors between fields and couplings")
        return new(J, h)
    end
end

const Marginal = OffsetArrays.OffsetArray{Float64, 2, Array{Float64,2}}
Marginal(N::Int) = OffsetArray(rand(2,N+2), 0:1, 0:N+1)

# We adopt the following convention:
# xᵢ = 1 → match,
# xᵢ = 0 → gap,
# nᵢ = 0 if xᵢ = 0 and no match has been yet found
# N + 1 > nᵢ > 0, xᵢ = 1 for all the matched cases
# N + 1 > nᵢ > 0, xᵢ = 0 for all the internal gaps
# nᵢ = N + 1 if xᵢ = 0 for last gaps

struct PBF
    L::Int # Potts-model length
    N::Int # Sequence-to-align length
    q::Int # Alphabet (20 amino acids + "-") gap = 21.
    P::Vector{Marginal} # P[i] → p[nᵢ, xᵢ] where nᵢ∈ [0,N+1] , xᵢ ∈ [0,1] ∀ i ∈ [1..L]
    B::Vector{Marginal}
    F::Vector{Marginal}
    newP::Marginal
    newB::Marginal
    newF::Marginal
    forward::Marginal
    central::Marginal
    bckward::Marginal
end

function PBF(L::Int, N::Int, q::Int)

    fwds = [Marginal(N) for i in 1:L]
    bcks = [Marginal(N) for i in 1:L]
    PBF(L,N,q,
        [Marginal(N) for i in 1:L],
        bcks,
        fwds,
        Marginal(N),
        Marginal(N),
        Marginal(N),
        Marginal(N),
        Marginal(N),
        Marginal(N))
end

function PBF(x::Jh,seq::String)
    q,L = size(x.h)
    N = length(seq)
    return PBF(L,N,q)
end

function PBF(x::Jh, seq::Seq)
    q,L = size(x.h)
    N = length(seq.intseq)
    return PBF(L,N,q)
end

struct AllVar
    pbf::PBF # marginal
    jh::Jh   # coupling
    seq::Seq # sequence to align
    alg::Alg # algorithm parameters.
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module ParaxialAnalysis

using StaticArrays: promote_tuple_eltype
using ..OpticSim
using ..OpticSim.Geometry
using ..OpticSim.Repeat
using StaticArrays
using LinearAlgebra
import LazySets
using Unitful
using Unitful.DefaultSymbols

include("LensletAssembly.jl")
include("HeadEyeModel/HeadEyeModel.jl")

end #module
export ParaxialAnalysis


text:
function Base.download(client::ZoteroClient, doc::Document; dl=false)
    request_json(client, "GET", joinpath("items", doc.key))
end

function Base.download(client::ZoteroClient, col::Collection; dl=false)
    request_json(client, "GET", joinpath("collections", col.key))
end

function download_item(client::ZoteroClient, key::String, dl=false, kwargs...)
    download(client, ParentDoc(; key=key); dl=dl, kwargs...)
end

function download_collection(client::ZoteroClient, key::String, dl=false, kwargs...)
    download(client, ParentDoc(;key=key); dl=dl, kwargs...)
end

function collect_items!(client::ZoteroClient, col::Collection; refresh=false, kwargs...)
    if isempty(col.docs) || refresh
        dicts = request_json(client, "GET", joinpath("collections", col.key, "items"); kwargs)
        docs = dict_to_doc.(dicts)
        org_docs = organize(docs)
        col.docs = org_docs
    end
    return col.docs
end


text:
# Seedlink with command-line stations
config_file = path*"/SampleFiles/seedlink.conf"
sta = ["CC.SEP", "UW.HDW"]
pat = ["?????.D"; "?????.D"]
trl = ".??.???.D"
sta_matrix = String.(reshape(split(sta[2],'.'), 1,2))

function SL_wait(ta::Array{Union{Task,Nothing}, 1}, t_interval::Int64)
  t = 0
  while t < 300
    if any([!istaskdone(t) for t in ta])
      sleep(t_interval)
      t += t_interval
    elseif t ≥ 60
      println("      one or more queries incomplete after 60 s; skipping test.")
      for i = 1:4
        ta[i] = Nothing
        GC.gc()
      end
      break
    else
      tf1 = fetch(ta[1])
      tf2 = fetch(ta[2])
      tf3 = fetch(ta[3])
      tf4 = fetch(ta[4])
      @test tf1 == tf2 == tf3 == tf4
      break
    end
  end
  return nothing
end

printstyled("  SeedLink\n", color=:light_green, bold=true)
printstyled("  (SeedLink tests require 4-6 minutes)\n", color=:green)

# has_stream
printstyled("    has_stream\n", color=:light_green)
ta = Array{Union{Task,Nothing}, 1}(undef, 4)
ta[1] = @async has_stream(sta, u="rtserve.iris.washington.edu")[2]
ta[2] = @async has_stream(sta, pat, u="rtserve.iris.washington.edu", d='.')[2]
ta[3] = @async has_stream(join(sta, ','), u="rtserve.iris.washington.edu")[2]
ta[4] = @async has_stream(sta_matrix, u="rtserve.iris.washington.edu")[1]
SL_wait(ta, 1)

# has_stream
printstyled("    has_sta\n", color=:light_green)
ta[1] = @async has_sta(sta[1], u="rtserve.iris.washington.edu")[1]
ta[2] = @async has_sta(sta[1]*trl, u="rtserve.iris.washington.edu")[1]
ta[3] = @async has_sta(sta, u="rtserve.iris.washington.edu")[1]
ta[4] = @async has_sta(parse_charr(sta), u="rtserve.iris.washington.edu")[1]
SL_wait(ta, 1)

# Attempting to produce errors
printstyled("    produce expected errors and warnings\n", color=:light_green)
redirect_stdout(out) do
  S1 = SeisData()

  @test_throws ErrorException seedlink!(S1, [sta[1]], ["*****.X"])

  S2 = seedlink([sta[1]], pat, x_on_err=false)
  write(S2.c[1], "BYE\r")
  close(S2.c[1])
  @test_throws ErrorException seedlink!(S2, [replace(sta[1], "SEP" => "XOX")], ["?????.D"])

  S3 = seedlink([replace(sta[1], "SEP" => "XOX")], ["*****.X"], x_on_err=false)
  write(S3.c[1], "BYE\r")
  close(S3.c[1])
end

# DATA mode
printstyled("    DATA mode\n", color=:light_green)
printstyled("      link 1: command-line station list\n", color=:light_green)
T = SeisData()
redirect_stdout(out) do
  seedlink!(T, sta, mode="DATA", refresh=9.9, kai=7.0, v=1)
end

printstyled("      link 2: station file\n", color=:light_green)
redirect_stdout(out) do
  seedlink!(T, config_file, mode="DATA", refresh=13.3, v=3)
end
wait_on_data!(T, tmax=50.0)

# FETCH mode (indistinguishable from DATA mode for most users)
printstyled("    FETCH mode\n", color=:light_green)
V = seedlink("GE.ISP..BH?.D", refresh=10.0, mode="FETCH", v=1)
printstyled("      link initialized\n", color=:light_green)
wait_on_data!(V, tmax=50.0)

# SeedLink time mode (more complicated)
printstyled("    TIME mode\n", color=:light_green)

# To ensure precise timing, we'll pass d0 and d1 as strings
st = 0.0
en = 60.0
dt = en-st
(d0,d1) = parsetimewin(st,en)

U = SeisData()
seedlink!(U, sta, mode="TIME", refresh=10.0, s=d0, t=d1, w=true)
printstyled("      first link initialized\n", color=:light_green)

# Seedlink with a config file
seedlink!(U, config_file, refresh=10.0, mode="TIME", s=d0, t=d1)
printstyled("      second link initialized\n", color=:light_green)

# Seedlink with a config string
redirect_stdout(out) do
  seedlink!(U, "CC.VALT..???, UW.ELK..EHZ", mode="TIME", refresh=10.0, kai=19.0, s=d0, t=d1, v=3)
end
printstyled("      third link initialized\n", color=:light_green)
wait_on_data!(U, tmax=50.0)


text:
@testset "core" begin

  # Accessor functions
  img = fill(1.0, (4,4))
  img[1:2,1:2] .= 2.0
  img[1:2,3:4] .= 3.0
  seg = fast_scanning(img, 0.5)

  @test labels_map(seg) == seg.image_indexmap
  @test segment_labels(seg) == seg.segment_labels
  @test segment_mean(seg) == seg.segment_means
  @test all([segment_mean(seg, i) == seg.segment_means[i] for i in keys(seg.segment_means)])
  @test segment_pixel_count(seg) == seg.segment_pixel_count
  @test all([segment_pixel_count(seg, i) == seg.segment_pixel_count[i] for i in keys(seg.segment_pixel_count)])

  img = [1 1.1 0.9 10 10.2 9.8; 1 1 1.1 4 3.8 4]
  r = fuzzy_cmeans(img, 2, 2.0)

  @test segment_labels(r) == [1,2]
  @test segment_mean(r, 1) == r.centers[:,1]
  @test segment_mean(r) == Dict(1=>r.centers[:,1], 2=>r.centers[:,2])
  @test all([segment_pixel_count(r,i) == 6 for i in 1:2])
  @test segment_pixel_count(r) == Dict(1=>6, 2=>6)

  # RAG
  img = fill(1.0, (10,10))
  img[:, 5:10] .= 2.0
  seg = fast_scanning(img, 0.5)
  g, vm = region_adjacency_graph(seg, (i,j)->sum(abs2, seg.segment_means[i]-seg.segment_means[j]))

  expectedg = SimpleWeightedGraph(2)
  add_edge!(expectedg, 1, 2, 1.0)
  expectedvm = Dict(1=>1, 2=>2)

  @test g == expectedg
  @test vm == expectedvm

  img = fill(1.0, (10,10))
  img[3:8,3:8] .= 2.0
  img[5:7,5:7] .= 3.0
  seg = fast_scanning(img, 0.5)
  g, vm = region_adjacency_graph(seg, (i,j)->1)

  expectedg = SimpleWeightedGraph(3)
  add_edge!(expectedg,1,2)
  add_edge!(expectedg,2,3)
  expectedvm = Dict(1=>1, 2=>2, 3=>3)

  @test g == expectedg
  @test vm == expectedvm

  # remove_segment
  img = fill(1.0, (10,10))
  img[1:4,:] .= 2.0
  img[:,5:10] .= 4.0
  seg = fast_scanning(img, 0.5)
  new_seg = remove_segment(seg, 1, (i,j)->(-seg.segment_pixel_count[j]))

  expected = fill(3, (10,10))
  expected[5:10,1:4] .= 2
  expected_labels = [2,3]
  expected_means = Dict(2=>1.0, 3=>68/19)
  expected_count = Dict(2=>24, 3=>76)

  @test all(label->(label in expected_labels), new_seg.segment_labels)
  @test all(label->(label in new_seg.segment_labels), expected_labels)
  @test expected_count == new_seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ new_seg.segment_means[label]), new_seg.segment_labels)
  @test new_seg.image_indexmap == expected

  new_seg = remove_segment(seg, 1, (i,j)->sum(abs2, seg.segment_means[i]-seg.segment_means[j]))

  expected = fill(3, (10,10))
  expected[:,1:4] .= 2
  expected_means = Dict(2=>1.4, 3=>4.0)
  expected_count = Dict(2=>40, 3=>60)

  @test all(label->(label in expected_labels), new_seg.segment_labels)
  @test all(label->(label in new_seg.segment_labels), expected_labels)
  @test expected_count == new_seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ new_seg.segment_means[label]), new_seg.segment_labels)
  @test new_seg.image_indexmap == expected

  remove_segment!(seg, 1, (i,j)->sum(abs2, seg.segment_means[i]-seg.segment_means[j]))

  @test all(label->(label in expected_labels), seg.segment_labels)
  @test all(label->(label in seg.segment_labels), expected_labels)
  @test expected_count == seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ seg.segment_means[label]), seg.segment_labels)
  @test seg.image_indexmap == expected

  # prune_segments
  img = fill(1.0, (10,10))
  img[3,3] = 2.0
  img[5,5] = 2.0
  img[8:9,7:8] .= 3.0
  seg = fast_scanning(img, 0.5)
  new_seg = prune_segments(seg, l->(seg.segment_pixel_count[l] < 2), (i,j)->sum(abs2, seg.segment_means[i]-seg.segment_means[j]))

  expected = fill(1, (10,10))
  expected[8:9,7:8] .= 4
  expected_labels = [1,4]
  expected_means = Dict(4=>3.0, 1=>49/48)
  expected_count = Dict(4=>4, 1=>96)

  @test all(label->(label in expected_labels), new_seg.segment_labels)
  @test all(label->(label in new_seg.segment_labels), expected_labels)
  @test expected_count == new_seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ new_seg.segment_means[label]), new_seg.segment_labels)
  @test new_seg.image_indexmap == expected

  new_seg = prune_segments(seg, [2,3], (i,j)->sum(abs2, seg.segment_means[i]-seg.segment_means[j]))

  @test all(label->(label in expected_labels), new_seg.segment_labels)
  @test all(label->(label in new_seg.segment_labels), expected_labels)
  @test expected_count == new_seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ new_seg.segment_means[label]), new_seg.segment_labels)
  @test new_seg.image_indexmap == expected

  # region_adjacency_graph
  img = fill(1.0, (10,10))
  img[4:6, 4:6] .= 0

  weight_fn(I, J) = exp(-abs(I.second - J.second)) / ((I.first[1]-J.first[1])^2 + (I.first[2]-J.first[2])^2)
  G, vertex2cartesian = region_adjacency_graph(img, weight_fn, 2)

  @test nv(G) == 100
  @test G.weights[45,46] == 1
  @test G.weights[43,44] == exp(-1)

  #test that all edges weights equal weights computed using weight_fn
  @test all(map(edge-> edge.weight == weight_fn(vertex2cartesian[edge.src]=>img[vertex2cartesian[edge.src]], vertex2cartesian[edge.dst]=>img[vertex2cartesian[edge.dst]]), edges(G)))

  #test that all edges are between nodes less than r pixels away
  @test all(map(edge-> max(abs(vertex2cartesian[edge.src][1] - vertex2cartesian[edge.dst][1]), abs(vertex2cartesian[edge.src][2] - vertex2cartesian[edge.dst][2])) <= 2 , edges(G)))

  @test region_adjacency_graph(img, weight_fn, 2) == region_adjacency_graph(img, weight_fn, CartesianIndex(2,2))

  G, vertex2cartesian = region_adjacency_graph(img, weight_fn, CartesianIndex(1,2))
  @test G.weights[45,47] == 0
  @test G.weights[45,65] != 0

end


text:
## Metsa-Oy Forest Products Supply Chain
## Edward J. Xu <edxu96@outlook.com>
## Jan 24th, 2020

function optim_3(alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm)

  println("Optimize model 3.")
  mod = Model(with_optimizer(Gurobi.Optimizer, Presolve=0, OutputFlag=0))

  ## Here-and-now decision variables in the first year.
  @variable(mod, h_t[1:6]>=0, integer=true)
  @variable(mod, y_i[1:5]>=0)
  @variable(mod, y_j[1:2]>=0)
  @variable(mod, y_paper>=0)
  @variable(mod, z_ik[1:5, 1:4]>=0)
  @variable(mod, z_jk[1:2, 1:4]>=0)
  @variable(mod, z_paper_k[1:4]>=0)

  ## Capacities for the second year are here-and-now decision variables.
  @variable(mod, 0 <= x_saw_nm[1:4, 2:3] <= r_saw * nu_saw)
  @variable(mod, 0 <= x_plywood_nm[1:4, 2:3] <= r_plywood * nu_plywood)
  @variable(mod, x_nmj[1:4, 2:3, 1:2] >= 0)
  @constraint(mod, [n=1:4, m=2:3, j=1:2], x_nmj[n, m, j] <= r_j[j] * nu_j[j])
  @variable(mod, 0 <= x_paper_nm[1:4, 2:3] <= r_paper * nu_paper)

  @constraint(mod, [n=2:4], x_saw_nm[1, 2] == x_saw_nm[n, 2])
  @constraint(mod, [n=2:4], x_plywood_nm[1, 2] == x_plywood_nm[n, 2])
  @constraint(mod, [n=2:4, j=1:2], x_nmj[1, 2, j] == x_nmj[n, 2, j])
  @constraint(mod, [n=2:4], x_paper_nm[1, 2] == x_paper_nm[n, 2])

  ## Wait-and-see decisions variables in the second and third year.
  @variable(mod, h_nmt[1:4, 2:3, 1:6]>=0, integer=true)
  @variable(mod, y_nmi[1:4, 2:3, 1:5]>=0)
  @variable(mod, y_nmj[1:4, 2:3, 1:2]>=0)
  @variable(mod, y_paper_nm[1:4, 2:3]>=0)
  @variable(mod, z_nmik[1:4, 2:3, 1:5, 1:4]>=0)
  @variable(mod, z_nmjk[1:4, 2:3, 1:2, 1:4]>=0)
  @variable(mod, z_paper_nmk[1:4, 2:3, 1:4]>=0)

  @objective(mod, Max,
    ## Objective fucntions regarding here-and-now decision variables.
    -sum(alpha_t[t] * h_t[t] + beta_t[t] * h_t[t]^2 for t = 1:6) +
    -sum(c_i[i] * y_i[i] for i=1:5) +
    -sum(c_j[j] * y_j[j] for j=1:2) - c_paper * y_paper +
    sum(alpha_t[t] * (
      h_t[t] - sum(a_i[i] * y_i[i] for i in i_t1[t])
      ) for t=1:3) +
    sum(alpha_t[t] * (
      h_t[t] + sum(b_i[i] * y_i[i] for i in i_t2[t-3]) -
      sum(a_j[j] * y_j[j] for j in j_t2[t-3]) -
        a_paper_t[t-3] * y_paper) for t=4:6) +
    sum(e_i[i] * y_i[i] for i=1:5) * p_fuel +
    sum(gamma_ik[i, k] * z_ik[i, k] - delta_ik[i, k] * z_ik[i, k]^2
      for i=1:5, k=1:4) +
    sum(gamma_jk[j, k] * z_jk[j, k] - delta_jk[j, k] * z_jk[j, k]^2
      for j=1:2, k=1:4) +
    sum(gamma_paper_k[k] * z_paper_k[k] -
      delta_paper_k[k] * z_paper_k[k]^2
      for k=1:4) +

    ## here-and-now capacity in the second year.
    ## those in the first scenario are used because they are the same.
    - (o_saw * (x_saw_nm[1, 2] - r_saw) +
      o_plywood * (x_plywood_nm[1, 2] - r_plywood) +
      sum(o_j[j] * (x_nmj[1, 2, j] - r_j[j]) for j=1:2) +
      o_paper * (x_paper_nm[1, 2] - r_paper)) +

    ## Objective fucntions regarding wait-and-see decision variables.
    sum(
      pi_n[n] * (
        sum(
          sigma^(m-1) * (
            ## 1, cost of timber procurement
            -sum(alpha_t[t] * h_nmt[n, m, t] + beta_t[t] * h_nmt[n, m, t]^2
              for t = 1:6) +
            ## 2, cost of wood production
            -sum(c_i[i] * y_nmi[n, m, i] for i=1:5) +
            ## 3, cost of pulp and paper production
            -sum(c_j[j] * y_nmj[n, m, j] for j=1:2) +
            -c_paper * y_paper_nm[n, m] +
            ## 4, profit of left timbers selling
            sum(alpha_t[t] * (h_nmt[n, m, t] - sum(a_i[i] * y_nmi[n, m, i]
              for i in i_t1[t])) for t=1:3) +
            sum(alpha_t[t] * (h_nmt[n, m, t] + sum(b_i[i] * y_nmi[n, m, i]
              for i in i_t2[t-3]) - sum(a_j[j] * y_nmj[n, m, j] for j in
              j_t2[t-3]) - a_paper_t[t-3] * y_paper_nm[n, m]) for t=4:6) +
            ## 5, profit of fuel wood selling
            p_fuel * sum(e_i[i] * y_nmi[n, m, i] for i=1:5) +
            ## 6, profit of wood selling
            sum(omega_i[i]^(m-1) * rho_nm[n, m] * gamma_ik[i, k] *
              z_nmik[n, m, i, k] - omega_i[i]^(m-1) * rho_nm[n, m] *
              delta_ik[i, k] * z_nmik[n, m, i, k]^2 for i=1:5, k=1:4) +
            ## 7, profit of pulp selling
            sum(omega_j[j]^(m-1) * rho_nm[n, m] * gamma_jk[j, k] *
              z_nmjk[n, m, j, k] - omega_j[j]^(m-1) * rho_nm[n, m] *
              delta_jk[j, k] * z_nmjk[n, m, j, k]^2 for j=1:2, k=1:4) +
            ## 8, profit of paper selling
            sum(omega_paper^(m-1) * rho_nm[n, m] * gamma_paper_k[k] *
              z_paper_nmk[n, m, k] - omega_paper^(m-1) * rho_nm[n, m] *
              delta_paper_k[k] * z_paper_nmk[n, m, k]^2 for k=1:4)
            ) for m = 2:3
          ) +

        ## Cost of capacity expansion in the third year
        - sigma * (
          o_saw * (x_saw_nm[n, 3] - x_saw_nm[1, 2]) +
          o_plywood * (x_plywood_nm[n, 3] - x_plywood_nm[1, 2]) +
          sum(o_j[j] * (x_nmj[n, 3, j] - x_nmj[1, 2, j]) for j=1:2) +
          o_paper * (x_paper_nm[n, 3] - x_paper_nm[n, 2])
          )

        ) for n=1:4
      )
    )

  ## Constraints for here-and-now decision variables
  @constraint(mod, [t=1:3],
    h_t[t] >= sum(a_i[i] * y_i[i] for i in i_t1[t]))
  @constraint(mod, [t=4:6],
    h_t[t] + sum(b_i[i] * y_i[i] for i in i_t2[t-3]) >=
    sum(a_j[j] * y_j[j] for j in j_t2[t-3]) + a_paper_t[t-3] * y_paper)
  @constraint(mod, [j=1:2],
    y_j[j] >= b_paper_j[j] * y_paper)
  @constraint(mod, [i=1:5],
    y_i[i] >= sum(z_ik[i, k] for k=1:4))
  @constraint(mod, [j=1:2],
    y_j[j] - b_paper_j[j] * y_paper >= sum(z_jk[j, k] for k=1:4))
  @constraint(mod, y_paper >= sum(z_paper_k[k] for k=1:4))
  @constraint(mod, sum(y_i[i] for i=1:3) <= r_saw)
  @constraint(mod, sum(y_i[i] for i=4:5) <= r_plywood)
  @constraint(mod, [j=1:2], y_j[j] <= r_j[j])
  @constraint(mod, y_paper <= r_paper)

  ## Constraints for wait-and-see decision variables.
  @constraint(mod, [n=1:4, t=1:3, m=2:3],
    h_nmt[n, m, t] >= sum(a_i[i] * y_nmi[n, m, i] for i in i_t1[t]))
  @constraint(mod, [n=1:4, t=4:6, m=2:3],
    h_nmt[n, m, t] + sum(b_i[i] * y_nmi[n, m, i] for i in i_t2[t-3]) >=
    sum(a_j[j] * y_nmj[n, m, j] for j in j_t2[t-3]) +
    a_paper_t[t-3] * y_paper_nm[n, m])
  @constraint(mod, [n=1:4, j=1:2, m=2:3],
    y_nmj[n, m, j] >= b_paper_j[j] * y_paper_nm[n, m])
  @constraint(mod, [n=1:4, i=1:5, m=2:3],
    y_nmi[n, m, i] >= sum(z_nmik[n, m, i, k] for k=1:4))
  @constraint(mod, [n=1:4, j=1:2, m=2:3], y_nmj[n, m, j] -
    b_paper_j[j] * y_paper_nm[n, m] >= sum(z_nmjk[n, m, j, k] for k=1:4))
  @constraint(mod, [n=1:4, m=2:3],
    y_paper_nm[n, m] >= sum(z_paper_nmk[n, m, k] for k=1:4))
  @constraint(mod, [n=1:4, m=2:3],
    sum(y_nmi[n, m, i] for i=1:3) <= x_saw_nm[n, m])
  @constraint(mod, [n=1:4, m=2:3],
    sum(y_nmi[n, m, i] for i=4:5) <= x_plywood_nm[n, m])
  @constraint(mod, [n=1:4, j=1:2, m=2:3], y_nmj[n, m, j] <= x_nmj[n, m, j])
  @constraint(mod, [n=1:4, m=2:3], y_paper_nm[n, m] <= x_paper_nm[n, m])

  ## relation between capacities
  @constraint(mod, [n=1:4], x_saw_nm[n, 2] <= x_saw_nm[n, 3])
  @constraint(mod, [n=1:4], x_plywood_nm[n, 2] <= x_plywood_nm[n, 3])
  @constraint(mod, [n=1:4, j=1:2], x_nmj[n, 2, j] <= x_nmj[n, 3, j])
  @constraint(mod, [n=1:4], x_paper_nm[n, 2] <= x_paper_nm[n, 3])

  @constraint(mod, r_saw <= x_saw_nm[1, 2])
  @constraint(mod, r_plywood <= x_plywood_nm[1, 2])
  @constraint(mod, [j=1:2], r_j[j] <= x_nmj[1, 2, j])
  @constraint(mod, r_paper <= x_paper_nm[1, 2])

  optimize!(mod)

  ##
  result_h_t = value_vec(h_t)
  result_y_i = value_vec(y_i)
  result_y_j = value_vec(y_j)
  result_y_paper = value(y_paper)
  result_z_ik = value_mat(z_ik)
  result_z_jk = value_mat(z_jk)
  result_z_paper_k = value_vec(z_paper_k)

  result_h_nmt = value_mat3_extra1(h_nmt)
  result_y_nmi = value_mat3_extra1(y_nmi)
  result_y_nmj = value_mat3_extra1(y_nmj)
  result_y_paper_nm = value_mat_extra1(y_paper_nm)
  result_z_nmik = value_mat4_extra1(z_nmik)
  result_z_nmjk = value_mat4_extra1(z_nmjk)
  result_z_paper_nmk = value_mat3_extra1(z_paper_nmk)

  result_x_saw_nm = value_mat_extra1(x_saw_nm)
  result_x_plywood_nm = value_mat_extra1(x_plywood_nm)
  result_x_nmj = value_mat3_extra1(x_nmj)
  result_x_paper_nm = value_mat_extra1(x_paper_nm)

  result_obj = objective_value(mod)
  println(result_obj)

  print_result_3(
    result_h_t, result_y_i, result_y_j, result_y_paper, result_z_ik,
    result_z_jk, result_z_paper_k, result_h_nmt, result_y_nmi, result_y_nmj,
    result_y_paper_nm, result_z_nmik, result_z_nmjk, result_z_paper_nmk,
    alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm, result_x_saw_nm,
    result_x_plywood_nm, result_x_nmj, result_x_paper_nm)

  # return result_h_t, result_y_i, result_y_j, result_y_paper, result_z_ik,
  #   result_z_jk, result_z_paper_k, result_h_nmt, result_y_nmi, result_y_nmj,
  #   result_y_paper_nm, result_z_nmik, result_z_nmjk, result_z_paper_nmk,
  #   alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
  #   c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
  #   gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
  #   omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
  #   o_j, o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm, result_x_saw_nm,
  #   result_x_plywood_nm, result_x_nmj, result_x_paper_nm
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "fts"
version = v"1.2.7"

# Collection of sources required to build fts
sources = [
    ArchiveSource("https://github.com/pullmoll/musl-fts/archive/v$version.zip",
                  "3b1fe92f1d8cb98488d1bce2ad078cd815f10a0fad0c03caf30229c9318f300b"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/musl-fts-*/
./bootstrap.sh
CFLAGS="-fPIC" ./configure --prefix=${prefix} --host=${target}
make -j${nproc}
make install
"""

# Select Unix platforms
platforms = [p for p in supported_platforms() if Sys.islinux(p) && libc(p) == "musl"]

# The products that we will ensure are always built
products = [
    LibraryProduct("libfts", :libfts),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
struct NGramHMM
        grams::Int
        dict::AhoCorasickAutomaton{Ti}
        words::Vector{String}
        usr_words::Int
        tags::Vector{String}
        pmp::Dict{String,Int}
        hpr::Vector{Tv}
        h2h::Matrix{Tv}
        h2v::Vector{Dict{Int, Tv}}
        INF::Vector{Tv}
        usrh2v::Vector{Dict{Int, Tv}}
        usrINF::Vector{Tv}
        usrw::Tv
end

function NGramHMM(grams::Int, corpus, dictpath::String, usrw::Float64)
        @assert grams >= 2
        ispath(dictpath) || @warn "dictpath is not a valid path"
        @assert usrw < 1.

        tags, words = KongYiji.postags(corpus), KongYiji.words(corpus)
        np, pmp = length(tags), KongYiji.dict(tags)
        B = np + 1
        nf = (B ^ (grams - 1)) - 1
        hpr, h2h, h2v, INF = fill(Tv(0), np), fill(Tv(0), (nf, np)), [Dict{Int, Tv}() for _ in 1:np], fill(Tv(0), np)
        usrh2v, usrINF = [Dict{Int, Tv}() for _ in 1:np], fill(Tv(0), np)
        for sent in sents(corpus)
                ns = length(sent)
                for i = 1:ns
                        ip, iw = sent[i]
                        if i == 1; hpr[ip] += 1;
                        else
                                from, pow = 0, 1
                                for j = i - 1:-1:max(1, i - grams + 1); from += sent[j][1] * pow; pow *= B; end
                                h2h[from, ip] += 1
                        end
                        h2v[ip][iw] = get(h2v[ip], iw, 0) + 1
                end
        end

        wmp = KongYiji.dict(words)

        if ispath(dictpath)
                cd(() -> begin
                                #@show pwd()
                                for tag in tags
                                        if !ispath(tag); continue; end
                                        ip = get(pmp, tag, 0)
                                        for file in readdir(tag)
                                                for word in eachline(joinpath(tag, file))
                                                        iw = getid!(wmp, word)
                                                        if ip > 0; usrh2v[ip][iw] = get(usrh2v[ip], iw, 0) + 1; end
                                                end
                                        end
                                end
                        end, dictpath)
        end

        dict = AhoCorasickAutomaton{Ti}(wmp)
        usr_words = length(wmp) - length(words)
        words = KongYiji.vec(wmp)
        normalize!(hpr)
        normalize!(h2h)
        normalize!(h2v, INF)
        log!(hpr)
        log!(h2h)
        normalize!(usrh2v, usrINF)

        NGramHMM(grams, dict, words, usr_words, tags, pmp, hpr, h2h, h2v, INF, usrh2v, usrINF, usrw)
end

NGramHMM() = load(KongYiji.dir("ngramhmm.jld2"))["m"]

function (m::NGramHMM)(xs::Vector{String}; recover=true, withpos=false)
        nc_max = mapreduce(ncodeunits, max, xs)
        nf, np = size(m.h2h)
        dp = fill(Tv(-Inf), (nf, nc_max + 1))
        pr = fill((0, 1), (nf, nc_max + 1))
        ma = fill(false, nc_max + 1)
        return [m(x, dp, pr, ma, recover, withpos) for x in xs]
end

function (m::NGramHMM)(x::String; o...)
        return m([x]; o...)[1]
end

function (m::NGramHMM)(x::String, dp::Matrix{Tv}, pr::Matrix{Tuple{Int,Int}}, ma::Vector{Bool}, recover::Bool, withpos::Bool)::Vector{String} where {T}

        origin = x
        x = mask(x)
        chrs = codeunits(x) #todo slow?
        vtxs = collect(eachmatch(m.dict, chrs))
        sort!(vtxs)
        nv, nc, np = length.((vtxs, chrs, m.hpr))
        for i = 1:nc + 1; dp[:, i] .= -oo end
        pv = 1
        dp[1:np, 1] .= m.hpr
        prei = 1
        B, M = np + 1, (np + 1) ^ (m.grams - 2)
        ma[1] = true; ma[2:nc + 1] .= false;
        nf = size(m.h2h, 1)

        for i in 1:nc + 1

                if ma[i]; prei = i; end

                while pv <= nv && vtxs[pv].s < i; pv = pv + 1; end

                go = i == nc + 1 || pv <= nv && vtxs[pv].s == i
                if !go; continue; end

                if !ma[i]
                        j = prei
                        fromL, fromR = 1, np
                        for g = 1:m.grams - 1
                                if dp[fromL, j] != -oo
                                        for from = fromL:fromR
                                                jp = from % B
                                                if jp == 0; continue; end
                                                for ip = 1:np
                                                        to = from % M * B + ip
                                                        maybe = dp[from, j] + m.h2h[from, ip] + logph2v(m, jp, 0)
                                                        if dp[to, i] < maybe
                                                                dp[to, i] = maybe
                                                                pr[to, i] = (from, j)
                                                        end
                                                end
                                        end
                                end
                                fromL = fromL * B + 1
                                fromR = fromR * B + np
                        end
                        #prei = i
                        ma[i] = true
                end

                while pv <= nv && vtxs[pv].s == i
                        vtx = vtxs[pv]
                        j = i + length(vtx)
                        ma[j] = true
                        fromL, fromR = 1, np
                        for g = 1:m.grams - 1
                                if dp[fromL, i] != -oo
                                        for from = fromL:fromR
                                                 ip = from % B
                                                 if ip == 0; continue; end
                                                 for jp = 1:np
                                                         to = from % M * B + jp
                                                         #if ip == 0; @show from, B, fromL, fromR, g; end
                                                         maybe = dp[from, i] + m.h2h[from, jp] + logph2v(m, ip, vtx.i)
                                                         if dp[to, j] < maybe
                                                                 dp[to, j] = maybe
                                                                 pr[to, j] = (from, i)
                                                         end
                                                 end
                                        end
                                end
                                #=
                                else
                                        @assert all(isequal(-oo), dp[fromL:fromR, i])
                                end
                                =#
                                fromL = fromL * B + 1
                                fromR = fromR * B + np
                        end
                        pv = pv + 1
                end
        end

        #@show m.h2h
        #@show maximum(dp[:, nc + 1])
        #@show dp
        #@show pr
        #@show vtxs
        v = (argmax(dp[:, nc + 1]), nc + 1)
        words = fill("", 0)
        while v[2] != 1
                pv = pr[v...]
                word = x[pv[2]:prevind(x, v[2])]
                push!(words, word)
                v = pv
        end
        reverse!(words)
        if recover; words = demask(words, origin); end
        ret = words

        if withpos
                v = (argmax(dp[:, nc + 1]), nc + 1)
                postags = fill("", 0)
                while v[2] != 1
                        pv = pr[v...]
                        postag = m.tags[v[1] % np]
                        push!(postags, postag)
                        v = pv
                end
                reverse!(postags)
                ret = [join(p, "/") for p in zip(words, postags)]
        end

        #@show origin
        #@show ret
        return ret
end

logph2v(m::NGramHMM, p, v) = ((1. - m.usrw) * get(m.h2v[p], v, m.INF[p]) + m.usrw * get(m.usrh2v[p], v, m.usrINF[p])) |> log

function ==(a::NGramHMM, b::NGramHMM)
        return all(fname -> getfield(a, fname) == getfield(b, fname), fieldnames(NGramHMM))
end


text:
# ------------------------------------------ #
#  Import functions defined in Rhone2020.jl  #
# ------------------------------------------ #

using Dates, MonteCarloMeasurements;

# ------------------------------------------ #
#    Load the data files into dictionaries   #
# ------------------------------------------ #

ctd145, ctd309, ctd265 = R.load_DCX22_CTD();
e_cond = 200*0.025 # error of conductivity in μS/cm, 2.5% of range FS
e_p = rhow*g*0.0005 # error of pressure, 0.05% FS; has to be multiplied with max. water column, e.g. *100 for ctd309-100mH20

# ------------------------------------------ #
#       Pick start/end points of traces      #
# ------------------------------------------ #

# plot conductivity to pick start and end points of individual injection signals
#using PyPlot
#pygui(true) # interactive window for figure
#plot(ctd309["1008"][:cond])
#plot(ctd265["1008"][:cond])

# manually determined start and end points of individual injection signals
# one column per injection
# 1st row: first index of ctd309
# 2nd row: last index of ctd309
# 3rd row: first index of ctd265
# 4th row: last index of ctd265
# 5th row: first index of ctd145
# 6th row: last index of ctd145

indices = Dict("0808" => ([16882 18408 19638 21934],  # traces 1 and 6-9: supersaturation
                          [17071 19079 19741 22000],
                          [17088 19355 19726 21992],
                          [17339 19474 19869 22080]),
               "0908" => ([8313 10072 11810 13638 15494 17204 18709 21049 22606 24239 26059 27914], # injection 13 too late for ctd265
                          [8386 10120 11900 13710 15560 17270 18784 21130 22682 24300 26140 27990],
                          [8339 10097 11833 13657 15514 17224 18732 21073 22630 24263 26086 27944],
                          [8440 10180 11950 13750 15590 17310 18815 21160 22716 24350 26170 28037]),
               "1008" => ([11829 15747 19984],
                          [11930 15795 20010],
                          [12070 15800 20040],
                          [12320 15880 20090]),
               "1108" => ([5307 7100 8213 9386 10902 11964 13122 18468 19575 25237],
                          [5337 7124 8246 9446 10984 12016 13152 18498 19604 25261],
                          [5329 7139 8273 9459 10960 12006 13164 18519 19631 25290],
                          [5367 7184 8330 9532 11050 12085 13213 18568 19693 25339]),
               "1208" => ([4134 5529 7557 8772 10018 25506],
                          [4181 5574 7592 8811 10050 25562],
                          [4148 5556 7592 8804 10052 25544],
                          [4218 5622 7662 8866 10104 25641]),
               "1308" => ([5437 7555 14943],
                          [5533 7606 14991],
                          [5441 7566 15017],
                          [5549 7641 15110]),
               "1408" => ([4642 5841 7764 9688 19782],
                          [4692 5893 7817 9729 19827],
                          [4649 5857 7805 9725 19822],
                          [4733 5955 7894 9830 19906]),
               "1708" => ([10072 12345 23422],
                          [10149 12437 23563],
                          [10101 12369 23444],
                          [10220 12514 23661]),
               "1808" => ([2416 21321],
                          [2530 21407],
                          [2421 21359],
                          [2540 21526]),
               "1908" => ([1732 2288 27037 28252],
                          [1925 2424 27187 28377],
                          [1740 2301 27091 28286],
                          [1939 2529 27318 28527]),
               "2008" => ([9960],[25900],[15160],[25900],[6213],[10571]), # not usable for calculations
               "2108" => ([722 2554 4313 6102 7894 9744 11763 13359 15143 16950 18758 20472 22224 23230 24193 26416 27983], # 309 and 265: end point = start point of next trace
                          [2554 4313 6102 7894 9744 11763 13359 15143 16950 18758 20472 22224 23230 24193 26416 27983 length(ctd309["2108"][:cond])],
                          [841 2670 4422 6195 7995 9879 11932 13499 15296 17101 18908 20634 22378 23389 24336 26564 28118],
                          [2670 4422 6195 7995 9879 11932 13499 15296 17101 18908 20634 22378 23389 24336 26564 28118 length(ctd265["2108"][:cond])],
                          [622 2446 4223 6021 7821 9627 11656 13233 15031 16830 18629 20363 22104 23122 24080 26307 27893],
                          [688 2512 4310 6100 7893 9760 11770 13355 15150 16943 18737 20461 22208 23228 24158 26400 27980]));

#plot to check if indices are ok
#using PyPlot
#pygui(true)
#date = "2108"
#figure()
#for tr in 1:length(indices[date][1])
#    plot(ctd309[date][:cond][indices[date][1][tr]:indices[date][2][tr]])
#end
#figure()
#for tr in 1:length(indices[date][1])
#    plot(ctd265[date][:cond][indices[date][3][tr]:indices[date][4][tr]])
#end

# ------------------------------------------ #
#            Indices of injections           #
# ------------------------------------------ #

# time of salt injection (translated to indices of ctd arrays)
t0 = Dict("0808" => ([16806 18186 19606 21906]),
          "0908" => ([8256 10026 11766 13596 15456 17166 18666 21006 22566 24196 26016 27876]),
          "1008" => ([11586 15441 19701]),
          "1108" => ([5291 7081 8191 9301 10816 11896 13066 18436 19521 25216]),
          "1208" => ([4126 5521 7546 8761 10006 25501]),
          "1308" => ([5431 7546 14926]),
          "1408" => ([4631 5831 7756 9676 19771]),
          "1708" => ([10066 12346 23411]),
          "1808" => ([2411 21316]),
          "1908" => ([1726 2281 27031 28246]),
          "2008" => ([4441]),
          "2108" => ([601 2416 4201 6001 7801 9601 11626 13201 15001 16801 18606 20341 22081 23101 24061 26281 27871]));
e_dt = 3 # estimated error for time difference in s between watch and ctd clock; dt error between two ctds is 0.5 s

# ------------------------------------------ #
#                Tracer masses               #
# ------------------------------------------ #

# masses of injected salt in kg
mass = Dict("0808" => [0.1;0.1;0.1;0.1],
            "0908" => [0.05;0.05;0.05;0.05;0.05;0.05;0.05;0.05;0.05;0.05;0.05;0.05],
            "1008" => [0.05;0.05;0.05],
            "1108" => [0.05;0.05;0.07;0.07;0.1;0.1;0.1;0.1;0.1;0.1],
            "1208" => [0.05;0.1;0.1;0.1;0.1;0.06],
            "1308" => [0.1;0.1;0.1],
            "1408" => [0.1;0.1;0.1;0.1;0.1],
            "1708" => [0.1;0.1;0.1],
            "1808" => [0.1;0.1],
            "1908" => [0.1;0.1;0.1;0.1],
            "2008" => [0.02],
            "2108" => [0.025;0.025;0.025;0.025;0.05;0.05;0.025;0.025;0.025;0.025;0.025;0.025;0.025;0.05;0.025;0.025;0.025]
            );

# repeatedly measured 100g with analogue suspended scale (as in the field) and compared it to kitchen scale:
m = [99.7,100.3,100.1,99.5,99.9,99.7,101.1,99.7,99.1,98.7,99.8,99.3,100.3,99.8,98.6,100.1,99.9,99.5,100.1,99.2,98.7,100.6];
e_m = 0.002 # error of mass in kg (small influence)

# ------------------------------------------ #
#              Depths of CTDs                #
# ------------------------------------------ #

e_dz = 1.0 # estimated error of vertical distance in m

# depths of ctds in m for each day [ctd309, ctd265, ctd145]
# upper 100m: add 0.04*(depth-100m),; 104m measured under tension
# lower 100m: add 0.01*depth; 101m measured under tension

ctd309["0808"][:depth], ctd265["0808"][:depth] = Particles.(n_partcl, Normal.((62.4, 164.6), e_dz))
ctd309["0908"][:depth], ctd265["0908"][:depth] = Particles.(n_partcl, Normal.((121.1, 171.7), e_dz))
ctd309["1008"][:depth], ctd265["1008"][:depth] = Particles.(n_partcl, Normal.((120.9, 171.7), e_dz))
ctd309["1108"][:depth], ctd265["1108"][:depth] = Particles.(n_partcl, Normal.((121.0, 171.7), e_dz))
ctd309["1208"][:depth], ctd265["1208"][:depth] = Particles.(n_partcl, Normal.((121.0, 171.7), e_dz))
ctd309["1308"][:depth], ctd265["1308"][:depth] = Particles.(n_partcl, Normal.((131.6, 171.9), e_dz))
ctd309["1408"][:depth], ctd265["1408"][:depth] = Particles.(n_partcl, Normal.((131.3, 171.7), e_dz))
ctd309["1708"][:depth], ctd265["1708"][:depth] = Particles.(n_partcl, Normal.((131.3, 171.7), e_dz))
ctd309["1808"][:depth], ctd265["1808"][:depth] = Particles.(n_partcl, Normal.((131.3, 171.7), e_dz))
ctd309["1908"][:depth], ctd265["1908"][:depth] = Particles.(n_partcl, Normal.((141.4, 171.7), e_dz))
ctd145["2008"][:depth], ctd309["2008"][:depth], ctd265["2008"][:depth] = Particles.(n_partcl, Normal.((10.4, 88.4, 189.9), e_dz))
ctd145["2108"][:depth], ctd309["2108"][:depth], ctd265["2108"][:depth] = Particles.(n_partcl, Normal.((10.4, 88.4, 189.9), e_dz))

# ------------------------------------------ #
#                Calibrations                #
# ------------------------------------------ #

# get functions converting conductivities to concentrations as well as shift of temperature sensors from 0°C
include(joinpath(@__DIR__,"calibration.jl"))
e_T = 0.05 # °C, error of temperature measurements introduced by calibration (otherwise it would be 0.1°C for PT sensor)

# ------------------------------------------ #
#           Subtract air pressure            #
# ------------------------------------------ #

stage = R.load_stage_pressure() # pressure sensor installed upstream, contains air pressure; sampling interval 2min
keller_p = R.load_DCX22_pressure() # pressure sensor installed in moulins over night; not used
for key in [collect(keys(indices));"1608"]
    if haskey(ctd309,key)
        local idx_stage309 = indexin(DateTime.(ctd309[key][:t]),DateTime.(stage[:t]))
        idx_stage309 = idx_stage309[idx_stage309.!==nothing]
        local idx_ctd309 = indexin(DateTime.(stage[:t]),DateTime.(ctd309[key][:t]))
        idx_ctd309 = idx_ctd309[idx_ctd309.!==nothing]
        local idx_stage265 = indexin(DateTime.(ctd265[key][:t]),DateTime.(stage[:t]))
        idx_stage265 = idx_stage265[idx_stage265.!==nothing]
        local idx_ctd265 = indexin(DateTime.(stage[:t]),DateTime.(ctd265[key][:t]))
        idx_ctd265 = idx_ctd265[idx_ctd265.!==nothing]
        ctd309[key][:dpress] = ctd309[key][:press][idx_ctd309]-stage[:airpress][idx_stage309]
        ctd309[key][:t_dpress] = ctd309[key][:t][idx_ctd309]
        ctd265[key][:dpress] = ctd265[key][:press][idx_ctd265]-stage[:airpress][idx_stage265]
        ctd265[key][:t_dpress] = ctd265[key][:t][idx_ctd265]
        ctd265[key][:idx_stage] = idx_stage265
    end
    if haskey(ctd145,key)
        local idx_stage145 = indexin(DateTime.(ctd145[key][:t]),DateTime.(stage[:t]))
        idx_stage145 = idx_stage145[idx_stage145.!==nothing]
        local idx_ctd145 = indexin(DateTime.(stage[:t]),DateTime.(ctd145[key][:t]))
        idx_ctd145 = idx_ctd145[idx_ctd145.!==nothing]
        ctd145[key][:dpress] = ctd145[key][:press][idx_ctd145]-stage[:airpress][idx_stage145]
        ctd145[key][:t_dpress] = ctd145[key][:t][idx_ctd145]
    end
    if haskey(keller_p,key) || key == "1608"
        local dt = keller_p[key][:t][1]-stage[:t]
        local i1_stage = findmin(abs.(dt))
        local iend_stage = min(i1_stage[2]+length(keller_p[key][:t])-1, length(stage[:t]))
        local iend_kellerp = iend_stage-i1_stage[2]+1
        keller_p[key][:dpress] = keller_p[key][:press][1:iend_kellerp]-stage[:airpress][i1_stage[2]:iend_stage]
        keller_p[key][:t_dpress] = keller_p[key][:t][1:iend_kellerp]
    end
end
stage[:dpress] = stage[:press]-stage[:airpress];


text:
module ConstantOptimizationModule

using LineSearches: LineSearches
using Optim: Optim
import ..CoreModule: CONST_TYPE, Node, Options, Dataset
import ..UtilsModule: get_birth_order
import ..EquationUtilsModule: get_constants, set_constants, count_constants
import ..LossFunctionsModule: score_func, eval_loss
import ..PopMemberModule: PopMember

# Proxy function for optimization
function opt_func(
    x::Vector{CONST_TYPE}, dataset::Dataset{T}, baseline::T, tree::Node, options::Options
)::T where {T<:Real}
    set_constants(tree, x)
    # TODO(mcranmer): This should use score_func batching.
    loss = eval_loss(tree, dataset, options)
    return loss
end

# Use Nelder-Mead to optimize the constants in an equation
function optimize_constants(
    dataset::Dataset{T}, baseline::T, member::PopMember, options::Options
)::Tuple{PopMember,Float64} where {T<:Real}
    nconst = count_constants(member.tree)
    num_evals = 0.0
    if nconst == 0
        return (member, 0.0)
    end
    x0 = get_constants(member.tree)
    f(x::Vector{CONST_TYPE})::T = opt_func(x, dataset, baseline, member.tree, options)
    if nconst == 1
        algorithm = Optim.Newton(; linesearch=LineSearches.BackTracking())
    else
        if options.optimizer_algorithm == "NelderMead"
            algorithm = Optim.NelderMead(; linesearch=LineSearches.BackTracking())
        elseif options.optimizer_algorithm == "BFGS"
            algorithm = Optim.BFGS(; linesearch=LineSearches.BackTracking())#order=3))
        else
            error("Optimization function not implemented.")
        end
    end
    result = Optim.optimize(
        f, x0, algorithm, Optim.Options(; iterations=options.optimizer_iterations)
    )
    num_evals += result.f_calls
    # Try other initial conditions:
    for i in 1:(options.optimizer_nrestarts)
        new_start =
            x0 .* (
                convert(CONST_TYPE, 1) .+
                convert(CONST_TYPE, 1//2) * randn(CONST_TYPE, size(x0, 1))
            )
        tmpresult = Optim.optimize(
            f,
            new_start,
            algorithm,
            Optim.Options(; iterations=options.optimizer_iterations),
        )
        num_evals += tmpresult.f_calls

        if tmpresult.minimum < result.minimum
            result = tmpresult
        end
    end

    if Optim.converged(result)
        set_constants(member.tree, result.minimizer)
        member.score, member.loss = score_func(dataset, baseline, member.tree, options)
        num_evals += 1
        member.birth = get_birth_order(; deterministic=options.deterministic)
    else
        set_constants(member.tree, x0)
    end
    return member, num_evals
end

end


text:
include("common.jl")

mmesh = Mesh(10, 10, 1/10, degree=1)
NT = 500
Δt = 2.0/NT
ηmax = 1
ηmin = 0.5

bdedge = bcedge((x1, y1, x2, y2)->(x1>1-1e-5 && x2>1-1e-5), mmesh)
bdnode = bcnode((x1, y1)->(y1<1e-5), mmesh)

invη = 50 * constant(ones(get_ngauss(mmesh)))
μ = 1.0
λ = 1.0

Forces = zeros(NT, 2mmesh.ndof)
T = eval_f_on_boundary_edge((x,y)->0.1, bdedge, mmesh)
T = [-T T]
rhs = compute_fem_traction_term(T, bdedge, mmesh)
for i = 1:NT
  Forces[i, :] = rhs
end
Forces = constant(Forces)

U, Sigma, Varepsilon = viscoelasticity_solver(invη, μ, λ, Forces, mmesh)

sess = Session(); init(sess)
Uval, Sigmaval, Varepsilonval = run(sess, [U, Sigma, Varepsilon])
p = visualize_displacement(Uval[1:5:end,:], mmesh)
saveanim(p, "space_u.gif")


text:
@doc raw"""
    AbstractMetric

Abstract type for the pseudo-Riemannian metric tensor ``g``, a family of smoothly
varying inner products on the tangent space. See [`inner`](@ref).

# Functor

    (metric::Metric)(M::AbstractManifold)
    (metric::Metric)(M::MetricManifold)

Generate the `MetricManifold` that wraps the manifold `M` with given `metric`.
This works for both a variable containing the metric as well as a subtype `T<:AbstractMetric`,
where a zero parameter constructor `T()` is availabe.
If `M` is already a metric manifold, the inner manifold with the new `metric` is returned.
"""
abstract type AbstractMetric end

"""
    IsMetricManifold <: AbstractTrait

Specify that a certain decorated Manifold is a metric manifold in the sence that it provides
explicit metric properties, extending/changing the default metric properties of a manifold.
"""
struct IsMetricManifold <: AbstractTrait end

"""
    IsDefaultMetric{G<:AbstractMetric}

Specify that a certain [`AbstractMetric`](@ref) is the default metric for a manifold.
This way the corresponding [`MetricManifold`](@ref) falls back to the default methods
of the manifold it decorates.
"""
struct IsDefaultMetric{G<:AbstractMetric} <: AbstractTrait
    metric::G
end
parent_trait(::IsDefaultMetric) = IsMetricManifold()

# piping syntax for decoration
(metric::AbstractMetric)(M::AbstractManifold) = MetricManifold(M, metric)
(::Type{T})(M::AbstractManifold) where {T<:AbstractMetric} = MetricManifold(M, T())

"""
    MetricManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractMetric} <: AbstractDecoratorManifold{𝔽}

Equip a [`AbstractManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.AbstractManifold) explicitly with a [`AbstractMetric`](@ref) `G`.

For a Metric AbstractManifold, by default, assumes, that you implement the linear form
from [`local_metric`](@ref) in order to evaluate the exponential map.

If the corresponding [`AbstractMetric`](@ref) `G` yields closed form formulae for e.g.
the exponential map and this is implemented directly (without solving the ode),
you can of course still implement that directly.

# Constructor

    MetricManifold(M, G)

Generate the [`AbstractManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.AbstractManifold) `M` as a manifold with the [`AbstractMetric`](@ref) `G`.
"""
struct MetricManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractMetric} <:
       AbstractDecoratorManifold{𝔽}
    manifold::M
    metric::G
end

function active_traits(f, M::MetricManifold, args...)
    return merge_traits(
        is_default_metric(M.manifold, M.metric) ? IsDefaultMetric(M.metric) : EmptyTrait(),
        IsMetricManifold(),
        active_traits(f, M.manifold, args...),
        is_metric_function(f) ? EmptyTrait() : IsExplicitDecorator(),
    )
end
# remetricise instead of double-decorating
(metric::AbstractMetric)(M::MetricManifold) = MetricManifold(M.manifold, metric)
(::Type{T})(M::MetricManifold) where {T<:AbstractMetric} = MetricManifold(M.manifold, T())

@doc raw"""
    RiemannianMetric <: AbstractMetric

Abstract type for Riemannian metrics, a family of positive definite inner
products. The positive definite property means that for ``X  ∈ T_p \mathcal M``, the
inner product ``g(X, X) > 0`` whenever ``X`` is not the zero vector.
"""
abstract type RiemannianMetric <: AbstractMetric end

decorated_manifold(M::MetricManifold) = M.manifold

get_embedding(M::MetricManifold) = get_embedding(M.manifold)

@doc raw"""
    change_metric(M::AbstractcManifold, G2::AbstractMetric, p, X)

On the [`AbstractManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.AbstractManifold) `M` with implicitly given metric ``g_1``
and a second [`AbstractMetric`](@ref) ``g_2`` this function performs a change of metric in the
sense that it returns the tangent vector ``Z=BX`` such that the linear map ``B`` fulfills

````math
g_2(Y_1,Y_2) = g_1(BY_1,BY_2) \quad \text{for all } Y_1, Y_2 ∈ T_p\mathcal M.
````

If both metrics are given in their [`local_metric`](@ref) (symmetric positive defintie) matrix
representations ``G_1 = C_1C_1^{\mathrm{H}}`` and ``G_2 = C_2C_2^{\mathrm{H}}``, where ``C_1,C_2`` denote their respective
Cholesky factors, then solving ``C_2C_2^{\mathrm{H}} = G_2 = B^{\mathrm{H}}G_1B = B^{\mathrm{H}}C_1C_1^{\mathrm{H}}B`` yields ``B = (C_1 \backslash C_2)^{\mathrm{H}}``,
where ``\cdot^{\mathrm{H}}`` denotes the conjugate transpose.

This function returns `Z = BX`.

# Examples

    change_metric(Sphere(2), EuclideanMetric(), p, X)

Since the metric in ``T_p\mathbb S^2`` is the Euclidean metric from the embedding restricted to ``T_p\mathbb S^2``, this just returns `X`

    change_metric(SymmetricPOsitiveDefinite(3), EuclideanMetric, p, X)

Here, the default metric in ``\mathcal P(3)`` is the [`LinearAffineMetric`](@ref) and the transformation can be computed as ``B=p``
"""
change_metric(::AbstractManifold, ::AbstractMetric, ::Any, ::Any)

function change_metric(M::AbstractManifold, G::AbstractMetric, p, X)
    Y = allocate_result(M, change_metric, X, p) # this way we allocate a tangent
    return change_metric!(M, Y, G, p, X)
end
function change_metric!(
    ::T,
    M::AbstractDecoratorManifold,
    Y,
    ::G,
    p,
    X,
) where {G<:AbstractMetric,T<:TraitList{<:IsDefaultMetric{<:G}}}
    return copyto!(M, Y, p, X)
end
function change_metric!(M::MetricManifold, Y, G::AbstractMetric, p, X)
    M.metric === G && return copyto!(M, Y, p, X) # no metric change
    # TODO: For local metric, inverse_local metric, det_local_metric: Introduce a default basis?
    B = DefaultOrthogonalBasis()
    G1 = local_metric(M, p, B)
    G2 = local_metric(G(M), p, B)
    x = get_coordinates(M, p, X, B)
    C1 = cholesky(G1).L
    C2 = cholesky(G2).L
    z = (C1 \ C2)'x
    return get_vector!(M, Y, p, z, B)
end

@trait_function change_metric(M::AbstractDecoratorManifold, G::AbstractMetric, X, p)
@trait_function change_metric!(M::AbstractDecoratorManifold, Y, G::AbstractMetric, X, p)

@doc raw"""
    change_representer(M::AbstractManifold, G2::AbstractMetric, p, X)

Convert the representer `X` of a linear function (in other words a cotangent vector at `p`)
in the tangent space at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.AbstractManifold) `M` given with respect to the
[`AbstractMetric`](@ref) `G2` into the representer with respect to the (implicit) metric of `M`.

In order to convert `X` into the representer with respect to the (implicitly given) metric ``g_1`` of `M`,
we have to find the conversion function ``c: T_p\mathcal M \to T_p\mathcal M`` such that

```math
    g_2(X,Y) = g_1(c(X),Y)
```

If both metrics are given in their [`local_metric`](@ref) (symmetric positive defintie) matrix
representations ``G_1`` and ``G_2`` and ``x,y`` are the local coordinates with respect to
the same basis of the tangent space, the equation reads

```math
   x^{\mathrm{H}}G_2y = c(x)^{\mathrm{H}}G_1 y \quad \text{for all } y \in ℝ^d,
```
where ``\cdot^{\mathrm{H}}`` denotes the conjugate transpose.
We obtain ``c(X) = (G_1\backslash G_2)^{\mathrm{H}}X``

For example `X` could be the gradient ``\operatorname{grad}f`` of a real-valued function
``f: \mathcal M \to ℝ``, i.e.

```math
    g_2(X,Y) = Df(p)[Y] \quad \text{for all } Y ∈ T_p\mathcal M.
```

and we would change the Riesz representer `X` to the representer with respect to the metric ``g_1``.

# Examples

    change_representer(Sphere(2), EuclideanMetric(), p, X)

Since the metric in ``T_p\mathbb S^2`` is the Euclidean metric from the embedding restricted to ``T_p\mathbb S^2``, this just returns `X`

    change_representer(SymmetricPositiveDefinite(3), EuclideanMetric(), p, X)

Here, the default metric in ``\mathcal P(3)`` is the [`LinearAffineMetric`](@ref) and the transformation can be computed as ``pXp``
"""
change_representer(::AbstractManifold, ::AbstractMetric, ::Any, ::Any)

function change_representer(M::AbstractManifold, G::AbstractMetric, p, X)
    Y = allocate_result(M, change_representer, X, p) # this way we allocate a tangent
    return change_representer!(M, Y, G, p, X)
end

@trait_function change_representer(M::AbstractDecoratorManifold, G::AbstractMetric, X, p)
@trait_function change_representer!(
    M::AbstractDecoratorManifold,
    Y,
    G::AbstractMetric,
    X,
    p,
)

# Default fallback II: Default metric (not yet hit, check subtyping?)
function change_representer!(
    ::T,
    M::AbstractDecoratorManifold,
    Y,
    ::G,
    p,
    X,
) where {G<:AbstractMetric,T<:TraitList{<:IsDefaultMetric{<:G}}}
    return copyto!(M, Y, p, X)
end
# Default fallback II: compute in local metric representations
function change_representer!(M::AbstractManifold, Y, G::AbstractMetric, p, X)
    M.metric === G && return copyto!(M, Y, p, X) # no metric change
    # TODO: For local metric, inverse_local metric, det_local_metric: Introduce a default basis?
    B = DefaultOrthogonalBasis()
    G1 = local_metric(M, p, B)
    G2 = local_metric(G(M), p, B)
    x = get_coordinates(M, p, X, B)
    z = (G1 \ G2)'x
    return get_vector!(M, Y, p, z, B)
end

"""
    connection(::MetricManifold)

Return the [`LeviCivitaConnection`](@ref) for a metric manifold.
"""
connection(::MetricManifold) = LeviCivitaConnection()

default_retraction_method(M::MetricManifold) = default_retraction_method(M.manifold)

@doc raw"""
    det_local_metric(M::AbstractManifold, p, B::AbstractBasis)

Return the determinant of local matrix representation of the metric tensor ``g``, i.e. of the
matrix ``G(p)`` representing the metric in the tangent space at ``p`` with as a matrix.

See also [`local_metric`](@ref)
"""
function det_local_metric(M::AbstractManifold, p, B::AbstractBasis)
    return det(local_metric(M, p, B))
end
@trait_function det_local_metric(M::AbstractDecoratorManifold, p, B::AbstractBasis)

function exp!(::TraitList{IsMetricManifold}, M::AbstractDecoratorManifold, q, p, X)
    return retract!(
        M,
        q,
        p,
        X,
        ODEExponentialRetraction(ManifoldsBase.default_retraction_method(M)),
    )
end

"""
    einstein_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_badefault_differential_backendckend())

Compute the Einstein tensor of the manifold `M` at the point `p`, see [https://en.wikipedia.org/wiki/Einstein_tensor](https://en.wikipedia.org/wiki/Einstein_tensor)
"""
function einstein_tensor(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=default_differential_backend(),
)
    Ric = ricci_tensor(M, p, B; backend=backend)
    g = local_metric(M, p, B)
    Ginv = inverse_local_metric(M, p, B)
    S = sum(Ginv .* Ric)
    G = Ric - g .* S / 2
    return G
end
@trait_function einstein_tensor(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    flat(N::MetricManifold{M,G}, p, X::TFVector)

Compute the musical isomorphism to transform the tangent vector `X` from the
[`AbstractManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.AbstractManifold) `M` equipped with [`AbstractMetric`](@ref) `G` to a cotangent by
computing

````math
X^♭= G_p X,
````
where ``G_p`` is the local matrix representation of `G`, see [`local_metric`](@ref)
"""
flat(::MetricManifold, ::Any, ::TFVector)

function flat!(
    ::TraitList{IsMetricManifold},
    M::AbstractDecoratorManifold,
    ξ::CoTFVector,
    p,
    X::TFVector,
)
    g = local_metric(M, p, ξ.basis)
    copyto!(ξ.data, g * X.data)
    return ξ
end
function flat!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    ξ::CoTFVector,
    p,
    X::TFVector,
) where {𝔽,TM<:AbstractManifold,G<:AbstractMetric}
    flat!(M.manifold, ξ, p, X)
    return ξ
end

function get_basis(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    B::AbstractBasis,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return get_basis(M.manifold, p, B)
end

function get_coordinates(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    X,
    B::AbstractBasis,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return get_coordinates(M.manifold, p, X, B)
end
function get_coordinates!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    Y,
    p,
    X,
    B::AbstractBasis,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return get_coordinates!(M.manifold, Y, p, X, B)
end

function get_vector(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    c,
    B::AbstractBasis,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return get_vector(M.manifold, p, c, B)
end
function get_vector!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    Y,
    p,
    c,
    B::AbstractBasis,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return get_vector!(M.manifold, Y, p, c, B)
end

@doc raw"""
    inverse_local_metric(M::AbstractcManifold{𝔽}, p, B::AbstractBasis)

Return the local matrix representation of the inverse metric (cometric) tensor
of the tangent space at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.AbstractManifold) `M` with respect
to the [`AbstractBasis`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/bases.html#ManifoldsBase.AbstractBasis) basis `B`.

The metric tensor (see [`local_metric`](@ref)) is usually denoted by ``G = (g_{ij}) ∈ 𝔽^{d×d}``,
where ``d`` is the dimension of the manifold.

Then the inverse local metric is denoted by ``G^{-1} = g^{ij}``.
"""
inverse_local_metric(::AbstractManifold, ::Any, ::AbstractBasis)
function inverse_local_metric(M::AbstractManifold, p, B::AbstractBasis)
    return inv(local_metric(M, p, B))
end
@trait_function inverse_local_metric(M::AbstractDecoratorManifold, p, B::AbstractBasis)

function Base.convert(::Type{MetricManifold{𝔽,MT,GT}}, M::MT) where {𝔽,MT,GT}
    return _convert_with_default(M, GT, Val(is_default_metric(M, GT())))
end

function _convert_with_default(
    M::MT,
    T::Type{<:AbstractMetric},
    ::Val{true},
) where {MT<:AbstractManifold}
    return MetricManifold(M, T())
end
function _convert_with_default(
    M::MT,
    T::Type{<:AbstractMetric},
    ::Val{false},
) where {MT<:AbstractManifold}
    return error(
        "Can not convert $(M) to a MetricManifold{$(MT),$(T)}, since $(T) is not the default metric.",
    )
end

function exp(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    X,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return exp(M.manifold, p, X)
end
function exp!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    q,
    p,
    X,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return exp!(M.manifold, q, p, X)
end

injectivity_radius(M::MetricManifold) = injectivity_radius(M.manifold)
function injectivity_radius(M::MetricManifold, m::AbstractRetractionMethod)
    return injectivity_radius(M.manifold, m)
end

@doc raw"""
    inner(N::MetricManifold{M,G}, p, X, Y)

Compute the inner product of `X` and `Y` from the tangent space at `p` on the
[`AbstractManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.AbstractManifold) `M` using the [`AbstractMetric`](@ref) `G`.
If `M` has `G` as its [`IsDefaultMetric`](@ref) trait,
this is done using `inner(M, p, X, Y)`, otherwise the [`local_metric`](@ref)`(M, p)` is employed as

````math
g_p(X, Y) = ⟨X, G_p Y⟩,
````
where ``G_p`` is the loal matrix representation of the [`AbstractMetric`](@ref) `G`.
"""
inner(::MetricManifold, ::Any, ::Any, ::Any)

function inner(
    ::TraitList{IsMetricManifold},
    M::AbstractDecoratorManifold,
    p,
    X::TFVector,
    Y::TFVector,
)
    X.basis === Y.basis ||
        error("calculating inner product of vectors from different bases is not supported")
    return dot(X.data, local_metric(M, p, X.basis) * Y.data)
end
function inner(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    X,
    Y,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return inner(M.manifold, p, X, Y)
end

"""
    is_default_metric(M::AbstractManifold, G::AbstractMetric)

returns whether an [`AbstractMetric`](@ref) is the default metric on the manifold `M` or not.
This can be set by defining this function, or setting the [`IsDefaultMetric`](@ref) trait for an
[`AbstractDecoratorManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/decorator.html#ManifoldsBase.AbstractDecoratorManifold).
"""
is_default_metric(M::AbstractManifold, G::AbstractMetric)

@trait_function is_default_metric(M::AbstractDecoratorManifold, G::AbstractMetric)
function is_default_metric(
    ::TraitList{IsDefaultMetric{G}},
    ::AbstractDecoratorManifold,
    ::G,
) where {G<:AbstractMetric}
    return true
end
is_default_metric(M::MetricManifold) = is_default_metric(M.manifold, M.metric)
is_default_metric(::AbstractManifold, ::AbstractMetric) = false

function is_point(
    ::TraitList{IsMetricManifold},
    M::MetricManifold{𝔽,TM,G},
    p,
    te=false;
    kwargs...,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return is_point(M.manifold, p, te; kwargs...)
end

function is_vector(
    ::TraitList{IsMetricManifold},
    M::MetricManifold{𝔽,TM,G},
    p,
    X,
    te=false,
    cbp=true;
    kwargs...,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return is_vector(M.manifold, p, X, te, cbp; kwargs...)
end

@doc raw"""
    local_metric(M::AbstractManifold{𝔽}, p, B::AbstractBasis)

Return the local matrix representation at the point `p` of the metric tensor ``g`` with
respect to the [`AbstractBasis`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/bases.html#ManifoldsBase.AbstractBasis) `B` on the [`AbstractManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.AbstractManifold) `M`.
Let ``d``denote the dimension of the manifold and $b_1,\ldots,b_d$ the basis vectors.
Then the local matrix representation is a matrix ``G\in 𝔽^{n\times n}`` whose entries are
given by ``g_{ij} = g_p(b_i,b_j), i,j\in\{1,…,d\}``.

This yields the property for two tangent vectors (using Einstein summation convention)
``X = X^ib_i, Y=Y^ib_i \in T_p\mathcal M`` we get ``g_p(X, Y) = g_{ij} X^i Y^j``.
"""
local_metric(::AbstractManifold, ::Any, ::AbstractBasis)
@trait_function local_metric(M::AbstractDecoratorManifold, p, B::AbstractBasis; kwargs...)

function local_metric(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    B::AbstractBasis,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return local_metric(M.manifold, p, B)
end

@doc raw"""
    local_metric_jacobian(
        M::AbstractManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend,
    )

Get partial derivatives of the local metric of `M` at `p` in basis `B` with respect to the
coordinates of `p`, ``\frac{∂}{∂ p^k} g_{ij} = g_{ij,k}``. The
dimensions of the resulting multi-dimensional array are ordered ``(i,j,k)``.
"""
local_metric_jacobian(::AbstractManifold, ::Any, B::AbstractBasis, ::AbstractDiffBackend)
function local_metric_jacobian(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=default_differential_backend(),
)
    n = size(p, 1)
    ∂g = reshape(_jacobian(q -> local_metric(M, q, B), p, backend), n, n, n)
    return ∂g
end
@trait_function local_metric_jacobian(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    log(N::MetricManifold{M,G}, p, q)

Copute the logarithmic map on the [`AbstractManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.AbstractManifold) `M` equipped with the [`AbstractMetric`](@ref) `G`.

If the metric was declared the default metric using the [`IsDefaultMetric`](@ref) trait or [`is_default_metric`](@ref), this method
falls back to `log(M,p,q)`. Otherwise, you have to provide an implementation for the non-default
[`AbstractMetric`](@ref) `G` metric within its [`MetricManifold`](@ref)`{M,G}`.
"""
log(::MetricManifold, ::Any...)

function log(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    q,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return log(M.manifold, p, q)
end
function log!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    X,
    p,
    q,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return log!(M.manifold, X, p, q)
end

@doc raw"""
    log_local_metric_density(M::AbstractManifold, p, B::AbstractBasis)

Return the natural logarithm of the metric density ``ρ`` of `M` at `p`, which
is given by ``ρ = \log \sqrt{|\det [g_{ij}]|}`` for the metric tensor expressed in basis `B`.
"""
log_local_metric_density(::AbstractManifold, ::Any, ::AbstractBasis)
function log_local_metric_density(M::AbstractManifold, p, B::AbstractBasis)
    return log(abs(det_local_metric(M, p, B))) / 2
end
@trait_function log_local_metric_density(M::AbstractDecoratorManifold, p, B::AbstractBasis)

manifold_dimension(M::MetricManifold) = manifold_dimension(M.manifold)

@doc raw"""
    metric(M::MetricManifold)

Get the metric ``g`` of the manifold `M`.
"""
metric(::MetricManifold)

function metric(M::MetricManifold)
    return M.metric
end

function norm(::TraitList{IsMetricManifold}, M::AbstractDecoratorManifold, p, X::TFVector)
    return sqrt(dot(X.data, local_metric(M, p, X.basis) * X.data))
end

function parallel_transport_to(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    X,
    q,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return parallel_transport_to(M.manifold, p, X, q)
end
function parallel_transport_to!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    Y,
    p,
    X,
    q,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return parallel_transport_to!(M.manifold, Y, p, X, q)
end

function project(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return project(M.manifold, p)
end
function project!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    q,
    p,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return project!(M.manifold, q, p)
end
function project(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    X,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return project(M.manifold, p, X)
end
function project!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    Y,
    p,
    X,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return project!(M.manifold, Y, p, X)
end

representation_size(M::MetricManifold) = representation_size(M.manifold)

@doc raw"""
    ricci_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend())

Compute the Ricci scalar curvature of the manifold `M` at the point `p` using basis `B`.
The curvature is computed as the trace of the Ricci curvature tensor with respect to
the metric, that is ``R=g^{ij}R_{ij}`` where ``R`` is the scalar Ricci curvature at `p`,
``g^{ij}`` is the inverse local metric (see [`inverse_local_metric`](@ref)) at `p` and
``R_{ij}`` is the Riccie curvature tensor, see [`ricci_tensor`](@ref). Both the tensor and
inverse local metric are expressed in local coordinates defined by `B`, and the formula
uses the Einstein summation convention.
"""
ricci_curvature(::AbstractManifold, ::Any, ::AbstractBasis)
function ricci_curvature(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=default_differential_backend(),
)
    Ginv = inverse_local_metric(M, p, B)
    Ric = ricci_tensor(M, p, B; backend=backend)
    S = sum(Ginv .* Ric)
    return S
end
@trait_function ricci_curvature(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    sharp(N::MetricManifold{M,G}, p, ξ::CoTFVector)

Compute the musical isomorphism to transform the cotangent vector `ξ` from the
[`AbstractManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.AbstractManifold) `M` equipped with [`AbstractMetric`](@ref) `G` to a tangent by
computing

````math
ξ^♯ = G_p^{-1} ξ,
````
where ``G_p`` is the local matrix representation of `G`, i.e. one employs
[`inverse_local_metric`](@ref) here to obtain ``G_p^{-1}``.
"""
sharp(::MetricManifold, ::Any, ::CoTFVector)

function sharp!(
    ::TraitList{IsMetricManifold},
    M::AbstractDecoratorManifold,
    X::TFVector,
    p,
    ξ::CoTFVector,
)
    Ginv = inverse_local_metric(M, p, X.basis)
    copyto!(X.data, Ginv * ξ.data)
    return X
end
function sharp!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    X::TFVector,
    p,
    ξ::CoTFVector,
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    sharp!(M.manifold, X, p, ξ)
    return X
end

function Base.show(io::IO, M::MetricManifold)
    return print(io, "MetricManifold($(M.manifold), $(M.metric))")
end
function Base.show(io::IO, i::IsDefaultMetric)
    return print(io, "IsDefaultMetric($(i.metric))")
end

function vector_transport_along(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    X,
    c::AbstractVector,
    m::AbstractVectorTransportMethod=default_vector_transport_method(M),
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return vector_transport_along(M.manifold, p, X, c, m)
end
function vector_transport_along!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    Y,
    p,
    X,
    c::AbstractVector,
    m::AbstractVectorTransportMethod=default_vector_transport_method(M),
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return vector_transport_along!(M.manifold, Y, p, X, c, m)
end

function vector_transport_direction(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    X,
    d,
    m::AbstractVectorTransportMethod=default_vector_transport_method(M),
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return vector_transport_direction(M.manifold, p, X, d, m)
end
function vector_transport_direction!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    Y,
    p,
    X,
    d,
    m::AbstractVectorTransportMethod=default_vector_transport_method(M),
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return vector_transport_direction!(M.manifold, Y, p, X, d, m)
end

function vector_transport_to(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    p,
    X,
    q,
    m::AbstractVectorTransportMethod=default_vector_transport_method(M),
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return vector_transport_to(M.manifold, p, X, q, m)
end
function vector_transport_to!(
    ::TraitList{IsDefaultMetric{G}},
    M::MetricManifold{𝔽,TM,G},
    Y,
    p,
    X,
    q,
    m::AbstractVectorTransportMethod=default_vector_transport_method(M),
) where {𝔽,G<:AbstractMetric,TM<:AbstractManifold}
    return vector_transport_to!(M.manifold, Y, p, X, q, m)
end

zero_vector(M::MetricManifold, p) = zero_vector(M.manifold, p)
zero_vector!(M::MetricManifold, X, p) = zero_vector!(M.manifold, X, p)

is_metric_function(::Any) = false
for mf in [
    change_metric,
    change_metric!,
    change_representer,
    change_representer!,
    christoffel_symbols_first,
    christoffel_symbols_second,
    christoffel_symbols_second_jacobian,
    det_local_metric,
    einstein_tensor,
    exp,
    exp!,
    flat!,
    gaussian_curvature,
    get_basis,
    get_coordinates,
    get_coordinates!,
    get_vector,
    get_vector!,
    get_vectors,
    inner,
    inverse_local_metric,
    inverse_retract,
    inverse_retract!,
    local_metric,
    local_metric_jacobian,
    log,
    log!,
    log_local_metric_density,
    mean,
    mean!,
    median,
    median!,
    mid_point,
    norm,
    parallel_transport_along,
    parallel_transport_along!,
    parallel_transport_direction,
    parallel_transport_direction!,
    parallel_transport_to,
    parallel_transport_to!,
    retract,
    retract!,
    ricci_curvature,
    ricci_tensor,
    riemann_tensor,
    sharp!,
    vector_transport_along,
    vector_transport_along!,
    vector_transport_direction,
    vector_transport_direction!,
    vector_transport_to,
    vector_transport_to!,
]
    @eval is_metric_function(::typeof($mf)) = true
end


text:
function user_input(prompt)
  print(prompt)
  answer = readline()
  chomp(answer)
end


text:
using Gadfly
using DataFrames
using VML

include("../utility.jl")

function main()
    N = [1,2,3,4,5,6,7,8,9,
         10,20,30,40,50,60,70,80,90,100,200,300,400,500,750,
         1_000,2_500,5_000,7_500,10_000,100_000,1_000_000]

    df = DataFrame(N=Int64[],ALLOC=Float64[],VJULIA=Float64[],BJULIA=Float64[],VML=Float64[])

    @inbounds for i in 1:length(N)
        T = Float64
        n = N[i]
        a = rand(n)
        dest = zeros(n)

        sqrt(a)

        nrep = 100
        ncycle = 10_000 ÷ n + 1

        # n < 100 ? gc_enable(false) : gc_enable(true)

        # gc()
        talloc = 0.0
        for j in 1:nrep
            talloc += @elapsed for k in 1:ncycle Vector{T}(n) end
        end
        talloc = talloc / nrep / ncycle / n * 1e9

        # gc()
        tvjulia = 0.0
        for j in 1:nrep
            tvjulia += @elapsed for k in 1:ncycle sqrt(a) end
        end
        tvjulia = tvjulia / nrep / ncycle / n * 1e9

        # gc()
        tbjulia = 0.0
        # for j in 1:nrep
        #     tbjulia += @elapsed for k in 1:ncycle broadcast!(sqrt,dest,a) end
        # end
        # tbjulia = tbjulia / nrep / ncycle / n * 1e9

        # gc()
        tvml = 0.0
        # for j in 1:nrep
        #     tvml += @elapsed for k in 1:ncycle VML.sqrt!(dest,a) end
        # end
        # tvml = tvml / nrep / ncycle / n * 1e9

        # n < 100 ? gc_enable(false) : gc_enable(true)

        push!(df,[n,talloc, tvjulia, tbjulia, tvml])

    end
    df
end

df = main();
df[:diff] = df[:VJULIA]-df[:ALLOC]
println(df)

layers = Layer[]

bench_alloc(),bench_jvect(),bench_jbroadcast!(),bench_mkl(),bench_mkl!(),bench_jvect2()

df = DataFrame() ; df[:N], df[:CPU] = bench_alloc()
push!(layers,layer(df, x="N", y="CPU", Geom.line, Theme(default_color=color("blue")))[1])

df = DataFrame() ; df[:N], df[:CPU] = bench_jvect()
push!(layers,layer(df, x="N", y="CPU", Geom.line, Theme(default_color=color("yellow")))[1])

df = DataFrame() ; df[:N], df[:CPU] = bench_jbroadcast!()
push!(layers,layer(df, x="N", y="CPU", Geom.line, Theme(default_color=color("green")))[1])

df = DataFrame() ; df[:N], df[:CPU] = bench_mkl()
push!(layers,layer(df, x="N", y="CPU", Geom.line, Theme(default_color=color("orange")))[1])

df = DataFrame() ; df[:N], df[:CPU] = bench_mkl!()
push!(layers,layer(df, x="N", y="CPU", Geom.line, Theme(default_color=color("red")))[1])

df = DataFrame() ; df[:N], df[:CPU] = bench_jvect2()
push!(layers,layer(df, x="N", y="CPU", Geom.line, Theme(default_color=color("black")))[1])

println(df)

p = Gadfly.plot(layers,
                # Scale.y_log10,
                Scale.x_log10,
                Guide.xlabel("n-element vector"),
                Guide.ylabel("CPU time in nsec/element"),
                Guide.title("CPU time with n elements"),
                Guide.manual_color_key("",
                ["ALLOC", "Base.sqrt","broadcast!","VML.sqrt","VML.sqrt!","ibench"],
                ["blue","yellow","green","orange","red","black"])
                )

path = Pkg.dir("MKL") * "/benchmark/sqrt_cpu_regarding_n/"
draw(PNG(path*"sqrt_map!.png", 20cm, 20cm), p)


text:
base_dir = string(dirname(@__FILE__))

@testset "read_data" begin
    include(joinpath(base_dir,"data_5bus_pu.jl"));
    include(joinpath(base_dir,"data_14bus_pu.jl"))
end