text:
using CSV, DataFrames, StanMCMCChains, MCMCChains

ProjDir = @__DIR__
cd(ProjDir)

d = CSV.read("Kline.csv", delim=';');
size(d) # Should be 10x5

# New col log_pop, set log() for population data
d[:log_pop] = map((x) -> log(x), d[:population]);
d[:society] = 1:10;

first(d, 5)

m12_6sl_model = "
  data {
    int N;
    int T[N];
    int N_societies;
    int society[N];
    int P[N];
  }
  parameters {
    real alpha;
    vector[N_societies] a_society;
    real bp;
    real<lower=0> sigma_society;
  }
  model {
    vector[N] mu;
    target += normal_lpdf(alpha | 0, 10);
    target += normal_lpdf(bp | 0, 1);
    target += cauchy_lpdf(sigma_society | 0, 1);
    target += normal_lpdf(a_society | 0, sigma_society);
    for(i in 1:N) mu[i] = alpha + a_society[society[i]] + bp * log(P[i]);
    target += poisson_log_lpmf(T | mu);
  }
  generated quantities {
    vector[N] log_lik;
    {
    vector[N] mu;
    for(i in 1:N) {
      mu[i] = alpha + a_society[society[i]] + bp * log(P[i]);
      log_lik[i] = poisson_log_lpmf(T[i] | mu[i]);
    }
    }
  }
";

# Define the Stanmodel and set the output format to :mcmcchain.

stanmodel = Stanmodel(name="m12.6.2sl_model",  num_samples=4000,
model=m12_6sl_model);

# Input data for cmdstan

m12_6sl_data = Dict("N" => size(d, 1), "T" => d[:total_tools],
"N_societies" => 10, "society" => d[:society], "P" => d[:population]);

# Sample using cmdstan

rc, a3d, cnames = stan(stanmodel, m12_6sl_data, ProjDir,
diagnostics=false, summary=false, CmdStanDir=CMDSTAN_HOME);

# Describe the draws

cnames1 = cnames
cnames1[9:18] =  ["a_society[$i]" for i in 1:10]
cnames1[21:30] =  ["log_lik[$i]" for i in 1:10]

pi = filter(p -> length(p) > 2 && p[end-1:end] == "__", cnames1)
p = filter(p -> !(p in  pi), cnames1)
p1 = ["alpha", "bp", "sigma_society"]
p2 = filter(p -> !(p in p1), p)

m12_6sl = MCMCChains.Chains(a3d,
  cnames1,
  Dict(
    :parameters => p1,
    :pooled => cnames1[9:18],
    :log_lik => cnames1[21:30],
    :internals => pi
  )
)

write(joinpath(ProjDir, "sections_m12_6sl.jls"), m12_6sl)
open(joinpath(ProjDir, "sections_m12_6sl.txt"), "w") do io
  describe(io, m12_6sl, section=:pooled);
end

describe(m12_6sl)

describe(m12_6sl, section=:pooled)

describe(m12_6sl, section=:log_lik)


text:
"""
HeunConfluentLocalExpansion
    To Compute such a function one needs to decide the indices
    at each particlar point. There are 8 possible choices for this
    we must create an interface that allows people to choose which
    expansion they want
"""
struct HeunConfluentLocalExpansion
    p::Complex{Float64}; α::Complex{Float64}; γ::Complex{Float64};
    δ::Complex{Float64}; σ::Complex{Float64};
    acoeffs::Array{Complex{Float64},1}
    rcoeffs::Array{Complex{Float64},1}
end

function ComputeSeries(p,α,γ,δ,σ)
    αₖ(k) = -(k + 1)*(k + γ);
    βₖ(k) = k*(k - 4*p + γ + δ - 1);
    γₖ(k) = 4*p*(k + α - 1);
    an(k) = βₙ(k)/αₙ(k) ;
    bn(k) = γₙ(k)/αₙ(k);
    ComputeSeriesFromab(an,bn; N=250)
end

ext: jl


text:
# Autogenerated wrapper script for GAP_jll for aarch64-linux-gnu-cxx11-julia_version+1.7.0
export gap, libgap

using GMP_jll
using Readline_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("GAP")
JLLWrappers.@declare_executable_product(gap)
JLLWrappers.@declare_library_product(libgap, "libgap.so.8")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll, Readline_jll, Zlib_jll)
    JLLWrappers.@init_executable_product(
        gap,
        "bin/gap",
    )

    JLLWrappers.@init_library_product(
        libgap,
        "lib/libgap.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()

    sym = dlsym(libgap_handle, :GAP_InitJuliaMemoryInterface)
    ccall(sym, Nothing, (Any, Ptr{Nothing}), @__MODULE__, C_NULL)

end  # __init__()


text:
@testset "Continuous space" begin

  # Basic model initialization
  space1 = ContinuousSpace((1, 1), 0.1; periodic = true)
  space2 = ContinuousSpace((1, 1), 0.1; periodic = false)
  space3 = ContinuousSpace((2, 1), 0.1; periodic = true)
  space4 = ContinuousSpace((1, 2), 0.1; periodic = true)
  space5 = ContinuousSpace((2, 2), 0.1; periodic = true)
  space6 = ContinuousSpace((1, ), 0.1; periodic = true)
  space7 = ContinuousSpace((1, 1, 1), 0.1; periodic = true)
  space8 = ContinuousSpace((1, 1, 1), 0.1; periodic = false)

  @test eltype(space1) == Float64
  @test length(space1.dims) == 2
  @test length(space2.dims) == 2
  @test length(space3.dims) == 2
  @test length(space4.dims) == 2
  @test length(space5.dims) == 2
  @test length(space6.dims) == 1
  @test length(space7.dims) == 3
  @test length(space8.dims) == 3

  @test_throws ArgumentError ContinuousSpace((-1,1), 0.1) # Cannot have negative extent
  @test_throws MethodError ContinuousSpace([1,1], 0.1) # Must be a tuple
  @test_throws MethodError ContinuousSpace(("one",1.0), 0.1) # Must contain reals

  @test space1.dims == ContinuousSpace((1,1)).dims
  @test space1.extent == ContinuousSpace((1,1)).extent
  model1 = ABM(Agent6, space1)
  model2 = ABM(Agent6, space2)

  @test nagents(model1) == 0
  @test model1.agents == Dict()

  # add_agent! with no existing agent (the agent is created)
  pos = (0.5, 0.5)
  vel = (0.2, 0.1)
  dia = 0.01
  add_agent!(pos, model2, vel, dia)
  @test collect(keys(model2.agents)) == [1]
  @test model2[1].pos == (0.5, 0.5)

  # move_agent! without provided update_vel! function
  move_agent!(model2[1], model2)
  @test model2[1].pos == (0.7, 0.6)

  kill_agent!(model2.agents[1], model2)
  @test length(model2.agents) == 0

  # add_agent! with an existing agent
  agent = Agent6(2, pos, vel, dia)
  add_agent!(agent, model2)
  @test Agents.defvel(agent, model2) === nothing
  @test collect(keys(model2.agents)) == [2]

  # agents within some range are found correctly
  agent2 = model2[2]
  new_pos = min.(agent2.pos .+ 0.0051, model2.space.extent.-0.0001)
  agent3 = Agent6(3, new_pos, vel, dia)
  add_agent_pos!(agent3, model2)
  n_ids = collect(nearby_ids(agent2, model2, agent2.weight, exact=true))
  @test n_ids == [3]
  n_ids = collect(nearby_ids(agent2, model2, agent2.weight/2, exact=true))
  @test n_ids == []

  # test that it finds both
  n_ids = collect(nearby_ids(agent2.pos, model2, agent2.weight; exact=true))
  @test sort!(n_ids) == [2, 3]

end

mutable struct AgentU1 <: AbstractAgent
  id::Int
  pos::NTuple{2,Float64}
  vel::NTuple{2,Float64}
end

mutable struct AgentU2 <: AbstractAgent
  id::Int
  pos::NTuple{2,Float64}
  vel::NTuple{2,Float64}
end

function ignore_six(model::ABM)
  [a.id for a in allagents(model) if !(typeof(a) <: Agent6)]
end

@testset "Interacting pairs" begin
  space = ContinuousSpace((10, 10), 0.2, periodic = false)
  model = ABM(Agent6, space; scheduler = model -> sort!(collect(keys(model.agents));rev=true))
  pos = [
    (7.074386436066224, 4.963014649338054)
    (5.831962448496828, 4.926297135685473)
    (5.122087781793935, 5.300031210394806)
    (3.9715633336430156, 4.8106570045816675)
  ]
  for i in 1:4
    add_agent_pos!(Agent6(i+2, pos[i], (0.0, 0.0), 0), model)
  end
  pairs = interacting_pairs(model, 2.0, :all).pairs
  @test length(pairs) == 5
  @test (3, 6) ∉ pairs

  space2 = ContinuousSpace((10, 10), 0.1,periodic = false)
  model2 = ABM(Agent6, space2; scheduler = Schedulers.by_id)
  for i in 1:4
    add_agent_pos!(Agent6(i, pos[i], (0.0, 0.0), 0), model2)
  end
  pairs = interacting_pairs(model2, 2.0, :nearest).pairs
  @test length(pairs) == 1
  pairs = interacting_pairs(model2, 2.5, :all).pairs
  @test length(pairs) == 5
  @test (1, 4) ∉ pairs

  space3 = ContinuousSpace((10,10), 1.0, periodic = false)
  model3 = ABM(Union{Agent6, AgentU1, AgentU2}, space3; warn = false)
  for i in 1:10
    add_agent_pos!(Agent6(i, (i/10, i/10), (0.0, 0.0), 0), model3)
  end
  for i in 11:20
    add_agent_pos!(AgentU1(i, (i/10-1, 0.5), (0.0, 0.0)), model3)
  end
  for i in 21:30
    add_agent_pos!(AgentU2(i, (0.45, i/10-2), (0.0, 0.0)), model3)
  end
  pairs = interacting_pairs(model3, 0.1, :types).pairs
  @test length(pairs) == 7
  for (a,b) in pairs
      @test typeof(model3[a]) !== typeof(model3[b])
  end
  @test (3, 6) ∉ pairs

  # Test that we have at least some Agent6's in this match
  @test any(typeof(model3[a]) <: Agent6 || typeof(model3[b]) <: Agent6 for (a,b) in pairs)
  pairs = interacting_pairs(model3, 0.2, :types; scheduler = ignore_six).pairs
  @test length(pairs) == 12
  # No Agent6's when using the ignore_six scheduler
  @test all(!(typeof(model3[a]) <: Agent6) && !(typeof(model3[b]) <: Agent6) for (a,b) in pairs)

  # Fix #288
  space = ContinuousSpace((1,1), 0.1; periodic = true)
  model = ABM(Agent6, space)
  pos = [(0.01, 0.01),(0.2,0.2),(0.5,0.5)]
  for i in pos
    add_agent!(i,model,(0.0,0.0),1.0)
  end
  pairs = collect(interacting_pairs(model, 0.29, :all))
  @test length(pairs) == 1
  (a,b) = first(pairs)
  @test (a.id, b.id) == (1,2)
  # Before the #288 fix, this would return (2,3) as a pair
  # which has a euclidean distance of 0.42
  pairs = collect(interacting_pairs(model, .3, :all))
  @test length(pairs) == 1
  (a,b) = first(pairs)
  @test (a.id, b.id) == (1,2)
end

@testset "nearest neighbor" begin
  space = ContinuousSpace((1,1), 0.1; periodic = true)
  model = ABM(Agent9, space)
  pos = [(0.01, 0.01),(0.2, 0.01),(0.2, 0.2),(0.5, 0.5)]
  for i in pos
    add_agent!(i,model,(0.0,0.0),nothing)
  end

  for agent in allagents(model)
    agent.f1 = nearest_neighbor(agent, model, sqrt(2)).id
  end

  @test model.agents[1].f1 == 2
  @test model.agents[2].f1 == 1
  @test model.agents[3].f1 == 2
  @test model.agents[4].f1 == 3

  @test size(collect(positions(model))) == (10, 10)
  @test ids_in_position((3, 3), model) == [3]
  @test sort!(collect(nearby_ids(model[3].pos, model, 0.2))) == [2, 3]
  @test Agents.distance_from_cell_center((0.22, 0.22), model) ≈ 0.04242640687
end


text:
__precompile__()

"""
The Conda module provides access to the [conda](http://conda.pydata.org/) packages
manager. Its main purpose is to be used as a BinDeps provider, to install binary
dependencies of other Julia packages.

The main functions in Conda are:

- `Conda.add(package)`: install a package;
- `Conda.rm(package)`: remove (uninstall) a package;
- `Conda.update()`: update all installed packages to the latest version;
- `Conda.list()`: list all installed packages.
- `Conda.add_channel(channel)`: add a channel to the list of channels;
- `Conda.channels()`: get the current list of channels;
- `Conda.rm_channel(channel)`: remove a channel from the list of channels;

To use Conda as a binary provider for BinDeps, the `Conda.Manager` type is proposed. A
small example looks like this:

```julia
# Declare dependency
using BinDeps
@BinDeps.setup
netcdf = library_dependency("netcdf", aliases = ["libnetcdf","libnetcdf4"])

using Conda
#  Use alternative conda channel.
Conda.add_channel("my_channel")
provides(Conda.Manager, "libnetcdf", netcdf)
```
"""
module Conda
using Compat
import Compat.String
using JSON

"Prefix for installation of all the packages."
const PREFIX = abspath(dirname(@__FILE__), "..", "deps", "usr")

"Prefix for the executable files installed with the packages"
const BINDIR = is_windows() ? joinpath(PREFIX, "Library", "bin") : joinpath(PREFIX, "bin")

"Prefix for the python scripts. On UNIX, this is the same than Conda.BINDIR"
const SCRIPTDIR = is_windows() ? joinpath(PREFIX, "Scripts") : BINDIR

"Prefix where the `python` command lives"
const PYTHONDIR = is_windows() ? PREFIX : BINDIR

"Path to the `conda` binary"
const conda = joinpath(SCRIPTDIR, "conda")

"Path to the condarc file"
const CONDARC = joinpath(PREFIX, "condarc-julia")

"""
Use a cleaned up environment for the command `cmd`.

Any environment variable starting by CONDA will interact with the run.
"""
function _set_conda_env(cmd)
    env = copy(ENV)
    to_remove = AbstractString[]
    for var in keys(env)
        if startswith(var, "CONDA")
            push!(to_remove, var)
        end
    end
    for var in to_remove
        pop!(env, var)
    end
    env["CONDARC"] = CONDARC
    setenv(cmd, env)
end

"Get the miniconda installer URL."
function _installer_url()
    res = "http://repo.continuum.io/miniconda/Miniconda-latest-"
    if is_apple()
        res *= "MacOSX"
    elseif is_linux()
        res *= "Linux"
    elseif is_windows()
        res *= "Windows"
    else
        error("Unsuported OS.")
    end
    res *= Sys.WORD_SIZE == 64 ? "-x86_64" : "-x86"
    res *= is_windows() ? ".exe" : ".sh"
    return res
end

is_windows() && include("outlook.jl")

"Install miniconda if it hasn't been installed yet; _install_conda(true) installs Conda even if it has already been installed."
function _install_conda(force=false)
    if is_windows()
          if is_outlook_running()
              error("""\n
              Outlook is running, and running the Miniconda installer will crash it.
              Please quit Outlook and then restart the installation.

              For more information, see:
                    https://github.com/Luthaf/Conda.jl/issues/15
                    https://github.com/conda/conda/issues/1084
              """)
          end
    end

    if force || !(is_windows() ? isfile(Conda.conda * ".exe") : isfile(Conda.conda))
        # Ensure PREFIX exists
        mkpath(PREFIX)
        info("Downloading miniconda installer ...")
        if is_unix()
            installer = joinpath(PREFIX, "installer.sh")
        end
        if is_windows()
            installer = joinpath(PREFIX, "installer.exe")
        end
        download(_installer_url(), installer)

        info("Installing miniconda ...")
        if is_unix()
            chmod(installer, 33261)  # 33261 corresponds to 755 mode of the 'chmod' program
            run(`$installer -b -f -p $PREFIX`)
        end
        if is_windows()
            if VERSION >= v"0.5.0-dev+8873" # Julia PR #13780
                run(Cmd(`$installer /S /AddToPath=0 /RegisterPython=0 /D=$PREFIX`, windows_verbatim=true))
            else
                # Workaround a bug in command-line argument parsing, see
                # https://github.com/Luthaf/Conda.jl/issues/17
                if match(r" {2,}", "") != nothing
                    error("The installer will fail when the path=\"$PREFIX\" contains two consecutive spaces")
                end
                run(`$installer /S /AddToPath=0 /RegisterPython=0 $(split("/D=$PREFIX"))`)
            end
        end
        Conda.add_channel("defaults")
    end
end

"Install a new package."
function add(pkg::AbstractString)
    _install_conda()
    run(_set_conda_env(`$conda install -y $pkg`))
end

"Uninstall a package."
function rm(pkg::AbstractString)
    _install_conda()
    run(_set_conda_env(`$conda remove -y $pkg`))
end

"Update all installed packages."
function update()
    _install_conda()
    for package in _installed_packages()
        run(_set_conda_env(`$conda update -y $package`))
    end
end

"List all installed packages as an dict of tuples with (version_number, fullname)."
function  _installed_packages_dict()
    _install_conda()
    package_dict = Dict{Compat.UTF8String, Tuple{VersionNumber, Compat.UTF8String}}()
    for line in eachline(`$(Conda.conda) list --export`)
        line = chomp(line)
        if !startswith(line, "#")
            name, version, build_string = split(line, "=")
            # As julia do not accepts xx.yy.zz.rr version number the last part is removed.
            # see issue https://github.com/JuliaLang/julia/issues/7282 a maximum of three levels is inserted
            version_number = join(split(version,".")[1:min(3,end)],".")
            try
                package_dict[name] = (convert(VersionNumber, version_number), line)
            catch
                package_dict[name] = (v"9999.9999.9999", line)
                warn("Failed parsing string: \"$(version_number)\" to a version number. Please open an issue!")
            end
        end
    end
    return package_dict
end

"List all installed packages as an array."
_installed_packages() = keys(_installed_packages_dict())

"List all installed packages to standard output."
function list()
    _install_conda()
    run(_set_conda_env(`$conda list`))
end

"Get the exact version of a package."
function version(name::AbstractString)
    _install_conda()
    packages = JSON.parse(readstring(`$conda list --json`))
    for package in packages
        if startswith(package, name) || ismatch(Regex("::$name"), package)
            return package
        end
    end
    error("Could not find the $name package")
end

"Search packages for a string"
function search(package::AbstractString)
    _install_conda()
    return collect(keys(JSON.parse(readstring(_set_conda_env(`$conda search $package --json`)))))
end

"Search a specific version of a package"
function search(package::AbstractString, ver::AbstractString)
    _install_conda()
    ret=JSON.parse(readstring(_set_conda_env(`$conda search $package --json`)))
    out = Compat.ASCIIString[]
    for k in keys(ret)
      for i in 1:length(ret[k])
        ret[k][i]["version"]==ver && push!(out,k)
      end
    end
    out
end

"Check if a given package exists."
function exists(package::AbstractString)
    if contains(package,"==")
      pkg,ver=split(package,"==")  # Remove version if provided
      return pkg in search(pkg,ver)
    else
      if package in search(package)
        # Found exactly this package
        return true
      else
        return false
      end
    end
end

"Get the list of channels used to search packages"
function channels()
    _install_conda()
    ret=JSON.parse(readstring(_set_conda_env(`$conda config --get channels --json`)))
    if haskey(ret["get"], "channels")
        return collect(Compat.ASCIIString, ret["get"]["channels"])
    else
        return Compat.ASCIIString[]
    end
end

"Add a channel to the list of channels"
function add_channel(channel::Compat.String)
    _install_conda()
    run(_set_conda_env(`$conda config --add channels $channel --force`))
end

"Remove a channel from the list of channels"
function rm_channel(channel::Compat.String)
    _install_conda()
    run(_set_conda_env(`$conda config --remove channels $channel --force`))
end

include("bindeps_conda.jl")

end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using Test
using FMIExport

include("model_description.jl")
include("export_simple.jl")


text:
# utility functions

function make_blocks(nobs, nblocks)
    maxbs, reminder = divrem(nobs, nblocks)

    res = fill(maxbs, nblocks)
    if reminder > 0
        res[1:reminder] = maxbs + 1
    end
    res
end

# the truncated newton method for matrix inversion
# adapted from https://en.wikipedia.org/wiki/Conjugate_gradient_method
# solves Ax = b for x, overwrites x
# stops if the error is < ɛ or after reaching max_iter
function truncated_newton!{T}(A::Matrix{T}, b::Vector{T},
                              x::Vector{T}, ɛ::T, max_iter::Int)
    r = b - A*x
    p = deepcopy(r)
    Ap = deepcopy(r)
    rsold = dot(r, r)

    n = length(r)

    for i in 1:max_iter
        # Ap[:] = A * p
        A_mul_B!(Ap, A, p)
        α = rsold / dot(p, Ap)
        # x += α * p
        BLAS.axpy!(α, p, x)
        # r -= α * Ap
        BLAS.axpy!(-α, Ap, r)
        rsnew = dot(r, r)
        rsnew < ɛ &&  break
        β = rsnew / rsold
        # p[:] = r + β * p
        BLAS.scal!(n, β, p, 1)
        BLAS.axpy!(1, r, p)
        rsold = rsnew
    end
    return x
end


text:
using BudanTables
using Test

@testset "BudanTables.jl" begin
    # Write your tests here.
end


text:
using DiffEqFlux, OrdinaryDiffEq, DiffEqSensitivity
using CUDA, Test, Zygote, Random, LinearAlgebra

CUDA.allowscalar(false)

H = CuArray(rand(Float32, 2, 2))
ann = FastChain(FastDense(1, 4, tanh))
p = initial_params(ann)

function func(x, p, t)
    ann([t],p)[1]*H*x
end

x0 = CuArray(rand(Float32, 2))
x1 = CuArray(rand(Float32, 2))

prob = ODEProblem(func, x0, (0.0f0, 1.0f0))

function evolve(p)
    solve(prob, Tsit5(), p=p, save_start=false,
          save_everystep=false, abstol=1e-4, reltol=1e-4,
          sensealg=QuadratureAdjoint(autojacvec=ZygoteVJP())).u[1]
end

function cost(p)
    x = evolve(p)
    c = sum(abs,x - x1)
    #println(c)
    c
end

grad = Zygote.gradient(cost,p)[1]
@test !iszero(grad[1])
@test iszero(grad[2:4])
@test !iszero(grad[5])
@test iszero(grad[6:end])

###
# https://github.com/SciML/DiffEqFlux.jl/issues/632
###

rng = MersenneTwister(1234)
m = 32
n = 16
Z = randn(rng, Float32, (n,m)) |> gpu
𝒯 = 2.0
Δτ = 0.1
ca_init = [zeros(1) ; ones(m)] |> gpu

function f(ca, Z, t)
  a = ca[2:end]

  a_unit = a / sum(a)
  w_unit = Z*a_unit
  Ka_unit = Z'*w_unit
  z_unit = dot(abs.(Ka_unit), a_unit)
  aKa_over_z = a .* Ka_unit / z_unit
  [sum(aKa_over_z) / m; -abs.(aKa_over_z)] |> gpu
end

function c(Z)
  prob = ODEProblem(f, ca_init, (0.,𝒯), Z, saveat=Δτ)
  sol = solve(prob, Tsit5(), sensealg=BacksolveAdjoint(), saveat=Δτ)
  sum(last(sol.u))
end

println("forward:", c(Z))
println("backward: ", Zygote.gradient(c, Z))


text:
module ExtraTrees

export extraTrees
export predict

abstract BinaryTree

type Node <: BinaryTree
  left::BinaryTree
  right::BinaryTree
  col::Int32
  cut::Float64
end

type Leaf <: BinaryTree
  value::Float32
end

type RegressionET
  trees::Vector{BinaryTree}
  ntry::Int64
  nodesize::Int64
  nrandomcuts::Int64
end

type Sampler
  x::Vector{Int64}
  nextIndex::Int64
  Sampler(x) = new(x, 1)
end

function hasNext(s::Sampler)
  return s.nextIndex <= length(s.x)
end

function reset!(s::Sampler)
  s.nextIndex = 1
end

function next!(s::Sampler)
  i = rand( 0:(length(s.x) - s.nextIndex) )
  tmp = s.x[s.nextIndex + i]
  s.x[s.nextIndex + i] = s.x[s.nextIndex]
  s.x[s.nextIndex]     = tmp
  s.nextIndex += 1
  return tmp
end

type RegressionData{T}
  x::T
  y::Vector{Float64}
  y2::Vector{Float64}
  N::Int

  RegressionData(x,y) = new(x, y, y.^2, length(y))
end

type VarScorer
  ysum::Float64
  y2sum::Float64
  n::Int64
end

function init!(scorer::VarScorer, data::RegressionData, ids)
  scorer.ysum  = 0.0
  scorer.y2sum = 0.0
  scorer.n     = length(ids)
  for id = ids
    scorer.ysum  += data.y[id]
    scorer.y2sum += data.y2[id]
  end
end

function init!(scorer::VarScorer, y, y2, ids)
  scorer.ysum  = sum(y)
  scorer.y2sum = sum(y2)
  scorer.n     = length(ids)
end

function isConstant(scorer::VarScorer)
  scorer.y2sum / scorer.n - (scorer.ysum / scorer.n) ^ 2 < 1e-7
end

## return variance for the given cut
function getScore(scorer::VarScorer, data::RegressionData, ids, col::Int64, cut::Float64)
  ## finding columns var
  ysumLeft  = 0.0
  y2sumLeft = 0.0
  countLeft = 0

  for id in ids
    if data.x[id, col] < cut
      ysumLeft  += data.y[id]
      y2sumLeft += data.y2[id]
      countLeft += 1
    end
  end
  countRight = length(ids) - countLeft
  varLeft    = y2sumLeft / countLeft - (ysumLeft / countLeft)^2
  varRight   = (scorer.y2sum - y2sumLeft) / countRight - ((scorer.ysum - ysumLeft) / countRight) ^ 2
  var = countLeft * varLeft + countRight * varRight
  return var
end

function extrema(x::Matrix{Float64}, ids, col::Int64)
  min = Inf
  max = -Inf
  for id = ids
    xv = x[id,col]
    if xv < min
      min = xv
    end
    if xv > max
      max = xv
    end
  end
  (min,max)
end

function trainTree(et::RegressionET,
                   data::RegressionData{SparseMatrixCSC{Bool,Int64}},
                   ids,
                   sampler::Sampler,
                   scorer::VarScorer)
  if length(ids) < et.nodesize
    return Leaf(mean(data.y[ids]))
  end
  Nsamples = size(data.x, 1)
  sids = sparsevec(ids, 1, Nsamples)'
  y    = sparsevec(ids, data.y[ids], Nsamples)'
  y2   = sparsevec(ids, data.y2[ids], Nsamples)'
  ## TODO
  init!(scorer, y, y2, ids)
end

function trainTree(et::RegressionET, data::RegressionData, ids, sampler::Sampler, scorer::VarScorer)
  if length(ids) < et.nodesize
    return Leaf(mean(data.y[ids]))
  end

  ## initialize scorer for current node
  init!(scorer, data, ids)

  if isConstant(scorer)
    return Leaf(scorer.ysum / length(ids))
  end

  bestVar = Inf
  bestCol = -1
  bestCut = 0.0

  ## random columns
  reset!(sampler)

  ## check K columns
  K = 0
  while hasNext(sampler)
    col = next!(sampler)
    range = extrema(data.x, ids, col)
    if range[1] == range[2]
      continue
    end
    for ncut = 1:et.nrandomcuts
      ## random cut
      cut = range[1] + (range[2] - range[1]) * rand()

      var = getScore(scorer, data, ids, col, cut)
      if var < bestVar
        bestVar = var
        bestCol = col
        bestCut = cut
      end
    end
    K += 1
    if K >= et.ntry
      break
    end
  end
  ## splitting samples
  idsLeft  = Int64[]
  idsRight = Int64[]
  for i = ids
    if data.x[i, bestCol] < bestCut
      push!(idsLeft, i)
    else
      push!(idsRight, i)
    end
  end

  leftNode  = trainTree(et, data, idsLeft, sampler, scorer)
  rightNode = trainTree(et, data, idsRight, sampler, scorer)

  return Node(leftNode, rightNode, bestCol, bestCut)
end

function extraTrees(x, y::Vector{Float64};
                            ntree     = 500,
                            ntry      = int( max(1, floor(size(x,1) / 3)) ),
                            nrandomcuts = 1,
                            nodesize  = 5)
  data  = RegressionData{typeof(x)}(x, y)
  trees = BinaryTree[]
  et    = RegressionET(trees, ntry, nodesize, nrandomcuts)
  scorer  = VarScorer(0.0, 0.0, 0)
  sampler = Sampler([1:size(x,2)])
  for i in 1:ntree
    push!(trees, trainTree(et, data, 1:length(y), sampler, scorer) )
  end
  return et
end

function predict(tree::Leaf, x::Vector{Float64})
  return tree.value
end

function predict(tree::Node, x::Vector{Float64})
  if x[tree.col] < tree.cut
    return predict(tree.left, x)
  else
    return predict(tree.right, x)
  end
end

function predict(et::RegressionET, x::Vector{Float64})
  sum = 0.0
  for tree = et.trees
    sum += predict(tree, x)
  end
  return sum / length(et.trees)
end

function predict(et::RegressionET, Xnew::Matrix{Float64})
  map(
    i -> predict(et, vec(Xnew[i,:])),
    1:size(Xnew, 1)
  )
end

end # module


text:
const ERR_TRAIN = ArgumentError("`IterationControl.train!` needs at "*
                                "least two arguments. ")

const ERR_NEEDS_LOSS = ArgumentError(
    "Encountered a control that needs losses but no losses found. ")

const ERR_NEEDS_TRAINING_LOSSES = ArgumentError(
    "Encountered a control that needs training losses but no training "*
    "losses found. ")

function train!(model, controls...; verbosity::Int=1)

    isempty(controls) && throw(ERR_TRAIN)

    control = composite(controls...)

    # before training:
    verbosity > 1 && @info "Using these controls: $(flat(control)). "

    # first training event:
    n = 1 # counts control cycles
    state = update!(control, model, verbosity, n)
    finished = done(control, state)

    # checks that model supports control:
    if needs_loss(control) && loss(model) === nothing
        throw(ERR_NEEDS_LOSS)
    end
    if needs_training_losses(control) && training_losses(model) === nothing
        throw(ERR_NEEDS_TRAINING_LOSSES)
    end

    # subsequent training events:
    while !finished
        n += 1
        state = update!(control, model, verbosity, n, state)
        finished = done(control, state)
    end

    # reporting final loss and training loss if available:
    _loss = IterationControl.loss(model)
    _training_losses = IterationControl.training_losses(model)
    if verbosity > 0
        _loss isa Nothing || @info "final loss: $_loss"
        _training_losses isa Nothing || isempty(_training_losses) ||
            @info "final training loss: $(_training_losses[end])"
        verbosity > 1 && @info "total control cycles: $n"
    end

    # finalization:
    return takedown(control, verbosity, state)
end


text:
for f in (:getindex, :view, :dotview)
    @eval begin
        # Int and CartesianIndex forward to the parent
        @propagate_inbounds Base.$f(l::LookupArray, i::Union{Int,CartesianIndex}) =
            Base.$f(parent(l), i)
        # AbstractArray and Colon the lookup is rebuilt around a new parent
        @propagate_inbounds Base.$f(l::LookupArray, i::Union{AbstractArray,Colon}) =
            rebuild(l; data=Base.$f(parent(l), i))
        # Selector gets processed with `selectindices`
        @propagate_inbounds Base.$f(l::LookupArray, i::SelectorOrInterval) = Base.$f(l, selectindices(l, i))
        # Everything else (like custom indexing from other packages) passes through to the parent
        @propagate_inbounds function Base.$f(l::LookupArray, i)
            Base.$f(parent(l), i)
        end
    end
end


text:
# Autogenerated wrapper script for OpenFAST_jll for powerpc64le-linux-gnu-libgfortran4
export aerodyn_driver, beamdyn_driver, dwm_driver_wind_farm, feam_driver, hydrodyn_driver, inflowwind_driver, libaeroacoustics, libaerodyn14lib, libaerodynlib, libafinfolib, libawaelib, libbeamdynlib, libelastodynlib, libextptfm_mckflib, libfeamlib, libfoamfastlib, libfvwlib, libhydrodynlib, libicedynlib, libicefloelib, libifwlib, libmapcpplib, libmaplib, libmoordynlib, libnwtclibs, libopenfast_postlib, libopenfast_prelib, libopenfastlib, libopenfoamtypeslib, liborcaflexlib, libscdataexlib, libscdataexttypeslib, libscfastlib, libsctypeslib, libservodynlib, libsubdynlib, libuaaerolib, libversioninfolib, libwdlib, moordyn_driver, openfast, orca_driver, servodyn_driver, subdyn_driver, turbsim, unsteadyaero_driver

using OpenBLAS32_jll
using LAPACK_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("OpenFAST")
JLLWrappers.@declare_library_product(libaeroacoustics, "libaeroacoustics.so")
JLLWrappers.@declare_library_product(libaerodyn14lib, "libaerodyn14lib.so")
JLLWrappers.@declare_library_product(libaerodynlib, "libaerodynlib.so")
JLLWrappers.@declare_library_product(libafinfolib, "libafinfolib.so")
JLLWrappers.@declare_library_product(libawaelib, "libawaelib.so")
JLLWrappers.@declare_library_product(libbeamdynlib, "libbeamdynlib.so")
JLLWrappers.@declare_library_product(libelastodynlib, "libelastodynlib.so")
JLLWrappers.@declare_library_product(libextptfm_mckflib, "libextptfm_mckflib.so")
JLLWrappers.@declare_library_product(libfeamlib, "libfeamlib.so")
JLLWrappers.@declare_library_product(libfoamfastlib, "libfoamfastlib.so")
JLLWrappers.@declare_library_product(libfvwlib, "libfvwlib.so")
JLLWrappers.@declare_library_product(libhydrodynlib, "libhydrodynlib.so")
JLLWrappers.@declare_library_product(libicedynlib, "libicedynlib.so")
JLLWrappers.@declare_library_product(libicefloelib, "libicefloelib.so")
JLLWrappers.@declare_library_product(libifwlib, "libifwlib.so")
JLLWrappers.@declare_library_product(libmapcpplib, "libmapcpplib.so")
JLLWrappers.@declare_library_product(libmaplib, "libmaplib.so")
JLLWrappers.@declare_library_product(libmoordynlib, "libmoordynlib.so")
JLLWrappers.@declare_library_product(libnwtclibs, "libnwtclibs.so")
JLLWrappers.@declare_library_product(libopenfast_postlib, "libopenfast_postlib.so")
JLLWrappers.@declare_library_product(libopenfast_prelib, "libopenfast_prelib.so")
JLLWrappers.@declare_library_product(libopenfastlib, "libopenfastlib.so")
JLLWrappers.@declare_library_product(libopenfoamtypeslib, "libopenfoamtypeslib.so")
JLLWrappers.@declare_library_product(liborcaflexlib, "liborcaflexlib.so")
JLLWrappers.@declare_library_product(libscdataexlib, "libscdataexlib.so")
JLLWrappers.@declare_library_product(libscdataexttypeslib, "libscdataextypeslib.so")
JLLWrappers.@declare_library_product(libscfastlib, "libscfastlib.so")
JLLWrappers.@declare_library_product(libsctypeslib, "libsctypeslib.so")
JLLWrappers.@declare_library_product(libservodynlib, "libservodynlib.so")
JLLWrappers.@declare_library_product(libsubdynlib, "libsubdynlib.so")
JLLWrappers.@declare_library_product(libuaaerolib, "libuaaerolib.so")
JLLWrappers.@declare_library_product(libversioninfolib, "libversioninfolib.so")
JLLWrappers.@declare_library_product(libwdlib, "libwdlib.so")
JLLWrappers.@declare_executable_product(aerodyn_driver)
JLLWrappers.@declare_executable_product(beamdyn_driver)
JLLWrappers.@declare_executable_product(dwm_driver_wind_farm)
JLLWrappers.@declare_executable_product(feam_driver)
JLLWrappers.@declare_executable_product(hydrodyn_driver)
JLLWrappers.@declare_executable_product(inflowwind_driver)
JLLWrappers.@declare_executable_product(moordyn_driver)
JLLWrappers.@declare_executable_product(openfast)
JLLWrappers.@declare_executable_product(orca_driver)
JLLWrappers.@declare_executable_product(servodyn_driver)
JLLWrappers.@declare_executable_product(subdyn_driver)
JLLWrappers.@declare_executable_product(turbsim)
JLLWrappers.@declare_executable_product(unsteadyaero_driver)
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, LAPACK_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libaeroacoustics,
        "lib/libaeroacoustics.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libaerodyn14lib,
        "lib/libaerodyn14lib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libaerodynlib,
        "lib/libaerodynlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libafinfolib,
        "lib/libafinfolib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libawaelib,
        "lib/libawaelib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbeamdynlib,
        "lib/libbeamdynlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libelastodynlib,
        "lib/libelastodynlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libextptfm_mckflib,
        "lib/libextptfm_mckflib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfeamlib,
        "lib/libfeamlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfoamfastlib,
        "lib/libfoamfastlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfvwlib,
        "lib/libfvwlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libhydrodynlib,
        "lib/libhydrodynlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicedynlib,
        "lib/libicedynlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicefloelib,
        "lib/libicefloelib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libifwlib,
        "lib/libifwlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmapcpplib,
        "lib/libmapcpplib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmaplib,
        "lib/libmaplib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmoordynlib,
        "lib/libmoordynlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnwtclibs,
        "lib/libnwtclibs.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenfast_postlib,
        "lib/libopenfast_postlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenfast_prelib,
        "lib/libopenfast_prelib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenfastlib,
        "lib/libopenfastlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenfoamtypeslib,
        "lib/libopenfoamtypeslib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liborcaflexlib,
        "lib/liborcaflexlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscdataexlib,
        "lib/libscdataexlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscdataexttypeslib,
        "lib/libscdataextypeslib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscfastlib,
        "lib/libscfastlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsctypeslib,
        "lib/libsctypeslib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libservodynlib,
        "lib/libservodynlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsubdynlib,
        "lib/libsubdynlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libuaaerolib,
        "lib/libuaaerolib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libversioninfolib,
        "lib/libversioninfolib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwdlib,
        "lib/libwdlib.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        aerodyn_driver,
        "bin/aerodyn_driver",
    )

    JLLWrappers.@init_executable_product(
        beamdyn_driver,
        "bin/beamdyn_driver",
    )

    JLLWrappers.@init_executable_product(
        dwm_driver_wind_farm,
        "bin/dwm_driver_wind_farm",
    )

    JLLWrappers.@init_executable_product(
        feam_driver,
        "bin/feam_driver",
    )

    JLLWrappers.@init_executable_product(
        hydrodyn_driver,
        "bin/hydrodyn_driver",
    )

    JLLWrappers.@init_executable_product(
        inflowwind_driver,
        "bin/inflowwind_driver",
    )

    JLLWrappers.@init_executable_product(
        moordyn_driver,
        "bin/moordyn_driver",
    )

    JLLWrappers.@init_executable_product(
        openfast,
        "bin/openfast",
    )

    JLLWrappers.@init_executable_product(
        orca_driver,
        "bin/orca_driver",
    )

    JLLWrappers.@init_executable_product(
        servodyn_driver,
        "bin/servodyn_driver",
    )

    JLLWrappers.@init_executable_product(
        subdyn_driver,
        "bin/subdyn_driver",
    )

    JLLWrappers.@init_executable_product(
        turbsim,
        "bin/turbsim",
    )

    JLLWrappers.@init_executable_product(
        unsteadyaero_driver,
        "bin/unsteadyaero_driver",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text: # TODO: gradients


text:
"""
intck(x,y) computes the difference between two date objects and return number of day. Any of x and y can be missing value.
"""
intck(x::Missing,y::Missing) = missing
intck(x::Missing,y) = missing
intck(x,y::Missing) = missing
intck(x::Date,y::Date) = (x-y).value

"""
rescale(x,minx,maxx,minval,maxval) rescales x to run from minval and maxval, given x originaly runs from minx to maxx.
"""
function rescale(x,minx,maxx,minval,maxval)
    -(-maxx*minval+minx*maxval)/(maxx-minx)+(-minval+maxval)*x/(maxx-minx)
end
rescale(x::Missing,minx,maxx,minval,maxval) = missing
rescale(x::Vector,minx,maxx,minval,maxval) = rescale.(x,minx,maxx,minval,maxval)
rescale(x,minx,maxx) = rescale(x,minx,maxx,0.0,1.0)

"""
stdze(x) standardizes an array. It return missing for missing data points.
"""
function stdze(x)
    all(ismissing,x) && return x
    meandata = mean(skipmissing(x))
    vardata = var(skipmissing(x))
    (x .- meandata) ./ sqrt(vardata)
end

"""
lag(x,k) Creates a lag-k of the provided array x. The output will be an array the
same size as x (the input array), and the its type will be Union{Missing, T} where T is the type of input.
"""
function lag(x,k)
    res=zeros(Union{typeof(x[1]),Missing},length(x))
    for i in 1:k
        res[i] = missing
    end
    for i in k+1:length(x)
        res[i] = x[i-k]
    end
    res
end

lag(x)=lag(x,1)

"""
lead(x,k) Creates a lead-k of the provided array x. The output will be an array the
same size as x (the input array), and the its type will be Union{Missing, T} where T is the type of input.
"""
function lead(x,k)
    res=zeros(Union{typeof(x[1]),Missing},length(x))
    for i in 1:length(x)-k
        res[i] = x[i+k]
    end
    for i in length(x)-k+1:length(x)
        res[i] = missing
    end
    res
end
lead(x)=lead(x,1)

"""
dttodate(x) converts SAS or STATA dates (which is the number of days after 01-01-1960) to a Julia Date object.
dttodate(DataFrame,cols) converts the given columns to Date object.
"""
dttodate(x::Missing) = missing
dttodate(x) = Date(1960,1,1)+Day(x)
dttodate(x::Date) = x
function dttodate!(df::DataFrame,cols)
    for i in cols
        df[!,i] = dttodate.(df[!,i])
    end
end

# modifying some Base functions to suit for working with data with missing values

Base.maximum(x::AbstractArray{Missing,1})=missing
function Base.maximum(x::AbstractArray{Union{T,Missing},1}) where {T <: Number}
    res=typemin(T)
    cnt_miss=0
    for i in 1:length(x)
        if !ismissing(x[i])
            if x[i]>res
                res=x[i]
            end
        else
            cnt_miss+=1
        end
    end
    if cnt_miss==length(x)
        res=missing
    end
    res
end
Base.maximum(x::AbstractArray{Union{Date,Missing},1})=maximum(skipmissing(x))

Base.minimum(x::AbstractArray{Missing,1})=missing
function Base.minimum(x::AbstractArray{Union{T,Missing},1}) where {T <: Number}
    res=typemax(T)
    cnt_miss=0
    for i in 1:length(x)
        if !ismissing(x[i])
            if res>x[i]
                res=x[i]
            end
        else
            cnt_miss+=1
        end
    end
    if cnt_miss==length(x)
        res=missing
    end
    res
end
Base.minimum(x::AbstractArray{Union{Date,Missing},1})=minimum(skipmissing(x))

Base.sum(x::AbstractArray{Missing,1})=missing
function Base.sum(x::AbstractArray{Union{T,Missing},1}) where {T <: Number}
    res=zero(T)
    cnt_miss=0
    for i in 1:length(x)
        if !ismissing(x[i])
            res+=x[i]
        else
            cnt_miss+=1
        end
    end
    if cnt_miss==length(x)
        res=missing
    end
    res
end

Statistics.mean(x::AbstractArray{Missing,1})=missing
function Statistics.mean(x::AbstractArray{Union{T,Missing},1}) where {T <: Number}
    res=0.0
    cnt_nonmiss=0
    for i in 1:length(x)
        if !ismissing(x[i])
            res+=x[i]
            cnt_nonmiss+=1
        end
    end
    if cnt_nonmiss>0
        res=res/cnt_nonmiss
    else
        res=missing
    end
    res
end

Statistics.var(x::AbstractArray{Missing,1},df=true)=missing
function Statistics.var(x::AbstractArray{Union{T,Missing},1},df=true) where {T <: Number}
    res=0.0
    ss=0.0
    sval=0.0
    cnt_nonmiss=0
    for i in 1:length(x)
        if !ismissing(x[i])
            sval+=x[i]
            ss+=x[i] * x[i]
            cnt_nonmiss+=1
        end
    end
    if cnt_nonmiss>1
        res=ss/cnt_nonmiss - (sval/cnt_nonmiss)*(sval/cnt_nonmiss)
        if df
            res=(cnt_nonmiss/(cnt_nonmiss-1))*res
        end
    elseif cnt_nonmiss==1
        res=0.0
    else
        res=missing
    end
    res
end
Statistics.std(x::AbstractArray{Missing,1},df=true)=missing
function Statistics.std(x::AbstractArray{Union{T,Missing},1},df=true) where {T <: Number}
    sqrt(var(x,df))
end

function Statistics.median(v::AbstractArray{T,1}) where T
    isempty(v) && throw(ArgumentError("median of an empty array is undefined, $(repr(v))"))
    eltype(v)>:Missing && all(ismissing, v) && return missing
    (eltype(v)<:AbstractFloat || eltype(v)>:AbstractFloat) && any(isnan, v) && return convert(eltype(v), NaN)
    any(ismissing,v) ? v2=collect(skipmissing(v)) : v2=v
    inds = axes(v2, 1)
    n = length(inds)
    mid = div(first(inds)+last(inds),2)
    if isodd(n)
        return middle(partialsort!(v2,mid))
    else
        m = partialsort!(v2, mid:mid+1)
        return middle(m[1], m[2])
    end
end

function Statistics.quantile(x::Array{T,1},v) where T
    all(ismissing,x) && return missing
    quantile(skipmissing(x),v)
end


text:
@testset "context" begin

ctx = CuCurrentContext()
dev = CuCurrentDevice()

let ctx2 = CuContext(dev)
    @test ctx2 == CuCurrentContext()    # ctor implicitly pushes
    activate(ctx)
    @test ctx == CuCurrentContext()

    @test CuDevice(ctx2) == dev

    CUDA.unsafe_destroy!(ctx2)
end

let global_ctx2 = nothing
    CuContext(dev) do ctx2
        @test ctx2 == CuCurrentContext()
        @test ctx != ctx2
        global_ctx2 = ctx2
    end
    @test !CUDA.isvalid(global_ctx2)
    @test ctx == CuCurrentContext()

    @test CuDevice(ctx) == dev
    @test CuCurrentDevice() == dev
    device_synchronize()
end

end

# FIXME: these trample over our globally-managed context

# @testset "primary context" begin

# pctx = CuPrimaryContext(device())

# @test !isactive(pctx)
# unsafe_reset!(pctx)
# @test !isactive(pctx)

# @test flags(pctx) == 0
# setflags!(pctx, CUDA.CTX_SCHED_BLOCKING_SYNC)
# @test flags(pctx) == CUDA.CTX_SCHED_BLOCKING_SYNC

# let global_ctx = nothing
#     CuContext(pctx) do ctx
#         @test CUDA.isvalid(ctx)
#         @test isactive(pctx)
#         global_ctx = ctx
#     end
#     @test !isactive(pctx)
#     @test !CUDA.isvalid(global_ctx)
# end

# CuContext(pctx) do ctx
#     @test CUDA.isvalid(ctx)
#     @test isactive(pctx)

#     unsafe_reset!(pctx)

#     @test !isactive(pctx)
#     @test !CUDA.isvalid(ctx)
# end

# let
#     @test !isactive(pctx)

#     ctx1 = CuContext(pctx)
#     @test isactive(pctx)
#     @test CUDA.isvalid(ctx1)

#     unsafe_reset!(pctx)
#     @test !isactive(pctx)
#     @test !CUDA.isvalid(ctx1)
#     CUDA.valid_contexts

#     ctx2 = CuContext(pctx)
#     @test isactive(pctx)
#     @test !CUDA.isvalid(ctx1)
#     @test CUDA.isvalid(ctx2)

#     unsafe_reset!(pctx)
# end

# end

@testset "cache config" begin

config = cache_config()

cache_config!(CUDA.FUNC_CACHE_PREFER_L1)
@test cache_config() == CUDA.FUNC_CACHE_PREFER_L1

cache_config!(config)

end

@testset "shmem config" begin

config = shmem_config()

shmem_config!(CUDA.SHARED_MEM_CONFIG_FOUR_BYTE_BANK_SIZE)
@test shmem_config() == CUDA.SHARED_MEM_CONFIG_FOUR_BYTE_BANK_SIZE

shmem_config!(config)

end

@testset "limits" begin

lim = limit(CUDA.LIMIT_DEV_RUNTIME_SYNC_DEPTH)

lim += 1
limit!(CUDA.LIMIT_DEV_RUNTIME_SYNC_DEPTH, lim)
@test lim == limit(CUDA.LIMIT_DEV_RUNTIME_SYNC_DEPTH)

limit!(CUDA.LIMIT_DEV_RUNTIME_SYNC_DEPTH, lim)

end


text:
using DelimitedFiles
using Printf
using Plots
using Statistics

"""
    imtz(t; clims=(-0.5,0.5), cmap=:pu_or, title="", anom=false, L=2*pi)

Plots a vertical rbsolve section at step `t`.

# Optional keywords
- `clims::Tuple{Real,Real}`: limits for the colormap
- `cmap::Symbol`: the colormap to use
- `title::String`: title of the plot
- `anom::Bool`: compute horizontal anomaly
- `L::Real`: X domain size
"""
function imtz(t; clims=(-0.5,0.5), cmap=:pu_or, title="", anom=false, L=2*pi)
    coord=readdlm("coord")
    (nz,nc)=size(coord)
    nz=convert(Int64,(nz-1)/2)+1
    coord=coord[1:nz,4]
    nz=nz-2
    filename=@sprintf("t%07d.zzz", t)
    tt=readdlm(filename)
    n=length(tt); n=n÷nz
    zz=coord[2:(nz+1)]
    xx=(0:(n-1))/n*L
    tt=reshape(tt, n, nz)
    if anom
       tt1=tt.-ones(n,1)*mean(tt,dims=1)
    else
       tt1=tt
    end
    plot(xx, zz.+0.5, tt1', st=:heatmap, c=cmap, clims=clims)
    xlabel!("x")
    ylabel!("z")
    title!(title)
end

max_stars_repo_name: jhardenberg/rbsolve

max_issues_repo_name: jhardenberg/rbsolve

max_forks_repo_name: jhardenberg/rbsolve


text:
using Test

include("alphametics.jl")

@testset "puzzle with three letters" begin
    @test solve("I + BB == ILL") == Dict(
        'I' => 1,
        'B' => 9,
        'L' => 0
    )
end

@testset "solution must have unique value for each letter" begin
    @test solve("A == B") === nothing
end

@testset "leading zero solution is invalid" begin
    @test solve("ACA + DD == BD") === nothing
end

@testset "puzzle with two digits final carry" begin
    @test solve("A + A + A + A + A + A + A + A + A + A + A + B == BCC") == Dict(
        'A' => 9,
        'B' => 1,
        'C' => 0
    )
end

@testset "puzzle with four letters" begin
    @test solve("AS + A == MOM") == Dict(
        'A' => 9,
        'S' => 2,
        'M' => 1,
        'O' => 0
    )
end

@testset "puzzle with six letters" begin
    @test solve("NO + NO + TOO == LATE") == Dict(
        'N' => 7,
        'O' => 4,
        'T' => 9,
        'L' => 1,
        'A' => 0,
        'E' => 2
    )
end

@testset "puzzle with seven letters" begin
    @test solve("HE + SEES + THE == LIGHT") == Dict(
        'E' => 4,
        'G' => 2,
        'H' => 5,
        'I' => 0,
        'L' => 1,
        'S' => 9,
        'T' => 7
    )
end

@testset "puzzle with eight letters" begin
    @test solve("SEND + MORE == MONEY") == Dict(
        'S' => 9,
        'E' => 5,
        'N' => 6,
        'D' => 7,
        'M' => 1,
        'O' => 0,
        'R' => 8,
        'Y' => 2
    )
end

@testset "puzzle with ten letters" begin
    @test solve("AND + A + STRONG + OFFENSE + AS + A + GOOD == DEFENSE") == Dict(
        'A' => 5,
        'D' => 3,
        'E' => 4,
        'F' => 7,
        'G' => 8,
        'N' => 0,
        'O' => 2,
        'R' => 1,
        'S' => 6,
        'T' => 9
    )
end

@testset "puzzle with ten letters and 199 addends" begin
    @test solve("THIS + A + FIRE + THEREFORE + FOR + ALL + HISTORIES + I + TELL + A + TALE + THAT + FALSIFIES + ITS + TITLE + TIS + A + LIE + THE + TALE + OF + THE + LAST + FIRE + HORSES + LATE + AFTER + THE + FIRST + FATHERS + FORESEE + THE + HORRORS + THE + LAST + FREE + TROLL + TERRIFIES + THE + HORSES + OF + FIRE + THE + TROLL + RESTS + AT + THE + HOLE + OF + LOSSES + IT + IS + THERE + THAT + SHE + STORES + ROLES + OF + LEATHERS + AFTER + SHE + SATISFIES + HER + HATE + OFF + THOSE + FEARS + A + TASTE + RISES + AS + SHE + HEARS + THE + LEAST + FAR + HORSE + THOSE + FAST + HORSES + THAT + FIRST + HEAR + THE + TROLL + FLEE + OFF + TO + THE + FOREST + THE + HORSES + THAT + ALERTS + RAISE + THE + STARES + OF + THE + OTHERS + AS + THE + TROLL + ASSAILS + AT + THE + TOTAL + SHIFT + HER + TEETH + TEAR + HOOF + OFF + TORSO + AS + THE + LAST + HORSE + FORFEITS + ITS + LIFE + THE + FIRST + FATHERS + HEAR + OF + THE + HORRORS + THEIR + FEARS + THAT + THE + FIRES + FOR + THEIR + FEASTS + ARREST + AS + THE + FIRST + FATHERS + RESETTLE + THE + LAST + OF + THE + FIRE + HORSES + THE + LAST + TROLL + HARASSES + THE + FOREST + HEART + FREE + AT + LAST + OF + THE + LAST + TROLL + ALL + OFFER + THEIR + FIRE + HEAT + TO + THE + ASSISTERS + FAR + OFF + THE + TROLL + FASTS + ITS + LIFE + SHORTER + AS + STARS + RISE + THE + HORSES + REST + SAFE + AFTER + ALL + SHARE + HOT + FISH + AS + THEIR + AFFILIATES + TAILOR + A + ROOFS + FOR + THEIR + SAFE == FORTRESSES") == Dict(
        'A' => 1,
        'E' => 0,
        'F' => 5,
        'H' => 8,
        'I' => 7,
        'L' => 2,
        'O' => 6,
        'R' => 3,
        'S' => 4,
        'T' => 9
    )
end

max_stars_repo_name: stepinski/julia

max_issues_repo_name: stepinski/julia

max_forks_repo_name: stepinski/julia


text:
#==============================================================================#
#
#          Structs defining buffers used for in-place computations
#
#==============================================================================#

"""
    AbstractBuffer{T} <: AbstractArray{T,1}

Types inheriting from `AbstractBuffer` define buffers used for various in-place
computations. Each subtype MUST:
- have a field `data` which is a one-dimensional vector with data
- provide a constructor of the form: NAME{eltype}(dim1, dim2, ...)
and it CAN:
- have multiple fields that provide views or reshaped views to `data`
"""
abstract type AbstractBuffer{T} <: AbstractArray{T,1} end

const _AB{T} = AbstractBuffer{T} where T # an alias

Base.size(var::_AB) = size(var.data)
Base.eltype(var::_AB{T}) where T = T

Base.getindex(var::_AB, i::Int) = var.data[i]
Base.getindex(var::_AB, I::Vararg{Int,N}) where N = var.data[I...]
Base.getindex(var::_AB, ::Colon) = var.data[:]
Base.getindex(var::_AB, kr::AbstractRange) = var.data[kr]

Base.setindex!(var::_AB, v, i::Int) = (var.data[i] = v)
Base.setindex!(var::_AB, v, I::Vararg{Int,N}) where N = (var.data[I...] = v)
Base.setindex!(var::_AB, v, ::Colon) = (var.data[:] .= v)
Base.setindex!(var::_AB, v, kr::AbstractRange) = (var.data[kr] .= v)

"""
    remove_curly(::Type{K}) where K

Utility function that removes all type-specifiers listed in the curly brackets.

# Examples
```julia-repl
julia> remove_curly(Array{Float64,1})
Array
```
"""
@generated function remove_curly(::Type{K}) where K
    name_without_curly = Meta.parse(string(K)).args[1]
    :( $name_without_curly )
end

"""
    get_curly(::Type{K}) where K

Utility function that returns a tuple with all type-specifiers listed in the
curly brackets.

# Examples
```julia-repl
julia> remove_curly(Array{Float64,1})
(:Float64, 1)
```
"""
@generated function get_curly(::Type{K}) where K
    curly_arguments = Meta.parse(string(K)).args[2:end]
    :( Tuple($curly_arguments) )
end

"""
    dimensions(::Type{K}) where K<:_AB

Utility function for structs inheriting from AbstractBuffer. Returns the
dimensions required for constructing the internal views of the data.
"""
dimensions(::Type{K}) where K<:_AB = get_curly(K)[2].args

@generated function Base.similar(::K) where K<:_AB{T} where T
    pure_name = Meta.parse(string((remove_curly(K))))
    constructor = Expr(
        :call,
        Expr(
            :curly,
            :($pure_name),
            T,
        )
    )
    append!(constructor.args, [dimensions(K)...])
    :(
        $constructor
    )
end

@generated function Base.similar(::K, ::Type{T}) where {K<:_AB,T}
    pure_name = Meta.parse(string((remove_curly(K))))
    constructor = Expr(
        :call,
        Expr(
            :curly,
            :($pure_name),
            T,
        )
    )
    append!(constructor.args, [dimensions(K)...])
    :(
        $constructor
    )
end

"""
    struct StandardEulerBuffer{T,D,Tb,Tσ,Tdw} <: AbstractBuffer{T}
        data::Vector{T}
        b::Tb
        σ::Tσ
        dW::Tdw
    end

Standard buffer for the Euler-Maruyama simulations. The data is stored in
`data` and `b`, `σ`, `dW` provide appropriately reshaped views to the
corresponding segments of `data`. `T` is the DataType of each data element
and `D` are dimensions needed to create the views.
"""
struct StandardEulerBuffer{T,D,Tb,Tσ,Tdw} <: AbstractBuffer{T}
    data::Vector{T}
    b::Tb
    σ::Tσ
    dW::Tdw

    function StandardEulerBuffer{T}(dim_process::Number, dim_wiener::Number) where T
        D, m = dim_process, dim_wiener
        data = zeros(T, D + D*m + m)

        b = view(data, 1:D)
        σ = reshape( view(data, (D+1):(D+D*m)), (D,m) )
        dW = view(data, (D+D*m+1):(D+D*m+m))
        new{T,(D,m),typeof(b),typeof(σ),typeof(dW)}(data, b, σ, dW)
    end
end

"""
    struct LinearDiffBuffer{T,D,Tb,TB,Tσ,Tdw} <: AbstractBuffer{T}
        data::Vector{T}
        b::Tb
        B::TB
        σ::Tσ
        dW::Tdw
    end

A buffer for Euler-Maruyama simulations of linear diffusions. Almost the same
as `StandardEulerBuffer`, but contains additional space for an intermediate
construction of a matrix `B` (and a corrsponding view).
"""
struct LinearDiffBuffer{T,D,Tb,TB,Tσ,Tdw} <: AbstractBuffer{T}
    data::Vector{T}
    b::Tb
    B::TB
    σ::Tσ
    dW::Tdw

    function LinearDiffBuffer{T}(dim_process::Number, dim_wiener::Number) where T
        D, m = dim_process, dim_wiener
        data = zeros(T, D + D*D + D*m + m)

        b = view(data, 1:D)
        B = reshape( view(data, (D+1):(D+D*D)), (D,D) )
        σ = reshape( view(data, (D+D*D+1):(D+D*D+D*m)), (D,m) )
        dW = view(data, (D+D*D+D*m+1):(D+D*D+D*m+m))
        new{T,(D,m),typeof(b),typeof(B),typeof(σ),typeof(dW)}(data, b, B, σ, dW)
    end
end


text:
module lcmtesttypes

export lcm_test_type_1, lcm_test_type_2, lcm_test_type_3, polynomial_t, polynomial_matrix_t
export hard_coded_example

using LCMCore, StaticArrays
using Random: randstring

mutable struct lcm_test_type_1 <: LCMType
    a::Int16
    blength::Int32
    b::Vector{Int64}
    c::SVector{3, Int32}
end

@lcmtypesetup(lcm_test_type_1,
    b => (blength, )
)

function Base.:(==)(x::lcm_test_type_1, y::lcm_test_type_1)
    x.a == y.a || return false
    x.blength == y.blength || return false
    x.b == y.b || return false
    x.c == y.c || return false
    true
end

function Base.rand(::Type{lcm_test_type_1})
    blength = rand(Int32(0) : Int32(10))
    lcm_test_type_1(rand(Int16), blength, rand(Int64, blength), rand(SVector{3, Int32}))
end

function hard_coded_example(::Type{lcm_test_type_1})
    ret = lcm_test_type_1()
    ret.a = 1
    ret.blength = 4
    ret.b = [4, 5, -6, 7]
    ret.c = [10, -1, 5]
    ret
end

mutable struct lcm_test_type_2 <: LCMType
    dlength::Int32
    f_inner_length::Int32
    a::Bool
    b::UInt8
    c::lcm_test_type_1
    d::Vector{lcm_test_type_1}
    e::SVector{3, lcm_test_type_1}
    f::Matrix{Int64}
end

@lcmtypesetup(lcm_test_type_2,
    d => (dlength, ),
    f => (3, f_inner_length)
)

function Base.:(==)(x::lcm_test_type_2, y::lcm_test_type_2)
    x.dlength == y.dlength || return false
    x.f_inner_length == y.f_inner_length || return false
    x.a == y.a || return false
    x.b == y.b || return false
    x.c == y.c || return false
    x.d == y.d || return false
    x.e == y.e || return false
    x.f == y.f || return false
    true
end

function Base.rand(::Type{lcm_test_type_2})
    dlength = rand(Int32(0) : Int32(10))
    f_inner_length = rand(Int32(0) : Int32(10))
    a = rand(Bool)
    b = rand(UInt8)
    c = rand(lcm_test_type_1)
    d = [rand(lcm_test_type_1) for i = 1 : dlength]
    e = SVector{3}([rand(lcm_test_type_1) for i = 1 : 3])
    f = rand(Int64, 3, f_inner_length)
    lcm_test_type_2(dlength, f_inner_length, a, b, c, d, e, f)
end

function hard_coded_example(::Type{lcm_test_type_2})
    ret = lcm_test_type_2()
    ret.dlength = 2
    ret.f_inner_length = 2
    ret.a = false
    ret.b = 6
    ret.c = hard_coded_example(lcm_test_type_1)
    ret.d = [hard_coded_example(lcm_test_type_1), hard_coded_example(lcm_test_type_1)]
    ret.d[1].a = 2
    ret.e = [hard_coded_example(lcm_test_type_1), hard_coded_example(lcm_test_type_1), hard_coded_example(lcm_test_type_1)]
    ret.e[3].c = [5, 3, 8]
    ret.f = [1 2; 3 4; 5 7]
    ret
end

mutable struct lcm_test_type_3 <: LCMType
    a::String
    blength::Int32
    b::Vector{String}
    c::SVector{2, String}
    d::Float64
end

@lcmtypesetup(lcm_test_type_3,
    b => (blength, )
)

function Base.:(==)(x::lcm_test_type_3, y::lcm_test_type_3)
    x.a == y.a || return false
    x.blength == y.blength || return false
    x.b == y.b || return false
    x.c == y.c || return false
    x.d == y.d || return false
    true
end

function Base.rand(::Type{lcm_test_type_3})
    a = randstring()
    blength = rand(Int32(0) : Int32(10))
    b = [randstring() for i = 1 : blength]
    c = SVector{2}([randstring() for i = 1 : 2])
    d = rand()
    lcm_test_type_3(a, blength, b, c, d)
end

function hard_coded_example(::Type{lcm_test_type_3})
    ret = lcm_test_type_3()
    ret.a = "abcd"
    ret.blength = 3
    ret.b = ["x*4f", "4^G32", "4"]
    ret.c = ["xyz", "wxy"]
    ret.d = 2.5
    ret
end

mutable struct polynomial_t <: LCMType
    timestamp::Int64
    num_coefficients::Int32
    coefficients::Vector{Float64}
end

@lcmtypesetup(polynomial_t,
    coefficients => (num_coefficients, )
)

function Base.:(==)(x::polynomial_t, y::polynomial_t)
    x.timestamp == y.timestamp || return false
    x.num_coefficients == y.num_coefficients || return false
    x.coefficients == y.coefficients || return false
    true
end

function Base.rand(::Type{polynomial_t})
    timestamp = rand(Int64)
    num_coefficients = rand(Int32(0) : Int32(10))
    coefficients = rand(num_coefficients)
    polynomial_t(timestamp, num_coefficients, coefficients)
end

function hard_coded_example(::Type{polynomial_t})
    polynomial_t(1234, 4, [5.0, 0.0, 1.0, 3.0])
end

mutable struct polynomial_matrix_t <: LCMType
    timestamp::Int64
    rows::Int32
    cols::Int32
    polynomials::Matrix{polynomial_t}
end

@lcmtypesetup(polynomial_matrix_t,
    polynomials => (rows, cols)
)

function Base.:(==)(x::polynomial_matrix_t, y::polynomial_matrix_t)
    x.timestamp == y.timestamp || return false
    x.rows == y.rows || return false
    x.cols == y.cols || return false
    x.polynomials == y.polynomials || return false
    true
end

function Base.rand(::Type{polynomial_matrix_t})
    timestamp = rand(Int64)
    rows = rand(Int32(0) : Int32(10))
    cols = rand(Int32(0) : Int32(10))
    polynomials = [rand(polynomial_t) for row = 1 : rows, col = 1 : cols]
    polynomial_matrix_t(timestamp, rows, cols, polynomials)
end

function hard_coded_example(::Type{polynomial_matrix_t})
    polynomials = Matrix{polynomial_t}(undef, 2, 1)
    polynomials[1] = polynomial_t(1234, 4, [5.0, 0.0, 1.0, 3.0])
    polynomials[2] = polynomial_t(1234, 4, [10.0, 0.0, 2.0, 6.0])
    polynomial_matrix_t(1234, 2, 1, polynomials)
end

end # module


text:
This is atest

# if x > y
#     return bob
# end

# for x in 1:10
#     println(x)
# end

max_stars_repo_name: haberdashPI/vscode-selection-utilities

max_issues_repo_name: haberdashPI/vscode-selection-utilities

max_forks_repo_name: haberdashPI/vscode-selection-utilities


text:
using Compat
@compat using Seismic
using Base.Test

d,e = SeisLinearEvents()
ddec = SeisDecimate(d,perc=50)

dmwni = SeisMWNI(ddec,dt=0.004,fmax=100,Niter_external=10,Niter_internal=3)

# test that quality factor is greater than 10 Decibels
quality_factor = 10*log10(norm(d[:],2)/norm(dmwni[:]-d[:],2))
println("Quality factor = ",quality_factor)
@test quality_factor > 5


text:
"""
```
t = find_threshold(Intermodes(), histogram, edges; maxiter = 8000)
```

Under the assumption that the histogram is bimodal the histogram is smoothed
using a length-3 mean filter until two modes remain. The threshold is then set
to the average value of the two modes.

# Output

Returns a real number `t` in `edges`. The `edges` parameter represents an
`AbstractRange` which specifies the intervals associated with the histogram bins.

# Details

If after `maxiter` iterations the smoothed histogram is still not bimodal then
the algorithm will fall back to using the `UnimodalRosin` method to select
a threshold.

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

## `maxiter`

An `Int` that specifies the maximum number of smoothing iterations. If left
unspecified a default value of 8000 is used.

# Example

Compute the threshold for the "cameraman" image in the `TestImages` package.

```julia

using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img,256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(Intermodes(), counts[1:end], edges)
```

# Reference

1. C. A. Glasbey, “An Analysis of Histogram-Based Thresholding Algorithms,” CVGIP: Graphical Models and Image Processing, vol. 55, no. 6, pp. 532–537, Nov. 1993. [doi:10.1006/cgip.1993.1040](https://doi.org/10.1006/cgip.1993.1040)
"""
function find_threshold(algorithm::Intermodes, histogram::AbstractArray, edges::AbstractRange)
    bimodal_histogram = smooth_histogram(histogram, 1000)
    indices = find_maxima_indices(bimodal_histogram)
    if length(indices) != 2
        @warn "Failed to find two modes. Falling back to `UnimodalRosin` method."
        return find_threshold(UnimodalRosin(), bimodal_histogram, edges)
    else
        index = round(Int, (indices[1] + indices[2]) / 2)
        return edges[index]
    end
end


text:
using Test
using DataFrames

# see https://discourse.julialang.org/t/writing-tests-in-vs-code-workflow-autocomplete-and-tooltips/57488
# see https://github.com/julia-vscode/julia-vscode/issues/800
if isdefined(@__MODULE__, :LanguageServer)
    # invoked by VS Code
    include("../src/OilData.jl")
    using .OilData
else
    # invoked during test
    using OilData
end

# Avoid include common.jl multiple times
if !@isdefined local_6cfddc95df03456583406d57fd5963ac
    local_6cfddc95df03456583406d57fd5963ac = true
end
if local_6cfddc95df03456583406d57fd5963ac
    include("common.jl")
end

function test_downsample()
    values = [1, 1.1, 2, 2.2, 2.5, 3, 4, 4.3, 5, 5.5, 6, 7, 8]
    @test downsample(values,2) == [1, 3, 5, 7, 8]
    @test downsample(values,3) == [1, 4, 7, 8]
end

function test_find_duplication()
    values = [1, 1.5, 2, 2, 2, 2, 3]
    dup_inds, dup_vals = find_duplication(values)
    @test dup_inds == [4,5,6]
    @test dup_vals == [2,2,2]
end

function test_find_integers()
    values = [0, 0.5, 1.0, 1.2, 2, 2.5, 3, 3.8, 4]
    inds,ints = find_integer_times(values)
    @test inds == [1,3,5,7,9]
    @test ints == [0, 1, 2, 3, 4]

    values = [0, 0.5, 2, 3, 3, 4]
    @test_throws ErrorException find_integer_times(values)
end

function test_find_indices()
    # indices are order by the values array,
    # not the filter
    values = [0, 0.5, 1.0, 1.2, 2, 2.5, 3, 3.8, 4]
    filter = [1, 3]
    inds = find_indices(filter, values)
    @test inds == [3,7]

    values = [0, 0.5, 1.0, 1.2, 2, 2.5, 3, 3.8, 4]
    filter = [3, 1]
    inds = find_indices(filter, values)
    @test inds == [3,7]

    # Reverse the order of values
    values = reverse([0, 0.5, 1.0, 1.2, 2, 2.5, 3, 3.8, 4])
    filter = [1, 4]
    inds = find_indices(filter, values)
    @test inds == [1,7]

    values = reverse([0, 0.5, 1.0, 1.2, 2, 2.5, 3, 3.8, 4])
    filter = [4, 1]
    inds = find_indices(filter, values)
    @test inds == [1,7]

    # Duplicated values
    values = [0, 1, 0.5, 2, 3, 3, 3.8, 4]
    filter = [1, 3]
    @test_throws ErrorException find_indices(filter, values)
end

function test_find_last_duplication()
    @test 1 == find_last_duplication(1, [1, 1.5, 2, 2, 2, 2, 3]) # 1
    @test 6 == find_last_duplication(3, [1, 1.5, 2, 2, 2, 2, 3]) # 2
    @test 7 == find_last_duplication(7, [1, 1.5, 2, 2, 2, 2, 3]) # 3
end

function test_repair_times()
    times = [1, 1.5, 2, 2, 2, 2, 3]
    @test [1, 1.5, 2, 2.25, 2.50, 2.75, 3] == repair_times(times)
end

function test_smooth_production()
    df = DataFrame( # Pair constructor
        "TIME" => [
            0,

            0.5,
            1.0,

            1.2,
            2,

            2.5,
            3,

            3.8,
            4,
        ],
        "WOPR" => [
            0,
            5,
            10,
            12,
            20,
            25,
            30,
            38,
            40,
        ],
    )

    rdf = smooth_production(df, "TIME", ["WOPR"])
    @test rdf[!,"TIME"] == [0, 1, 2, 3, 4]
    @test rdf[!,"WOPR"] == [0.0, 7.5, 18.4, 27.5, 38.4]
end

@testset "utility/smooth" begin
    test_downsample()
    test_find_duplication()
    test_find_integers()
    test_find_indices()
    test_find_last_duplication()
    test_repair_times()
    test_smooth_production()
end

nothing


text:
"""
    floatmod(x, y)

This is the equivalent of `x // y` for floating point.
"""
floatmod(x, y) = x - y * floor(x / y)

"""
    fibonacci_points(n)

Constructs lat-long points on a sphere from a Fibonacci sequence.

There are 2n + 1 points on the sphere. It returns a `(2,2n+1)` grid of lat-long.

These points are a good start on an evenly-spaced set of points on the
surface of a sphere. This algorithm is from González [^1]. The usual way
to use these points is to make them a little more even by adjustment
with a successive over-relaxation algorithm. That algorithm moves them
a little at each step due to a repulsive force among neighbors.
Then you make a Voronoi on the sphere and create triangles from this.
The Voronoi will be almost all hexagons, but there will be a few polygons
with 5 and 7 sides.

[^1]: González, Á. (2010). Measurement of areas on a sphere using Fibonacci
and latitude–longitude lattices. Mathematical Geosciences, 42(1), 49.
"""
function fibonacci_points(n)
    ϕ = (1 + sqrt(5)) / 2
    points = zeros(Float64, 2, 2n + 1)
    for i in -n:n
        points[1, i + n + 1] = asin(2i / (2n + 1)) * 180 / π
        lon = floatmod(i, ϕ) * 360 / ϕ
        if lon < -180
            points[2, i + n + 1] = lon + 360
        elseif lon >= 180
            points[2, i + n + 1] = lon - 360
        else
            points[2, i + n + 1] = lon
        end
    end
    points
end


text:
using PCL
using Cxx

PT = pcl.PointXYZRGB
kumamon = pcl.PointCloud{PT}("kumamon.pcd")
@show length(kumamon)

kumamon_zfiltered = pcl.PointCloud{PT}()
kumamon_xfiltered = pcl.PointCloud{PT}()
kumamon_yfiltered = pcl.PointCloud{PT}()

zpass = pcl.PassThrough{pcl.PointXYZRGB}()
pcl.setInputCloud(zpass, kumamon)
pcl.setFilterFieldName(zpass, "z")
pcl.setFilterLimits(zpass, 0.0, 1.5)
filter(zpass, kumamon_zfiltered)
@show length(kumamon_zfiltered)

xpass = pcl.PassThrough{pcl.PointXYZRGB}()
pcl.setInputCloud(xpass, kumamon_zfiltered)
pcl.setFilterFieldName(xpass, "x")
pcl.setFilterLimits(xpass, -0.4, 0.4)
filter(xpass, kumamon_xfiltered)
@show length(kumamon_xfiltered)

ypass = pcl.PassThrough{pcl.PointXYZRGB}()
pcl.setInputCloud(ypass, kumamon_xfiltered)
pcl.setFilterFieldName(ypass, "y")
pcl.setFilterLimits(ypass, 0, 0.23)
filter(ypass, kumamon_yfiltered)
@show length(kumamon_yfiltered)


text:
using Plots, Manifolds, ManifoldsBase, ManifoldDiffEq, Documenter, PyPlot
ENV["GKSwstype"] = "100"

makedocs(
    # for development, we disable prettyurls
    format = Documenter.HTML(prettyurls = false, assets = ["assets/favicon.ico"]),
    modules = [ManifoldDiffEq],
    authors = "Seth Axen, Mateusz Baran, Ronny Bergmann, and contributors.",
    sitename = "ManifoldDiffEq.jl",
    pages = [
        "Home" => "index.md",
        "Lie group action solvers" => "lie_group_solvers.md",
        "Frozen coefficients solvers" => "frozen_coefficients.md",
        "Error estimation" => "error_estimation.md",
        "Notation" => "notation.md",
        "Examples" => "examples.md",
    ],
)
deploydocs(
    repo = "github.com/JuliaManifolds/ManifoldDiffEq.jl.git",
    push_preview = true,
    devbranch = "main",
)


text:
abstract type BasicModel end
abstract type ModelResults end

### Mean Forecast Model
mutable struct MeanForecast <: BasicModel
    y::Vector{AbstractFloat}
    h::Integer
    level::Union{Nothing, Vector{AbstractFloat}}
    function MeanForecast(y = []; h = 5, level = [80, 95]) # Constructor
        if length(y) <= 0
            error("The input array is empty")
        end
        if h <= 0
            error("The h-value must be positive")
        end

        cleanPredInterval_(level)
        return new(y, h, level)
    end
end

function predict(model::MeanForecast)
    y = model.y
    h = model.h
    level = model.level
    data_size = length(y)

    # Compute mean, concatenate it onto the original data
    data_mean = mean(y)
    mean_vector = ones(h + data_size) .* data_mean
    fittedvalues = mean_vector
    residuals = y - fittedvalues[1:data_size] # compute residuals

    ### Coverage Probability Estimate, will calculate for 80 and 95 for now.
    c_80 = 1.28
    c_95 = 1.96

    K = 0 # number of parameters estimated in the forecasting method

    lower = zeros(length(level), h)
    upper = lower

    # Compute Prediction Intervals
    if h != 0
        # Calculate one-step prediction intervals
        one_step = sqrt(1 / (data_size - K) * (transpose(residuals) * residuals))
        σ = one_step
    end
    if h > 1
        # Calculate multi-step prediction intervals
        multi_step = fill(one_step * sqrt(1 + 1/data_size), h)
        multi_step[1] = one_step
        σ = multi_step
    end

    # Get the last h values
    yhat = fittedvalues[data_size + 1 : data_size + h]

    # Calculate prediction intervals, store in upper / lower matrices
    lower[1, :] = yhat - c_80 * σ
    lower[2, :] = yhat - c_95 * σ

    upper[1, :] = yhat + c_80 * σ
    upper[2, :] = yhat + c_95 * σ

    results = BasicMethodResults(model, fittedvalues, lower, upper)
    return results
end

function fit(model::MeanForecast)
    return model
end

### Naive Forecast Model
mutable struct NaiveForecast <: BasicModel
    y::Vector{AbstractFloat}
    h::Integer
    level::Union{Nothing, Vector{AbstractFloat}}
    function NaiveForecast(y = []; h = 5, level = [80, 95]) # Constructor
        if length(y) <= 0
            error("The input array is empty")
        end
        if h <= 0
            error("The h-value must be positive")
        end
        cleanPredInterval_(level)
        return new(y, h, level)
    end
end

function predict(model::NaiveForecast)
    y = model.y
    h = model.h
    level = model.level
    data_size = length(y)

    # Forecasts are equal to last value of observations
    lastValue = last(y)
    fitted_chunk = fill(lastValue, (h))
    fittedvalues = vcat(y, fitted_chunk)

    # Compute Prediction Intervals
    lower = zeros(length(level), h)
    upper = lower
    for i in 1:length(level)
        dist = Normal(0.0, 1.0)
        qq = quantile(dist, 0.5 * (1 + level[i]/100))
        lower[i, :] = fitted_chunk - qq * fitted_chunk
        upper[i, :] = fitted_chunk + qq * fitted_chunk
    end

    results = BasicMethodResults(model, fittedvalues, lower, upper)
    return results
end

function fit(model::NaiveForecast)
    return model
end

### Seasonal Naive Model
mutable struct SeasonalNaiveForecast <: BasicModel
    y::Vector{AbstractFloat}
    h::Integer
    m::Integer
    level::Union{Nothing, Vector{AbstractFloat}}
    function SeasonalNaiveForecast(y = []; h = 5, m = 12, level = [80, 95]) # Constructor
        if length(y) <= 0
            error("The input array is empty")
        end
        if h <= 0
            error("The h-value must be positive")
        end
        if m <= 0
            error("The m-value must be positive")
        end
        cleanPredInterval_(level)
        return new(y, h, m, level)
    end
end

function fit(model::SeasonalNaiveForecast)
    y = model.y
    h = model.h
    m = model.m
    level = model.level
    data_size = length(y)

    # Forecasts are equal to last value of observations
    k = trunc(Integer, (h-1)/m)
    fitted_chunk = y[data_size + 1 - m * (k+1) : data_size + h - m * (k+1)]
    fittedvalues = vcat(y, fitted_chunk)

    # Compute Prediction Intervals
    lower = zeros(length(level), h)
    upper = lower
    for i in 1:length(level)
        dist = Normal(0.0, 1.0)
        qq = quantile(dist, 0.5 * (1 + level[i]/100))
        lower[i, :] = fitted_chunk - qq * fitted_chunk
        upper[i, :] = fitted_chunk + qq * fitted_chunk
    end

    results = BasicMethodResults(model, fittedvalues, lower, upper)
    return results
end

# Ensures Prediction Intervals are numerical and on the interval (0, 100)
function cleanPredInterval_(level)
    for i in level
        if !isa(i, Number)
            print(isa(i, Number))
            error("Ensure that all level inputs are numeric")
        elseif i < 0 || i > 100
            error("Ensure that all confidence interval inputs are in the interval [0, 100]")
        end
    end
end

# Results Object - Common across all "Basic Method" models
mutable struct BasicMethodResults <: ModelResults # Stores results from fit(MeanForecast)
    model::BasicModel
    fittedvalues::Vector{Float64} # Vector of fitted values
    lower::Matrix{Float64} # Lower bound of prediction intervals
    upper::Matrix{Float64} # Upper bound of prediction intervals
    function BasicMethodResults(model, fittedvalues, lower, upper)
        return new(model, fittedvalues, lower, upper)
    end
end


text:
using LinearMaps
using IterativeSolvers
using ForwardDiff
using ProgressMeter

"""
Compute the independent-particle susceptibility. Will blow up for large systems.
Drop all non-diagonal terms with (f(εn)-f(εm))/(εn-εm) factor less than `droptol`.
"""
function compute_χ0(ham; droptol=0, temperature=ham.basis.model.temperature)
    # We're after χ0(r,r') such that δρ = ∫ χ0(r,r') δV(r') dr'
    # where (up to normalizations)
    # ρ = ∑_nk f(εnk - εF) |ψnk|^2
    # ∑_nk f(εnk - εF) = N_el
    # Everything is summed on k so we omit it for notational simplicity

    # We differentiate wrt a variation δV of the external potential
    # δρ = ∑_n (f'n δεn |ψn|^2 + 2Re fn ψn* δψn - f'n δεF |ψn|^2
    # with fn = f(εnk - εF), f'n = f'(εnk - εF)
    # δN_el = 0 = ∑_n f'n (δεn - δεF)

    # Now we use from first order perturbation theory
    # δεn = <ψn|δV|ψn>
    # δψn = ∑_{m != n} <ψm|δV|ψn> |ψm> / (εn-εm)

    # for δεF we get, with DOS = -∑_n f'n and LDOS = -∑_n f'n |ψn|^2
    # δεF = 1/DOS ∫ δV(r) LDOS(r)dr

    # for δρ we note ρnm = ψn* ψm, and we get
    # δρ = LDOS δεF + ∑_n f'n <ρn|δV> ρn + ∑_{n,m != n} 2Re fn ρnm <ρmn|δV> / (εn-εm)
    # δρ = LDOS δεF + ∑_n f'n <ρn|δV> ρn + ∑_{n,m != n} (fn-fm)/(εn-εm) ρnm <ρnm|δV>
    # The last two terms merge with the convention that (f(x)-f(x))/(x-x) = f'(x) into
    # δρ = LDOS δεF + ∑_{n,m} (fn-fm) ρnm <ρmn|δV> / (εn-εm)
    # Therefore the kernel is LDOS(r) LDOS(r') / DOS + ∑_{n,m} (fn-fm)/(εn-εm) ρnm(r) ρmn(r')
    basis = ham.basis
    model = basis.model
    fft_size = basis.fft_size
    @assert model.spin_polarization == :none
    filled_occ = filled_occupation(model)
    dVol = basis.model.unit_cell_volume / prod(basis.fft_size)

    EVs = [eigen(Hermitian(Array(Hk))) for Hk in ham.blocks]
    Es = [EV.values for EV in EVs]
    Vs = [EV.vectors for EV in EVs]
    occ, εF = find_occupation(basis, Es, temperature=temperature)

    χ0 = zeros(eltype(basis), prod(fft_size), prod(fft_size))
    for ik = 1:length(basis.kpoints)
        if length(basis.ksymops[ik]) != 1
            error("Kpoint symmetry not supported")
        end
        N = length(G_vectors(basis.kpoints[ik]))
        @assert N < 10_000
        E = Es[ik]
        V = Vs[ik]
        Vr = cat(G_to_r.(Ref(basis), Ref(basis.kpoints[ik]), eachcol(V))..., dims=4)
        Vr = reshape(Vr, prod(fft_size), N)
        @showprogress "Computing χ0 for kpoint $ik/$(length(basis.kpoints)) ..." for m = 1:N, n = 1:N
            enred = (E[n] - εF) / temperature
            @assert occ[ik][n] ≈ filled_occ * Smearing.occupation(model.smearing, enred)
            ddiff = Smearing.occupation_divided_difference
            ratio = filled_occ * ddiff(model.smearing, E[m], E[n], εF, temperature)
            (n != m) && (abs(ratio) < droptol) && continue
            # dVol because inner products have a dVol so that |f> becomes <dVol f|
            # can take the real part here because the nm term is complex conjugate of mn
            # TODO optimize this a bit... use symmetry nm, reduce allocs, etc.
            factor = basis.kweights[ik] * ratio * dVol
            @views χ0 .+= factor .* real(conj((Vr[:, m] .* Vr[:, m]'))
                                            .*   (Vr[:, n] .* Vr[:, n]'))
        end
    end

    # Add variation wrt εF
    if temperature > 0
        ldos = vec(LDOS(εF, basis, Es, Vs))
        dos = DOS(εF, basis, Es)
        χ0 .+= (ldos .* ldos') .* dVol ./ dos
    end

    χ0
end

"""
Returns the change in density δρ for a given δV. Drop all non-diagonal terms with
(f(εn)-f(εm))/(εn-εm) factor less than `droptol`. If `sternheimer_contribution`
is false, only compute excitations inside the provided orbitals.
"""
function apply_χ0(ham, δV, ψ, εF, eigenvalues; droptol=0,
                  sternheimer_contribution=true, temperature=ham.basis.model.temperature)
    if droptol > 0 && sternheimer_contribution == true
        error("Droptol cannot be positive if sternheimer contribution is to be computed.")
    end
    # δρ = ∑_nk (f'n δεn |ψn|^2 + 2Re fn ψn* δψn - f'n δεF |ψn|^2
    basis = ham.basis
    T = eltype(basis)
    model = basis.model
    fft_size = basis.fft_size
    @assert model.spin_polarization == :none
    filled_occ = filled_occupation(model)
    dVol = basis.model.unit_cell_volume / prod(basis.fft_size)

    δρ = zero(δV)
    for ik = 1:length(basis.kpoints)
        if length(basis.ksymops[ik]) != 1
            error("Kpoint symmetry not supported")
        end
        for n = 1:size(ψ[ik], 2)
            ψnk = @view ψ[ik][:, n]
            ψnk_real = G_to_r(basis, basis.kpoints[ik], ψnk)
            εnk = eigenvalues[ik][n]

            # 2Re fn ψn* δψn
            # we split δψn into its component on the computed and uncomputed states:
            # δψn = P ψn + Q ψn
            # we compute Pψn explicitly by sum over states, and
            # Q δψn by solving the Sternheimer equation
            # (H-εn) Q δψn = -Q δV ψn
            # where Q = sum_n |ψn><ψn|

            # explicit contributions
            for m = 1:size(ψ[ik], 2)
                εmk = eigenvalues[ik][m]
                ddiff = Smearing.occupation_divided_difference
                ratio = filled_occ * ddiff(model.smearing, εmk, εnk, εF, temperature)
                (n != m) && (abs(ratio) < droptol) && continue
                ψmk = @view ψ[ik][:, m]
                ψmk_real = G_to_r(basis, basis.kpoints[ik], ψmk)
                # ∑_{n,m != n} (fn-fm)/(εn-εm) ρnm <ρmn|δV>
                ρnm = conj(ψnk_real) .* ψmk_real
                weight = dVol*dot(ρnm, δV)
                δρ .+= real(basis.kweights[ik] .* ratio .* weight .* ρnm)
            end

            # Sternheimer contributions. TODO add preconditioning here
            !(sternheimer_contribution) && continue
            fnk = filled_occ * Smearing.occupation(model.smearing, (εnk-εF) / temperature)
            abs(fnk) < eps(T) && continue
            # compute δψn by solving Q (H-εn) Q δψn = -Q δV ψn
            # we err on the side of caution here by applying Q a lot,
            # there are optimizations to be made here
            Q(ϕ) = ϕ - ψ[ik] * (ψ[ik]' * ϕ)
            rhs = - Q(r_to_G(basis, basis.kpoints[ik], δV .* ψnk_real))
            function QHQ(ϕ)
                Qϕ = Q(ϕ)
                Q(ham.blocks[ik] * Qϕ - εnk * Qϕ)
            end
            J = LinearMap{eltype(ψ[ik])}(QHQ, size(ham.blocks[ik], 1))
            δψnk = cg(J, rhs)
            δψnk_real = G_to_r(basis, basis.kpoints[ik], δψnk)
            δρ .+= 2 .* fnk .* basis.kweights[ik] .* real(conj(ψnk_real) .* δψnk_real)
        end
    end

    # Add variation wrt εF
    if temperature > 0
        ldos = LDOS(εF, basis, eigenvalues, ψ, temperature=temperature)
        dos = DOS(εF, basis, eigenvalues, temperature=temperature)
        δρ .+= ldos .* dot(ldos, δV) .* dVol ./ dos
    end

    δρ
end


text:
# empirical_line_lists.jl
# generate empirical mask, based on Eric Ford's RvSpectML.jl/examples/expres_analyze_line_by_line.jl, adapted to NEID data by Alex Wise

#TODO: Make sure barycentric shifting is working in template making (i thought its only 20 m/s but its actualy closer to 1 km/s)

""" Generate an empirical mask by finding spectral lines, fiting them with gaussians with a linear offset, and filtering the fitted line params

Inputs (must be defined in param.jl):
   inst::Module - instrument used to collect the data. Works with: EXPRES

   data_path::String - string specifying the path to the data

Optional inputs:

   times_to_use::UnitRange{Int64} - indices of the data to use out of the result of inst.make_manifest(). #TODO: change this to take a list of datetimes in a universal format to mitigate risk of index mismatch

   isSolarData::Bool - boolean specifying whether or not the spectra are taken of the Sun

   output_dir::String - path from the current directory to the directory where the empricial mask should be saved (Default: "outputs/masks/")

   pipeline_plan::PipelinePlan - pipeline tracker for RvSpectML

   verbose:Bool - determines whether descriptions of each step are printed

Steps:
   1) add functions/scripts to load data: EXPRES, NEID, HARPS-N, solarDataOrNighttimeDataBoolean param, times_to_use param, make sure data are normalized
   1.5) make a function to generate times_to_use for NEID solar data, using +/- X hours of solar noon (see 1.1 in to-do list.txt)
   2) make a template spectrum using the loaded data (see :template section below)
   3) find/fit lines in template spectrum (:fit lines section below) - FUTURE - add option to fit doublets / preset list of lines
   4) select lines based on fits (already done below probs)

Outputs:
   Empirical Mask::DataFrame: empirically-fitted and filtered RV mask including (specify fields)
"""
function generateEmpiricalMask( ; output_dir::String="", pipeline_plan::PipelinePlan = PipelinePlan(), verbose::Bool=true)

   #isSolarData=false
   #output_dir = "outputs/masks/"
   #pipeline_plan = PipelinePlan()
   #verbose = true
   reset_all_needs!(pipeline_plan) #TODO: revise module PipelinePlan or my usage of the module so this line is not needed.

   if verbose println("# Reading in customized parameters from param.jl.")  end
   paths_to_search_for_param = [pwd(),joinpath(pwd(),"inputs"),joinpath(pkgdir(RvLineList),"inputs")]
   eval(read_data_paths(paths_to_search=paths_to_search_for_param, filename="param.jl"))

   global norm_type
   if norm_type != :continuum
      println("Warning: norm_type=:continuum is required by generateEmpiricalMask() but found norm_type=:", String(norm_type), ". Setting global norm_type=:continuum in order to continue.")
      norm_type = :continuum
   end

   if need_to(pipeline_plan,:read_spectra)
      if inst == EXPRES
         all_spectra = EXPRES_read_spectra(data_path)
      elseif inst == NEID
         all_spectra = combine_NEID_daily_obs(get_NEID_best_days(startDate=Date(2021,01,01), endDate=Date(2021,09,30), nBest=100))
      else
         print("Error: spectra failed to load; inst not supported.")
      end
      dont_need_to!(pipeline_plan,:read_spectra)
   end

   #make sure EXPRES spectra are normalized
   if typeof(get_inst(all_spectra)) <: AnyEXPRES   RvLineList.continuum_normalize_spectra!(all_spectra)   end

   #extract orders_to_use from all_spectra. We use remove_bad_chunks=false because neid has some nans due to a bad detector column
   order_list_timeseries = extract_orders(all_spectra, pipeline_plan, orders_to_use=orders_to_use, remove_bad_chunks=false, recalc=true);

   #remove nans from data before making template spectra - TODO: turn this into an interpolation or modify template construction to allow nans
   for i in 1:length(order_list_timeseries)
     for j in 1:length(order_list_timeseries[1])
       order_list_timeseries[i][j].flux[isnan.(order_list_timeseries[i][j].flux)] .= 1.0
     end
   end

   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)
      GC.gc()   # run garbage collector for deallocated memory
      #map(i->order_list_timeseries.metadata[i][:rv_est] = 0.0, 1:length(order_list_timeseries) )
      # Smothing is broken with GP interpolation.  Need to fix.  In mean time, here's a Sinc interpolation workaround
      @time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = RvSpectML.make_template_spectra(order_list_timeseries, smooth_factor=2.0) #, alg=:Sinc)
      #@time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = RvSpectML.make_template_spectra(order_list_timeseries, alg=:Sinc)
      if save_data(pipeline_plan, :template)
         #using JLD2, FileIO
         save(joinpath(output_dir,fits_target_str * "_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 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), ones(Int64,length(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_logy = RvSpectML.LineFinder.find_lines_in_chunklist(cl, plan=RvSpectML.LineFinder.LineFinderPlan(line_width=line_width_50,min_deriv2=0.5, use_logλ=true, use_logflux=true), verbose=false)  # TODO: Automate threshold for finding a line
      lines_in_template = RvSpectML.LineFinder.find_lines_in_chunklist(cl, plan=RvSpectML.LineFinder.LineFinderPlan(line_width=line_width_50,min_deriv2=0.5, use_logλ=true, use_logflux=false), verbose=false)  # 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)
         CSV.write(joinpath(output_dir,"linefinder",fits_target_str * "_linefinder_lines.csv"), lines_in_template )
         CSV.write(joinpath(output_dir,"linefinder",fits_target_str * "_linefinder_line_fits.csv"), fits_to_lines )
      end

      dont_need_to!(pipeline_plan,:fit_lines);
   end

   select_line_fits_with_good_depth_width_slope(fits_to_lines, parse(Float64,quant) / 100, verbose=verbose, output_dir=output_dir)

end #end function generateEmpiricalMask()

#read spectra code for EXPRES data
#TODO: move this somewhere else?
function EXPRES_read_spectra(data_path::String; verbose::Bool=false)
   if verbose println("# Finding what data files are avaliable.")  end
   if isfile(joinpath(pipeline_output_path_afw5465,"manifest.csv"))
      if verbose println("# Reading in manifest from manifest.csv") end
      df_files  = CSV.read(joinpath(pipeline_output_path_afw5465,"manifest.csv"), DataFrame)
      @assert size(df_files,1) >= 1
      @assert hasproperty(df_files,:Filename)
      @assert hasproperty(df_files,:bjd)
   else
      if verbose println("# Generating manifest file manifest.csv") end
      df_files = inst.make_manifest(data_path, inst)
      CSV.write(joinpath(pipeline_output_path_afw5465,"manifest.csv"), df_files)
   end
   df_files_use = df_files |> @take(max_spectra_to_use) |> DataFrame
   if verbose println("# Reading in ", size(df_files_use,1), " FITS files.")  end
   @time all_spectra = map(inst.read_data,eachrow(df_files_use))
end

#author: Eric Ford
#adapted from: RvSpectML.jl/examples/expres_analyze_line_by_line.jl
function select_line_fits_with_good_depth_width_slope(line_fits_df::DataFrame, quantile_threshold::Real; verbose::Bool = false, output_dir::String = "" )
   fit_distrib = line_fits_df |> @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=first(_.line_id),  frac_converged=mean(_.fit_converged)  } ) |>
            @filter(_.frac_converged == 1.0 ) |> DataFrame

   std_depth_treshold = quantile(fit_distrib.std_depth,quantile_threshold)
   median_σ_width_treshold = 2000 # quantile(sqrt.(fit_distrib.median_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,1-quantile_threshold)
   std_σ_width_treshold = quantile(sqrt.(fit_distrib.std_σ²)./fit_distrib.median_λc.*RvSpectML.speed_of_light_mps,quantile_threshold)
   std_b_treshold = quantile(fit_distrib.std_b,quantile_threshold)
   std_a_treshold = quantile(fit_distrib.std_a,quantile_threshold)
   good_lines_alt = fit_distrib |>
      @filter( 0.05 <= _.median_depth <= 0.95 ) |>
      #@filter( sqrt.(_.median_σ²)./_.median_λc.*RvSpectML.speed_of_light_mps >= median_σ_width_treshold ) |>
      @filter( sqrt.(_.std_σ²)./_.median_λc.*RvSpectML.speed_of_light_mps <= std_σ_width_treshold ) |>
      @filter( _.std_b < std_b_treshold) |>

      DataFrame
   if verbose
      println("# Found ", size(good_lines_alt,1), " good lines (std_depth_width_slope), rejected ", size(fit_distrib,1)-size(good_lines_alt,1), " lines.")
   end
   if length(output_dir) > 0
      val_str = Printf.@sprintf("%1.2f",quantile_threshold)
      airVacString = get_airVacString(inst)
      CSV.write(joinpath(output_dir,fits_target_str * "_good_lines_fit_quant=" * val_str * airVacString * ".csv"), good_lines_alt )
   end
   return good_lines_alt
end

#generate function that finds 100 days between jan 1, 2021 and sept 30, 2021, with highest number of num_rvs.good according to summary_1.csv
function get_NEID_best_days(; startDate::Date=Date(2021,01,01), endDate::Date=Date(2021,09,30), nBest::Int64=100)
   @assert endDate > startDate
   #summary_1 = CSV.read("/home/awise/data/neid/solar/summary_1.csv", DataFrame)
   summary_1 = CSV.read(joinpath(pipeline_output_path_ebf11,"summary_1.csv"), DataFrame)
   summary_1_in_date_range = summary_1[(summary_1[!,"obs_date.string"] .>= startDate) .& (summary_1[!,"obs_date.string"] .<= endDate),:]
   indices = sort(reverse(sortperm(summary_1_in_date_range[!,"num_rvs.good"]))[1:nBest])
   NEID_best_days = summary_1_in_date_range[indices,:]
   println(string("Taking all NEID days with at least ",minimum(NEID_best_days."num_rvs.good")," good RVs."))
   return NEID_best_days."obs_date.string"
end

#combine NEID clean daily mean spectra from daily_ccfs files into a vector of type Spectra2DBasic
function combine_NEID_daily_obs(dates::Vector{Date})
   date_dirs = reduce.(joinpath,split.(Dates.format.(dates,ISODateFormat),"-"))
   obs_list = joinpath.(pipeline_output_path_ebf11,date_dirs,"daily_ccfs_1.jld2")
   manifest_list = joinpath.(pipeline_output_path_ebf11,date_dirs,"manifest.csv")
   #manifest_list = joinpath.(pipeline_output_path_ebf11,date_dirs,"manifest.csv")
   #reduce(vcat,CSV.read.(manifest_list,DataFrame))
   @time obs_lambda_flux_var = [getindex.(Ref(load(fn)),["mean_lambda", "mean_clean_flux_continuum_normalized", "mean_clean_var"]) for fn in obs_list]
   daily_rvs = zeros(length(manifest_list))
   for i in 1:length(manifest_list)
      manifest_i = CSV.read(manifest_list[i], DataFrame)
      index_minimum_hour_angle = argmin(abs.(manifest_i[:,"hour_angle"]))
      daily_rvs[i] = manifest_i[index_minimum_hour_angle,"ssbz"] * speed_of_light_mps
   end
   obs_lambda = [obs[1] for obs in obs_lambda_flux_var]
   obs_flux = [obs[2] for obs in obs_lambda_flux_var]
   obs_var = [obs[3] for obs in obs_lambda_flux_var]
   obs_metadata = [Dict{Symbol,Any}(:date=>dates[i], :rv_est=>daily_rvs[i], :bjd=>datetime2julian(DateTime(Dates.year(dates[i]), Dates.month(dates[i]), Dates.day(dates[i]), 12))) for i in 1:length(dates)]
   all_spectra = map(obs -> Spectra2DBasic(obs[1], obs[2], obs[3], NEID2D(), metadata=obs[4]), eachrow(hcat(obs_lambda,obs_flux,obs_var,obs_metadata))) #I think there's a better way to do this, perhaps by using obs_lambda_flux_var direclty
   return all_spectra
end

max_stars_repo_name: alexander-wise/RvLineList

max_issues_repo_name: alexander-wise/RvLineList

max_forks_repo_name: alexander-wise/RvLineList


text:
# Exercise 3-1
# repeatlyrics()

function repeatlyrics()
    printlyrics()
    printlyrics()
end

function printlyrics()
    println("I'm a lumberjack, and I'm okay.")
    println("I sleep all night and I work all day.")
end

# Function definition needs to be before the call if calling in global scope.
# If calling from within another function, it can be before.

# Exercise 3-2
function rightjustify(s)
    println(" "^(70-length(s)) * s)
end
# rightjustify("monty")

# Exercise 3-3
function printtwice(bruce)
    println(bruce)
    println(bruce)
end

function dotwiceold(f)
    f()
    f()
end

function dotwice(f, val)
    f(val)
    f(val)
end

function printspam()
    println("spam")
end

# println("dotwice(printtwice, \"spam\")")
# dotwice(printtwice, "spam")

function dofour(f, val)
    dotwice(f, val)
    dotwice(f, val)
end

# println("dofour(printtwice, \"spam\")")
# dofour(printtwice, "spam")

# Exercise 3-4
function printmiddlerow()
    println("|", " "^9, "|", " "^9, "|")
end

function printgrid()
    println("+ - - - - + - - - - +")
    dotwiceold(printmiddlerow)
    dotwiceold(printmiddlerow)
    println("+ - - - - + - - - - +")
    dotwiceold(printmiddlerow)
    dotwiceold(printmiddlerow)
    println("+ - - - - + - - - - +")
end

# printgrid()

function printsmallmiddlerow()
    println("|", " "^5, "|")
end

function printsmallgrid()
    println("+ - - +")
    dotwiceold(printsmallmiddlerow)
    println("+ - - +")
end

printsmallgrid()