text:
# Julia wrapper for header: /usr/include/scip/scip_nonlinear.h
# Automatically generated using Clang.jl wrap_c

function SCIPaddSquareLinearization(scip, sqrcoef, refpoint, isint, lincoef, linconstant, success)
    ccall((:SCIPaddSquareLinearization, libscip), Cvoid, (Ptr{SCIP_}, Cdouble, Cdouble, UInt32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}), scip, sqrcoef, refpoint, isint, lincoef, linconstant, success)
end

function SCIPaddSquareSecant(scip, sqrcoef, lb, ub, refpoint, lincoef, linconstant, success)
    ccall((:SCIPaddSquareSecant, libscip), Cvoid, (Ptr{SCIP_}, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}), scip, sqrcoef, lb, ub, refpoint, lincoef, linconstant, success)
end

function SCIPaddBilinLinearization(scip, bilincoef, refpointx, refpointy, lincoefx, lincoefy, linconstant, success)
    ccall((:SCIPaddBilinLinearization, libscip), Cvoid, (Ptr{SCIP_}, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}), scip, bilincoef, refpointx, refpointy, lincoefx, lincoefy, linconstant, success)
end

function SCIPaddBilinMcCormick(scip, bilincoef, lbx, ubx, refpointx, lby, uby, refpointy, overestimate, lincoefx, lincoefy, linconstant, success)
    ccall((:SCIPaddBilinMcCormick, libscip), Cvoid, (Ptr{SCIP_}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, UInt32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}), scip, bilincoef, lbx, ubx, refpointx, lby, uby, refpointy, overestimate, lincoefx, lincoefy, linconstant, success)
end

function SCIPcomputeBilinEnvelope1(scip, bilincoef, lbx, ubx, refpointx, lby, uby, refpointy, overestimate, xcoef, ycoef, constant, lincoefx, lincoefy, linconstant, success)
    ccall((:SCIPcomputeBilinEnvelope1, libscip), Cvoid, (Ptr{SCIP_}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, UInt32, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}), scip, bilincoef, lbx, ubx, refpointx, lby, uby, refpointy, overestimate, xcoef, ycoef, constant, lincoefx, lincoefy, linconstant, success)
end

function SCIPcomputeBilinEnvelope2(scip, bilincoef, lbx, ubx, refpointx, lby, uby, refpointy, overestimate, alpha1, beta1, gamma1, alpha2, beta2, gamma2, lincoefx, lincoefy, linconstant, success)
    ccall((:SCIPcomputeBilinEnvelope2, libscip), Cvoid, (Ptr{SCIP_}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, UInt32, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{UInt32}), scip, bilincoef, lbx, ubx, refpointx, lby, uby, refpointy, overestimate, alpha1, beta1, gamma1, alpha2, beta2, gamma2, lincoefx, lincoefy, linconstant, success)
end

function SCIPcreateNlpiProb(scip, nlpi, nlrows, nnlrows, nlpiprob, var2idx, nlscore, cutoffbound, setobj, onlyconvex)
    ccall((:SCIPcreateNlpiProb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NLPI}, Ptr{Ptr{SCIP_NLROW}}, Cint, Ptr{SCIP_NLPIPROBLEM}, Ptr{SCIP_HASHMAP}, Ptr{Cdouble}, Cdouble, UInt32, UInt32), scip, nlpi, nlrows, nnlrows, nlpiprob, var2idx, nlscore, cutoffbound, setobj, onlyconvex)
end

function SCIPupdateNlpiProb(scip, nlpi, nlpiprob, var2nlpiidx, nlpivars, nlpinvars, cutoffbound)
    ccall((:SCIPupdateNlpiProb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NLPI}, Ptr{SCIP_NLPIPROBLEM}, Ptr{SCIP_HASHMAP}, Ptr{Ptr{SCIP_VAR}}, Cint, Cdouble), scip, nlpi, nlpiprob, var2nlpiidx, nlpivars, nlpinvars, cutoffbound)
end

function SCIPaddNlpiProbRows(scip, nlpi, nlpiprob, var2idx, rows, nrows)
    ccall((:SCIPaddNlpiProbRows, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NLPI}, Ptr{SCIP_NLPIPROBLEM}, Ptr{SCIP_HASHMAP}, Ptr{Ptr{SCIP_ROW}}, Cint), scip, nlpi, nlpiprob, var2idx, rows, nrows)
end


text:
include("zmodel.jl")
_dx = 2.5
_dy = 1.0
_d = Array(Float32,512,512)

dlimit(a::Int,b::Int,c::Int)=(a<b ? b : a>c ? c : a)

function dview{T}(data::Array{T,3}, img::Int; temp::Matrix{T}=_d, dx=_dx, dy=_dy)
    z = mod1(img-Z0,Z)
    dx = round(Int64,dx * z)
    dy = round(Int64,dy * z)
    X,Y = size(temp)
    @inbounds for x=1:X, y=1:Y
        xx = dlimit(x+dx, 1, X)
        yy = dlimit(y+dy, 1, Y)
        temp[x,y] = data[xx,yy,img]
    end
    zview(temp)
end

function dfix{T}(data::Array{T,3}; temp::Matrix{T}=_d, dx=_dx, dy=_dy, out::Array{T,3}=similar(data))
    @assert dx>=0 && dy>=0
    copy!(out,data)
    X,Y,I = size(data)
    for i=1:I
        z = mod1(i-Z0,Z)
        x = round(Int64,dx * z)
        y = round(Int64,dy * z)
        copy!(sub(out,1:X-x,1:Y-y,i),sub(data,x+1:X,y+1:Y,i))
    end
    return out
end

ext: jl

max_stars_repo_name: JuliaTagBot/NeuronTracking

max_issues_repo_name: JuliaTagBot/NeuronTracking

max_forks_repo_name: JuliaTagBot/NeuronTracking


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using JuMP, LinearAlgebra, Geodesy, Cbc, Clp, ProgressBars, Printf, DataStructures

function _get_default_milp_optimizer()
    return optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
end

function _get_default_lp_optimizer()
    return optimizer_with_attributes(Clp.Optimizer, "LogLevel" => 0)
end

function _print_graph_stats(instance::Instance, graph::Graph)::Nothing
    @info @sprintf("    %12d time periods", instance.time)
    @info @sprintf("    %12d process nodes", length(graph.process_nodes))
    @info @sprintf("    %12d shipping nodes (plant)", length(graph.plant_shipping_nodes))
    @info @sprintf(
        "    %12d shipping nodes (collection)",
        length(graph.collection_shipping_nodes)
    )
    @info @sprintf("    %12d arcs", length(graph.arcs))
    return
end

function solve(
    instance::Instance;
    optimizer = nothing,
    output = nothing,
    marginal_costs = true,
    return_model = false,
)

    milp_optimizer = lp_optimizer = optimizer
    if optimizer == nothing
        milp_optimizer = _get_default_milp_optimizer()
        lp_optimizer = _get_default_lp_optimizer()
    end

    @info "Building graph..."
    graph = RELOG.build_graph(instance)
    _print_graph_stats(instance, graph)

    @info "Building optimization model..."
    model = RELOG.build_model(instance, graph, milp_optimizer)

    @info "Optimizing MILP..."
    JuMP.optimize!(model)

    if !has_values(model)
        error("No solution available")
    end

    if marginal_costs
        @info "Re-optimizing with integer variables fixed..."
        all_vars = JuMP.all_variables(model)
        vals = OrderedDict(var => JuMP.value(var) for var in all_vars)
        JuMP.set_optimizer(model, lp_optimizer)
        for var in all_vars
            if JuMP.is_binary(var)
                JuMP.unset_binary(var)
                JuMP.fix(var, vals[var])
            end
        end
        JuMP.optimize!(model)
    end

    @info "Extracting solution..."
    solution = get_solution(model, marginal_costs = marginal_costs)

    if output != nothing
        write(solution, output)
    end

    if return_model
        return solution, model
    else
        return solution
    end
end

function solve(filename::AbstractString; heuristic = false, kwargs...)
    @info "Reading $filename..."
    instance = RELOG.parsefile(filename)
    if heuristic && instance.time > 1
        @info "Solving single-period version..."
        compressed = _compress(instance)
        csol = solve(compressed; output = nothing, marginal_costs = false, kwargs...)
        @info "Filtering candidate locations..."
        selected_pairs = []
        for (plant_name, plant_dict) in csol["Plants"]
            for (location_name, location_dict) in plant_dict
                push!(selected_pairs, (plant_name, location_name))
            end
        end
        filtered_plants = []
        for p in instance.plants
            if (p.plant_name, p.location_name) in selected_pairs
                push!(filtered_plants, p)
            end
        end
        instance.plants = filtered_plants
        @info "Solving original version..."
    end
    sol = solve(instance; kwargs...)
    return sol
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
using SafeTestsets

@safetestset prefixsearch_single = "test single value" begin
    using Fleck: PrefixSearchTree, choose, sum
	t = PrefixSearchTree([3])
	@test sum(t) == 3
	@test choose(t, 2)[1] == 1
	push!(t, [(1, 4)])
	@test sum(t) == 4
	@test choose(t, 3.3)[1] == 1
end

@safetestset prefixsearch_double = "test double value" begin
    using Fleck: PrefixSearchTree, choose, sum
    t = PrefixSearchTree([3, 1])
	@test sum(t) == 4
	@test choose(t, 2)[1] == 1
	@test choose(t, 3)[1] == 2
	@test choose(t, 3.7)[1] == 2
	push!(t, [(1, 4)])
	v = [(2, 1), (3.3, 1), (4.1, 2)]
	for (guess, result) in v
		@test choose(t, guess)[1] == result
	end
	push!(t, [(2, 2)])
	v = [(2, 1), (3.3, 1), (5.1, 2)]
	for (guess, result) in v
		@test choose(t, guess)[1] == result
	end
end

@safetestset prefixsearch_three = "test three values" begin
    using Fleck: PrefixSearchTree, choose, sum
    t = PrefixSearchTree([3, 1, 2])
    @test sum(t) == 6
    @test choose(t, 2.1)[1] == 1
    @test choose(t, 3)[1] == 2
    @test choose(t, 3.1)[1] == 2
    @test choose(t, 5)[1] == 3
    push!(t, ((2, 2), (3, 3)))
    v = [(2, 1), (3, 2), (4.8, 2), (5.1, 3), (7.9, 3)]
    for (guess, result) in v
        @test choose(t, guess)[1] == result
    end
end

@safetestset prefixsearch_three_floats = "three floats" begin
    using Fleck: PrefixSearchTree, choose, sum
    t = PrefixSearchTree([3.5, 1.5, 2.5])
    @test abs(sum(t) - 7.5) < 1e-9
    @test choose(t, 2.1)[1] == 1
    @test choose(t, 3.5)[1] == 2
    @test choose(t, 3.7)[1] == 2
    @test choose(t, 5.5)[1] == 3
    push!(t, ((2, 2.5), (3, 3.5)))
    v = [(2, 1),(3.5, 2),(4.8, 2),(6.1, 3),(7.9, 3)]
    for (guess, result) in v
        @test choose(t, guess)[1] == result
    end
end

@safetestset prefixsearch_four_floats = "four floats" begin
    using Fleck: PrefixSearchTree, choose, sum
    vals = [3, 1, 2, 4]
    t = PrefixSearchTree(vals)
    @test sum(t) == 10
    v = [(2, 1), (3, 2),(4.1, 3), (9.9, 4)]
    for (guess, result) in v
        @test choose(t, guess)[1] == result
    end
end

@safetestset prefixsearch_five = "five values" begin
    using Fleck: PrefixSearchTree, choose, sum
    vals = [3, 0, 2, 4, 1]
    t = PrefixSearchTree(vals)
    @test sum(t) == 10
    v = [(2, 1), (3, 3), (4.1, 3), (9.9, 5)]
    for (guess, result) in v
        @test choose(t, guess)[1] == result
    end
    push!(t, [(2, 7), (3, 3)])
    @test sum(t) == 18
    v = [(17.7, 5)]
    for (guess, result) in v
        @test choose(t, guess)[1] == result
    end
end


text:
function manhattan_distance(set1::Union{Vector, Array}, set2::Union{Vector, Array})

    @assert length(set1) == length(set2) "Sets don't have the same length."

    sum = zero(eltype(set1))

    for i in 1:length(set1)
        sum += abs(set2[i] - set1[i])
    end

    return sum
end


text:
# my_julia_main.jl
# include("my_code.jl")
using Sockets
using Dates
using Base.Filesystem
using Libdl

#hsa121345
# # Base.MainInclude.
# ghDLLRawProc = Libdl.dlopen("RawProc.dll")
Base.include(Main, "RawFileofflineProc/BaseLib.jl")

Base.MainInclude.include( "./RawFileofflineProc/ReadFiberData.jl")
Base.MainInclude.include( "./AlarmAnalyssys/AlarmPro/C01-RAW3AlgEntryRT.jl")
fileRawDataChannel = nothing

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    # @warn "hello"
   # return 0;
    global  fileRawDataChannel
    # ghDLLRawProc = Libdl.dlopen("RawProc.dll")
    Main.dllLoad( ) # baseLib.jl
    # Name = ghDLLRawProc
    # # baseGlbPath= replace(baseGlbPath,"\\"=>"/")
    # evs ="ghDLLRawProc = Libdl.dlopen(\"RawProc.dll\")"
    # reps = Base.Meta.parse( evs )
    # Core.eval( Main,reps )

    fileRawDataChannel = Channel(300)
    curpath = pwd()

    # println("::::--",curpath ," $ghDLLRawProc, $(Main.ghDLLRawProc)")
    # println("}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}")
    cfgfile = joinpath(curpath,"../Jl-Aux/glbCfg/baseCfg.jl")
    Main.getCfg( cfgfile )

    cfgObj = Main.cfgObj
    println()
    needFlowControl1 = 0
    try
    needFlowControl1 = ENV["NeedFlowControl"]
    catch
    end
    if needFlowControl1 == 0
        needFlowControl = false
    else
        needFlowControl = true
    end

    @async ClearExpiredFiles(cfgObj.RAW1MaxSpace, cfgObj.RAW2MaxSpace, cfgObj.RAW3MaxSpace, cfgObj.RAW4MaxSpace,cfgObj.Unit,cfgObj.folderA , cfgObj.folderB)
    while true
        try
            println("++++++++++++Init begin ")
            Main.Init(cfgObj.needSaveFile,cfgObj.needJuliaSend ,cfgObj.Logfolder,  cfgObj.folderA, cfgObj.folderB,needFlowControl)
            println("++++++++++++Init Ok")
            Main.InitIPAddressSvr(cfgObj.IpStr,cfgObj.Port ) ;
            println("++++++++++++InitIPAddressSvr Ok")
            Main.StartupMainThread()
            println("++++++++++++StartupMainThread Ok")
        catch(e)
            println("Exception: -------Init Svr Error :",e)
            break
        end

        try
            Main.SegParaAlarmProc.RAW3AlgEntryC01(true)
        catch(e)
            println("Exception: --------- Socket Receive data Error :",e)
            break
        end

        sleep(20)
    end
    return 0;
end
# exit()

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
using .PlotlyJS
import .PlotlyJS: plot

# NB: This file expects to be loaded *after* /src/compat/plotlyjs_wignerseitz.jl

const KPATH_COL = Ref("rgb(95,39,205)") # "nasu purple" (https://flatuicolors.com/palette/ca)

# ---------------------------------------------------------------------------------------- #

function plot(kp::KPath{3}, layout::Layout=DEFAULT_PLOTLY_LAYOUT_3D;
              config::PlotConfig = PlotConfig(responsive=true, displaylogo=false))
    setting(kp) !== CARTESIAN && (kp = cartesianize(kp))

    tpaths = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(paths(kp)))
    for (i,path) in enumerate(paths(kp))
        kvs = [points(kp)[lab] for lab in path]
        tpaths[i] = PlotlyJS.scatter3d(
            x=getindex.(kvs, 1), y=getindex.(kvs, 2), z=getindex.(kvs, 3),
            hoverinfo = "none",
            mode="lines", line=attr(color=KPATH_COL[], width=8))
    end
    tpoints = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(points(kp)))
    for (i,(lab, kv)) in enumerate(points(kp))
        tpoints[i] = PlotlyJS.scatter3d(
            x=kv[1:1], y=kv[2:2], z=kv[3:3],
            mode="marker", hovertext=string(lab)*" = 2π"*string(round.(kv/(2π), digits=3)),
            hoverinfo="text",
            marker=attr(color=KPATH_COL[], size=6, line=attr(color="white", width=1)))
    end
    return plot(vcat(tpoints, tpaths), layout; config=config)
end

# ---------------------------------------------------------------------------------------- #
function plot(kp::KPath{2}, layout::Layout=DEFAULT_PLOTLY_LAYOUT_2D;
              config::PlotConfig = PlotConfig(responsive=true, displaylogo=false))
    setting(kp) !== CARTESIAN && (kp = cartesianize(kp))

    tpaths = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(paths(kp)))
    for (i,path) in enumerate(paths(kp))
        kvs = [points(kp)[lab] for lab in path]
        tpaths[i] = PlotlyJS.scatter(
            x=getindex.(kvs, 1), y=getindex.(kvs, 2),
            hoverinfo = "none",
            mode="lines", line=attr(color=KPATH_COL[], width=5))
    end
    tpoints = Vector{GenericTrace{Dict{Symbol,Any}}}(undef, length(points(kp)))
    for (i,(lab, kv)) in enumerate(points(kp))
        tpoints[i] = PlotlyJS.scatter(
            x=kv[1:1], y=kv[2:2],
            mode="marker", hovertext=string(lab)*" = 2π"*string(round.(kv/(2π), digits=3)),
            hoverinfo="text",
            marker=attr(color=KPATH_COL[], size=8, line=attr(color="white", width=1)))
    end
    return plot(vcat(tpoints, tpaths), layout; config=config)
end
# ---------------------------------------------------------------------------------------- #

function plot(c::Cell{D}, kp::KPath{D},
              layout::Layout=(D==3 ? DEFAULT_PLOTLY_LAYOUT_3D : DEFAULT_PLOTLY_LAYOUT_2D);
              config::PlotConfig = PlotConfig(responsive=true, displaylogo=false)
              ) where D

    D ∉ (2,3) && error("must be 2D or 3D Cell and KPath")
    # ::Cell
    Pᶜ  = plot(c, layout)
    tsᶜ = Pᶜ.plot.data
    layout′ = Pᶜ.plot.layout # grab modified layout

    # get rid of "inside" axis lines; can be identified based on color
    filter!(tsᶜ) do t
        l = get(t, :line, nothing)
        l == nothing && return true
        c = get(l, :color, nothing)
        c == nothing && return true
        return !(c == BASIS_LIGHT_COL || c == AXIS_LIGHT_COL)
    end

    # ::KPath
    Pᵏᵖ  = plot(kp, layout′)
    tsᵏᵖ = Pᵏᵖ.plot.data

    # combine traces and plot
    ts = vcat(tsᶜ, tsᵏᵖ)
    return PlotlyJS.plot(ts, layout′; config=config)
end
# ---------------------------------------------------------------------------------------- #


text:
"""
   @lens struct Example
      field_name_a
   end

This macro will generate for each field of the struct, a `Lens` of the type
`Lens{Example,:field_name}` where
`getr(Lens{Example,:field_name}(), Example(value))` returns `value` and
`setr(Lens{Example,:field_name}(), Example(value), another_value)` returns `Example(another_value)`.
"""
macro lens(ex)
    lens(ex)
end

function lens(ex)
    if ex.head == :struct       # check to see if the expression is a struct
        if ex.args[1] == false      # check to see if it is immutable
            struct_name = get_name(ex.args[2])    # get the name of the new struct
            field_names = Symbol[]      # collect up the names of the fields
            for arg in ex.args[3].args
                if arg isa Expr
                    if arg.head == :(::)
                        push!(field_names, arg.args[1])
                    end
                end
            end
            #for each field name, make expression that build the getr and setter for the type
            final_quote = Expr(:block, esc(ex), (map(field_names) do name
                get = """begin function Lunettes.getr(::Lens{$struct_name, :$name}, a)
                        a.$name
                    end \n; """
                set = "function Lunettes.setr(::Lens{$struct_name, :$name}, a, c) \n " * " $struct_name( "
                for fname in field_names
                    if name == fname
                        set = set * " c, "
                    else
                        set = set * " a.$fname, "
                    end
                end
                set *=" );  end end"
                out = get * set
                out
                Meta.parse(out) |> esc
            end)..., nothing )
            final_quote
        else
            error("Expected a `struct` expression, but received a `mutable struct` instead.")
        end
    else
        error("Expected a `struct` expression.")
    end
end

function get_name(s::Symbol)
   s
end

function get_name(ex::Expr)
   ex.args[1]
end


text:
valid = 0

open("day02/day02.in") do io
    global valid
    for line in readlines(io)
        range_all, target_char, test_string = split(line, " ")
        target_char = target_char[1]
        a = parse(Int32, split(range_all, "-")[1])
        b = parse(Int32, split(range_all, "-")[2])

        if (target_char == test_string[a]) ⊻ (target_char == test_string[b])
            valid += 1
        end
    end
end

println(valid)

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
#!/usr/bin/env julia

# Code written entirely by Eric Alcaide.ç
# Testing the SearchMethods.jl package

# Start Test Script
using SearchMethods
using Base.Test

# Run tests
# An implementation of a Tree data structure for test purposes.
tree = Dict(
    1 => [2, 3, 4, 5],
	2 => [6, 7, 8],
	3 => [9, 6],
	4 => [11],
	5 => [12, 13],
	6 => [],
	7 => [14],
	8 => [],
	9 => [],
	10 => [],
	11 => [15],
	12 => [],
	13 => [16, 17, 18],
	14 => [19, 20],
	15 => [21, 22],
	16 => [23],
	17 => [],
	18 => [24, 25],
	19 => [],
	20 => [],
	21 => [],
	22 => [],
	23 => [26],
	24 => [],
	25 => [],
	26 => [27, 28, 29],
	27 => [],
	28 => [],
	29 => [],
 )

# Sample expand function
function expand_node(node, tree)
	return tree[node]
end

# Sample goal_check function - will return true
function goal_check_true(node)
	if node == 28 return true
	else return false
	end
end

# Sample goal_check function - will return false
function goal_check_false(node)
	if node == 30 return true
	else return false
	end
end

# Sample cost function
function cost(node, neighbor)
	return abs(node-neighbor)
end

# Sample heuristic function.
function heuristic(node, goal)
    return (abs(node-goal)/2)^2
end

tic()
# TESTS FOR UNINFORMED SEARCH ALGORITHMS
println("UNINFORMED METHODS")
# Sample test for the Breadth Frirst Search algorithm.
println("BFS Test")
@time @test if SearchMethods.bfs(1, tree, goal_check_true, expand_node) != false return true else false end
@time @test if SearchMethods.bfs(1, tree, goal_check_false, expand_node) == false return true else false end
# Sample test for the Depth First Search algorithm.
println("DFS Test")
@time @test if SearchMethods.dfs(1, tree, goal_check_true, expand_node) != false return true else false end
@time @test if SearchMethods.dfs(1, tree, goal_check_false, expand_node) == false return true else false end
# Sample test for the Depth Limited Search algorithm
println("DLS Test")
@time @test if SearchMethods.dls(1, tree, goal_check_true, expand_node, 9) != false return true else false end
@time @test if SearchMethods.dls(1, tree, goal_check_false, expand_node, 9) == false return true else false end
# Sample test for the IDS algorithm
# Sample test for the Iterative Deepening Search algorithm
println("IDS Test")
@time @test if SearchMethods.ids(1, tree, goal_check_true, expand_node, 2) != false return true else false end
@time @test if SearchMethods.ids(1, tree, goal_check_false, expand_node, 2) == false return true else false end
# Sample test for the Uniform Cost Search algorithm
println("UCS Test")
@time @test if SearchMethods.ucs(1, tree, goal_check_true, expand_node, cost) != false return true else false end
@time @test if SearchMethods.ucs(1, tree, goal_check_false, expand_node, cost) == false return true else false end
 # TESTS FOR INFORMED SEARCH ALGORITHMS
println("\nINFORMED METHODS")
# Sample test for the Greedy Best-First Search algorithm.
println("GBS Test")
@time @test if SearchMethods.greedy(1, tree, goal_check_true, expand_node, heuristic) != false return true else false end
@time @test if SearchMethods.greedy(1, tree, goal_check_false, expand_node, heuristic) == false return true else false end
# Sample test for the A* Search algorithm.
println("A*S Test")
@time @test if SearchMethods.a_star(1, tree, goal_check_true, expand_node, cost, heuristic) != false return true else false end
@time @test if SearchMethods.a_star(1, tree, goal_check_false, expand_node, cost, heuristic) == false return true else false end

toc()


text:
using Revise
using Plots
using DynamicalSystems
using DifferentialEquations
using Basins

### First version using high-level function.
b=0.1665
#b=0.22
ds = Systems.thomas_cyclical(b = b)
#integ=integrator(ds, reltol=1e-9)
xg=range(-6.,6.,length=200)
yg=range(-6.,6.,length=200)
pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9)

@time bsn = Basins.basins_map2D(xg, yg, pmap)

plot(xg,yg,bsn.basin',seriestype=:heatmap)


text:
import Base: setindex!

setindex!(x::TS, v, ::Colon) = x.values[:] .= v
setindex!(x::TS, v, ::Colon, ::Colon) = x[:] = v
setindex!(x::TS, v, ::Colon, c) = x.values[:, findcols(c,x)] .= v
setindex!(x::TS, v, r, ::Colon) = x.values[findrows(r,x), :] .= v
setindex!(x::TS, v, r, c) = x.values[findrows(r,x), findcols(c,x)] .= v

setindex!(x::TS, v, r::R) where {R<:TimeType} = x.values[findrows(r,x), :] .= v
setindex!(x::TS, v, r::R) where {R<:AbstractVector{<:TimeType}} = x.values[findrows(r,x), :] .= v
setindex!(x::TS, v, r::R) where {R<:AbstractString} = x.values[findrows(r,x), :] .= v
setindex!(x::TS, v, r::R) where {R<:Int} = x.values[findrows(r,x), :] .= v
setindex!(x::TS, v, r::R) where {R<:AbstractVector{<:Integer}} = x.values[findrows(r,x), :] .= v

setindex!(x::TS, v, c::C) where {C<:Symbol} = x.values[:, findcols(c,x)] .= v
setindex!(x::TS, v, c::C) where {C<:AbstractVector{Symbol}} = x.values[:, findcols(c,x)] .= v


text:
# ---
# title: 535. Encode and Decode TinyURL
# id: problem535
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Math
# link: <https://leetcode.com/problems/encode-and-decode-tinyurl/description/>
# hidden: true
# ---
#
# > Note: This is a companion problem to the [System
# > Design](https://leetcode.com/discuss/interview-question/system-design/)
# > problem: [Design TinyURL](https://leetcode.com/discuss/interview-
# > question/124658/Design-a-URL-Shortener-\(-TinyURL-\)-System/).
#
# TinyURL is a URL shortening service where you enter a URL such as
# `https://leetcode.comhttps://leetcode.com/problems/design-tinyurl` and it returns a short URL such
# as `http://tinyurl.com/4e9iAk`.
#
# Design the `encode` and `decode` methods for the TinyURL service. There is no
# restriction on how your encode/decode algorithm should work. You just need to
# ensure that a URL can be encoded to a tiny URL and the tiny URL can be decoded
# to the original URL.
#
#
## @lc code=start
using LeetCode

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


text:
###
### Series as linear functionals on polynomials.
###
### Hankel matrices associated to series.
###
### Bernard Mourrain
###

export series, hankel, hankelbasis

#-----------------------------------------------------------------------
"""
```
hankel(σ::Series{C,M}, L1::Vector{M}, L2::Vector{M}) -> Array{C,2}
```
Hankel matrix of ``σ`` with the rows indexed by the list of
polynomials L1 and the columns by L2. The entries are the dot product for ``σ``
of the corresponding elements in L1 and L2.

Example
-------
```jldoctest
julia> L =[1, x1, x2, x1^2, x1*x2, x2^2]

julia> H = hankel(s,L,L)
6x6 Array{Float64,2}:
  4.0   5.0   7.0    5.0  11.0  13.0
  5.0   5.0  11.0   -1.0  17.0  23.0
  7.0  11.0  13.0   17.0  23.0  25.0
  5.0  -1.0  17.0  -31.0  23.0  41.0
 11.0  17.0  23.0   23.0  41.0  47.0
 13.0  23.0  25.0   41.0  47.0  49.0
```
"""
function hankel(sigma::Series{C,M}, L1::AbstractVector, L2::AbstractVector) where {C,M}
   H = fill(zero(C), (length(L1), length(L2)));
   for i in 1:length(L1)
      for j in 1:length(L2)
         H[i,j] = dot(sigma, L1[i]*L2[j])
      end
   end
   H
end

#------------------------------------------------------------------------
"""
```
series(H::Matrix{C}, L1::Vector{M}, L2::Vector{M}) -> Series{C,M}
```
Compute the series associated to the Hankel matrix H, with rows (resp. columns)
indexed by the array of monomials L1 (resp. L2).
"""
function series(H::Array{C,2} , L1::Vector{M}, L2::Vector{M}) where {C, M}
    res = zero(Series{C,M})
    i = 1
    for m1  in L1
        j=1
        for m2 in L2
            m = m1*m2
            if res[m] == zero(C)
	        res[m] = H[i,j]
	    end
            j+=1
        end
        i+=1
    end
   res
end

#----------------------------------------------------------------------
"""
Generate the table of Hankel matrices ``H_α`` associated to the monomials ``x^α`` generating the space of Hankel matrices indexed by L1 (for the rows) and L2 (for the columns).
"""
function hankelbasis(L1,L2)
    m1 = length(L1)
    m2 = length(L2)
    table = Dict( L1[1]*L2[1] => fill(0.,m1,m2))
    for i in 1:m1
        for j in 1:m2
            M = get(table,L1[i]*L2[j], 0)
            if M != 0
                M[i,j]=1.0
            else
                M = fill(0.,m1, m2)
                M[i,j]=1.0
                table[L1[i]*L2[j]] = M
            end
        end
    end
    table
end


text:
using Pkg|
 Pkg.activate(".")

verbose = true
if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectML
include("shared/scripts.jl")
 if verbose   println("# Loading other packages")    end
 using DataFrames, Query, Statistics, Dates

include("read_neid_solar_data_20190918.jl")

if need_to(pipeline_plan,:extract_orders)
   if verbose println("# Extracting order list timeseries from spectra.")  end
   @assert !need_to(pipeline_plan,:read_spectra)
   order_list_timeseries = RvSpectML.make_order_list_timeseries(all_spectra)
   order_list_timeseries = RvSpectML.filter_bad_chunks(order_list_timeseries,verbose=true)
   #RvSpectML.normalize_spectra!(order_list_timeseries,all_spectra);
   dont_need_to!(pipeline_plan,:extract_orders);
 end
 order_list_timeseries

if need_to(pipeline_plan,:read_line_list)
   if verbose println("# Reading line list for CCF: ", linelist_for_ccf_filename, ".")  end
   lambda_range_with_good_data = get_λ_range(all_spectra)
   espresso_filename = joinpath(pkgdir(RvSpectML),"data","masks",linelist_for_ccf_filename)
   espresso_df = RvSpectML.read_linelist_espresso(espresso_filename)
   line_list_df = NEID.filter_line_list(espresso_df,first(all_spectra).inst)
   dont_need_to!(pipeline_plan,:read_line_list);
 end
 #line_list_df

line_list_no_tellurics_df = line_list_df
dont_need_to!(pipeline_plan,:clean_line_list_tellurics)
# TODO: Filter line list for tellurics when using non-EXPRES data
#=
#need_to!(pipeline_plan,:clean_line_list_tellurics)
if need_to(pipeline_plan,:clean_line_list_tellurics)
   if verbose println("# Removing lines with telluric contamination.")  end    # Currently only works for EXPRES data
   @assert !need_to(pipeline_plan,:read_line_list)
   @assert !need_to(pipeline_plan,:read_spectra)
   line_list_no_tellurics_df  = make_clean_line_list_from_tellurics_expres(line_list_df, all_spectra, Δv_to_avoid_tellurics = 30.0e3) #14000.0)
   dont_need_to!(pipeline_plan,:clean_line_list_tellurics);
 end
 #line_list_no_tellurics_df
=#

#need_to!(pipeline_plan,:ccf_total)
if need_to(pipeline_plan,:ccf_total)
   if verbose println("# Computing CCF.")  end
   @assert !need_to(pipeline_plan,:extract_orders)
   @assert !need_to(pipeline_plan,:clean_line_list_tellurics)
   mask_shape = RvSpectML.CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=tophap_ccf_mask_scale_factor)
   #line_list = RvSpectML.CCF.BasicLineList(line_list_df.lambda, line_list_df.weight)
   line_list = RvSpectML.CCF.BasicLineList(line_list_no_tellurics_df.lambda, line_list_no_tellurics_df.weight)
   ccf_plan = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=ccf_mid_velocity, range_no_mask_change=RvSpectMLBase.max_bc)
   v_grid = RvSpectML.CCF.calc_ccf_v_grid(ccf_plan)
   @time ccfs = RvSpectML.CCF.calc_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan)
   mask_shape_expr = RvSpectML.CCF.GaussianCCFMask(order_list_timeseries.inst, scale_factor=9)
   # Warning:  CCF with a shape other than a tophat is still experimental
   ccf_plan_expr = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape_expr, line_list=line_list, midpoint=ccf_mid_velocity)
   @time ccfs_expr = RvSpectML.CCF.calc_ccf_chunklist_timeseries_expr(order_list_timeseries, ccf_plan_expr) #, verbose=true)
   println("# Ratio of max(ccfs_expr)/max(ccfs) = ", mean(maximum(ccfs_expr,dims=1)./maximum(ccfs,dims=1)) )
   #=
   mask_shape_expr2 = RvSpectML.CCF.CosCCFMask(order_list_timeseries.inst, scale_factor=18)   # Why does such a large value help so much?
   ccf_plan_expr2 = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape_expr2, line_list=line_list, midpoint=ccf_mid_velocity)
   @time ccfs_expr2 = RvSpectML.CCF.calc_ccf_chunklist_timeseries_expr(order_list_timeseries, ccf_plan_expr2) #, verbose=true)
   =#
   if save_data(pipeline_plan, :ccf_total)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_ccfs.csv"),Tables.table(ccfs',header=Symbol.(v_grid)))
      CSV.write(joinpath(output_dir,target_subdir * "_ccfs_expr.csv"),Tables.table(ccfs_expr',header=Symbol.(v_grid)))
   end
   dont_need_to!(pipeline_plan,:ccf_total)
 end

if make_plot(pipeline_plan, :ccf_total)
   using Plots
   t_idx = 20
   plt = plot(v_grid,ccfs[:,t_idx]./maximum(ccfs[:,t_idx],dims=1),label=:none)
   scatter!(plt,v_grid,ccfs_expr[:,t_idx]./maximum(ccfs_expr[:,t_idx],dims=1),markersize=1.2,label=:none)
   xlabel!("v (m/s)")
   ylabel!("CCF")
   if save_plot(pipeline_plan,:ccf_total)   savefig(plt,joinpath(output_dir, target_subdir * "_ccf_sum.png"))   end
   display(plt)
end

if make_plot(pipeline_plan, :ccf_total)
   using Plots
   zvals = ccfs./maximum(ccfs,dims=1).-mean(ccfs./maximum(ccfs,dims=1),dims=2)
   colorscale = cgrad(:balance)
   plt = heatmap(v_grid,collect(1:size(ccfs,2)),zvals', c=colorscale, clims=(-maximum(abs.(zvals)),maximum(abs.(zvals))) )
   xlabel!("v (m/s)")
   ylabel!("Observation #")
   title!("CCF(v,t)-<CCF>(v) vs time")
   if save_plot(pipeline_plan,:ccf_total)   savefig(plt,joinpath(output_dir, target_subdir * "_ccf_sum_vs_time_heatmap.png"))   end
   display(plt)
end

#need_to!(pipeline_plan, :rvs_ccf_total)
if need_to(pipeline_plan, :rvs_ccf_total)
   if verbose println("# Measuring RVs from CCF.")  end
   @assert !need_to(pipeline_plan,:ccf_total)
   fit_gaussian_to_ccf = RVFromCCF.MeasureRvFromCCFGaussian()
   #fit_quadratic_to_ccf = RVFromCCF.MeasureRvFromCCFQuadratic()
   rvs_ccf_gauss = RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs,alg=fit_gaussian_to_ccf)
   rvs_ccf_gauss2 = RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs_expr,alg=fit_gaussian_to_ccf)
   #rvs_ccf_gauss3 = RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs_expr2,alg=fit_gaussian_to_ccf)
   println("# RMS of RVs:          Tophat ", std(rvs_ccf_gauss), "   Gaussian ", std(rvs_ccf_gauss))
   rms_rv_binned = bin_rvs_consecutive(rvs_ccf_gauss,4)
   rms_rv_expr_binned = bin_rvs_consecutive(rvs_ccf_gauss2,4)
   println("# RMS of binned RVs:  Tophat ", std(rms_rv_binned), "   Gaussian ", std(rms_rv_expr_binned))
   #rms_rv_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf_gauss)
   #rms_rv_expr_within_night = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf_gauss2)
   #println("# RMS within night RVs: Tophat ",rms_rv_within_night,  "   Gaussian ",rms_rv_expr_within_night)
   # Store estimated RVs in metadata
   map(i->order_list_timeseries.metadata[i][:rv_est] = rvs_ccf_gauss[i]-mean(rvs_ccf_gauss), 1:length(order_list_timeseries) )
   if save_data(pipeline_plan, :rvs_ccf_total)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_rvs_ccf.csv"),DataFrame("Time [MJD]"=>order_list_timeseries.times,"CCF RV [m/s]"=>rvs_ccf_gauss))
   end
   dont_need_to!(pipeline_plan, :rvs_ccf_total)
end

if make_plot(pipeline_plan, :rvs_ccf_total)
   using Plots
   rvs_ccf_gauss .-= mean(rvs_ccf_gauss)
   rvs_ccf_gauss2 .-= mean(rvs_ccf_gauss2)
   plt = scatter(rvs_ccf_gauss,markersize=3,label="RVs CCF Tophat", legend=:bottomright)
   scatter!(plt,rvs_ccf_gauss2,markersize=3,label="RVs CCF Gaussian")
   ylabel!("v (m/s)")
   xlabel!("Time (#)")
   if save_plot(pipeline_plan,:rvs_ccf_total)   savefig(plt,joinpath(output_dir, target_subdir * "_rvs_ccf_sum.png"))   end
   display(plt)

   #=
   using CSV
   df_yale_resutls = CSV.read(joinpath(homedir(),"Data/EXPRES/inputs/101501/101501_activity.csv"),DataFrame)
   rvs_yale_ccf = df_yale_resutls["CCF RV [m/s]"]
   rvs_yale_ccf .-= mean(rvs_yale_ccf)
   rvs_yale_cbc = df_yale_resutls["CBC RV [m/s]"]
   rvs_yale_cbc .-= mean(rvs_yale_cbc)
   scatter!(rvs_yale_cbc,markersize=2,label="Yale CBC")
   scatter!(rvs_yale_ccf,markersize=2,label="Yale CCF")
   =#
   #diff = rvs_yale_ccf.-rvs_yale_cbc
   #diff = rvs_ccf_gauss.-rvs_yale_cbc
   #diff = rvs_ccf_gauss2.-rvs_yale_cbc
   #diff = rvs_ccf_gauss.-rvs_yale_ccf
   #diff = rvs_ccf_gauss2.-rvs_yale_ccf
end

if make_plot(pipeline_plan, :rvs_ccf_total)
   diff = rvs_ccf_gauss.-rvs_ccf_gauss2
   println(std(diff))
   plt = scatter(order_list_timeseries.times,diff,markersize=4,label="Delta RV")
   ylabel!("Δv (m/s) (Two mask shapes)")
   xlabel!("Time (d)")
   if save_plot(pipeline_plan,:rvs_ccf_total)   savefig(plt,joinpath(output_dir, target_subdir * "_rvs_ccf_sum.png"))   end
   display(plt)
end

#need_to!(pipeline_plan, :ccf_orders)
 if need_to(pipeline_plan, :ccf_orders)  # Compute order CCF's & measure RVs
   tstart = now()    # Compute CCFs for each order
   order_ccfs = RvSpectML.CCF.calc_order_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan)
   println("# Order CCFs runtime: ", now()-tstart)

   if save_data(pipeline_plan, :ccf_orders)
      using CSV
      inst = EXPRES.EXPRES2D()
      for (i, order) in orders_to_use
         if !(sum(order_ccfs[:,i,:]) > 0)   continue    end
         local t = Tables.table( order_ccfs[:,i,:]', header=Symbol.(v_grid) )
         CSV.write(joinpath(output_dir,target_subdir * "_ccf_order=" * string(order) * ".csv"),t)
      end
   end
   dont_need_to!(pipeline_plan, :ccf_orders);
end

if make_plot(pipeline_plan, :ccf_orders)
   # order ccfs averaged over time
   obs = 1:size(order_ccfs,3)
   #plot(v_grid,reshape(sum(order_ccfs[:,ord,:],dims=3)./maximum(sum(order_ccfs[:,ord,:],dims=3),dims=1), size(order_ccfs,1),size(order_ccfs[:,ord,:],2)), label=:none)
   # resudiuals of order ccfs to time and chunk averaged ccf
   #plt = plot(v_grid,reshape(sum(order_ccfs[:,ord,:],dims=3)./maximum(sum(order_ccfs[:,ord,:],dims=3),dims=1),size(order_ccfs,1),size(order_ccfs[:,ord,:],2)).-
   #            reshape(sum(order_ccfs[:,ord,:],dims=(2,3))./maximum(sum(order_ccfs[:,ord,:],dims=(2,3))),size(order_ccfs,1) ), label=:none)
   order_labels = map(c->order_list_timeseries.chunk_list[1].data[c].λ.indices[2], 1:size(order_ccfs,2))
   zvals = reshape(sum(order_ccfs[:,:,obs],dims=3)./maximum(sum(order_ccfs[:,:,obs],dims=3),dims=1),size(order_ccfs,1),size(order_ccfs[:,:,obs],2)).-
   reshape(sum(order_ccfs[:,:,obs],dims=(2,3))./maximum(sum(order_ccfs[:,:,obs],dims=(2,3))),size(order_ccfs,1) )
   plt = heatmap(v_grid,order_labels, zvals',c = cgrad(:balance), clims=(-maximum(abs.(zvals)),maximum(abs.(zvals))) )

   xlabel!("v (m/s)")
   ylabel!("Order ID")
   title!("CCF-<CCF> for obs ID=" * string(obs))
   if save_plot(pipeline_plan,:ccf_orders)   savefig(plt,joinpath(output_dir, target_subdir * "_ccf_orders.png"))   end
   display(plt)
end

if make_plot(pipeline_plan, :ccf_orders)
   ord = 3 # TODO:  I think this order illustrates the issue Alex is going to fix with lines dropping off at the boundary of orders.  Check his fix solves the issue.
   #plot(v_grid,reshape(sum(order_ccfs[:,ord,:],dims=3)./maximum(sum(order_ccfs[:,ord,:],dims=3),dims=1), size(order_ccfs,1),size(order_ccfs[:,ord,:],2)), label=:none)
   # resudiuals of order ccfs to time and chunk averaged ccf
   #plt = plot(v_grid,reshape(sum(order_ccfs[:,ord,:],dims=3)./maximum(sum(order_ccfs[:,ord,:],dims=3),dims=1),size(order_ccfs,1),size(order_ccfs[:,ord,:],2)).-
   #            reshape(sum(order_ccfs[:,ord,:],dims=(2,3))./maximum(sum(order_ccfs[:,ord,:],dims=(2,3))),size(order_ccfs,1) ), label=:none)
   zvals = reshape(sum(order_ccfs[:,ord,:],dims=3)./maximum(sum(order_ccfs[:,ord,:],dims=3),dims=1),size(order_ccfs,1),size(order_ccfs[:,ord,:],2)).-
   reshape(sum(order_ccfs[:,ord,:],dims=(2,3))./maximum(sum(order_ccfs[:,ord,:],dims=(2,3))),size(order_ccfs,1) )
   plt = heatmap(v_grid,1:size(order_ccfs,3), zvals',c = cgrad(:balance), clims=(-maximum(abs.(zvals)),maximum(abs.(zvals))) )
   xlabel!("v (m/s)")
   ylabel!("Observation ID")
   title!("CCF-<CCF> for order=" * string(ord))
   if save_plot(pipeline_plan,:ccf_orders)   savefig(plt,joinpath(output_dir, target_subdir * "_ccf_orders.png"))   end
   display(plt)
end

RvSpectML.Pipeline.make_plot!(pipeline_plan,:scalpels)
if make_plot(pipeline_plan, :scalpels)
   @assert !need_to(pipeline_plan, :rvs_ccf_total)
   @assert !need_to(pipeline_plan, :ccf_total)
   include("../scripts/plots/scalpels_code.jl")
   plt = Scalpels.make_plots_scalpels(rvs_ccf_gauss, ccfs, max_num_basis=4, v_grid=v_grid, times=order_list_timeseries.times, output_path="examples/output/figures")
   display(plt)
end

#need_to!(pipeline_plan, :template)
if need_to(pipeline_plan, :template)  # Compute order CCF's & measure RVs
   if verbose println("# Making template spectra.")  end
   @assert !need_to(pipeline_plan,:extract_orders)
   @assert !need_to(pipeline_plan,:rvs_ccf_total)
   @time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = RvSpectML.make_template_spectra(order_list_timeseries)
   if save_data(pipeline_plan, :template)
      #using CSV
      #CSV.write(joinpath(output_dir,target_subdir * "_template.csv"),DataFrame("λ"=>spectral_orders_matrix.λ,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
      using JLD2, FileIO
      save(joinpath(output_dir,target_subdir * "_matrix.jld2"), Dict("λ"=>spectral_orders_matrix.λ,"spectra"=>spectral_orders_matrix.flux,"var_spectra"=>spectral_orders_matrix.var,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
   end
   dont_need_to!(pipeline_plan, :template);
end

if make_plot(pipeline_plan,:template)
   chunkid = 10
   idx = spectral_orders_matrix.chunk_map[chunkid]
   plt = scatter(spectral_orders_matrix.λ[idx],(f_mean[idx].-1.0)./maximum(abs.((f_mean[idx].-1.0))),markersize=1.0,label=:none)
   scatter!(plt,spectral_orders_matrix.λ[idx],deriv[idx]./maximum(abs.(deriv[idx])),markersize=1.0,label=:none)
   scatter!(plt,spectral_orders_matrix.λ[idx],deriv2[idx]./maximum(abs.(deriv2[idx])),markersize=1.0,label=:none)
   xlabel!("λ(Å)")
   ylabel!("f(λ), f'(λ), f''(λ), all standardized")
   title!("Template spectrum for chunk " * string(chunkid) )
end

#need_to!(pipeline_plan, :dcpca)
if need_to(pipeline_plan, :dcpca)
   if verbose println("# Performing Doppler constrained PCA analysis.")  end
   @assert !need_to(pipeline_plan,:rvs_ccf_total)
   @assert !need_to(pipeline_plan,:template)
   using MultivariateStats
   Δrvs_dcpca = calc_rvs_from_taylor_expansion(spectral_orders_matrix; mean=f_mean, deriv=deriv).rv
   dcpca_out, M = RvSpectML.DCPCA.doppler_constrained_pca(spectral_orders_matrix.flux, deriv, Δrvs_dcpca)
   frac_var_unexplained = 1.0.-cumsum(principalvars(M))./tvar(M)
   num_basis = length(frac_var_unexplained)
   println("# Fraction of variance unexplained = ", frac_var_unexplained[1:min(5,num_basis)])
   if save_data(pipeline_plan, :dcpca)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_dcpca_basis.csv"),  Tables.table(M.proj) )
      CSV.write(joinpath(output_dir,target_subdir * "_dcpca_scores.csv"), Tables.table(dcpca_out) )
   end
   dont_need_to!(pipeline_plan, :dcpca);
end

if make_plot(pipeline_plan, :dcpca)  # Ploting results from DCPCA
  include("../scripts/plots/dcpca.jl")
  # Set parameters for plotting analysis
  plt = scatter(frac_var_unexplained, xlabel="Number of PCs", ylabel="Frac Variance Unexplained")
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir, target_subdir * "_dcpca_frac_var.png"))   end
  display(plt)
end

if make_plot(pipeline_plan, :dcpca)
  plt_order = 10
  plt = RvSpectML.plot_basis_vectors(order_grids, f_mean, deriv, M.proj, idx_plt = spectral_orders_matrix.chunk_map[plt_order], num_basis=2 )
  #xlims!(5761.5,5766)
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir, target_subdir * "_dcpca_basis.png"))   end
  display(plt)
end

if make_plot(pipeline_plan, :dcpca)
  plt = RvSpectML.plot_basis_scores(order_list_timeseries.times, rvs_ccf_gauss, dcpca_out, num_basis=min(5,num_basis) )
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir, target_subdir * "_dcpca_scores.png"))   end
  display(plt)
end

if make_plot(pipeline_plan, :dcpca)
  plt = RvSpectML.plot_basis_scores_cor( rvs_ccf_gauss, dcpca_out, num_basis=min(5,num_basis))
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir, target_subdir * "_dcpca_scores_cor.png"))   end
  display(plt)
end

#need_to!(pipeline_plan,:fit_lines)
if need_to(pipeline_plan,:fit_lines)
   if verbose println("# Performing fresh search for lines in template spectra.")  end
   cl = ChunkList(map(grid->ChunkOfSpectrum(spectral_orders_matrix.λ,f_mean, var_mean, grid), spectral_orders_matrix.chunk_map))
   #= # We're done with the spectral_orders_matrix, so we can release the memory now
   spectral_orders_matrix = nothing
   GC.gc()
   need_to!(pipeline_plan,:template)
   =#
   lines_in_template = RvSpectML.LineFinder.find_lines_in_chunklist(cl, plan=RvSpectML.LineFinder.LineFinderPlan(min_deriv2=3))  # TODO: Automate threshold for finding a line

   if verbose println("# Finding above lines in all spectra.")  end
   @time fits_to_lines = RvSpectML.LineFinder.fit_all_lines_in_chunklist_timeseries(order_list_timeseries, lines_in_template )

   if save_data(pipeline_plan,:fit_lines)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_linefinder_lines.csv"), lines_in_template )
      CSV.write(joinpath(output_dir,target_subdir * "_linefinder_line_fits.csv"), fits_to_lines )
      #CSV.write(joinpath(output_dir,target_subdir * "_linefinder_line_fits_clean.csv"), lines_to_try )
   end
   dont_need_to!(pipeline_plan,:fit_lines);
end

#need_to!(pipeline_plan,:clean_line_list_blends)
if need_to(pipeline_plan,:clean_line_list_blends)
   # Exploratory data analysis of the distribution of line properties over time to figure out how to select "clean" lines
   @assert !need_to(pipeline_plan,:fit_lines)
   fit_distrib = fits_to_lines |> @groupby(_.line_id) |>
            @map( { median_a=median(_.fit_a), median_b=median(_.fit_b), median_depth=median(_.fit_depth), median_σ²=median(_.fit_σ²), median_λc=median(_.fit_λc),
                   std_a=std(_.fit_a), std_b=std(_.fit_b), std_depth=std(_.fit_depth), std_σ²=std(_.fit_σ²), std_λc=std(_.fit_λc), line_id=_.line_id, frac_converged=mean(_.fit_converged)  } ) |>
            @filter(_.frac_converged == 1.0 ) |> DataFrame
                  # , min_telluric_model_all_obs=minimum(_.min_telluric_model_this_obs),
   good_lines = fit_distrib |>
            #@filter(_.min_telluric_model_all_obs == 1.0 ) |> # Already done
            @filter( _.median_depth > 0.1 ) |>
            @filter( _.median_σ² < 0.015) |>
            @filter( _.std_σ² < 0.004) |>
            @filter(_.std_b < 0.06) |>
            @filter( _.std_depth < 0.006 ) |>
            #@filter(_.std_a < 0.025) |>             # Unncessary?
            #@filter( -0.25 < _.median_b < 0.25) |>  # Unncessary?
            DataFrame
      println("# Found ", size(lines_in_template,1), " lines, including ",size(good_lines,1), " good lines, rejected ", size(lines_in_template,1)-size(good_lines,1), " lines.")
      bad_lines = fit_distrib |> @filter(_.std_λc >= 0.002 ) |> DataFrame
      good_lines_high_scatter = good_lines |> @filter(_.std_λc >= 0.002 ) |> DataFrame
      if size(bad_lines,1) >= 1
         println("# ", size(bad_lines,1), " lines have large λc scatter, including ", size(good_lines_high_scatter,1), " good lines.")
      end
      lines_to_try = lines_in_template[first.(good_lines[!,:line_id]),:]
      dont_need_to!(pipeline_plan,:clean_line_list_blends);
 end
 fit_distrib

#=if make_plot(pipeline_plan,:fit_lines) || true # Need to add key to pipeline_plan
scatter(good_lines.median_depth,good_lines.std_λc,markersize=2,label=:none)
scatter(good_lines.std_depth,good_lines.std_λc,markersize=2,label=:none)
scatter(good_lines.std_σ²,good_lines.std_λc,markersize=2,label=:none)
scatter(good_lines.median_σ²,good_lines.std_λc,markersize=2,label=:none)

scatter(good_lines.median_b,good_lines.std_λc,markersize=2,label=:none)
scatter(good_lines.std_b,good_lines.std_λc,markersize=2,label=:none)

scatter(good_lines.median_b,good_lines.std_λc,markersize=2,label=:none)
scatter(good_lines.std_b,good_lines.std_λc,markersize=2,label=:none)
end
=#

# Reset steps of pipeline_plan to rerun with new linelist.
need_to!(pipeline_plan, :ccf_total);
 need_to!(pipeline_plan, :rvs_ccf_total);
 need_to!(pipeline_plan, :ccf_orders);
 need_to!(pipeline_plan, :rvs_ccf_orders);
 need_to!(pipeline_plan, :template);
 need_to!(pipeline_plan, :dcpca);

#need_to!(pipeline_plan,:ccf_total)
if need_to(pipeline_plan,:ccf_total)
   if verbose println("# Computing CCFs with new line list.")  end
   @assert !need_to(pipeline_plan,:extract_orders)
   @assert !need_to(pipeline_plan,:clean_line_list_blends)
   #mask_shape = RvSpectML.CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=tophap_ccf_mask_scale_factor)
   perm = sortperm(lines_to_try.fit_λc)
   line_list2 = RvSpectML.CCF.BasicLineList(lines_to_try.fit_λc[perm], lines_to_try.fit_depth[perm] )
   #mask_shape = RvSpectML.CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=tophap_ccf_mask_scale_factor*3)
   mask_shape = RvSpectML.CCF.GaussianCCFMask(order_list_timeseries.inst, scale_factor=4)
   ccf_plan2 = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape, step=100.0, line_list=line_list2, midpoint=0.0)
   v_grid2 = RvSpectML.CCF.calc_ccf_v_grid(ccf_plan2)
   @time ccfs2 = RvSpectML.CCF.calc_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan2)
   # Write CCFs to file
   if save_data(pipeline_plan,:ccf_total)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_ccfs2.csv"),Tables.table(ccfs2',header=Symbol.(v_grid)))
   end
   dont_need_to!(pipeline_plan,:ccf_total)
end

if make_plot(pipeline_plan, :ccf_total)
   using Plots
   t_idx = 20
   using Plots
   plt = plot(v_grid,ccfs[:,t_idx]./maximum(ccfs[:,t_idx],dims=1),label=:none)
   scatter!(plt,v_grid2,ccfs2[:,t_idx]./maximum(ccfs2[:,t_idx],dims=1),markersize=1.2,label=:none)
   xlabel!("v (m/s)")
   ylabel!("CCF")
   if save_plot(pipeline_plan,:ccf_total)   savefig(plt,joinpath(output_dir, target_subdir * "_ccf2_sum.png"))   end
   display(plt)
end

#need_to!(pipeline_plan,:rvs_ccf_total)
if need_to(pipeline_plan, :rvs_ccf_total)
   if verbose println("# Measuring RVs from CCF.")  end
   @assert !need_to(pipeline_plan,:ccf_total)
   rvs_ccf_gauss2 = RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid2,ccfs2,alg=fit_gaussian_to_ccf)
   # Store estimated RVs in metadata
   println("# RMS of RVs:           Orig lines ", std(rvs_ccf_gauss),  "    Cleaned lines ", std(rvs_ccf_gauss2))
   rms_rv_binned = bin_rvs_consecutive(rvs_ccf_gauss,4)
   rms_rv_binned2 = bin_rvs_consecutive(rvs_ccf_gauss2,4)
   println("# RMS of binned RVs:  Orig lines ", std(rms_rv_binned2), "   Cleanded lines ", std(rms_rv_expr_binned2))
   #rms_rv_nightly2 = bin_rvs_nightly(times=order_list_timeseries.times,rvs=rvs_ccf_gauss2)
   #println("# RMS of nightly RVs:   Orig lines ", std(rms_rv_nightly), "    Cleaned lines ", std(rms_rv_nightly2))
   #rms_rv_within_night2 = rms_rvs_within_night(times=order_list_timeseries.times,rvs=rvs_ccf_gauss2)
   #println("# RMS within night RVs: Orig lines ",rms_rv_within_night,  "   Cleaned liens ",rms_rv_within_night2)
   map(i->order_list_timeseries.metadata[i][:rv_est] = rvs_ccf_gauss2[i]-mean(rvs_ccf_gauss2), 1:length(order_list_timeseries) )
   if save_data(pipeline_plan, :rvs_ccf_total)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_rvs_ccf2.csv"),DataFrame("Time [MJD]"=>order_list_timeseries.times,"CCF RV [m/s]"=>rvs_ccf_gauss))
   end
   dont_need_to!(pipeline_plan, :rvs_ccf_total);
end

#need_to!(pipeline_plan, :template)
if need_to(pipeline_plan, :template)
   @assert !need_to(pipeline_plan,:rvs_ccf_total)
   @assert !need_to(pipeline_plan,:clean_line_list_blends)
   # Revert to velocities before cleaning for now
   map(i->order_list_timeseries.metadata[i][:rv_est] = rvs_ccf_gauss[i]-mean(rvs_ccf_gauss), 1:length(order_list_timeseries) )
   #chunk_list_df2 = lines_to_try |> @select(:fit_min_λ,:fit_max_λ) |> @rename(:fit_min_λ=>:lambda_lo, :fit_max_λ=>:lambda_hi) |> DataFrame
   expand_chunk_factor = 4
   chunk_list_df2 = lines_to_try |> @select(:fit_λc,:fit_min_λ,:fit_max_λ) |> @map({:lambda_lo=>_.fit_λc-expand_chunk_factor*(_.fit_λc-_.fit_min_λ), :lambda_hi=>_.fit_λc+expand_chunk_factor*(_.fit_max_λ-_.fit_λc)}) |> DataFrame
   chunk_list_timeseries2 = RvSpectML.make_chunk_list_timeseries(all_spectra,chunk_list_df2)

   # Check that no NaN's included
   #(chunk_list_timeseries2, chunk_list_df2) = RvSpectML.filter_bad_chunks(chunk_list_timeseries2,chunk_list_df2)
   chunk_list_timeseries2 = RvSpectML.filter_bad_chunks(chunk_list_timeseries2)
   println(size(chunk_list_df2), " vs ", num_chunks(chunk_list_timeseries2) )

   if verbose println("# Making template spectra.")  end
   @time ( spectral_orders_matrix2, f_mean2, var_mean2, deriv_2, deriv2_2, order_grids2 )  = RvSpectML.make_template_spectra(chunk_list_timeseries2)

   if save_data(pipeline_plan, :template)
      #using CSV
      #CSV.write(joinpath(output_dir,target_subdir * "_template2.csv"),DataFrame("λ"=>spectral_orders_matrix2.λ,"flux_template"=>f_mean2,"var"=>var_mean2, "dfluxdlnλ_template"=>deriv_2,"d²fluxdlnλ²_template"=>deriv2_2))
      using JLD2, FileIO
      save(joinpath(output_dir,target_subdir * "_matrix2.jld2"), Dict("λ"=>spectral_orders_matrix2.λ,"spectra"=>spectral_orders_matrix2.flux,"var_spectra"=>spectral_orders_matrix2.var,"flux_template"=>f_mean2,"var"=>var_mean2, "dfluxdlnλ_template"=>deriv_2,"d²fluxdlnλ²_template"=>deriv2_2))
   end
   dont_need_to!(pipeline_plan, :template)
end

#rvs_dcpca = calc_rvs_from_taylor_expansion(spectral_orders_matrix; mean=f_mean, deriv=deriv)
#std(rvs_dcpca.rv)

#std(rvs_dcpca2.rv)

#need_to!(pipeline_plan, :dcpca)
if need_to(pipeline_plan, :dcpca)
   @assert !need_to(pipeline_plan,:rvs_ccf_total)
   @assert !need_to(pipeline_plan,:template)
   if verbose println("# Performing Doppler constrained PCA analysis.")  end
   using MultivariateStats
   Δrvs_dcpca2 = calc_rvs_from_taylor_expansion(spectral_orders_matrix2; mean=f_mean2, deriv=deriv_2).rv
   dcpca2_out, M2 = RvSpectML.DCPCA.doppler_constrained_pca(spectral_orders_matrix2.flux, deriv_2, Δrvs_dcpca2)
   frac_var_unexplained2 = 1.0.-cumsum(principalvars(M2))./tvar(M2)
   num_basis = length(frac_var_unexplained2)
   println("# Fraction of variance unexplained (orig) = ", frac_var_unexplained[1:min(5,num_basis)])
   println("# Fraction of variance unexplained (clean) = ", frac_var_unexplained2[1:min(5,num_basis)])
   if save_data(pipeline_plan, :dcpca)
      using CSV
      CSV.write(joinpath(output_dir,target_subdir * "_dcpca_basis2.csv"),  Tables.table(M2.proj) )
      CSV.write(joinpath(output_dir,target_subdir * "_dcpca_scores2.csv"), Tables.table(dcpca2_out) )
   end
   dont_need_to!(pipeline_plan, :dcpca)
end

if make_plot(pipeline_plan, :dcpca)  # Ploting results from DCPCA
  using Plots
  plt = scatter(frac_var_unexplained2, xlabel="Number of PCs", ylabel="Frac Variance Unexplained")
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir, target_subdir * "_dcpca2_frac_var.png"))   end
  display(plt)
end

if make_plot(pipeline_plan, :dcpca)
   plt_line = 10
   plt = RvSpectML.plot_basis_vectors(order_grids2, f_mean2, deriv_2, M2.proj, idx_plt = spectral_orders_matrix2.chunk_map[plt_line], num_basis=min(4,num_basis), label=plt_line)
   if save_plot(pipeline_plan,:dcpca)
     savefig(plt,joinpath(output_dir, target_subdir * "_dcpca2_basis.png"))
   end
   suspicous_lines = Int[]
   anim = @animate for plt_line ∈ 1:length(order_grids2)
      local idx = spectral_orders_matrix2.chunk_map[plt_line]
      argmin_pixel = argmin(f_mean2[idx])
      if 0.45*length(idx) <= argmin_pixel <= 0.55*length(idx)
        RvSpectML.plot_basis_vectors(order_grids2, f_mean2, deriv_2, M2.proj, idx_plt = spectral_orders_matrix2.chunk_map[plt_line], num_basis=min(4,num_basis), label=plt_line  )
     else
        push!(suspicous_lines,plt_line)
     end
   end
   gif(anim, joinpath(output_dir, target_subdir * "_dcpca2_basis.gif"), fps = 25)
  display(plt)
end

if make_plot(pipeline_plan, :dcpca)
  plt = RvSpectML.plot_basis_scores(order_list_timeseries.times, rvs_ccf_gauss2, dcpca2_out, num_basis=min(4,num_basis) )
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir, target_subdir * "_dcpca2_scores.png"))   end
  display(plt)

end

if make_plot(pipeline_plan, :dcpca)
  plt = RvSpectML.plot_basis_scores_cor( rvs_ccf_gauss2, dcpca2_out, num_basis=min(4,num_basis))
  if save_plot(pipeline_plan,:dcpca)   savefig(plt,joinpath(output_dir, target_subdir * "_dcpca2_score_cor.png"))   end
  display(plt)
end


text:
## loading prover

include("load_cheaplogic.jl")

include("lps.jl")

println("lps is loaded")

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
# test Odo, Obsv2

using IncrementalInference
using Statistics
using Test

N=100
fg = initfg()

doors = zeros(1,1)
pd = kde!(doors,[3.0])
pd = resample(pd,N);
bws = getBW(pd)[:,1]
doors2 = getPoints(pd);

@testset "test evaluation of pose pose constraint" begin

v1 = addVariable!(fg,:x1, ContinuousScalar,N=N)
f1  = addFactor!(fg,[v1], Prior( pd )) #, samplefnc=getSample

# tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addVariable!(fg,:x2, ContinuousScalar, N=N)
odoc = LinearConditional(Normal(50.0,2.0)) # Odo(50.0*ones(1,1),2.0*ones(1,1),[1.0])
f2 = addFactor!(fg, [:x1; :x2], odoc ) #, samplefnc=getSample

# @test isInitialized(fg, :x1)

pts = approxConv(fg, :x1x2f1, :x2)
# pts = evalFactor2(fg, f2, v2.label)
@show Statistics.mean(pts,dims=2)
@test norm(Statistics.mean(pts,dims=2)-[50.0]) < 15.0

tree, smt, hist = solveTree!(fg)
# ensureAllInitialized!(fg)
# tree = wipeBuildNewTree!(fg, drawpdf=false)
# inferOverTree!(fg, tree)

@test norm(Statistics.mean(getVal(fg, :x2),dims=2)-[50.0]) < 15.0

end

# using RoMEPlotting
#
# plotKDE( getVertKDE(fg,:x2) )
# plotKDE(kde!(pts))

#


text:
# Module for tilecoding functionality
include("utilities/TileCoding.jl")

"""
    TileCoder(num_tilings, num_tiles, num_features, num_ints)
Tile coder for coding all features together.
"""
mutable struct TileCoder <: AbstractFeatureCreator
    # Main Arguments
    num_tilings::Int64
    num_tiles::Int64
    num_features::Int64
    num_ints::Int64

    # Optional Arguments
    wrap::Bool
    wrapwidths::Float64

    iht::TileCoding.IHT
    TileCoder(num_tilings, num_tiles, num_features, num_ints=1; wrap=false, wrapwidths=0.0) =
        new(num_tilings, num_tiles, num_features, num_ints, wrap, 0.0, TileCoding.IHT(num_tilings*(num_tiles+1)^num_features * num_ints))
end

function create_features(fc::TileCoder, s; ints=[], readonly=false)
    if fc.wrap
        return 1 .+ TileCoding.tileswrap!(fc.iht, fc.num_tilings, s.*fc.num_tiles, fc.wrapwidths, ints, readonly)
    else
        return 1 .+ TileCoding.tiles!(fc.iht, fc.num_tilings, s.*fc.num_tiles, ints, readonly)
    end
end

feature_size(fc::TileCoder) = fc.num_tilings

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
import .CuArrays: CuArray
import .CuArrays.CUDAdrv: CuPtr, synchronize
import .CuArrays.CUDAdrv.Mem: DeviceBuffer

function Base.cconvert(::Type{MPIPtr}, buf::CuArray{T}) where T
    Base.cconvert(CuPtr{T}, buf) # returns DeviceBuffer
end

# CuArrays <= v1.3
function Base.unsafe_convert(::Type{MPIPtr}, buf::DeviceBuffer)
    reinterpret(MPIPtr, buf.ptr)
end
# CuArrays > v1.3
function Base.unsafe_convert(::Type{MPIPtr}, X::CuArray{T}) where T
    reinterpret(MPIPtr, Base.unsafe_convert(CuPtr{T}, X))
end
# only need to define this for strided arrays: all others can be handled by generic machinery
function Base.unsafe_convert(::Type{MPIPtr}, V::SubArray{T,N,P,I,true}) where {T,N,P<:CuArray,I}
    X = parent(V)
    pX = Base.unsafe_convert(CuPtr{T}, X)
    pV = pX + ((V.offset1 + V.stride1) - first(LinearIndices(X)))*sizeof(T)
    return reinterpret(MPIPtr, pV)
end

function Buffer(arr::CuArray)
    Buffer(arr, Cint(length(arr)), Datatype(eltype(arr)))
end


text:
using DataFrames
using CSV
using Random
using Dates

mutable struct Trace
    seed::Any
    algorithm::AbstractAlgorithm
    individual::Instance
    rows::DataFrames.DataFrame
    optimum::NamedTuple{(:individual, :iteration),Tuple{Any,Integer}}

    function Trace(
        algorithm::AbstractAlgorithm,
        individual::Instance;
        seed::Integer = ceil(Integer, (time() * 10e6) % typemax(Int)),
        rows::DataFrames.DataFrame = DataFrames.DataFrame(),
        optimum = (individual = missing, iteration = -1),
    )
        new(seed, algorithm, individual, rows, optimum)
    end
end

function deepparameters(object)
    allfields = Dict()
    for fieldname in fieldnames(typeof(object))
        if typeof(getfield(object, fieldname)) <: Real
            push!(allfields, (fieldname => getfield(object, fieldname)))
        else
            subfield = getfield(object, fieldname)
            push!(allfields, (fieldname => string(getfield(object, fieldname))))
            for key_value in deepparameters(subfield)
                push!(
                    allfields,
                    (
                        Symbol(string(fieldname) * "_" * string(key_value[1])) =>
                            key_value[2]
                    ),
                )
            end
        end
    end
    allfields
end

function realparameters(algorithm::T) where {T<:AbstractAlgorithm}
    realfields = []
    for fieldname in fieldnames(typeof(algorithm))
        if typeof(getfield(algorithm, fieldname)) <: Real
            push!(realfields, fieldname)
        end
    end
    realfields
end

function storeresult(trace::Trace)
    date_format = Dates.format(now(), "yyyymmddHHMMSS_s")
    id = trace.seed

    result = Dict(
        "date" => date_format,
        "algorithm" => string(typeof(trace.algorithm)),
        "size" => length(trace.individual.problem),
        "iteration" => trace.optimum[2],
        "seed" => trace.seed,
    )

    for fieldname in realparameters(trace.algorithm)
        push!(result, (string(fieldname) => getfield(trace.algorithm, fieldname)))
    end

    result_df = DataFrame(result)

    csvfilename = string(typeof(trace.algorithm))
    writeheader = false
    rlock = Threads.Condition()
    lock(rlock)
    if !ispath("results")
        mkdir("results")
    end
    if !ispath("results/$csvfilename.csv")
        touch("results/$csvfilename.csv")
        writeheader = true
    end
    CSV.write("results/$csvfilename.csv", result_df, append = true, header = writeheader)
    unlock(rlock)
end

function info(trace::Trace, text, data)
    @debug text data...
end

function optimumfound(trace::Trace, x::Instance, iteration::Integer)
    storeresult(trace)
end

function record(db::Trace, x::Instance, iteration::Integer, isoptimum::Bool = false)
    tuple = (iteration = iteration, fitness = fitness(x))
    info(db, "New row is inserted.", tuple)
    push!(db.rows, (iteration = iteration, fitness = fitness(x)))
    db.optimum = (individual = x.individual, iteration = iteration)
    isoptimum && optimumfound(db, x, iteration)
end

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 837951f4-68ec-11eb-048d-6b27b9055428
using Pkg, DrWatson

# ╔═╡ 91489170-68ec-11eb-3611-150d5febf223
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanQuap
	using StatisticalRethinking
end

# ╔═╡ 7046ea12-68ec-11eb-18ae-a724ca79a21f
md" ## Clip-07-32-34s.jl"

# ╔═╡ d603e8be-8381-11eb-1c62-8d87639fc38b
versioninfo()

# ╔═╡ ace5ed4c-68ec-11eb-05ed-1f1ff2d72e1b
begin
	df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);
	scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])
	data = (N=size(df, 1), D=df.Divorce_s, A=df.MedianAgeMarriage_s,
		M=df.Marriage_s)
end;

# ╔═╡ 4c47f312-68ed-11eb-37de-ff782118d258
stan5_1 = "
data {
	int < lower = 1 > N; // Sample size
	vector[N] D; // Outcome
	vector[N] A; // Predictor
}
parameters {
	real a; // Intercept
	real bA; // Slope (regression coefficients)
	real < lower = 0 > sigma;    // Error SD
}
transformed parameters {
	vector[N] mu;               // mu is a vector
	for (i in 1:N)
		mu[i] = a + bA * A[i];
}
model {
	a ~ normal(0, 0.2);         //Priors
	bA ~ normal(0, 0.5);
	sigma ~ exponential(1);
	D ~ normal(mu , sigma);     // Likelihood
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

# ╔═╡ 4c482dfa-68ed-11eb-2db9-a54ab816beb7
begin
	m5_1s = SampleModel("m5.1s", stan5_1)
	rc5_1s = stan_sample(m5_1s; data)
	if success(rc5_1s)
		post5_1s_df = read_samples(m5_1s, :dataframe)
		PRECIS(post5_1s_df[:, [:a, :bA, :sigma]])
	end
end

# ╔═╡ 8f107504-68ed-11eb-002b-372884f598b0
stan5_2 = "
data {
	int N;
	vector[N] D;
	vector[N] M;
}
parameters {
	real a;
	real bM;
	real<lower=0> sigma;
}
transformed parameters {
	vector[N] mu;
	for (i in 1:N)
		mu[i]= a + bM * M[i];

}
model {
	a ~ normal( 0 , 0.2 );
	bM ~ normal( 0 , 0.5 );
	sigma ~ exponential( 1 );
	D ~ normal( mu , sigma );
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

# ╔═╡ cbb9ef56-68ed-11eb-387b-252ff6dc5b21
begin
	m5_2s = SampleModel("m5.2s", stan5_2);
	rc5_2s = stan_sample(m5_2s; data)
	if success(rc5_2s)
		post5_2s_df = read_samples(m5_2s, :dataframe)
		PRECIS(post5_2s_df[:, [:a, :bM, :sigma]])
	end
end

# ╔═╡ 71a48f7a-68ee-11eb-1b0e-331393e95551
stan5_3 = "
data {
  int N;
  vector[N] D;
  vector[N] M;
  vector[N] A;
}
parameters {
  real a;
  real bA;
  real bM;
  real<lower=0> sigma;
}
transformed parameters {
	vector[N] mu;
	mu = a + + bA * A + bM * M;
}
model {
  a ~ normal( 0 , 0.2 );
  bA ~ normal( 0 , 0.5 );
  bM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  D ~ normal( mu , sigma );
}
generated quantities{
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

# ╔═╡ b5b487ec-68ee-11eb-3ee7-0d98367028ff
begin
	init = (a = 0.0, bA = -1.0, bM = 0.0, sigma = 1.0)
	q5_3s, m5_3s, o5_3s = stan_quap("m5.3s", stan5_3; data, init)

	if !isnothing(m5_3s)
		post5_3s_df = read_samples(m5_3s, :dataframe)
		PRECIS(post5_3s_df[:, [:a, :bA, :bM, :sigma]])
	end
end

# ╔═╡ 7310fa28-7132-11eb-32dc-e9aae12b4f9e
md" ### Compare the model cofficients."

# ╔═╡ 84d0ef98-7132-11eb-1d64-3d911674a004
begin
	plot_models([m5_1s, m5_2s, m5_3s], [:a, :bA, :bM, :sigma])
end

# ╔═╡ 8501f844-69b3-11eb-24f8-d3b87b745be7
if !isnothing(m5_3s)
	waic(m5_3s)
end

# ╔═╡ c2ba84b0-68f2-11eb-3c5a-17e03541dcb1
begin
    b5_3s = post5_3s_df[:, [:a, :bA, :bM, :sigma]]
    mu5_3s = b5_3s.a .+ b5_3s.bM * df.Marriage_s' +
		b5_3s.bA * df.MedianAgeMarriage_s'
	log_lik5_3s = logpdf.(Normal.(mu5_3s, post5_3s_df.sigma),  df.Divorce_s')
	waic(log_lik5_3s)
end

# ╔═╡ 9ee34a3a-6d4e-11eb-1e6e-b191bea5527b
    df_waic = compare([m5_1s, m5_2s, m5_3s], :waic)

# ╔═╡ b05f0250-7009-11eb-0370-4969eb2fd5ce
plot_models([m5_1s, m5_2s, m5_3s], :waic)

# ╔═╡ 85a11cf4-6d40-11eb-3424-5b873aa5bc88
md"With quap():
```
      PSIS    SE   dPSIS   dSE   pPSIS   weight
 m5.1 127.6 14.69   0.0    NA     4.7    0.71
 m5.3 129.4 15.10   1.8   0.90    5.9    0.29
 m5.2 140.6 11.21  13.1  10.82    3.8    0.00
```

or, with ulam():
```
       PSIS    SE dPSIS  dSE pPSIS weight
m5.1u 126.0 12.83   0.0   NA   3.7   0.67
m5.3u 127.4 12.75   1.4 0.75   4.7   0.33
m5.2u 139.5  9.95  13.6 9.33   3.0   0.00
```
"

# ╔═╡ c6955180-6d40-11eb-0723-e55c69316f63
df_psis = compare([m5_1s, m5_2s, m5_3s], :psis)

# ╔═╡ 2c4b5ffc-68f4-11eb-21d6-7d2db834d3f8
begin
	loo5_1s, loos5_1s, pk5_1s = psisloo(m5_1s)
	pk_plot(pk5_1s)
	annotate!([(13 + 1, pk5_1s[13] + 0.02, Plots.text(df[13, :Loc],
		6, :red, :right))])
end

# ╔═╡ 3e956c3c-6ee6-11eb-2aff-d73d8c0d9dcd
waic(m5_1s)

# ╔═╡ 468ae700-6ee6-11eb-289d-8dabcdca1b7d
exp(std(loos5_1s))

# ╔═╡ aedb425c-68f4-11eb-319d-a958c2aacde1
begin
	loo5_3s, loos5_3s, pk5_3s = psisloo(m5_3s)
	pk_plot(pk5_3s)
	annotate!([(13 + 1, pk5_3s[13] + 0.02, Plots.text(df[13, :Loc],
		6, :red, :right))])
end

# ╔═╡ f3afc7c2-6d68-11eb-10a6-65ab80be96a6
waic_5_3s_pw = waic(m5_3s; pointwise=true)

# ╔═╡ bb70ce3a-6956-11eb-26ad-e19934980e11
waic_5_3s_pw.penalty[13]

# ╔═╡ eb7c7c14-6956-11eb-3c98-3b4ee223bccc
begin
	scatter(pk5_3s, waic_5_3s_pw.penalty,
		xlab="PSIS Pareto k", ylab="WAIC penalty", leg=false)
	vline!([0.5])
	annotate!([([pk5_3s[13]], [waic_5_3s_pw.penalty[13] + 0.02],
		Plots.text(df[13, :Loc], 6, :red, :right))])
	annotate!([([pk5_3s[20]], [waic_5_3s_pw.penalty[20] + 0.02],
		Plots.text(df[20, :Loc], 6, :red, :right))])

end

# ╔═╡ c3169ed8-694e-11eb-0616-97d25d8b56e1
md" ## End of clip-07-32-34s.jl"

# ╔═╡ Cell order:
# ╟─7046ea12-68ec-11eb-18ae-a724ca79a21f
# ╠═837951f4-68ec-11eb-048d-6b27b9055428
# ╠═91489170-68ec-11eb-3611-150d5febf223
# ╠═d603e8be-8381-11eb-1c62-8d87639fc38b
# ╠═ace5ed4c-68ec-11eb-05ed-1f1ff2d72e1b
# ╠═4c47f312-68ed-11eb-37de-ff782118d258
# ╠═4c482dfa-68ed-11eb-2db9-a54ab816beb7
# ╠═8f107504-68ed-11eb-002b-372884f598b0
# ╠═cbb9ef56-68ed-11eb-387b-252ff6dc5b21
# ╠═71a48f7a-68ee-11eb-1b0e-331393e95551
# ╠═b5b487ec-68ee-11eb-3ee7-0d98367028ff
# ╟─7310fa28-7132-11eb-32dc-e9aae12b4f9e
# ╠═84d0ef98-7132-11eb-1d64-3d911674a004
# ╠═8501f844-69b3-11eb-24f8-d3b87b745be7
# ╠═c2ba84b0-68f2-11eb-3c5a-17e03541dcb1
# ╠═9ee34a3a-6d4e-11eb-1e6e-b191bea5527b
# ╠═b05f0250-7009-11eb-0370-4969eb2fd5ce
# ╟─85a11cf4-6d40-11eb-3424-5b873aa5bc88
# ╠═c6955180-6d40-11eb-0723-e55c69316f63
# ╠═2c4b5ffc-68f4-11eb-21d6-7d2db834d3f8
# ╠═3e956c3c-6ee6-11eb-2aff-d73d8c0d9dcd
# ╠═468ae700-6ee6-11eb-289d-8dabcdca1b7d
# ╠═aedb425c-68f4-11eb-319d-a958c2aacde1
# ╠═f3afc7c2-6d68-11eb-10a6-65ab80be96a6
# ╠═bb70ce3a-6956-11eb-26ad-e19934980e11
# ╠═eb7c7c14-6956-11eb-3c98-3b4ee223bccc
# ╟─c3169ed8-694e-11eb-0616-97d25d8b56e1


text:
using LinearAlgebra, StaticArrays
import ProgressMeter
using DynamicalSystemsBase: MinimalDiscreteIntegrator

export lyapunovspectrum, lyapunov, local_growth_rates

#####################################################################################
#                               Lyapunov Spectum                                    #
#####################################################################################
"""
    lyapunovspectrum(ds::DynamicalSystem, N [, k::Int | Q0]; kwargs...) -> λs

Calculate the spectrum of Lyapunov exponents [^Lyapunov1992] of `ds` by applying
a QR-decomposition on the parallelepiped matrix `N` times. Return the
spectrum sorted from maximum to minimum.

The third argument `k` is optional, and dictates how many lyapunov exponents
to calculate (defaults to `dimension(ds)`).
Instead of passing an integer `k` you can pass
a pre-initialized matrix `Q0` whose columns are initial deviation vectors (then
`k = size(Q0)[2]`).

See also [`lyapunov`](@ref), [`local_growth_rates`](@ref).

## Keyword Arguments
* `u0 = get_state(ds)` : State to start from.
* `Ttr = 0` : Extra "transient" time to evolve the system before application of the
  algorithm. Should be `Int` for discrete systems. Both the system and the
  deviation vectors are evolved for this time.
* `Δt = 1` : Time of individual evolutions
  between successive orthonormalization steps. For continuous systems this is approximate.
* `show_progress = false` : Display a progress bar of the process.
* `diffeq` is a `NamedTuple` (or `Dict`) of keyword arguments propagated into
  `init` of DifferentialEquations.jl.
  See [`trajectory`](@ref) for examples. Only valid for continuous systems.

## Description
The method we employ is "H2" of [^Geist1990], originally stated in [^Benettin1980]. The deviation vectors
defining a `D`-dimensional parallepiped in tangent space
are evolved using the tangent dynamics of the system.
A QR-decomposition at each step yields the local growth rate for each dimension
of the parallepiped. The growth rates are
then averaged over `N` successive steps, yielding the lyapunov exponent spectrum
(at each step the parallepiped is re-normalized).

## Performance Notes
This function uses a [`tangent_integrator`](@ref).
For loops over initial conditions and/or
parameter values one should use the low level method that accepts
an integrator, and `reinit!` it to new initial conditions.
See the "advanced documentation" for info on the integrator object.
The low level method is
```julia
lyapunovspectrum(tinteg, N, Δt::Real, Ttr::Real)
```

If you want to obtain the convergence timeseries of the Lyapunov spectrum,
use the method
```julia
ChaosTools.lyapunovspectrum_convergence(tinteg, N, Δt, Ttr)
```
(not exported).

[^Lyapunov1992]: A. M. Lyapunov, *The General Problem of the Stability of Motion*, Taylor & Francis (1992)

[^Geist1990]: K. Geist *et al.*, Progr. Theor. Phys. **83**, pp 875 (1990)

[^Benettin1980]: G. Benettin *et al.*, Meccanica **15**, pp 9-20 & 21-30 (1980)
"""
lyapunovspectrum(ds::DS, N, k::Int = dimension(ds); kwargs...) =
lyapunovspectrum(ds, N, orthonormal(dimension(ds), k); kwargs...)

function lyapunovspectrum(ds::DS{IIP, S, D}, N, Q0::AbstractMatrix;
        Ttr::Real = 0, Δt::Real = 1, u0 = get_state(ds), show_progress = false,
        diffeq = NamedTuple(), kwargs...
    ) where {IIP, S, D}

    if !isempty(kwargs)
        @warn DIFFEQ_DEP_WARN
        diffeq = NamedTuple(kwargs)
    end

    if typeof(ds) <: DDS
        @assert typeof(Ttr) == Int
        integ = tangent_integrator(ds, Q0; u0)
    else
        integ = tangent_integrator(ds, Q0; u0, diffeq)
    end
    λ = lyapunovspectrum(integ, N, Δt, Ttr, show_progress)
    return λ
end

function lyapunovspectrum(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress = false)
    if show_progress
        progress = ProgressMeter.Progress(N; desc = "Lyapunov Spectrum: ", dt = 1.0)
    end
    B = copy(get_deviations(integ)) # for use in buffer
    if Ttr > 0
        t0 = integ.t
        while integ.t < t0 + Ttr
            step!(integ, Δt)
            Q, R = _buffered_qr(B, get_deviations(integ))
            set_deviations!(integ, Q)
        end
    end

    k = size(get_deviations(integ))[2]
    λ = zeros(stateeltype(integ), k)
    t0 = integ.t
    for i in 1:N
        step!(integ, Δt)
        Q, R = _buffered_qr(B, get_deviations(integ))
        for j in 1:k
            @inbounds λ[j] += log(abs(R[j,j]))
        end
        set_deviations!(integ, Q)
        show_progress && ProgressMeter.update!(progress, i)
    end
    λ ./= (integ.t - t0)
    return λ
end

# For out-of-place systems, this is just standard QR decomposition.
# For in-place systems, this is a more performant buffered version.
function _buffered_qr(B::SMatrix, Y) # Y are the deviations
    Q, R = LinearAlgebra.qr(Y)
    return Q, R
end
function _buffered_qr(B::Matrix, Y) # Y are the deviations
    B .= Y
    Q, R = LinearAlgebra.qr!(B)
    return Q, R
end

lyapunovspectrum(ds::DiscreteDynamicalSystem{IIP, T, 1}, a...; kw...) where {IIP, T} =
error("For discrete 1D systems, only method with state type = number is implemented.")

function lyapunovspectrum(ds::DDS{false, T, 1}, N; Ttr = 0) where {T}
    x = get_state(ds); f = ds.f
    p = ds.p; t0 = ds.t0
    if Ttr > 0
        for i in t0:(Ttr+t0)
            x = f(x, p, i)
        end
    end
    λ = zero(T)
    for i in (t0+Ttr):(t0+Ttr+N)
        x = f(x, p, i)
        λ += log(abs(ds.jacobian(x, p, i)))
    end
    return λ/N
end

#####################################################################################
#                           Lyapunov Spectrum Convergence                           #
#####################################################################################
"""
    lyapunovspectrum_convergence(ds::DynamicalSystem, N [, k::Int | Q0]; kwargs...) -> λs, t

Computes the convergence of the Lyapunov exponents, calculated exactly as in `lyapunovspectrum`.
(See [`lyapunovspectrum`](@ref) for information on the parameters and implementation).
Returns `λs`, a vector of `N` vectors, wherein vector of index `i` contains the Lyapunov
spectrum evaluated up to the time `t[i]`. The final vector is equal to the returned vector
of `lyapunovspectrum`.

## Performance Notes
This function uses a [`tangent_integrator`](@ref).
For loops over initial conditions and/or
parameter values one should use the low level method that accepts
an integrator, and `reinit!` it to new initial conditions.
See the "advanced documentation" for info on the integrator object.
The low level method is
```julia
lyapunovspectrum_convergence(tinteg, N, Δt::Real, Ttr::Real)
```
"""
lyapunovspectrum_convergence(ds::DS, N, k::Int = dimension(ds); kwargs...) =
lyapunovspectrum_convergence(ds, N, orthonormal(dimension(ds), k); kwargs...)

function lyapunovspectrum_convergence(ds::DS{IIP, S, D}, N, Q0::AbstractMatrix;
        Ttr::Real = 0, Δt::Real = 1, u0 = get_state(ds), show_progress = false, diffeq...
    ) where {IIP, S, D}

    if typeof(ds) <: DDS
        @assert typeof(Ttr) == Int
        integ = tangent_integrator(ds, Q0; u0)
    else
        integ = tangent_integrator(ds, Q0; u0, diffeq...)
    end
    λ = lyapunovspectrum_convergence(integ, N, Δt, Ttr, show_progress)
    return λ
end

function lyapunovspectrum_convergence(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress=false)
    if show_progress
        progress = ProgressMeter.Progress(N; desc = "Lyapunov Spectrum: ", dt = 1.0)
    end
    B = copy(get_deviations(integ)) # for use in buffer
    if Ttr > 0
        t0 = integ.t
        while integ.t < t0 + Ttr
            step!(integ, Δt)
            Q, R = _buffered_qr(B, get_deviations(integ))
            set_deviations!(integ, Q)
        end
    end

    k = size(get_deviations(integ))[2]
    T = stateeltype(integ)
    t0 = integ.t; t = zeros(T, N); t[1] = t0
    λs = [zeros(T, k) for i in 1:N];

    for i in 2:N
        step!(integ, Δt)
        Q, R = _buffered_qr(B, get_deviations(integ))
        for j in 1:k
            @inbounds λs[i][j] = λs[i-1][j] + log(abs(R[j,j]))
        end
        t[i] = integ.t
        set_deviations!(integ, Q)
        show_progress && ProgressMeter.update!(progress, i)
    end
    popfirst!(λs); popfirst!(t)
    for j in eachindex(t)
        λs[j] ./= (t[j] - t0)
    end
    return λs, t
end

lyapunovspectrum_convergence(ds::DiscreteDynamicalSystem{IIP, T, 1}, a...; kw...) where {IIP, T} =
error("For discrete 1D systems, only method with state type = number is implemented.")

function lyapunovspectrum_convergence(ds::DDS{false, T, 1}, N; Ttr = 0) where {T}
    x = get_state(ds); f = ds.f
    p = ds.p; t0 = ds.t0
    if Ttr > 0
        for i in t0:(Ttr+t0)
            x = f(x, p, i)
        end
    end
    t = (t0+Ttr):(t0+Ttr+N)
    λs = zeros(T, length(t));
    for idx = 2:length(t)
        i = t[idx]
        x = f(x, p, i)
        @inbounds λs[idx] = λs[idx-1] + log(abs(ds.jacobian(x, p, i)))
    end
    t = collect(t)
    popfirst!(λs); popfirst!(t)
    λs ./= t .- t0
    return λs, t
end

#####################################################################################
#                           Maximum Lyapunov Exponent                               #
#####################################################################################
inittest_default(D) = (state1, d0) -> state1 .+ d0/sqrt(D)

"""
    lyapunov(ds::DynamicalSystem, Τ; kwargs...) -> λ

Calculate the maximum Lyapunov exponent `λ` using a method due to Benettin [^Benettin1976],
which simply
evolves two neighboring trajectories (one called "given" and one called "test")
while constantly rescaling the test one.
`T`  denotes the total time of evolution (should be `Int` for discrete systems).

See also [`lyapunovspectrum`](@ref), [`local_growth_rates`](@ref).

## Keyword Arguments
* `u0 = get_state(ds)` : Initial condition.
* `Ttr = 0` : Extra "transient" time to evolve the trajectories before
  starting to measure the expontent. Should be `Int` for discrete systems.
* `d0 = 1e-9` : Initial & rescaling distance between the two neighboring trajectories.
* `upper_threshold = 1e-6` : Upper distance threshold for rescaling.
* `lower_threshold = 1e-12` : Lower distance threshold for rescaling (in order to
   be able to detect negative exponents).
* `Δt = 1` : Time of evolution between each check of
  distance exceeding the thresholds. For continuous
  systems this is approximate.
* `inittest = (u1, d0) -> u1 .+ d0/sqrt(D)` :
  A function that given `(u1, d0)`
  initializes the test state with distance
  `d0` from the given state `u1` (`D` is the dimension
  of the system). This function can be used when you want to avoid
  the test state appearing in a region of the phase-space where it would have
  e.g. different energy or escape to infinity.
* `diffeq` is a `NamedTuple` (or `Dict`) of keyword arguments propagated into
  `init` of DifferentialEquations.jl.
  See [`trajectory`](@ref) for examples. Only valid for continuous systems.

## Description
Two neighboring trajectories with initial distance `d0` are evolved in time.
At time ``t_i`` their distance ``d(t_i)`` either exceeds the `upper_threshold`,
or is lower than `lower_threshold`, which initializes
a rescaling of the test trajectory back to having distance `d0` from
the given one, while the rescaling keeps the difference vector along the maximal
expansion/contraction direction: `` u_2 \\to u_1+(u_2−u_1)/(d(t_i)/d_0)``.

The maximum
Lyapunov exponent is the average of the time-local Lyapunov exponents
```math
\\lambda = \\frac{1}{t_{n} - t_0}\\sum_{i=1}^{n}
\\ln\\left( a_i \\right),\\quad a_i = \\frac{d(t_{i})}{d_0}.
```

## Performance Notes
This function uses a [`parallel_integrator`](@ref).
For loops over initial conditions and/or
parameter values one should use the low level method that accepts
an integrator, and `reinit!` it to new initial conditions.
See the "advanced documentation" for info on the integrator object.
The low level method is
```
lyapunov(pinteg, T, Ttr, Δt, d0, ut, lt)
```

[^Benettin1976]: G. Benettin *et al.*, Phys. Rev. A **14**, pp 2338 (1976)
"""
function lyapunov(ds::DS, T;
        u0 = get_state(ds),
        Ttr = 0,
        d0=1e-9,
        upper_threshold = 1e-6,
        lower_threshold = 1e-12,
        inittest = inittest_default(dimension(ds)),
        Δt = 1,
        diffeq = NamedTuple(), kwargs...
    )

    if !isempty(kwargs)
        @warn DIFFEQ_DEP_WARN
        diffeq = NamedTuple(kwargs)
    end

    ST = stateeltype(ds)
    lower_threshold ≤ d0 ≤ upper_threshold || throw(ArgumentError(
    "d0 must be between thresholds!"))
    D = dimension(ds)
    if typeof(ds) <: DDS
        pinteg = parallel_integrator(ds, [deepcopy(u0), inittest(u0, d0)])
    else
        pinteg = parallel_integrator(ds, [deepcopy(u0), inittest(u0, d0)]; diffeq)
    end
    λ::ST = lyapunov(pinteg, T, Ttr, Δt, d0, upper_threshold, lower_threshold)
    return λ
end

function lyapunov(pinteg, T, Ttr, Δt, d0, ut, lt)
    # transient
    t0 = pinteg.t
    while pinteg.t < t0 + Ttr
        step!(pinteg, Δt)
        d = λdist(pinteg)
        lt ≤ d ≤ ut || rescale!(pinteg, d/d0)
    end

    t0 = pinteg.t
    d = λdist(pinteg)
    d == 0 && error("Initial distance between states is zero!!!")
    rescale!(pinteg, d/d0)
    λ = zero(d)
    while pinteg.t < t0 + T
        d = λdist(pinteg)
        #evolve until rescaling:
        while lt ≤ d ≤ ut
            step!(pinteg, Δt)
            d = λdist(pinteg)
            pinteg.t ≥ t0 + T && break
        end
        # local lyapunov exponent is simply the relative distance of the trajectories
        a = d/d0
        λ += log(a)
        rescale!(pinteg, a)
    end
    # Do final rescale, in case no other happened
    d = λdist(pinteg)
    a = d/d0
    λ += log(a)
    return λ/(pinteg.t - t0)
end

lyapunov(ds::DDS{false, T, 1}, N; Ttr = 0) where {T} = lyapunovspectrum(ds, N; Ttr = Ttr)

################ Helper functions that allow a single definition ######################
function λdist(integ::AbstractODEIntegrator{Alg, IIP, M}) where {Alg, IIP, M<:Matrix}
    d = 0.0
    for i in 1:size(integ.u)[1]
        d += (integ.u[i, 1] - integ.u[i, 2])^2
    end
    return sqrt(d)
    # return norm(view(integ.u, :, 1) .- view(integ.u, :, 2))
end
# No-annotation case is with Vectors
function λdist(integ)
    @inbounds s = zero(eltype(integ.u[1]))
    @inbounds for k in 1:length(integ.u[1])
        x = (integ.u[1][k] - integ.u[2][k])
        s += x*x
    end
    return sqrt(s)
end
function λdist(integ::MinimalDiscreteIntegrator{true, Vector{S}}) where {S<:SVector}
    return norm(integ.u[1] - integ.u[2])
end
function λdist(integ::AbstractODEIntegrator{Alg, IIP, Vector{S}}) where {Alg, IIP, S<:SVector}
    return norm(integ.u[1] - integ.u[2])
end

# Rescales:
function rescale!(integ::MinimalDiscreteIntegrator{true, Vector{S}}, a) where {S<:SVector}
    integ.u[2] = integ.u[1] + (integ.u[2] - integ.u[1])/a
    u_modified!(integ, true)
end
function rescale!(integ::MinimalDiscreteIntegrator{true, Vector{S}}, a) where {S<:Vector}
    @. integ.u[2] = integ.u[1] + (integ.u[2] - integ.u[1])/a
    u_modified!(integ, true)
end
function rescale!(integ::AbstractODEIntegrator{Alg, IIP, M}, a) where {Alg, IIP, M<:Matrix}
    for i in 1:size(integ.u)[1]
        integ.u[i, 2] = integ.u[i,1] + (integ.u[i,2] - integ.u[i,1])/a
    end
    u_modified!(integ, true)
end
function rescale!(integ::AbstractODEIntegrator{Alg, IIP, Vector{S}}, a) where {Alg, IIP, S<:Vector}
    @. integ.u[2] = integ.u[1] + (integ.u[2] - integ.u[1])/a
    u_modified!(integ, true)
end
function rescale!(integ::AbstractODEIntegrator{Alg, IIP, Vector{S}}, a) where {Alg, IIP, S<:SVector}
    integ.u[2] = integ.u[1] + (integ.u[2] - integ.u[1])/a
    u_modified!(integ, true)
end

#####################################################################################
#                              Local Growth Rates                                   #
#####################################################################################
"""
    local_growth_rates(ds, points::Dataset; S=100, Δt=5, kwargs...) → λlocal
Compute the exponential local growth rate(s) of perturbations of the dynamical system
`ds` for initial conditions given in `points`. For each initial condition `u ∈ points`,
`S` total perturbations are created and evolved for time `Δt`. The exponential local growth
rate is defined simply by `log(g/g0)/Δt` with `g0` the initial pertrubation size
and `g` the size after `Δt`. Thus, `λlocal` is a matrix of size `(length(points), S)`.

This function is a modification of [`lyapunov`](@ref). It uses the full nonlinear dynamics
to evolve the perturbations, but does not do any re-scaling, thus allowing
probing state and time dependence of perturbation growth. The actual growth
is given by `exp(λlocal * Δt)`.

The output of this function is sometimes referred as "Nonlinear Local Lyapunov Exponent".

## Keyword Arguments
* `perturbation`: If given, it should be a function `perturbation(ds, u, j)` that
  outputs a pertrubation vector (preferrably `SVector`) given the system, current initial
  condition `u` and the counter `j ∈ 1:S`. If not given, a random perturbation is
  generated with norm given by the keyword `e = 1e-6`.
* `diffeq` is a `NamedTuple` (or `Dict`) of keyword arguments propagated into
  `init` of DifferentialEquations.jl.
  See [`trajectory`](@ref) for examples. Only valid for continuous systems.
"""
function local_growth_rates(ds::DynamicalSystem, points;
        S = 100, Δt = 5, e = 1e-6,
        perturbation = (ds, u, j) -> _random_Q0(ds, u, j, e),
        diffeq = NamedTuple(), kwargs...
    )

    if !isempty(kwargs)
        @warn DIFFEQ_DEP_WARN
        diffeq = NamedTuple(kwargs)
    end

    λlocal = zeros(length(points), S)
    Q0 = perturbation(ds, points[1], 1)
    states = [points[1], points[1] .+ Q0]
    pinteg = parallel_integrator(ds, states; diffeq)

    for (i, u) in enumerate(points)
        for j in 1:S
            Q0 = perturbation(ds, u, j)
            states[1] = u
            states[2] = states[1] .+ Q0
            g0 = norm(Q0)
            reinit!(pinteg, states)
            step!(pinteg, Δt, true)
            g = norm(get_state(pinteg, 2) .- get_state(pinteg, 1))
            λlocal[i, j] = log(g/g0) / (pinteg.t - pinteg.t0)
        end
    end
    return λlocal
end

function _random_Q0(ds, u, j, e)
    D, T = dimension(ds), eltype(ds)
    Q0 = randn(Random.GLOBAL_RNG, SVector{D, T})
    Q0 = e * Q0 / norm(Q0)
end

#####################################################################################
#                    Numerical Lyapunov (from reconstruction)                       #
#####################################################################################
using Neighborhood, StaticArrays
using Distances: Metric, Cityblock, Euclidean
export NeighborNumber, WithinRange
export lyapunov_from_data
export Cityblock, Euclidean
# Everything in this section is based on Ulrich Parlitz [1]

"""
    lyapunov_from_data(R::Dataset, ks;  refstates, w, distance, ntype)
Return `E = [E(k) for k ∈ ks]`, where `E(k)` is the average logarithmic distance
between states of a neighborhood that are evolved in time for `k` steps
(`k` must be integer). The slope of `E` vs `k` approximate the maximum Lyapunov exponent,
see below.
Typically `R` is the result of delay coordinates of a single timeseries.

## Keyword Arguments

* `refstates = 1:(length(R) - ks[end])` : Vector of indices
  that notes which
  states of the reconstruction should be used as "reference states", which means
  that the algorithm is applied for all state indices contained in `refstates`.
* `w::Int = 1` : The [Theiler window](@ref).
* `ntype = NeighborNumber(1)` : The neighborhood type. Either [`NeighborNumber`](@ref)
  or [`WithinRange`](@ref). See [Neighborhoods](@ref) for more info.
* `distance::Metric = Cityblock()` : The distance function used in the
  logarithmic distance of nearby states. The allowed distances are `Cityblock()`
  and `Euclidean()`. See below for more info. The metric for finding neighbors is
  always the Euclidean one.

## Description
If the dataset exhibits exponential divergence of nearby states, then it should hold
```math
E(k) \\approx \\lambda\\cdot k \\cdot \\Delta t + E(0)
```
for a *well defined region* in the `k` axis, where ``\\lambda`` is
the approximated maximum Lyapunov exponent. ``\\Delta t`` is the time between samples in the
original timeseries. You can use [`linear_region`](@ref) with arguments `(ks .* Δt, E)` to
identify the slope (= ``\\lambda``) immediatelly, assuming you
have choosen sufficiently good `ks` such that the linear scaling region is bigger
than the saturated region.

The algorithm used in this function is due to Parlitz[^Skokos2016], which itself
expands upon Kantz [^Kantz1994]. In sort, for
each reference state a neighborhood is evaluated. Then, for each point in this
neighborhood, the logarithmic distance between reference state and neighborhood
state(s) is calculated as the "time" index `k` increases. The average of the above over
all neighborhood states over all reference states is the returned result.

If the `Metric` is `Euclidean()` then use the Euclidean distance of the
full `D`-dimensional points (distance ``d_E`` in ref.[^Skokos2016]).
If however the `Metric` is `Cityblock()`, calculate
the absolute distance of *only the first elements* of the `m+k` and `n+k` points
of `R` (distance ``d_F`` in ref.[^Skokos2016], useful when `R` comes from delay embedding).

[^Skokos2016]: Skokos, C. H. *et al.*, *Chaos Detection and Predictability* - Chapter 1 (section 1.3.2), Lecture Notes in Physics **915**, Springer (2016)

[^Kantz1994]: Kantz, H., Phys. Lett. A **185**, pp 77–87 (1994)
"""
function lyapunov_from_data(
        R::AbstractDataset{D, T}, ks;
        refstates = 1:(length(R) - ks[end]),
        w = 1,
        distance = Cityblock(),
        ntype = NeighborNumber(1),
    ) where {D, T}
    Ek = lyapunov_from_data(R, ks, refstates, Theiler(w), distance, ntype)
end

function lyapunov_from_data(
        R::AbstractDataset{D, T},
        ks::AbstractVector{Int},
        ℜ::AbstractVector{Int},
        theiler,
        distance::Metric,
        ntype::SearchType
    ) where {D, T}

    # ℜ = \Re<tab> = set of indices that have the points that one finds neighbors.
    # n belongs in ℜ and R[n] is the "reference state".
    # Thus, ℜ contains all the reference states the algorithm will iterate over.
    # ℜ is not estimated. it is given by the user. Most common choice:
    # ℜ = 1:(length(R) - ks[end])

    # ⩅(n) = \Cup<tab> = neighborhood of reference state n
    # which is evaluated for each n and for the given neighborhood type

    timethres = length(R) - ks[end]
    if maximum(ℜ) > timethres
        erstr = "Maximum index of reference states is > length(R) - ks[end] "
        erstr*= "and the algorithm cannot be performed on it. You have to choose "
        erstr*= "reference state indices of at most up to length(R) - ks[end]."
        throw(ArgumentError(erstr))
    end
    E = zeros(T, length(ks))
    E_n, E_m = copy(E), copy(E)
    tree = KDTree(R, Euclidean())
    skippedm = 0; skippedn = 0

    for n in ℜ
        # The ⋓(n) can be evaluated on the spot instead of being pre-calculated
        # for all reference states. Precalculating is faster, but allocates more memory.
        # Since ⋓[n] doesn't depend on `k` one can then interchange the loops:
        # Instead of k being the outermost loop, it becomes the innermost loop!
        point = R[n]
        ⋓ = isearch(tree, point, ntype, theiler(n))
        for m in ⋓
            # If `m` is nearer to the end of the timeseries than k allows
            # is it completely skipped (and length(⋓) reduced).
            if m > timethres
                skippedm += 1
                continue
            end
            for (j, k) in enumerate(ks) #ks should be small (of order 10 to 100 MAX)
                E_m[j] = log(delay_distance(distance, R, m, n, k))
            end
            E_n .+= E_m # no need to reset E_m
        end
        if skippedm >= length(⋓)
            skippedn += 1
            skippedm = 0
            continue # be sure to continue if no valid point!
        end
        E .+= E_n ./ (length(⋓) - skippedm)
        skippedm = 0
        fill!(E_n, zero(T)) #reset distances for n-th reference state
    end

    if skippedn >= length(ℜ)
        ers = "Skipped number of points ≥ length(R)...\n"
        ers*= "Could happen because all the neighbors fall within the Theiler "
        ers*= "window. Fix: increase neighborhood size."
        error(ers)
    end
    E ./= (length(ℜ) - skippedn)
end

@inline function delay_distance(::Cityblock, R, m, n, k)
    @inbounds abs(R[m+k][1] - R[n+k][1])
end

@inline function delay_distance(::Euclidean, R, m, n, k)
    @inbounds norm(R[m+k] - R[n+k])
end


text:
const state = gensym("state")

macro addr(expr::Expr, addr)
    if expr.head != :call
        error("syntax error in @addr at $(expr)")
    end
    fn = esc(expr.args[1])
    args = map(esc, expr.args[2:end])
    Expr(:call, :addr, esc(state), fn, Expr(:tuple, args...), esc(addr))
end

macro splice(expr::Expr)
    if expr.head != :call
        error("syntax error in @splice at $(expr)")
    end
    invocation = expr.args[1]
    args = esc(Expr(:tuple, expr.args[2:end]...))
    Expr(:call, :splice, esc(state), esc(invocation), args)
end

macro read(addr)
    Expr(:call, :read, esc(state), esc(addr))
end

marked_for_ad(arg::Symbol) = false

marked_for_ad(arg::Expr) = (arg.head == :macrocall && arg.args[1] == Symbol("@ad"))

strip_marked_for_ad(arg::Symbol) = arg

function strip_marked_for_ad(arg::Expr)
    if (arg.head == :macrocall && arg.args[1] == Symbol("@ad"))
		if length(arg.args) == 3 && isa(arg.args[2], LineNumberNode)
			arg.args[3]
		elseif length(arg.args) == 2
			arg.args[2]
		else
            error("Syntax error at $arg")
        end
    else
        arg
    end
end

##############################
# utility for examining code #
##############################

function strip_lineinfo(expr::Expr)
    @assert !(expr.head == :line)
    new_args = []
    for arg in expr.args
        if (isa(arg, Expr) && arg.head == :line) || isa(arg, LineNumberNode)
        elseif isa(arg, Expr) && arg.head == :block
            stripped = strip_lineinfo(arg)
            append!(new_args, stripped.args)
        else
            push!(new_args, strip_lineinfo(arg))
        end
    end
    Expr(expr.head, new_args...)
end

function strip_lineinfo(expr)
    expr
end

export strip_lineinfo

export @addr
export @read
export @splice

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen