text:
#=
See - http://ape-package.ird.fr/misc/BitLevelCodingScheme.html for the bit-level
coding scheme :

Nucleotide	    IUPAC code	 Bit-level code	  Value
---------------------------------------------------
A	                  A	      10001000	       136
G	                  G	      01001000	        72
C	                  C	      00101000	        40
T	                  T	      00011000	        24
A or G	              R	      11000000	       192
A or C	              M	      10100000	       160
A or T	              W	      10010000	       144
G or C	              S	      01100000	        96
G or T	              K	      01010000	        80
C or T	              Y	      00110000	        48
A or G or C	          V	      11100000	       224
A or C or T	          H	      10110000	       176
A or G or T	          D	      11010000	       208
G or C or T	          B	      01110000	       112
A or G or C or T      N	      11110000	       240
Alignment gap        (–)	  00000100	         4
Unknown character    (?)	  00000010	         2
---------------------------------------------------

This has been adapated (last two rows have changed) to
allow matches between alignment gaps / ? and any nucleotide:

Nucleotide	    IUPAC code	 Bit-level code	  Value
---------------------------------------------------
A	                  A	      10001000	       136
G	                  G	      01001000	        72
C	                  C	      00101000	        40
T	                  T	      00011000	        24
A or G	              R	      11000000	       192
A or C	              M	      10100000	       160
A or T	              W	      10010000	       144
G or C	              S	      01100000	        96
G or T	              K	      01010000	        80
C or T	              Y	      00110000	        48
A or G or C	          V	      11100000	       224
A or C or T	          H	      10110000	       176
A or G or T	          D	      11010000	       208
G or C or T	          B	      01110000	       112
A or G or C or T      N	      11110000	       240
Alignment gap        (–)	  11110100	       244  ###### CHANGED
Unknown character    (?)	  11110010	       242  ###### CHANGED
---------------------------------------------------

Function	               C code	                      Value returned
------------------------------------------------------------------------------------------
KnownBase(a)	           a & 8	                      8 if a is known surely
IsAdenine(a)	           a == 136	                      1 if a is adenine
IsGuanine(a)	           a == 72	                      1 if a is guanine
IsCytosine(a)       	   a == 40	                      1 if a is cytosine
IsThymine(a)	           a == 24	                      1 if a is thymine
IsPurine(a)	               a & 55	                      0 if a is a purine
IsPyrimidine(a)	           a & 199	                      0 if a is a pyrimidine
DifferentBase(a, b) 	   (a & b) < 16	                  1 if a and b are different surely
SameBase(a, b)	           KnownBase(a) && a == b	      1 if a and b are the same surely
-------------------------------------------------------------------------------------------
=#

# byte_dict is a map from the Char representation of IUPAC nucleotide
# codes to uint8 values according to the (modified) bit-level coding scheme
# for nucleotides developed by Emmanuel Paradis:
# http://ape-package.ird.fr/misc/BitLevelCodingScheme.html
function make_byte_dict()
    byte_dict = Dict{Char, UInt8}('A' => 136,
                                    'G' => 72,
                                    'C' => 40,
                                    'T' => 24,
                                    'R' => 192,
                                    'M' => 160,
                                    'W' => 144,
                                    'S' => 96,
                                    'K' => 80,
                                    'Y' => 48,
                                    'V' => 224,
                                    'H' => 176,
                                    'D' => 208,
                                    'B' => 112,
                                    'N' => 240,
                                    '-' => 244,
                                    '?' => 242)

    return byte_dict
end

#=
byte_array is a map from the byte/uint8 representation of IUPAC characters
as they would appear in ASCII text (e.g. a fasta file) to the E. Paradis
bit-level coding scheme. E.g. 'A', whose byte value is 65, is element 65
in the array below, which maps to '0x88'/136/10001000 - which is the encoding
for A, and so on for the other nucleotides. The highest value is Y, which is 89
(and encodes 0x30 - 48/00110000). This makes for faster look up than converting
the file from bytes to strings/chars then looking up the corresponding byte
in a dict.
=#
function make_byte_array()
    byte_array = [0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00, #10
                  0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00, #20
                  0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00, #30
                  0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00, #40
                  0x00,0x00,0x00,0x00,0xf4,0x00,0x00,0x00,0x00,0x00, #50
                  0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00, #60
                  0x00,0x00,0xf2,0x00,0x88,0x70,0x28,0xd0,0x00,0x00, #70
                  0x48,0xb0,0x00,0x00,0x50,0x00,0xa0,0xf0,0x00,0x00, #80
                  0x00,0xc0,0x60,0x18,0x00,0xe0,0x90,0x00,0x30,0x00] #90
    return byte_array
end

# nuc_dict maps from uint8 representation of nucleotides to their IUPAC code
function make_nuc_dict()
    nuc_dict = Dict{UInt8, Char}(136 => 'A',
                                 72 => 'G',
                                 40 => 'C',
                                 24 => 'T',
                                 192 => 'R',
                                 160 => 'M',
                                 144 => 'W',
                                 96 => 'S',
                                 80 => 'K',
                                 48 => 'Y',
                                 224 => 'V',
                                 176 => 'H',
                                 208 => 'D',
                                 112 => 'B',
                                 240 => 'N',
                                 244 => '-',
                                 242 => '?')

    return nuc_dict
end

function make_score_array()
    A = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #20
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #40
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #60
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #80
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #100
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #120
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #140
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #160
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #180
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #200
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #220
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #240
         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] #260

    A[136] = 12
    A[72] = 12
    A[40] = 12
    A[24] = 12
    A[192] = 6
    A[160] = 6
    A[144] = 6
    A[96] = 6
    A[80] = 6
    A[48] = 6
    A[224] = 4
    A[176] = 4
    A[208] = 4
    A[112] = 4
    A[240] = 3
    A[244] = 3
    A[242] = 3

    return A
end

# score_dict is a map from uint8 values for IUPAC nucleotide codes to an
# integer for how unambiguous they are. The score is caculated as:
# 12 * 1/possible real nucleotides. E.g. an 'A' / 136::uint8 scores 12, but an
# 'N' (A, C, G or T) / 240::uint8 scores 3
function make_score_dict()
    score_dict = Dict{UInt8, Int64}(136 => 12,
                                     72 => 12,
                                     40 => 12,
                                     24 => 12,
                                     192 => 6,
                                     160 => 6,
                                     144 => 6,
                                     96 => 6,
                                     80 => 6,
                                     48 => 6,
                                     224 => 4,
                                     176 => 4,
                                     208 => 4,
                                     112 => 4,
                                     240 => 3,
                                     244 => 3,
                                     242 => 3)

    return score_dict
end

max_stars_repo_name: cov-ert/julialign

max_issues_repo_name: cov-ert/julialign

max_forks_repo_name: cov-ert/julialign


text:
module FoldsDagger

export DaggerEx, foldx_dagger, transduce_dagger

using Dagger: DArray, delayed, distribute
using SplittablesBase: amount, halve
using Referenceables: ReferenceableArray, referenceable
using Transducers:
    Executor,
    Map,
    Reduced,
    Transducer,
    Transducers,
    combine,
    complete,
    foldl_nocomplete,
    reduced,
    start,
    unreduced

# TODO: Don't import internals from Transducers:
using Transducers:
    DefaultInit,
    DefaultInitOf,
    EmptyResultError,
    IdentityTransducer,
    maybe_usesimd,
    restack,
    retransform

include("core.jl")
include("darray.jl")

end

ext: jl


text:
# module TestBasics

using Test, Dates, Temporal

@testset "Initialization" begin
    @test Temporal.autocol(1:27)[end] == :AA
    @test Temporal.autoidx(2) == [today()-Day(1), today()]
    @test isa(X, TS)
    @test isa(Y, TS)
    @test size(X) == size(Y) == (N, K)
    @test X.values == Y.values == data
    @test X.index == collect(dates)
    @test Y.index == collect(times)
    @test size(TS(rand(N,K))) == (N,K)
    @test TS(X.values, X.index, ['A','B','C','D']) == X
    @test TS(X.values, X.index, ["A","B","C","D"]) == X
    @test TS(X.values[:,1], X.index, 'A') == X[:,1]
    @test TS(X.values[:,1], X.index, "A") == X[:,1]
    @test TS(X.values[1,:]', X.index[1], X.fields) == X[1]
    @test size(TS()) == (0,0)
    # underscores in colnames (issue 13)
    u = TS(rand(N), X.index, "under_score")
    @test u.fields[1] == :under_score
    show(stdout, TS(rand(252,4)))
    show(stdout, TS(rand(252,4))[1,:])
    print("\n")
    # iterator protocol
    for (t, x) in X
        @test t in X.index
        @test size(x,1) == size(X,2)
    end
end

@testset "Rename" begin
    @testset "inplace" begin
        @testset "Pair" begin
            ts = TS(data, times, [:a, :b, :c, :d])
            rename!(ts, :a => :A)
            @test ts.fields == [:A, :b, :c, :d]
        end

        @testset "Several Pairs" begin
            ts = TS(data, times, [:a, :b, :c, :d])
            rename!(ts, :a => :A, :d => :D)
            @test ts.fields == [:A, :b, :c, :D]
        end

        @testset "Dict" begin
            ts = TS(data, times, [:a, :b, :c, :d])
            d = Dict(:a => :A, :d => :D)
            rename!(ts, d...)
            @test ts.fields == [:A, :b, :c, :D]
        end

        @testset "Lambda Function" begin
            ts = TS(data, times, [:a, :b, :c, :d])
            f = colname -> Symbol(uppercase(string(colname)))
            rename!(f, ts)
            @test ts.fields == [:A, :B, :C, :D]
        end

        @testset "Function with composition" begin
            ts = TS(data, times, [:a, :b, :c, :d])
            f = Symbol ∘ uppercase ∘ string
            rename!(f, ts)
            @test ts.fields == [:A, :B, :C, :D]
        end

        @testset "Function with do block" begin
            ts = TS(data, times, [:a, :b, :c, :d])
            rename!(ts) do x
                x |> string |> uppercase |> Symbol
            end
            @test ts.fields == [:A, :B, :C, :D]
        end

        @testset "Function with automatic String/Symbol" begin
            ts = TS(data, times, [:a, :b, :c, :d])
            rename!(uppercase, ts, String)
            @test ts.fields == [:A, :B, :C, :D]
        end
    end

    @testset "not in place" begin
        @testset "Pair" begin
            ts1 = TS(data, times, [:a, :b, :c, :d])
            ts2 = rename(ts1, :a => :A)
            @test ts1.fields == [:a, :b, :c, :d]
            @test ts2.fields == [:A, :b, :c, :d]
        end

        @testset "Several Pairs" begin
            ts1 = TS(data, times, [:a, :b, :c, :d])
            ts2 = rename(ts1, :a => :A, :d => :D)
            @test ts1.fields == [:a, :b, :c, :d]
            @test ts2.fields == [:A, :b, :c, :D]
        end

        @testset "Lambda Function" begin
            ts1 = TS(data, times, [:a, :b, :c, :d])
            f = colname -> Symbol(uppercase(string(colname)))
            ts2 = rename(f, ts1)
            @test ts1.fields == [:a, :b, :c, :d]
            @test ts2.fields == [:A, :B, :C, :D]
        end

        @testset "Function with automatic String/Symbol" begin
            ts1 = TS(data, times, [:a, :b, :c, :d])
            ts2 = rename(uppercase, ts1, String)
            @test ts1.fields == [:a, :b, :c, :d]
            @test ts2.fields == [:A, :B, :C, :D]
        end
    end
end

# end


text:
module Chart

using DataFrames #this needs to be here, Thyme or Scallop not good enough
using Compose
using Gadfly

export xchart,
       priceplot,
       returnsplot

include("xchart.jl")
include("priceplot.jl")
include("returnsplot.jl")

end


text:
# Could be called ActualNewton...
struct Newton <: NearlyNewtonMethod end


text:
using LinearAlgebra
using SparseArrays

using SymbolicWedderburn
using DynamicPolynomials

include(joinpath(dirname(@__DIR__), "examples", "action_polynomials.jl"))
include(joinpath(dirname(@__DIR__), "examples", "dihedral.jl"))
include(joinpath(dirname(@__DIR__), "examples", "util.jl"))

using JuMP
using SCS

const OPTIMIZER = optimizer_with_attributes(
    SCS.Optimizer,
    "acceleration_lookback" => 10,
    "max_iters" => 3_000,
    "alpha" => 1.2,
    "eps" => 1e-6,
    "linear_solver" => SCS.DirectSolver,
)

@polyvar x y
const robinson_form = x^6 + y^6 - x^4 * y^2 - y^4 * x^2 - x^4 - y^4 - x^2 - y^2 + 3x^2 * y^2 + 1

struct DihedralAction <: OnMonomials end

SymbolicWedderburn.coeff_type(::DihedralAction) = Float64
function SymbolicWedderburn.action(::DihedralAction, el::DihedralElement, mono::AbstractMonomial)
    if iseven(el.reflection + el.id)
        var_x, var_y = x, y
    else
        var_x, var_y = y, x
    end
    sign_x = 1 <= el.id <= 2 ? -1 : 1
    sign_y = 2 <= el.id ? -1 : 1
    return mono([x, y] => [sign_x * var_x, sign_y * var_y])
end

@testset "Dihedral Action" begin
    G = DihedralGroup(4)
    @test all(
        SymbolicWedderburn.action(DihedralAction(), g, robinson_form) == robinson_form
        for g in DihedralGroup(4)
    )

    T = Float64

    wedderburn = WedderburnDecomposition(
        T,
        G,
        DihedralAction(),
        monomials([x,y], 0:DynamicPolynomials.maxdegree(robinson_form)),
        monomials([x,y], 0:(DynamicPolynomials.maxdegree(robinson_form) ÷ 2))
    )

    M = let bh = monomials([x,y], 0:(DynamicPolynomials.maxdegree(robinson_form) ÷ 2))
        [SymbolicWedderburn.basis(wedderburn)[x * y] for x in bh, y in bh]
    end

    m = let m = JuMP.Model(OPTIMIZER)
        JuMP.@variable m t
        JuMP.@objective m Max t
        psds = [
            JuMP.@variable(m, [1:d, 1:d] in PSDCone())
            for d in size.(direct_summands(wedderburn), 1)
        ]

        # preallocating
        M_orb = similar(M, T)
        # these could be preallocated as well:
        # Mπs = zeros.(T, size.(psds))
        # tmps = SymbolicWedderburn._temps(wedderburn)

        C = DynamicPolynomials.coefficients(robinson_form - t, basis(wedderburn))

        for iv in invariant_vectors(wedderburn)
            c = dot(C, iv)
            M_orb = invariant_constraint!(M_orb, M, iv)
            # Mπs = SymbolicWedderburn.diagonalize!(Mπs, M_orb, wedderburn, tmps)
            Mπs = SymbolicWedderburn.diagonalize(M_orb, wedderburn)

            JuMP.@constraint m sum(
                dot(Mπ, Pπ) for (Mπ, Pπ) in zip(Mπs, psds) if !iszero(Mπ)
            ) == c
        end
        m
    end

    optimize!(m)

    @test isapprox(value(m[:t]), -3825 / 4096, rtol = 1e-4)
    status = termination_status(m)
    @test status ∈ (MOI.OPTIMAL, MOI.ALMOST_OPTIMAL)
end


text:
abstract type AbstractDeviceFormulation end

mutable struct DeviceModel{D <: PSY.Device,
                           B <: PSI.AbstractDeviceFormulation}
    device::Type{D}
    formulation::Type{B}
end

mutable struct InitialCondition{T <: Union{PJ.ParameterRef, Float64}}
    device::PSY.Device
    value::T
end


text:
#mainly just exercies the basic constructors

n = 8

serialComm = SerialComm{Int32, Bool, Int16}()
srcMap = BlockMap(n, n, serialComm)
desMap = BlockMap(n, n, serialComm)

function basicTest(impor)
    if isa(impor, Import)
        data = impor.importData
    else
        # basic exports are about the same anyways
        data = impor.exportData
    end
    @test srcMap == data.source
    @test desMap == data.target
    @test n == data.numSameIDs
    @test isa(data.distributor, Distributor{Int32, Bool, Int16})
    @test [] == data.permuteToLIDs
    @test [] == data.permuteFromLIDs
    @test [] == data.remoteLIDs
    @test [] == data.exportLIDs
    @test [] == data.exportPIDs
    @test true == data.isLocallyComplete
end

# basic import
basicTest(Import(srcMap, desMap))
basicTest(Import(srcMap, desMap, nothing))

# basic export
basicTest(Export(srcMap, desMap))
basicTest(Export(srcMap, desMap, nothing))

max_stars_repo_path: test/Import-Export Tests.jl

max_issues_repo_path: test/Import-Export Tests.jl

max_forks_repo_path: test/Import-Export Tests.jl


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

name = "boost"
version = v"1.71.0"

# Collection of sources required to build boost
sources = [
    ArchiveSource("https://dl.bintray.com/boostorg/release/$(version)/source/boost_$(version.major)_$(version.minor)_$(version.patch).tar.bz2",
                  "d73a8da01e8bf8c7eda40b4c84915071a8c8a0df4a6734537ddde4a8580524ee"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/boost*/

CXX=$CXX_FOR_BUILD ./bootstrap.sh --prefix=$prefix --without-libraries=python --with-toolset=gcc

# Patch adapted from
# https://svnweb.freebsd.org/ports/head/devel/boost-libs/files/patch-boost_math_tools_config.hpp?revision=439932&view=markup
# to be able to build long double math libraries
atomic_patch -p1 ../patches/boost_math_tools_config_hpp.patch

rm project-config.jam
toolset=gcc
targetos=linux
extraargs=

# BinaryBuilderBase compiler wrappers don't like it when we use -march=ANYTHING.
# So we patch that out. However, we cannot just insert an empty string; so we instead
# add another "harmless" option; we choose `-Wall` which is already passed anyway
sed -i "s/-march=i686/-Wall/g" tools/build/src/tools/gcc.*

if [[ $target == *apple* ]]; then
    targetos=darwin
    toolset=darwin-6.0
    extraargs="binary-format=mach-o link=shared"
    echo "using darwin : 6.0 : $CXX : <cxxflags>-stdlib=libc++ <linkflags>-stdlib=libc++ ;" > project-config.jam
elif [[ $target == x86_64*mingw* ]]; then
    targetos=windows
    extraargs="address-model=64 binary-format=pe abi=ms link=shared"
elif [[ $target == i686*mingw* ]]; then
    targetos=windows
    extraargs="address-model=32 binary-format=pe abi=ms link=shared"
elif [[ $target == *freebsd* ]]; then
    targetos=freebsd
    toolset=clang-6.0
    extraargs="address-model=64 link=shared"
    echo "using clang : 6.0 : $CXX : <linkflags>\\"$LDFLAGS\\" ;" > project-config.jam
fi
./b2 -j${nproc} toolset=$toolset target-os=$targetos $extraargs variant=release --prefix=$prefix --without-python --layout=system install

install_license LICENSE_1_0.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libboost_atomic", :libboost_atomic),
    LibraryProduct("libboost_chrono", :libboost_chrono),
    LibraryProduct("libboost_container", :libboost_container),
    LibraryProduct("libboost_context", :libboost_context),
    LibraryProduct("libboost_contract", :libboost_contract),
    LibraryProduct("libboost_coroutine", :libboost_coroutine),
    LibraryProduct("libboost_date_time", :libboost_date_time),
    LibraryProduct("libboost_filesystem", :libboost_filesystem),
    LibraryProduct("libboost_graph", :libboost_graph),
    LibraryProduct("libboost_iostreams", :libboost_iostreams),
    # boost_locale segfaults on windows, see https://github.com/benlorenz/boostBuilder/issues/2
    #LibraryProduct("libboost_locale", :libboost_locale),
    LibraryProduct("libboost_log", :libboost_log),
    LibraryProduct("libboost_log_setup", :libboost_log_setup),
    LibraryProduct("libboost_math_c99", :libboost_math_c99),
    LibraryProduct("libboost_math_c99l", :libboost_math_c99l),
    LibraryProduct("libboost_math_c99f", :libboost_math_c99f),
    LibraryProduct("libboost_math_tr1", :libboost_math_tr1),
    LibraryProduct("libboost_math_tr1f", :libboost_math_tr1f),
    LibraryProduct("libboost_math_tr1l", :libboost_math_tr1l),
    LibraryProduct("libboost_prg_exec_monitor", :libboost_prg_exec_monitor),
    LibraryProduct("libboost_program_options", :libboost_program_options),
    LibraryProduct("libboost_random", :libboost_random),
    LibraryProduct("libboost_regex", :libboost_regex),
    LibraryProduct("libboost_serialization", :libboost_serialization),
    LibraryProduct("libboost_stacktrace_basic", :libboost_stacktrace_basic),
    LibraryProduct("libboost_stacktrace_noop", :libboost_stacktrace_noop),
    LibraryProduct("libboost_system", :libboost_system),
    LibraryProduct("libboost_thread", :libboost_thread),
    LibraryProduct("libboost_timer", :libboost_timer),
    LibraryProduct("libboost_type_erasure", :libboost_type_erasure),
    LibraryProduct("libboost_unit_test_framework", :libboost_unit_test_framework),
    LibraryProduct("libboost_wave", :libboost_wave),
    LibraryProduct("libboost_wserialization", :libboost_wserialization),
]

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

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using Compat
import Compat.view

# Simple static kd-trees.

abstract KDNode

immutable KDTree{T <: AbstractFloat}
	# A matrix of n, m-dimensional observations
	xs::AbstractMatrix{T}

	# Permutaiton of the data used by the structure (to avoid sorting or
	# otherwise modifying xs)
	perm::Vector{Int}

	root::KDNode

	verts::Set{Vector{T}}

	# Top-level bounding box
	bounds::Matrix{T}
end

# Construct a kd-tree
#
# Args:
#   xs: Data organized in the tree.
#   leaf_size_factor: Stop spliting if a node contains
#       fewer than leaf_size_factor * n elements.
#   leaf_diameter_factor: Stop spliting if a node's bounding
#       hypercube has diameter less that leaf_diameter_factor
#       times the diameter of the root node's bounding hypercube.
#
# Returns:
#   A KDTree object
#
function KDTree{T <: AbstractFloat}(xs::AbstractMatrix{T},
	                                leaf_size_factor=0.05,
	                                leaf_diameter_factor=0.0)
	n, m = size(xs)
	perm = collect(1:n)

	bounds = Array(T, (2, m))
	for j in 1:m
		col = xs[:,j]
		bounds[1, j] = minimum(col)
		bounds[2, j] = maximum(col)
	end
	diam = diameter(bounds)

	leaf_size_cutoff = ceil(Integer, (leaf_size_factor * n))
	leaf_diameter_cutoff = leaf_diameter_factor * diam
	verts = Set{Vector{T}}()

	# Add verticies defined by the bounds
	for vert in product([bounds[:,j] for j in 1:m]...)
		push!(verts, T[vert...])
	end

	root = build_kdtree(xs, view(perm, 1:n), bounds,
		                leaf_size_cutoff, leaf_diameter_cutoff, verts)

	KDTree(xs, collect(1:n), root, verts, bounds)
end

immutable KDLeafNode <: KDNode
end

immutable KDInternalNode{T <: AbstractFloat} <: KDNode
	j::Int # dimension on which the data is split
	med::T # median value where the split occours
	leftnode::KDNode
	rightnode::KDNode
end

# Compute the diamater of a hypercube defined by bounds
#
# Args:
#   bounds: A 2 by n matrix where bounds[1,i] gives the lower bound in
#           the ith dimension and bonuds[2,j] the upper.
#
# Returns:
#   Computed diameter
#
function diameter(bounds::Matrix)
	euclidean(vec(bounds[1,:]), vec(bounds[2,:]))
end

# Recursively build a kd-tree
#
# Args:
#   xs: Data being orginized.
#   perm: Permutation of the data, used to avoid
#       directly sorting or modifying xs.
#   bounds: Bounding hypercube of the node.
#   leaf_size_cutoff: stop spliting on nodes with more than
#       this many values.
#   leaf_diameter_cutoff: stop splitting on nodes with less
#        than this diameter.
#   verts: current set of vertexes
#
# Modifies:
#   perm, verts
#
# Returns:
#   Either a KDLeafNode or a KDInternalNode
#
function build_kdtree{T}(xs::AbstractMatrix{T},
	                     perm::SubArray,
	                     bounds::Matrix{T},
	                     leaf_size_cutoff::Int,
	                     leaf_diameter_cutoff::T,
	                     verts::Set{Vector{T}})
	n, m = size(xs)

	if length(perm) <= leaf_size_cutoff || diameter(bounds) <= leaf_diameter_cutoff
		return KDLeafNode()
	end

	# split on the dimension with the largest spread
	j = 1
	maxspread = 0
	for k in 1:m
		xmin = Inf
		xmax = -Inf
		for i in perm
			xmin = min(xmin, xs[i, k])
			xmax = max(xmax, xs[i, k])
		end
		if xmax - xmin > maxspread
			maxspread = xmax - xmin
			j = k
		end
	end

	# find the median and partition
	if isodd(length(perm))
		mid = div(length(perm), 2)
		select!(perm, mid, by=i -> xs[i, j])
		med = xs[perm[mid], j]
		mid1 = mid
		mid2 = mid + 1
	else
		mid1 = div(length(perm), 2)
		mid2 = mid1 + 1
		select!(perm, mid1, by=i -> xs[i, j])
		select!(perm, mid2, by=i -> xs[i, j])
		med = (xs[perm[mid1], j] + xs[perm[mid2], j]) / 2
	end

	leftbounds = copy(bounds)
	leftbounds[2, j] = med
	leftnode = build_kdtree(xs, view(perm, 1:mid1), bounds,
		                    leaf_size_cutoff, leaf_diameter_cutoff, verts)

	rightbounds = copy(bounds)
	rightbounds[1, j] = med
	rightnode = build_kdtree(xs, view(perm, mid2:length(perm)), bounds,
		                     leaf_size_cutoff, leaf_diameter_cutoff, verts)

	coords = Array(Array, m)
	for i in 1:m
		if i == j
			coords[i] = [med]
		else
			coords[i] = bounds[:, i]
		end
	end

	for vert in product(coords...)
		push!(verts, T[vert...])
	end

	KDInternalNode{T}(j, med, leftnode, rightnode)
end

# Given a bounding hypecube, return its verticies
function bounds_verts(bounds::Matrix)
	collect(product([bounds[:, i] for i in 1:size(bounds, 2)]...))
end

# Traverse the tree to the bottom and return the verticies of
# the leaf node's bounding hypercube.
function traverse{T}(kdtree::KDTree{T}, xs::AbstractVector{T})
	m = size(kdtree.bounds, 2)

	if length(xs) != m
		error("$(m)-dimensional kd-tree searched with a length $(length(x)) vector.")
	end

	for j in 1:m
		if xs[j] < kdtree.bounds[1, j] || xs[j] > kdtree.bounds[2, j]
			error(
				"""
				Loess cannot perform extrapolation. Predict can only be applied
				to points within the bounding hypercube of the data used to train
				the model.
				""")
		end
	end

	bounds = copy(kdtree.bounds)
	node = kdtree.root
	while !isa(node, KDLeafNode)
		if xs[node.j] <= node.med
			bounds[2, node.j] = node.med
			node = node.leftnode
		else
			bounds[1, node.j] = node.med
			node = node.rightnode
		end
	end

	bounds_verts(bounds)
end


text:
abstract type AbstractParticle{T <: AbstractFloat} end

abstract type AbstractParticleActivator end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
nbits(U) = round(Int, log2(size(U,1)))

@generated function bits(n, ::Type{Val{N}}) where N
  :($([:(Bool((n >> $(i-1)) & 0x01)) for i = N:-1:1]...),)
end

@generated function int(bits::NTuple{N}) where N
  reduce((x, b) -> :(($x<<1)+$b), [:(bits[$i]) for i = 1:N]; init = 0)
end

int(n::Bool) = int((n,))

function classical(f, ::Type{Val{N}}) where N
  U = zeros(Int, 2^N, 2^N)
  for i = 0:2^N-1
    U[int(f(bits(i, Val{N})...))+1, i+1] = 1
  end
  return U
end

classical(f, N) = classical(f, Val{N})

permutation(is::NTuple{N}) where N =
  classical((x...) -> map(i -> x[i], is), Val{N})

pad(U, l, r) = foldr(⊗, [repeated(I, l)..., U, repeated(I, r)...])

headfirst(n, is...) = permutation((is..., setdiff(1:n, is)...))

function measure(n, m)
  M = zeros(Int, 2^n, 2^(n+m))
  for i = 0:2^n-1, j = 0:2^m-1
    mn = i<<m+j
    M[i+1, mn+1] = 1
  end
  return M
end


text:
#=
downsample.jl
timing tests
=#

import MIRT
using BenchmarkTools: @btime

function downsample3_time()
	x = ones(Float32, 100, 200, 40)
	@btime MIRT.downsample3_perm($x, (2, 2, 2)) # fastest
	@btime MIRT.downsample3_loop($x, [2, 2, 2], T=Float32) # 2x slower
	nothing
end

#= ir28 with 1.4.2
  3.3 ms (63 allocations: 8.01 MiB)
  5.4 ms (200011 allocations: 9.54 MiB)
=#

downsample3_time()


text:
function my_handler(myendpoints::Dict{Int64,Array{NamedTuple{(:route, :toexec, :urlparams),Tuple{Union{String,Regex},Function,Union{Nothing,Dict{Int64,String}}}},1}},tonumber::Dict{String,Char},formats::Dict{String,Function},cleanurl::Function,constantheaders::Array{Pair{String,String},1})

	function myqueryparams(input::SubString{String})::Dict{String,String}
		salida = Dict{String,String}()
        for e in split(input,"&")
        	items = split(e,"=")
        	if (length(items)==2)
        		salida[items[1]]=items[2]
            end
        end
        return salida
	end

	function urlparams(input::String,params::Dict{Int64,String})::Dict{String,String}
		salida = Dict{String,String}()

        for (index, value) in enumerate(split(input,"/"))
        	if haskey(params, index)
        		salida[params[index]]=value
            end
        end #845.574 ns (16 allocations: 768 bytes)

        return salida
	end

	function processrequest(request::HTTP.Request)

		data = split(request.target,"?")
		nvalues= 0
	    body = ""
	    myquery= Dict()

	    mycleanurl= cleanurl(data[1])
	    if(length(data[1])>1)
		  nvalues = length(split(mycleanurl,"/"))
		end

	    if(!eof(IOBuffer(HTTP.payload(request))))
	    	body = get(formats, HTTP.header(request,"Content-Type"), formats["*/*"] )(String(request.body))
        end

        if (length(data)>1)
            myquery = myqueryparams(data[2])
        end

	    return (
			query = myquery
			,body = body
			,headers = request.headers
			,searchroute = parse(Int64, string(tonumber[request.method] , nvalues) )
			,url = mycleanurl
	    )

    end

	function search(value)
	    function run(item)
			return occursin(item.route,value)
		end
		()-> (run)
	end

	function service(myrequest)
		myresponse = HTTP

		myparams = Dict{String,String}()
		result = get(myendpoints, myrequest.searchroute, 1)

		if (typeof(result)!== Int64)
			f = iterate(Iterators.filter(search(myrequest.url).run, result))

			if f !== nothing

				if f[1].urlparams !== nothing
					myparams = urlparams(myrequest.url, f[1].urlparams)
				end

				return f[1].toexec(
			   	(
			   		 query = myrequest.query
			   		,body = myrequest.body
			   		,headers = myrequest.headers
			   		,params= myparams
			   	)
				,myresponse
			    )

		   end
		end

		return myendpoints[0][1].toexec(myrequest,myresponse)
	end

	function handler(request::HTTP.Request)
              myresponse = service(processrequest(request))
              map( x -> HTTP.setheader(myresponse,x ), constantheaders)
              return myresponse
	end

    ()->(handler)
end

	function useCORS(;AllowOrigins::String = "*", AllowHeaders::String = "Origin, Content-Type, Accept", AllowMethods::String = "GET,POST,PUT,DELETE", MaxAge::String = "178000")
		headersalways(["Access-Control-Allow-Origin" => AllowOrigins])
		CORSenabled[1] = (request, HTTP) -> begin

        return HTTP.Response(200
        	, HTTP.mkheaders(["Content-Type" => "text/plain"
			,"Access-Control-Allow-Methods" => AllowMethods
        	,"Access-Control-Allow-Headers" => AllowHeaders
        	,"Access-Control-Max-Age" => MaxAge])
        	, body="")
        end

	end

	function headersalways(values::Array{Pair{String,String},1})
	  map( x -> push!(constantheaders, x), values)
	  @info(values)
	end

	function notfound(text::String)
	  if occursin(".html", text)
	    notfound_message= File(text)
	    addnotfound(notfound_message,myendpoints)
	  else
	    addnotfound(text,myendpoints)
	  end
	end

	function start( ;host::String = "127.0.0.1", port::Int64 = 8086, verbose::Bool = false, sslconfig::Union{MbedTLS.SSLConfig, Nothing}=nothing)
		Handler = my_handler(myendpoints,tonumber,formats,cleanurl,constantheaders)
		if ':' in host
			HTTP.serve(Handler.handler, Sockets.IPv6(host), port, verbose=verbose, sslconfig=sslconfig)
		end

		HTTP.serve(Handler.handler, Sockets.IPv4(host), port, verbose=verbose, sslconfig=sslconfig)
	end


text:
function Base.show(io::IO, bls::BLSPeriodogram{T,S}) where {T,S}
    N = length(bls.power)
    print(io, "$N-element BLSPeriodogram{$T,$S}")
end

function Base.show(io::IO, ::MIME"text/plain", bls::BLSPeriodogram)
    p = params(bls)
    println(io, "BLSPeriodogram\n==============")
    println(io, "input dim: ", length(bls.t))
    println(io, "output dim: ", length(bls.power))
    println(io, "period range: ", range_str(extrema(bls.periods)))
    println(io, "duration range: ", range_str(extrema(bls.duration_in)))
    println(io, "objective: ", bls.objective)
    println(io, "\nparameters\n----------")
    println(io, "index: ", p.index)
    println(io, "period: ", p.period)
    println(io, "duration: ", p.duration)
    println(io, "t0: ", p.t0)
    println(io, "depth: ", p.depth, " ± ", p.depth_err)
    println(io, "snr: ", p.snr)
    print(io, "log-likelihood: ", p.loglike)
end

range_str((min, max)) = "$min - $max"


text:
# # Polarizability by linear response
#
# We compute the polarizability of a Helium atom. The polarizability
# is defined as the change in dipole moment
# ```math
# \mu = \int r ρ(r) dr
# ```
# with respect to a small uniform electric field ``E = -x``.
#
# We compute this in two ways: first by finite differences (applying a
# finite electric field), then by linear response. Note that DFTK is
# not really adapted to isolated atoms because it uses periodic
# boundary conditions. Nevertheless we can simply embed the Helium
# atom in a large enough box (although this is computationally wasteful).
#
# As in other tests, this is not fully converged, convergence
# parameters were simply selected for fast execution on CI,

using DFTK
using LinearAlgebra

a = 10.
lattice = a * I(3)  # cube of ``a`` bohrs
He = ElementPsp(:He, psp=load_psp("hgh/lda/He-q2"))
atoms = [He => [[1/2; 1/2; 1/2]]]  # Helium at the center of the box

kgrid = [1, 1, 1]  # no k-point sampling for an isolated system
Ecut = 30
tol = 1e-8

## dipole moment of a given density (assuming the current geometry)
function dipole(basis, ρ)
    rr = [a * (r[1] - 1/2) for r in r_vectors(basis)]
    sum(rr .* ρ) * basis.dvol
end;

# ## Using finite differences
# We first compute the polarizability by finite differences.
# First compute the dipole moment at rest:
model = model_LDA(lattice, atoms; symmetries=false)
basis = PlaneWaveBasis(model; Ecut, kgrid)
res   = self_consistent_field(basis, tol=tol)
μref  = dipole(basis, res.ρ)

# Then in a small uniform field:
ε = .01
model_ε = model_LDA(lattice, atoms; extra_terms=[ExternalFromReal(r -> -ε * (r[1] - a/2))],
                    symmetries=false)
basis_ε = PlaneWaveBasis(model_ε; Ecut, kgrid)
res_ε = self_consistent_field(basis_ε, tol=tol)
με = dipole(basis_ε, res_ε.ρ)

#-
polarizability = (με - μref) / ε

println("Reference dipole:  $μref")
println("Displaced dipole:  $με")
println("Polarizability :   $polarizability")

# The result on more converged grids is very close to published results.
# For example [DOI 10.1039/C8CP03569E](https://doi.org/10.1039/C8CP03569E)
# quotes **1.65** with LSDA and **1.38** with CCSD(T).

# ## Using linear response
# Now we use linear response to compute this analytically; we refer to standard
# textbooks for the formalism. In the following, ``\chi_0`` is the
# independent-particle polarizability, and ``K`` the
# Hartree-exchange-correlation kernel. We denote with ``\delta V_{\rm ext}`` an external
# perturbing potential (like in this case the uniform electric field). Then:
# ```math
# \delta\rho = \chi_0 \delta V = \chi_0 (\delta V_{\rm ext} + K \delta\rho),
# ```
# which implies
# ```math
# \delta\rho = (1-\chi_0 K)^{-1} \chi_0 \delta V_{\rm ext}.
# ```
# From this we identify the polarizability operator to be ``\chi = (1-\chi_0 K)^{-1} \chi_0``.
# Numerically, we apply ``\chi`` to ``\delta V = -x`` by solving a linear equation
# (the Dyson equation) iteratively.

using KrylovKit

## Apply (1- χ0 K)
function dielectric_operator(δρ)
    δV = apply_kernel(basis, δρ; ρ=res.ρ)
    χ0δV = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, δV)
    δρ - χ0δV
end

## δVext is the potential from a uniform field interacting with the dielectric dipole
## of the density.
δVext = [-a * (r[1] - 1/2) for r in r_vectors(basis)]
δVext = cat(δVext; dims=4)

## Apply χ0 once to get non-interacting dipole
δρ_nointeract = apply_χ0(res.ham, res.ψ, res.εF, res.eigenvalues, δVext)

## Solve Dyson equation to get interacting dipole
δρ = linsolve(dielectric_operator, δρ_nointeract, verbosity=3)[1]

println("Non-interacting polarizability: $(dipole(basis, δρ_nointeract))")
println("Interacting polarizability:     $(dipole(basis, δρ))")

# As expected, the interacting polarizability matches the finite difference
# result. The non-interacting polarizability is higher.


text:
import Term.renderables: Renderable, RenderableText, AbstractRenderable

@testset "\e[34mRenderables - Renderable" begin
    # simply test the creation and type of Renderables
    @test Renderable("text") isa RenderableText

    r1 = Renderable("asdasda")
    @test Renderable(r1) isa AbstractRenderable
    @test Renderable(r1.segments[1]) isa AbstractRenderable
end

@testset "\e[34mRenderables - RenderableText" begin
    @test RenderableText("sadasdasdasdasdas") isa RenderableText

    r = RenderableText("a"^100; width = 25)
    @test r.measure.w == 25
    @test r.measure.h == 4

    r = RenderableText("a"^100, "red bold"; width = 25)
    @test r.measure.w == 25
    @test r.measure.h == 4
end


text:
using ODMXMLTools
using Documenter

makedocs(
        modules = [ODMXMLTools],
        sitename = "ODMXMLTools.jl",
        authors = "Vladimir Arnautov",
        pages = [
            "Home" => "index.md",
            ],
        )

deploydocs(repo = "github.com/PharmCat/ODMXMLTools.jl.git", devbranch = "main", forcepush = true
)


text:
# =========================================================================== #
# Activation Functions
# =========================================================================== #

abstract type AbstractActivationFunction
end

struct TanhActivationFunction <: AbstractActivationFunction
end

function (act_fcn::TanhActivationFunction)(x)
    return tanh(x)
end

function differentiate(act_fcn::TanhActivationFunction, x)
    return one(typeof(x)) - act_fcn(x)^2
end

function eval(act_fcn::TanhActivationFunction, x)
    return act_fcn(x)
end

struct ReluActivationFunction <: AbstractActivationFunction
end

function (act_fcn::ReluActivationFunction)(x)
    T = typeof(x)
    z = zero(T)
    return (x < z) * z + (x >= z) * x
end

function differentiate(act_fcn::ReluActivationFunction, x)
    T = typeof(x)
    z = zero(T)
    return (x < z) * z + (x >= z) * one(T)
end

function eval(act_fcn::ReluActivationFunction, x)
    return act_fcn(x)
end

struct IdentityActivationFunction <: AbstractActivationFunction
end

function (act_fcn::IdentityActivationFunction)(x)
    return x
end

function differentiate(act_fcn::IdentityActivationFunction, x)
    return one(typeof(x))
end

function eval(act_fcn::IdentityActivationFunction, x)
    return act_fcn(x)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
export Arduino, sendCommand

abstract type Arduino <: Device end

@mustimplement cmdStart(ard::Arduino)
@mustimplement cmdEnd(ard::Arduino)
@mustimplement cmdDelim(ard::Arduino)
@mustimplement serialDevice(ard::Arduino)

function sendCommand(ard::Arduino, cmdString::String)
  cmd = cmdStart(ard) * cmdString * cmdEnd(ard) * cmdDelim(ard)
  return query(serialDevice(ard), cmd)
end

Base.@kwdef struct SimpleArduino <: Arduino
  commandStart::String = "!"
  commandEnd::String = "*"
  delim::String = "#"
  sd::SerialDevice
end
cmdStart(ard::SimpleArduino) = ard.commandStart
cmdEnd(ard::SimpleArduino) = ard.commandEnd
cmdDelim(ard::SimpleArduino) = ard.delim
serialDevice(ard::SimpleArduino) = ard.sd


text:
module ImarisWriter

using Libdl
using ImarisWriter_jll
#using Lz4_jll, Zlib_jll, HDF5_jll

# Need the following DLLs
# 1. bpImarisWriter96.dll
# 2. lz4.dll
# 3. zlib.dll
# 4. hdf5.dll

# Write your package code here.

include("imariswriter_common.jl")
include("imariswriter_api.jl")

#const bpImarisWriter96_path = joinpath( @__DIR__(), "..", "bpImarisWriter96.dll")
#const bpImarisWriter96 = :bpImarisWriter96
const bpImarisWriter96_path = ImarisWriter_jll.libbpImarisWriter96_path
const bpImarisWriter96 = ImarisWriter_jll.libbpImarisWriter96
const bpImarisWriter96_h = ImarisWriter_jll.libbpImarisWriter96_handle

function __init__()
    #global bpImarisWriter96_h
    #bpImarisWriter96_h =  dlopen( bpImarisWriter96_path )
    @info bpImarisWriter96_path
end

end


text:
# coding: utf-8

# In[1]:

#heap sort is the alternative version of selection sort
#heap sort is implemented on a binary heap instead of a list
#well, u can argue binary heap is expressed in list
#binary heap is a complete binary tree
#for each node it has no more than two children
#the tree is as far left as possible
#for this binary heap, we keep minimum at root
#similar to selection sort
#we only need to get the smallest number for each round
#and our list keeps shrinking until we are left with one element for final round
#however, the traversal is kinda different from selection sort
#what we do is to compare two children with parent node
#we keep the smallest among the family of three as parent node for each node
#until we reach the root and we successfully keep the smallest at root
#next step is to remove the root and we find the minimum for the new tree
#we would do it recursively in this script
#otherwise we would have to use loops in loops

#in this version we use min heap
#if u intend to use max heap,check out geeks for geeks
# https://www.geeksforgeeks.org/heap-sort/

# In[2]:

function heap_sort(arr)

    #denote layer as the height of the tree
    layer=0
    maxlen=length(arr)

    #when we are left with one element
    #that is the base case for recursion
    if maxlen<2

        return arr

    end

    #this part is to find out how many layers we got for binary heap
    #note that layer starts from 1 in julia
    while maxlen>2^layer

        maxlen-=2^layer
        layer+=1

    end

    #as layer starts from 1
    #we have to use layer+1 to make sure each layer has been traveled
    #until we reach the base case, layer 1, the root
    for i in layer+1:-1:1

        #a special property of binary heap
        #the beginning of each layer equals to 2^(i-1)
        #each layer contains 2^(i-1) elements
        for j in 2^(i-1):2^i-1

            #a special property of binary heap as well
            #the index of left child is always even number
            #and for the right child, always odd number
            right=j*2+1
            left=j*2

            #we use try function in case the node is the leaf
            #and if there is no left child
            #there wont be right child
            #according to the far left principle of binary heap
            #we compare the parent with two children nodes
            #we only keep the smallest as parent node
            try

                if arr[left]<arr[j]

                    arr[left],arr[j]=arr[j],arr[left]

                end

                try

                    if arr[right]<arr[j]

                        arr[right],arr[j]=arr[j],arr[right]

                    end

                    if arr[right]<arr[left]

                        arr[right],arr[left]=arr[left],arr[right]

                    end

                catch e

                    if !(isa(e,BoundsError))

                        throw(e)

                    end

                end

            catch e

                if !(isa(e,BoundsError))

                    throw(e)

                end

            end

        end

    end

    #recursively, we shrink the size of the list
    #by removing its root which is the smallest element for the current tree
    arr[2:end]=heap_sort(arr[2:end])

    return arr

end

# In[3]:

#heap sort s time complexity is o(nlogn)
#however, this script is written in recursive function
#it may not achieve the same speed as iterations

# In[4]:

for _ in 1:100

    test_arr=rand(100)

    if !(heap_sort(test_arr)==sort(test_arr))
        printstyled("Erreur",color=:red)
    end

end

max_stars_repo_path: heap sort.jl

max_stars_repo_name: je-suis-tm/search-and-sort

max_issues_repo_path: heap sort.jl

max_issues_repo_name: HafizSabia/search-and-sort

max_forks_repo_path: heap sort.jl

max_forks_repo_name: HafizSabia/search-and-sort


text:
# SimLynx/src/variable.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

include("statistics.jl")
include("history.jl")

"""
A variable automatically maintains the history and statistics of its value over
time.
"""
mutable struct Variable{T<:Real}
    value::T
    prev_update::Float64
    history::Union{History{T}, Nothing}
    stats::Union{Stats{T}, Nothing}
    function Variable{T}(value::T;
                         data=:accumulate,
                         history=false,
                         stats=true) where {T<:Real}
        if data == :none
            return new(value, 0.0, nothing, nothing)
        elseif data == :tally
            return new(value, 0.0,
                       history ? TalliedHistory{T}(value) : nothing,
                       stats ? TalliedStats{T}(value) : nothing)
        elseif data == :accumulate
            return new(value, 0.0,
                       history ? AccumulatedHistory{T}(value) : nothing,
                       stats ? AccumulatedStats{T}(value) : nothing)
        else
            throw(ArgumentError("expected :none, :tally, or :accumulate, got $data"))
        end
    end
    function Variable{T}(;
                         data=:accumulate,
                         history=false,
                         stats=true) where {T<:Real}
        if data == :none
            return new(zero(T), 0.0, nothing, nothing)
        elseif data == :tally
            return new(zero(T), 0.0,
                       history ? TalliedHistory{T}() : nothing,
                       stats ? TalliedStats{T}() : nothing)
        elseif data == :accumulate
            return new(zero(T), 0.0,
                       history ? AccumulatedHistory{T}() : nothing,
                       stats ? AccumulatedStats{T}() : nothing)
        else
            throw(ArgumentError("expected :none, :tally, or :accumulate, got $data"))
        end
    end
end

function Base.getproperty(var::Variable{T}, name::Symbol) where T<:Real
    value = getfield(var, :value)
    # If they just want the value, return it.
    if name == :value
        return value
    end
    prev_update = getfield(var, :prev_update)
    # If they just want the previous update time, return it.
    if name == :prev_update
        return prev_update
    end
    # Otherwise, synchronize the variable and return the requested field.
    history = getfield(var, :history)
    stats = getfield(var, :stats)
    duration = current_time() - prev_update
    if duration > 0.0
        if isa(history, AccumulatedHistory)
            update!(history, value, duration)
        end
        if isa(stats, AccumulatedStats)
            update!(stats, value, duration)
        end
        setfield!(var, :prev_update, current_time())
    end
    if name == :history
        return history
    elseif name == :stats
        return stats
    else
        throw(ArgumentError("expected a Variable field, given $name"))
    end
end

"Return the current value of a variable. (deprecated)"
# get(var::Variable) = var.value

function Base.setproperty!(var::Variable{T}, name::Symbol, value::T) where T<:Real
    if name!= :value
        throw(ArgumentError("expected Variable field :value, given $name"))
    end
    # Synchronize the variable
    old_value = getfield(var, :value)
    prev_update = getfield(var, :prev_update)
    history = getfield(var, :history)
    stats = getfield(var, :stats)
    duration = current_time() - prev_update
    if duration > 0.0
        if isa(history, AccumulatedHistory)
            update!(history, old_value, duration)
        end
        if isa(stats, AccumulatedStats)
            update!(stats, old_value, duration)
        end
        setfield!(var, :prev_update, current_time())
    end
    # Update the variable's value
    setfield!(var, :value, value)
    if isa(history, AccumulatedHistory)
        update!(history, var.value, 0.0)
    end
    if isa(history, TalliedHistory)
        update!(history, var.value)
    end
    if isa(stats, AccumulatedStats)
        update!(stats, var.value, 0.0)
    end
    if isa(stats, TalliedStats)
        update!(stats, var.value)
    end
    return nothing
end

"Set the value of a variable. (deprecated)"
# set!(var::Variable, value) = var.value = value

increment!(var::Variable, x) = var.value += x
decrement!(var::Variable, x) = var.value -= x

function print_stats(var::Variable;
                     title::String="Statistics")
    stats = var.stats
    if !isnothing(stats)
        print_stats(stats, title=title)
    end
    return nothing
end

function plot_history(var::Variable;
                      file::Union{String, Nothing}=nothing,
                      title="History")
    history = var.history
    if !isnothing(history)
        plot_history(history, file=file, title=title)
    end
    return nothing
end


text:
"""
    cost,i1,i2 = fastdtw(seq1,seq2,dist,radius)

Perform dynamic-time warping using the FastDTW algorithm to measure the distance between two sequences. Note that regular DTW often performs better than FastDTW https://arxiv.org/abs/2003.11246

Computes FastDTW approximation to the DTW, described in Salvador & Chan,
Intelligent Data Analysis (2007).

See also [`dtw`](@ref), [`dtw_cost`](@ref), [`dtwnn`](@ref).
"""
function fastdtw(
        seq1::AbstractArray,
        seq2::AbstractArray,
        dist::Union{SemiMetric, Function},
        radius::Int,
    )

    MinSize = max(radius + 2, 10)
    N1 = lastlength(seq1)
    N2 = lastlength(seq2)
    if N1 <= MinSize || N2 <= MinSize
        return (dtw(seq1, seq2, dist))
    end

    # Call recursively on a pair of sequences half this length
    compressed1 = compress2(seq1)
    compressed2 = compress2(seq2)
    _cost, lowrescol, lowresrow = fastdtw(compressed1, compressed2, dist, radius)

    # Now resample that path to the finer resolution, find the correct
    # window around it, and get the DTW given that window.
    hirescol, hiresrow = expandpath(lowrescol, lowresrow, N1, N2)
    idx2min, idx2max = computewindow(hirescol, hiresrow, radius)
    cost1, newcol, newrow = dtw(seq1, seq2, dist, idx2min, idx2max)
end

# Given a path through low-res space, generate an approximate path
# through high-res space. It should have dimension Ncol x Nrow

@inbounds function expandpath(lowrescol, lowresrow, Ncol, Nrow)
    @assert div(Ncol+1,2) == lowrescol[end]
    @assert div(Nrow+1,2) == lowresrow[end]
    Np = length(lowrescol)
    @assert Np == length(lowresrow)

    hirescol = zeros(eltype(lowrescol), 2*Np)
    hiresrow = zeros(eltype(lowresrow), 2*Np)
    hirescol[1] = hiresrow[1] = c = r = 1
    for i=1:Np-1
        # Select plan according to the next move in lowres path.
        if lowrescol[i+1] == lowrescol[i]  # Next move is up
            r += 1
            hirescol[2*i] = c
            hiresrow[2*i] = r
            r += 1
            hirescol[2*i+1] = c
            hiresrow[2*i+1] = r

        elseif lowresrow[i+1] == lowresrow[i] # Next move is sideways
            c += 1
            hirescol[2*i] = c
            hiresrow[2*i] = r
            c += 1
            hirescol[2*i+1] = c
            hiresrow[2*i+1] = r

        else  # Next move is diagonal.
            c += 1; r += 1
            hirescol[2*i] = c
            hiresrow[2*i] = r
            c += 1; r += 1
            hirescol[2*i+1] = c
            hiresrow[2*i+1] = r
        end
    end
    hirescol[end] = Ncol
    hiresrow[end] = Nrow
    # When expanding to an odd numbered size, it's possible to repeat
    # the last step.  Fix that:
    if hirescol[end]==hirescol[end-1] && hiresrow[end]==hiresrow[end-1]
        hirescol = hirescol[1:end-1]
        hiresrow = hiresrow[1:end-1]
    end
    hirescol, hiresrow
end

# yshort = compress2(y)
#   Returns a shortened time series that is half the length of the input sequence.
#   The length of the compressed sequence is always even.
function compress2(seq::AbstractArray)
    # Navg = div(length(seq), 2)
    n = lastlength(seq)
    evenseq = 0.5*(seq[!,1:2:n-1]+seq[!,2:2:n])
    if lastlength(seq)%2 == 1
        return cat(evenseq, seq[!, n], dims=ndims(seq))
    end
    evenseq
end

# Given the lists of (col,row) indices for the optimal path, compute a "window"
# around that path of the given radius.
# Returns (rowmin, rowmax), each a vector of length pathcols[end], representing
# for each column, the minimum and maximum row numbers used in that column.

@inbounds function computewindow(pathcols, pathrows, radius)
    Np = length(pathcols)
    @assert Np == length(pathrows)
    Ncol = pathcols[end]
    Nrow = pathrows[end]

    # Find the min/max row at each column in the path.
    pathmin = zeros(Int, Ncol)
    pathmax = zeros(Int, Ncol)
    for i=1:Np
        c,r = pathcols[i], pathrows[i]
        pathmax[c] = r
        if pathmin[c] == 0
            pathmin[c] = r
        end
    end

    # The window in each column for "radius" r starts at the pathmin
    # of the rth-previous column and ends at the pathmax of the
    # rth-next column, plus (in each case) the radius.
    if radius < Ncol-1 && radius < Nrow-1
        rowmin = vcat(fill(1,radius), pathmin[1:end-radius] .- radius)
        rowmax = vcat(pathmax[radius+1:end] .+ radius, fill(Nrow,radius))

        # Window values must be in the range [1:Nrow].
        for c=1:Ncol
            if rowmin[c]<1; rowmin[c]=1; end
            if rowmax[c]>Nrow; rowmax[c]=Nrow; end
        end
    else
        rowmin = fill(1,Ncol)
        rowmax = fill(Nrow,Ncol)
    end
    rowmin, rowmax
end


text:
# copypasta from kemal for fasta and fastq i/o

# i/o functions use file parsing from BioSequences package

### FASTA files

"""
    read_fasta_records(filename)

Read .fasta file contents.
"""
function read_fasta_records(filename)
    stream = open(FASTA.Reader, filename)
    records = FASTA.Record[]
    for entry in stream
        push!(records, entry)
    end
    return records
end

"""
    read_fasta(filename; seqtype=String)

Read .fasta file contents, parse, and return sequences as type `seqtype`
"""
function read_fasta(filename; seqtype=String)
    records = read_fasta_records(filename)
    return seqtype[FASTA.sequence(seqtype, r) for r in records]
end

"""
    read_fasta_with_names(filename; seqtype=String)

Read .fasta file contents, parse, and return names and sequences as type `seqtype`.
"""
function read_fasta_with_names(filename; seqtype=String)
    records = read_fasta_records(filename)
    return [FASTA.identifier(r) for r in records], seqtype[FASTA.sequence(seqtype, r) for r in records]
end

"""
    read_fasta_with_names_in_other_order(filename; seqtype=String)

Read .fasta file contents, parse, and return sequences as type `seqtype` with names.
Now that's what I call convenience.
"""
function read_fasta_with_names_in_other_order(filename::String; seqtype=String)
    names, seqs = read_fasta_with_names(filename, seqtype=seqtype)
    return seqs, names
end

"""
    write_fasta(filename::String, seqs; names = String[])

Write given `seqs` and optional `names` to a .fasta file with given filepath.
"""
function write_fasta(filename::String, seqs; names = String[])
    if length(names) > 0 && length(names) != length(seqs)
        error("number of sequences does not match number of names")
    end
    if length(names) == 0
        names = ["seq_$i" for i in 1:length(seqs)]
    end
    stream = open(FASTA.Writer, filename)
    for (name, seq) in zip(names, seqs)
        write(stream, FASTA.Record(name, seq))
    end
    close(stream)
end

### FASTQ files

"""
    read_fastq_records(filename)

Read .fastq file contents.
"""
function read_fastq_records(filename)
    stream = open(FASTQ.Reader, filename)
    records = FASTQ.Record[]
    for record in stream
        if any([q < 0 for q in record.quality])
            error("$(record.name) in $filename contains negative phred values")
        end
        push!(records, record)
    end
    return records
end

"""
    read_fastq(filename; seqtype=String)

Read .fastq file contents, parse, and return sequences as `seqtype` type, phreds, and names.
"""
function read_fastq(filename; seqtype=String, min_length=nothing, max_length=nothing, err_rate=nothing)
    records = read_fastq_records(filename)
    seqs = seqtype[]
    phreds = Vector{Phred}[]
    names = String[]
    for record in records
		if err_rate != nothing && mean(phred_to_p(Array{Int8,1}(FASTQ.quality(record, :sanger)))) > err_rate
			continue
		end
		if min_length != nothing && length(FASTQ.sequence(seqtype, record)) < min_length
			continue
		end
		if max_length != nothing && length(FASTQ.sequence(seqtype, record)) > max_length
			continue
        end
		push!(seqs, FASTQ.sequence(seqtype, record))
        push!(phreds, FASTQ.quality(record, :sanger))
        push!(names, FASTQ.identifier(record))
    end
    return seqs, phreds, names
end

"""
    write_fastq(filename, seqs, phreds::Vector{Vector{Phred}};
                     names=String[], DNASeqType = false)

Write given sequences, phreds, names to .fastq file with given file path.
If `names` not provided, gives names 'seq_1', etc.
"""
function write_fastq(filename, seqs, phreds::Vector{Vector{Phred}};
                     names=String[], DNASeqType = false)
    if !DNASeqType
        seqs = [DNASequence(s) for s in seqs]
    end
    stream = open(FASTQ.Writer, filename)
    i = 0
    if length(names) != length(seqs)
        names = [string("seq_", i) for i in 1:length(seqs)]
    end
    for (s, q, n) in zip(seqs, phreds, names)
        i += 1
        write(stream, FASTQ.Record(n, s, q))
    end
    close(stream)
end


text:
module MadNLPKrylov

import MadNLP:
    @kwdef, Logger, @debug, @warn, @error,
    AbstractOptions, AbstractIterator, set_options!, @sprintf,
    solve_refine!, mul!, ldiv!, size, StrideOneVector
import IterativeSolvers:
    FastHessenberg, ArnoldiDecomp, Residual, init!, init_residual!, expand!, Identity,
    orthogonalize_and_normalize!, update_residual!, gmres_iterable!, GMRESIterable, converged,
    ModifiedGramSchmidt

@kwdef mutable struct Options <: AbstractOptions
    krylov_restart::Int = 5
    krylov_max_iter::Int = 10
    krylov_tol::Float64 = 1e-10
    krylov_acceptable_tol::Float64 = 1e-5
end

struct VirtualMatrix
    mul!::Function
    m::Int
    n::Int
end
mul!(y,A::VirtualMatrix,x) = A.mul!(y,x)
size(A::VirtualMatrix,i) = (A.m,A.n)[i]

struct VirtualPreconditioner
    ldiv!::Function
end
ldiv!(Pl::VirtualPreconditioner,x::StrideOneVector{Float64}) = Pl.ldiv!(x)

mutable struct Solver <: AbstractIterator
    g::Union{Nothing,GMRESIterable}
    res::Vector{Float64}
    opt::Options
    logger::Logger
end

function Solver(res::Vector{Float64},_mul!,_ldiv!;
                opt=Options(),logger=Logger(),
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),kwargs...)
    !isempty(kwargs) && (for (key,val) in kwargs; option_dict[key]=val; end)
    set_options!(opt,option_dict)

    g = GMRESIterable(VirtualPreconditioner(_ldiv!),
                      Identity(),Float64[],Float64[],res,
                      ArnoldiDecomp(VirtualMatrix(_mul!,length(res),length(res)),opt.krylov_restart, Float64),
                      Residual(opt.krylov_restart, Float64),
                      0,opt.krylov_restart,1,
                      opt.krylov_max_iter,opt.krylov_tol,0.,
                      ModifiedGramSchmidt())

    return Solver(g,res,opt,logger)
end

function solve_refine!(x::StridedVector{Float64},
                       is::Solver,
                       b::AbstractVector{Float64})
    @debug(is.logger,"Iterator initiated")
    is.res.=0
    x.=0
    gmres_iterable_update!(is.g,x,b)

    noprogress = 0
    oldres = Inf
    iter = 0
    @debug(is.logger,"iter ||res||")
    @debug(is.logger,@sprintf("%4i %6.2e",iter,is.g.residual.current))
    for (~,res) in enumerate(is.g)
        iter += 1
        mod(iter,10)==0 && @debug(is.logger,"iter ||res||")
        @debug(is.logger,@sprintf("%4i %6.2e",iter,res))
        # res >= oldres && res > is.opt.krylov_acceptable_tol && (noprogress += 1)
        # oldres = res
        # if noprogress >= 3
        #     @debug(is.logger,@sprintf(
        #         "Iterative solver terminated with %4i refinement steps and residual = %6.2e",
        #         iter,res))
        #     return :Singular
        # end
    end
    @debug(is.logger,@sprintf(
        "Iterative solver terminated with %4i refinement steps and residual = %6.2e",iter,is.g.residual.current))

    return (is.g.residual.current < is.opt.krylov_acceptable_tol ? :Solved : :Failed)
end

function gmres_iterable_update!(g,x,b)
    g.x=x
    g.b=b
    g.mv_products = 0
    g.residual.current = init!(g.arnoldi,g.x,g.b,g.Pl,g.Ax,initially_zero = true)
    init_residual!(g.residual, g.residual.current)
    g.k=1
    g.β=g.residual.current
end

end # module


text:
#Schaback convexity preserving rational spline interpolation
using LinearAlgebra
using LeastSquaresOptim
import AQFED.Math: normpdf, inv, ClosedTransformation

export SchabackRationalSpline, makeConvexSchabackRationalSpline, FirstDerivativeBoundary, SecondDerivativeBoundary, FitResult,  fitConvexSchabackRationalSpline

struct SchabackRationalSpline{U<:Number,T<:Real}
    x::Vector{U}
    y::Vector{T}
    h::Vector{T}
    m::Vector{T}
    SchabackRationalSpline(U, T, n::Int) = new{U,T}(zeros(U, n), zeros(T, n), zeros(T, n), zeros(T, n))
    SchabackRationalSpline(x::Vector{U}, y::Vector{T}, h::Vector{T}, m::Vector{T}) where {U<:Real,T<:Real} =
        new{U,T}(x, y, h, m)
end

abstract type PPBoundary end
struct FirstDerivativeBoundary{T} <: PPBoundary
    value::T
end
struct SecondDerivativeBoundary{T} <: PPBoundary
    value::T
end

function initLeftBoundary(b::FirstDerivativeBoundary, delta, d, h, m)
    d[1] = delta[1] - b.value
    m[1] = sqrt(2 * d[1] / (h[1] * m[2])) * sign(d[1])
end

function initLeftBoundary(b::SecondDerivativeBoundary, delta, d, h, m)
    d[1] = b.value
    m[1] = d[1]
end

function initRightBoundary(b::FirstDerivativeBoundary, delta, d, h, m)
    d[end] = b.value - delta[end]
    m[end] = sqrt(2 * d[end] / (h[end] * m[end-1])) * sign(d[end])
end

function initRightBoundary(b::SecondDerivativeBoundary, delta, d, h, m)
    d[end] = b.value
    m[end] = d[end]
end

abstract type SecondDerivativeGuess end
struct Schaback{T} <: SecondDerivativeGuess end
struct Normal{T} <: SecondDerivativeGuess
    μ::T
    σ::T
end

function initSecondDerivative(guessType::Schaback, x, d, h, m)
    for j = 2:length(m)-1
        m[j] = (2 * d[j]) / (h[j] + h[j-1])
    end
end

function initSecondDerivative(guessType::Normal, x, d, h, m)
    for j = 2:length(m)-1
        m[j] = normpdf((x[j] - guessType.μ) / guessType.σ) / guessType.σ
    end
end

struct FitResult
    measure
    iterations
    x
    fx
    details
end
#Newton optimization
function fitConvexSchabackRationalSpline(
    x0::AbstractArray{U},
    y0::AbstractArray{T},
    weights::AbstractArray{T},
    leftBoundary::PPBoundary,
    rightBoundary::PPBoundary;
    penalty = 0.0,
    guessType::SecondDerivativeGuess = Schaback{Float64},
)::Tuple{SchabackRationalSpline{U,T}, FitResult} where {T,U}
    x = copy(x0)
    y = copy(y0)
    h = x[2:end] - x[1:end-1]
    delta = (y[2:end] - y[1:end-1]) ./ h
    d = zeros(T, length(x))
    m = zeros(T, length(x))

    n = length(delta)
    leftBoundaryOverride = typeof(leftBoundary) == SecondDerivativeBoundary{T} ? leftBoundary : SecondDerivativeBoundary((eps(T)))
    rightBoundaryOverride = typeof(rightBoundary) == SecondDerivativeBoundary{T} ? rightBoundary : SecondDerivativeBoundary((eps(T)))
    initLeftBoundary(leftBoundaryOverride, delta, d, h, m)
    initRightBoundary(rightBoundaryOverride, delta, d, h, m)

    sign1 = sign(d[1])
    for i = 2:n
        d[i] = delta[i] - delta[i-1] #corresponds to D[j]
        if sign(d[i]) != sign1
            d[i] = sign1 * eps(one(d[1]))
        end
    end
    #starting value
    initSecondDerivative(guessType, x, d, h, m)
    @. m = max(m, sqrt(eps(T))) #a too small estimate may create instabilities
    @. m = sign(m) * (sign(m) * m^(1 / 3))
    #tridiagonal system
    lower = zeros(T, n)
    middle = zeros(T, n + 1)
    upper = zeros(T, n)
    b = zeros(T, n + 1)
    for j = 2:n
        lower[j-1] = h[j] / h[j-1]
        upper[j] = one(T)
        middle[j] = -(upper[j] + lower[j-1])
    end
    #boundary for linear system can only be second derivative
    middle[1] = one(T)
    upper[1] = zero(T)
    middle[end] = one(T)
    lower[end] = zero(T)
    tri = Tridiagonal(lower, middle, upper)
    transform = ExpMinTransformation(eps(T)^(1/3))
    iter = 0
    #m longer if first derivative boundary ?!?
    function obj!(fvec, x::AbstractArray{W}) where {W}
        #c0 and cn are exact, for first der boundary, use eqn 5.144 5.145
        #to determine line i=1 and i=n-1 => u appears on line 2
        ml = zeros(W,n+1)
        b = zeros(W, n + 1)
        @. ml[2:end-1] = transform(x)
        ml[1] = m[1]
        ml[end] = m[end]
        b[1] = y[1]
        b[n+1] = y[n+1]
        for j = 2:n
            b[j] = h[j]*(h[j-1] * ml[j-1] + h[j] * ml[j+1] ) * ml[j]^2  / 2
        end
        yHat = tri \ b
        @. fvec[1:n-1] = weights[2:end-1] * (yHat[2:end-1] - y[2:end-1]) #end points are exact
        if penalty > 0
            #Strain on M
            for j = 1:n-1
                 s2,s1,s0 = ml[j+2]^3 , ml[j+1]^3 ,  ml[j]^3
                 g2 = 2 * ((s2 - s1) / h[j+1] - (s1 - s0) / h[j]) / (h[j+1] + h[j])
                 #g2 = (s2 - s1) /  h[j+1]
                 fvec[j+n-1] = g2 * sqrt(h[j])  * penalty #strain is sometimes unstable, second derivative is better
            end
        end
        iter += 1
        fvec
    end
    c0 = @. inv(transform, m[2:end-1])
    outlen = n - 1
    if penalty > 0
        outlen += n - 1
    end
    #optimization not great for very small option prices due to machine epsilon issues in Schaback eval combined with large inverse vega weights
    fit = optimize!(
        LeastSquaresProblem(
            x = c0,
            f! = obj!,
            autodiff = :forward,
            # g! = jac!, #useful to debug issue with ForwardDiff NaNs
            output_length = outlen,
        ),
        LevenbergMarquardt();
        iterations = 1024*4,
        x_tol = T(1e-3),
        f_tol = eps(T),
        g_tol = eps(T)
    )
    c0 = fit.minimizer
    fvec = zeros(T, outlen)
    obj!(fvec, c0)
    @. y[2:end-1] += fvec[1:n-1] / weights[2:end-1]
    @. m[2:end-1] = transform(c0)
    #println(iter, " Schaback fit ", fit, fvec) #obj(fit.x))  #fit.f
    return SchabackRationalSpline(x, y, h, m), FitResult(fit.ssr, iter, c0, fvec, fit)
end

#Fixed point iteration
function makeConvexSchabackRationalSpline(
    x0::AbstractArray{U},
    y0::AbstractArray{T},
    leftBoundary::PPBoundary,
    rightBoundary::PPBoundary;
    iterations::Int = 8,
) where {T,U}
    x = copy(x0)
    y = copy(y0)
    h = x[2:end] - x[1:end-1]
    delta = (y[2:end] - y[1:end-1]) ./ h
    d = zeros(T, length(x))
    m = zeros(T, length(x))

    n = length(delta)
    initLeftBoundary(leftBoundary, delta, d, h, m)
    initRightBoundary(rightBoundary, delta, d, h, m)

    sign1 = sign(d[1])
    for i = 2:n
        d[i] = delta[i] - delta[i-1] #corresponds to D[j]
        if sign(d[i]) != sign1
            d[i] = sign1 * eps(one(d[1]))
        end
    end
    #starting value
    for j = 2:n
        m[j] = (2 * d[j]) / (h[j] + h[j-1])
    end
    @. m = sign(m) * (sign(m) * m^(1 / 3))
    #Gauss Seidel iterations
    for loop = 1:iterations
        initLeftBoundary(leftBoundary, delta, d, h, m)
        initRightBoundary(rightBoundary, delta, d, h, m)
        for j = 2:n
            m[j] = sqrt(2 * d[j] / (h[j-1] * m[j-1] + h[j] * m[j+1])) * sign(d[j])
        end
    end
    return SchabackRationalSpline(x, y, h, m)
end

(p::SchabackRationalSpline{U,T})(x::U) where {U,T} = evaluate(p, x)
Base.length(p::SchabackRationalSpline) = length(p.x)
Base.broadcastable(p::SchabackRationalSpline) = Ref(p)

function evaluateRational(s::SchabackRationalSpline{U,T}, index::Int, z::U) where {U,T}
    hi = (s.x[index+1] - s.x[index])
    h = (s.x[index+1] - z)
    m11 = s.m[index+1]
    if abs(m11) <= zero(U)
        return s.y[index+1] - h * ((s.y[index+1] - s.y[index]) / hi)
    end
    m12 = m11 * m11
    m13 = m12 * m11
    mr = (s.m[index] - m11) / (hi * s.m[index])
    return s.y[index+1] - h * ((s.y[index+1] - s.y[index]) / hi + hi / 2 * s.m[index] * m12) +
           (m13 * h^2) / (2 * (1 - h * mr))
end

function evaluateRationalDerivative(s::SchabackRationalSpline{U,T}, index::Int, z::U) where {U,T}
    hi = (s.x[index+1] - s.x[index])
    h = (s.x[index+1] - z)
    m11 = s.m[index+1]
    if abs(m11) <= zero(U)
        return (s.y[index+1] - s.y[index]) / hi
    end
    m12 = m11 * m11
    m13 = m12 * m11
    mr = (s.m[index] - m11)
    mhi = (hi * s.m[index])
    return (s.y[index+1] - s.y[index]) / hi +
           hi / 2 * s.m[index] * m12 +
           m13 / 2 * h * mhi * (-2 - mr * h / (mhi - mr * h)) / (mhi - mr * h)
end

function evaluateRationalSecondDerivative(s::SchabackRationalSpline{U,T}, index::Int, z::U) where {U,T}
    hi = (s.x[index+1] - s.x[index])
    h = (s.x[index+1] - z)
    m11 = s.m[index+1]
    if abs(m11) <= zero(U)
        return zero(T)
    end
    mr = (s.m[index] - m11)
    mhi = (hi * s.m[index])
    return m11^3 * ( mhi + mr * h * mhi * (2 +  mr * h / (mhi - mr * h)) / (mhi - mr * h)) / (mhi - mr * h)
end

function locateLowerIndex(self::SchabackRationalSpline{T,U}, z::T) where {T,U}
    i = searchsortedfirst(self.x, z)  # x[i-1]<z<=x[i]
    if i >= length(self.x)
        return length(self.x) - 1
    elseif i > 1
        return i - 1
    end
end

function evaluate(self::SchabackRationalSpline{T,U}, z::T) where {T,U}
    i = locateLowerIndex(self, z)
    return evaluateRational(self, i, z)
end

function evaluateDerivative(self::SchabackRationalSpline{T,U}, z::T) where {T,U}
    i = locateLowerIndex(self, z)
    return evaluateRationalDerivative(self, i, z)
end

function evaluateSecondDerivative(self::SchabackRationalSpline{T,U}, z::T) where {T,U}
    i = locateLowerIndex(self, z)
    return evaluateRationalSecondDerivative(self, i, z)
end


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

using Test

include("StudentData.jl")
using .StudentData

# Binary Search - Returns the index of the element if found,
# otherwise returns -1
function binarySearch(arr::Vector, elem;
    desc = true,
    compareFn = desc ? (x, y) -> x > y : (x, y) -> x < y,
    startindex = 1)
    mid = ceil(Int, length(arr) / 2)
    if mid == 0
        -1
    elseif arr[mid] == elem
        startindex + mid - 1
    elseif compareFn(arr[mid], elem)
        binarySearch(arr[mid+1:end], elem,
            desc = desc, compareFn = compareFn, startindex = startindex + mid)
    else
        binarySearch(arr[1:mid-1], elem,
            desc = desc, compareFn = compareFn, startindex = startindex)
    end
end

@testset "Binary search" begin
    @test binarySearch([10, 9, 8, 7, 6, 5, 4, 3, 2, 1], 7) == 4
    @test binarySearch([10, 9, 8, 7, 6, 5, 4, 3, 2, 1], 10) == 1
    @test binarySearch([10, 9, 8, 7, 6, 5, 4, 3, 2, 1], 1) == 10
    @test binarySearch([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 7, desc = false) == 7
    @test binarySearch([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 10, desc = false) == 10
    @test binarySearch([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 1, desc = false) == 1
    @test binarySearch(sortedStudents, Student("mark", "s2020127", 7.9),
        compareFn = studentCompareFn) == 3
end

max_stars_repo_name: codekomali/JuliaPractice

max_issues_repo_name: codekomali/JuliaPractice

max_forks_repo_name: codekomali/JuliaPractice


text:
module LatticeQCD
    using Requires

    #include("./SLMC/logdet.jl")

    include("./rhmc/AlgRemez.jl")
    include("./rhmc/rhmc.jl")

    include("./gaugefields/SUN_generator.jl")
    include("./output/verboseprint.jl")
    include("./fermions/cgmethod.jl")

    include("./autostaples/wilsonloops.jl")

    include("./system/system_parameters.jl")
    #include("parallel.jl")
    #include("site.jl")
    include("./system/rand.jl")
    include("./actions/actions.jl")
    include("./gaugefields/gaugefields.jl")

    #include("gaugefields.jl")
    include("./fermions/fermionfields.jl")
    include("./liealgebra/liealgebrafields.jl")

    include("./rationalapprox/rationalapprox.jl")

    include("./fermions/clover.jl")

    include("./fermions/diracoperator.jl")
    include("./fermions/misc.jl")

    include("./output/io.jl")
    include("./output/ildg_format.jl")

    include("./system/LTK_universe.jl")
    include("./gaugefields/smearing.jl")

    include("./output/print_config.jl")

    #include("cg.jl")

    include("./measurements/measurements.jl")
    include("./heatbath/heatbath.jl")
    include("./md/md.jl")
    include("./system/wizard.jl")

    include("./SLMC/SLMC.jl")

    include("./system/mainrun.jl")
    include("./output/analyze.jl")

    function __init__()
        @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
            include("./demo/demo.jl")
            import .Demo:demo
            export demo

            import .Analyze:plot_plaquette,plot_polyakov,plot_plaq_and_poly
            export plot_plaquette,plot_polyakov,plot_plaq_and_poly
        end
    end

    import .LTK_universe:Universe,show_parameters,make_WdagWmatrix,calc_Action,set_β!,set_βs!
    import .Actions:Setup_Gauge_action,Setup_Fermi_action,GaugeActionParam_autogenerator
    import .Measurements:calc_plaquette,measure_correlator,Measurement,calc_polyakovloop,measure_chiral_cond,calc_topological_charge,
                measurements,Measurement_set
    import  .MD:md_initialize!,MD_parameters_standard,md!,metropolis_update!,construct_MD_parameters
    import .System_parameters:Params,print_parameters,parameterloading,Params_set,make_parameters#,parameterloading2
    import .Print_config:write_config
    import .Smearing:gradientflow!
    import .ILDG_format:ILDG,load_gaugefield
    import .Heatbath:heatbath!
    import .Wilsonloops:make_plaq
    import .IOmodule:saveU,loadU,loadU!
    import .Wizard:run_wizard
    import .Mainrun:run_LQCD
    import .RationalApprox:calc_exactvalue,calc_Anϕ,calc_det
    #,run_LQCD!

    import .Analyze:analyze,get_plaquette,get_polyakov,get_plaquette_average,get_polyakov_average,get_trjs

    #import .Fermionfields:make_WdagWmatrix

    export Setup_Gauge_action,Setup_Fermi_action,GaugeActionParam_autogenerator
    export Universe,set_β!,set_βs!
    export calc_plaquette,calc_polyakovloop,calc_topological_charge
    export md_initialize!,MD_parameters_standard,md!,metropolis_update!,construct_MD_parameters
    export show_parameters
    export Params,print_parameters,parameterloading,Params_set,make_parameters#,parameterloading2
    export measure_correlator,measure_chiral_cond,Measurement,measurements,Measurement_set
    export gradientflow!
    export ILDG,load_gaugefield
    export make_WdagWmatrix
    export heatbath!
    export make_plaq
    export calc_Action
    export calc_topological_charge
    export saveU,loadU,loadU!
    export run_LQCD,run_LQCD!

    export write_config
    export run_wizard
    export analyze,get_plaquette,get_polyakov,get_plaquette_average,get_polyakov_average,get_trjs

end


text:
"""
    NormInfinityBridge{T}

The `NormInfinityCone` is representable with LP constraints, since
``t \\ge \\max_i \\lvert x_i \\rvert`` if and only if
``t \\ge x_i`` and ``t \\ge -x_i`` for all ``i``.
"""
struct NormInfinityBridge{T, F, G} <: SetMapBridge{T, MOI.Nonnegatives, MOI.NormInfinityCone, F, G}
    constraint::CI{F, MOI.Nonnegatives}
end
function concrete_bridge_type(::Type{<:NormInfinityBridge{T}}, G::Type{<:MOI.AbstractVectorFunction}, ::Type{MOI.NormInfinityCone}) where T
    F = MOIU.promote_operation(+, T, G, G)
    return NormInfinityBridge{T, F, G}
end

function map_set(::Type{<:NormInfinityBridge}, set::MOI.NormInfinityCone)
    return MOI.Nonnegatives(2 * (MOI.dimension(set) - 1))
end
function inverse_map_set(::Type{<:NormInfinityBridge}, set::MOI.Nonnegatives)
    return MOI.NormInfinityCone(div(MOI.dimension(set), 2) + 1)
end

function map_function(::Type{<:NormInfinityBridge{T}}, func) where T
    scalars = MOIU.eachscalar(func)
    t = scalars[1]
    lb = scalars[2:end]
    ub = MOIU.operate(-, T, lb)
    f_new = MOIU.operate(vcat, T, ub, lb)
    for i in 1:(2 * (length(scalars) - 1))
        MOIU.operate_output_index!(+, T, i, f_new, t)
    end
    return f_new
end
function inverse_map_function(::Type{<:NormInfinityBridge{T}}, func) where T
    scalars = MOIU.eachscalar(func)
    t = MOIU.operate!(/, T, sum(scalars), T(length(scalars)))
    d = div(length(scalars), 2)
    x = MOIU.operate!(/, T, MOIU.operate!(-, T, scalars[(d + 1):end], scalars[1:d]), T(2))
    return MOIU.operate(vcat, T, t, x)
end
# Given a_i is dual on t - x_i >= 0 and b_i is dual on t + x_i >= 0,
# the dual on (t, x) in NormInfinityCone is (u, v) in NormOneCone, where
# v_i = -a_i + b_i and u = sum(a) + sum(b).
function adjoint_map_function(::Type{<:NormInfinityBridge}, func)
    scalars = MOIU.eachscalar(func)
    t = sum(scalars)
    d = div(length(scalars), 2)
    x = (scalars[(d + 1):end] - scalars[1:d])
    return vcat(t, x)
end
function inverse_adjoint_map_function(::Type{<:NormInfinityBridge{T}}, func::AbstractVector{T}) where T
    # This is used by `MOI.ConstraintDualStart`.
    # The result should belong to `MOI.Nonnegatives` and the sum of the elements should
    # be `t`.
    # Only one of dual of the upper and lower bound should be active so one of the two
    # duals is zero. We know which one since they should be nonnegative.
    # Then if `t = sum abs(x_i)`, we will indeed have only one of them being zero.
    t = func[1]
    y = func[2:end]
    lb = [x > 0 ?  x : zero(x) for x in y]
    ub = [x < 0 ? -x : zero(x) for x in y]
    x = [ub; lb]
    return x .+ (t - sum(x)) / length(x)
end

"""
    NormOneBridge{T}

The `NormOneCone` is representable with LP constraints, since
``t \\ge \\sum_i \\lvert x_i \\rvert`` if and only if there exists a vector y such that
``t \\ge \\sum_i y_i`` and ``y_i \\ge x_i``, ``y_i \\ge -x_i`` for all ``i``.
"""
struct NormOneBridge{T, F, G, H} <: AbstractBridge
    y::Vector{MOI.VariableIndex}
    ge_index::CI{F, MOI.GreaterThan{T}}
    nn_index::CI{G, MOI.Nonnegatives}
end
function bridge_constraint(::Type{NormOneBridge{T, F, G, H}}, model::MOI.ModelLike, f::H, s::MOI.NormOneCone) where {T, F, G, H}
    f_scalars = MOIU.eachscalar(f)
    d = MOI.dimension(s)
    y = MOI.add_variables(model, d - 1)
    ge_index = MOIU.normalize_and_add_constraint(model, MOIU.operate(-, T, f_scalars[1], MOIU.operate(sum, T, y)), MOI.GreaterThan(zero(T)), allow_modify_function=true)
    lb = f_scalars[2:d]
    ub = MOIU.operate(-, T, lb)
    lb = MOIU.operate!(+, T, lb, MOI.VectorOfVariables(y))
    ub = MOIU.operate!(+, T, ub, MOI.VectorOfVariables(y))
    f_new = MOIU.operate(vcat, T, ub, lb)
    nn_index = MOI.add_constraint(model, f_new, MOI.Nonnegatives(2d - 2))
    return NormOneBridge{T, F, G, H}(y, ge_index, nn_index)
end

MOI.supports_constraint(::Type{NormOneBridge{T}}, ::Type{<:MOI.AbstractVectorFunction}, ::Type{MOI.NormOneCone}) where T = true
MOIB.added_constrained_variable_types(::Type{<:NormOneBridge}) = Tuple{DataType}[]
MOIB.added_constraint_types(::Type{NormOneBridge{T, F, G, H}}) where {T, F, G, H} = [(F, MOI.GreaterThan{T}), (G, MOI.Nonnegatives)]
function concrete_bridge_type(::Type{<:NormOneBridge{T}}, H::Type{<:MOI.AbstractVectorFunction}, ::Type{MOI.NormOneCone}) where T
    S = MOIU.scalar_type(H)
    F = MOIU.promote_operation(+, T, S, S)
    G = MOIU.promote_operation(+, T, H, H)
    return NormOneBridge{T, F, G, H}
end

# Attributes, Bridge acting as a model
MOI.get(b::NormOneBridge, ::MOI.NumberOfVariables) = length(b.y)
MOI.get(b::NormOneBridge, ::MOI.ListOfVariableIndices) = b.y
MOI.get(b::NormOneBridge{T, F, G, H}, ::MOI.NumberOfConstraints{F, MOI.GreaterThan{T}}) where {T, F, G, H} = 1
MOI.get(b::NormOneBridge{T, F, G, H}, ::MOI.NumberOfConstraints{G, MOI.Nonnegatives}) where {T, F, G, H} = 1
MOI.get(b::NormOneBridge{T, F, G, H}, ::MOI.ListOfConstraintIndices{F, MOI.GreaterThan{T}}) where {T, F, G, H} = [b.ge_index]
MOI.get(b::NormOneBridge{T, F, G, H}, ::MOI.ListOfConstraintIndices{G, MOI.Nonnegatives}) where {T, F, G, H} = [b.nn_index]

# References
function MOI.delete(model::MOI.ModelLike, c::NormOneBridge)
    MOI.delete(model, c.nn_index)
    MOI.delete(model, c.ge_index)
    MOI.delete(model, c.y)
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintFunction, c::NormOneBridge{T, F, G, H}) where {T, F, G, H}
    ge_func = MOI.get(model, MOI.ConstraintFunction(), c.ge_index)
    nn_func = MOIU.eachscalar(MOI.get(model, MOI.ConstraintFunction(), c.nn_index))
    t = MOIU.operate!(+, T, ge_func, MOIU.operate!(/, T, sum(nn_func), T(2)))
    d = div(length(nn_func), 2)
    x = MOIU.operate!(/, T, MOIU.operate!(-, T, nn_func[(d + 1):end], nn_func[1:d]), T(2))
    return MOIU.convert_approx(H, MOIU.remove_variable(MOIU.operate(vcat, T, t, x), c.y))
end
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintSet, c::NormOneBridge)
    dim = 1 + div(MOI.dimension(MOI.get(model, MOI.ConstraintSet(), c.nn_index)), 2)
    return MOI.NormOneCone(dim)
end

function MOI.supports(
    ::MOI.ModelLike,
    ::Union{MOI.ConstraintPrimalStart, MOI.ConstraintDualStart},
    ::Type{<:NormOneBridge})

    return true
end
function MOI.set(model::MOI.ModelLike, attr::MOI.ConstraintPrimalStart,
                 bridge::NormOneBridge{T}, value) where T
    x_value = value[1 .+ (1:length(bridge.y))]
    y_value = abs.(x_value)
    for i in eachindex(bridge.y)
        MOI.set(model, MOI.VariablePrimalStart(), bridge.y[i], y_value[i])
    end
    MOI.set(model, attr, bridge.nn_index, [y_value - x_value; y_value + x_value])
    MOI.set(model, attr, bridge.ge_index, value[1] - reduce(+, y_value, init=zero(T)))
    return
end
function MOI.get(model::MOI.ModelLike,
                 attr::Union{MOI.ConstraintPrimal, MOI.ConstraintPrimalStart},
                 bridge::NormOneBridge)
    ge_primal = MOI.get(model, attr, bridge.ge_index)
    nn_primal = MOI.get(model, attr, bridge.nn_index)
    t = ge_primal + sum(nn_primal) / 2
    d = length(bridge.y)
    x = (nn_primal[(d + 1):end] - nn_primal[1:d]) / 2
    return vcat(t, x)
end
# Given a_i is dual on y_i - x_i >= 0 and b_i is dual on y_i + x_i >= 0 and c is dual on t - sum(y) >= 0,
# the dual on (t, x) in NormOneCone is (u, v) in NormInfinityCone, where
# v_i = -a_i + b_i and u = c.
function MOI.get(model::MOI.ModelLike,
                 attr::Union{MOI.ConstraintDual, MOI.ConstraintDualStart},
                 bridge::NormOneBridge)
    t = MOI.get(model, attr, bridge.ge_index)
    nn_dual = MOI.get(model, attr, bridge.nn_index)
    d = length(bridge.y)
    x = nn_dual[(d + 1):end] - nn_dual[1:d]
    return vcat(t, x)
end
# value[1 + i] = nn_dual[d + i] - nn_dual[i]
# and `nn_dual` is nonnegative. By complementarity slackness, only one of each
# `nn_dual` can be nonzero (except if `x = 0`) so we can set
# depending on the sense of `value[1 + i]`.
function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintDualStart,
                 bridge::NormOneBridge, value)
    t = MOI.set(model, MOI.ConstraintDualStart(), bridge.ge_index, value[1])
    d = length(bridge.y)
    nn_dual = zeros(eltype(value), 2d)
    for i in eachindex(bridge.y)
        v = value[1 + i]
        if v < 0
            nn_dual[i] = -v
        else
            nn_dual[d + i] = v
        end
    end
    MOI.set(model, MOI.ConstraintDualStart(), bridge.nn_index, nn_dual)
    return
end


text: using Test, QuHamiltonian

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
@testset "Neighbor List Unit Tests" begin
    atoms = [
        Atom(:Ar, (@SVector [0.5, 0.1, 0.3])u"bohr"),
        Atom(:Ar, (@SVector [0.7, 2.5, -2.1])u"bohr"),
        Atom(:Ar, (@SVector [1.1, -3.1, 7.3])u"bohr")
    ]
    box = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]u"bohr"
    boundary_conditions = [Periodic(), DirichletZero(), Periodic()]
    system = FlexibleSystem(atoms, box, boundary_conditions)

    rcutoff = 5.0
    nnlist = neighborlist(system, rcutoff)

    true_j = [
        Int[2, 3],
        Int[],
        Int[]
    ]

    true_R = [
        Float64[norm([-0.2, -2.4, 0.4]), norm([0.4, 3.2, 0.0])],
        Float64[],
        Float64[]
    ]

    true_r = [
        Vector{Float64}[[-0.2, -2.4, 0.4], [0.4, 3.2, 0.0]],
        Vector{Float64}[],
        Vector{Float64}[]
    ]

    @test length(nnlist) == length(system)
    @test nnlist.j == true_j
    @test all(nnlist.R .≈ true_R)
    @test all(all(ri .≈ true_ri) for (ri, true_ri) in zip(nnlist.r, true_r))
end


text:
using LinearAlgebra

"""
    DBAResult(cost,converged,iterations,cost_trace)

Holds results of a DTW Barycenter Averaging (DBA) fit.
"""
mutable struct DBAResult
    cost::Float64
    converged::Bool
    iterations::Int
    cost_trace::Vector{Float64}
end

"""
    avgseq, results = dba(sequences, [dist=SqEuclidean()]; kwargs...)

Perfoms DTW Barycenter Averaging (DBA) given a collection of `sequences`
and the current estimate of the average sequence.

Example usage:

    x = Sequence([1., 2., 2., 3., 3., 4.])
    y = Sequence([1., 3., 4.])
    z = Sequence([1., 2., 2., 4.])
    avg,result = dba([x,y,z])
"""
function dba(
        sequences::AbstractVector{T},
        method::DTWMethod = ClassicDTW(),
        dist::SemiMetric = SqEuclidean();
        init_center::T = rand(sequences),
        iterations::Int = 1000,
        rtol::Float64 = 1e-5,
        store_trace::Bool = false,
        show_progress::Bool = true,
        i2min::AbstractVector =[],
        i2max::AbstractVector = []
    ) where {T<:Sequence}

    # method for computing dtw
    dtwdist = DTWDistance(method,dist)

    # initialize dbavg as a random sample from the dataset
    nseq = length(sequences)
    dbavg = deepcopy(init_center)

    # storage for each iteration
    newavg = Sequence(zeros(size(dbavg.val)))
    counts = zeros(Int, length(dbavg))

    # variables storing optimization progress
    converged = false
    iter = 0
    cost,newcost = Inf,Inf
    cost_trace = Float64[]

    # display optimization progress
    if show_progress
        p = ProgressMeter.ProgressThresh(rtol)
    end

    ## main loop ##
    while !converged && iter < iterations

        # do an iteration of dba
        newcost = dba_iteration!(newavg, dbavg, counts, sequences, dtwdist; i2min=i2min,i2max=i2max)
        iter += 1

        # store history of cost while optimizing (optional)
        store_trace && push!(cost_trace, newcost)

        # check convergence
        Δ = (cost-newcost)/newcost
        if Δ < rtol
            converged = true
        else
            # update estimate
            cost = newcost
            dbavg = deepcopy(newavg)
        end

        # update progress bar
        if show_progress
            ProgressMeter.update!(p, Δ; showvalues =[(:iteration,iter),
                                                 (Symbol("max iteration"),iterations),
                                                 (:cost,cost)])
        end
    end

    return newavg, DBAResult(newcost,converged,iter,cost_trace)
end

"""
    newavg, cost = dba_iteration(dbavg, sequences, dist)

Performs one iteration of DTW Barycenter Averaging (DBA) given a collection of
`sequences` and the current estimate of the average sequence, `dbavg`. Returns
an updated estimate, and the cost/loss of the previous estimate
"""
function dba_iteration!(
        newavg::T,
        oldavg::T,
        counts::Array{Int,1},
        sequences::AbstractVector{T},
        d::DTWDistance;
        i2min::AbstractVector=[],
        i2max::AbstractVector=[]
    ) where {T<:Sequence}

    # sum of dtw dist of all sequences to center
    total_cost = 0.0

    # store stats for barycenter averages
    rmul!(counts,0)
    rmul!(newavg,0)

    # main ploop
    for seq in sequences
        # time warp signal versus average
        # if one of the two is empty, use unconstrained window. If both are nonempty, but not the same length, distpath will throw error
        if isempty(i2min) && isempty(i2max)
          cost, i1, i2 = distpath(d, oldavg, seq)
        else
          cost, i1, i2 = distpath(d, oldavg, seq, i2min, i2max)
        end
        total_cost += cost

        # store stats for barycentric average
        for j=1:length(i2)
            counts[i1[j]] += 1
            newavg[i1[j]] += seq[i2[j]]
        end
    end

    # compute average and return total cost
    for i in eachindex(newavg)
        newavg[i] = newavg[i]/counts[i]
    end

    return total_cost
end

# Wrapper for AbstractArray of one-dimensional time series.
function dba( s::AbstractArray, args...; kwargs... )
    dba(_sequentize(s), args...; kwargs...)
end

# weirdly enought, this works for the dtw_dba_miniexample,  but does not work for dtw_dbaclust
#@generated function _sequentize{T,N}(s::AbstractArray{T,N})
#    :( Sequence[ Sequence(@ncall($N, view, s, n-> n==$N ? i : Colon())) for i = 1:size(s,2) ] )
#end

# This function replaces the compile-time function which doesn't work in v0.6
# Compile-time function fixed, generated functions to not support untyped generators
function _sequentize(s::AbstractArray)
    [Sequence(s[:,i]) for i=1:size(s,2)]
end


text:
function handle_callbacks!(integrator)
  discrete_callbacks = integrator.opts.callback.discrete_callbacks
  continuous_callbacks = integrator.opts.callback.continuous_callbacks
  atleast_one_callback = false

  continuous_modified = false
  discrete_modified = false
  saved_in_cb = false
  if !(typeof(continuous_callbacks)<:Tuple{})
    time,upcrossing,event_occured,idx,counter =
              find_first_continuous_callback(integrator,continuous_callbacks...)
    if event_occured
      integrator.event_last_time = true
      continuous_modified,saved_in_cb = apply_callback!(integrator,continuous_callbacks[idx],time,upcrossing)
    else
      integrator.event_last_time = false
    end
  end
  if !(typeof(discrete_callbacks)<:Tuple{})
    discrete_modified,saved_in_cb = apply_discrete_callback!(integrator,discrete_callbacks...)
  end
  if !saved_in_cb
    savevalues!(integrator)
  end

  integrator.u_modified = continuous_modified || discrete_modified
  if integrator.u_modified
    handle_callback_modifiers!(integrator)
  end
  integrator.u_modified = false
end

function DiffEqBase.savevalues!(integrator::AbstractSundialsIntegrator,force_save=false)
    uType = eltype(integrator.sol.u)
    while !isempty(integrator.opts.saveat) &&
        integrator.tdir*top(integrator.opts.saveat) < integrator.tdir*first(integrator.tout)

        curt = pop!(integrator.opts.saveat)
        tmp = integrator(curt)
        save_value!(integrator.sol.u,tmp,uType,integrator.sizeu,Val{false})
        push!(integrator.sol.t,curt)
        if integrator.opts.dense
          tmp = integrator(curt,Val{1})
          save_value!(integrator.sol.interp.du,tmp,uType,integrator.sizeu,Val{false})
        end
    end

    if integrator.opts.save_everystep || force_save
        save_value!(integrator.sol.u,integrator.u,uType,integrator.sizeu)
        push!(integrator.sol.t, integrator.t)
        if integrator.opts.dense
          tmp = integrator(integrator.t,Val{1})
          save_value!(integrator.sol.interp.du,tmp,uType,integrator.sizeu)
        end
    end
end

function save_value!(save_array,val,::Type{T},sizeu,
                     make_copy::Type{Val{bool}}=Val{true}) where {T <: Number,bool}
    push!(save_array,first(val))
end
function save_value!(save_array,val,::Type{T},sizeu,
                     make_copy::Type{Val{bool}}=Val{true}) where {T <: Vector,bool}
    bool ? save = copy(val) : save = val
    push!(save_array,save)
end
function save_value!(save_array,val,::Type{T},sizeu,
                     make_copy::Type{Val{bool}}=Val{true}) where {T <: Array,bool}
    bool ? save = copy(val) : save = val
    push!(save_array,reshape(save,sizeu))
end

function handle_callback_modifiers!(integrator::CVODEIntegrator)
  CVodeReInit(integrator.mem,integrator.t,integrator.u)
end

function handle_callback_modifiers!(integrator::IDAIntegrator)
  IDAReInit(integrator.mem,integrator.t,integrator.u,integrator.du)
end

@inline function DiffEqBase.add_tstop!(integrator::AbstractSundialsIntegrator,t)
  t < integrator.t && error("Tried to add a tstop that is behind the current time. This is strictly forbidden")
  push!(integrator.opts.tstops,t)
end

function DiffEqBase.add_saveat!(integrator::AbstractSundialsIntegrator,t)
  integrator.tdir * (t - integrator.t) < 0 && error("Tried to add a saveat that is behind the current time. This is strictly forbidden")
  push!(integrator.opts.saveat,t)
end

@inline DiffEqBase.get_tmp_cache(integrator::AbstractSundialsIntegrator) = (integrator.tmp,)

@inline function DiffEqBase.u_modified!(integrator::AbstractSundialsIntegrator,bool::Bool)
  integrator.u_modified = bool
end

@inline function DiffEqBase.terminate!(integrator::AbstractSundialsIntegrator)
  integrator.opts.tstops.valtree = typeof(integrator.opts.tstops.valtree)()
end

@inline function DiffEqBase.get_du(integrator::CVODEIntegrator)
  integrator(integrator.t,Val{1})
end

@inline function DiffEqBase.get_du!(out,integrator::CVODEIntegrator)
  integrator(out,integrator.t,Val{1})
end

@inline function DiffEqBase.get_du(integrator::IDAIntegrator)
  reshape(integrator.du,integrator.sizedu)
end

@inline function DiffEqBase.get_du!(out,integrator::IDAIntegrator)
  out .= reshape(integrator.du,integrator.sizedu)
end


text:
# ---
# title: 297. Serialize and Deserialize Binary Tree
# id: problem297
# author: zhwang
# date: 2022-03-11
# difficulty: Hard
# categories: Tree, Design
# link: <https://leetcode.com/problems/serialize-and-deserialize-binary-tree/description/>
# hidden: true
# ---
#
# Serialization is the process of converting a data structure or object into a
# sequence of bits so that it can be stored in a file or memory buffer, or
# transmitted across a network connection link to be reconstructed later in the
# same or another computer environment.
#
# Design an algorithm to serialize and deserialize a binary tree. There is no
# restriction on how your serialization/deserialization algorithm should work.
# You just need to ensure that a binary tree can be serialized to a string and
# this string can be deserialized to the original tree structure.
#
# **Clarification:** The input/output format is the same as [how LeetCode
# serializes a binary tree](https://leetcode.com/faq/#binary-tree). You do not necessarily need to
# follow this format, so please be creative and come up with different
# approaches yourself.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/15/serdeser.jpg)
#
#
#
#     Input: root = [1,2,3,null,null,4,5]
#     Output: [1,2,3,null,null,4,5]
#
#
# **Example 2:**
#
#
#
#     Input: root = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: root = [1]
#     Output: [1]
#
#
# **Example 4:**
#
#
#
#     Input: root = [1,2]
#     Output: [1,2]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[0, 104]`.
#   * `-1000 <= Node.val <= 1000`
#
#
## @lc code=start
using LeetCode

serialize(::Nothing) = "[]"
function serialize(root::TreeNode{Int})::String
    res, queue, hasnew = String[], Union{TreeNode,Nothing}[root], true
    while hasnew
        hasnew = false
        for _ in 1:length(queue)
            node = popfirst!(queue)
            if !isnothing(node)
                push!(res, string(node.val))
                push!(queue, node.left)
                push!(queue, node.right)
                all(isnothing.([node.left, node.right])) || (hasnew = true)
            else
                push!(res, "null")
            end
        end
    end
    return "[" * join(res, ',') * "]"
end

function deserialize(data::String)::Union{TreeNode,Nothing}
    data == "[]" && return nothing
    vals = split(data[2:(end - 1)], ',')
    root, n = TreeNode(parse(Int, vals[1])), length(vals)
    queue, m = [root], 2
    while m <= n
        for _ in 1:length(queue) ## new nodes
            node = popfirst!(queue)
            if vals[m] != "null"
                node.left = TreeNode(parse(Int, vals[m]))
                push!(queue, node.left)
            end
            if vals[m + 1] != "null"
                node.right = TreeNode(parse(Int, vals[m + 1]))
                push!(queue, node.right)
            end
            m += 2
        end
    end
    return root
end

## @lc code=end