text:
# Sioux Falls network data
# http://www.bgu.ac.il/~bargera/tntp/

#Link travel time = free flow time * ( 1 + b * (flow/capacity)^Power ).
#Link generalized cost = Link travel time + toll_factor * toll + distance_factor * distance

# Traffic Assignment Data structure
mutable struct TA_Data
    network_name::String

    number_of_zones::Int
    number_of_nodes::Int
    first_thru_node::Int
    number_of_links::Int

    init_node::Array{Int,1}
    term_node::Array{Int,1}
    capacity::Array{Float64,1}
    link_length::Array{Float64,1}
    free_flow_time::Array{Float64,1}
    b::Array{Float64,1}
    power::Array{Float64,1}
    speed_limit::Array{Float64,1}
    toll::Array{Float64,1}
    link_type::Array{Int64,1}

    total_od_flow::Float64

    travel_demand::Array{Float64,2}
    od_pairs::Array{Tuple{Int64,Int64},1}

    toll_factor::Float64
    distance_factor::Float64

    best_objective::Float64
end

function net_dataframe(td::TA_Data)
    df = DataFrame()
    df.init_node = td.init_node
    df.term_node = td.term_node
    df.capacity = td.capacity
    df.link_length = td.link_length
    df.free_flow_time = td.free_flow_time
    df.b = td.b
    df.power = td.power
    df.speed_limit = td.speed_limit
    df.toll = td.toll
    df.link_type = td.link_type
    return df
end

search_sc(s,c) = something(findfirst(isequal(c), s), 0)

function download_tntp(force_download=false)
  ta_root_dir = joinpath(dirname(dirname(@__FILE__)))
  dest = ta_root_dir
  data_dir = joinpath(dest, "TransportationNetworks-$(TNTP_SHA)")

  # Download
  if !isdir(data_dir) || force_download
    if isdir(data_dir)
      rm(data_dir)
    end
    file = joinpath(ta_root_dir, "tntp.zip")
    dl = download("https://github.com/bstabler/TransportationNetworks/archive/$(TNTP_SHA).zip", file)
    run(unpack_cmd(file, dest, ".zip", ""))
    rm(file)
  end

  return data_dir
end

function read_ta_network(network_name)
  tntp_dir = download_tntp()
  network_dir = joinpath(tntp_dir, network_name)

  @assert ispath(network_dir)

  network_data_file = ""
  trip_table_file = ""

  for f in readdir(network_dir)
    if occursin(".zip", lowercase(f))
      zipfile = joinpath(network_dir, f)
      run(unpack_cmd(zipfile, network_dir, ".zip", ""))
      rm(zipfile)
    end
  end

  for f in readdir(network_dir)
    if occursin("_net", lowercase(f)) && occursin(".tntp", lowercase(f))
      network_data_file = joinpath(network_dir, f)
  elseif occursin("_trips", lowercase(f)) && occursin(".tntp", lowercase(f))
      trip_table_file = joinpath(network_dir, f)
    end
  end

  @assert network_data_file != ""
  @assert trip_table_file != ""

  return network_data_file, trip_table_file
end

function load_ta_network(network_name; best_objective=-1.0, toll_factor=0.0, distance_factor=0.0)
  network_data_file, trip_table_file = read_ta_network(network_name)

  load_ta_network(network_name, network_data_file, trip_table_file, best_objective=best_objective, toll_factor=toll_factor, distance_factor=distance_factor)
end

function load_ta_network(network_name, network_data_file, trip_table_file; best_objective=-1.0, toll_factor=0.0, distance_factor=0.0)

    @assert ispath(network_data_file)
    @assert ispath(trip_table_file)

    ##################################################
    # Network Data
    ##################################################

    number_of_zones = 0
    number_of_links = 0
    number_of_nodes = 0
    first_thru_node = 0

    n = open(network_data_file, "r")

    while (line=readline(n)) != ""
        if occursin("<NUMBER OF ZONES>", line)
            number_of_zones = parse(Int, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<NUMBER OF NODES>", line)
            number_of_nodes = parse(Int, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<FIRST THRU NODE>", line)
            first_thru_node = parse(Int, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<NUMBER OF LINKS>", line)
            number_of_links = parse(Int, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<END OF METADATA>", line)
            break
        end
    end

    @assert number_of_links > 0

    init_node = Array{Int64}(undef, number_of_links)
    term_node = Array{Int64}(undef, number_of_links)
    capacity = zeros(number_of_links)
    link_length = zeros(number_of_links)
    free_flow_time = zeros(number_of_links)
    b = zeros(number_of_links)
    power = zeros(number_of_links)
    speed_limit = zeros(number_of_links)
    toll = zeros(number_of_links)
    link_type = Array{Int64}(undef, number_of_links)

    idx = 1
    while !eof(n)
      line = readline(n)
        if occursin("~", line) || line == ""
            continue
        end

        if occursin(";", line)
            line = strip(line, [' ', '\n', ';'])
            line = replace(line, ";" => "")

            numbers = split(line)
            init_node[idx] = parse(Int64, numbers[1])
            term_node[idx] = parse(Int64, numbers[2])
            capacity[idx] = parse(Float64, numbers[3])
            link_length[idx] = parse(Float64, numbers[4])
            free_flow_time[idx] = parse(Float64, numbers[5])
            b[idx] = parse(Float64, numbers[6])
            power[idx] = parse(Float64, numbers[7])
            speed_limit[idx] = parse(Float64, numbers[8])
            toll[idx] = parse(Float64, numbers[9])
            link_type[idx] = Int(round(parse(Float64, numbers[10])))

            idx = idx + 1
        end
    end

    ##################################################
    # Trip Table
    ##################################################

    number_of_zones_trip = 0
    total_od_flow = 0

    f = open(trip_table_file, "r")

    while (line=readline(f)) != ""
        if occursin("<NUMBER OF ZONES>", line)
            number_of_zones_trip = parse(Int, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<TOTAL OD FLOW>", line)
            total_od_flow = parse(Float64, line[ search_sc(line, '>')+1 : end ] )
        elseif occursin("<END OF METADATA>", line)
            break
        end
    end

    @assert number_of_zones_trip == number_of_zones # Check if number_of_zone is same in both txt files
    @assert total_od_flow > 0

    travel_demand = zeros(number_of_zones, number_of_zones)
    od_pairs = Array{Tuple{Int64, Int64}}(undef, 0)

    origin = -1

    while !eof(f)
        line = readline(f)

        if line == ""
            origin = -1
            continue
        elseif occursin("Origin", line)
            origin = parse(Int, split(line)[2] )
        elseif occursin(";", line)
            pairs = split(line, ";")
            for i=1:size(pairs)[1]
                if occursin(":", pairs[i])
                    pair = split(pairs[i], ":")
                    destination = parse(Int64, strip(pair[1]) )
                    od_flow = parse(Float64, strip(pair[2]) )

                    # println("origin=$origin, destination=$destination, flow=$od_flow")

                    travel_demand[origin, destination] = od_flow
                    push!(od_pairs, (origin, destination))
                end
            end
        end
    end

    # Preparing data to return
    ta_data = TA_Data(
        network_name,
        number_of_zones,
        number_of_nodes,
        first_thru_node,
        number_of_links,
        init_node,
        term_node,
        capacity,
        link_length,
        free_flow_time,
        b,
        power,
        speed_limit,
        toll,
        link_type,
        total_od_flow,
        travel_demand,
        od_pairs,
        toll_factor,
        distance_factor,
        best_objective)

    return ta_data

end # end of load_network function

function read_ta_summary(network_data_file)
  @assert ispath(network_data_file)

  number_of_zones = 0
  number_of_links = 0
  number_of_nodes = 0
  first_thru_node = 0

  n = open(network_data_file, "r")

  search_sc(s,c) = something(findfirst(isequal(c), s), 0)

  while (line=readline(n)) != ""
      if occursin("<NUMBER OF ZONES>", line)
          number_of_zones = parse(Int, line[ search_sc(line, '>')+1 : end ] )
      elseif occursin("<NUMBER OF NODES>", line)
          number_of_nodes = parse(Int, line[ search_sc(line, '>')+1 : end ] )
      elseif occursin("<FIRST THRU NODE>", line)
          first_thru_node = parse(Int, line[ search_sc(line, '>')+1 : end ] )
      elseif occursin("<NUMBER OF LINKS>", line)
          number_of_links = parse(Int, line[ search_sc(line, '>')+1 : end ] )
      elseif occursin("<END OF METADATA>", line)
          break
      end
  end

  return number_of_zones, number_of_links, number_of_nodes
end

function summarize_ta_data(;markdown=false)
  data_dir = download_tntp()

  # Test
  df = DataFrame(Network = String[], Zones = Int[], Links = Int[], Nodes = Int[])
  dic = OrderedDict()
  for d in readdir(data_dir)
    if isdir(joinpath(data_dir, d))
      try
        network_data_file, trip_table_file = read_ta_network(d)
        number_of_zones, number_of_links, number_of_nodes = read_ta_summary(network_data_file)
        dic[d] = (number_of_zones, number_of_links, number_of_nodes)
        push!(df, (d, number_of_zones, number_of_links, number_of_nodes))
      catch

      end

    end
  end

  max_len = [0, 0, 0, 0]
  for net in sort(dic)
    max_len[1] = max( length(net[1]), max_len[1] )
    max_len[2] = max( length(digits(net[2][1])), max_len[2] )
    max_len[3] = max( length(digits(net[2][2])), max_len[3] )
    max_len[4] = max( length(digits(net[2][3])), max_len[4] )
  end

  function format(mlen, val)
    len = 0
    str = ""
    if isa(val, String)
      len = length(val)
      str = val * " "^(mlen - len + 1 )
    elseif isa(val, Number)
      len = length(digits(val))
      str = " "^(mlen - len + 1) * string(val)
    end
    return str
  end

  if markdown
      println("-"^(17+sum(max_len)))
      println("| $(format(max_len[1],"Network")) | $(format(max_len[2],"Zones")) | $(format(max_len[3],"Links")) | $(format(max_len[4],"Nodes")) |")
      println("| $(format(max_len[1],":---")) | $(format(max_len[2],"---:")) | $(format(max_len[3],"---:")) | $(format(max_len[4],"---:")) |")
      for net in sort(dic)
        println("| $(format(max_len[1],net[1])) | $(format(max_len[2],net[2][1])) | $(format(max_len[3],net[2][2])) | $(format(max_len[4],net[2][3])) |")
      end
      println("-"^(17+sum(max_len)))
  end

  if !markdown
      return df
  end
end # end of summarize_ta_data


text:
using ExponentialAction, LinearAlgebra, SparseArrays, Test

Tsingle = (Float32, ComplexF32)
Tdouble = (Float64, ComplexF64)

@testset "expv" begin
    n = 10
    @testset "expv(t::$Tt, A::$MT{$TA}, B::Array{$TB,$(length(Bdims2)+1)}), tscale=$tscale, shift=$shift" for
        MT in (Matrix, Diagonal, Bidiagonal, Tridiagonal),
        Tset in (Tsingle, Tdouble),
        Bdims2 in ((), (4,)),
        Tt in Tset,
        TA in Tset,
        TB in Tset,
        tscale in Tt.((0.1, 1, 10)),
        shift in (true, false)

        t = tscale * randn(Tt)
        if MT <: Bidiagonal
            A = MT(randn(TA, n, n), :U)
        else
            A = MT(randn(TA, n, n))
        end
        B = randn(TB, n, Bdims2...)
        T = Base.promote_eltype(t, A, B)
        rT = real(T)
        @inferred expv(t, A, B; shift=shift)
        @test expv(t, A, B; shift=shift) ≈ exp(Matrix(t * A)) * B
        @test eltype(expv(t, A, B; shift=shift)) === T
    end

    @testset "expv(t::$Tt, A::SparseMatrixCSC{$TA}, B::Array{$TB,$(length(Bdims2)+1)}), tscale=$tscale, shift=$shift" for
        Tset in (Tsingle, Tdouble),
        Bdims2 in ((), (4,)),
        Tt in Tset,
        TA in Tset,
        TB in Tset,
        tscale in Tt.((0.1, 1, 10)),
        shift in (true, false)

        t = tscale * randn(Tt)
        A = VERSION ≥ v"1.1" ? sprandn(TA, n, n, 1 / n) : TA.(sprandn(n, n, 1 / n))
        B = randn(TB, n, Bdims2...)
        T = Base.promote_eltype(t, A, B)
        rT = real(T)
        @inferred expv(t, A, B; shift=shift)
        @test expv(t, A, B; shift=shift) ≈ exp(Matrix(t * A)) * B
        @test eltype(expv(t, A, B; shift=shift)) === T
    end

    @testset "errors if tolerance too low" begin
        t = randn()
        A = randn(10, 10)
        B = randn(10, 10)

        @test_throws DomainError expv(t, A, B; tol=eps(Float64) / 100)
        @test_throws DomainError expv(big(t), big.(A), big.(B))
    end
end

ext: jl


text:
# Implementation of Smith-Waterman algorithm for Local Alignment of strings
# Nov 21, 2017
# Match = 1; mis-match = -1; gop = {gap opening penalty};
# gep = {gap extension penalty}
# Modified from sw_align.jl to adapt for processing of SRS record strings

# Create scoring matrix and path matrix
function sw_align(seq_a::String, seq_b::String)
    seq_a = "^$seq_a"
    seq_b = "^$seq_b"
    match    = 1.0
    mismatch = -1.0
    gop      = -1.0
    #gep      = -0.5

    col      = length(seq_a)
    row      = length(seq_b)
    matrix   = fill(0.0, (row,col))
    path     = fill("N",(row,col))

    pathvals = ["-","|","M"]
    for i = 2:row
        for j = 2:col
            scores = []
            push!(scores,matrix[i,j-1] + gop)
            push!(scores,matrix[i-1,j] + gop)
            if seq_a[j] == seq_b[i]
                push!(scores, matrix[i-1,j-1] + match)
            else
                push!(scores, matrix[i-1,j-1] + mismatch)
            end
            val,ind = findmax(scores)
            if val < 0
                matrix[i,j] = 0
                #path[i,j] = pathvals[ind]
            else
                matrix[i,j] = val
                path[i,j] = pathvals[ind]
            end

        end
    end
    return matrix, path
end

# Traceback
function traceback(matrix::Array{Float64,2},path::Array{String,2},seq_a::String,seq_b::String)
    seq_a = "^$seq_a"
    seq_b = "^$seq_b"

    irow = size(matrix)[1]
    jcol = size(matrix)[2]

    align1 = []
    align2 = []

    maxval = 0
    ind_x  = 1
    ind_y  = 1
    for i=1:irow
        for j=1:jcol
            if matrix[i,j] > maxval
                maxval = matrix[i,j]
                ind_x  = i
                ind_y  = j
            end
        end
    end
    irow = ind_x
    jcol = ind_y

    check = false
    i = 0
    while !check
        if irow == 1 && jcol == 1
            #push!(align1,seq_a[jcol])
            #push!(align2,seq_b[irow])
            break
        elseif path[irow,jcol] == "M"
            push!(align1,seq_a[jcol])
            push!(align2,seq_b[irow])
            irow -= 1
            jcol -= 1
        elseif path[irow,jcol] == "-"
            push!(align1,seq_a[jcol])
            push!(align2,"-")
            jcol -= 1
            i = i + 1
            #println("\+\+")
        elseif path[irow,jcol] == "|"
            push!(align1,"-")
            push!(align2,seq_b[irow])
            irow -= 1
            #println("\-\-")
        elseif path[irow,jcol] == "N"
            break
        end
    end
    if maxval > 1.0
        if align1[1] == align2[1] && align1[1] == ' '
            deleteat!(align1,1)
            deleteat!(align2,1)
            ind_x = ind_x - 1
            #println("\#")
        end
        if align1[end] == align2[end] && align1[end] == ' '
            pop!(align1)
            pop!(align2)
            #ind_x = ind_x + 1
            #println("\$")
        end
        m = collect((m.match for m = eachmatch(r"\-", join(align2, ""))))
        return align1[end:-1:1], align2[end:-1:1], maxval, ind_x-1, ind_y, seq_b[ind_x-length(align2)+i+1:ind_x]
    else
        m = collect((m.match for m = eachmatch(r"\-", join(align2, ""))))
        return align1[end:-1:1], align2[end:-1:1], maxval, ind_x-1, ind_y, seq_b[ind_x-length(align2)+i+1:ind_x]
    end
end

#seq_a = "Thitmine"
#seq_b = "Thiamine is an amino acid"
#seq_b = "thiamine"

#seq_a = uppercase(seq_a)
#seq_b = uppercase(seq_b)
#matrix, path = sw_align(seq_a,seq_b)
#algn_a, algn_b, scr, indx, indy, mstr = traceback(matrix,path,seq_a,seq_b)

max_stars_repo_name: bcbi/SolrPlantAPI

max_issues_repo_name: bcbi/SolrPlantAPI

max_forks_repo_name: bcbi/SolrPlantAPI


text:
"""
    pca_reduction(wv::WordVectors, rdim=7, outdim=size(wv.vectors,1); [do_pca=true])

Post-processes word embeddings `wv` by removing the first `rdim` PCA components
from the word vectors and also reduces the dimensionality to `outdim` through
a subsequent PCA transform, if `do_pca=true`.

# Arguments
  * `wv::WordVectors` the word embeddings
  * `rdim::Int` the number of PCA components to remove from the data
     (default 7)
  * `outdim::Int` the output dimensionality of the data after the PCA
     dimensionality reduction; it is performed only if `do_pca=true`
     and the default value is the same as that of the input embeddings
     i.e. no reduction

# Keyword arguments
  * `do_pca::Bool` whether to perform a PCA transform of the
     post-processed data (default `true`)

# References:
  * [Vikas Raunak "Simple and effective dimensionality reduction for
     word embeddings", NIPS 2017 Workshop](https://arxiv.org/abs/1708.03629)
"""
function pca_reduction(wv::WordVectors{S,T,H},
                       rdim::Int=7,
                       outdim::Int=size(wv.vectors,1);
                       do_pca::Bool=true
                      ) where {S<:AbstractString, T<:Real, H<:Integer}

    # Perform first post-processing
    X = _pca_postprocessing(wv.vectors, rdim)

    # Do PCA and post-process again
    if do_pca
        outdim = clamp(outdim, 1, size(X,1))
        pratio = ifelse(size(wv.vectors,1)==outdim, 1.0, 0.99)
        M = fit(PCA, X, maxoutdim=outdim, pratio=pratio)
        X = transform(M, X)
        X = _pca_postprocessing(X, rdim)
    end

    return WordVectors{S,T,H}(wv.vocab, X, wv.vocab_hash)
end

function _pca_postprocessing(X::AbstractMatrix{T}, rdim::Int=7) where {T<:AbstractFloat}
    # Subtract the mean
    @debug "Subtracting the mean..."
    X .-= mean(X, dims=2)

    # Compute the first d PCA components
    @debug "Computing PCA components..."
    m, n = size(X)
    rdim = clamp(rdim, 1, m)
    M = fit(PCA, X, pratio=1.0, mean=0)
    M = __handle_pca_dimensions(M, m)
    Xd = transform(M, X)
    Xdv = [Xd[:,i]*Xd[:,i]' for i in 1:rdim]

    # Eliminate top d components
    @debug "Eliminating the top $rdim components..."
    Xout = similar(X)
    @inbounds @simd for i in 1:n
        Xout[:,i] = X[:,i] .- mapreduce(x->x*X[:,i], +, Xdv)
    end
    return Xout
end

# Introduces 0 components into PCA transform to force the number
# of components be equal to the number of dimensions explicitly
# specified
function __handle_pca_dimensions(M::PCA{T}, m) where {T}
    pcadim = length(M.prinvars)
    if pcadim < m
        proj = zeros(T, m, m)
        proj[:, 1:pcadim] .+= M.proj
        prinvars = zeros(T, m)
        prinvars[1:pcadim] .+= M.prinvars
        M = PCA(M.mean, proj, prinvars, M.tprinvar, M.tvar)
    end
    return M
end

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã ``2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Linearized Derivation tree (aka parse tree or syntax tree).  Wraps DerivationTrees.
Default is DFS order.
"""
module LinearDerivTrees

export Grammar
export LDTParams, LDTActions, LinearDerivTree, AbstractTreeOrder, DepthFirst, BreadthFirst
export initialize!, step!, play!, actionspace, isdone, iscomplete, get_derivtree, get_expr
export rand_with_retry!, maxlength, current_opennode

using Reexport
using ExprSearch
@reexport using DerivationTrees
using GrammaticalEvolution
using DataStructures
using RLESUtils, MemPools
using JLD

import GrammaticalEvolution.Grammar
import DerivationTrees: initialize!, actionspace, iscomplete, get_expr, pretty_string, get_sym
import ExprSearch: get_derivtree
import Base: length, push!, convert, rand!, getindex, empty!, copy!
import AbstractTrees.print_tree

const STACKSIZE = 50 #initial stack allocation size, default value 1024 allocates too much unused memory

abstract AbstractTreeOrder

immutable DepthFirst <: AbstractTreeOrder end
immutable BreadthFirst <: AbstractTreeOrder end

immutable LDTParams
    grammar::Grammar
    maxsteps::Int64
    order::AbstractTreeOrder
end
LDTParams(grammar::Grammar, maxsteps::Int64) = LDTParams(grammar, maxsteps, DepthFirst())

type LDTActions
    actions::Vector{Int64}
end
LDTActions() = LDTActions(Int64[])
LDTActions(ind::ExampleIndividual) = LDTActions(ind.genome+1)

type LinearDerivTree
    params::LDTParams
    derivtree::DerivationTree
    opennodes::Stack
    actions::LDTActions
end

function LinearDerivTree(params::LDTParams; nodepool::Union{Void,MemPool}=nothing)
    dt_params = DerivTreeParams(params.grammar)
    derivtree = (nodepool == nothing) ?
        DerivationTree(dt_params) :
        DerivationTree(dt_params, nodepool)
    stack = Stack(DerivTreeNode, STACKSIZE)
    actions = LDTActions()
    LinearDerivTree(params, derivtree, stack, actions)
end

type LinearDerivTreeSerial
    params::LDTParams
    derivtree::DerivationTree
    actions::LDTActions
end

function JLD.writeas(t::LinearDerivTree)
    LinearDerivTreeSerial(t.params, t.derivtree, t.actions)
end

function convert(::Type{LinearDerivTree}, t::LinearDerivTreeSerial)
    stack = Stack(DerivTreeNode, STACKSIZE)
    LinearDerivTree(t.params, t.derivtree, stack, t.actions)
end

push!(actions::LDTActions, a::Int64) = push!(actions.actions, a)
length(actions::LDTActions) = length(actions.actions)
getindex(actions::LDTActions, ids) = getindex(actions.actions, ids)

function copy!(dst::LDTActions, src::LDTActions)
    resize!(dst.actions, length(src.actions))
    copy!(dst.actions, src.actions)
end

convert(::Type{Array}, a::LDTActions) = copy(a.actions)
empty!(a::LDTActions) = empty!(a.actions)

get_expr(tree::LinearDerivTree) = get_expr(tree.derivtree)
empty!(stack::Stack) = empty!(stack.store)

function copy!(dst::LinearDerivTree, src::LinearDerivTree)
    dst.params = src.params
    copy!(dst.derivtree, src.derivtree)
    dst.opennodes = deepcopy(src.opennodes)
    copy!(dst.actions, src.actions)
    dst
end

get_derivtree(tree::LinearDerivTree) = tree.derivtree

function pretty_string(tree::LinearDerivTree, fmt::Format, capitalize::Bool=false)
    pretty_string(tree.derivtree, fmt, capitalize)
end

"""
Peek at top of opennode stack without popping
"""
current_opennode(tree::LinearDerivTree) = first(tree.opennodes)

function initialize!(tree::LinearDerivTree)
  root = initialize!(tree.derivtree)
  empty!(tree.opennodes)
  push!(tree.opennodes, root)
  empty!(tree.actions)
  process_non_decisions!(tree)
end

function step!(tree::LinearDerivTree, a::Int64)
  opennodes = tree.opennodes
  if isempty(opennodes)
    return #we're done
  end
  as = actionspace(tree)
  a = ((a-1) % length(as)) + 1 #1-indexed
  push!(tree.actions, a)
  node = pop!(opennodes)
  children = expand_node!(tree.derivtree, node, a)
  add_nodes!(tree.params.order, opennodes, children)
  process_non_decisions!(tree)
end

play!(tree::LinearDerivTree, x) = play!(tree, LDTActions(x))

function play!(tree::LinearDerivTree, actions::LDTActions)
  initialize!(tree)
  n = 0
  while !isdone(tree) && n < length(actions)
    n += 1
    step!(tree, actions[n])
  end
  return n #number of actions actually used
end

#iterator functions
#start(tree::LinearDerivTree) = start(tree.root)
#next(tree::LinearDerivTree, stack::Stack) = next(tree.root, stack)
#done(tree::LinearDerivTree, stack::Stack) = done(tree.root, stack)
#
##iterate through nodes in depth-first search order
#function start(startnode::DerivTreeNode)
  #stack = Stack(DerivTreeNode, STACKSIZE)
  #push!(stack, startnode)
  #return stack
#end
#
#function next(startnode::DerivTreeNode, stack::Stack)
  #node = pop!(stack)
  #for i = length(node.children):-1:1
    #push!(stack, node.children[i])
  #end
  #return node, stack
#end
#
#done(startnode::DerivTreeNode, stack::Stack) = isempty(stack)

function iscomplete(tree::LinearDerivTree)
    @assert length(tree.opennodes) == tree.derivtree.nopen #sanity check only
    @assert isempty(tree.opennodes) == iscomplete(tree.derivtree) #sanity check only
    iscomplete(tree.derivtree)
end

function isdone(tree::LinearDerivTree)
  return iscomplete(tree) || length(tree.actions) >= tree.params.maxsteps
end

function process_non_decisions!(tree::LinearDerivTree)
  opennodes = tree.opennodes
  while !isempty(opennodes) && !isa(top(opennodes).rule, DecisionRule)
    node = pop!(opennodes)
    children = expand_node!(tree.derivtree, node)
    add_nodes!(tree.params.order, opennodes, children)
  end
end

function add_nodes!(::DepthFirst, opennodes::Stack, nodes_to_add::Vector{DerivTreeNode})
    for i = length(nodes_to_add):-1:1
        push!(opennodes, nodes_to_add[i])
    end
end

function add_nodes!(::BreadthFirst, opennodes::Stack, nodes_to_add::Vector{DerivTreeNode})
    for i = length(nodes_to_add):-1:1
        push!(opennodes, nodes_to_add[i])
    end
end

"""
Range of valid action numbers at next open node
"""
function actionspace(tree::LinearDerivTree)
  if isempty(tree.opennodes)
    return 0:-1 #special -1 action for terminal state
  end
  node = top(tree.opennodes)
  return actionspace(node)
end

"""
Generate a random tree.
Returns true if completed successfully, false otherwise
"""
function rand!(tree::LinearDerivTree)
    initialize!(tree)
    while !isdone(tree)
        as = actionspace(tree)
        step!(tree, rand(as))
    end
    return iscomplete(tree)
end
"""
Generate a random tree up to a maximum number of retries to get a complete tree.
"""
function rand_with_retry!(tree::LinearDerivTree, retries::Int64=5)
    while retries > 0
        compl = rand!(tree)
        if compl
            return true
        end
        retries -= 1
    end
    return false
end

get_sym(tree::LinearDerivTree) = isempty(tree.opennodes) ? :() : get_sym(top(tree.opennodes))

#interface to AbstractTrees.jl
function print_tree(io::IO, tree::LinearDerivTree, args...; kwargs...)
    print_tree(io, tree.derivtree, args...; kwargs...)
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
export dropna, selectkeys, selectvalues

"""
`select(t::Table, which::Selection)`

Select all or a subset of columns, or a single column from the table.

`Selection` is a type union of many types that can select from a table. It can be:

1. `Integer` -- returns the column at this position.
2. `Symbol` -- returns the column with this name.
3. `Pair{Selection => Function}` -- selects and maps a function over the selection, returns the result.
4. `AbstractArray` -- returns the array itself. This must be the same length as the table.
5. `Tuple` of `Selection` -- returns a table containing a column for every selector in the tuple. The tuple may also contain the type `Pair{Symbol, Selection}`, which the selection a name. The most useful form of this when introducing a new column.

# Examples:

Selection with `Integer` -- returns the column at this position.

```jldoctest select
julia> tbl = table([0.01, 0.05], [2,1], [3,4], names=[:t, :x, :y], pkey=:t)
Table with 2 rows, 3 columns:
t     x  y
──────────
0.01  2  3
0.05  1  4

julia> select(tbl, 2)
2-element Array{Int64,1}:
 2
 1

```

Selection with `Symbol` -- returns the column with this name.

```jldoctest select
julia> select(tbl, :t)
2-element Array{Float64,1}:
 0.01
 0.05

```

Selection with `Pair{Selection => Function}` -- selects some columns and maps a function over it, then returns the mapped column.

```jldoctest select
julia> select(tbl, :t=>t->1/t)
2-element Array{Float64,1}:
 100.0
  20.0

```

Selection with `AbstractArray` -- returns the array itself.

```jldoctest select
julia> select(tbl, [3,4])
2-element Array{Int64,1}:
 3
 4

```
Selection with `Tuple`-- returns a table containing a column for every selector in the tuple.

```jldoctest select
julia> select(tbl, (2,1))
Table with 2 rows, 2 columns:
x  t
───────
2  0.01
1  0.05

julia> vx = select(tbl, (:x, :t)=>p->p.x/p.t)
2-element Array{Float64,1}:
 200.0
  20.0

julia> select(tbl, (:x,:t=>-))
Table with 2 rows, 2 columns:
x  t
────────
1  -0.05
2  -0.01
```

Note that since `tbl` was initialized with `t` as the primary key column, selections that retain the
key column will retain its status as a key. The same applies when multiple key columns are selected.

Selection with a custom array in the tuple will cause the name of the columns to be removed and replaced with integers.

```jldoctest select
julia> select(tbl, (:x, :t, [3,4]))
Table with 2 rows, 3 columns:
1  2     3
──────────
2  0.01  3
1  0.05  4
```

This is because the third column's name is unknown. In general if a column's name cannot be determined, then selection
returns an iterable of tuples rather than named tuples. In other words, it strips column names.

To specify a new name to a custom column, you can use `Symbol => Selection` selector.

```jldoctest select
julia> select(tbl, (:x,:t,:z=>[3,4]))
Table with 2 rows, 3 columns:
x  t     z
──────────
2  0.01  3
1  0.05  4

julia> select(tbl, (:x, :t, :minust=>:t=>-))
Table with 2 rows, 3 columns:
x  t     minust
───────────────
2  0.01  -0.01
1  0.05  -0.05

julia> select(tbl, (:x, :t, :vx=>(:x,:t)=>p->p.x/p.t))
Table with 2 rows, 3 columns:
x  t     vx
──────────────
2  0.01  200.0
1  0.05  20.0

```
"""
function Base.select(t::AbstractIndexedTable, which)
    ColDict(t)[which]
end

# optimization
@inline function Base.select(t::NextTable, which::Union{Symbol, Int})
    getfield(columns(t), which)
end

function selectkeys(x::NDSparse, which; kwargs...)
    ndsparse(rows(keys(x), which), values(x); kwargs...)
end

function selectvalues(x::NDSparse, which; presorted=true, copy=false, kwargs...)
    ndsparse(keys(x), rows(values(x), which); presorted=presorted, copy=copy, kwargs...)
end

"""
`reindex(t::Table, by[, select])`

Reindex `t` by columns selected in `by`.
Keeps columns selected by `select` as non-indexed columns.
By default all columns not mentioned in `by` are kept.

Use [`selectkeys`](@ref) to reindex and NDSparse object.

```jldoctest reindex
julia> t = table([2,1],[1,3],[4,5], names=[:x,:y,:z], pkey=(1,2))

julia> reindex(t, (:y, :z))
Table with 2 rows, 3 columns:
y  z  x
───────
1  4  2
3  5  1

julia> pkeynames(t)
(:y, :z)

julia> reindex(t, (:w=>[4,5], :z))
Table with 2 rows, 4 columns:
w  z  x  y
──────────
4  5  1  3
5  4  2  1

julia> pkeynames(t)
(:w, :z)

```
"""
function reindex end

function reindex(T::Type, t, by, select; kwargs...)
    if isa(by, SpecialSelector)
        return reindex(T, t, lowerselection(t, by), select; kwargs...)
    end
    if !isa(by, Tuple)
        return reindex(T, t, (by,), select; kwargs...)
    end
    if T <: NextTable && !isa(select, Tuple) && !isa(select, SpecialSelector)
        return reindex(T, t, by, (select,); kwargs...)
    end
    perm = sortpermby(t, by)
    if isa(perm, Base.OneTo)
        convert(T, rows(t, by), rows(t, select); presorted=true, copy=false, kwargs...)
    else
        convert(T, rows(t, by)[perm], rows(t, select)[perm]; presorted=true, copy=false, kwargs...)
    end
end

function reindex(t::NextTable, by=pkeynames(t), select=excludecols(t, by); kwargs...)
    reindex(collectiontype(t), t, by, select; kwargs...)
end

function reindex(t::NDSparse, by=pkeynames(t), select=valuenames(t); kwargs...)
    reindex(collectiontype(t), t, by, select; kwargs...)
end

canonname(t, x::Symbol) = x
canonname(t, x::Int) = colnames(t)[colindex(t, x)]

"""
`map(f, t::Table; select)`

Apply `f` to every row in `t`. `select` selects fields
passed to `f`.

Returns a new table if `f` returns a tuple or named tuple.
If not, returns a vector.

# Examples

```jldoctest map
julia> t = table([0.01, 0.05], [1,2], [3,4], names=[:t, :x, :y])
Table with 2 rows, 3 columns:
t     x  y
──────────
0.01  1  3
0.05  2  4

julia> manh = map(row->row.x + row.y, t)
2-element Array{Int64,1}:
 4
 6

julia> polar = map(p->@NT(r=hypot(p.x + p.y), θ=atan2(p.y, p.x)), t)
Table with 2 rows, 2 columns:
r    θ
────────────
4.0  1.24905
6.0  1.10715

```

`select` argument selects a subset of columns while iterating.

```jldoctest map

julia> vx = map(row->row.x/row.t, t, select=(:t,:x)) # row only cotains t and x
2-element Array{Float64,1}:
 100.0
  40.0

julia> map(sin, polar, select=:θ)
2-element Array{Float64,1}:
 0.948683
 0.894427

```
"""
function map(f, t::AbstractIndexedTable; select=nothing) end

function map(f, t::Dataset; select=nothing, copy=false, kwargs...)
    if isa(f, Tup) && select===nothing
        select = colnames(t)
    elseif select === nothing
        select = valuenames(t)
    end

    x = map_rows(f, rows(t, select))
    isa(x, Columns) ? table(x; copy=false, kwargs...) : x
end

function _nonna(t::Union{Columns, NextTable}, by=(colnames(t)...))
    indxs = [1:length(t);]
    if !isa(by, Tuple)
        by = (by,)
    end
    bycols = columns(t, by)
    d = ColDict(t)
    for (key, c) in zip(by, bycols)
        x = rows(t, c)
        filt_by_col!(!isnull, x, indxs)
        if isa(x, Array{<:DataValue})
            y = Array{eltype(eltype(x))}(length(x))
            y[indxs] = map(get, x[indxs])
            x = y
        elseif isa(x, DataValueArray)
            x = x.values # unsafe unwrap
        end
        d[key] = x
    end
    (d[], indxs)
end

"""
`dropna(t[, select])`

Drop rows which contain NA values.

```jldoctest dropna
julia> t = table([0.1, 0.5, NA,0.7], [2,NA,4,5], [NA,6,NA,7],
                  names=[:t,:x,:y])
Table with 4 rows, 3 columns:
t    x    y
─────────────
0.1  2    #NA
0.5  #NA  6
#NA  4    #NA
0.7  5    7

julia> dropna(t)
Table with 1 rows, 3 columns:
t    x  y
─────────
0.7  5  7
```
Optionally `select` can be speicified to limit columns to look for NAs in.

```jldoctest dropna

julia> dropna(t, :y)
Table with 2 rows, 3 columns:
t    x    y
───────────
0.5  #NA  6
0.7  5    7

julia> t1 = dropna(t, (:t, :x))
Table with 2 rows, 3 columns:
t    x  y
───────────
0.1  2  #NA
0.7  5  7
```

Any columns whose NA rows have been dropped will be converted
to non-na array type. In our last example, columns `t` and `x`
got converted from `Array{DataValue{Int}}` to `Array{Int}`.
Similarly if the vectors are of type `DataValueArray{T}`
(default for `loadtable`) they will be converted to `Array{T}`.
```julia
julia> typeof(column(dropna(t,:x), :x))
Array{Int64,1}
```
"""
function dropna(t::Dataset, by=(colnames(t)...))
    subtable(_nonna(t, by)...)
end

filt_by_col!(f, col, indxs) = filter!(i->f(col[i]), indxs)

Base.@deprecate select(arr::NDSparse, conditions::Pair...) filter(conditions, arr)
Base.@deprecate select(arr::NDSparse, which::DimName...; agg=nothing) selectkeys(arr, which; agg=agg)

"""
`filter(pred, t::Union{NextTable, NDSparse}; select)`

Filter rows in `t` according to `pred`. `select` choses the fields that act as input to `pred`.

`pred` can be:

- A function - selected structs or values are passed to this function
- A tuple of `column => function` pairs: applies to each named column the corresponding function, keeps only rows where all such conditions are satisfied.

By default, `filter` iterates a table a row at a time:
```jldoctest filter
julia> t = table(["a","b","c"], [0.01, 0.05, 0.07], [2,1,0],
                 names=[:n, :t, :x])
Table with 3 rows, 3 columns:
n    t     x
────────────
"a"  0.01  2
"b"  0.05  1
"c"  0.07  0

julia> filter(p->p.x/p.t < 100, t) # whole row
Table with 2 rows, 3 columns:
n    t     x
────────────
"b"  0.05  1
"c"  0.07  0

```

By default, `filter` iterates by values of an `NDSparse`:

```jldoctest filter
julia> x = ndsparse(@NT(n=["a","b","c"], t=[0.01, 0.05, 0.07]), [2,1,0])
2-d NDSparse with 3 values (Int64):
n    t    │
──────────┼──
"a"  0.01 │ 2
"b"  0.05 │ 1
"c"  0.07 │ 0

julia> filter(y->y<2, x)
2-d NDSparse with 2 values (Int64):
n    t    │
──────────┼──
"b"  0.05 │ 1
"c"  0.07 │ 0
```

If select is specified. (See [Selection convention](@ref select)) then, the selected values will be iterated instead.

```jldoctest filter
julia> filter(iseven, t, select=:x)
Table with 2 rows, 3 columns:
n    t     x
────────────
"a"  0.01  2
"c"  0.07  0

julia> filter(p->p.x/p.t < 100, t, select=(:x,:t))
Table with 2 rows, 3 columns:
n    t     x
────────────
"b"  0.05  1
"c"  0.07  0
```

`select` works similarly for `NDSparse`:
```jldoctest filter
julia> filter(p->p[2]/p[1] < 100, x, select=(:t, 3))
2-d NDSparse with 2 values (Int64):
n    t    │
──────────┼──
"b"  0.05 │ 1
"c"  0.07 │ 0
```
Here 3 represents the third column, which is the values, `p` is a tuple of `t` field and the value.

Filtering by many single columns can be done by passing a tuple of `column_name => function` pairs.

```jldoctest filter
julia> filter((:x=>iseven, :t=>a->a>0.01), t)
Table with 1 rows, 3 columns:
n    t     x
────────────
"c"  0.07  0

julia> filter((3=>iseven, :t=>a->a>0.01), x) # NDSparse
2-d NDSparse with 1 values (Int64):
n    t    │
──────────┼──
"c"  0.07 │ 0

```

"""
function Base.filter(fn, t::Dataset; select=valuenames(t))
    x = rows(t, select)
    indxs = find(map(fn, x))
    subtable(t, indxs, presorted=true)
end

function Base.filter(pred::Tuple, t::Dataset; select=nothing)
    indxs = [1:length(t);]
    x = select === nothing ? t : rows(t, select)
    for p in pred
        if isa(p, Pair)
            c, f = p
            filt_by_col!(f, rows(x, c), indxs)
        else
            filt_by_col!(p, x, indxs)
        end
    end
    subtable(t, indxs, presorted=true)
end

function Base.filter(pred::Pair, t::Dataset; select=nothing)
    filter((pred,), t, select=select)
end

# We discard names of fields in a named tuple. keeps it consistent
# with map and reduce, we don't select using those
function Base.filter(pred::NamedTuple, t::Dataset; select=nothing)
    filter(astuple(pred), t, select=select)
end


text:
# A few supported manifolds

export
  # general manifolds
  Euclidean,
  Circle,
  ℝ,
  Euclid,
  Euclid2,
  Euclid3,
  Euclid4,
  SE2_Manifold,
  SE3_Manifold,

  # special quirks that still need to be fixed
  Circle1
  # Circular,

export
  coords,
  uncoords,
  getPointsManifold

#

const _AMP_CIRCLE = Manifolds.Circle()

# this is just wrong and needs to be fixed
const Circle1 = Circle()

const Euclid =  TranslationGroup(1)
const EuclideanManifold = Euclid

const Euclid2 = TranslationGroup(2)
const Euclid3 = TranslationGroup(3)
const Euclid4 = TranslationGroup(4)

# TODO if not easy simplification exists, then just deprecate this
const SE2_Manifold = SpecialEuclidean(2)
const SE3_Manifold = SpecialEuclidean(3)

#


text:
using Test
using Bijectors, ForwardDiff, LinearAlgebra
using Random: seed!

seed!(1)

@testset "InvertibleBatchNorm" begin
    x = randn(2, 20)
    bn = InvertibleBatchNorm(2)

    @test inv(inv(bn)) == bn
    @test inv(bn)(bn(x)) ≈ x
    @test (inv(bn) ∘ bn)(x) ≈ x
    @test_throws ErrorException forward(bn, randn(10,2))
    @test logabsdetjac(inv(bn), bn(x)) ≈ - logabsdetjac(bn, x)

    y, ladj = forward(bn, x)
    @test log(abs(det(ForwardDiff.jacobian(bn, x)))) ≈ sum(ladj)
    @test log(abs(det(ForwardDiff.jacobian(inv(bn), y)))) ≈ sum(logabsdetjac(inv(bn), y))

    test_functor(bn, (b = bn.b, logs = bn.logs))
end

@testset "PlanarLayer" begin
    for i in 1:4
        flow = PlanarLayer(2)
        z = randn(2, 20)
        forward_diff = log(abs(det(ForwardDiff.jacobian(t -> flow(t), z))))
        our_method = sum(forward(flow, z).logabsdetjac)

        @test our_method ≈ forward_diff
        @test inv(flow)(flow(z)) ≈ z
        @test (inv(flow) ∘ flow)(z) ≈ z
    end

    w = ones(10)
    u = zeros(10)
    b = 1.0
    flow = PlanarLayer(w, u, b)
    z = ones(10, 100)
    @test inv(flow)(flow(z)) ≈ z

    test_functor(flow, (w = w, u = u, b = b))
    test_functor(inv(flow), (orig = flow,))

    @testset "find_alpha" begin
        for wt_y in (-20.3, -3, -3//2, 0.0, 5, 29//4, 12.3)
            # the root finding algorithm assumes wt_u_hat ≥ -1 (satisfied for the flow)
            # |wt_u_hat| < eps checks that empty brackets are handled correctly
            # https://github.com/TuringLang/Bijectors.jl/issues/204
            for wt_u_hat in (-1, -1//2, -1e-20, 0, 1e-20, 3, 11//3, 17.2)
                for b in (-19.3, -8//3, -1, 0.0, 1//2, 3, 4.3)
                    # find α that solves wt_y = α + wt_u_hat * tanh(α + b)
                    α = @inferred(Bijectors.find_alpha(wt_y, wt_u_hat, b))

                    # check if α is an approximate solution to the considered equation
                    # have to set atol if wt_y is zero (otherwise only equality is checked)
                    @test wt_y ≈ α + wt_u_hat * tanh(α + b) atol=iszero(wt_y) ? 1e-14 : 0.0
                end
            end
        end

        # floating point issues
        # https://github.com/TuringLang/Bijectors.jl/issues/204
        wt_y = 0.8845640339582252
        wt_u_hat = 0.8296950433716855
        b = -1e8
        @test Bijectors.find_alpha(wt_y, wt_u_hat, b) ≈ wt_y + wt_u_hat
    end
end

@testset "RadialLayer" begin
    for i in 1:4
        flow = RadialLayer(2)
        z = randn(2, 20)
        forward_diff = log(abs(det(ForwardDiff.jacobian(t -> flow(t), z))))
        our_method = sum(forward(flow, z).logabsdetjac)

        @test our_method ≈ forward_diff
        @test inv(flow)(flow(z)) ≈ z rtol=0.2
        @test (inv(flow) ∘ flow)(z) ≈ z rtol=0.2
    end

    α_ = 1.0
    β = 1.0
    z_0 = zeros(10)
    z = ones(10, 100)
    flow = RadialLayer(α_, β, z_0)
    @test inv(flow)(flow(z)) ≈ z

    test_functor(flow, (α_ = α_, β = β, z_0 = z_0))
    test_functor(inv(flow), (orig = flow,))
end

@testset "Flows" begin
    d = MvNormal(zeros(2), I)
    b = PlanarLayer(2)
    flow = transformed(d, b)  # <= Radial flow

    y = rand(flow)
    @test logpdf(flow, y) != 0.0

    x = rand(d)
    y = flow.transform(x)
    res = forward(flow.transform, x)
    lp = logpdf_forward(flow, x, res.logabsdetjac)

    @test res.rv ≈ y
    @test logpdf(flow, y) ≈ lp rtol=0.1

    # flow with unconstrained-to-constrained
    d1 = Beta()
    b1 = inv(bijector(d1))
    d2 = InverseGamma()
    b2 = inv(bijector(d2))

    x = rand(d) .+ 10
    y = b(x)

    sb = stack(b1, b1)
    @test all((sb ∘ b)(x) .≤ 1.0)

    sb = stack(b1, b2)
    cb = (sb ∘ b)
    y = cb(x)
    @test (0 ≤ y[1] ≤ 1.0) && (0 < y[2])
end


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

using Test
using Base.Threads

@test nthreadpools() == 2
@test threadpool() == :default
@test threadpool(2) == :interactive
dtask() = @test threadpool(current_task()) == :default
itask() = @test threadpool(current_task()) == :interactive
dt1 = @spawn dtask()
dt2 = @spawn :default dtask()
it = @spawn :interactive itask()
wait(dt1)
wait(dt2)
wait(it)

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
const TIVN = Tuple{Int, Union{Vector, Nothing}}

function read_input(input::Vector)::Matrix{Int}
  n = length(input)
  dist = fill(zero(Int), n, n)

  for (rix, row) ∈ input |> enumerate
    for (cjx, col) ∈ split(row, r"\s+") |> enumerate
      dist[rix, cjx] = parse(Int, col)
    end
  end

  dist
end

function max_path(d::Matrix{Int}; show_path=false)::TIVN
  """
  Start from penultimate row downto 1st row using max. from
  2 possible cells (from the row below)
  """
  for prow_ix ∈ size(d)[1] - 1:-1:1   # penultimate row index
    for cjx ∈ 1:prow_ix
      d[prow_ix, cjx] += max(d[prow_ix + 1, cjx], d[prow_ix + 1, cjx + 1])
    end
  end

  show_path ? (d[1, 1], find_path(d)) : (d[1, 1], nothing)
end

function find_path(d::Matrix{Int})::Vector
  path = [(1, 1)]

  for row ∈ 2:size(d)[1]
    push!(path, coord_max(d, path[end]))
  end

  path
end

function coord_max(d, coord)
  (x, y) = coord
  d[x + 1, y] > d[x + 1, y + 1] ? (x + 1, y) : (x + 1, y + 1)
end

function calc_max_path(input::AbstractString; show_path=false)::TIVN
  vstr = if occursin(".txt", input)
    # file
    open(input) do f
      readlines(f)
    end
  else
    split(input, "\n")
  end

  vstr .|>
    strip |>
    a -> filter(s -> length(s) > 0, a) .|>
    String |>
    read_input |>
    d -> max_path(d; show_path)
end

max_stars_repo_name: pascal-p/Euler

max_issues_repo_name: pascal-p/Euler

max_forks_repo_name: pascal-p/Euler


text:
module SAIMFitKit

    export  fit_SAIM,
            plot_fit,
            #Types
            SAIMOptics

    # using statement for external packages -
    using Plots
    using LsqFit
    using Images, FileIO
    using JLD
    using Statistics
    using CSV
    using DataFrames

    # include package codes -
    include("Types.jl")
    include("DataFit.jl")
    include("OneColor.jl")
    include("TwoColor.jl")
    include("ThreeColor.jl")
    include("Utilities.jl")

end # module


text:
mutable struct PrimitiveVars
    density::Float64
    pressure::Float64
    velocity::Float64
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
@recipe function f(p::VDPTagProblem)
    m = mdp(p)
    ratio --> :equal
    xlim --> (-5, 5)
    ylim --> (-5, 5)
    bs = m.barriers
    # @series begin
    #     lim = (-3.5,3.5)
    #     pts = linspace( -3, 3, 20)
    #     xys = [Vec2(x, y) for x in pts, y in pts]
    #     xs = [xy[1] for xy in xys]
    #     ys = [xy[2] for xy in xys]
    #     seriestype := quiver
    #     label := ""
    #     color --> :lightblue
    #     quiver := (x,y)->0.1*vdp_dynamics(m.mu, Vec2(x,y)),
    #     xs, ys
    # end
    if bs isa CardinalBarriers
        for dir in cardinals()
            ends = (bs.start*dir, (bs.start+bs.len)*dir)
            color := :black
            linewidth --> 4
            label --> ""
            @series [v[1] for v in ends], [v[2] for v in ends]
        end
    end
end

@recipe function f(pomdp::VDPTagPOMDP, h::POMDPHistory{TagState})
    ratio --> :equal
    xlim --> (-5, 5)
    ylim --> (-5, 5)
    @series mdp(pomdp), h
    @series begin
        label := "belief"
        belief_hist(h)[end]
    end
end

@recipe function f(p::VDPTagProblem, h::SimHistory)
    m = mdp(p)
    ratio --> :equal
    xlim --> (-5, 5)
    ylim --> (-5, 5)
    @series begin
        label := "path"
        x = [s.agent[1] for s in state_hist(h)[1:end-1]]
        y = [s.agent[2] for s in state_hist(h)[1:end-1]]
        x, y
    end
    @series begin
        a = action_hist(h)[end]
        if a isa TagAction && a.look
            color := :blue
        else
            color := :red
        end
        s = state_hist(h)[end-1]
        label := "current agent position"
        pts = Plots.partialcircle(0, 2*pi, 100, m.tag_radius)
        x, y = Plots.unzip(pts)
        x+s.agent[1], y+s.agent[2]
    end
    @series begin
        seriestype := :scatter
        label := "current target"
        pos = state_hist(h)[end-1].target
        color --> :orange
        [pos[1]], [pos[2]]
    end
    @series begin m end
end

@recipe function f(pc::ParticleCollection{TagState})
    seriestype := :scatter
    x = [p.target[1] for p in particles(pc)]
    y = [p.target[2] for p in particles(pc)]
    color --> :black
    markersize --> 0.1
    x, y
end

"Create a gif of a history and return the filename."
function gif(p::VDPTagProblem, h::SimHistory)

end

"Create a quiver plot of the equations and the barriers"
function Plots.quiver(p::VDPTagProblem)
    m = mdp(p)
    lim = (-3.5,3.5)
    pts = linspace( -3, 3, 16)
    xys = [Vec2(x, y) for x in pts, y in pts]
    xs = [xy[1] for xy in xys]
    ys = [xy[2] for xy in xys]
    quiver(xs, ys,
            quiver = (x,y)->0.1*vdp_dynamics(m.mu, Vec2(x,y)),
            color=:lightblue,
          )
    plot!(p) # to get barriers
    plot!(xlim=lim, ylim=lim)
end


text:
"""
    function cudaAdapt_(code)

Function that adapts vector declarations and math functions to CUDA.
"""
function cudaAdapt_(code::Expr)

    #Subs all math
    l = [i for i in names(Base.Math) if isdefined(CUDA,i) && i != :(^)]
    for i in l
        code = subs_(code,i,Meta.parse(string("CUDA.","$i")))
    end
    #Subs vectors
    for i in [:zeros,:ones]
        code = subs_(code,i,Meta.parse(string("CUDA.","$i")))
    end
    code = subs_(code,:Array,Meta.parse(string("CUDA.CuArray")))
    code = subs_(code,:Int,Meta.parse(string("Int32")))
    code = subs_(code,:Float64,Meta.parse(string("Float32")))

    return code
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
"""Theorems that we can prove with this package"""
module Theorems

include("Napoleon.jl")
include("Centroid.jl")
include("Ceva.jl")

end


text:
# Autogenerated wrapper script for marble_jll for x86_64-linux-gnu-cxx03
export astro, marble_qt, marblewidget

using Qt_jll
using Libglvnd_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("marble")
JLLWrappers.@declare_library_product(astro, "libastro.so.1")
JLLWrappers.@declare_executable_product(marble_qt)
JLLWrappers.@declare_library_product(marblewidget, "libmarblewidget-qt5.so.28")
function __init__()
    JLLWrappers.@generate_init_header(Qt_jll, Libglvnd_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        astro,
        "lib/libastro.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        marble_qt,
        "bin/marble-qt",
    )

    JLLWrappers.@init_library_product(
        marblewidget,
        "lib/libmarblewidget-qt5.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test
using DelimitedFiles

@testset "Find neighbors" begin

    include("../src/init.jl")

    @testset "Neighbors inside domain" begin
        actual = neighbors(4, 12, 3, 4)
        expected = Dict{String,Int64}("S" => 3, "W" => 1, "N" => 5, "E" => 7)
        @test (actual == expected)
    end;

    @testset "Neighbors near boundary" begin
        actual = neighbors(11, 12, 3, 4)
        expected = Dict{String,Int64}("S" => 10, "W" => 8, "N" => -1, "E" => -1)
        @test (actual == expected)
    end;

end;

@testset "End-to-end tests" begin

    reference_8x4_output = readdlm("test/resources/reference_8x4_output.dat")

    @testset "1 process, 8x4 grid" begin
        cmd = `mpirun -np 1 julia src/diffusion.jl test/resources/test1_params.yml`
        run(cmd)
        expected = reference_8x4_output
        actual = readdlm("test/resources/test1.dat")
        @test all(actual .≈ expected)
    end;

    @testset "1 process, 4x8 grid" begin
        cmd = `mpirun -np 1 julia src/diffusion.jl test/resources/test2_params.yml`
        run(cmd)
        nx = 4
        ny = 8
        expected = transpose(reshape(reference_8x4_output, nx + 2, ny + 2))
        actual = reshape(readdlm("test/resources/test2.dat"), ny + 2, nx + 2)
        @test all(actual .≈ expected)
    end;

    @testset "2x1 processes, 8x4 grid" begin
        cmd = `mpirun -np 2 julia src/diffusion.jl test/resources/test3_params.yml`
        run(cmd)
        expected = reference_8x4_output
        actual = readdlm("test/resources/test3.dat")
        @test all(actual .≈ expected)
    end;

    @testset "2x1 processes, 4x8 grid" begin
        cmd = `mpirun -np 2 julia src/diffusion.jl test/resources/test4_params.yml`
        run(cmd)
        nx = 4
        ny = 8
        expected = transpose(reshape(reference_8x4_output, nx + 2, ny + 2))
        actual = reshape(readdlm("test/resources/test4.dat"), ny + 2, nx + 2)
        @test all(actual .≈ expected)
    end;

    @testset "1x2 processes, 8x4 grid" begin
        cmd = `mpirun -np 2 julia src/diffusion.jl test/resources/test5_params.yml`
        run(cmd)
        expected = reference_8x4_output
        actual = readdlm("test/resources/test5.dat")
        @test all(actual .≈ expected)
    end;

    @testset "1x2 processes, 4x8 grid" begin
        cmd = `mpirun -np 2 julia src/diffusion.jl test/resources/test6_params.yml`
        run(cmd)
        nx = 4
        ny = 8
        expected = transpose(reshape(reference_8x4_output, nx + 2, ny + 2))
        actual = reshape(readdlm("test/resources/test6.dat"), ny + 2, nx + 2)
        @test all(actual .≈ expected)
    end;

end;


text:
module RingPolymers

include("ring_polymer.jl")
export nbeads

end # module


text:
# Note to devs. Use GLPK or Cbc for models with linear constraints and linear cost functions
# Use OSQP for models with quadratic cost function and linear constraints and ipopt otherwise
@testset "All PowerModels models construction" begin
    networks = [
        (PM.ACPPowerModel, fast_ipopt_optimizer),
        (PM.ACRPowerModel, fast_ipopt_optimizer),
        (PM.ACTPowerModel, fast_ipopt_optimizer),
        #(PM.IVRPowerModel, fast_ipopt_optimizer), #instantiate_ivp_expr_model not implemented
        (PM.DCPPowerModel, fast_ipopt_optimizer),
        (PM.DCMPPowerModel, fast_ipopt_optimizer),
        (PM.NFAPowerModel, fast_ipopt_optimizer),
        (PM.DCPLLPowerModel, fast_ipopt_optimizer),
        (PM.LPACCPowerModel, fast_ipopt_optimizer),
        (PM.SOCWRPowerModel, fast_ipopt_optimizer),
        (PM.SOCWRConicPowerModel, scs_solver),
        (PM.QCRMPowerModel, fast_ipopt_optimizer),
        (PM.QCLSPowerModel, fast_ipopt_optimizer),
        #(PM.SOCBFPowerModel, fast_ipopt_optimizer), # not implemented
        (PM.BFAPowerModel, fast_ipopt_optimizer),
        #(PM.SOCBFConicPowerModel, fast_ipopt_optimizer), # not implemented
        (PM.SDPWRMPowerModel, scs_solver),
        (PM.SparseSDPWRMPowerModel, scs_solver),
        (PTDFPowerModel, fast_ipopt_optimizer),
    ]
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    for (network, solver) in networks
        template = get_thermal_dispatch_template_network(
            NetworkModel(network; PTDF=PSY.PTDF(c_sys5)),
        )
        ps_model = DecisionModel(template, c_sys5; optimizer=solver)
        @test build!(ps_model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        @test ps_model.internal.container.pm !== nothing
        # TODO: Change test
        # @test :nodal_balance_active in keys(ps_model.internal.container.expressions)
    end
end

@testset "Network Copper Plate" begin
    template = get_thermal_dispatch_template_network(CopperPlatePowerModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    test_results = IdDict{System, Vector{Int}}(
        c_sys5 => [120, 0, 120, 120, 24],
        c_sys14 => [120, 0, 120, 120, 24],
        c_sys14_dc => [120, 0, 120, 120, 24],
    )
    constraint_keys = [PSI.ConstraintKey(CopperPlateBalanceConstraint, PSY.System)]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 240000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )

    for (ix, sys) in enumerate(systems)
        ps_model = DecisionModel(template, sys; optimizer=HiGHS_optimizer)

        @test build!(ps_model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_keys)
        moi_tests(
            ps_model,
            false,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(ps_model, [MOI.OPTIMAL], test_obj_values[sys], 10000)
    end
    template = get_thermal_dispatch_template_network(
        NetworkModel(CopperPlatePowerModel; use_slacks=true),
    )
    ps_model_re = DecisionModel(
        template,
        PSB.build_system(PSITestSystems, "c_sys5_re");
        optimizer=GLPK_optimizer,
    )
    @test build!(ps_model_re; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    psi_checksolve_test(ps_model_re, [MOI.OPTIMAL], 240000.0, 10000)
end

@testset "Network DC-PF with PTDF Model" begin
    template = get_thermal_dispatch_template_network(StandardPTDFModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_keys = [
        PSI.ConstraintKey(RateLimitConstraint, PSY.Line, "lb"),
        PSI.ConstraintKey(RateLimitConstraint, PSY.Line, "ub"),
        PSI.ConstraintKey(CopperPlateBalanceConstraint, PSY.System),
        PSI.ConstraintKey(NetworkFlowConstraint, PSY.Line),
    ]
    PTDF_ref = IdDict{System, PTDF}(
        c_sys5 => PTDF(c_sys5),
        c_sys14 => PTDF(c_sys14),
        c_sys14_dc => PTDF(c_sys14_dc),
    )
    test_results = IdDict{System, Vector{Int}}(
        c_sys5 => [264, 0, 264, 264, 168],
        c_sys14 => [600, 0, 600, 600, 504],
        c_sys14_dc => [600, 48, 552, 552, 456],
    )
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 340000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    for (ix, sys) in enumerate(systems)
        template = get_thermal_dispatch_template_network(
            NetworkModel(StandardPTDFModel; PTDF=PTDF_ref[sys]),
        )
        ps_model = DecisionModel(template, sys; optimizer=HiGHS_optimizer)

        @test build!(ps_model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_keys)
        moi_tests(
            ps_model,
            false,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(
            ps_model,
            [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL],
            test_obj_values[sys],
            10000,
        )
    end
    # PTDF input Error testing
    ps_model = DecisionModel(template, c_sys5; optimizer=GLPK_optimizer)
    @test build!(
        ps_model;
        console_level=Logging.AboveMaxLevel,  # Ignore expected errors.
        output_dir=mktempdir(cleanup=true),
    ) == PSI.BuildStatus.FAILED
end

@testset "Sparse Network DC-PF with PTDFPowerModel" begin
    template = get_thermal_dispatch_template_network(PTDFPowerModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_keys = [
        PSI.ConstraintKey(RateLimitConstraint, PSY.Line, "lb"),
        PSI.ConstraintKey(RateLimitConstraint, PSY.Line, "ub"),
        PSI.ConstraintKey(CopperPlateBalanceConstraint, PSY.System),
        PSI.ConstraintKey(NetworkFlowConstraint, PSY.Line),
    ]
    test_results = IdDict{System, Vector{Int}}(
        c_sys5 => [264, 0, 264, 264, 168],
        c_sys14 => [600, 0, 600, 600, 504],
        c_sys14_dc => [600, 48, 552, 552, 456],
    )
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 340000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    for (ix, sys) in enumerate(systems)
        ps_model = DecisionModel(template, sys; optimizer=HiGHS_optimizer)

        @test build!(ps_model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_keys)
        moi_tests(
            ps_model,
            false,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(
            ps_model,
            [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL],
            test_obj_values[sys],
            10000,
        )
    end
end

@testset "Network DC lossless -PF network with PowerModels DCPlosslessForm" begin
    template = get_thermal_dispatch_template_network(DCPPowerModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_keys = [
        PSI.ConstraintKey(PSI.RateLimitConstraint, PSY.Line, "ub"),
        PSI.ConstraintKey(PSI.RateLimitConstraint, PSY.Line, "lb"),
        PSI.ConstraintKey(PSI.NodalBalanceActiveConstraint, PSY.Bus),
    ]
    test_results = IdDict{System, Vector{Int}}(
        c_sys5 => [384, 0, 408, 408, 288],
        c_sys14 => [936, 0, 1080, 1080, 840],
        c_sys14_dc => [984, 96, 984, 984, 840],
    )
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 342000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    for (ix, sys) in enumerate(systems)
        ps_model = DecisionModel(template, sys; optimizer=ipopt_optimizer)
        @test build!(ps_model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_keys)
        moi_tests(
            ps_model,
            false,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(
            ps_model,
            [MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
            test_obj_values[sys],
            1000,
        )
    end
end

@testset "Network Solve AC-PF PowerModels StandardACPModel" begin
    template = get_thermal_dispatch_template_network(ACPPowerModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    # Check for voltage and angle constraints
    constraint_keys = [
        PSI.ConstraintKey(RateLimitConstraintFromTo, PSY.Line),
        PSI.ConstraintKey(RateLimitConstraintToFrom, PSY.Line),
        PSI.ConstraintKey(PSI.NodalBalanceActiveConstraint, PSY.Bus),
        PSI.ConstraintKey(PSI.NodalBalanceReactiveConstraint, PSY.Bus),
    ]
    test_results = IdDict{System, Vector{Int}}(
        c_sys5 => [1056, 0, 384, 384, 264],
        c_sys14 => [2832, 0, 720, 720, 696],
        c_sys14_dc => [2832, 96, 672, 672, 744],
    )
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 340000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    for (ix, sys) in enumerate(systems)
        ps_model = DecisionModel(template, sys; optimizer=ipopt_optimizer)
        @test build!(ps_model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_keys)
        moi_tests(
            ps_model,
            false,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(
            ps_model,
            [MOI.TIME_LIMIT, MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
            test_obj_values[sys],
            10000,
        )
    end
end

# TODO: Add constraint tests for these models, other is redundant with first test
@testset "Network Solve AC-PF PowerModels NFAPowerModel" begin
    template = get_thermal_dispatch_template_network(NFAPowerModel)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    objfuncs = [GAEVF, GQEVF, GQEVF]
    constraint_keys = [PSI.ConstraintKey(PSI.NodalBalanceActiveConstraint, PSY.Bus)]
    test_results = Dict{System, Vector{Int}}(
        c_sys5 => [264, 0, 264, 264, 120],
        c_sys14 => [600, 0, 600, 600, 336],
        c_sys14_dc => [648, 96, 552, 552, 384],
    )
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 300000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    for (ix, sys) in enumerate(systems)
        ps_model = DecisionModel(template, sys; optimizer=ipopt_optimizer)
        @test build!(ps_model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_keys)
        moi_tests(
            ps_model,
            false,
            test_results[sys][1],
            test_results[sys][2],
            test_results[sys][3],
            test_results[sys][4],
            test_results[sys][5],
            false,
        )
        psi_checkobjfun_test(ps_model, objfuncs[ix])
        psi_checksolve_test(
            ps_model,
            [MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
            test_obj_values[sys],
            10000,
        )
    end
end

# TODO: Add constraint tests for these models, other is redundant with first test
@testset "Other Network AC PowerModels models" begin
    networks = [#ACPPowerModel, Already tested
        ACRPowerModel,
        ACTPowerModel,
    ]
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    # TODO: add model specific constraints to this list. Voltages, etc.
    constraint_keys = [
        PSI.ConstraintKey(PSI.NodalBalanceActiveConstraint, PSY.Bus),
        PSI.ConstraintKey(PSI.NodalBalanceReactiveConstraint, PSY.Bus),
    ]
    ACR_test_results = Dict{System, Vector{Int}}(
        c_sys5 => [1056, 0, 240, 240, 264],
        c_sys14 => [2832, 0, 240, 240, 696],
        c_sys14_dc => [2832, 96, 240, 240, 744],
    )
    ACT_test_results = Dict{System, Vector{Int}}(
        c_sys5 => [1344, 0, 384, 384, 840],
        c_sys14 => [3792, 0, 720, 720, 2616],
        c_sys14_dc => [3696, 96, 672, 672, 2472],
    )
    test_results = Dict(zip(networks, [ACR_test_results, ACT_test_results]))
    for network in networks, sys in systems
        template = get_thermal_dispatch_template_network(network)
        ps_model = DecisionModel(template, sys; optimizer=fast_ipopt_optimizer)
        @test build!(ps_model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_keys)
        moi_tests(
            ps_model,
            false,
            test_results[network][sys][1],
            test_results[network][sys][2],
            test_results[network][sys][3],
            test_results[network][sys][4],
            test_results[network][sys][5],
            false,
        )
        @test ps_model.internal.container.pm !== nothing
    end
end

# TODO: Add constraint tests for these models, other is redundant with first test
@testset "Network DC-PF PowerModels quadratic loss approximations models" begin
    networks = [DCPLLPowerModel, LPACCPowerModel]
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    c_sys14_dc = PSB.build_system(PSITestSystems, "c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    # TODO: add model specific constraints to this list. Bi-directional flows etc
    constraint_keys = [PSI.ConstraintKey(PSI.NodalBalanceActiveConstraint, PSY.Bus)]
    test_obj_values = IdDict{System, Float64}(
        c_sys5 => 340000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )
    DCPLL_test_results = Dict{System, Vector{Int}}(
        c_sys5 => [528, 0, 408, 408, 288],
        c_sys14 => [1416, 0, 1080, 1080, 840],
        c_sys14_dc => [1416, 96, 984, 984, 840],
    )
    LPACC_test_results = Dict{System, Vector{Int}}(
        c_sys5 => [1200, 0, 384, 384, 840],
        c_sys14 => [3312, 0, 720, 720, 2616],
        c_sys14_dc => [3264, 96, 672, 672, 2472],
    )
    test_results = Dict(zip(networks, [DCPLL_test_results, LPACC_test_results]))
    for network in networks, (ix, sys) in enumerate(systems)
        template = get_thermal_dispatch_template_network(network)
        ps_model = DecisionModel(template, sys; optimizer=ipopt_optimizer)
        @test build!(ps_model; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        psi_constraint_test(ps_model, constraint_keys)
        moi_tests(
            ps_model,
            false,
            test_results[network][sys][1],
            test_results[network][sys][2],
            test_results[network][sys][3],
            test_results[network][sys][4],
            test_results[network][sys][5],
            false,
        )
        @test ps_model.internal.container.pm !== nothing
        psi_checksolve_test(
            ps_model,
            [MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
            test_obj_values[sys],
            10000,
        )
    end
end

@testset "Network Unsupported Power Model Formulations" begin
    for network in PSI.UNSUPPORTED_POWERMODELS
        template = get_thermal_dispatch_template_network(network)
        ps_model = DecisionModel(
            template,
            PSB.build_system(PSITestSystems, "c_sys5");
            optimizer=ipopt_optimizer,
        )
        @test build!(
            ps_model;
            console_level=Logging.AboveMaxLevel,  # Ignore expected errors.
            output_dir=mktempdir(cleanup=true),
        ) == PSI.BuildStatus.FAILED
    end
end


text:
reduce_kvecs(vec::Vector,tol) = vec

function reduce_kvecs(vecs::Vector{Vector{T}},tol::T) where T<:AbstractFloat
    all_inds = collect(eachindex(vecs))
    vecs = map(vecs) do w
        ind_ins = findall([norm(v - w) < tol for v in vecs[all_inds]])
        inds = all_inds[ind_ins]
        deleteat!(all_inds,ind_ins)
        isempty(inds) ? [zero(T),-one(T)] :  mean(vecs[inds])
    end
    vecs = deleteat!(vecs, findall([w[2] < -sqrt(tol) for w in vecs]))

    return vecs
end


text: accuracy(pred,y) = mean(onecold(pred |>  cpu).==onecold(y |> cpu))


text:
initial_temperature(model::SaturationModel,p) = 400.0

function temperature_impl(mt::SingleSatP,model::T,p) where T<:SaturationModel
    f0(t) = pressure_impl(QuickStates.sat_t(),model,t)-p
    df0(t) = ForwardDiff.derivative(f0,t)
    t0 = initial_temperature(model,p)
    t = Roots.find_zero(f0, t0)
    return t
end

function temperature_impl(mt::SingleΦP,model::SaturationModel,p)
    f0(t) = temperature_impl(QuickStates.sat_t(),model,t)-p
    df0(t) = ForwardDiff.derivative(f0,t)
    t0 = initial_temperature(model,p)
    t = Roots.find_zero((f0, df0), t0, Roots.Newton())
    return t
end

function pressure_impl(mt::SingleΦT,model::SaturationModel,p)
    return pressure_impl(QuickStates.sat_t(),model,t)
end

function pressure(model::SaturationModel,st::ThermodynamicState,unit=u"Pa")
    return pressure(state_type(st),model,st,unit)
end
function pressure(mt::SingleSatT,model::SaturationModel,st::ThermodynamicState,unit)
    t = temperature(FromState(),st)
    p = pressure_impl(mt,model,t)
    return convert_unit(u"Pa",unit,p)
end

function temperature(model::SaturationModel,st::ThermodynamicState,unit=u"K")
    return temperature(state_type(st),model,st,unit)
end
function temperature(mt::SingleSatP,model::SaturationModel,st::ThermodynamicState,unit)
    p = pressure(FromState(),st)
    t = temperature_impl(mt,model,p)
    return convert_unit(u"K",unit,t)
end

function critical_sat_interpolation(model::SaturationModel,p)
    #ln(pr) = h*(1-1/tr)
    #1/(1-ln(pr)/h)*tc = t
    tc = temperature(model,CriticalPoint())
    pc = pressure(model,CriticalPoint())
    t7 = 0.7*tc
    p7 = pressure_impl(QuickStates.sat_t(),model,t7)

    h = 2.3333333333333335*log(pc/p7)
    return 1/(1-log(p/pc)/h)*tc
end

#for multiple values


text:
module Get_kurtosis

export get_kurtosis

using Distributions, Statistics, Dierckx, SeisIO

"""
    get_kurtosis(data::SeisChannel,kurtsis_tw_sparse::Float64; timewinlength::Float64=60)

    compute kurtosis at each timewindow

# Input:
    - `data::SeisData`    : SeisData from SeisIO
    - `kurtosis_tw_sparse::Float64` : time length of span for kurtosis time window
    - `timewinlength::Float64`  : time window to calculate kurtosis
    kurtosis evaluation following Baillard et al.(2013)
"""
function get_kurtosis(data::SeisChannel, kurtosis_tw_sparse::Float64; timewinlength::Float64=60.0)

    #convert window lengths from seconds to samples
    TimeWin = trunc(Int,timewinlength * data.fs)
    data.misc["kurtosis"] = fast_kurtosis_series(data.x, TimeWin, kurtosis_tw_sparse)

    return data

end

"""
    fast_kurtosis_series(v::RealArray, TimeWin::Int64)

    fast compute kurtosis series at each timewindow

# Input:
    - `v::RealArray`    : SeisData from SeisIO
    - `N::Int64`  : time window length to calculate kurtosis

    kurtosis evaluation following Baillard et al.(2013)
"""
function fast_kurtosis_series(v::Array{Float64, 1}, TN::Int64, kurtosis_tw_sparse::Float64)

    kurt = zeros(length(v))
    n = length(v)
    kurt_grid = 1:n

    KTWSparse = kurtosis_tw_sparse

    if n < TN error("Kurtosis time window is larger than data length. Decrease time window.") end

    cm2 = 0.0  # empirical 2nd centered moment (variance)
    cm4 = 0.0  # empirical 4th centered moment

    # !DEPRECATED!: 1. compute mean value at each time window by numerical sequence
    # 2. use mapreduce to sum values

    # # first term
    # Trace = @views v[1:TN]
    # z2 = zeros(TN)
    # m0 = mean(Trace)
    #
    # cm2 = Statistics.varm(Trace, m0, corrected=false)
    # cm4 = fourthmoment(Trace, m0, corrected=false)
    #
    # # fill first part with kurtosis at TN
    # kurt[1:TN] .= (cm4 / (cm2 * cm2)) - 3.0
    #
    # @simd for k = TN:n-1
    #
    #     diff1 = @inbounds @views (v[k-TN+1] - v[k+1])/TN
    #     m1 = m0 - diff1
    #     Trace = @views v[k-TN+2:k+1]
    #     cm2 = Statistics.varm(Trace, m1, corrected=false)
    #     cm4 = fourthmoment(Trace, m1, corrected=false) #sum(xi - m)^4 / N
    #     kurt[k+1] = (cm4 / (cm2 * cm2)) - 3.0
    #     m0 = m1
    # end

    # first term
    Trace = @views v[1:TN]
    z2 = zeros(TN)
    m0 = mean(Trace)

    cm2 = Statistics.varm(Trace, m0, corrected=false)
    cm4 = fourthmoment(Trace, m0, corrected=false)

    # fill first part with kurtosis at TN
    kurt[1:TN] .= (cm4 / (cm2 * cm2)) - 3.0

    @simd for k = TN:n-1

        diff1 = @inbounds @views (v[k-TN+1] - v[k+1])/TN
        m1 = m0 - diff1
        Trace = @views v[k-TN+2:k+1]
        cm2 = Statistics.varm(Trace, m1, corrected=false)
        cm4 = fourthmoment(Trace, m1, corrected=false) #sum(xi - m)^4 / N
        kurt[k+1] = (cm4 / (cm2 * cm2)) - 3.0
        m0 = m1
    end

    return kurt

end

#---following functions are modified from Statistics.jl---#

centralizedabs4fun(m) = x -> abs2.(abs2.(x - m))
centralize_sumabs4(A::AbstractArray, m) =
    mapreduce(centralizedabs4fun(m), +, A)
centralize_sumabs4(A::AbstractArray, m, ifirst::Int, ilast::Int) =
    Base.mapreduce_impl(centralizedabs4fun(m), +, A, ifirst, ilast)

function centralize_sumabs4!(R::AbstractArray{S}, A::AbstractArray, means::AbstractArray) where S
    # following the implementation of _mapreducedim! at base/reducedim.jl
    lsiz = Base.check_reducedims(R,A)
    isempty(R) || fill!(R, zero(S))
    isempty(A) && return R

    if Base.has_fast_linear_indexing(A) && lsiz > 16 && !has_offset_axes(R, means)
        nslices = div(length(A), lsiz)
        ibase = first(LinearIndices(A))-1
        for i = 1:nslices
            @inbounds R[i] = centralize_sumabs4(A, means[i], ibase+1, ibase+lsiz)
            ibase += lsiz
        end
        return R
    end
    indsAt, indsRt = Base.safe_tail(axes(A)), Base.safe_tail(axes(R)) # handle d=1 manually
    keep, Idefault = Broadcast.shapeindexer(indsRt)
    if Base.reducedim1(R, A)
        i1 = first(Base.axes1(R))
        @inbounds for IA in CartesianIndices(indsAt)
            IR = Broadcast.newindex(IA, keep, Idefault)
            r = R[i1,IR]
            m = means[i1,IR]
            @simd for i in axes(A, 1)
                r += abs2(abs2(A[i,IA] - m))
            end
            R[i1,IR] = r
        end
    else
        @inbounds for IA in CartesianIndices(indsAt)
            IR = Broadcast.newindex(IA, keep, Idefault)
            @simd for i in axes(A, 1)
                R[i,IR] += abs2(abs2(A[i,IA] - means[i,IR]))
            end
        end
    end
    return R
end

function fourthmoment!(R::AbstractArray{S}, A::AbstractArray, m::AbstractArray; corrected::Bool=true) where S
    if isempty(A)
        fill!(R, convert(S, NaN))
    else
        rn = div(length(A), length(R)) - Int(corrected)
        centralize_sumabs4!(R, A, m)
        R .= R .* (1 // rn)
    end
    return R
end

"""
    fourthmoment(v, m; dims, corrected::Bool=true)
Compute the fourthmoment of a collection `v` with known mean(s) `m`,
optionally over the given dimensions. `m` may contain means for each dimension of
`v`. If `corrected` is `true`, then the sum is scaled with `n-1`,
whereas the sum is scaled with `n` if `corrected` is `false` where `n = length(v)`.
!!! note
    If array contains `NaN` or [`missing`](@ref) values, the result is also
    `NaN` or `missing` (`missing` takes precedence if array contains both).
    Use the [`skipmissing`](@ref) function to omit `missing` entries and compute the
    variance of non-missing values.
"""
fourthmoment(A::AbstractArray, m::AbstractArray; corrected::Bool=true, dims=:) = _fourthmoment(A, m, corrected, dims)

_fourthmoment(A::AbstractArray{T}, m, corrected::Bool, region) where {T} =
    fourthmoment!(Base.reducedim_init(t -> abs2(t)/2, +, A, region), A, m; corrected=corrected)

fourthmoment(A::AbstractArray, m; corrected::Bool=true) = _fourthmoment(A, m, corrected, :)

function _fourthmoment(A::AbstractArray{T}, m, corrected::Bool, ::Colon) where T
    n = length(A)
    n == 0 && return oftype((abs2(zero(T)) + abs2(zero(T)))/2, NaN)
    return centralize_sumabs4(A, m) / (n - Int(corrected))
end

end


text:
# Download some assets necessary for docs/testing not stored in the repo
import Downloads

const directory = joinpath(@__DIR__, "src", "examples")
mkpath(directory)

for (file, url) in [
        "periodic-rve.msh" => "https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/gh-pages/assets/periodic-rve.msh",
        "periodic-rve-coarse.msh" => "https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/gh-pages/assets/periodic-rve-coarse.msh",
    ]
    afile = joinpath(directory, file)
    if !isfile(afile)
        Downloads.download(url, afile)
    end
end


text:
module ErrorsModule

"""
```jldoctest
julia> a = 1
2

```

```jldoctest
```
"""
func(x) = x

end

using Documenter

makedocs(modules = [ErrorsModule])


text:
function goodbad(f, xs)
    good = []
    bad = []
    for x in xs
        try
            push!(good, x => f(x))
        catch e
            push!(bad, x => e)
        end
    end
    (good, bad)
end

"""
naive tester to separate the ones that lower and those that dont
"""
function test_suite()
    models = semantic()
    f = x -> ODESystem(readSBML(x))
    goodbad(f, models)
end

function lower_one(fn, df; verbose=false)
    k = 0
    n_dvs = 0
    n_ps = 0
    err = ""
    try
        ml = readSBML(fn)
        k = 1
        sys = ODESystem(ml)
        n_dvs = length(states(sys))
        n_ps = length(parameters(sys))
        k = 2
        prob  = ODEProblem(ml, (0, 1000.0))
        k = 3
        sol = solve(prob, TRBDF2(), dtmax=0.5; force_dtmin=true, unstable_check=unstable_check = (dt,u,p,t) -> any(isnan, u))
        k = 4
    catch e
        verbose && @info fn => e
        err = string(e)
        if length(err) > 1000 # cutoff since I got ArgumentError: row size (9088174) too large
            err = err[1:1000]
        end
    finally
        push!(df, (fn, k, n_dvs, n_ps, err))
        verbose && printstyled("$fn done with a code $k\n"; color=:green)
    end
end

function lower_fns(fns; write_fn=nothing)
    df = DataFrame(file=String[], retcode=Int[], n_dvs=Int[], n_ps=Int[], error=String[])
    # @sync Threads.@threads
    for fn in fns
        lower_one(fn, df)
    end
    write_fn !== nothing && CSV.write("logs/$(write_fn)", df)
    df
end


text:
using Logging
function ginicoeff(Dat::AbstractMatrix;dim=1,nosamplecorr=false)

    # Negative values or one-element series (not admitted)
    IDXnan = isnan.(Dat)
    IDX    = any(Dat.<0,dims=dim) .| (sum(.!IDXnan,dims=dim).<2)
    if dim == 1
        Dat[:,IDX[:]] .= 0
    else
        Dat[IDX[:],:] .= 0
    end

    if any(IDX)
        @warn "Check IDX for negative values or one-element series " * string(IDX)
    end

    # Total numel
    totNum = sum(.!IDXnan,dims=dim)

    # Replace NaNs
    Dat[IDXnan] .= 0

    # Sort In
    Dat = sort(Dat,dims=dim)

    # Calculate frequencies for each series
    freq = reverse(cumsum(ones(size(Dat)),dims=dim),dims=dim)

    # Totals
    tot = sum(Dat,dims=dim)

    # Gini's coefficient
    coeff = totNum .+ 1 .- 2 .*(sum(Dat.*freq,dims=dim)./tot)

    # Sample correction
    if nosamplecorr
        coeff = coeff./totNum
    else
        coeff = coeff./(totNum .-1 )
    end

    return coeff
end

ginicoeff(Dat::AbstractVector;nosamplecorr=false) = ginicoeff(reshape(Dat,length(Dat),1),nosamplecorr=false)[1]


text:
"""
    run_with_timeout(
        command; log=stdout, timeout = 5*60, name = "",
        wait_time = 1, verbose = true
    )

Runs `command` and pipes all output to `log`. The process will be terminated after
`timeout` seconds without any output. `name` describes the process in log messages,
and `verbose` determines whether meta-logs ("process started" etc.) will be printed.
"""
function run_with_timeout(
        command; log=stdout, timeout = 40*60, name = "",
        wait_time = 1, verbose = true, kill_timeout = 60,
        max_timeout = 3*60*60
    )
    print_interval = 60/wait_time
    print_in = print_interval

    out_io = IOBuffer()
    err_io = IOBuffer()
    out_file, err_file = "", ""
    if VERSION < v"1.1"
        out_file, out_io = mktemp()
        err_file, err_io = mktemp()
    end
    pipe = pipeline(command, stdout = out_io, stderr = err_io)
    process = run(pipe, wait = false)
    if VERSION < v"1.1"
        out_io = open(out_file)
        err_io = open(err_file)
    end
    job_start = timeout_start = time()
    task = @async begin
        logfallback = false
        io = try
            log isa String ? open(log, "w") : log
        catch err
            @error "Error opening logfile, falling back to stdout" error=err
            logfallback = true
            stdout
        end
        try
            tstart = time()
            verbose && @info("starting $name")
            while process_running(process)
                elapsed = time() - timeout_start
                total = time() - job_start
                if elapsed > timeout || total > max_timeout
                    verbose && @info("Terminating $name")
                    kill(process)
                    # Handle scenarios where SIGTERM is blocked/ignored/handled by the process
                    start_time = time()
                    while process_running(process)
                        if time() - start_time > kill_timeout
                            verbose && @info("Killing $name")
                            kill(process, Base.SIGKILL)
                        end
                        sleep(5)
                    end
                    break
                end
                errstr, outstr = readstr_buffer.((out_io, err_io))
                is_silent = length(errstr) == 0 && length(outstr) == 0
                isempty(outstr) || println(io, outstr)
                isempty(errstr) || println(io, errstr)

                # if something printed reset timeout
                if !is_silent
                    timeout_start = time()
                end

                sleep(wait_time)
                print_in -= 1

                if print_in <= 0 && verbose
                    print(".")
                    print_in = print_interval
                end
            end

            verbose && println()
            verbose && @info("$name completed in $(round(time() - tstart, digits=1)) seconds")
        catch err
            @error "Error while running $(name) with timeout." error=err
        finally
            errstr, outstr = readstr_buffer.((out_io, err_io))
            isempty(outstr) || println(io, outstr)
            isempty(errstr) || println(io, errstr)

            flush(io)
            if log isa String && !logfallback
                close(io)
            end
        end
    end
    return process, task
end

function readstr_buffer(x::IOStream)
    return read(x, String)
end

function readstr_buffer(x::Base.GenericIOBuffer{Array{UInt8,1}})
    return String(take!(x))
end


text:
include((@__DIR__)*"/../preamble.jl")

include("reflecting_model/model_def.jl")
include("default_params.jl")

# look at time to exit of the interval [0,1]
# use the bounded fluid queue, but change the boundary to be at 1
model = BoundedFluidQueue(
    model.T[1:2,1:2],rates(model)[1:2],[0.0 1.0],[1.0 0.0],1.0,
)

dg_model(order) = build_discretised_model(DGMesh{t},model,(1.0/3)/((order+1)/2),2)
# order1_model(order) = build_discretised_model(DGMesh{t},model,1.0/3/order,1)
# qbdrap_model(order) = build_discretised_model(FRAPMesh{t},model,1.0/3,order)

d0_map_point_mass(dq) = (interior_point_mass(eps(),1,dq).coeffs)
d0_map_exp(dq) = (SFMDistribution((x,i)->(i∈(1:2))*exp(-x)/(1-exp(-model.b))/2,dq).coeffs)

# change d0_map here to do the make_approximations for ohter initial confitions
# remember to change write directory below too
models = make_approximations(orders,approx_types,d0_map_exp,(args...)->args[1])

for k1 in keys(models)
    for k2 in keys(models[k1])
        models[k1][k2]["B"][1,:] = spzeros(1,size(models[k1][k2]["B"],2))
        models[k1][k2]["B"][end,:] = spzeros(1,size(models[k1][k2]["B"],2))
        droptol!(models[k1][k2]["B"],eps()*1000)
    end
end

pth = mkpath((@__DIR__)*"/hitting_times/data/exp")

h = 0.005/3
t_vec = 0.0:h:10.0
hitting_times_cdf = zeros(length(t_vec),3)
hitting_times_cdf[:,1] = t_vec
for k1 in keys(models)
    for k2 in keys(models[k1])
        for (c,t) in enumerate(t_vec)
            models[k1][k2]["coeffs_mapped"] =
                DiscretisedFluidQueues._integrate(
                    models[k1][k2]["coeffs_mapped"],
                    models[k1][k2]["B"],
                    h,
                    StableRK4(h),
                    x->GeneralisedMUSCL.fun(x,GeneralisedMUSCL.generate_params(models[k1][k2]["dq"])...),
                )
            hitting_times_cdf[c,2:3] = models[k1][k2]["coeffs_mapped"][[1;end]]
        end
        CSV.write(
            pth*"/order_"*string(k1)*"_model_"*string(k2)*".csv",
            DataFrame(
                t=hitting_times_cdf[:,1],
                phase_2=hitting_times_cdf[:,2],
                phase_1=hitting_times_cdf[:,3],
            )
        )
    end
end


text:
# Docu for profiling Cuda Code in Julia: https://cuda.juliagpu.org/stable/development/profiling/

#using BenchmarkTools
using CUDA

a = rand(1024,1024,100);

ad = CUDA.rand(1024,1024,100);

@time sin.(a)

sin.(ad);

println("Finished")


text:
using Test
using Aqua
using Tables
using Transducers

using Distributed
nprocs() > 1 || addprocs(1)

@everywhere using TrackingTimers

@time @testset "Basics + `show`" begin
    t = TrackingTimer()
    @test sprint(show, t) == "TrackingTimer(…)"
    @test occursin("No entries.", sprint(show, MIME"text/plain"(), t))

    TrackingTimers.@timeit t "timeit" sleep(1)
    @test occursin("timeit", sprint(show, MIME"text/plain"(), t))
    TrackingTimers.synchronize!(t) # in case it evaluates the RHS first
    @test Tables.rows(t) == t.results
    @test length(t.results) == 1
    @test t.results[1].pid == 1
end

@time @testset "Distributed" begin
    t = TrackingTimer()

    @everywhere slp(s) = (sleep(s / 100); s)
    slp_t = t(slp)

    slp_t(1)
    @test length(Tables.rows(t)) == 1
    slp_t(2)
    @test length(Tables.rows(t)) == 2

    @test pmap(slp_t, 1:10) == 1:10
    @test length(t.results) == 2 # haven't synchronized yet
    @test length(Tables.rows(t)) == 12 # have synchronized

    pids = Tables.getcolumn(Tables.columns(t), :pid)
    @test pids isa AbstractVector{Int}
    @test 2 ∈ pids
end

@time @testset "Threaded" begin
    t = TrackingTimer()

    sqrt_t = t(sqrt)
    result = @sync [Threads.@spawn sqrt_t(i) for i in 1:10]
    @test fetch.(result) == sqrt.(1:10)
    TrackingTimers.synchronize!(t)
    @test length(t.results) == 10

    tids = Tables.getcolumn(Tables.columns(t), :thread_id)
    @test tids isa AbstractVector{Int}
    @test 2 ∈ tids
end

@time @testset "Transducers" begin
    t = TrackingTimer()
    xs = 1:1000
    sin_t = t(sin)
    t_result = foldxt(+, Map(sin_t), xs)
    tids = Tables.getcolumn(Tables.columns(t), :thread_id)

    if Threads.nthreads() > 1
        @test length(unique(tids)) >= 2
    end
    @test length(Tables.rows(t)) == length(xs)

    d_result = foldxd(+, Map(sin_t), xs)
    pids = Tables.getcolumn(Tables.columns(t), :pid)
    @test procs() ⊆ pids

    TrackingTimers.synchronize!(t)
    @test length(Tables.rows(t)) == 2 * length(xs) == length(t.results)
end

@time @testset "Tables interface" begin
    t = TrackingTimer()
    xs = 1:1000
    sin_t = t(sin)
    foldxt(+, Map(sin_t), xs)

    @test Tables.istable(t)
    @test Tables.rowaccess(typeof(t))
    @test !Tables.columnaccess(typeof(t))
    TrackingTimers.synchronize!(t) # in case it evaluates the RHS first
    @test Tables.rows(t) == t.results

    col_tbl = Tables.columns(t)
    @test Tables.columnnames(col_tbl) ==
          (:name, :time, :gctime, :n_allocs, :bytes, :thread_id, :pid)
    @test col_tbl.name isa AbstractVector{String}
    @test col_tbl.time isa AbstractVector{Float64}
    @test Tables.getcolumn(col_tbl, :time) == col_tbl.time
end

@time @testset "Aqua tests" begin
    Aqua.test_all(TrackingTimers; ambiguities=false)
end

max_stars_repo_name: ericphanson/TrackingTimers

max_issues_repo_name: ericphanson/TrackingTimers

max_forks_repo_name: ericphanson/TrackingTimers


text:
struct WebSocketClient <: WebSocketHandler
    endpoint:: String
    client:: WSClient
    callbacks:: Dict{UInt64, Function}

    # https://coincheck.com/ja/documents/exchange/api#websocket-overview
    WebSocketClient() = new("wss://ws-api.coincheck.com", WSClient(), Dict())
end
export WebSocketClient

on_text(handler:: WebSocketClient, s:: String) = begin
    try
        for x in handler.callbacks
            x[2](JSON.parse(s))
        end
    catch ex
        # TODO error-handling
        println(ex)
    end
end
function connect(client:: WebSocketClient)
    # Connect to Coincheck server
    wsconnect(client.client, URI(client.endpoint), client)
end
function consume!(callback:: Function, client:: WebSocketClient)
    while true
        id = rand(0:UInt64(1) << 63)
        if !haskey(client.callbacks, id)
            client.callbacks[id] = callback
            return id
        end
    end
end
function delete!(id:: UInt64, client:: WebSocketClient)
    if haskey(client.callbacks, id)
        Base.delete!(client.callbacks, id)
    end
end
function subscribe(client:: WebSocketClient, channels)
    # Make requests
    for channel = collect(channels)
        json = JSON.json(Dict("type" => "subscribe", "channel" => channel))
        send_text(client.client, json)
    end
end


text:
module jInvSeismic
include("Utils/Utils.jl")
include("BasicFWI/BasicFWI.jl")
include("EikonalInv/EikonalInv.jl")
include("FWI/FWI.jl")
end


text:
module MultistartOptimization

export MinimizationProblem, LocationValue, NLoptLocalMethod, local_minimization, TikTak,
    multistart_minimization

using ArgCheck: @argcheck
using Base.Threads: @spawn, fetch
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF
using NLopt: NLopt
using Parameters: @unpack
using Sobol: SobolSeq, Sobol

####
#### structures for the problem and results
####

"""
$(TYPEDEF)

Wrapper for a minimization problem.

# Fields

$(FIELDS)
"""
struct MinimizationProblem{F,T}
    "The function to be minimized."
    objective::F
    "Lower bounds (a vector of real numbers)."
    lower_bounds::T
    "Upper bounds (a vector of real numbers)."
    upper_bounds::T
    # FIXME constructor checks
end

"""
$(TYPEDEF)

A location-value pair.

# Fields

$(FIELDS)
"""
struct LocationValue{T <: AbstractVector{<:Real}, S <: Real}
    "Location (a vector of real numbers)."
    location::T
    "The value of the objective at `location`."
    value::S
end

####
#### internal utilities
####

"""
$(SIGNATURES)

Evaluate and return points of an `N`-element Sobol sequence.

An effort is made to parallelize the code using `Threads` when available.
"""
function sobol_starting_points(minimization_problem::MinimizationProblem, N::Integer)
    @unpack objective, lower_bounds, upper_bounds = minimization_problem
    s = SobolSeq(lower_bounds, upper_bounds)
    skip(s, N)                  # better uniformity
    map(fetch, map(x -> @spawn(LocationValue(x, objective(x))), Iterators.take(s, N)))
end

"""
$(SIGNATURES)

Helper function to keep the `N` points with the lowest `value`.
"""
function _keep_lowest(xs, N)
    @argcheck 1 ≤ N ≤ length(xs)
    partialsort(xs, 1:N, by = p -> p.value)
end

####
#### local minimization
####

Base.@kwdef struct NLoptLocalMethod
    algorithm::NLopt.Algorithm
    xtol_abs::Float64 = 1e-8
    xtol_rel::Float64 = 1e-8
    maxeval::Int = 0
    maxtime::Float64 = 0.0
end

"""
Solve `minimization_problem` using `local_method`, starting from `x`. Return a
`LocationValue`.
"""
function local_minimization(opt, x)
    optf, optx, ret = NLopt.optimize(opt, x)
    # FIXME ret is ignored, save it? is it useful?
    LocationValue(optx, optf)
end

####
#### multistart minimization
####

struct TikTak
    quasirandom_N::Int
    initial_N::Int
    θ_min::Float64
    θ_max::Float64
    θ_pow::Float64
end

"""
$(SIGNATURES)

The “TikTak” multistart method, as described in *Arnoud, Guvenen, and Kleineberg (2019)*.

This implements the *multistart* part, can be called with arbitrary local methods, see
[`multistart_minimization`](@ref).

# Arguments

- `quasirandom_N`: the number of quasirandom points for the first pass (using a Sobol
  sequence).

# Keyword arguments

- `keep_ratio`: the fraction of best quasirandom points which are kept

- `θ_min` and `θ_max` clamp the weight parameter, `θ_pow` determines the power it is raised
  to.

The defaults are from the paper cited above.
"""
function TikTak(quasirandom_N; keep_ratio = 0.1, θ_min = 0.1, θ_max = 0.995, θ_pow = 0.5)
    @argcheck 0 < keep_ratio ≤ 1
    TikTak(quasirandom_N, ceil(keep_ratio * quasirandom_N), θ_min, θ_max, θ_pow)
end

function _weight_parameter(t::TikTak, i)
    @unpack initial_N, θ_min, θ_max, θ_pow = t
    clamp((i / initial_N)^θ_pow, θ_min, θ_max)
end

"""
$(SIGNATURES)

Solve `minimization_problem` by using `local_method` within `multistart_method`.
"""
function multistart_minimization(multistart_method::TikTak, local_method,
                                 minimization_problem)
    @unpack quasirandom_N, initial_N, θ_min, θ_max, θ_pow = multistart_method
    quasirandom_points = sobol_starting_points(minimization_problem, quasirandom_N)
    initial_points = _keep_lowest(quasirandom_points, initial_N)
    function _step(visited_minimum, (i, initial_point))
        θ = _weight_parameter(multistart_method, i)
        x = @. (1 - θ) * initial_point.location + θ * visited_minimum.location
        local_minimum = local_minimization(local_method, minimization_problem, x)
        local_minimum.value < visited_minimum.value ? local_minimum : visited_minimum
    end
    foldl(_step, enumerate(initial_points); init = first(initial_points))
end

end # module