text:
export structure

"""
	phase_structure(datatype::DataType, args...)
Takes a DataType (such as `Int8`) and a series of integers to return
a sorted Tuple of those integers converted to that DataType. i.e. takes
a series of alleles and returns a genotype. Returns `missing` if args are
`missing`. Used internally in PopGen.structure file reader.
**Example**
```
phase_structure(Int8, 1,2,3,4,3,4,6,1)
(1, 1, 2, 3, 3, 4, 4, 6)
phase_structure(Int16, missing, missing)
missing
```
"""
function phase_structure(datatype::DataType, args...)
    all(ismissing.(args)) && return missing
    return Tuple(datatype.(sort([args...])))
end

"""
    structure(infile::String; kwargs...)
Load a Structure format file into memory as a PopData object.
- `infile::String` : path to Structure file

### Keyword Arguments
- `extracols::Integer`: how many additional optional columns there are beyond Stucture's POPDATA the reader needs to ignore (default: `0`)
    - these include POPFLAG, LOCDATA, or anything else you might have added
- `extrarows::Integer` : how many additional optional rows there are beyond the first row of locus names (default: `0`)
- `missingval::String`  : the value used to identify missing values in the data (default: `"-9"`)
- `silent::Bool`   : whether to print file information during import (default: `false`)
- `allow_monomorphic::Bool` : whether to keep monomorphic loci in the dataset (default: `false`)
- `faststructure::Bool`: whether the file is fastStructure format (default: `false`)

### File must follow this Structure format:
- the file is `tab` or `space` delimited **but not both**
- first row is locus names separated by the delimiter
    - leading/trailing whitespaces are tolerated
    - optional rows allowed **after** the locus names
- number of rows per sample = ploidy
    - e.g. if diploid, that sample would have 2 rows
    - multi-column variant not supported
    - all samples must have the same ploidy
- first data column is sample name
- second data column is population ID
    - optional columns allowed **after** the population ID (2nd) column
- remaining columns are the genotype for that individual for that locus

### Structure file example:
```
locus_1	locus_2	locus_3	locus_4	locus_5
walnut_01	1	-9	145	66	0	92
walnut_01	1	-9	-9	64	0	94
walnut_02	1	106	142	68	1	92
walnut_02	1	106	148	64	0	94
walnut_03	2	110	145	-9	0	92
walnut_03	2	110	148	66	1	-9
```

### fastStructure file format:
- the file is `tab` or `space` delimited **but not both**
- no first row of loci names
- number of rows per sample = ploidy
    - e.g. if diploid, that sample would have 2 rows
- first data column is sample name
- second data column is population ID
- remaining columns are the genotype for that individual for that locus
- usually, first 6 colums are empty (but not necessary)
- **no** extra rows or columns.
### fastStructure file example:
```
chestnut_01	1	-9	145	66	0	92
chestnut_01	1	-9	-9	64	0	94
chestnut_02	1	106	142	68	1	92
chestnut_02	1	106	148	64	0	94
chestnut_03	2	110	145	-9	0	92
chestnut_03	2	110	148	66	1	-9
```

## Example
```
walnuts = structure("juglans_nigra.str", extracols = 0, extrarows = 0)
```
"""
function structure(infile::String; silent::Bool = false, extracols::Int = 0, extrarows::Int = 0, allow_monomorphic::Bool = false, missingval::String = "-9", faststructure::Bool = false)
    # find the delimiter
    isfile(infile) || throw(ArgumentError("$infile not found."))
    first_row = strip(open(readline, infile))
    if occursin("\t", first_row) & occursin(" ", first_row)
        error("$infile contains both tab and space delimiters. Please format the file so it uses either one or the other.")
    elseif occursin("\t", first_row)
        delim = "\t"
    elseif occursin(" ", first_row)
        delim = " "
    else
        error("Please format $infile to be either tab or space delimited")
    end

    if faststructure == true
        data_row = 1
        #first_row = split(strip(open(readline, infile)), delim)
        n_loci = length(split(first_row, delim))
        locinames = ["marker_$i" for i in 1:n_loci-2]
    else
        data_row = 2 + extrarows
        locinames = Symbol.(split(first_row, delim))
    end
    # read in the file as a table
    geno_parse = CSV.read(
        infile,
        DataFrame,
        delim = delim,
        header = false,
        skipto = data_row,
        ignorerepeated = true,
        missingstring = [missingval],
        types = Dict(1=>String, 2=>String)
    )
    # ignore any extra columns
    if !iszero(extracols)
        geno_parse = geno_parse[!, Not(collect(3:2+n))]
    end

    # fix names, just in case
    rename!(geno_parse, append!([:name, :population], locinames))
    geno_parse.name .= replace.(geno_parse.name, "-" => "_")
    geno_parse = stack(geno_parse, Not([1,2]))
    rename!(geno_parse, [:name, :population, :locus, :genotype])

    # convert columns to PooledArrays
    geno_parse.name = PooledArray(geno_parse.name, compress = true)
    geno_parse.population = PooledArray(geno_parse.population, compress = true)
    geno_parse.locus = PooledArray(geno_parse.locus, compress = true)
    if !silent
        @info "\n $(truncatepath(abspath(infile)))\n data: loci = $(length(geno_parse.locus.pool)), samples = $(length(geno_parse.name.pool)), populations = $(length(geno_parse.population.pool))"
        println()
    end
    grp = groupby(geno_parse, [:name, :population, :locus])

    geno_df =
        try
            DataFrames.combine(grp, 4 => _SNP => :genotype)
        catch
            DataFrames.combine(grp, 4 => _MSat => :genotype)
        end
    pd_out = PopData(geno_df)
    !allow_monomorphic && dropmonomorphic!(pd_out)
    return pd_out
end

"""
    structure(data::PopData; filename::String, faststructure::Bool, delim::String)
Write a `PopData` object to a Stucture format file
- `data`: the `PopData` object you wish to write to a Structure file
### keyword arguments
- `filename`: a `String` of the output filename
- `delim` : a `String` of either `"tab"` or `"space"` indicating the delimiter (default: `"tab"`)
- `faststructure`: true/false of whether the output should be formatted for fastStructure (default: `false`)

```
cats = @nancycats;
fewer_cats = omit(cats, name = samplenames(cats)[1:10]);
structure(fewer_cats, filename = "filtered_nancycats.str", faststructure = true)
```
"""
function structure(data::PopData; filename::String, faststructure::Bool = false, delim::String = "tab")
    # index both dataframes
    genos_gdf = groupby(data.genodata, :name)
    meta_gdf = groupby(data.metadata.sampleinfo, :name)
    # get the sample names to iterate keys over

    outfile = open(filename, "w")

    # check delimiter
    if delim == "tab"
        dlm = "\t"
    elseif delim == "space"
        dlm = " "
    else
        throw(ArgumentError("Please choose from either \"tab\" (default) or \"space\" delimiters."))
    end

    faststructure == false && println(outfile, join(loci(data), dlm))

    # remap populations as integers
    pops = unique(data.sampleinfo.population)
    pop_mappings = Dict{String,Integer}()
    [pop_mappings[j] = i for (i,j) in enumerate(pops)]

    for sampl in data.sampleinfo.name
        ploid = meta_gdf[(name = sampl,)].ploidy |> first
        pop_id = meta_gdf[(name = sampl,)].population |> first

        # copy so as to not overwrite
        genos = copy(genos_gdf[(name = sampl,)].genotype)

        # replace missing with -9's
        miss_idx = findall(ismissing, genos)
        for i in miss_idx
            genos[i] = phase_structure(Int8, fill(-9, ploid)...)
        end

        # write the alleles to the file
        for allele in 1:ploid
            allele_row = join(getindex.(genos, allele), dlm)
            println(outfile, sampl, dlm, pop_mappings[pop_id], dlm, allele_row)
        end
    end
    close(outfile)
end


text:
using Combinatorics

meandiff(a::Vector{T}, b::Vector{T}) where T <: Real = mean(a) - mean(b)

function bifurcate(a::AbstractVector, sel::Vector{T}) where T <: Integer
    x         = a[sel]
    asel      = trues(length(a))
    asel[sel] = false
    y         = a[asel]
    return x, y
end

function permutation_test(treated::Vector{T}, control::Vector{T}) where T <: Real
    effect0 = meandiff(treated, control)
    pool    = vcat(treated, control)
    tlen    = length(treated)
    plen    = length(pool)
    better = worse = 0
    for subset in combinations(1:plen, tlen)
        t, c = bifurcate(pool, subset)
        if effect0 < meandiff(t, c)
            better += 1
        else
            worse += 1
        end
    end
    return better, worse
end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ---
# title: 805. Split Array With Same Average
# id: problem805
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math
# link: <https://leetcode.com/problems/split-array-with-same-average/description/>
# hidden: true
# ---
#
# In a given integer array A, we must move every element of A to either list B
# or list C. (B and C initially start empty.)
#
# Return true if and only if after such a move, it is possible that the average
# value of B is equal to the average value of C, and B and C are both non-empty.
#
#
#
#     **Example :**
#     Input:
#     [1,2,3,4,5,6,7,8]
#     Output: true
#     Explanation: We can split the array into [1,4,5,8] and [2,3,6,7], and both of them have the average of 4.5.
#
#
# **Note:**
#
#   * The length of `A` will be in the range [1, 30].
#   * `A[i]` will be in the range of `[0, 10000]`.
#
#
#
#
## @lc code=start
using LeetCode

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


text:
export inv_rotation

################################################################################
#                              Inverse Rotations
################################################################################

"""
    inv_rotation(R)

Compute the inverse rotation of `R`, which can be:

* A direction cosina matrix (`DCM`);
* An Euler angle and axis (`EulerAngleAxis`);
* A set of Euler anlges (`EulerAngles`); or
* A quaternion (`Quaternion`).

The output will have the same type as `R` (DCM or quaternion).

# Remarks

If `R` is a DCM, than its transpose is computed instead of its inverse to reduce
the computational burden. The both are equal if the DCM has unit norm. This must
be verified by the user.

If `R` is a quaternion, than its conjugate is computed instead of its inverse to
reduce the computational burden. The both are equal if the quaternion has unit
norm. This must be verified by the used.

# Example

```julia-repl
julia> D = angle_to_dcm(+pi/3,+pi/4,+pi/5,:ZYX);

julia> inv_rotation(D)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  0.353553  -0.492816  0.795068
  0.612372   0.764452  0.201527
 -0.707107   0.415627  0.572061

julia> ea = EulerAngleAxis(30*pi/180, [1;0;0]);

julia> inv_rotation(ea)
EulerAngleAxis{Float64}:
  Euler angle:   0.5236 rad ( 30.0000 deg)
   Euler axis: [ -1.0000,  -0.0000,  -0.0000]

julia> Θ = EulerAngles(-pi/3, -pi/2, -pi, :YXZ);

julia> inv_rotation(Θ)
EulerAngles{Float64}:
  R(Z):   3.1416 rad ( 180.0000 deg)
  R(X):   1.5708 rad (  90.0000 deg)
  R(Y):   1.0472 rad (  60.0000 deg)

julia> q = angle_to_quat(+pi/3,+pi/4,+pi/5,:ZYX);

julia> inv_rotation(q)
Quaternion{Float64}:
  + 0.8200711519756747 - 0.06526868310243991.i - 0.45794027732580056.j - 0.336918398289752.k
```

"""
@inline inv_rotation(D::DCM) = D'
@inline inv_rotation(ea::EulerAngleAxis) = inv(ea)
@inline inv_rotation(Θ::EulerAngles) = inv(Θ)
@inline inv_rotation(q::Quaternion) = conj(q)


text:
export raytrace

# ----------------- #
# Blinn Phong Model #
# ----------------- #

"""
    light(s::Object, origin, direction, dist, lgt::Light, eye_pos, scene, obj_num, bounce)

Implements the Blinn Phong rendering algorithm. This function is merely for
internal usage and should in no case be called by the user. This function is
quite general and supports user defined **Objects**. For support of custom
Objects have a look at the examples.

!!! warning
    Don't try to use this function by itself. But if you are a person who likes
    to ignore warnings look into the way [`raytrace`](@ref) calls this.
"""
function light(s::Object, origin, direction, dist, lgt::Light, eye_pos,
               scene, obj_num, bounce)
    pt = origin + direction * dist
    normal = get_normal(s, pt, direction)
    dir_light, intensity = get_shading_info(lgt, pt)
    dir_origin = normalize(eye_pos - pt)
    nudged = pt + normal * 0.0001f0 # Nudged to miss itself

    # Shadow
    light_distances = map(x -> intersect(x, nudged, dir_light), scene)
    seelight = fseelight(obj_num, light_distances)

    # Ambient
    color = get_color(s, pt, Val(:ambient))

    # Lambert Shading (diffuse)
    visibility = max.(dot(normal, dir_light), 0.0f0)
    color += ((get_color(s, pt, Val(:diffuse)) * intensity) * visibility) * seelight

    # Reflection
    if bounce < 2
        rayD = normalize(direction - normal * 2.0f0 * dot(direction, normal))
        color += raytrace(nudged, rayD, scene, lgt, eye_pos, bounce + 1) * reflection(s)
    end

    # Blinn-Phong shading (specular)
    phong = dot(normal, normalize(dir_light + dir_origin))
    color += (get_color(s, pt, Val(:specular)) *
              (clamp.(phong, 0.0f0, 1.0f0) .^ specular_exponent(s))) * seelight

    return color
end

"""
    raytrace(origin::Vec3, direction::Vec3, scene::Vector, lgt::Light,
             eye_pos::Vec3, bounce::Int)
    raytrace(origin::Vec3, direction::Vec3, scene, lgt::Vector{Light},
             eye_pos::Vec3, bounce::Int)
    raytrace(origin::Vec3, direction::Vec3, scene::BoundingVolumeHierarchy,
             lgt::Light, eye_pos::Vec3, bounce::Int)
    raytrace(;origin = nothing, direction = nothing, scene = nothing,
              light_source = nothing, global_illumination = false)

Computes the color contribution to every pixel by tracing every single ray.
Internally it calls the `light` function which implements Blinn Phong Rendering
and adds up the color contribution for each object.

The `eye_pos` is simply the `origin` when called by the user. However, the origin
keeps changing across the recursive calls to this function and hence it is
necessary to keep track of the `eye_pos` separately.

The `bounce` parameter allows the configuration of global illumination. To turn off
global illumination set the `bounce` parameter to `>= 2`. As expected rendering is
much faster if global illumination is off but at the same time is much less photorealistic.

!!! note
    The support for multiple lights is primitive as we loop over the lights.
    Even though it is done in a parallel fashion, it is not the best way to
    do so. Nevertheless it exists just for the sake of experimentation.

!!! note
    None of the parameters (except global_illumination) in the keyword argument version of
    [`raytrace`](@ref) is optional. It is just present for convenience.
"""
raytrace(;origin = nothing, direction = nothing, scene = nothing,
          light_source = nothing, global_illumination = false) =
    raytrace(origin, direction, scene, light_source, origin, global_illumination ? 0 : 2)

function raytrace(origin::Vec3, direction::Vec3, scene::Vector,
                  lgt::Light, eye_pos::Vec3, bounce::Int = 0)
    distances = map(x -> intersect(x, origin, direction), scene)

    dist_reshaped = reducehcat(distances)
    nearest = map(idx -> minimum(dist_reshaped[idx, :]), 1:size(dist_reshaped, 1))

    h = .!isinf.(nearest)

    color = Vec3(0.0f0)

    for (c, (s, d)) in enumerate(zip(scene, distances))
        hit = map((x, y, z) -> ifelse(y == z, x, zero(x)), h, d, nearest)
        if any(hit)
            dc = extract(hit, d)
            originc = extract(hit, origin)
            dirc = extract(hit, direction)
            cc = light(s, originc, dirc, dc, lgt, eye_pos, scene, c, bounce)
            color += place(cc, hit)
        end
    end

    return color
end

# FIXME: This is a temporary solution. Long term solution is to support
#        a light vector in the `light` function itself.
function raytrace(origin::Vec3, direction::Vec3, scene,
                  lgt::Vector{L}, eye_pos::Vec3, bounce::Int = 0) where {L<:Light}
    colors = pmap(x -> raytrace(origin, direction, scene, x, eye_pos, bounce), lgt)
    return sum(colors)
end

# ------------------------ #
# General Helper Functions #
# ------------------------ #

"""
    fseelight(n::Int, light_distances)

Checks if the ``n^{th}`` object in the scene list can see the light source.
"""
function fseelight(n, light_distances)
    ldist = reducehcat(light_distances)
    seelight = map(idx -> minimum(ldist[idx, :]) == ldist[idx,n], 1:size(ldist, 1))
    return seelight
end

# fseelight(n, light_distances) = map((x...) -> min(x...) == x[n], light_distances...)

"""
    reducehcat(x)

This is simply `reduce(hcat(x))`. The version of Zygote we are currently using
can't differentiate through this function. So we define a custom adjoint for this.
"""
reducehcat(x) = reduce(hcat, x)


text:
using Documenter
using POMDPSimulators

makedocs(
    modules = [POMDPSimulators],
    format = Documenter.HTML(),
    sitename = "POMDPSimulators.jl",
    pages = ["index.md",
             "which.md",
             "rollout.md",
             "parallel.md",
             "history_recorder.md",
             "histories.md",
             "stepthrough.md",
             "display.md",
             "sim.md"]
)

deploydocs(
    repo = "github.com/JuliaPOMDP/POMDPSimulators.jl.git",
)


text:
export term, @sym

const DOMAIN_TYPES = Dict()

include("integer.jl")
include("real.jl")
include("complex.jl")
include("number.jl")

const SymbolicType = Union{SymInteger, SymReal, SymComplex, SymNumber}

term(x) = x
track(f, x::T, xs::T...) where {T <: SymbolicType} = T(Term(f, collect(Any, map(term, (x, xs...)))))
track(f, xs...) = track(f, promote(xs...)...)

for T in [SymInteger, SymReal, SymComplex, SymNumber]
    @eval Base.:(^)(x::$T, n::Number) = track(^, x, n)
    @eval Base.:(^)(x::$T, n::Integer) = track(^, x, n)
    @eval Base.:(^)(x::$T, n::Real) = track(^, x, n)
    @eval Base.:(^)(x::$T, n::AbstractFloat) = track(^, x, n)
end

macro sym(ex::Expr)
    symm(ex)
end

function symm(ex::Expr)
    if ex.head === :call
        ts, xs = parse_sym(ex)
    else
        ex.head === :tuple || throw(Meta.ParseError("expect a tuple, e.g @sym x, y, z"))

        xs = Expr(:tuple)
        ts = Expr(:tuple)
        for each in ex.args
            name, e = parse_sym(each)
            push!(xs.args, e)
            push!(ts.args, name)
        end
    end
    return :($(esc(ts)) = $xs)
end

function parse_sym(name::Symbol)
    name, :(SymNumber($(Meta.quot(name))))
end

function parse_sym(ex::Expr)
    ex.head === :call || throw(Meta.ParseError("expect sym in domain, got $ex"))
    if ex.args[1] === :in # domain
        name, domain = ex.args[2], ex.args[3]
        if haskey(DOMAIN_TYPES, domain)
            T = DOMAIN_TYPES[domain]
            return name, :($T($(Meta.quot(name))))
        else
            throw(Meta.ParseError("unsupported domain $(domain)"))
        end
    else
        throw(Meta.ParseError("invalid syntax $ex"))
    end
end


text:
abstract type FunctionalImage <: Imaging end

const FMRI{P,S,T,A,Rl,Ri,Ei,E} = FunctionalMagneticResonanceImaging{P,S,T,A,Rl,Ri,Ei,E}
datatypefolder(::FMRI) = "func"

struct TaskFMRI{P,S,T,A,Rl,Ri,Ei,E} <: FMRI{P,S,T,A,Rl,Ri,Ei,E}
    participant::P
    seslabel::S
    tasklabel::T
    acqlabel::A
    reclabel::Rl
    runindex::Ri
    echoindex::Ei
    ext::E
end

struct RestingStateFMRI{P,S,T,A,Rl,Ri,Ei,E} <: FMRI{P,S,T,A,Rl,Ri,Ei,E}
    participant::P
    seslabel::S
    tasklabel::T
    acqlabel::A
    reclabel::Rl
    runindex::Ri
    echoindex::Ei
    ext::E
end

"""
    tasklabel(bidsp::FMRI) -> "_task-<label>"

Required
"""
tasklabel(bidsp::FMRI) = "_task-$(bidsp.tasklabel)"

acqlabel(bidsp::FMRI{P,S,A,L,Rl,Ri}) where {P,S,A,L,Rl,Ri} = "_acq-$(bidsp.seslabel)"
acqlabel(bidsp::FMRI{P,S,Nothing,L,RL,Ri}) where {P,S,L,Ri} = ""

seslabel(bidsp::FMRI{P,S,A,L,Rl,Ri}) where {P,S,A,L,Rl,Ri} = "_ses-$(bidsp.seslabel)"
seslabel(bidsp::FMRI{P,Nothing,A,L,RL,Ri}) where {P,A,L,Rl,Ri} = ""

recindex(bidsp::FMRI{P,S,A,L,Rl,Ri}) where {P,S,A,L,Rl,Ri} = "_rec-$(bidsp.reclabel)"
recindex(bidsp::FMRI{P,S,A,L,Nothing,Ri}) where {P,S,A,L,Ri} = ""

runindex(bidsp::FMRI{P,S,A,L,Rl,Ri}) where {P,S,A,L,Rl,Ri} = "_run-$(bidsp.runindex)"
runindex(bidsp::FMRI{P,S,A,L,Rl,Nothing}) where {P,S,A,L,Rl} = ""

"""
    echoindex(::FMRI) -> "_echo-<index>"
"""
echoindex(bidsp::FMRI{P,S,A,L,Rl,Ri,Ei}) where {P,S,A,L,Rl,Ri,Ei} = "_echo-$(bidsp.echoindex)"
echoindex(bidsp::FMRI{P,S,A,L,Rl,Ri,Nothing}) where {P,S,A,L,Rl,Ri} = ""

function filename(bidsp::FMRI)
    string(participant(bidsp),
           seslabel(bidsp),
           tasklabel(bidsp),
           acqlabel(bidsp),
           reclabel(bidsP),
           runindex(bidsp),
           echoindex(bidsp),
           "_bold",
           join(string.(".", extensions(bidsp)))
          )
end

struct MultiBandFMRI{P,S,T,A,Rl,Ri,Ei,E} <: FMRI{P,S,T,A,Rl,Ri,Ei,E}
    participant::P
    seslabel::S
    tasklabel::T
    acqlabel::A
    reclabel::Rl
    runindex::Ri
    echoindex::Ei
    ext::E
end

function filename(bidsp::MultiBandFMRI)
    string(participant(bidsp),
           seslabel(bidsp),
           tasklabel(bidsp),
           acqlabel(bidsp),
           reclabel(bidsP),
           runindex(bidsp),
           echoindex(bidsp),
           "_sbref",
           join(string.(".", extensions(bidsp)))
          )
end


text:
using MPI, MPIArrays
MPI.Init()
rank = MPI.Comm_rank(MPI.COMM_WORLD)
x = MPIArray{Int64}(5, 5)
y = MPIArray{Int64}(5, 5)
z = CyclicMPIArray(Int64, 9, 9,  proc_grids=(2, 2), blocksizes=(2, 2))

forlocalpart!(x->fill!(x, rank), z)

# this is error, because z.win is nothing
# show(stdout, "text/plain", z)

rma!(z) do
    if rank == 0
        show(stdout, "text/plain", z)
        println()
    end
end

MPI.Finalize()


text:
function ner_datadep_register()
    register(DataDep("NER Model Weights",
        """
        The weights for NER Sequence Labelling Model.
        """,
        "https://github.com/JuliaText/TextAnalysis.jl/releases/download/v0.6.0/ner_weights.tar.xz",
        "6290353b66c9bdbb794ddcb6063ab52c30145d3918f2f115f19e21fa994282e6",
        post_fetch_method = function(fn)
            unpack(fn)
            dir = "weights"
            innerfiles = readdir(dir)
            mv.(joinpath.(dir, innerfiles), innerfiles)
            rm(dir)
        end
    ))

    register(DataDep("NER Model Dicts",
        """
        The character and words dict for NER Sequence Labelling Model.
        """,
        "https://github.com/JuliaText/TextAnalysis.jl/releases/download/v0.6.0/ner_dicts.tar.xz",
        "40cfa37da216b990eb9c257aa7994e34d7a7a59d69b2506c6f39120f2688dc11",
        post_fetch_method = function(fn)
            unpack(fn)
            dir = "model_dicts"
            innerfiles = readdir(dir)
            mv.(joinpath.(dir, innerfiles), innerfiles)
            rm(dir)
        end
    ))
end


text:
__precompile__()

module SimpleGradientDescent

import MathProgBase
import MathProgBase: setwarmstart!, optimize!, status, getsolution
import ForwardDiff
import StaticArrays: SVector, MVector

export DifferentiableFunction,
       NaiveSolver,
       getmodel,
       setwarmstart!,
       optimize!,
       getsolution,
       status

include("differentiablefunction.jl")
include("solver.jl")

end # module


text:
module CDCS

using SparseArrays
using MATLAB

export cdcs

# The fields should be integer values but the type should be `Float64` to work
# with CDCS
mutable struct Cone
    f::Float64 # number of free primal variables / linear dual equality constraints
    l::Float64 # length of LP cone
    q::Vector{Float64} # list of second-order cone dimension
    s::Vector{Float64} # list of semidefinite constraints
end
Cone(f::Real, l::Real) = Cone(f, l, Float64[], Float64[])

to_vec(x::Vector) = x
to_vec(x::Float64) = [x]
function to_vec(x::AbstractMatrix)
    @assert isempty(x) || isone(size(x, 2))
    return vec(x)
end

# Solve the primal/dual pair
# min c'x,      max b'y
# s.t. Ax = b,   c - A'x ∈ K
#       x ∈ K
function cdcs(A::Union{Matrix{Float64}, SparseMatrixCSC{Float64}},
              b::Vector{Float64}, c::Vector{Float64},
              K::Cone = Cone(0, size(A, 2));
              kws...)
    @assert size(A, 1) == length(c)
    @assert size(A, 2) == length(b)
    @assert size(A, 1) == K.f + K.l + sum(K.q) + sum(K.s.^2)
    pars = Dict{String, Any}(string(key) => value for (key, value) in kws)
    # There are 4 output arguments: x, y, z and info so we use `4` below
    x, y, z, info = mxcall(:cdcs, 4, A, b, c, K, pars)
    return to_vec(x), to_vec(y), to_vec(z), info
end

include("MOI_wrapper.jl")

end # module


text:
using ClassicalOrthogonalPolynomials, ContinuumArrays, BandedMatrices, LazyArrays, ForwardDiff, Test
import LazyArrays: rowsupport, colsupport

@testset "Ultraspherical" begin
    @testset "Transforms" begin
        U = Ultraspherical(1)
        x = axes(U,1)
        Un = U[:,Base.OneTo(5)]
        @test factorize(Un) isa ContinuumArrays.TransformFactorization
        @test (Un \ x) ≈ [0,0.5,0,0,0]
        @test (U * (U \ exp.(x)))[0.1] ≈ exp(0.1)
    end

    @testset "Operators" begin
        @testset "Lowering" begin
            λ = 1
            wC1 = Weighted(Ultraspherical(λ))
            wC2 = Weighted(Ultraspherical(λ+1))
            L = wC1 \ wC2
            @test L isa BandedMatrix
            @test bandwidths(L) == (2,0)
            @test L[1,1] ≈ 2λ*(2λ+1)/(4λ*(λ+1))
            @test L[3,1] ≈ -2/(4λ*(λ+1))
            u = [randn(10); zeros(∞)]
            @test wC1[0.1,:]'*(L*u) ≈ wC2[0.1,:]'*u
        end

        @testset "Weighted Derivative" begin
            T = Chebyshev()
            wC1 = Weighted(Ultraspherical(1))
            wC2 = Weighted(Ultraspherical(2))
            x = axes(wC2,1)
            D = Derivative(x)

            @test wC1 \ (D*wC2) isa BandedMatrix
            @test (wC1 \ (D*wC2))[1:5,1:5] == BandedMatrix(-1 => [-1.5,-4,-7.5,-12])

            u = wC2 * [randn(5); zeros(∞)]
            @test (D*u)[0.1] ≈ ((D*T) * (T\u))[0.1]
        end

        @testset "Interrelationships" begin
            @testset "Chebyshev–Ultrashperical" begin
                T = ChebyshevT()
                U = ChebyshevU()
                C = Ultraspherical(2)
                D = Derivative(axes(T,1))

                @test C\C === pinv(C)*C === Eye(∞)
                D₀ = U\(D*T)
                D₁ = C\(D*U)
                @test D₁ isa BandedMatrix
                @test (D₁*D₀)[1:10,1:10] == diagm(2 => 4:2:18)
                @test D₁*D₀ isa MulMatrix
                @test bandwidths(D₁*D₀) == (-2,2)

                S₁ = (C\U)[1:10,1:10]
                @test S₁ isa BandedMatrix{Float64}
                @test S₁ == diagm(0 => 1 ./ (1:10), 2=> -(1 ./ (3:10)))

                @test (U\C)[1:10,1:10] ≈ inv((C\U)[1:10,1:10])
                @test (T\C)[1:10,1:10] ≈ inv((C\T)[1:10,1:10])
                @test bandwidths(U\C) == bandwidths(T\C) == (0,∞)
                @test colsupport(U\C,5) == colsupport(T\C,5) == 1:5
                @test rowsupport(U\C,5) == rowsupport(T\C,5) == 5:∞
            end
            @testset "Legendre" begin
                @test Ultraspherical(0.5) \ (UltrasphericalWeight(0.0) .* Ultraspherical(0.5)) == Eye(∞)
                @test Legendre() \ (UltrasphericalWeight(0.0) .* Ultraspherical(0.5)) == Eye(∞)
                @test (Legendre() \ Ultraspherical(1.5))[1:10,1:10] ≈ inv(Ultraspherical(1.5) \ Legendre())[1:10,1:10]
            end
        end

        @testset "Conversion" begin
            R = Ultraspherical(3.5) \ Ultraspherical(0.5)
            c = [[1,2,3,4,5]; zeros(∞)]
            @test Ultraspherical(3.5)[0.1,:]' * (R * c) ≈ Ultraspherical(0.5)[0.1,:]' * c
            Ri = Ultraspherical(0.5) \ Ultraspherical(3.5)
            @test Ri[1:10,1:10] ≈ inv(R[1:10,1:10])
            @test Ultraspherical(0.5)[0.1,:]' * (Ri * c) ≈ Ultraspherical(3.5)[0.1,:]' * c
        end
    end

    @testset "test on functions" begin
        T = Chebyshev()
        U = Ultraspherical(1)
        D = Derivative(axes(T,1))
        f = T*Vcat(randn(10), Zeros(∞))
        @test (U*(U\f)).args[1] isa Ultraspherical
        @test (U*(U\f))[0.1] ≈ f[0.1]
        @test (D*f)[0.1] ≈ ForwardDiff.derivative(x -> (ChebyshevT{eltype(x)}()*f.args[2])[x],0.1)
    end

    @testset "Evaluation" begin
        C = Ultraspherical(2)
        @test @inferred(C[0.1,Base.OneTo(0)]) == Float64[]
        @test @inferred(C[0.1,Base.OneTo(1)]) == [1.0]
        @test @inferred(C[0.1,Base.OneTo(2)]) == [1.0,0.4]
        @test @inferred(C[0.1,Base.OneTo(3)]) == [1.0,0.4,-1.88]
    end

    @testset "special syntax" begin
        @test ultrasphericalc.(0:5, 2, 0.3) == Ultraspherical(2)[0.3, 1:6]
    end

    @testset "Ultraspherical vs Chebyshev and Jacobi" begin
        @test Ultraspherical(1) == ChebyshevU()
        @test ChebyshevU() == Ultraspherical(1)
        @test Ultraspherical(0) ≠ ChebyshevT()
        @test ChebyshevT() ≠ Ultraspherical(0)
        @test Ultraspherical(1) ≠ Jacobi(1/2,1/2)
        @test Jacobi(1/2,1/2) ≠ Ultraspherical(1)
        @test Ultraspherical(1/2) == Jacobi(0,0)
        @test Ultraspherical(1/2) == Legendre()
        @test Jacobi(0,0) == Ultraspherical(1/2)
        @test Legendre() == Ultraspherical(1/2)

        @test Ultraspherical(1/2) \ (JacobiWeight(0,0) .* Jacobi(0,0)) isa Diagonal
        @test (JacobiWeight(0,0) .* Jacobi(0,0)) \ Ultraspherical(1/2) isa Diagonal
    end

    @testset "D^2 * mapped" begin
        T = chebyshevt(0..1)
        C = ultraspherical(2,0..1)
        r = axes(T,1)
        D = Derivative(r)
        D₂ = C \ (D^2 * T)
        # r²D₂ = C \ (r.^2 .* (D^2 * T))

        c = [randn(100); zeros(∞)]
        @test C[0.1,:]'*(D₂ * c) ≈ 4*(Derivative(axes(ChebyshevT(),1))^2 * (ChebyshevT() * c))[2*0.1-1]
        @test_broken C[0.1,:]'*(r²D₂ * c) ≈ 0.1^2 * C[0.1,:]'*(D₂ * c)
    end

    @testset "BigFloat" begin
        U = Ultraspherical{BigFloat}(1)
        T = ChebyshevT{BigFloat}()
        x = axes(U,1)
        D = Derivative(x)

        @test Weighted(T) \ (D * Weighted(U)) isa BandedMatrix{BigFloat}

        C³ = Ultraspherical{BigFloat}(3)
        c = [1; 2; 3; zeros(BigFloat,∞)]
        @test C³[big(1)/10,:]'*(C³ \ U) * c ≈ U[big(1)/10,:]'c
    end
end


text:
"""
The base type of a fractional integral sense, in FractionalCalculus.jl, all of the fractional integral senses belong to ```FracIntAlg```
"""
abstract type FracIntAlg end

"""
Riemann-Liouville sense fractional integral algorithms

Note this two algorithms belong to direct compute, precise are ensured, but maybe cause more memory allocation and take more compilation time.
"""
abstract type RLInt <: FracIntAlg end

"""
Using the direct mathematic expression:

```math
(J^αf)(t)=\\frac{1}{\\Gamma(α)} \\int_0^t(t-τ)^{α-1}f(τ)dτ
```

By using [QuadGK](https://github.com/JuliaMath/QuadGK.jl) calculate the integration and obtain the value.
"""
struct RL_Direct <: RLInt end

"""
With first derivative known, we can use the Riemann Liouville sense to obtain the fractional integral more effcient.
"""
struct RL_Direct_First_Diff_Known <: RLInt end

"""
@article{LI20113352,
title = {Numerical approaches to fractional calculus and fractional ordinary differential equation},
author = {Changpin Li and An Chen and Junjie Ye},
}
"""

"""
Using piecewise linear interpolation:

```math
    y_n(t)=\\frac{t-t_{i+1}}{t_i-t_{i+1}}y(t_i)+\\frac{t-t_i}{t_{i+1}-t_i}y(t_{i+1})
```

Constitute the original function with the interpolation and implement the summation.
"""
struct RL_Piecewise <: RLInt end

"""
@book{oldham_spanier_1984,
title={The fractional calculus: Theory and applications of differentiation and integration to arbitrary order},
author={Oldham, Keith B. and Spanier, Jerome},
year={1984}}
"""

"""
Using the **Staircase approximation** to approximate the original function and implement the summation.
"""
struct RLInt_Approx <: RLInt end

"""
Deploying the [Linear Interpolation](https://en.wikipedia.org/wiki/Linear_interpolation) between ``f_{j+1}`` and ``f_j``, **RL_LinearInterp** method is more precise than **RLInt_Approx**.
"""
struct RL_LinearInterp <: RLInt end

"""
@article{2009,
title={Matrix approach to discrete fractional calculus II: Partial fractional differential equations},
DOI={10.1016/j.jcp.2009.01.014},
author={Podlubny, Igor and Chechkin, Aleksei and Skovranek, Tomas and Chen, YangQuan and Vinagre Jara, Blas M.},
}

"""

"""
Using [Triangular Strip Matrix](https://en.wikipedia.org/wiki/Triangle_strip) to discrete the integral.
"""
struct RLInt_Matrix <: RLInt end

####################################
###     Type defination done     ###
####################################

"""
Check if the format of nargins is correct
"""
function checks(f, α, start_point, end_point)
    α % 1 != 0 ? nothing : error("Decimal number only!")
    if isa(end_point, Number)
        start_point < end_point ? nothing : DomainError("Domain error! Start point must smaller than end point")
    elseif isa(end_point, AbstractArray)
        start_point < minimum(end_point) ? nothing : DomainError("Vector domain error! Start point must smaller than end point")
    else
        ErrorException("Please input correct point you wan tto compute")
    end

    #The fractional integral of number is relating with the end_point.
    if isa(f, Number)
        f == 0 ? 0 : 2*f*sqrt(end_point/pi)
    end
end

"""
# Riemann Liouville sense fractional integral

    fracint(f::Function, α, start_point, end_point, h, RL_Direct())

### Example:

```julia-repl
julia> fracint(x->x^5, 0.5, 0, 2.5, 1e-8, RL_Direct())
```

Riemann_Liouville fractional integral using complex step differentiation.
Returns a tuple (1.1639316474512205, 1.0183453796725215e-8), which contains the value of this derivative is 1.1639316474512205, and the error estimate is 1.0183453796725215e-8
"""
function fracint(f::Union{Function, Number}, α, start_point, end_point::Real, h, ::RL_Direct)
    checks(f, α, start_point, end_point)

    temp1 = f(start_point) .* (end_point-start_point) .^α
    #Use Complex differentiation to obtain the differentiation
    g(τ) = imag(f(τ .+ 1*im*h) ./ h) .* (end_point-τ) .^α
    temp2 = quadgk(g, start_point, end_point)
    result = (temp1 .+ temp2) ./gamma(α+1)
    return result
end

function fracint(f::Union{Function, Number}, α::Float64, start_point, end_point::AbstractArray, h, ::RL_Direct)::Vector
    ResultArray = Float64[]

    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracint(f, α, start_point, value, h, RL_Direct())[1])
    end
    return ResultArray
end

"""
# Riemann Liouville sense fractional integral with first diff known.

    fracint(f, fd, α, start_point, end_point, RL_Direct_First_Diff_Known())

### Example

```julia-repl
julia> fracint(x->x^5, x->5*x^4, 0.5, 0, 2.5, RL_Direct_First_Diff_Known())
```

With first order derivative known, we can directly use it in the computation of α order fraction integral
"""
function fracint(f::Function, fd::Function, α, start_point, end_point, ::RL_Direct_First_Diff_Known)
    checks(f, α, start_point, end_point)

    #The fractional integral of number is relating with the end_point.
    if typeof(f) <: Number
        if f == 0
            return 0
        else
            return 2*f*sqrt(end_point/pi)
        end
    end

    temp1 = f(start_point) .* (end_point-start_point) .^α
    g(τ) = fd(τ) .* (end_point-τ) .^α
    temp2 = quadgk(g, start_point, end_point)
    result = (temp1 .+ temp2) ./ gamma(α+1)
    return result
end

function fracint(f::Union{Function, Number}, α::Float64, start_point, end_point::AbstractArray, ::RL_Direct_First_Diff_Known)::Vector
    ResultArray = Float64[]

    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracint(f, α, start_point, value, RL_Direct_First_Diff_Known())[1])
    end
    return ResultArray
end

"""
# Riemann Liouville sense fractional integral using piecewise interpolation.

    fracint(f, α, end_point, h, RL_Piecewise())

### Example

```julia-repl
julia> fracint(x->x^5, 0.5, 2.5, 0.0001, RL_Piecewise())
```

By deploying Piecewise interpolation to approximate the original function, with small step size, this method is fast and take little memory allocation.
"""
function fracint(f::Union{Function, Number}, α::Float64, end_point, h, ::RL_Piecewise)::Float64
    checks(f, α, 0, end_point)

    #Init
    n=end_point/h
    result=0

    for i in range(0, n, step=1)
        result = result + W(i, n, α)*f(i*h)
    end

    result1 = result*h^α/gamma(α+2)
    return result1
end
function W(i, n, α)
    if i==0
        return n^α*(α+1-n)+(n-1)^(α+1)
    elseif i==n
        return 1
    else
        return (n-i-1)^(α+1)+(n-i+1)^(α+1)-2*(n-i)^(α+1)
    end
end

function fracint(f::Union{Function, Number}, α::Float64, end_point::AbstractArray, h, ::RL_Piecewise)::Vector
    ResultArray = Float64[]

    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracint(f, α, value, h, RL_Piecewise())[1])
    end
    return ResultArray
end

"""
# Riemann Liouville sense fractional integral approximation.

    fracint(f, α, end_point, h, RLInt_Approx())

### Example

```julia-repl
julia> fracint(x->x^5, 0.5, 2.5, 0.0001, RLInt_Approx())
```
"""
function fracint(f::Union{Function, Number}, α::Float64, end_point, h, ::RLInt_Approx)::Float64
    checks(f, α, 0, end_point)

    α = -α
    n = end_point/h
    result=0

    for i in range(0, n-1, step=1)
        result += (f(end_point-i*end_point/n)+f(end_point-(i+1)*end_point/n))/2*((i+1)^(-α)-i^(-α))
    end

    result1=result*end_point^(-α)*n^α/gamma(1-α)
    return result1
end

function fracint(f::Union{Function, Number}, α::Float64, end_point::AbstractArray, h, ::RLInt_Approx)::Vector
    ResultArray = Float64[]

    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracint(f, α, value, h, RLInt_Approx())[1])
    end
    return ResultArray
end

"""
# Riemann Liouville sense fractional integral using Linear interpolation.

    fracint(f, α, end_point, h, RL_LinearInterp())

### Example

```julia-repl
julia> fracint(x->x^5, 0.5, 2.5, 0.0001, RL_LinearInterp())
```

**RL_LinearInterp** is more complex but more precise.
"""
function fracint(f::Union{Function, Number}, α::Float64, end_point::Number, h, ::RL_LinearInterp)::Float64

    checks(f, α, 0, end_point)

    α = -α
    n = end_point/h
    result = 0

    for i in range(0, n-1, step=1)
        result+=((i+1)*f(end_point-i*end_point/n)-i*f(end_point-(i+1)*end_point/n))/(-α)*((i+1)^(-α)-i^(-α))+(f(end_point-(i+1)*end_point/n)-f(end_point-i*end_point/n))/(1-α)*((i+1)^(1-α)-i^(1-α))
    end

    result1=result*end_point^(-α)*n^α/gamma(-α)
    return result1
end

function fracint(f::Union{Function, Number}, α::Float64, end_point::AbstractArray, h, ::RL_LinearInterp)::Vector
    ResultArray = Float64[]

    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracint(f, α, value, h, RL_LinearInterp())[1])
    end
    return ResultArray
end

"""
# Riemann Liouville sense integral using Triangular Strip Matrix to discrete and compute.

    fracint(f, α, end_point, h, RLInt_Matrix())

Using Triangular Strip Matrix to approximate fractional integral.

### Example

```julia-repl
julia> fracint(x->x^5, 0.5, 2.5, 0.0001, RLInt_Matrix())
```

!!! info
    Triangular Strip Matrix method returns the derivative in the interval ``[0, T]`` in ```Vector```

!!! tip
    With the advancing Triangular Strip Matrix method, you can not only compute fractional integrals, integer order, higher order integral is also supported!!
Try to set α as an integer, arbitrary integer of course! I promise you would enjoy it😏
"""
function fracint(f, α::Number, end_point, h::Float64, ::RLInt_Matrix)
    N=Int64(end_point/h+1)
    tspan=collect(0:h:end_point)
    return J(N, α, h)*f.(tspan)
end

function omega(n, p)
    omega = zeros(n+1)

    omega[1]=1
    for i in range(1, n, step=1)
        omega[i+1]=(1-(p+1)/i)*omega[i]
    end

    return omega

end
function J(N, p, h::Float64)
    result=zeros(N, N)
    temp=omega(N, -p)

    for i in range(1, N, step=1)
        result[i, 1:i]=reverse(temp[1:i])
    end

    return h^p*result
end

## Macros for convenient computing.
"""
    @fracint(f, α, point)

Return the α-order integral of **f** at specific point.

```julia-repl
julia> @fracint(x->x, 0.5, 1)
0.7522525439593486
```
"""
macro fracint(f::Union{Function, Number}, α::Float64, point)
    return :(fracint($f, $α, $point, 0.0001, RLInt_Approx()))
end

"""
    @semifracint(f, point)

Return the semi-integral of **f** at specific point.

```julia-repl
julia> @semifracint(x->x, 1)
0.7522525439593486
```
"""
macro semifracint(f::Union{Function, Number}, point)
    return :(fracint($f, 0.5, $point, 0.0001, RLInt_Approx()))
end


text:
# Principal Component Analysis

"""
Linear Principal Component Analysis
"""
struct PCA{T<:Real} <: LinearDimensionalityReduction
    mean::AbstractVector{T}     # sample mean: of length d (mean can be empty, which indicates zero mean)
    proj::AbstractMatrix{T}     # projection matrix: of size d x p
    prinvars::AbstractVector{T} # principal variances: of length p
    tprinvar::T                 # total principal variance, i.e. sum(prinvars)
    tvar::T                     # total input variance
end

## constructor

function PCA(mean::AbstractVector{T}, proj::AbstractMatrix{T}, pvars::AbstractVector{T}, tvar::T) where {T<:Real}
    d, p = size(proj)
    (isempty(mean) || length(mean) == d) ||
        throw(DimensionMismatch("Dimensions of mean and projection matrix are inconsistent."))
    length(pvars) == p ||
        throw(DimensionMismatch("Dimensions of projection matrix and principal variables are inconsistent."))
    tpvar = sum(pvars)
    tpvar <= tvar || isapprox(tpvar,tvar) || throw(ArgumentError("principal variance cannot exceed total variance."))
    PCA(mean, proj, pvars, tpvar, tvar)
end

## properties
"""
    size(M)

Returns a tuple with the dimensions of input (the dimension of the observation space)
and output (the dimension of the principal subspace).
"""
size(M::PCA) = size(M.proj)

"""
    mean(M::PCA)

Returns the mean vector (of length `d`).
"""
mean(M::PCA) = fullmean(size(M.proj,1), M.mean)

"""
    projection(M::PCA)

Returns the projection matrix (of size `(d, p)`). Each column of the projection matrix corresponds to a principal component.
The principal components are arranged in descending order of the corresponding variances.
"""
projection(M::PCA) = M.proj

"""
    eigvecs(M::PCA)

Get the eigenvalues of the PCA model `M`.
"""
eigvecs(M::PCA) = projection(M)

"""
    principalvars(M::PCA)

Returns the variances of principal components.
"""
principalvars(M::PCA) = M.prinvars
principalvar(M::PCA, i::Int) = M.prinvars[i]

"""
    eigvals(M::PCA)

Get the eigenvalues of the PCA model `M`.
"""
eigvals(M::PCA) = principalvars(M)

"""
    tprincipalvar(M::PCA)

Returns the total variance of principal components, which is equal to `sum(principalvars(M))`.
"""
tprincipalvar(M::PCA) = M.tprinvar

"""
    tresidualvar(M::PCA)

Returns the total residual variance.
"""
tresidualvar(M::PCA) = M.tvar - M.tprinvar

"""
    var(M::PCA)

Returns the total observation variance, which is equal to `tprincipalvar(M) + tresidualvar(M)`.
"""
var(M::PCA) = M.tvar

"""
    r2(M::PCA)
    principalratio(M::PCA)

Returns the ratio of variance preserved in the principal subspace, which is equal to `tprincipalvar(M) / var(M)`.
"""
r2(M::PCA) = M.tprinvar / M.tvar
const principalratio = r2

"""
    loadings(M::PCA)

Returns model loadings, i.e. the weights for each original variable when calculating the principal component.
"""
loadings(M::PCA) = sqrt.(principalvars(M))' .* projection(M)

## use

"""
    predict(M::PCA, x::AbstractVecOrMat{<:Real})

Given a PCA model `M`, retur transform observations `x` into principal components space, as

\$\\mathbf{y} = \\mathbf{P}^T (\\mathbf{x} - \\boldsymbol{\\mu})\$

Here, `x` can be either a vector of length `d` or a matrix where each column is an observation,
and `\\mathbf{P}` is the projection matrix.
"""
predict(M::PCA, x::AbstractVecOrMat{T}) where {T<:Real} = transpose(M.proj) * centralize(x, M.mean)

"""
    reconstruct(M::PCA, y::AbstractVecOrMat{<:Real})

Given a PCA model `M`, returns a (approximately) reconstructed observations
from principal components space, as

\$\\tilde{\\mathbf{x}} = \\mathbf{P} \\mathbf{y} + \\boldsymbol{\\mu}\$

Here, `y` can be either a vector of length `p` or a matrix where each column
gives the principal components for an observation, and \$\\mathbf{P}\$ is the projection matrix.
"""
reconstruct(M::PCA, y::AbstractVecOrMat{T}) where {T<:Real} = decentralize(M.proj * y, M.mean)

## show & dump

function show(io::IO, M::PCA)
    idim, odim = size(M)
    print(io, "PCA(indim = $idim, outdim = $odim, principalratio = $(r2(M)))")
end

#### PCA Training

## auxiliary

const default_pca_pratio = 0.99

function check_pcaparams(d::Int, mean::AbstractVector, md::Int, pr::Real)
    isempty(mean) || length(mean) == d ||
        throw(DimensionMismatch("Incorrect length of mean."))
    md >= 1 || error("`maxoutdim` parameter must be a positive integer.")
    0.0 < pr <= 1.0 || throw(ArgumentError("principal ratio must be a positive real value ≤ 1.0."))
end

function choose_pcadim(v::AbstractVector{T}, ord::Vector{Int}, vsum::T, md::Int,
                       pr::Real) where {T<:Real}
    md = min(length(v), md)
    k = 1
    a = v[ord[1]]
    thres = vsum * convert(T, pr)
    while k < md && a < thres
        a += v[ord[k += 1]]
    end
    return k
end

## core algorithms
"""
    pcacov(C, mean; ...)

Compute and return a PCA model based on eigenvalue decomposition of a given covariance matrix `C`.

**Parameters:**
- `C`: The covariance matrix of the samples.
- `mean`: The mean vector of original samples, which can be a vector of length `d`,
           or an empty vector `Float64[]` indicating a zero mean.

*Note:* This function accepts two keyword arguments: `maxoutdim` and `pratio`.
"""
function pcacov(C::AbstractMatrix{T}, mean::AbstractVector{T};
                maxoutdim::Int=size(C,1),
                pratio::Real=default_pca_pratio) where {T<:Real}

    check_pcaparams(size(C,1), mean, maxoutdim, pratio)
    Eg = eigen(Symmetric(C))
    ev = Eg.values
    ord = sortperm(ev; rev=true)
    vsum = sum(ev)
    k = choose_pcadim(ev, ord, vsum, maxoutdim, pratio)
    v, P = extract_kv(Eg, ord, k)
    PCA(mean, P, v, vsum)
end

"""
    pcasvd(Z, mean, tw; ...)

Compute and return a PCA model based on singular value decomposition of a centralized sample matrix `Z`.

**Parameters:**
- `Z`: a matrix of centralized samples.
- `mean`: The mean vector of the **original** samples, which can be a vector of length `d`,
          or an empty vector `Float64[]` indicating a zero mean.
- `n`: a number of samples.

*Note:* This function accepts two keyword arguments: `maxoutdim` and `pratio`.
"""
function pcasvd(Z::AbstractMatrix{T}, mean::AbstractVector{T}, n::Real;
                maxoutdim::Int=min(size(Z)...),
                pratio::Real=default_pca_pratio) where {T<:Real}

    check_pcaparams(size(Z,1), mean, maxoutdim, pratio)
    Svd = svd(Z)
    v = Svd.S::Vector{T}
    U = Svd.U::Matrix{T}
    for i = 1:length(v)
        @inbounds v[i] = abs2(v[i]) / n
    end
    ord = sortperm(v; rev=true)
    vsum = sum(v)
    k = choose_pcadim(v, ord, vsum, maxoutdim, pratio)
    si = ord[1:k]
    PCA(mean, U[:,si], v[si], vsum)
end

## interface functions
"""
    fit(PCA, X; ...)

Perform PCA over the data given in a matrix `X`. Each column of `X` is an **observation**.

**Keyword arguments**

- `method`: The choice of methods:
    - `:auto`: use `:cov` when `d < n` or `:svd` otherwise (*default*).
    - `:cov`: based on covariance matrix decomposition.
    - `:svd`: based on SVD of the input data.
- `maxoutdim`: The output dimension, i.e. dimension of the transformed space (*min(d, nc-1)*)
- `pratio`: The ratio of variances preserved in the principal subspace (*0.99*)
- `mean`: The mean vector, which can be either of
    - `0`: the input data has already been centralized
    - `nothing`: this function will compute the mean (*default*)
    - a pre-computed mean vector

**Notes:**

- The output dimension `p` depends on both `maxoutdim` and `pratio`, as follows. Suppose
  the first `k` principal components preserve at least `pratio` of the total variance, while the
  first `k-1` preserves less than `pratio`, then the actual output dimension will be \$\\min(k, maxoutdim)\$.

- This function calls [`pcacov`](@ref) or [`pcasvd`](@ref) internally, depending on the choice of method.
"""
function fit(::Type{PCA}, X::AbstractMatrix{T};
             method::Symbol=:auto,
             maxoutdim::Int=size(X,1),
             pratio::Real=default_pca_pratio,
             mean=nothing) where {T<:Real}

    @assert !SparseArrays.issparse(X) "Use Kernel PCA for sparse arrays"

    d, n = size(X)

    # choose method
    if method == :auto
        method = d < n ? :cov : :svd
    end

    # process mean
    mv = preprocess_mean(X, mean)

    # delegate to core
    if method == :cov
        C = covm(X, isempty(mv) ? 0 : mv, 2)
        M = pcacov(C, mv; maxoutdim=maxoutdim, pratio=pratio)
    elseif method == :svd
        Z = centralize(X, mv)
        M = pcasvd(Z, mv, n; maxoutdim=maxoutdim, pratio=pratio)
    else
        throw(ArgumentError("Invalid method name $(method)"))
    end

    return M::PCA
end


text:
function nnComputeOutputs(nn,datax;dIdx=1:size(datax,2))
  out = zeros(nn.nNeurons[end],length(dIdx));
  for d=dIdx
    out[:,d] = nnForward(nn,datax[:,d]);
  end

  return out;
end

max_stars_repo_name: krometis/SimpleNeuralNet

max_issues_repo_name: krometis/SimpleNeuralNet

max_forks_repo_name: krometis/SimpleNeuralNet


text:
using Test
using TestSetExtensions
using Logging

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

using .RationalFunctionFields

# Ask me about this
AA = RationalFunctionFields.AbstractAlgebra
Sing = RationalFunctionFields.Singular
Nemo = RationalFunctionFields.Nemo
Primes = RationalFunctionFields.Primes

@info "Testing started"

@testset "All the tests" begin

    @includetests ARGS

end

@info "All tests OK"


text:
const WorldCache = Dict{Tuple{Int, Int}, Int}
const MAXWORLDTYPE = 3

"""
Looks for a tile and generates it if not existing.
"""
function getTile(x::Int, y::Int, db=getDB())
    data = DBInterface.execute(db, """
        SELECT * FROM worldmap WHERE x = $x AND y = $y
    """) |> DataFrame
    try
        return data.type[1]
    catch
        @debug "Tile doesn't exist, generating."
        r = rand(0:MAXWORLDTYPE)
        DBInterface.execute(db, """
        INSERT INTO worldmap (x, y, type)
        VALUES ($x, $y, $r)
        """)
        return getTile(x, y)
    end
end
"""
Caches the tile.
There is no force option, as tiles don't change.
"""
function getTile(x::Int, y::Int, cache::WorldCache)
    if (x, y) in cache
        cache[(x, y)]
    else
        v = getTile(x, y)
        cache[(x, y)] = v
        v
    end
end
getTile(pos::Position, args...) = getTile(pos.x, pos.y, args...)


text:
# Native execution support

export
    @cuda,
    CuCodegenContext

# NOTE: keep this in sync with the architectures supported by NVPTX
#       (see lib/Target/NVPTX/NVPTXGenSubtargetInfo.inc)
const architectures = [
    (v"2.0", "sm_20"),
    (v"2.1", "sm_21"),
    (v"3.0", "sm_30"),
    (v"3.5", "sm_35"),
    (v"5.0", "sm_50") ]

# TODO: allow code generation without actual device/ctx?
# TODO: require passing the codegen context to @cuda, so we can have multiple
#       contexts active, generating for and executing on multiple GPUs
type CuCodegenContext
    ctx::CuContext
    dev::CuDevice
    triple::String
    arch::String

    CuCodegenContext(ctx::CuContext) = CuCodegenContext(ctx, device(ctx))

    function CuCodegenContext(ctx::CuContext, dev::CuDevice)
        # Determine the triple
        if haskey(ENV, "CUDA_FORCE_GPU_TRIPLE")
            triple = ENV["CUDA_FORCE_GPU_TRIPLE"]
        else
            # TODO: detect 64/32
            triple = "nvptx64-nvidia-cuda"
        end

        # Determine the architecture
        if haskey(ENV, "CUDA_FORCE_GPU_ARCH")
            arch = ENV["CUDA_FORCE_GPU_ARCH"]
        else
            cap = capability(dev)
            if cap < architectures[1][1]
                error("No support for SM < $(architectures[1][1])")
            end
            for i = 2:length(architectures)
                if cap < architectures[i][1]
                    arch = architectures[i-1][2]
                    break
                end
            end
        end

        ccall(:jl_init_codegen_ptx, Void, (String, String), triple, arch)

        new(ctx, dev, triple, arch)
    end
end

function destroy(ctx::CuCodegenContext)
    # TODO: we don't look at ctx, but we should, in order to prevent the user
    #       from keeping a ctx after destruction, or destroying one context
    #       while expecting to destroy another (cfr. other destroy methods,
    #       asserting obj==active_obj)
    ccall(:jl_destroy_codegen_ptx, Void, ())
end

#
# macros/functions for native Julia-CUDA processing
#

immutable TypeConst{val} end
value{val}(::Type{TypeConst{val}}) = val

macro cuda(config::Expr, callexpr::Expr)
    # Sanity checks
    if config.head != :tuple || !(2 <= length(config.args) <= 3)
        error("first argument to @cuda should be a tuple (gridDim, blockDim, [shmem])")
    end
    if length(config.args) == 2
        push!(config.args, :0)
    end
    if callexpr.head != :call
        error("second argument to @cuda should be a function call")
    end
    kernel_func_sym = callexpr.args[1]
    if !isa(kernel_func_sym, Symbol)
        # NOTE: cannot support Module.Kernel calls, because TypeConst only works
        #       on symbols (why?). If it allowed Expr's or even Strings (through
        #       string(expr)), we could specialise the stagedfunction on that.
        error("only simple function calls are supported")
    end
    if search(string(kernel_func_sym), "kernel_").start != 1
        error("kernel function should start with \"kernel_\"")
    end

    # HACK: wrap the function symbol in a type, so we can specialize on it in
    #       the staged function
    kernel_func_const = TypeConst{kernel_func_sym}()
    esc(Expr(:call, CUDA.generate_launch, config, kernel_func_const,
             callexpr.args[2:end]...))
end

# TODO: is this necessary? Just check the method cache?
func_cache = Dict{(Symbol, Tuple), CuFunction}()

immutable ArgRef
    typ::Type
    ref::Union(Symbol, Expr)
end

# TODO: nested stagedfunction?

# Read the type and symbol information for each argument of the stagedfunction
function read_arguments(argspec::Tuple)
    args = Array(ArgRef, length(argspec))

    for i in 1:length(argspec)
        args[i] = ArgRef(argspec[i], :( argspec[$i] ))

        if args[i].typ <: DevicePtr{Void} || isbits(args[i].typ)
            # these will be passed as-is
        elseif args[i].typ <: CuArray
            # known conversion to CuDeviceArray
        elseif args[i].typ <: CuManaged
            # should be fine -- will be checked in manage_arguments()
        else
            # TODO: warn optionally?
            # TODO: also display variable name, if possible?
            @warn("you passed an unmanaged argument -- assuming input/output (costly!)")
            args[i] = ArgRef(CuInOut{args[i].typ}, :( CuInOut($(args[i].ref)) ))
        end
    end

    return args
end

# Replace CuManaged-wrapped arguments by the underlying argument it manages,
# managing data up- and downloading along the way
function manage_arguments(args::Array{ArgRef})
    setup = Expr[]
    managed_args = Array(ArgRef, length(args))
    teardown = Expr[]

    # TODO: accessing args[i].ref reconstructs the managed container
    #       avoid by stripping Cu*() and accessing the underlying var?

    for i in 1:length(args)
        if args[i].typ <: CuManaged
            input = (args[i].typ <: CuIn) || (args[i].typ <: CuInOut)
            output = (args[i].typ <: CuOut) || (args[i].typ <: CuInOut)

            managed_type = eltype(args[i].typ)
            if managed_type <: Array
                # set-up
                managed_var = gensym("arg$(i)")
                if input
                    # TODO: CuArray shouldn't auto-upload! Use to_device or smth
                    #       maybe upload/download(CuArray)
                    push!(setup, :( $managed_var = CuArray($(args[i].ref).data) ))
                else
                    # create without initializing
                    push!(setup, :( $managed_var =
                        CuArray($(eltype(managed_type)),
                                size($(args[i].ref).data)) ))
                end

                # TODO: N=1 -- how to / does it support higher dimensions?
                #       does this even make sense?
                managed_args[i] = ArgRef(CuArray{eltype(managed_type), 1},
                                         managed_var)

                # tear-down
                if output
                    data_var = gensym("arg$(i)_data")
                    append!(teardown, [
                        :( $data_var = to_host($managed_var) ),
                        :( copy!($(args[i].ref).data, $data_var) ) ])
                end
                push!(teardown, :(free($managed_var)))
            elseif isbits(managed_type)
                error("managed bits types are not supported -- use an array instead")
            else
                # TODO: support more types, for example CuOut(status::Int)
                error("invalid managed type -- cannot handle $managed_type")
            end
        else
            managed_args[i] = args[i]
        end
    end

    return (setup, managed_args, teardown)
end

# Convert the arguments to be used in combination with @cucall
function convert_arguments(args::Array{ArgRef})
    converted_args = Array(ArgRef, length(args))

    for i in 1:length(args)
        if args[i].typ <: DevicePtr || isbits(args[i].typ)
            # pass these as-is
            converted_args[i] = args[i]
        elseif args[i].typ <: CuArray
            # TODO: pass a CuDeviceArray-wrapped pointer,
            #       so our type matches the actual argument?
            # NOTE: when wrapping ptr in a CuDeviceArray,
            #       the type would now be Void because of Ptr{Void}
            #       (no constructor accepting a pointer?)
            converted_args[i] = ArgRef(CuDeviceArray{eltype(args[i].typ)},
                                      :( $(args[i].ref).ptr ))
        else
            error("invalid argument type -- cannot handle $(args[i].typ)")
        end

    end

    return converted_args
end

# Construct the necessary argument conversions for launching a PTX kernel
# with given Julia arguments
stagedfunction generate_launch(config::(CuDim, CuDim, Int),
                               func_const::TypeConst, argspec::Any...)
    exprs = Expr(:block)

    # Process the arguments
    args = read_arguments(argspec)
    (managing_setup, managed_args, managing_teardown) = manage_arguments(args)
    kernel_args = convert_arguments(managed_args)

    # Compile the function
    kernel_specsig = tuple([arg.typ for arg in kernel_args]...)
    kernel_func_sym = value(func_const)
    if haskey(func_cache, (kernel_func_sym, kernel_specsig))
        ptx_func = func_cache[kernel_func_sym, kernel_specsig]
    else
        # trigger function compilation
        kernel_func = eval(:($(current_module()).$kernel_func_sym))
        try
            precompile(kernel_func, kernel_specsig)
        catch err
            print("\n\n\n*** Compilation failed ***\n\n")
            # this is most likely caused by some boxing issue, so dump the ASTs
            # to help identifying the offending variable
            print("-- lowered AST --\n\n", code_lowered(kernel_func, kernel_specsig), "\n\n")
            #print("-- typed AST --\n\n", code_typed(kernel_func, kernel_specsig), "\n\n")
            throw(err)
        end

        # check if the function actually exists, by mimicking code_llvm()
        # (precompile silently ignores invalid func/specsig combinations)
        # TODO: is there a more clean way to check this?
        kernel_llvm = ccall(:jl_dump_function, Any,
              (Any,Any,Bool,Bool),
              kernel_func, kernel_specsig, false, false)::ByteString
        if kernel_llvm == ""
            error("no method found for $kernel_func$kernel_specsig")
        end
        @debug("LLVM function IR: $(kernel_llvm)")

        # trigger module compilation
        module_ptx = ccall(:jl_to_ptx, Any, ())::String
        @debug("PTX module contents: $(module_ptx)")

        # create CUDA module
        ptx_mod = try
            # TODO: what with other kernel calls? entirely new module? or just
            # add them?
            CuModule(module_ptx)
        catch err
            if isa(err, CuError) && err.code == 209
                # CUDA_ERROR_NO_BINARY_FOR_GPU (#209) usually indicates the PTX
                # code was invalid, so try to assembly using "ptxas" manually in
                # order to get some more information
                try
                    readall(`ptxas`)
                    (path, io) = mktemp()
                    print(io, module_ptx)
                    close(io)
                    run(`ptxas --gpu-name=$(get(cgctx).arch) $path`)
                    rm(path)
                end
            end
            throw(err)
        end

        # Get internal function name
        # FIXME: just get new module / clean slate for every CUDA
        # module, with a predestined function name for the kernel? But
        # then what about type specialization?
        ptx_func_name = ccall(:jl_dump_function_name, Any, (Any, Any),
                              kernel_func, kernel_specsig)
        @debug("PTX function name: $(ptx_func_name)")

        # Get CUDA function object
        ptx_func = CuFunction(ptx_mod, ptx_func_name)

        # Cache result to avoid unnecessary compilation
        func_cache[(kernel_func_sym, kernel_specsig)] = ptx_func
    end

    # Throw everything together and generate the final runtime call
    append!(exprs.args, managing_setup)
    kernel_arg_expr = Expr(:tuple, [arg.ref for arg in kernel_args]...)
    push!(exprs.args, :( CUDA.exec(config, $ptx_func, $kernel_arg_expr) ))
    append!(exprs.args, managing_teardown)

    exprs
end

# Perform the run-time API calls launching the kernel
function exec(config::(CuDim, CuDim, Int), ptx_func::CuFunction, args::(Any...))
    grid  = config[1]
    block = config[2]
    shared_bytes = config[3]

    # Launch the kernel
    launch(ptx_func, grid, block, args, shmem_bytes=shared_bytes)
end


text:
# =======================================
# Functionality to project the flowpipe
# =======================================
import LazySets.Approximations: project

# add a "time" variable by taking the cartesian product of the flowpipe ℱ with each time lapse
function add_time(ℱ::Vector{ReachSet{Zonotope{Float64}}})
    ℱ_with_time = Vector{ReachSet{Zonotope{Float64}}}(undef, length(ℱ))
    @inbounds for i in eachindex(ℱ)
        t0, t1 = ℱ[i].t_start, ℱ[i].t_end
        radius = (t1 - t0)/2.0
        Xi = ℱ[i].X × Zonotope([t0 + radius], hcat(radius))
        Xi = convert(Zonotope, Xi)
        ℱ_with_time[i] = ReachSet{Zonotope{Float64}}(Xi, t0, t1)
    end
    return ℱ_with_time
end

function project(sol::ReachSolution{Zonotope{Float64}})
    N = length(sol.Xk)  # number of reach sets
    n = dim(first(sol.Xk).X) # state space dimension
    options = copy(sol.options)
    πℱ = Vector{ReachSet{Zonotope{Float64}}}(undef, N) # preallocated projected reachsets
    πvars = sol.options[:plot_vars] # variables for plotting
    @assert length(πvars) == 2

    if 0 ∈ πvars
        # add the time variable to the flowpipe (assuming it's not already
        # part of the model)
        ℱ = add_time(sol.Xk)
        n += 1
        options[:n] += 1 # TODO : remove when option is removed
        πvars = copy(πvars)
        πvars[first(indexin(0, πvars))] = n # time index is added in the end
    else
        ℱ = sol.Xk
    end

    M = sparse([1, 2], πvars, [1.0, 1.0], 2, n)
    for i in eachindex(ℱ)
        t0, t1 = ℱ[i].t_start, ℱ[i].t_end
        πℱ_i = linear_map(M, ℱ[i].X)
        πℱ_i = Zonotope(πℱ_i.center, πℱ_i.generators)
        πℱ_i = reduce_order(πℱ_i, options[:max_order])
        πℱ[i] = ReachSet{Zonotope{Float64}}(πℱ_i, t0, t1)
    end
    return ReachSolution(πℱ, options)
end


text:
importall Bukdu
import Base: show

type MarkdownController <: ApplicationController
end

function index(c::MarkdownController)
    code = c[:params]["code"]
    render(Markdown, """
```
julia> $code
$(eval(parse(code)))
```
""")
end

layout(::Layout, body) = "$body"
show(::MarkdownController) = render(Markdown/Layout, "`cool`")

Router() do
    get("/:code", MarkdownController, index)
    get("/mark/down", MarkdownController, show)
end

using Base.Test
conn = (Router)(get, "/mark/down")
@test 200 == conn.status
@test "
cool
" == conn.resp_body

conn = (Router)(get, "/1+2")
@test 200 == conn.status
@test """
julia> 1+2\n3
""" == conn.resp_body

conn = (Router)(get, "/1+2")
@test 200 == conn.status
@test """
julia> 1+2\n3
""" == conn.resp_body

conn = (Router)(get, "/")
@test 200 == conn.status
@test """
julia> \nnothing
""" == conn.resp_body

Logger.have_color(false)
let oldout = STDERR
   rdout, wrout = redirect_stdout()

conn = (Router)(get, "/undefined_variable")
@test 400 == conn.status
@test "bad request UndefVarError(:undefined_variable)" == conn.resp_body

   reader = @async readstring(rdout)
   redirect_stdout(oldout)
   close(wrout)

@test startswith(wait(reader), "ERROR")
end


text:
# render a block of a RegressionTable
function render(io::IO, block::Array{String, 2}, colWidths::Vector{Int64}, align::String, settings::RenderSettings = asciiSettings; isHeader::Bool = false)

    c = size(block,2)

    if length(colWidths) != c
        error("colWidths has invalid length.")
    end
    if length(align) != c
        error("align string has invalid length.")
    end

    # print the whole thing
    for row = 1:size(block,1)
        s = settings.linestart
        for col = 1:c
            # if the string is too long, truncate it
            # (this sometimes happens with column headers)
            # if length(block[row,col]) > colWidths[col]
            #     printstring = block[row,col][1:colWidths[col]]
            # else
            # end
            printstring = block[row,col]

            if align[col] == 'l'
                s = s * rpad(printstring,colWidths[col])
            elseif align[col] == 'r'
                s = s * lpad(printstring,colWidths[col])
            elseif align[col] == 'c'
                l = iseven(colWidths[col]-length(printstring)) ? Int64((colWidths[col]-length(printstring))/2) : Int64((colWidths[col]-length(printstring)+1)/2)
                r = iseven(colWidths[col]-length(printstring)) ? Int64((colWidths[col]-length(printstring))/2) : Int64((colWidths[col]-length(printstring)-1)/2)
                # if the printstring is too long, so be it
                l = max(l,0)
                r = max(r,0)
                s = s * (" " ^ l) * printstring * (" " ^ r)
            end
            if col < c
                s = s * (isHeader ? settings.headercolsep : settings.colsep)
            end
        end
        s = s * settings.linebreak
        println(io, s)
    end

end

# render a whole table
function render(io::IO, tab::RegressionTable, align::String, settings::RenderSettings)

    c = columns(tab)

    if length(align) != c
        error("align string has invalid length.")
    end

    # construct, but not print, the header
    colWidths = column_widths(tab, align)
    adjusted = true

    global header_vec = []

    while adjusted
      header_vec = map(1:size(tab.header, 1)) do i
        h = header(tab, tab.header[i:i,:], settings, colWidths)
        @unpack adjusted, h = adjust_widths!(colWidths, h, settings)
        @unpack headerArray, hr = headerrule(h, settings)

        (h=h, hr=hr, headerArray=headerArray)
      end
    end
    # START RENDERING

    # header
    println(io, settings.header(c,align))

    # headerString
    # if tab.headerString == ""
    #     # don't print anything
    # elseif length(tab.headerString) == 1
    #     println(io, tab.headerString ^ (sum(colWidths) + (columns(tab)-1)*length(settings.colsep)))
    # else
    #     println(io, tab.headerString)
    # end

    # print toprule
    if length(settings.toprule)==1
        # one character, extend over the whole line
        println(io, settings.toprule ^ (sum(colWidths) + (columns(tab)-1)*length(settings.colsep))  )
    else
        println(io, settings.toprule)
    end

    # header
    map(header_vec) do head
        @unpack headerArray, h, hr = head
        render(io, headerArray, h.headerWidths, ("c" ^ size(headerArray,2)), settings, isHeader=true)
        if hr.print_headerrule_separately
            println(io, hr.headerrule)
        end
    end

    # bodies
    for b = 1:size(tab.bodies,1)

        render(io, tab.bodies[b], colWidths, align, settings)

        # if we're not at the last block, print the midrule
        if b < size(tab.bodies,1)
            if length(settings.midrule)==1
                # one character, extend over the whole line
                println(io, settings.midrule ^ (sum(colWidths) + (columns(tab)-1)*length(settings.colsep))  )
            else
                println(io, settings.midrule)
            end
        end
    end

    # print bottomrule
    if length(settings.bottomrule)==1
        # one character, extend over the whole line
        println(io, settings.bottomrule ^ (sum(colWidths) + (columns(tab)-1)*length(settings.colsep))  )
    else
        println(io, settings.bottomrule)
    end

    # footerString
    # if tab.footerString == ""
    #     # don't print anything
    # elseif length(tab.footerString) == 1
    #     println(io, tab.footerString ^ (sum(colWidths) + (columns(tab)-1)*length(settings.colsep)))
    # else
    #     println(io, tab.footerString)
    # end

    # footer
    println(io, settings.footer(c,align))

end


text:
populationtracker = [
    SomaticEvolution.ModuleTracker(
        [1, 2, 3, 4, 4, 3, 4, 3, 4],
        [0.0, 16.892494137270344, 68.68837372525104, 73.9360343780634, 95.06744895741765, 166.45694063272188, 186.08040859501872, 218.77527717255302, 227.24560639149834],
        Cell[Cell([2, 3, 9, 15], 1), Cell([2, 4], 1), Cell([2, 3, 10], 1), Cell([2, 3, 9, 16], 1)],
        SomaticEvolution.CloneTracker[],
        1,
        0
    )
    SomaticEvolution.ModuleTracker(
        [1, 2, 3],
        [166.45694063272188, 215.96281393553738, 216.40410240545967],
        Cell[Cell([1, 6, 8, 11, 13], 1), Cell([1, 6, 8, 12], 1), Cell([1, 6, 8, 11, 14], 1)],
        SomaticEvolution.CloneTracker[],
        2,
        1
    )
    SomaticEvolution.ModuleTracker(
        [1],
        [218.77527717255302],
        Cell[Cell(Int64[1, 5], 1)],
        SomaticEvolution.CloneTracker[],
        3,
        1
    )
]

μ = 2
clonalmutations = 3
rng = MersenneTwister(12)

@testset "multilevel" begin
    #gets list of all mutations in the population
    @test Set(SomaticEvolution.get_mutationlist(populationtracker)) ==
        Set([1, 2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16])
    @test Set(SomaticEvolution.get_mutationlist(populationtracker[1])) ==
        Set([2, 3, 9, 15, 4, 10, 16])
    mutationlist = SomaticEvolution.get_mutationlist(populationtracker)

    #check mean is approx correct for poisson distributed mutations
    expandedmutationids = SomaticEvolution.get_expandedmutationids(3, collect(1:100000), clonalmutations, rng, fixedmu=false)
    @test  mean(reduce(vcat,map(length,values(expandedmutationids)))) ≈ 3 atol=0.01

    module1 = deepcopy(populationtracker[2])
    expandedmutationids = Dict(1=>[4], 2=>[5,6], 3=>[7], 6=>[8,9,10], 8=>[11], 11=>[], 12=>[], 13=>[12,13,14], 14=>[15,16])
    SomaticEvolution.expandmutations!(module1, expandedmutationids, clonalmutations)
    @test sort(clonal_mutation_ids(module1)) == [1, 2, 3, 4, 8, 9, 10, 11]
    @test sort(module1.cells[1].mutations) == [1, 2, 3, 4, 8, 9, 10, 11, 12, 13, 14]
    @test sort(module1.cells[2].mutations) == [1, 2, 3, 4, 8, 9, 10, 11]
    @test sort(module1.cells[3].mutations) == [1, 2, 3, 4, 8, 9, 10, 11, 15, 16]

    expandedmutationids = SomaticEvolution.get_expandedmutationids(μ, mutationlist, clonalmutations, rng, fixedmu=true)
    @test Set(mutationlist) == keys(expandedmutationids)
    @test all(reduce(vcat,map(length,values(expandedmutationids))) .== 2)
    SomaticEvolution.processresults!(populationtracker, μ, clonalmutations, rng, fixedmu=true)
    @test clonal_mutations(populationtracker[1]) == 2 * 1 + 3
    @test length(populationtracker[1].cells[1].mutations) == 2 * 1 + 3 + 2 * 3
    @test length(populationtracker[1].cells[2].mutations) == 2 * 1 + 3 + 2 * 1
    @test length(populationtracker[3].cells[1].mutations) == 2 * 2 + 3

end

SomaticEvolution.ModuleTracker(
    [1, 2, 3, 4, 4, 3, 4, 3, 4],
    [0.0, 16.892494137270344, 68.68837372525104, 73.9360343780634, 95.06744895741765, 166.45694063272188, 186.08040859501872, 218.77527717255302, 227.24560639149834],
    Cell[Cell([2, 3, 9, 15], 1), Cell([2, 4], 1), Cell([2, 3, 10], 1), Cell([2, 3, 9, 16], 1)],
    SomaticEvolution.CloneTracker[],
    1,
    0
)

@testset "singlesim" begin

end


text:
#%%

using Plots
gr()
# %%
xs=-pi:0.1:pi
anim = @animate for i=1:100
    plot(xs,sin.(i*xs));
end
gif(anim, "anim_fps15.gif", fps = 15)
gif(anim, "anim_fps30.gif", fps = 30)

# %%
@gif for i=1:100
    xs=-pi:0.1:pi
    plot(xs,sin.(i*xs));
end every 10

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
World Geodetic Coordinate System of 1984 (WGS 84)
Standardized coordinate system for Earth
Global ellipsoidal reference surface
"""
const WGS84  = OpenStreetMapX.Ellipsoid(a = 6378137.0, f_inv = 298.257223563)
const OSGB36 = OpenStreetMapX.Ellipsoid(a = 6377563.396, b = 6356256.909)
const NAD27  = OpenStreetMapX.Ellipsoid(a = 6378206.4,   b = 6356583.8)

#########################
### s ###
#########################
"""
Point Translator gets longitude
"""
getX(lla::OpenStreetMapX.LLA) = lla.lon
"""
Point Translator gets lattitude
"""
getY(lla::OpenStreetMapX.LLA) = lla.lat
"""
Point Translator gets altitude
"""
getZ(lla::OpenStreetMapX.LLA) = lla.alt
"""
Point Translator gets enu `east` value
"""
getX(enu::OpenStreetMapX.ENU) = enu.east
"""
Point Translator gets enu `north` value
"""
getY(enu::OpenStreetMapX.ENU) = enu.north
"""
Point Translator gets `up` value
"""
getZ(enu::OpenStreetMapX.ENU) = enu.up

"""
    distance(a::ENU, b::ENU)

Calculates a distance between two points `a` and `b`
"""
distance(a::ENU, b::ENU) = OpenStreetMapX.distance(a.east, a.north, a.up,
                                    b.east, b.north, b.up)

"""
    distance(a::ECEF, b::ECEF)

Calculates a distance between two points `a` and `b`
"""
distance(a::ECEF, b::ECEF) = OpenStreetMapX.distance(a.x, a.y, a.z,
                                      b.x, b.y, b.z)

function distance(x1, y1, z1, x2, y2, z2)
    return sqrt((x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2)
end


text:
export OSMPlot, osmplot, osmplot!

##########################################################################################
# OSMPlot recipe
##########################################################################################

"""
Define OSMPlot plotting function with some attribute defaults.

## Arguments

`osm::LightOSM.OSMGraph`

## Keyword arguments

`graphplotkwargs::NamedTuple = (; )` : All kwargs are passed on to `GraphMakie.graphplot!`.
    All kwargs that work with graphplot will also work here (see [GraphMakie docs](https://juliaplots.org/GraphMakie.jl/stable/#The-graphplot-Recipe) for reference).
    Extending the defaults can be done by providing `graphplotkwargs = (; kwargs...)`.
`hide_elabels::Bool = false` : Show or hide edge labels.
`hide_nlabels::Bool = true` : Show or hide node labels.
`buildings::Union{Dict{Integer, LightOSM.Building}, Nothing} = nothing` : Buildings polygons
    are plotted if this is not nothing.
`buildingskwargs::NamedTuple = (; )` : All kwargs are passed on to `Makie.poly` and will
    overwrite the default plotting behaviour
`inspect_nodes::Bool = false` : Enables/disables inspection of OpenStreetMap nodes.
`inspect_edges::Bool = true` : Enables/disables inspection of OpenStreetMap ways.
"""
@recipe(OSMPlot, osm) do scene
    Attributes(
        # general
        graphplotkwargs = NamedTuple(),
        hide_elabels = true,
        hide_nlabels = true,
        buildings = nothing,
        buildingskwargs = NamedTuple(),

        # inspection
        inspect_nodes = false,
        inspect_edges = true,

        # internal
        sorted_edges = [],
        index_to_way = Dict(),
    )
end

function Makie.plot!(osmplot::OSMPlot{<:Tuple{<:OSMGraph}})
    osm = osmplot.osm[]

    # Need to make a PR to have an `index_to_way` dict in OSMGraph. For now we create
    # a sorted list of edges (source and destination node ids) as well as a dict that maps
    # edge index to its related way. We need to use edges(::AbstractGraph) to preserve edge order.
    osmplot.sorted_edges = collect([e.src, e.dst] for e in edges(osm.graph))
    osmplot.index_to_way = Dict(
        zip(1:Graphs.ne(osm.graph),
            (osm.highways[way] for way in
             (osm.edge_to_highway[edge] for edge in
              ([osm.index_to_node[s], osm.index_to_node[d]] for (s, d) in osmplot.sorted_edges[])))
        )
    )

    # Node positions
    # OSMGraph.node_coordinates is in lat/lon format. Reversing it provides lon/lat
    # which then creates standard north-oriented maps.
    node_pos = Point2.(reverse.(osm.node_coordinates))

    # OSMMakie defaults (see defaults.jl for details)
    edge_defaults = set_edge_defaults(osmplot)
    node_defaults = set_node_defaults(osmplot, edge_defaults.edge_width, edge_defaults.edge_color)

    # If user provided buildings, plot them as polys below ways layer
    if !isnothing(osmplot.buildings[])
        building_polys = @lift(get_building_polys($(osmplot.buildings)))
        bp = poly!(osmplot, building_polys;
            color = BUILDINGSCOLORS,
            strokecolor = colorant"#444", strokewidth = 0.5,
            osmplot.buildingskwargs...)
        bp.inspectable[] = false # Disable building inspection for now
    end

    # Create the ways layer as a graphplot
    # User-provided graphplotkwargs will overwrite defaults
    gp = graphplot!(osmplot, osm.graph;
        layout = _ -> node_pos,
        osmplot.graphplotkwargs...,
        node_defaults...,
        edge_defaults...
    )

    # Setup inspection on mouse hover
    # Tracks inspect_edges for ways
    gp.plots[1].plots[1].inspectable = lift(identity, osmplot.inspect_edges)
    # Always disabled for one-way arrows
    gp.plots[2].inspectable[] = false
    # Tracks inspect_nodes for nodes
    gp.plots[3].inspectable = lift(identity, osmplot.inspect_nodes)

    return osmplot
end


text:
using MonoMath
using Test
using Unicode

@testset "MonoMath" begin
  # Check that all codepoints used in the samples are valid Unicode
  for (i, s) in enumerate(MonoMath.text_samples)
    @testset "Sample $i" begin
      for (j, c) in enumerate(s)
        @testset "Codepoint $j ($c = $(repr(codepoint(c))))" begin
          @test Unicode.isassigned(c)
        end
      end
    end
  end
end

max_stars_repo_name: waldyrious/mono-math

max_issues_repo_name: waldyrious/mono-math

max_forks_repo_name: waldyrious/mono-math


text:
export LambdaActor, lambda

"""
    LambdaActor{D, N, E, C}(on_next::N, on_error::E, on_complete::C) where D

Lambda actor wraps `on_next`, `on_error`, `on_complete` callbacks for data, error and complete events.
Should not be used explicitly, use [`lambda`](@ref) creation operator instead.

# Constructor arguments
- `on_next`: Callback for data event. Optional. Default is `nothing`.
- `on_error`: Callback for error event. Optional. Default is `nothing`.
- `on_complete`: Callback for complete event. Optional. Default is `nothing`.

See also: [`Actor`](@ref), [`lambda`](@ref)
"""
struct LambdaActor{D, NextCallback, ErrorCallback, CompleteCallback} <: Actor{D}
    on_next     :: NextCallback
    on_error    :: ErrorCallback
    on_complete :: CompleteCallback
end

function on_next!(actor::LambdaActor{D}, data::D) where D
    if actor.on_next !== nothing
        actor.on_next(data)
    end
end

function on_error!(actor::LambdaActor, err)
    if actor.on_error !== nothing
        actor.on_error(err)
    end
end

function on_complete!(actor::LambdaActor)
    if actor.on_complete !== nothing
        actor.on_complete()
    end
end

struct LambdaActorFactory{NextCallback, ErrorCallback, CompleteCallback} <: AbstractActorFactory
    on_next     :: NextCallback
    on_error    :: ErrorCallback
    on_complete :: CompleteCallback
end

function create_actor(::Type{L}, factory::LambdaActorFactory{N, E, C}) where { L, N, E, C }
    return LambdaActor{L, N, E, C}(factory.on_next, factory.on_error, factory.on_complete)
end

"""
    lambda(; on_next = nothing, on_error = nothing, on_complete = nothing)
    lambda(::Type{T}; on_next = nothing, on_error = nothing, on_complete = nothing) where T

Creation operator for the 'LambdaActor' actor.

# Examples

```jldoctest
using Rocket

actor = lambda(Int; on_next = (d) -> println(d))
actor isa LambdaActor{Int}

# output
true
```

See also: [`LambdaActor`](@ref), [`AbstractActor`](@ref)
"""
lambda(; on_next::N = nothing, on_error::E = nothing, on_complete::C = nothing)          where { N <: Union{Nothing, Function}, E <: Union{Nothing, Function}, C <: Union{Nothing, Function} }    = LambdaActorFactory{N, E, C}(on_next, on_error, on_complete)
lambda(::Type{T}; on_next::N = nothing, on_error::E = nothing, on_complete::C = nothing) where { T, N <: Union{Nothing, Function}, E <: Union{Nothing, Function}, C <: Union{Nothing, Function} } = LambdaActor{T, N, E, C}(on_next, on_error, on_complete)


text:
using Test, AdmittanceModels, LinearAlgebra

@testset "TransmissionLine with terminations" begin
    Z0 = 50.0
    tline = TransmissionLine(["0", "1"], 1e8, Z0, 5e-3, δ=50e-6)
    ω = 2π * collect(range(4, stop=6, length=100)) * 1e9
    pso = PSOModel(tline)
    bbox0 = canonical_gauge(Blackbox(ω, pso))
    bbox1 = Blackbox(ω, tline)
    @test isapprox(admittance_matrices(bbox0), admittance_matrices(bbox1), atol=1e-5)
    @test !isapprox(admittance_matrices(bbox0), admittance_matrices(bbox1), atol=1e-6)
    # now terminate it with Z0
    resistor = ParallelComponent("1", 0, 1/Z0, 0)
    casc = Cascade([tline, resistor], [open_ports_except => ["0"]])
    S = [x[1,1] for x in scattering_matrices(Blackbox(ω, PSOModel(casc)), [Z0])]
    @test all([abs(x) < 1e-4 for x in S])
    S = [x[1,1] for x in scattering_matrices(Blackbox(ω, casc), [Z0])]
    @test all([abs(x) < 1e-15 for x in S])
    # now terminate it with Z0/2
    z0 = Z0/2
    resistor = ParallelComponent("1", 0, 1/z0, 0)
    casc = Cascade([tline, resistor], [open_ports_except => ["0"]])
    S = [x[1,1] for x in scattering_matrices(Blackbox(ω, PSOModel(casc)), [Z0])]
    correct_S = ((z0 - Z0)/(z0 + Z0)) * exp.(-2im * ω * tline.len/tline.propagation_speed)
    @test isapprox(S, correct_S, rtol=1e-4)
    S = [x[1,1] for x in scattering_matrices(Blackbox(ω, casc), [Z0])]
    @test isapprox(S, correct_S, rtol=1e-14)
end

@testset "reflection λ/4 resonator" begin
    Z0 = 50.0
    tline = TransmissionLine(["open", "coupler", "short"], 1e8, Z0, 5e-3, locations=[1e-3], δ=50e-6)
    capacitor = SeriesComponent("in", "coupler", 0, 0, 10e-15)
    casc = Cascade([tline, capacitor], [short_ports => ["short"], open_ports_except => ["in"]])
    ω = 2π * collect(range(4, stop=6, length=500)) * 1e9
    bbox0 = canonical_gauge(Blackbox(ω, PSOModel(casc)))
    bbox1 = Blackbox(ω, casc)
    circ = cascade_and_unite(Circuit(tline, "tline"), Circuit(capacitor))
    g = AdmittanceModels.ground
    circ = unite_vertices(circ, g, "short")
    pso = PSOModel(circ, [("in", g)], ["in"])
    bbox2 = canonical_gauge(Blackbox(ω, pso))
    @test isapprox(admittance_matrices(bbox0), admittance_matrices(bbox1), atol=1e-3)
    @test isapprox(admittance_matrices(bbox0), admittance_matrices(bbox2), atol=1e-11)
    @test !isapprox(admittance_matrices(bbox0), admittance_matrices(bbox1), atol=1e-4)
end


text:
#   Unit "hilbert" of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019-2022,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements the Hilbert transform and analytic signal
#   estimations for signals of arbitrary length, using FFTW.

"""
```julia
   (1)
   function analyticsignal( X  :: Union{Vector{T}, Matrix{T}},
                            wl :: Int     = size(X, 1);
                        nonlinear :: Bool    =  false,
                        planner   :: Planner =  getplanner,
                        ⏩       :: Bool    =  true) where T<:Real

   (2)
   function analyticsignal( 𝐗      :: Vector{Matrix{T}},
                            wl     :: Int;
                        nonlinear  ::  Bool      =  false,
                        planner    ::  Planner   =  getplanner,
                        ⏩        ::  Bool      =  true) where T<:Real
```

(1)

Compute the analytic signal(AS) of vector `X` or of
all column vectors of matrix `X` via the FFT and iFFT procedure,
as explained in Marple(1999).
If `wl`=size(X, 1) (default), use the standard method passing to the FFT
and iFFT all samples in `X` altogether, whereas if `wl`<size(X, 1) a
sliding-windows method is used (see below).

Return the analytic signal `𝑌`, a complex vector if `X` is a vector
or a complex matrix holding in its columns the analytic signal of the
columns of `X` if `X` is a matrix. `𝑌` has the same number of samples (rows)
as `X`.

Contrarely to what is done in the [DSP](https://github.com/JuliaDSP/DSP.jl)
package, the DC level of the signal is removed, therefore,
if the input signal features a non-zero DC level,
the real part of the AS will be equal to the input signal with the
DC level removed. The imaginary part of `𝑌` is the Hilbert transform
of such no-DC `X`.

The sliding-windows AS allows an efficient estimation of the AS for vectors and
matrices when they are very large; it proceeds
computing the AS on 50% sliding overlapping windows and forming the AS
by retaining the central half of each window. The number of points effectively
used to obtain the final estimation is ``wl``÷2 (integer division).
`wl` must be even for using this estimation. This procedure
produces edge effects, thus the first and last ``wl÷4`` samples of the AS
estimation should be discarded. In practice, one requires the AS of a
larger data segment and trims at least ``wl÷4`` samples at the beginning and
end of the estimation. This is done automatically
by the [`TFanalyticsignal`](@ref) function.

Also, the sliding-windows AS method creates small discontinuities at sample
``wl``÷4 and then every ``wl``÷2 samples,
therefore ``wl`` should be chosen as large as possible.

!!! note "Nota Bene"

    In order to avoid FFT computation of very long epochs,
    if `wl` > 2^14, then `wl` is set to 2^10. Below this limit, as long as
    the computations are feasable, use the standard method. If you absolutely
    need to use the sliding-windows method, see [window length in FFTW](@ref)
    for setting efficiently argument `wl`.

    The input signal should be previously band-pass or high-pass filtered
    so as not to contain frequency components below the first discrete
    Fourier frequency obtained with windows of `wl` samples, that is,
    below sr/wl Hz, where sr is the sampling rate.

**Optional Keyword Arguments**:

`nonlinear`, if true, the analytic signal is normalized so that its amplitude
is ``1.0`` at all points. This allow non-linear univariate and bivariate estimations
(see [timefrequencyuni.jl](@ref) and [timefrequencybi.jl](@ref)).

`planner` is an instance of the [`Planner`](@ref) object, holding the forward
and backward FFTW plans used to compute the FFTs and the iFFTs.
By default the planner is computed, but it can be passed as an
argumet here if it is pre-computed. This is interesting if the
`analyticsignal` function is to be invoked repeatedly.

if `⏩` is true, the method is run in multi-threaded mode across the series
in `X` if the number of series is at least twice the number of threads Julia
is instructed to use. See [Threads](@ref).

(2)

Compute the analytic signal for all ``k``
multivariate data matrices given as a vector of matrices `𝐗`.
Return a vector of matrices hodling the corresponding analytic signals
as in method (1). The FFT and iFFT plans are computed only once.
The ``k`` matrices in `𝐗` may have different number of columns (i.e., different
number of series) and different number of rows (samples).
However, the number of rows must be larger than `wl` for all of them.

If `⏩` is true, this method run in multi-threaded mode across the
matrices in `𝐗` if the number of matrices is at least twice
the number of threads Julia is instructed to use, otherwise it
tries to run each analytic signal estimation in multi-threaded mode
as per method (1). See [Threads](@ref).

This function is called by the following functions operating
on time-frequency reprsentations: [`TFanalyticsignal`](@ref),
[`TFamplitude`](@ref), [`TFphase`](@ref), [`meanAmplitude`](@ref),
[`concentration`](@ref), [`meanDirection`](@ref), [`comodulation`](@ref),
[`coherence`](@ref).

**References**
Marple S.L. (1999)
Computing the Discrete-Time Analytic Signal via FFT.
IEEE Transactions on Signal Processing 47(9), 2600-3.

**Examples**:
```julia
using FourierAnalysis, FFTW, LinearAlgebra, Statistics, Plots, DSP
t=128; lab=["x", "real(y)", "imag(y)"]

# Analytic signal of one vector
x=sinusoidal(10, 2, 128, t, π/2; DC=10) # cosine
y=analyticsignal(x)
# note that real(y) is x without the DC level, i.e., x=real(y)+DC
plot([x, real(y), imag(y)]; labels=lab)

# make a check
s=sinusoidal(10, 2, 128, t, 0) # sine
norm(s-imag(y)) # should be zero

# Analytic Signal by DSP.jl
y2=hilbert(x)
norm(s-imag(y2)) # should be zero
# DSP.jl does not remove the DC level
# thus x=real(y2) in this case
plot([x, real(y2), imag(y2)]; labels=lab)

# Analytic signal of multiple vectors
x=hcat(x, sinusoidal(10, 3, 128, t, π/2; DC=10))
y=analyticsignal(x)

# sliding-windows analytic signal of one vector
# (note edge effects)
x=sinusoidal(10, 2, 128, t*4, π/2; DC=0)
y=analyticsignal(x, t)
plot([x, real(y), imag(y)]; labels=lab)

# sliding-windows analytic signal of multiple vectors
x=hcat(x, sinusoidal(10, 3, 128, t*4, π/2; DC=0))
y=analyticsignal(x, t)
```
"""
function analyticsignal( X  :: Union{Vector{T}, Matrix{T}},
                         wl :: Int     = size(X, 1);
                  nonlinear :: Bool    =  false,
                  planner   :: Planner =  getplanner,
                  ⏩       :: Bool    =  true) where T<:Real

   # get all needed paramters (see below the _paramHT! function)
   X_, 𝚙, i𝚙, e, tₐₗₗ, n, wl½, two_wl⁻¹, cT, f, g, ζ = _paramHT!(X, wl, planner)

   𝑌=zeros(cT, tₐₗₗ, n) # Preallocate memory for the Analytic Signal of X

   # Add to `𝑌` the analytic signal of a vector `X_` or of the column vectors of
   # matrix `X_`, computed on the epoch(s) of duration `t` samples starting at
   # sample `at`. It does not return anything.
   # ARGUMENTS:
   # `n` (Int) is the number of columns in `X_`; it must be 1 if `X_` is a vector.
   # `wl½` = wl÷2
   # `two_wl⁻¹` = 2/wl
   # `e` (Int) the number of sliding windows
   # `𝚙` is the forward FFTW plan performing the FFT.
   # `i𝚙` is the backward FFTW plan performing the iFFT (Hilbert transform).
   # `ζ`, =zeros(cT, wl-wl½_), a zero-vector appended to the iFFT vectors, since the second half of the FFT is not computed
   # `f`, the lower limit (in samples) of the central half of each Hilbert transform to cumulate with respect to the FFT window
   # `g`, the upper limit (in samples) of the central half of each Hilbert transform to cumulate with respect to the FFT window
   function _analyticsignal!(at, 𝚗)
      #if 𝚗==1 println(at:at+wl-1, " ", f, " ", g, " ", at+f-1:at+g-1) end ###
      y = 𝚙*X_[at:at+wl-1, 𝚗]        # FFT
      @inbounds begin                # transform
         y[1]=0.
         for i=2:wl½ y[i] *= two_wl⁻¹ end
      end
      # iFFT and cumulate
      e==1 ? 𝑌[:, 𝚗]+=i𝚙*vcat(y, ζ) : 𝑌[at+f-1:at+g-1, 𝚗]+=(i𝚙*vcat(y, ζ))[f:g]
   end

   # add the Analytic Signal
   as!(𝚎, 𝚗)=_analyticsignal!((𝚎*wl½)+1, 𝚗)
   _thread(⏩, n) ? (for 𝚎=0:e-1 @threads for 𝚗=1:n as!(𝚎, 𝚗) end end) : (for 𝚗=1:n, 𝚎=0:e-1 as!(𝚎, 𝚗) end)

   # return a vector if `X` is a vector, a matrix if `X` is a matrix.
   # for the sliding-windows version eliminate the first and last
   # t½ samples that have been previously padded.
   𝑌_ = n==1 ? (e>1 ? 𝑌[:][wl½+1:end-wl½] : 𝑌[:]) : (e>1 ? 𝑌[wl½+1:end-wl½, 1:n] : 𝑌)
   if nonlinear @inbounds for i in eachindex(𝑌_) 𝑌_[i]/=abs(𝑌_[i]) end end
   return 𝑌_
end

function analyticsignal( 𝐗      :: Vector{Matrix{T}},
                         wl     :: Int;
                     nonlinear  ::  Bool      =  false,
                     planner    ::  Planner   =  getplanner,
                     ⏩        ::  Bool      =  true) where T<:Real

   cT, k=typeof(Complex(T(0))), length(𝐗)

   planner ≠ getplanner ? plan=planner : plan=Planner(plan_estimate, -1.0, wl, T, true)

   𝒀=Vector{Matrix{cT}}(undef, k)

   # function to add the analytic signal
   as(i, ⏩)=analyticsignal(𝐗[i], wl; nonlinear=nonlinear, planner=plan, ⏩=⏩)
   _thread(⏩, k) ? (@threads for i=1:k 𝒀[i]=as(i, false) end) : (for i=1:k 𝒀[i]=as(i, ⏩) end)

   return 𝒀
end

# Internal function: prepare all parameters that are needed for
# Analytic Signal and sliding-windows Analytic Signal estimations
#
# ARGUMENTS:
# `X` is the input data vector or matrix (n time-series)
# `wl` (Int) is the epoch length for FFT
# `planner` a pre-computed FFT `Planner` or `getPlanner` to compute it
#
# RETURN
# `X_`: if a standard anaytic signal is requested (t==size(X, 1)) then X_=X
#       else X_ is the data in `X` with wl½ zeros prepended and appended.
# `𝚙` the forward FFTW plan. It is computed if argument `planner`=`getplanner`.
# `i𝚙` the backward FFTW plan. It is computed if argument `planner`=`getplanner`.
# `e` (Int) the number of sliding windows, 1 for the standard AS method.
# `tₐₗₗ` (Int) number of samples total. For sliding-windows anaytic signal,
#     this is not size(X, 1) as zeros are padded.
# `n` (Int) 1 if `X` is a Vector, the number of columns of `X` otherwise.
# `wl½` `t`÷2 computed by the right aritmetic bit-wise operation `wl`>>1.
# `two_wl⁻¹` = 2/wl.
# `cT` the complex type corresponding to type `T` in function definition
#     (e.g., ComplexF64 if T=Float64).
# `f`, the lower limit (in samples) of the central half of each
#     Hilbert transform to cumulate with respect to the FFT window.
# `g`, the upper limit (in samples) of the central half of each
#     Hilbert transform to cumulate with respect to the FFT window.
# `ζ`, =zeros(cT, wl-wl½_), a zero-vector append to the iFFT vectors,
#     since the second half of the FFT is not computed
# NB if wl > 2^14 then t is set to 2^10. This affects the default
#     behavior of all AS functions.
function _paramHT!( X       :: AbstractArray{T},
                    wl      :: Int,
                    planner :: Planner) where T<:Real

    size(X, 1) < wl && @error 📌*", the number of samples in input matrix is smaller than the desired window length."
    wl > 2^14 ? wl = 2^10 : nothing
    n, wl½, two_wl⁻¹ = size(X, 2), (wl>>1), T(2/wl)
    if size(X, 1)>wl
        X_=[zeros(T, wl½, n); X; zeros(T, wl½, n)] else X_=X
        isodd(wl) && @error 📌*", for Welch-like Analytic Signal estimation `wl` must be even."
    end
    tₐₗₗ = size(X_, 1) # for sliding-windows anaytic signal, this is not size(X, 1)
    cT = typeof(Complex(T(0)))
    if planner ≠ getplanner
       𝚙=planner.p
       i𝚙=planner.ip
    else
       𝚙=plan_rfft(zeros(T, wl), flags=plan_estimate, timelimit=-1.0)
       i𝚙=plan_bfft(zeros(cT, wl), flags=plan_estimate, timelimit=-1.0)
    end
    f=wl½÷2+1 # lower limit of central region to copy
    g=f+wl½-1 # upper limit of central region to copy
    e = (tₐₗₗ-wl)÷wl½+1 # number of 50% overlapping epochs or 1 if tₐₗₗ=wl
    ζ=zeros(cT, wl-wl½-1)

    return X_, 𝚙, i𝚙, e, tₐₗₗ, n, wl½, two_wl⁻¹, cT, f, g, ζ
end


text:
mutable struct PrincetonDriver <: LaneFollowingDriver
    a::Float64 # predicted acceleration
    σ::Float64 # optional stdev on top of the model, set to zero or NaN for deterministic behavior
    k::Float64 # proportional constant for speed tracking [s⁻¹]
    v_des::Float64 # desired speed [m/s]

    function PrincetonDriver(;
        σ::Float64 = NaN,
        k::Float64 = 1.0,
        v_des::Float64 = 29.0,
        )

        retval = new()
        retval.a = NaN
        retval.σ = σ
        retval.k = k
        retval.v_des = v_des
        retval
    end
end
get_name(::PrincetonDriver) = "PrincetonDriver"
function set_desired_speed!(model::PrincetonDriver, v_des::Float64)
    model.v_des = v_des
    model
end
function track_longitudinal!(model::PrincetonDriver, v_ego::Float64, v_oth::Float64, headway::Float64)

    v_des = model.v_des
    k = model.k

    if !isnan(v_oth)
        v_des = min(v_oth*(1-exp(-k*headway/v_oth - 1)), v_des)
    end

    Δv = v_des - v_ego
    model.a = Δv*k # predicted accel to match target speed

    model
end
function Base.rand(model::PrincetonDriver)
    if isnan(model.σ) || model.σ ≤ 0.0
        LaneFollowingAccel(model.a)
    else
        LaneFollowingAccel(rand(Normal(model.a, model.σ)))
    end
end
function Distributions.pdf(model::PrincetonDriver, a::LaneFollowingAccel)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        pdf(Normal(model.a, model.σ), a.a)
    end
end
function Distributions.logpdf(model::PrincetonDriver, a::LaneFollowingAccel)
    if isnan(model.σ) || model.σ ≤ 0.0
        Inf
    else
        logpdf(Normal(model.a, model.σ), a.a)
    end
end

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
using Test
using GatherShot

@testset "select_tests finds simplest" begin
    @test GatherShot.select_tests(Bool[1]) == [1]
end

@testset "select_tests gets independent" begin
    outcomes = Bool[
        1 0 0;
        0 1 0;
        0 0 1
    ]
    @test GatherShot.select_tests(outcomes) == [1, 2, 3]
end

@testset "select_tests ignores non-dominant" begin
    outcomes = Bool[
        1 0 0;
        1 1 0;
        0 0 1
    ]
    @test GatherShot.select_tests(outcomes) == [2, 3]
end