text:
# --------------------- Functions for normalizing tabular datasets --------------------- */
"""
    normalize!(type::Type{Std}, data, mean, std; kwargs...)

Returns standardized data by the specified mean and standard deviation.

If `data` is an `AbstractArray`, in keyword argument `dims` dimensions can be specified
over which mean and standard deviation are computed.
"""
function normalize!(type::Type{Std}, data, mean, std; kwargs...)
    change!(data, mean, std)
    return
end

"""
    normalize!(type::Type{Std}, data; kwargs...)

Returns standardized data by the mean and variance of `data`.

If `data` is an `AbstractArray`, in keyword argument `dims` dimensions can be specified
over which mean and standard deviation are computed.
"""
function normalize!(type::Type{Std}, data; kwargs...)
    mean, std = meanstd(data; kwargs...)
    normalize!(type, data, mean, std; kwargs...)
    return
end

"""
    normalize!(type::Type{MinMax}, data, min, max; kwargs...)

Returns min-max scaled data based on specified minimum and maximum values.

If `data` is an `AbstractArray`, in keyword argument `dims` dimensions can be specified
over which (columns or rows) minimum and maximum are computed.
"""
function normalize!(type::Type{MinMax}, data, min, max; kwargs...)
    change!(data, min, max)
    return
end

"""
    normalize!(type::Type{MinMax}, data; kwargs...)

Returns min-max scaled data based on minimum and maximum values from `data`.

If `data` is an `AbstractArray`, in keyword argument `dims` dimensions can be specified
over which mean and standard deviation are computed.
"""
function normalize!(type::Type{MinMax}, data; kwargs...)
    min, max = minmax(data; kwargs...)
    normalize!(type, data, min, max; kwargs...)
    return
end

"""
    normalize!(type::Type{L2}, data, norm)

Returns normalized data by the specified `norm`.

If `data` is an `AbstractArray`, in keyword argument `dims` dimensions can be specified
over which norm is computed. For more see [`l2norm`](@ref)
"""

function normalize!(type::Type{L2}, data, norm; kwargs...)
    change!(data, norm)
    return
end

"""
    normalize!(type::Type{L2}, data; kwargs...)

Returns normalized data by the Euclidan norm of `data` elements.

If `data` is an `AbstractArray`, in keyword argument `dims` dimensions can be specified
over which mean and standard deviation are computed.
"""
function normalize!(type::Type{L2}, data; kwargs...)
    norm = l2norm(data; kwargs...)
    normalize!(type, data, norm; kwargs...)
    return
end

function normalize!(data; kwargs...)
    error("Type of normalization not specified, try Std, MinMax or L2.")
end

# ------------------ Util functions for normalizing tabular datasets ------------------- */
"""
    _select(data, selection)

    Returns data, range specified by `selection`.
"""
_select(data::AbstractArray, selection) = data[selection[1], selection[2]]

_select(data::AbstractDataFrame, selection) = data[!, selection]

"""
    meanstd(data; dims::Int=1)

Returns mean and standard deviation of columns or rows of given `data`.

Mean and standard deviation are returned in a row vector (`Tuple` of two `1xn Matrix`).

`data` can be an `AbstractArray` or `AbstractDataFrame`.
If `data` is an `AbstractArray`, a keyword argument `dims` can be provided to compute values
over dimensions.

This function uses `mean` and `std` from the package `Statistics`.
"""
function meanstd(data::AbstractArray; dims::Int=1)
    selected = _select(data, selectnumeric(data; dims=dims))
    mean(selected; dims=dims), std(selected; dims=dims)
end
function meanstd(data::AbstractDataFrame)
    selected = _select(data, selectnumeric(data))
    mean.(eachcol(selected)), std.(eachcol(selected))
end

"""
    minmax(data; dims::Int=1)

Returns minimum and maximum of columns or rows of given `data`.
If `data` is an `AbstractArray`, `dims` specifies whether it is over columns or rows.
"""
function minmax(data::AbstractArray; dims::Int=1)
    selected = _select(data, selectnumeric(data))
    minimum(selected; dims=dims), maximum(selected; dims=dims)
end

function minmax(data::AbstractDataFrame)
    selected = _select(data, selectnumeric(data))
    minimum.(eachcol(selected)), maximum.(eachcol(selected))
end

"""
    l2norm(data; dims::Int=1)

Returns norm of columns or rows of given `data`.

If `data` is an `AbstractArray`, in keyword argument `dims` dimensions can be specified
over which norm is computed. If `data` is an `AbstractDataFrame`, then norm is calculated
for columns.

This function uses `norm` from the package `LinearAlgebra`.
"""
function l2norm(data::AbstractArray; dims::Int=1)
    selected = _select(data, selectnumeric(data))
    if dims == 1    #by columns
        return makerow(norm.(eachcol(selected)))
    elseif dims == 2
        return norm.(eachrow(selected))
    else
        throw(ArgumentError("Wrong dimensions."))
    end
end

function l2norm(data::AbstractDataFrame)
    selected = _select(data, selectnumeric(data))
    return norm.(eachcol(selected))
end

"""
    selectnumeric(data::AbstractDataFrame)

Returns array of number of columns which contain values of type `Number`.
"""
function selectnumeric(data::AbstractDataFrame)
    n = Base.size(data, 2)
    nums = Int[]

    for i in 1:n-1
        if eltype(data[:,i]) <: Number
            append!(nums, i)
        end
    end

    return nums
end

"""
    selectnumeric(data::AbstractArray; dims::Int=1)

Returns tuple array of number of columns or rows which contain values of type `Number` and
other element of the tuple is `Colon` :.
"""
function selectnumeric(data::AbstractArray; dims::Int=1)
    nums = Int[]

    if dims == 1
        n = Base.size(data, 2)

        for i in 1:n
            if eltype(data[1,i]) <: Number && eltype(data[2,i]) <: Number
                append!(nums, i)
            end
        end

        return (:, nums)

    elseif dims == 2
        n = Base.size(data, 1)

        for i in 1:n
            if eltype(data[i,1]) <: Number && eltype(data[i,1]) <: Number
                append!(nums, i)
            end
        end

        return (nums,:)

    else
        throw.ArgumentError("Unsopported dimensions. (1 or 2 are accepted)")
    end
end

"""
    makerow(vec::Matrix)

For matrix permutes dimensions (makes a row vector from column vector). If `vec` is already
a row vector then does nothing.
"""
makerow(vec::AbstractArray) = Base.size(vec, 1) == 1 ? vec : permutedims(vec)

"""
    change(data, substracted, divided)

Changes `data` elements following a formula: `data .- substracted ./ devided`.
"""
function change!(data::AbstractArray, substracted, devided)
    selection = selectnumeric(data)
    data[selection[1], selection[2]] .= (data[selection[1], selection[2]] .- substracted) ./ devided
    return
end

function change!(data::AbstractDataFrame, substracted, devided)
    for (i, col) in enumerate(selectnumeric(data))
        data[!, col] = (data[:,col] .- substracted[i]) ./ devided[i]
    end

    return
end

"""
    change!(data, divided)

Changes `data` elements following a formula: `data ./ devided`.
"""
function change!(data::AbstractArray, devided)
    selection = selectnumeric(data)
    data[selection[1], selection[2]] .= data[selection[1], selection[2]] ./ devided
    return
end

function change!(data::AbstractDataFrame, devided)
    for (i, col) in enumerate(selectnumeric(data))
        data[!, col] = data[:,col] ./ devided[i]
    end

    return
end

# ------------------ Util functions for binarizing tabular datasets -------------------- */
"""
    classes(dataset)

Prints unique names of target values.
"""
classes(dataset::DatasetName) = sort(unique(labels(dataset)))

"""
    binarize(data, pos_labels)

Binarize data for specified positive labels. (Data must contain just labels no feature valeus.)
"""
binarize(data, pos_labels::AbstractString) = binarize(data, [pos_labels])
binarize(data, pos_labels) = [in(i, pos_labels) for i in data]

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
module Serializer

import SearchLight

function serialize end

function deserialize end

function serializables end # return Vector{Symbol}

include("../serializers/JsonSerializer.jl")

end

ext: jl


text:
#=
    Utility functions for testing the DataFlowTasks package.
=#
using DataFlowTasks: R,W,RW

# define a simple function to mimic a problem with fork-join parallelization.
# The execution starts with one node, spawns `n` independent nodes, and then
# joint them later at a last node. The  computation waits for the last node, and
# each block works for `s` seconds
function fork_join(n,s)
    A = rand(2n)
    @dspawn do_work(s) (A,) (RW,)
    for i in 1:n
        @dspawn do_work(s) (view(A,[i,i+n]),) (RW,)
    end
    res = @dspawn do_work(s) (A,) (R,)
    return res
end

function do_work(t)
    ti = time()
    while (time()-ti) < t end
    return
end


text:
# # Example: Dying philosopher

using Julio
using ContextManagers: @with, onfail

# Following function demonstrate the [Dining philosophers problem -
# Wikipedia](https://en.wikipedia.org/wiki/Dining_philosophers_problem) and a
# "quick-and-dirty" way to recover from the deadlock by killing a philosopher
# after a timeout.
#
# This example is based on
# [David Beazley - Die Threads - YouTube](https://www.youtube.com/watch?v=xOyJiN3yGfU)
# which demonstrated this with a threading library for Python called
# [Thredo](https://github.com/dabeaz/thredo).

function dying_philosopher(nphils = 5)
    sticks = [Julio.Lock() for _ in 1:nphils]
    Julio.withtaskgroup() do tg
        phils = [Julio.cancelscope() for _ in 1:nphils]
        for (n, scope) in pairs(phils)
            Julio.spawn!(tg) do
                open(scope) do
                    @with(sticks[n]) do
                        sleep(rand())
                        @with(
                            onfail() do
                                println("But what is death? $n")
                            end,
                            sticks[mod1(n + 1, nphils)]
                        ) do
                            println("eating $n")
                        end
                    end
                end
            end
        end
        Julio.spawn!(tg) do
            sleep(10)
            n = rand(1:nphils)
            println("Killing philosopher $n")
            Julio.cancel!(phils[n])
        end
    end
end

# Running the above function prints something like
#
# ```julia
# julia> dying_philosopher();
# Killing philosopher 4
# But what is death? 4
# eating 3
# eating 2
# eating 1
# eating 5
# ```
#
# Of course, this is not a "correct" solution to the dining philosophers problem
# *per se*. However, it is easy to accidentally write a program producing
# deadlocks and discovering this in an interactive session. It is important that
# we can reliably recover from a deadlock with `Ctrl+C` without breaking the
# entire Julia runtime so that we can debug the program. It helps writing and
# fixing concurrent programs.


text:
# Assume that the prism axis is orthogonal to the prism base.  This makes it much easier to
# implement surfpt_nearby by separating the axis and base dimensions (i.e., the closest
# point on the surface can be found independently in axis and base dimensions).
#
# Supporting skewed prisms is not impossible, but this will be done in the future if the
# demand is high.

# It is possible to combine Cylinder and Prism by defining a more general prism type that
# accepts any base shape.  We could use Shape{2} to create a 2D base shape and store it as
# a base shape.  However there is one problem.  We want to use
#
#   surfpt_nearby(x_base::SVector{2}, shape_base::Shape{2})
#
# to find the closest point in the base dimension (such a point is on the side of the prism),
# and compare the distance to it with the distance from x to the base plane along the axis
# dimension.  Then, if x is inside the prism for example, then this comparison will tell us
# which of the base and side is the closest surface from x.  However, this means that we
# have to take the difference between the output of the above function and x_base to figure
# out the distance in the base dimension.  Provided that for many Shape types ∆x is first
# calculated inside surfpt_nearby and x+∆x is returned as the found surface point, having to
# recover ∆x by subtraction is not only inefficient, but also suffers from loss of
# significant digits.
#
# Maybe we should have not implemented surfpt_nearby, but a different function that returns
# ∆x, because the surface point x + ∆x can be easily calculated.  Let's consider fixing this
# in the future.
#
# Hmmm...  Because the surface point is not calculated by adding ∆x to x for Ball and
# Ellipse, maybe implemeting a function returning ∆x is not a very good idea.  Then, it may
# not be a bad idea to calculate ∆x by subtracting x from the surface point.  Let's try to
# implement the general prism.

export Prism

mutable struct Prism{B<:Shape{2}} <: Shape{3,9}
    c::SVector{3,Float64}  # prism center
    b::B  # base shape described in prism coordinates (i.e, when translating prism, do not need to translate b)
    h2::Float64  # height * 0.5
    p::SMatrix{3,3,Float64,9}  # projection matrix to prism coordinates; must be orthonormal (see surfpt_nearby)
    Prism{B}(c,b,h2,p) where {B} = new(c,b,h2,p)  # suppress default outer constructor
end

Prism(c::SVector{3,<:Real},
      b::B,
      h::Real=Inf,
      axes::SMatrix{3,3,<:Real}=SMatrix{3,3,Float64}(I)  # columns are axes vectors: first two columns span prism base, and last column is prism axis
      ) where {B<:Shape{2}} =
    Prism{B}(c, b, 0.5h, inv(axes ./ sqrt.(sum(abs2,axes,dims=Val(1)))))

Prism(c::AbstractVector{<:Real}, b::Shape{2}, h::Real=Inf, axes::AbstractMatrix{<:Real}=Matrix{Float64}(I,length(c),length(c))) =
    Prism(SVector{3}(c), b, h, SMatrix{3,3}(axes))

Base.:(==)(s1::Prism, s2::Prism) = s1.c==s2.c && s1.b==s2.b && s1.h2==s2.h2 && s1.p==s2.p
Base.isapprox(s1::Prism, s2::Prism) = s1.c≈s2.c && s1.b≈s2.b && s1.h2≈s2.h2 && s1.p≈s2.p
Base.hash(s::Prism, h::UInt) = hash(s.c, hash(s.b, hash(s.h2, hash(s.p, hash(:Prism, h)))))

function level(x::SVector{3,<:Real}, s::Prism)
    y = s.p * (x - s.c)  # coordinates after projection
    ya = y[3]  # scalar: coordinate in axis dimension
    yb = y[SVector(1,2)]  # SVector{2}: coordinate in base dimensions

    return max(abs(ya)/s.h2 - 1.0, level(yb,s.b))
end

function surfpt_nearby(x::SVector{3,<:Real}, s::Prism)
    ax = inv(s.p)  # prism axes: columns are not only unit vectors, but also orthogonal

    y = s.p * (x - s.c)  # x in prism coordinates
    ya = y[3]  # scalar: coordinate in axis dimension
    yb = y[SVector(1,2)]  # SVector{2}: coordinates in base dimensions

    la = abs(ya)
    abs∆a = abs(s.h2 - la)  # scalar: distance between x and base point closest to x
    surfa = SVector(yb.data..., copysign(s.h2, ya))  # SVector{3}: coordinates of base point closest to x
    nouta = SVector(0.0, 0.0, copysign(1.0, ya))  # SVector{3}: outward direction normal at surfa
    onbnda = abs∆a ≤ Base.rtoldefault(Float64) * s.h2
    isouta = s.h2<la || onbnda

    surfb2, noutb2 = surfpt_nearby(yb, s.b)  # (SVector{2}, SVector{2}): side point closest to x and outward direction normal to side there
    abs∆b = norm(surfb2 - yb)  # scalar: distance between x and side point closest to x
    surfb = SVector(surfb2.data..., ya)  # SVector{3}: coordinates of side point closest to x
    noutb = SVector(noutb2.data..., 0.0)  # SVector{3}: outward direction normal to side surface at surfb
    basesize = abs.((-)(bounds(s.b)...))  # SVector{2}: size of bounding rectancle of base
    onbndb = abs∆b ≤ Base.rtoldefault(Float64) * max(basesize.data...)
    isoutb = yb∉s.b || onbndb

    if isouta && isoutb  # x outside in both axis and base dimensions
        surf = SVector(surfb[1], surfb[2], surfa[3])
        nout = (onbnda && onbndb) ? (noutb + nouta) : (y - surf)
        nout = norm(nout)==Inf ? isinf.(nout) .* sign.(nout) : normalize(nout)  # e.g., return [0,0,-1] for nout = [1,-2,-Inf]
    elseif !isouta && isoutb  # x outside in base dimensions, but inside prism in axis dimension
        (surf, nout) = (surfb, noutb)
    elseif isouta && !isoutb # x outside in axis dimension, but inside prism in base dimensions
        (surf, nout) = (surfa, nouta)
    else  # !isouta && !isoutb: x strictly inside prism
        (surf, nout) = (abs∆a ≤ abs∆b) ? (surfa, nouta) : (surfb, noutb)
    end

    return ax*(surf+s.c), ax*nout
end

function bounds(s::Prism)
    ax = inv(s.p)  # prism axes: columns are not only unit vectors, but also orthogonal
    a = ax[:,3]  # SVector{3}
    h2a = s.h2 * a

    l0, u0 = bounds_ctrcut(s)  # (SVector{3}, SVector{3})
    l1, u1 = l0+h2a, u0+h2a  # (SVector{3}, SVector{3})
    l2, u2 = l0-h2a, u0-h2a  # (SVector{3}, SVector{3})

    return min.(l1,l2)+s.c, max.(u1,u2)+s.c
end

include("cylinder.jl")
include("polygonal.jl")
include("sectoral.jl")


text:
module DoubletFinder

import Pkg
import PyCall

# export functions

export findV
export findZAbs
export getblocks1
export getblocks2
export sep_blocks
export three_sigmadet
export atomicdata
export rb_limits
export df_interest
export datacleanandtag
export iscontained
export getINFO
export possible_doublets_parallel
export possible_doublets
export Find

using Distributed
@everywhere using Distributed
using Pkg
@everywhere using Pkg
Pkg.activate(".")
@everywhere Pkg.activate(".")
Pkg.instantiate()
using CSV
@everywhere using CSV
using PyCall
@everywhere using PyCall
using Test
@everywhere using Test
using Statistics
@everywhere using Statistics
using DataFrames, Query, DataFramesMeta
@everywhere using DataFrames, Query, DataFramesMeta
using Images
@everywhere using Images
using LinearAlgebra
@everywhere using LinearAlgebra
using DSP
@everywhere using DSP
using Printf
@everywhere using Printf
ss = pyimport("scipy.signal")
@everywhere ss = pyimport("scipy.signal")
using Distributions
@everywhere using Distributions

@everywhere const c_light=299792.458; # speed of light in km/s
@everywhere const LYA = 1215.6; # Lyman alpha rest-wavelength in angstrom

@everywhere begin
    """
        findV(zAbs, zEm)

    "Computes the velocity at a redshift, zAbs, relative to the emission redshift, zEm."

    # Examples
    ```julia-repl
    julia> findV(3.12,3.12)
    0.0
    ```
    """
    function findV(zAbs::Float64, zEm::Float64)
        @assert length(zAbs) == length(zEm)
        v = ((((1.0+zEm)^2) - ((1.0+zAbs)^2))/(((1.0+zEm)^2) + ((1.0+zAbs)^2))) * c_light
        return v
    end
end

@everywhere begin
    """
        findZAbs(v, zEm)

    "Computes the redshift at a given velocity, v, relative to the emission redshift, zEm."

    # Examples
    ```julia-repl
    julia> findZAbs(0, 3.12)
    3.12
    ```
    """
    function findZAbs(v::Float64, zEm::Float64)
        @assert length(v) == length(zEm)
        temp = ((1.0+zEm)^2)*((c_light-v)/(c_light+v))
        zAbs = temp^0.5 - 1.0
        return zAbs
    end
end

@everywhere begin
    """
        getblocks1(vs)

    "Returns an array of tuple_arrays that are indices of the input vector, vs, which have values less than 1"

    """
    function getblocks1(vs::Array{Float64,1})
        @assert length(vs)>0
        blocks = Tuple{Int, Int}[]
        inblock, start = false, 0, 0
        for (i, v) in enumerate(vs)
            if inblock
                if v >= 1.0
                    push!(blocks, (start, i-1))
                    inblock = false
                end
            else
                if v < 1.0
                    start = i
                    inblock = true
                end
            end
        end
        inblock && push!(blocks, (start, length(vs)))
        blocks
    end
end

@everywhere begin
    """
        getblocks2(vs)

    "Returns an array of tuple_arrays that are indices of the input vector, vs, which have values less than 1. A more Pythonic implementation."

    """
    function getblocks2(vs::Array{Float64,1})
        @assert length(vs)>0
        t = [false; vs .< 1.0; false]
        dt = diff(t)
        f = findall(==(1), dt)
        l = findall(==(-1), dt) .- 1
        collect(zip(f, l))
    end
end

@everywhere begin
    """
        getblocks2(vs)

    "Returns an array of tuples, where each tuple is separated by more than 1 element"

    """
    function sep_blocks(blocks::Array{Tuple{Int64,Int64},1})
        n_blocks=Array{Tuple{Int64,Int64},1}(undef,length(blocks))
        n_n_blocks=1
        for i in 1:length(blocks)
            if blocks[i][2]-blocks[i][1]>1
                n_blocks[n_n_blocks]=blocks[i]
                n_n_blocks=n_n_blocks+1
            end
        end
        resize!(n_blocks,n_n_blocks-1)
        return n_blocks
    end
end

@everywhere begin
    """
        threesigmadet(df,tuple_array)

    Returns true if the trough (absorption line) is a 3 sigma detection. Takes input as df, the data frame containing the spectral data, and the `tuple_array` which identifies the indices corresponding to beginning and ending of the trough.

    """
    function threesigmadet(df::DataFrame,tuple_array::Tuple{Int64,Int64})
        @assert size(df,2)>0
        @assert length(tuple_array)>0

        val=minimum(df[:FLUX][tuple_array[1]:tuple_array[2]])
        loc_min=findall(x->x==val,df[:FLUX])[1]
        if (1-val)/df[:ERROR][loc_min][1] >= 3
            return true
        else
            return false
        end
    end
end

@everywhere begin
    """
        atomicdata(all_ions)

    Returns a tuple consisting of dictionary of species in `all_ions`. The file `atoms_UVES.csv` is read
    from the current path. It is a library containing information on commonly observed absorption lines in a quasar spectrum. It contains information such as rest-wavelength, oscillator strength, atomic mass. The file is ingested into DataFrame format for further manipulation. `all_ions` denotes the most commonly observed ions; it is an array string, e.g. ["NV","CIV","SiIV","OVI"]

    """
    function atomicdata(all_ions::Array{String,1}=["NV","CIV","SiIV","OVI","MgII","MgI","FeII","AlII","AlIII","MnII","CrII","ZnII","Lya","HI","NeVIII"])

        # The data file containing atomic information of interest
        path=dirname(pathof(DoubletFinder))
        path_atomicfile=joinpath(path,"atoms_UVES.csv")

        @assert isfile(path_atomicfile)
        df_atomic_data=CSV.read(path_atomicfile,header=["Ion","Rest_Wavelength","f","Gamma","Mass","Column6"],delim='\t')

        #delete unused column
        deletecols!(df_atomic_data, :Column6)

        #create a new column with only the Ion Identifier name
        df_atomic_data[:Identifier_Ion]=[split(i,r"[0-9]+")[1] for i in df_atomic_data[:Ion]]

        #create a new column with only the Wavelength Identifier name
        df_atomic_data[:Identifier_Wavelength]=[split(i,r"[a-zA-Z]+")[2] for i in df_atomic_data[:Ion]]

        #get the properties of all_ions from the atomic_data table
        species = Dict{String,Dict{String,Tuple{Float64,Float64}}}()
        for search_trans in all_ions
            #species[search_trans] = Dict()
            species[search_trans] = Dict{String,Tuple{Float64,Float64}}()
            for i in 1:1:length(df_atomic_data[1])
                atom = split(df_atomic_data[:Ion][i],r"[0-9]+")[1]
                tranwav = split(df_atomic_data[:Ion][i],r"[a-zA-Z]+")[2]
                current_oscillator = df_atomic_data[:f][i]
                if ( atom == search_trans )
                    species[atom][tranwav] = df_atomic_data[:Rest_Wavelength][i] , current_oscillator
                end
            end
        end

        for i in species
            species[i[1]]=sort(species[i[1]])
        end
        return species
    end
end

@everywhere begin
    """
        rb_limits(wavelength,zEm,species)

    Returns a tuple of blue limit and red limit of the spectrum. Takes as input a wavelength vector, emission redshift, and the dictionary of species.

    """
    function rb_limits(wavelength::Array{Float64,1},zEm::Float64,species::Dict{String,Dict{String,Tuple{Float64,Float64}}})
        # determining the blue limit
        @assert length(wavelength)>0
        @assert zEm>0
        if LYA * (1. + zEm) > wavelength[1]
           blue_limit = LYA * (1. + zEm)
        else
           blue_limit = wavelength[1]
        end
        #Emission redshift of quasar - 3000 km/s
        emission_limit = 1. + zEm - 0.01
        red_limit =  Dict{String,Float64}()
        for specie in species
             # Choose either the reddest_transition or the end of the spectrum
            reddest = maximum(species[specie[1]])[1]
            if ((species[specie[1]][reddest][1] * emission_limit) < wavelength[end])
                red_limit[specie[1]] = species[specie[1]][reddest][1] * emission_limit
            else
                red_limit[specie[1]] = wavelength[end]
            end
        end
        max_red = maximum(values(red_limit));
        return blue_limit, max_red
    end
end

@everywhere begin
    """
        df_interest(df,zEm,blue_limit,red_limit)

    Returns a DataFrame within the limits of interest. Takes as input the full DataFrame, emission redshift, and the limits.

    """
    function df_interest(df::DataFrame,zEm::Float64,blue_limit::Float64,max_red::Float64,
            species::Dict{String,Dict{String,Tuple{Float64,Float64}}})
        @assert max_red>blue_limit
        @assert zEm>0
        wavelengths=df[:WAVELENGTH]::Array{Float64,1}

        #dataframe of interest between the blue and red limits
        df_int=df[blue_limit.<= wavelengths .<=max_red,:]
        #adding the column for rest wavelength
        df_int[:Rest_Wavelength]=df_int[:WAVELENGTH]/(1.0+zEm);
        #create velocity columns appended to the dataframe for each of the species in the library
        for specie in species
            for transition in keys(specie[2])
                df_int[Symbol(join([specie[1],transition,"_Vel"]))]=c_light*(df_int[:Rest_Wavelength].^2 .- species[specie[1]][transition][1]^2)./(df_int[:Rest_Wavelength].^2 .+ species[specie[1]][transition][1]^2)
            end
        end
        return df_int
    end
end

@everywhere begin
    """
        datacleanandtag(df,telluriclines)

    Returns a DataFrame with telluric regions set to a value of 1 and data is cleaned for negative errors and fluxes.

    """
    function datacleanandtag(df::DataFrame,telluriclines::Array{Tuple{Float64,Float64},1}=[(9300.0,9630.0),(7594.0,7700.0),(6868.0,6932.0),(6277.0,6318.0)])
        # clean cosmic rays/bad high flux values
        df[:ERROR][df[:FLUX].>1.5].=0
        df[:FLUX][df[:FLUX].>1.5].=1
        for i in telluriclines
            # Avoid Telluric features, in order from strongest to weakest by setting those regions to a value of 1.
            df[:FLUX][(df[:WAVELENGTH].>i[1]) .& (df[:WAVELENGTH].<i[2])].=1;
        end
        # Remove negative errors
        df[:FLUX][df[:ERROR].<0] .= 1;
        df[:ERROR][df[:ERROR].<0] .= 0;
        # Remove negative flux
        df[:ERROR][df[:FLUX].<0] .= 0;
        df[:FLUX][df[:FLUX].<0] .= 0;
        # tag regions with flux less than 1
        df[:TAG]=df[:FLUX] .< 1.0;
        df[:MASK]=zeros(size(df[:TAG])[1]);

        df[:MASK][df[:TAG].==true] .=1;
        df[:MASK][df[:TAG].==false] .=0;
        return df
    end
end

@everywhere begin
    """
        iscontained(array,value)

    Returns true if the inputted array contains the value of interest within it.

    """
    function iscontained(array::Array{Float64,1},value::Float64)
        if minimum(array) <= value <= maximum(array)
            return true
        else
            return false
        end
    end
end

@everywhere begin
    """
        getINFO(df,zEm,blocks,indice,ion,transition,species_dictionary)

    Returns the indice, minimum velocity, maximum velocity, median redshift of the trough, an array of redshifts corresponding to components if a doublet is detected. Takes as input a dataframe, emission redshift, the tuple of troughs, index for the tuple describing a trough, ion of interest, transition of interest, and the `species_dictionary`
    """
    function getINFO(df::DataFrame,zEm::Float64,blocks::Array{Tuple{Int64,Int64},1},indice::Int64,ion::AbstractString,transition::String,species::Dict{String,Dict{String,Tuple{Float64,Float64}}}=atomicdata(df,zEm))

        d=df[blocks[indice][1]:blocks[indice][2],:]

    #     if length(d[:FLUX])>=5
    #         x=savitzkyGolay(Float64[ a for a in d[:FLUX] if !ismissing(a) ],5,1) #smoothing the array
    #     elseif length(d[:FLUX])>=3 & length(d[:FLUX])<5
    #         x=savitzkyGolay(Float64[ a for a in d[:FLUX] if !ismissing(a) ],3,1) #smoothing the array
    #     else
    #         x=d[:FLUX]
    #     end

        length(d[:FLUX])>=5 ? x=ss.savgol_filter(d[:FLUX],5,1) : x=ss.savgol_filter(d[:FLUX],3,1)

        xs=findlocalminima(x)
        yvals = [idx[1] for idx in xs]
        xvals=d[Symbol(join([ion,transition,"_Vel"]))][yvals]

        ele_list=Array{String}(undef,length(keys(species[ion])))

        n_ele_list=1

        for i in keys(species[ion])
            ele_list[n_ele_list]= i
            n_ele_list=n_ele_list+1
        end
        resize!(ele_list, n_ele_list-1)

        filter!(e->e != minimum(ele_list),ele_list)

        zlocs=Array{Float64}(undef, length(xvals)*length(ele_list))

        n_zlocs=1
        for j in xvals
            for ele in ele_list
                if length(findall(x->x==true,iscontained.([df[Symbol(join([ion,ele,"_Vel"]))][blocks[i][1]:blocks[i][2]] for i in 1:length(blocks)],j))) > 0
                    zlocs[n_zlocs]= findZAbs(-j, zEm)
                    n_zlocs=n_zlocs+1
                end
            end
        end
        resize!(zlocs, n_zlocs-1)

        if transition==minimum(species[ion])[1]
            return (indice,blocks[indice],minimum(d[Symbol(join([ion,transition,"_Vel"]))]),maximum(d[Symbol(join([ion,transition,"_Vel"]))]),
                Statistics.median(d[:WAVELENGTH])/species[ion][transition][1]-1,zlocs)
        else
            return (indice,blocks[indice],minimum(d[Symbol(join([ion,transition,"_Vel"]))]),maximum(d[Symbol(join([ion,transition,"_Vel"]))]),
                Statistics.median(d[:WAVELENGTH])/species[ion][transition][1]-1,NaN)
        end
    end
end

@everywhere begin
    """
        possible_doublets_parallel(df,zEm,blocks,species_dictionary,search_ions)

    Returns an array of tuples with the possible doublets and the other information of indice location, maximum and minimum velocities of the troughs, and component redshifts. Parallel version.

    """
    function possible_doublets_parallel(df::DataFrame,zEm::Float64,blocks::Array{Tuple{Int64,Int64},1},dict_atomicdata::Dict{String,Dict{String,Tuple{Float64,Float64}}},search_ions::Dict{String,Tuple{Int64,Int64}}=Dict("NV" => (-70000,10000), "CIV" => (-70000,10000), "SiIV" => (-70000,10000),"OVI" => (-70000,10000)))
        doublet_array=Dict{String,Array{Tuple{Float64,Tuple{Int64,Int64},Float64,Float64,Float64,Array{Float64,1}},1}}()
        for specie in keys(search_ions)
            doublet_array[specie]=pmap(y->getINFO(df,zEm,blocks,y,specie,minimum(keys(dict_atomicdata[specie])),dict_atomicdata),1:length(blocks),batch_size=Int(round(length(blocks)/nworkers())))
            doublet_array[specie] = [x for x in doublet_array[specie] if (x[3]>=search_ions[specie][1]) & (x[4]<=search_ions[specie][2]) & (length(x[6])!=0)]
        end
        return doublet_array
    end
end

@everywhere begin
    """
        possible_doublets(df,zEm,blocks,species_dictionary,search_ions)

    Returns an array of tuples with the possible doublets and the other information of indice location, maximum and minimum velocities of the troughs, and component redshifts. Serial version.

    """
    function possible_doublets(df::DataFrame,zEm::Float64,blocks::Array{Tuple{Int64,Int64},1},dict_atomicdata::Dict{String,Dict{String,Tuple{Float64,Float64}}},search_ions::Dict{String,Tuple{Int64,Int64}}=Dict("NV" => (-70000,10000), "CIV" => (-70000,10000), "SiIV" => (-70000,10000),"OVI" => (-70000,10000)))
        doublet_array=Dict{String,Array{Tuple{Float64,Tuple{Int64,Int64},Float64,Float64,Float64,Array{Float64,1}},1}}()
        for specie in keys(search_ions)
            doublet_array[specie]=map(y->getINFO(df,zEm,blocks,y,specie,minimum(keys(dict_atomicdata[specie])),dict_atomicdata),1:length(blocks))
            doublet_array[specie] = [x for x in doublet_array[specie] if (x[3]>=search_ions[specie][1]) & (x[4]<=search_ions[specie][2]) & (length(x[6])!=0)]
        end
        return doublet_array
    end
end

@everywhere begin
    """
    Find(file,zEm,search_ions)

    Returns an array of tuples with the possible doublets and the other information of indice location, maximum and minimum velocities of the troughs, and component redshifts.
    """
    function Find(file::String,zEm::Float64,search_ions::Dict{String,Tuple{Int64,Int64}})
        df_spectral_data = CSV.read(file,header=true,delim='\t',allowmissing=:none);
        # dictionary of species
        species_dict=atomicdata();
        # red and blue limits
        rb_limit=rb_limits(df_spectral_data[:WAVELENGTH],zEm,species_dict);
        # DataFrame of interest
        df_int=df_interest(df_spectral_data,zEm,rb_limit[1],rb_limit[2],species_dict);
        # Cleaned DataFrame of interest
        df_clean=datacleanandtag(df_int);
        # Locations of possible absorption troughs
        blocks=getblocks1(df_clean[:FLUX]);
        blocks_3s=sep_blocks(blocks)[pmap((x->threesigmadet(df_clean,x)),sep_blocks(blocks),batch_size=Int(round(length(sep_blocks(blocks))/nworkers())))];
        return possible_doublets_parallel(df_clean,zEm,blocks_3s,species_dict,search_ions)
    end
end

end #module


text:
push!(LOAD_PATH, "../src/")
using Documenter, SeqSpace

makedocs(
    sitename = "SeqSpace.jl",
    authors  = "Nicholas Noll, Madhav Mani, Boris Shraiman",
    modules  = [SeqSpace],
    pages    = [
        "Home" => "index.md",
        "Science" => [
            "sci/normalize.md",
            "sci/inference.md",
            "sci/autoencode.md",
        ],
        "Library" => [
            "lib/distance.md",
            "lib/generate.md",
            "lib/infer.md",
            "lib/io.md",
            "lib/manifold.md",
            "lib/model.md",
            "lib/normalize.md",
            "lib/pointcloud.md",
            "lib/queue.md",
            "lib/rank.md",
            "lib/scrna.md",
            "lib/util.md",
            "lib/voronoi.md",
        ],
        "Command Line" => [
            "cli/drosophila.md",
        ],
    ]
)

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
"""
 Tst
"""
function convertstringarray(s)
    if !ismissing(s)
        m = match(r"(?<=\[).*(?=\])", s)
        if isnothing(m)
            println(s)
        end
        m = replace(m.match, r"\s" => "")
        m = split(m, ",")
        return convert.(String, m)
    else
        return []
    end
end


text:
using TSXPlaygroundBenchmarks
SUITE = TSXPlaygroundBenchmarks.setup()


text:
using OrdinaryDiffEq, DiffEqDevTools, Test, Random
testTol = 0.2

Random.seed!(100)

# Test that the infrustructure works

f1 = (u,p,t) -> 2u
f2 = (u,p,t) -> 2u

prob = SplitODEProblem(f1,f2,1.0,(0.0,1.0))
sol = solve(prob,SplitEuler(),dt=1/10)
sol2 = solve(prob,Euler(),dt=1/10)
@test sol2[end] == sol[end]
@test sol2(0.345) == sol(0.345)

f3 = (u,p,t) -> 4u
prob2 = ODEProblem(f3,1.0,(0.0,1.0))
sol3 = solve(prob2,Euler(),dt=1/10)
@test sol3[end] == sol[end]
@test sol3(0.345) == sol(0.345)

u = rand(4,2)
f1 = (du,u,p,t) -> du.=2u
f2 = (du,u,p,t) -> du.=2u
prob = SplitODEProblem(f1,f2,u,(0.0,1.0))
sol = solve(prob,SplitEuler(),dt=1/10)
sol2 = solve(prob,Euler(),dt=1/10)

@test sol2[end] == sol[end]
@test sol2(0.345) == sol(0.345)

f3 = (du,u,p,t) -> du.=4u
prob2 = ODEProblem(f3,u,(0.0,1.0))
sol3 = solve(prob2,Euler(),dt=1/10)

@test sol3[end] == sol[end]
@test sol3(0.345) == sol(0.345)

println("OOP Tests")
# Now test only the first part
println("Testing only first part of Split ODE")

f1 = (u,p,t) -> 2u
f2 = (u,p,t) -> zero(u)

ff_split = SplitFunction(f1, f2; analytic=(u0,p,t)->exp(2t)*u0)
prob = SplitODEProblem(ff_split,1.0,(0.0,1.0))

sol = solve(prob,KenCarp3())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CFNLIRK3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

sol = solve(prob,KenCarp3(nlsolve=NLFunctional()))
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

sol = solve(prob,KenCarp4())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp4())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp5())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp5())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

sol = solve(prob,KenCarp47())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp47())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp58())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp58())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

# IMEXEuler
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEuler())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# IMEXEulerARK
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEulerARK())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# CNAB2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNAB2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# CNLF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNLF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF3
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF3())
@test_broken sim.𝒪est[:l∞] ≈ 3 atol=testTol

# SBDF4
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF4())
@test_broken sim.𝒪est[:l∞] ≈ 4 atol=testTol

# IRKC
dts = 1 .//2 .^(12:-1:8)
sim = test_convergence(dts,prob,IRKC())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# Now test only the second part
println("Testing only second part of Split ODE")

f1 = (u,p,t) -> zero(u)
f2 = (u,p,t) -> 2u

ff_split2 = SplitFunction(f1, f2; analytic=(u0,p,t)->exp(2t)*u0)
prob = SplitODEProblem(ff_split2,1.0,(0.0,1.0))

sol = solve(prob,KenCarp3())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CFNLIRK3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

sol = solve(prob,KenCarp4())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp4())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp5())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp5())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

sol = solve(prob,KenCarp47())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp47())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp58())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp58())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

# IMEXEuler
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEuler())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# IMEXEulerARK
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEulerARK())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# CNAB2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNAB2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# CNLF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNLF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF3
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF3())
@test_broken sim.𝒪est[:l∞] ≈ 3 atol=testTol

# SBDF4
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF4())
@test_broken sim.𝒪est[:l∞] ≈ 4 atol=testTol

# IRKC
dts = 1 .//2 .^(12:-1:8)
sim = test_convergence(dts,prob,IRKC())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# Test together
println("Testing both parts of Split ODE together")

f1 = (u,p,t) -> u
f2 = (u,p,t) -> 2u

ff_split3 = SplitFunction(f1, f2; analytic=(u0,p,t)->exp(3t)*u0)
prob = SplitODEProblem(ff_split3,1.0,(0.0,1.0))

sol = solve(prob,KenCarp3())
dts = 1 .//2 .^(12:-1:8)
sim = test_convergence(dts,prob,KenCarp3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

dts = 1 .//2 .^(12:-1:8)
sim = test_convergence(dts,prob,CFNLIRK3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

sol = solve(prob,KenCarp4())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp4())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp5())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp5())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

sol = solve(prob,KenCarp47())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp47())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp58())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp58())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

# IMEXEuler
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEuler())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# IMEXEulerARK
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEulerARK())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# CNAB2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNAB2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# CNLF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNLF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF3
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF3())
@test_broken sim.𝒪est[:l∞] ≈ 3 atol=testTol

# SBDF4
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF4())
@test_broken sim.𝒪est[:l∞] ≈ 4 atol=testTol

# IRKC
dts = 1 .//2 .^(12:-1:8)
sim = test_convergence(dts,prob,IRKC())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

println("IIP Tests")
# Now test only the first part
println("Testing only first part of Split ODE")

f1 = (du,u,p,t) -> du .= 2u
f2 = (du,u,p,t) -> du .= 0.0

ff_split4 = SplitFunction(f1, f2; analytic=(u0,p,t)->exp(2t)*u0)
prob = SplitODEProblem(ff_split4,rand(4,2),(0.0,1.0))

sol = solve(prob,KenCarp3())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CFNLIRK3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

sol = solve(prob,KenCarp4())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp4())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp5())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp5())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

sol = solve(prob,KenCarp47())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp47())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp58())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp58())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

# IMEXEuler
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEuler())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# IMEXEulerARK
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEulerARK())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# CNAB2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNAB2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# CNLF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNLF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF3
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF3())
@test_broken abs(sim.𝒪est[:l∞]-3) < testTol

# SBDF4
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF4())
@test_broken sim.𝒪est[:l∞] ≈ 4 atol=testTol

# IRKC
dts = 1 .//2 .^(12:-1:8)
sim = test_convergence(dts,prob,IRKC())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# Now test only the second part
println("Testing only second part of Split ODE")

f1 = (du,u,p,t) -> du.= 0.0
f2 = (du,u,p,t) -> du.= 2u

ff_split5 = SplitFunction(f1, f2; analytic=(u0,p,t)->exp(2t)*u0)
prob = SplitODEProblem(ff_split5,rand(4,2),(0.0,1.0))

sol = solve(prob,KenCarp3())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CFNLIRK3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

sol = solve(prob,KenCarp4())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp4())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp5())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp5())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

sol = solve(prob,KenCarp47())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp47())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp58())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp58())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

# IMEXEuler
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEuler())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# IMEXEulerARK
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEulerARK())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# CNAB2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNAB2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# CNLF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNLF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF3
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF3())
@test_broken sim.𝒪est[:l∞] ≈ 3 atol=testTol

# SBDF4
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF4())
@test_broken sim.𝒪est[:l∞] ≈ 4 atol=testTol

# IRKC
dts = 1 .//2 .^(12:-1:8)
sim = test_convergence(dts,prob,IRKC())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# Test together
println("Testing both parts of Split ODE together")

f1 = (du,u,p,t) -> du .= u
f2 = (du,u,p,t) -> du .= 2u

ff_split6 = SplitFunction(f1, f2; analytic=(u0,p,t)->exp(3t)*u0)
prob = SplitODEProblem(ff_split6,rand(4,2),(0.0,1.0))

sol = solve(prob,KenCarp3())
dts = 1 .//2 .^(12:-1:8)
sim = test_convergence(dts,prob,KenCarp3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

dts = 1 .//2 .^(12:-1:8)
sim = test_convergence(dts,prob,CFNLIRK3())
@test sim.𝒪est[:l∞] ≈ 3 atol=testTol

sol = solve(prob,KenCarp4())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp4())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp5())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp5())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

sol = solve(prob,KenCarp47())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp47())
@test sim.𝒪est[:l∞] ≈ 4 atol=testTol

sol = solve(prob,KenCarp58())
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,KenCarp58())
@test sim.𝒪est[:l∞] ≈ 5 atol=testTol

# IMEXEuler
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEuler())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# IMEXEulerARK
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,IMEXEulerARK())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol

# CNAB2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNAB2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# CNLF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,CNLF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF2
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF2())
@test sim.𝒪est[:l∞] ≈ 2 atol=testTol

# SBDF3
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF3())
@test_broken sim.𝒪est[:l∞] ≈ 3 atol=testTol

# SBDF4
dts = 1 .//2 .^(8:-1:4)
sim = test_convergence(dts,prob,SBDF4())
@test_broken sim.𝒪est[:l∞] ≈ 4 atol=testTol

# IRKC
dts = 1 .//2 .^(12:-1:8)
sim = test_convergence(dts,prob,IRKC())
@test sim.𝒪est[:l∞] ≈ 1 atol=testTol


text:
using PyCall
pushfirst!(PyVector(pyimport("sys")["path"]), "", "./core/", "./src/core/")
@pyimport plotwindow

module Context
    const pi = Base.pi
    const e = Base.ℯ
    const oo = Inf

    help() = Main.help()
    file() = Main.file()

    function get(x::Symbol)
        val = try
            Base.eval(Context, x)
        catch e
            if isa(e, UndefVarError)
                error("variable $(e.var) not defined.")
            end
        end
        val
    end

    function set(x::Main.Symbolic, value)
        if x ∈ [:pi, :e, :oo]
            error("$(x) is a constant and cannot be redefined.")
        end
        thunk = Main.createThunk(value)
        eval( :($x = $thunk) )
    end

    eval!(x) = error("eval! can only accept variables (you passed: $(typeof(x)))")
    function eval!(x::Symbol)
        value = eval( Context.get(x) )
        set(x, value)
    end

    eval(int::Main.Integral) = Main.evalIntegral(int)
    eval(t::Main.Thunk) = eval(t.expr)

    function eval(x::Symbol)
        val = try
            Base.eval(Context, x)
        catch e
            if isa(e, UndefVarError)
                error("variable $(e.var) not defined.")
            end
        end
        val
    end

    function eval(expr::Expr)
        if expr.head == :call && expr.args[1] == :eval!
            return eval!(expr)
        end

        eval_expr = copy(expr)
        if expr.head == :call
            if expr.args[1] == :plot
                return plot(expr.args[2:end]...)
            end
            evalThunk(x) = x
            evalThunk(t::Main.Thunk) = eval(t)
            args = map(evalThunk ∘ eval, expr.args[2:end])
            eval_expr.args = [expr.args[1], args...]
        end

        Base.eval(Context, eval_expr)
    end

    function plot(args...)
        n = length(args)
        if (n < 3 || n % 3 != 0)
            error("plot() requires a number of arguments divisible by 3 (you passed $n)")
        end

        plots = []
        data = [args[i:i+2] for i in 1:3:n]
        for (expr, a, b) in data
            a, b = eval(a), eval(b)
            step = (b - a) / 100
            xs, ys = @eval begin
                local func(x) = $expr
                local xs, ys = [], []
                for x in $a:$step:$b
                    try
                        y = func(x)
                        push!(xs, x)
                        push!(ys, y)
                    catch e
                        !isa(e, DomainError) && rethrow()
                    end
                end
                @assert length(xs) == length(ys)
                xs, ys
            end
            push!(plots, (xs, ys, "y = $expr"))
        end

        Main.plotwindow.createPlot(plots)
    end

    # Evals all eager parts of expression
    eagerCalc(x) = x
    function eagerCalc(ex::Expr)
        if ex.head == :call && length(ex.args) == 2
            func, arg = ex.args[1], ex.args[2]
            if func == :eval
                return eval(arg)
            elseif func == :eval!
                return eval!(arg)
            end
        end

        Expr(ex.head, map(eagerCalc, ex.args)...)
    end
end


text:
using DiffEqJump, DiffEqBase
# using BenchmarkTools
using Test, Graphs

Nsims        = 100
reltol       = 0.05
non_spatial_mean = [65.7395, 65.7395, 434.2605] #mean of 10,000 simulations

dim = 1
linear_size = 5
dims = Tuple(repeat([linear_size], dim))
num_nodes = prod(dims)
starting_site = trunc(Int,(linear_size^dim + 1)/2)
u0 = [500,500,0]
end_time = 10.0
diffusivity = 1.0

domain_size = 1.0 #μ-meter
mesh_size = domain_size/linear_size
# ABC model A + B <--> C
num_species = 3
reactstoch = [[1 => 1, 2 => 1],[3 => 1]]
netstoch = [[1 => -1, 2 => -1, 3 => 1],[1 => 1, 2 => 1, 3 => -1]]
rates = [0.1/mesh_size, 1.]
majumps = MassActionJump(rates, reactstoch, netstoch)

# spatial system setup
hopping_rate = diffusivity * (linear_size/domain_size)^2

# Starting state setup
starting_state = zeros(Int, length(u0), num_nodes)
starting_state[:,starting_site] .= u0

tspan = (0.0, end_time)
prob = DiscreteProblem(starting_state, tspan, rates)
hopping_constants = [hopping_rate for i in starting_state]
# algs = [NSM(), DirectCRDirect()]

function get_mean_end_state(jump_prob, Nsims)
    end_state = zeros(size(jump_prob.prob.u0))
    for i in 1:Nsims
        sol = solve(jump_prob, SSAStepper())
        end_state .+= sol.u[end]
    end
    end_state/Nsims
end

# testing
grids = [CartesianGridRej(dims), CartesianGridIter(dims), Graphs.grid(dims)]
jump_problems = JumpProblem[JumpProblem(prob, NSM(), majumps, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(false,false)) for grid in grids]
push!(jump_problems, JumpProblem(prob, DirectCRDirect(), majumps, hopping_constants=hopping_constants, spatial_system = grids[1], save_positions=(false,false)))
# setup flattenned jump prob
push!(jump_problems, JumpProblem(prob, NRM(), majumps, hopping_constants=hopping_constants, spatial_system = grids[1], save_positions=(false,false)))
# test
for spatial_jump_prob in jump_problems
    solution = solve(spatial_jump_prob, SSAStepper())
    mean_end_state = get_mean_end_state(spatial_jump_prob, Nsims)
    mean_end_state = reshape(mean_end_state, num_species, num_nodes)
    diff =  sum(mean_end_state, dims = 2) - non_spatial_mean
    for (i,d) in enumerate(diff)
        @test abs(d) < reltol*non_spatial_mean[i]
    end
end

#using non-spatial SSAs to get the mean
# non_spatial_rates = [0.1,1.0]
# reactstoch = [[1 => 1, 2 => 1],[3 => 1]]
# netstoch = [[1 => -1, 2 => -1, 3 => 1],[1 => 1, 2 => 1, 3 => -1]]
# majumps = MassActionJump(non_spatial_rates, reactstoch, netstoch)
# non_spatial_prob = DiscreteProblem(u0,(0.0,end_time), non_spatial_rates)
# jump_prob = JumpProblem(non_spatial_prob, Direct(), majumps)
# non_spatial_mean = get_mean_end_state(jump_prob, 10000)


text:
# SPDX-License-Identifier: X11
# 2020-12-07
# Day 1, Part 2

function main()
  input = parse.(Int, readlines())
  n = length(input)
  nums = falses(2020)

  for x ∈ input
    nums[x] = true
  end

  for i ∈ 1:n - 1
    for j ∈ i + 1:n
      x = 2020 - input[i] - input[j]
      if x > 0 && nums[x]
        println(input[i] * input[j] * x)
        return
      end
    end
  end
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
using Test
using NeXLParticle
using StaticArrays
using CoordinateTransformations

@testset "Translate" begin
    oldPts = [ SA[ 10.0 * rand(Float64,2)...] for _ in 1:20 ]
    θ, sc, off = π/3, SA[ 1.1 0.0; 0.0 0.9], SA[10.0,-10.0]
    tr = AffineMap( SA[ cos(θ) -sin(θ) ; sin(θ) cos(θ) ] * sc, off)
    newPts = (x->SA[ (tr(x) .+ 0.1*rand(2))...]).(oldPts)

    trp = translate(oldPts, newPts)
    diff = inv(trp).(newPts) .- oldPts
    for d in diff
        @test abs(d[1])<0.1
        @test abs(d[2])<0.1
    end
    diff = trp.(oldPts) .- newPts
    for d in diff
        @test abs(d[1])<0.1
        @test abs(d[2])<0.1
    end
end


text:
module WOA

using NCDatasets, Glob, CSV, Missings

if (pwd() in LOAD_PATH) == false
    push!(LOAD_PATH, pwd());
end

#import C2PO: gc_distance, rad2deg, deg2rad, histc, meshgrid, nan, findNaNmin, findNaNmax # this module stores additional function needed for analysis (e.g. gc_distance)

export npp_load_woa, load_woa, season, Profile, WOAdata

include("NPP/npp_path_setup.jl")

mutable struct Profile
    z::Array{Float64}
    temp::Array{Float64}
    salt::Array{Float64}
end

mutable struct WOAdata
    lat::Array{Float64}
    lon::Array{Float64}
    z::Array{Float64}
    temp::Array{Float64}
    salt::Array{Float64}
end

# 1: winter, 2: spring, 3: summer, 4: fall
const season = Dict([(1,"Winter"), (2,"Spring"), (3,"Summer"), (4,"Fall")]);
const seasontempfile = Dict([(1,"woa18_decav81B0_t13_04.nc"), (2,"woa18_decav81B0_t14_04.nc"), (3,"woa18_decav81B0_t15_04.nc"), (4,"woa18_decav81B0_t16_04.nc")]);
const seasonsaltfile = Dict([(1,"woa18_decav81B0_s13_04.nc"), (2,"woa18_decav81B0_s14_04.nc"), (3,"woa18_decav81B0_s15_04.nc"), (4,"woa18_decav81B0_s16_04.nc")]);

# this function loads the global WOA dataset, not necessary (overkill) for comparison with NPP CTD data.
function load_woa(i::Int64 = 3, woadir::String=string(ENV["HOME"], "/Research/data/WOA18/"))
    #seasonlocal = season;
    println(season[i])

    temproot = "https://data.nodc.noaa.gov/thredds/dodsC/ncei/woa/temperature/decav81B0/0.25/";
    saltroot = "https://data.nodc.noaa.gov/thredds/dodsC/ncei/woa/salinity/decav81B0/0.25/";

    #woatemppath = Glob.glob("woa*81B0_t*.nc", woadir);
    #woasaltpath = Glob.glob("woa*81B0_s*.nc", woadir);
    woatemppath = temproot * seasontempfile[i];
    woasaltpath = saltroot * seasonsaltfile[i];

    # load WOA18 location data for the specified season
    dstemp = Dataset(woatemppath);
    dssalt = Dataset(woasaltpath);
    latnc = dstemp["lat"];
    lonnc = dstemp["lon"];
    depthnc = dstemp["depth"];
    lon = Float64.(collect(Missings.replace(lonnc[:],NaN)));
    lat = Float64.(collect(Missings.replace(latnc[:],NaN)));
    z = -Float64.(collect(Missings.replace(depthnc[:],NaN)));
    #(LON,LAT) = meshgrid(lon,lat);

    tempnc = dstemp["t_an"];
    saltnc = dssalt["s_an"];
    temp = Float64.(collect(Missings.replace(tempnc[:],NaN)));
    salt = Float64.(collect(Missings.replace(saltnc[:],NaN)));
    woa = WOAdata(lat, lon, z, temp, salt);
    #return lat, lon, z, temp, salt;
    return woa;
end

# this function load the WOA data only at the closest WOA location for each NPP CTD station (requires internet connection)
function npp_load_woa(i::Int64 = 3, ctdlocationpath::String = ctdlocationpath)
    # load WOA data nearest to NPP CTD station for a specific season
    # i = 1: winter, 2: spring, 3: summer, 4: fall

    #woadir = string(ENV["HOME"], "/Research/data/WOA18/");
    #ctdlocationpath = string(ENV["HOME"], "/Research/NPP/data/CTD/NPP_ctd_locations.csv")

    #woa = load_woa(i,woadir); # only use this if want to load the global WOA dataset.
    temproot = "https://data.nodc.noaa.gov/thredds/dodsC/ncei/woa/temperature/decav81B0/0.25/";
    saltroot = "https://data.nodc.noaa.gov/thredds/dodsC/ncei/woa/salinity/decav81B0/0.25/";

    #woatemppath = Glob.glob("woa*81B0_t*.nc", woadir); # use this if a version of the WOA data is downloaded locally
    #woasaltpath = Glob.glob("woa*81B0_s*.nc", woadir);
    woatemppath = temproot * seasontempfile[i];
    woasaltpath = saltroot * seasonsaltfile[i];

    # load WOA18 location data for the specified season
    dstemp = Dataset(woatemppath);
    dssalt = Dataset(woasaltpath);
    latnc = dstemp["lat"];
    lonnc = dstemp["lon"];
    depthnc = dstemp["depth"];
    lon = Float64.(collect(Missings.replace(lonnc[:],NaN)));
    lat = Float64.(collect(Missings.replace(latnc[:],NaN)));
    z = -Float64.(collect(Missings.replace(depthnc[:],NaN)));
    #(LON,LAT) = meshgrid(lon,lat);

    tempnc = dstemp["t_an"];
    saltnc = dssalt["s_an"];
    #temp = Float64.(collect(Missings.replace(tempnc[:],NaN)));
    #salt = Float64.(collect(Missings.replace(saltnc[:],NaN)));
    #woa = WOAdata(lat, lon, z, temp, salt);

    #woa = load_woa(i, woadir);

    # location NPP CTD location data
    ctdloc = CSV.read(ctdlocationpath); # reading in the NPP CTD location data
    lonctd = convert(Vector,ctdloc.lon); # converting from dataframe to vector
    latctd = convert(Vector,ctdloc.lat);

    prof = Profile[];

    for j = 1:length(lonctd)
        println("Loading Station: " * string(j))
        #latind = findmin(abs.(woa.lat .- latctd[j]))[2];
        #lonind = findmin(abs.(woa.lon .- lonctd[j]))[2];;
        #push!(prof,Profile(woa.z, woa.temp[lonind,latind,:,1], woa.salt[lonind,latind,:,1]));
        latind = findmin(abs.(lat .- latctd[j]))[2];
        lonind = findmin(abs.(lon .- lonctd[j]))[2];

        temp = collect(Missings.replace(tempnc[lonind,latind,:,1],NaN));
        salt = collect(Missings.replace(saltnc[lonind,latind,:,1],NaN));
        push!(prof, Profile(z, temp, salt));
    end
    return prof
end

end #module


text:
module OPpoly
#
# Module OPpoly contains functions manipulating multivariate polynomials.
# A polynomial is represented as a vector of coefficients.  The indexing
# is arranged (taking a D=3,ord=4 example) as
# c000, c100, c010, c001, c200, c110, c101, c020, c011, c002, c300, c210,
# c201, c120, c111, c102, c030, c021, c012, c003, c400, c310, c301, c220,
# c211, c202, c130, c121, c112, c103, c040, c031, c022, c013, c004
#

export tensorToList, listToTensor, polyrec, polyeval, rateval,
       solveLSPolyCoeffs

# =======================================================================
function tensorToList(t)
#
# Convert between tensor (00,10,01,20,11,02,30,21,12,03,...) and list
# (1,2,3,...) indexing.
#
# Version:        Changes:
# --------        -------------
# 13.04.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 13.04.2021      Tested on sample problems up to D=4 and ord=3.
#
# Inputs:
# -------
# t               : D-by-N matrix of integers, with D the dimension.
#
# Outputs:
# --------
# L               : Length N vector of integers, containing the indices.

   D = size(t,1)
   N = size(t,2)

   L = ones(Int,N)
   for ip in 1:N
      ord = sum(t[1:D,ip])
      for id in 1:D
         De = D - id + 1
         L[ip] = L[ip] + binomial((ord-1)+De,De)
         ord = ord - t[id,ip]
      end
   end

   return L

end # tensorToList

# =======================================================================
function listToTensor(L, D)
#
# Convert between list (1,2,3,...) and tensor indexing
# (00,10,01,20,11,02,30,21,12,03,...).
#
# Version:        Changes:
# --------        -------------
# 13.04.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 13.04.2021      Verified that this inverts tensorToList.
#
# Inputs:
# -------
# L               : Index list, vector of length N.
# D               : Dimension, integer.
#
# Outputs:
# --------
# t               : D-by-N matrix of tensor indices.

   N = length(L)

   # The matrix A is consistent for any dimension D.
   A = zeros(Int,D,D)
   for id = 1:D
      A[id,id:D] = ones(Int,1,D-id+1)
   end

   b = zeros(Int,D,N)
   for ip in 1:N

      Le = L[ip]

      for id in 1:D

         De = D - id + 1

         ind = 1
         sav = 0
         flg = true
         while flg
            if (ind >= Le)
               flg = false
               Le = Le - sav
            else
               b[id,ip] = b[id,ip] + 1
               sav = ind
               ind = binomial(De+b[id,ip],De)
            end
         end

      end

   end

   t = zeros(Int,D,N)
   t[:,:] = A\b

   return t

end # listToTensor

# =======================================================================
function polyrec(x, ord; ptype=0)
#
# Recursively evaluate the normalized terms of a polynomial; that is,
# the value of each term when the coefficient is equal to 1.
#
# Version:        Changes:
# --------        -------------
# 22.04.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 22.04.2021      Hand-checked the first few terms of a D=2 case.
#
# Inputs:
# -------
# x               : Column vector of length D.
# ord             : The order of individual x components, up to which
#                   the polynomial is to be computed.
# ptype           : Specifies the type of polynomial.
#                   0: Nominal   x[n+1] = x x[n]
#                   1: Chebyshev T[n+1] = 2xT[n] - T[n-1]
#                   2: Legendre  P[n+1] = ((2n+1)xP[n] - nP[n-1])/(n+1)
#                   3: Hermite   H[n+1] = (x H[n] - sq(n)H[n-1])/sq(n+1)
#
# Outputs:
# --------
# psi             : A D-by-ord+1 matrix containing the terms' values
#                   at the given x vector.

   D = size(x,1)
   psi = zeros(D,ord+1)
   psi[:,1] = ones(D)

   if (ord > 0)
      if (ptype == 1)
         # Chebyshev
         for id in 1:D
            psi[id,2] = x[id]
            for n in 1:ord-1
               psi[id,n+2] = 2*x[id]*psi[id,n+1] - psi[id,n]
            end
         end
      elseif (ptype == 2)
         # Legendre
         for id in 1:D
            psi[id,2] = x[id]
            for n in 1:ord-1
               psi[id,n+2] = ((2*n+1)*x[id]*psi[id,n+1] - n*psi[id,n])/(n+1)
            end
         end
      elseif (ptype == 3)
         # Hermite
         for id in 1:D
            psi[id,2] = x[id]
            sqn = 1.0
            for n in 1:ord-1
               sqn1 = sqrt(n+1)
               psi[id,n+2] = (x[id]*psi[id,n+1] - sqn*psi[id,n])/sqn1
               sqn = sqn1  # Avoid recomputing sqrt's.
            end
         end
      else
         # Default to nominal.
         for id in 1:D
            psi[id,2] = x[id]
            for n = 1:ord-1
               psi[id,n+2] = x[id]*psi[id,n+1]
            end
         end
      end
   end # if ord > 0.

   return psi

end # polyrec

# =======================================================================
function polyeval(p, t, x; ptype=0)
#
# Evaluate a polynomial defined by coefficients p of the tensor product
# t at the coordinates x.  This is currently designed to be easy to use,
# not optimal in terms of speed and storage.
#
# Version:        Changes:
# --------        -------------
# 22.04.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 22.04.2021
#
# Inputs:
# -------
# p               : Vector of polynomial coefficients.  See the module
#                   header for the convention used for indexing.
# t               : Tensor specifying the x component orders in each
#                   polynomial term.  D-by-N matrix.
# x               : Column vector or matrix; each column contains an x
#                   vector of length D.
# ptype           : Specifies the type of polynomial.  See polyrec.
#
# Outputs:
# --------
# y               : A vector containing the outputs, one scalar output
#                   for each x vector.

   D = size(x,1)   # Dimension of x.
   M = size(x,2)   # Number of points.
   N = size(p,1)   # Number of coefficients.

   ord = maximum(sum(t,dims=1))

   y = zeros(M)

   for ip in 1:M

      # Compute the individual polynomial terms in a recursive way.
      if (D == 1)
         # This is required since the ":" notation is not permitted
         # if x is a scalar.
         psi = polyrec(x[ip],ord,ptype=ptype)
      else
         psi = polyrec(x[:,ip],ord,ptype=ptype)
      end

      for ic in 1:N
         val = 1.0
         for id = 1:D
            val = val*psi[id,t[id,ic]+1]
         end
         y[ip] = y[ip] + p[ic]*val
      end

   end

   return y

end # polyeval: version including explicit input of the sequence.

function polyeval(p, x; ptype=0)
#
# Evaluate a polynomial defined by coefficients p at the coordinates x.
# This is currently designed to be easy to use, not optimal in terms of
# speed and storage.
#
# This version generates the tensor t automatically, based on the
# assumed sequence of polynomials indicated in the module header.
#
# Version:        Changes:
# --------        -------------
# 22.04.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 22.04.2021
#

   D = size(x,1)   # Dimension of x.
   N = size(p,1)   # Number of coefficients.

   L = range(1,stop=N)
   t = listToTensor(L,D)

   y = polyeval(p,t,x,ptype=ptype)

   return y

end # polyeval: version with an assumed polynomial sequence.

# =======================================================================
function rateval(p, tp, q, tq, x; ptype=0)
#
# Evaluate a rational function defined by numerator coefficients p and
# denominator coefficients q at the coordinates x.  The order associated
# with each coefficient is given in the matrices tp,tq.
#
#
# Version:        Changes:
# --------        -------------
# 22.04.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 22.04.2021      Verified some example calculations.
#
# Inputs:
# -------
# p,q             : Coefficients of numerator and denominator polynomials.
# tp,tq           : Tensors specifying the x component orders in each
#                   polynomial term.  D-by-Np,Nq matrices.
# x               : Column vector or matrix; each column contains an x
#                   vector.
#
# Outputs:
# --------
# y               : A vector containing the outputs, one scalar output
#                   for each x vector.

   yp = polyeval(p,tp,x,ptype=ptype)
   yq = polyeval(q,tq,x,ptype=ptype)

   return (yp ./ yq)

end # rateval: explicit sequence.

function rateval(p, q, x; ptype=0)
#
# Evaluate a rational function defined by numerator coefficients p and
# denominator coefficients q at the coordinates x.
#
# Version:        Changes:
# --------        -------------
# 22.04.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 22.04.2021
#

   yp = polyeval(p,x,ptype=ptype)
   yq = polyeval(q,x,ptype=ptype)

   return (yp ./ yq)

end # rateval: implicit sequence.

#=

# =======================================================================
function polyderiv(p, t, k, x; cof=0)
#
# Evaluate the derivatives of a polynomial at specified points.  At the
# moment this supports only nominal polynomials (ptype = 0).  That is,
# if a column of t contains, say, (1, 0, 3, 2) then the polynomial term
# is (x1 x3^3 x4^2).
#
# Version:        Changes:
# --------        -------------
# 18.05.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 18.05.2021
#
# Inputs:
# -------
# p               : Vector of polynomial coefficients, length N.  See
#                   the module header for the convention used for
#                   indexing.
# t               : Tensor specifying the x component orders in each
#                   polynomial term.  D-by-N matrix.
# k               : An integer specifying which derivative to take:
#                   1 = x1, 2 = x2, ...
# x               : Column vector or matrix of dimension D-by-M; each
#                   column contains an x vector of length D.
#
# Outputs:
# --------
# pd, td          : Coefficients and tensor describing the derivative
#                   polynomial.

   ptype = 0       # Limited functionality in order to get something
                   # working fast.

   D = size(x,1)   # Dimension of x.
   M = size(x,2)   # Number of points.
   N = size(p,1)   # Number of coefficients.

   ord = maximum(sum(t,dims=1))

   dy = zeros(M,L)

   for ip in 1:M

      # Compute the individual polynomial terms in a recursive way.
      if (D == 1)
         # This is required since the ":" notation is not permitted
         # if x is a scalar.
         psi = polyrec(x[ip],ord,ptype=ptype)
      else
         psi = polyrec(x[:,ip],ord,ptype=ptype)
      end

      for ic in 1:N
         val = 1.0
         for id in 1:D
            val = val*psi[id,t[id,ic]+1]
         end
         y[ip] = y[ip] + p[ic]*val
      end

   end

   return y

end # polyderiv

=#

# =======================================================================
#function polyint()
#
# Evaluate
#
# Version:        Changes:
# --------        -------------
# 13.04.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 13.04.2021
#
# Inputs:
# -------
#
#
# Outputs:
# --------
#

#end # polyint

# ===================================================================
function solveLSPolyCoeffs(x, y, tp, tq; ptype=0, wgt=1.0)
#
# Given a list of samples, solves a (weighted) least-squares
# polynomial or rational-function fit.
#
# [Note, this should be augmented to accommodate complex coefficients.]
#
# Version:        Changes:
# --------        -------------
# 22.04.2021      Original code.
#
# Version:        Verification:
# --------        -------------
# 22.04.2021      Verified the fit to a known quadratic rational
#                 polynomial function in two dimensions.
#
# Inputs:
# -------
# x               : D-by-Ns vector of coordinates.
# y               : Length Ns vector of sample values.
# tp,tq           : Tensors specifying the x component orders in each
#                   polynomial term.  D-by-Np,Nq matrices.  Note that
#                   tq should always include one all-zeros term as its
#                   first entry.
# ptype           : Specifies the type of polynomial.  See polyrec.
# wgt             : A vector of weights, the square roots of the weight
#                   matrix diagonal.
#
# Outputs:
# --------
# p,q             : Vector of best-fit polynomial coefficients; p
#                   for the numerator polynomials, q for the
#                   denominator polynomials.

   D  = size(x,1)
   Ns = size(x,2)
   Np = size(tp,2)
   Nq = size(tq,2)
   Nc = Np + Nq

   # The maximum polynomial order in the individual variables.
   op = maximum(tp)
   oq = maximum(tq)
   omax = max(op,oq)

   if (omax == 0)
      # We are dealing with fitting a least-squares constant, not a
      # polynomial.
      A = ones(Ns,1) .* wgt
      AT = transpose(A)
      p = (AT*A)\(AT*(wgt .* y))
      q = 1.0
   else

      # The matrix in the equation A [p;q] = y.  The "-1" is for the
      # forced constant = 1 term in the denominator.
      A = ones(Ns,Nc-1)
      for isamp in 1:Ns

         # Compute the individual polynomial terms in a recursive way.
         psi = polyrec(x[:,isamp],omax,ptype=ptype)

         for ip in 1:Np
            for id in 1:D
               A[isamp,ip] = A[isamp,ip]*psi[id,tp[id,ip]+1]
            end
         end
         for iq in 1:Nq-1  # Nq-1: Neglect the constant 1.0 in the denominator.
            ind = Np + iq
            for id in 1:D
               A[isamp,ind] = A[isamp,ind]*psi[id,tq[id,iq+1]+1]
            end
            A[isamp,ind] = -A[isamp,ind]*y[isamp]
         end
      end

      A = A .* wgt

      AT = A'
      c = (AT*A)\(AT*(wgt .* y))
      p = c[1:Np]
      if (oq > 0)
         q = [1.0;c[Np+1:Nc-1]]
      else
         q = 1.0
      end

   end

   return p, q

end  # solveLSPolyCoeffs

end # module

max_stars_repo_name: SINTEF-Energy-Wind/OutliersPrimer

max_issues_repo_name: SINTEF-Energy-Wind/OutliersPrimer

max_forks_repo_name: SINTEF-Energy-Wind/OutliersPrimer


text:
using LinearAlgebra
using Plots
using IterativeSolvers

function func(y,t,ϵ)
    return -(1/ϵ)*(y-sin(t));
end

function actual_solution(t,ϵ)
    return (sin(t) - ϵ*cos(t))/(1+ϵ^2) + (1 + ϵ/(1+ϵ^2))*exp(-t/ϵ);
end

function compute_solution(t,ϵ)
    n = length(t);
    y = zeros(n);
    for i in 1:n
       y[i] = actual_solution(t[i],ϵ);
    end
    return y
end

function explicit_euler(t,y0,ϵ)
    n = length(t);
    y = zeros(n);
    y[1] = y0;
    dt = t[2]-t[1];
    for i in 2:n
       y[i] = y[i-1] + dt*func(y[i-1],t[i-1],ϵ);
    end
    return y;
end

function sor_implicit_euler(t,y0,ϵ)
    n = length(t);
    y = zeros(n);
    y[1] = y0;
    dt = t[2]-t[1];
    res = 100;
    τ = 1e-6;
    iter = 1e6;
    κ2 = (dt/ϵ);
    κ1 = 1/(1 + κ2);
    κ3 = 1/κ1;
    it = 0;
    while (it<iter && res>τ)
        for i in 2:n
           y[i] = κ1*(y[i-1] + κ2*sin(t[i]));
        end
        for i in 2:n
           res+= κ2*sin(t[i]) + y[i-1] - κ3*y[i];
        end
        res = sqrt(abs(res/n^2));
        it += 1;
    end
    return y;
end

function sor_trapezoidal_method(t,y0,ϵ)
    n = length(t);
    y = zeros(n);
    y[1] = y0;
    dt = t[2]-t[1];
    res = 100;
    τ = 1e-6;
    iter = 1e6;
    κ1 = (dt/(2*ϵ));
    κ2 = 1 - κ1;
    κ3 = 1 + κ1;
    κ4 = 1/κ3;
    it = 0;
    while (it<iter && res>τ)
        for i in 2:n
            y[i] = κ4*(κ1*(sin(t[i]) + sin(t[i-1])) + κ2*y[i-1]);
        end
        for i in 2:n
            res += κ1*(sin(t[i]) + sin(t[i-1])) + κ2*y[i-1] - κ3*y[i];
        end
        res = sqrt(abs(res/n^2));
        it += 1;
    end
    return y;
end

function norm(y1,y2)
    n = length(y1);
    norm = 0;
    for i = 1:n
        norm += abs(y1[i]-y2[i]);
    end
    norm = sqrt(norm/n*n);
    println("The norm is ",norm)
end

x0 = 0;
xl = pi;
N = 99;
y0 = 1;
x = LinRange(x0,xl,N)

ϵ = 1;
y_a = compute_solution(x,ϵ);
y = sor_trapezoidal_method(x,y0,ϵ);
plot(x,y)
p1 = plot!(x,y_a,title="Symplectic ϵ=1")
norm(y,y_a);

ϵ = 0.1;
y_a = compute_solution(x,ϵ);
y = sor_trapezoidal_method(x,y0,ϵ);
plot(x,y)
p2 = plot!(x,y_a,title="Symplectic  ϵ=0.1")
norm(y,y_a);

ϵ = 0.01;
y_a = compute_solution(x,ϵ);
y = sor_trapezoidal_method(x,y0,ϵ);
plot(x,y)
p3 = plot!(x,y_a,title="Symplectic  ϵ=0.01")
norm(y,y_a);

ϵ = 0.001;
y_a = compute_solution(x,ϵ);
y = sor_trapezoidal_method(x,y0,ϵ);
plot(x,y)
p4 = plot!(x,y_a,title="Symplectic  ϵ=0.001")
norm(y,y_a);

plot(p1,p2,p3,p4,layout=(2,2),legend=false)

x0 = 0;
xl = pi;
N = 99;
y0 = 1;
x = LinRange(x0,xl,N)

ϵ = 1;
y_a = compute_solution(x,ϵ);
y = explicit_euler(x,y0,ϵ);
plot(x,y)
p1 = plot!(x,y_a,title="Explcit Euler ϵ=1")
norm(y,y_a);

ϵ = .1;
y_a = compute_solution(x,ϵ);
y = explicit_euler(x,y0,ϵ);
plot(x,y)
p2 = plot!(x,y_a,title="Explcit Euler ϵ=0.1")
norm(y,y_a);

ϵ = .01;
y_a = compute_solution(x,ϵ);
y = explicit_euler(x,y0,ϵ);
plot(x,y)
p3 = plot!(x,y_a,title="Explcit Euler ϵ=0.01")
norm(y,y_a);

ϵ = .001;
y_a = compute_solution(x,ϵ);
y = explicit_euler(x,y0,ϵ);
plot(x,y)
p4 = plot!(x,y_a,title="Explcit Euler ϵ=0.001")
norm(y,y_a);

plot(p1,p2,p3,p4,layout=(2,2),legend=false)

x0 = 0;
xl = pi;
N = 99;
y0 = 1;
x = LinRange(x0,xl,N);

ϵ = 1;
y_a = compute_solution(x,ϵ);
y = sor_implicit_euler(x,y0,ϵ);
plot(x,y)
p1 = plot!(x,y_a,title="Implicit Euler ϵ=1")
norm(y,y_a);

ϵ = 0.1;
y_a = compute_solution(x,ϵ);
y = sor_implicit_euler(x,y0,ϵ);
plot(x,y)
p2 = plot!(x,y_a,title="Implicit Euler ϵ=0.1")
norm(y,y_a);

ϵ = 0.01;
y_a = compute_solution(x,ϵ);
y = sor_implicit_euler(x,y0,ϵ);
plot(x,y)
p3 = plot!(x,y_a,title="Implicit Euler ϵ=0.01")
norm(y,y_a);

ϵ = 0.001;
y_a = compute_solution(x,ϵ);
y = sor_implicit_euler(x,y0,ϵ);
plot(x,y)
p4 = plot!(x,y_a,title="Implicit Euler ϵ=0.001")
norm(y,y_a);

plot(p1,p2,p3,p4,layout=(2,2),legend=false)

max_stars_repo_name: dynamic-queries/PDEs

max_issues_repo_name: dynamic-queries/PDEs

max_forks_repo_name: dynamic-queries/PDEs


text:
# MC analysis on the speed of sound

using Statistics

MAX = Dict{Int, Vector{Float64}}()
MEAN = Dict{Int, Vector{Float64}}()

for mc in 1:10, ss in 420:20:420 #270:20:410

    !haskey(MAX, ss) ? (MAX[ss]=Float64[]) : nothing
    !haskey(MEAN, ss) ? (MEAN[ss]=Float64[]) : nothing

    intensityLeft = beamformBasic(leftwav, lefttemplate, soundSpeed=ss)
    intensityRight = beamformBasic(rightwav, righttemplate, soundSpeed=ss)

    pcLeft = intensityToCircularKDE(intensityLeft, SNRfloor=0.8, N=200)
    pcRight = intensityToCircularKDE(intensityRight, SNRfloor=0.8, N=200)

    fg = initfg()

    # add two beacons
    addVariable!(fg, :l0, Point2, labels=["BEACON"])
    addFactor!(fg, [:l0], Prior(MvNormal(l0, lcov)))

    addVariable!(fg, :l1, Point2, labels=["BEACON"])
    addFactor!(fg, [:l1], Prior(MvNormal(l1, lcov)))

    # add unknown pose location
    addVariable!(fg, :x0, Pose2, labels=["POSE"])
    addFactor!(fg, [:x0], PartialPrior(Normal(0.0, 0.01), (2,)) )

    # add two bearing only measurements
    addFactor!(fg, [:x0; :l0], Pose2Point2BearingRange(pcLeft, Rayleigh(1.0)))
    addFactor!(fg, [:x0; :l1], Pose2Point2BearingRange(pcRight, Rayleigh(1.0)))

    batchSolve!(fg, N=100)

    pts = getVal(fg, :x0)
    pts0 = pts[:,pts[1,:] .< 1.0]
    X00 = kde!(pts0[1,:])

    push!(MAX[ss], getKDEMax(X00)[1])
    push!(MEAN[ss], getKDEMean(X00)[1])

end

@save "/home/dehann/.julia/dev/Caesar/examples/tracking/acoustic/mc_final.jld2" MAX MEAN

for mm in MAX
  @show mm[1], Statistics.mean(mm[2]), std(mm[2])
  @show mm[1], Statistics.mean(MEAN[mm[1]]), std(MEAN[mm[1]])
end

using DataFrames

DF = DataFrame[]

for mm in keys(MAX)
  push!(DF, DataFrame(
    speed = mm,
    error = MEAN[mm],
    Legend = "mean"
  ) )
  # push!(DF, DataFrame(
  #   speed = mm,
  #   error = MAX[mm],
  #   Legend = "Max"
  # ) )
end

alldata = vcat(DF...)

pl = Gadfly.plot(alldata, x=:speed, y=:error, Geom.boxplot)

pl |> SVG("/home/dehann/.julia/dev/Caesar/examples/tracking/acoustic/soundSpeed.svg",10cm,6cm)

## Linear regression

using GLM

alldata

ols = lm(@formula(error ~ speed), alldata)

round.(stderror(ols), digits=5)

gg = (x) -> 0.00231747*x - 0.815986

0.815986/0.00231747

using Colors

ll = Gadfly.layer(gg, 200, 450, Theme(default_color=colorant"red"))

push!(pl.layers, ll[1])

pl

#


text:
using JuMP
using Cbc
model = Model(Cbc.Optimizer)
# set_optimizer_attribute(model, "cuts","off")

month_num = 8
month = 1:month_num

# paramerters
q = [5000] # product fix cost
p = [100] # cost per product
d = [400, 400, 800, 800, 1200, 1200, 1200, 1200] # demand of product by month
s_init = [200] # initial stock
h = [5] # stock cost per product month

@variable(model, x[j in month] >= 0)
@variable(model, s[j in month] >= 0)
@variable(model, y[j in month], Bin)

# 库存平衡
@constraint(model ,dem_sat1, s_init[1] + x[1] == d[1] + s[1])
@constraint(model, dem_sat[j in month[2:end]], s[j - 1] + x[j] == d[j] + s[j] )

# 最大产量约束
@constraint(model, vub[j in month], x[j] <= sum(d[k]  for k in j:month_num) * y[j])

# 目标函数
@expression(model ,cost, sum(p[1] * x[j] + q[1] * y[j] for j in month) + sum(h[1] * s[j] for j in month[1:end - 1]))
@objective(model, Min, cost)
optimize!(model)

# 显示结果
cost_val = value.(cost)
x_val = value.(x)
s_val = value.(s)
@show cost_val;
@show x_val;
@show s_val;

max_stars_repo_name: Behemoth-s/ProductPlanningByMIP

max_issues_repo_name: Behemoth-s/ProductPlanningByMIP

max_forks_repo_name: Behemoth-s/ProductPlanningByMIP


text:
using Distributions
using DataStructures

include("sample_sobol.jl")
include("analyze_sobol.jl")
include("test_functions/ishigami.jl")

# define the data
data = SobolData(
    params = OrderedDict(:x1 => Uniform(-3.14159265359, 3.14159265359),
        :x2 => Uniform(-3.14159265359, 3.14159265359),
        :x3 => Uniform(-3.14159265359, 3.14159265359)),
    N = 1000
)

# generate samples using Sobol sequence
samples = sample(data)

# run model (example)
Y = ishigami(samples)

# perform Sobol Analysis
analyze(data, Y)


text:
import Base.copy

"""
    ImpurityQuantities

Contains all quantities of a given channel, computed by DMFT
"""
struct ImpurityQuantities
    Γ::SharedArray{Complex{Float64},3}
    χ::SharedArray{Complex{Float64},3}
    χ_ω::SharedArray{Complex{Float64},1}
    χ_loc::Complex{Float64}
    usable_ω::AbstractArray{Int,1}
    tailCoeffs::AbstractArray{Float64,1}
end

"""
    NonLocalQuantities

Contains all non local quantities computed by the lDGA code
"""
mutable struct NonLocalQuantities{T1 <: Union{Complex{Float64}, Float64}, T2 <: Union{Complex{Float64}, Float64}}
    χ::SharedArray{T1,2}
    γ::SharedArray{T2,3}
    usable_ω::AbstractArray
    λ::Float64
end

Base.copy(x::T) where T <: Union{NonLocalQuantities, ImpurityQuantities} = T([deepcopy(getfield(x, k)) for k ∈ fieldnames(T)]...)

const ΓT = SharedArray{Complex{Float64},3}
const BubbleT = SharedArray{Complex{Float64},3}
const GνqT = SharedArray{Complex{Float64},2}
const qGridT = Array{Tuple{Int64,Int64,Int64},1}


text:
@testset "c\"typedef ...\"" begin
	@eval c``

	@eval c"""
	typedef struct {
	} __attribute__((packed)) EmptyStructTypedef;
	"""
	@test sizeof(c"EmptyStructTypedef") == 0

	@eval c"""
	typedef struct {
		char c;
	} __attribute__((packed)) CcharStructTypedef;
	"""
	@test sizeof(c"CcharStructTypedef") == sizeof(Cchar)
	@test :c in fieldnames(c"CcharStructTypedef")

	@eval c"""
	typedef union {
		int i;
	} __attribute__((packed)) CintUnionTypedef;
	"""
	@test sizeof(c"CintUnionTypedef") == sizeof(Cint)
	@test :i in fieldnames(c"CintUnionTypedef")

	@eval c"""
	typedef enum {
		CENUM_TYPEDEF_1,
		CENUM_TYPEDEF_2 = 1 << 8,
		CENUM_TYPEDEF_3 = CENUM_TYPEDEF_2,
	} __attribute__((packed)) CenumTypedef;
	"""
	@test sizeof(c"CenumTypedef") == 2
	@test c"CENUM_TYPEDEF_3" != c"CENUM_TYPEDEF_1"
	@test c"CENUM_TYPEDEF_3" == c"CENUM_TYPEDEF_2"

	@eval c"""
	typedef unsigned int CuintArrayTypedef[32];
	"""
	@test sizeof(c"CuintArrayTypedef") == sizeof(Cuint)*32

	@eval c"""
	typedef int CintTypedef;
	"""
	@test sizeof(c"CintTypedef") == sizeof(Cint)
end


text:
# Random sub-space ensemble module
module RandomSubspace

	using j4pr.ClassifierCombiner
	using MLDataPattern: getobs, datasubset, ObsDim
	export SubspaceEnsemble, randomsubspace_train, randomsubspace_exec

	"""
	Random subspace ensemble object (it describes the ensemble).

	# Fields
	 * `idx::Vector{Vector{Int}}` is a vector of variable index vectors, 1 vector/ensemble member
	 * `f_train::Function` is the training function
	 * `f_exec::Function` is the execution function
	 * `members::Vector` is a vector containing the trained classifiers
	 * `combiner::ClassifierCombiner.AbstractCombiner` the means by which ensemble outputs are combined. Non-trainable combiners should be used. If a trainable
	    combiner is to be used, it has to be trained outside the ensemble
	* `parallel_train::Bool` specifies whether to parallelize training or not
	* `parallel_execution::Bool` specifies whether to parallelize execution or not

	"""
	struct SubspaceEnsemble{F, T<:Vector{S} where S, U<:AbstractCombiner}
		idx::Vector{Vector{Int}}	# vector of vectors of variable indices
		f_train::Function		# training function
		f_exec::F			# execution function
		members::T			# ensemble members (trained)
		combiner::U			# combiner (only non-trainable combiners supported)
		parallel_train::Bool		# whether to parallelize ensemble training or not
		parallel_execution::Bool	# whether to parallelize ensemble execution or not
	end

	# Aliases
	const SubspaceEnsembleCell{F<:Void, T, U} = SubspaceEnsemble{F,T,U}		# no execution function, it is dynamically generated in training
	const SubspaceEnsembleClassic{F<:Function, T, U} = SubspaceEnsemble{F,T,U} 	# the execution function is specified

	# Printers
	Base.show(io::IO, m::SubspaceEnsembleCell) = print(io, "Sub-space ensemble, $(length(m.idx)) members, Cell version")
	Base.show(io::IO, m::SubspaceEnsembleClassic) = print(io, "Sub-space ensemble, $(length(m.idx)) members")

	"""
		subspace_indices(n, L, M, replace)

	Function that generates a `Vector{Vector{Int}` containing the `L` random subsets of variable indices for a given dataset of `n` variables using `M` variables.
	If `replace` is `true`, the same variable can appear in multiple subsets and the number of subsets has no upper bound.
	"""
	function subspace_indices(n::Int, L::Int, M::Int, replace::Bool)
		@assert n > 1 "The dataset has to have more than 1 variable."
		@assert L > 0 "The random subset ensemble has to have at least 1 member."
		@assert M > 0 "The number of random variables has to be at least 1. "

		# Generate V based on n,L,M
		if replace
			V=Vector{Vector{Int}}(L)
			@inbounds for i in 1:L
				V[i] = randperm(n)[1:M]
			end
			return V
		else
			@assert L*M <=n "Cannot split $n variables to an ensemble of $L elements with $M distinct variables each."
			ro = randperm(n) # random order
			V=Vector{Vector{Int}}(L)
			k=1
			@inbounds for i in 1:M:L*M
				V[k] = ro[i:i+M-1]
				k+=1
			end
			return V
		end
	end

	"""
		randomsubspace_train(X, y, L, M, f_train, combiner, replace; parallel_train=true, parallel_execution=false)

	Trains a random subspace classifier. The function runs `f_train` on subspaces of `X` and returns `RandomSubspace` object.

	# Arguments
	 * `X::AbstractMatrix` is the training data; columns are assumed to be samples, lines variables
	 * `y::AbstractVector` are the labels
	 * `L::Int` is the ensemble size
	 * `M::Int` is the number of variables for each ensemble member
	 * `f_train::Function` is the training function for the ensemble members
	 * `f_exec::Union{Function,Void}` is the execution function for the ensemble members
	 * `combiner::ClassifierCombiner.AbstractCombiner` is the combiner
	 * `replace::Bool` specifies whether the variable subsets can contain (`true`) of not (`false`) the same variable. If `true`, any number of subsets can be generated; if
	 `false` a finite maximum number can be generated but all subsets will be distinct

	# Keyword arguments
	 * `parallel_train::Bool=false` specifies whether to paralellize or not the training. Generally a good ideea except for lazy classifiers (e.g. kNN).
	 * `parallel_execution::Bool=false` specifies whether to paralellize or not the execution. Unless the execution is time consuming (e.g. Parzen classifier), it should be left to `false`.

	It is important to note that the training method calls internally `x->f_train((x,y))` where `x` are the subsets of `X` constructed using `L`, `M` and the total number of variables of `X`.
	Considering a generic training function `g` desired to be used in the ensemble, of the following forms:
	 * `g(X, y, arg1, arg2)`, one should provide as argument something as `f_train = (x)->g_train(x[1], x[2], arg1, arg2) # uses labels`
	 * `g(X, arg1, arg2)`, one should provide as argument something as `f_train = (x)->g_train(x[1], arg1, arg2) # ignores labels`
	 * `g(X, arg1, arg2)` where `X==(data,labels)`, one should provide as argument something as `f_train = (x)->g_train(x, arg1, arg2) # MLDataUtils container case`
	 * `g(args...)::CellFunU`, one can simply provide `g.f` (i.e. the training function of the untrained function cell) since it already contains the training arguments.

	 During execution, `f_exec(ensemble.members[i], Xᵢ)` will be called for all subspace subsets `Xᵢ` of `X` corresponding to the members of the ensemble hence `f_exec` must have
	 a signature `f_exec(M, X)` where `M` is represents the learned model and `X` represents the data.
	"""
	function randomsubspace_train(X::T where T<:AbstractMatrix, y::S where S<:AbstractVector, L::Int, M::Int, f_train::Function, f_exec::Union{Function,Void},
			       		combiner::S where S<:AbstractCombiner, replace::Bool;
					parallel_train::Bool=true, parallel_execution::Bool=false)

		# Calculate the number of variables
		n = size(X,1)

		# Calculate subspace indices
		idx = subspace_indices(n, L, M, replace)

		# Return data subsets (Vector of datasubsets of X) and V
		if parallel_train
			members = pmap(x->f_train((x,y)), (getobs(datasubset(X,i,ObsDim.First())) for i in idx) )
		else
			members = map(x->f_train((x,y)), (getobs(datasubset(X,i,ObsDim.First())) for i in idx) )
		end

		return SubspaceEnsemble(idx, f_train, f_exec, members, combiner, parallel_train, parallel_execution)
	end

	"""
		randomsubspace_exec(ensemble, X)

	Executes the random subspace ensemble.

	# Arguments
	 * `ensemble::RandomSubspace.SubspaceEnsemble` the subspace ensemble
	 * `X::AbstractMatrix` is the data; columns are assumed to be samples, lines variables
	"""
	function randomsubspace_exec(ensemble::SubspaceEnsembleClassic, X::T where T<:AbstractMatrix)
		if ensemble.parallel_execution
			out = @sync @parallel (vcat) for i in eachindex(ensemble.members)
				ensemble.f_exec(ensemble.members[i], getobs(datasubset(X, ensemble.idx[i], ObsDim.First())))
			end
		else
			out = [ensemble.f_exec(ensemble.members[i], getobs(datasubset(X, ensemble.idx[i], ObsDim.First()))) for i in eachindex(ensemble.members)]

		end

		return combiner_exec(ensemble.combiner, out)
	end

	# j4pr case: no execution function, the ensemble members contain it already ('f_exec' is not present at all)
	function randomsubspace_exec(ensemble::SubspaceEnsembleCell, X::T where T<:AbstractMatrix )
		if ensemble.parallel_execution
			out = @sync @parallel (vcat) for i in eachindex(ensemble.members)
				ensemble.members[i]( getobs(datasubset(X, ensemble.idx[i], ObsDim.First())) )
			end
		else
			out = [ensemble.members[i](getobs(datasubset(X, ensemble.idx[i],ObsDim.First()))) for i in eachindex(ensemble.members)]
		end

		return combiner_exec(ensemble.combiner, out)
	end
end

##########################
# FunctionCell Interface #
##########################
"""
	randomsubspace(f, L, M, combiner, replace; parallel_train=true, parallel_execution=false)

Creates an untrained random subspace ensemble.

# Arguments
 * `f::CellFunU` is a untrained cell corresponding to the classifer/transform that the ensemble will contain
 * `L::Int` is the ensemble size
 * `M::Int` is the number of variables for each ensemble member
 * `f_train::Function` is the training function
 * `combiner::ClassifierCombiner.AbstractCombiner` is the combiner
 * `replace::Bool` specifies whether the variable subsets can contain (`true`) or not (`false`) the same variable. If `true`, any number of subsets can be generated; if
 `false` a finite maximum number can be generated but all subsets will be distinct (default: `true`)

# Keyword arguments
 * `parallel_train::Bool=true` specifies whether to paralellize or not the training.
 * `parallel_execution::Bool=false` specifies whether to paralellize or not the execution.

For more information:
	[1] T.K. Ho "The Random Subspace Method for Constructing Decision Forests" 1998 IEEE Transactions on Pattern Analysis and Machine Intelligence 20(8): 832-844.
 	[2] L. Kuncheva et al. "Random Subspace Ensembles for fMRI Classification" IEEE Transactions on Medical Imaging 29 (2): 531-542.
	[3] L. Kuncheva "Combining Pattern Classifiers 2'nd Ed." 2014, ISBN 978-1-118-31523-1

Try `?j4pr.RandomSubspace.randomsubspace_train` and `?j4pr.RandomSubspace.randomsubspace_exec` and for more details.

# Examples
```
julia> using j4pr

julia> D=DataGenerator.iris() # Get the iris dataset
Iris Dataset, 150 obs, 4 vars, 1 target(s)/obs, 3 distinct values: "virginica"(50),"setosa"(50),"versicolor"(50)

julia> w=j4pr.knn(5,smooth=:ml) # 5-NN classifier, maximum likelihood posterior smoothing
5-NN classifier: smooth=ml, no I/O size information, untrained

julia> we=randomsubspace(w, 5, 2,j4pr.ClassifierCombiner.GeneralizedMeanCombiner(5,3,1.0)) # ensemble of 5 members, 2 variables; the combiner uses 3 classes
Sub-space ensemble (5-NN classifier: smooth=ml, 5 members × 2 vars), no I/O size information, untrained

julia> we_trained = D |>we # train ensemble
Sub-space ensemble (5-NN classifier: smooth=ml, 5 members × 2 vars), 4->3, trained

julia> +D |> we_trained # run ensemble on training data
3×150 Array{Float64,2}:
 1.0          1.0          1.0          1.0          1.0          1.0          1.0          …  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16
 2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16     2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  0.12
 2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16  2.22045e-16     1.0          1.0          1.0          1.0          1.0          1.0          0.88
```
"""
randomsubspace(f::F where F<:CellFunU, L::Int, M::Int,
	       combiner::ClassifierCombiner.AbstractCombiner=ClassifierCombiner.NoCombiner(),
	       replace::Bool=true; parallel_train::Bool=true, parallel_execution::Bool=false) =
	FunctionCell(randomsubspace, (f, L, M, combiner, replace), ModelProperties(), "Sub-space ensemble ($(f.tinfo), $L members × $M vars)";
			parallel_train=parallel_train, parallel_execution=parallel_execution) # untrained function cell

############################
# DataCell/Array Interface #
############################
"""
	randomsubspace(x, f, L, M, combiner, replace; parallel_train=true, parallel_execution=false)

Trains a random subspace ensemble using the data `x`.
"""
# Training
randomsubspace(x::T where T<:CellDataL, f::F where F<:CellFunU, L::Int, M::Int, combiner::ClassifierCombiner.AbstractCombiner=ClassifierCombiner.NoCombiner(),
	       replace::Bool=true; parallel_train::Bool=true, parallel_execution::Bool=false) =
	randomsubspace((getx!(x), gety(x)), f, L, M, combiner, replace; parallel_train=parallel_train, parallel_execution=parallel_execution)

randomsubspace(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, f::F where F<:CellFunU, L::Int, M::Int,
	combiner::ClassifierCombiner.AbstractCombiner=ClassifierCombiner.NoCombiner(), replace::Bool=true; parallel_train::Bool=true, parallel_execution::Bool=false) =
	randomsubspace((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), f, L, M, combiner, replace; parallel_train=parallel_train, parallel_execution=parallel_execution)

randomsubspace(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, f::F where F<:CellFunU, L::Int, M::Int,
	       combiner::ClassifierCombiner.AbstractCombiner=ClassifierCombiner.NoCombiner(), replace::Bool=true;
	       parallel_train::Bool=true, parallel_execution::Bool=false) =
begin

	@assert nobs(x[1]) == nobs(x[2]) "[randomsubspace] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Transform labels first
	enc = labelencn(x[2])
	yenc = label2ind.(x[2],enc)

	# Train model
	ensembledata = RandomSubspace.randomsubspace_train(x[1], yenc, L, M, f.f, nothing, combiner, replace;
						    parallel_train=parallel_train, parallel_execution=parallel_execution)

	# Build model properties
	modelprops = ModelProperties(nvars(x[1]), length(enc.label), enc)

	FunctionCell(randomsubspace, Model(ensembledata, modelprops),
	      "Sub-space ensemble ($(f.tinfo), $L members × $M vars)" ) # trained function cell
end

# Execution
randomsubspace(x::T where T<:CellData, model::Model{<:RandomSubspace.SubspaceEnsemble}) =
	datacell(randomsubspace(getx!(x), model), gety(x))
randomsubspace(x::T where T<:AbstractVector, model::Model{<:RandomSubspace.SubspaceEnsemble}) =
	randomsubspace(mat(x, LearnBase.ObsDim.Constant{2}()), model)
randomsubspace(x::T where T<:AbstractMatrix, model::Model{<:RandomSubspace.SubspaceEnsemble}) =
	RandomSubspace.randomsubspace_exec(model.data,x)