text:
using BenchmarkTools
using Dates
using JSON
include("../src/SSSTensor.jl")

log_root = "./logs/"
log_folder = log_root *"Contraction-logs/"

if !isdir(log_root)
    mkdir(log_root)
end
if !isdir(log_folder)
    mkdir(log_folder)
end

 scales = 5
 n_0 = 100000
 ord = 3
 tol = 1e-9
 sparsity = 10

root_name = log_folder*
  "Contraction-profiling_$(Dates.month(now()))-$(Dates.day(now()))-"*
   "$(Dates.year(now()))_$(Dates.hour(now())):$(Dates.minute(now()))"

#initialize dictionaries
Trials = BenchmarkGroup()

if false  #fulll contraction tests

Trials["matSST_contract_k_1"] = BenchmarkGroup()
Trials["dictSST_contract_k_1"] = BenchmarkGroup()
Trials["dictSST_contract"] = BenchmarkGroup()

for i = 1:scales
    n = n_0^i

    nnz = Int(floor(sparsity*n))
    m = ord -1

    x = rand(n)
    y = Array{Float64,1}(undef,n)

    indices = rand(1:n,nnz,ord)
    for i =1:nnz
        indices[i,:] = indices[i,sortperm(indices[i, :])]
    end
    vals = rand(nnz)

    A = ssten.SSSTensor([(indices[i,:],vals[i]) for i in 1:nnz],n)

#    Trials["matSST_contract_k_1"]["n:$(n)"]  = @benchmark ssten.contract_k_1($indices,$vals,$n, $x)
    Trials["dictSST_contract_k_1"]["n:$(n)"] =  @benchmark ssten.contract_k_1($A, $x)
#    Trials["dictSST_contract"]["n:$(n)"] =      @benchmark ssten.contract($A, $x, $m)

    for k in keys(Trials)
        println("Results for Trial: $(k) for n = $(n)")
        display(Trials[k]["n:$(n)"])
        print("\n")
    end

end

BenchmarkTools.save(root_name*".json", Trials)
end

if true  #edge wise contraction operations
ord_0 = 3

Trials["row_edge_contract_k_1"] = BenchmarkGroup()
Trials["tuple_contract_k_1"] = BenchmarkGroup()
Trials["tuple_contract"] = BenchmarkGroup()

n = n_0

for i = 1:scales

    ord = ord_0^i

    m = ord -1

    x = rand(n)
    y = zeros(n)

    indices = rand(1:n,ord)
    sort!(indices)
    val = rand()

    e = (indices,val)

    Trials["row_edge_contract_k_1"]["ord:$(ord)"]  = @benchmark ssten.contract_edge_k_1!($indices,$val,$x,$y)
    Trials["tuple_contract_k_1"]["ord:$(ord)"] =  @benchmark ssten.contract_edge_k_1($e,$x)
    Trials["tuple_contract"]["ord:$(ord)"] =      @benchmark ssten.contract_edge($e,$x,$m)

    for k in keys(Trials)
        println("Results for Trial: $(k) for ord = $(ord)")
        display(Trials[k]["ord:$(ord)"])
        print("\n")
    end

end

end

# unprofiled
#contract(A,x,m)
#contract(A::Array{N,k}, x::Array{M,1},m::Int64)
#contract_k_1(A::SSSTensor, x::Array{N,1})
#contract_multi(A::SSSTensor, Vs::Array{N,2})
#contract(A::SSSTensor,v::Array{N,1},u::Array{N,1})
#contract_k_1!(indices::Array{Int,2},nnz::Array{N,1}, x::Array{N,1},y::Array{N,1})


text:
using HTTP
using JSON
using RCall

export get_lat_lon, vgregion_lat_lon_bounds

#nominatim requires as user agent and will block you otherwise
HTTP.setuseragent!("Mozilla/5.0 (Windows NT 10.0; rv:68.0) Gecko/20100101 Firefox/68.0")

const NominatimSearchPrefix = "https://nominatim.openstreetmap.org/search?q="

# Nominatim will block you if you access more than once per second so keep track
# of time.
struct NominatimSearcher
    lasttime::Float64
    NominatimSearcher() = new(0.0)
end

const DefaultSearcher = NominatimSearcher()

"""
    get_lat_lon(address::String)

Returns (lat,lon) coordinate for address (typically in middle of street)
as well as additional data returned from Nominatim.
"""
function get_lat_lon(address::String, n::NominatimSearcher = DefaultSearcher)
    addressFormatted = join(split(address),"+")
    searchurl = NominatimSearchPrefix * addressFormatted * "&format=json&limit=1"

    if time() - n.lasttime < 1.0
        sleep(1.05 - (time() - n.lasttime))
        n.lasttime = time()
    end

    rawdata = HTTP.get(searchurl)

    metadata = JSON.parse(String(rawdata.body))[1]
    coordinates = parse(Float64, metadata["lat"]), parse(Float64, metadata["lon"])

    return coordinates, metadata
end

# Yttergränser för Västra Götaland från Google Maps 2020-04-22:
const VGRegion = Dict(
    :south => (57.140190, 13.127416), # Gislaved V
    :east  => (58.687169, 14.714468), # Karlsborg N
    :north => (59.260856, 12.208255), # Årjäng SO
    :west  => (58.887149, 10.963946)  # väster om Nordkoster
)

function vgregion_lat_lon_bounds(delta = 0.10)
    lats = map(first, values(VGRegion))
    lons = map(t -> t[2], values(VGRegion))

    # Boundingbox for ggmap given as lowerleftlon, lowerleftlat, upperightlon, upperrightlat
    return (minimum(lons)-delta, minimum(lats)-delta, maximum(lons)+delta, maximum(lats)+delta)
end

# google maptypes: roadmap, terrain, satellite, hybrid
# osm maptypes: none...

#get_stamenmap(bbox = c(left = -95.80204, bottom = 29.38048, right =
#-94.92313, top = 30.14344), zoom = 10, maptype = c("terrain",
#"terrain-background", "terrain-labels", "terrain-lines", "toner",
#"toner-2010", "toner-2011", "toner-background", "toner-hybrid",
#"toner-labels", "toner-lines", "toner-lite", "watercolor"),
#crop = TRUE, messaging = FALSE, urlonly = FALSE,
#color = c("color", "bw"), force = FALSE, where = tempdir(), ...)
function make_stamen_map(boundingbox, filename = "t.png";
    maptype = "toner-lite", color="bw")
    R"library(ggmap)"
    left, bottom, right, top = boundingbox
    R"bbox <- c(left = $left, bottom = $bottom, right = $right, top = $top)"
    R"zoomlevel <- calc_zoom(bbox)"
    R"mymap <- get_stamenmap(bbox, zoom=zoomlevel, maptype=$maptype, crop=FALSE, color=$color)"
    R"p <- ggmap(mymap)"
    R"ggsave($filename, p)"
    @info "Saved to file $filename"
    filename
end

ext: jl


text:
using RobotOSData
using Test

include("testmsgs.jl")
include("testbags.jl")


text:
using GeometricAlgebra
using Documenter

DocMeta.setdocmeta!(GeometricAlgebra, :DocTestSetup, quote
	using Revise, GeometricAlgebra
end; recursive=true)

doctest(GeometricAlgebra)


text:
import ColorTypes: RGB

## Code for PIGMENT type and its sons ##################################################

abstract type Pigment end

## Unfirom Pigment

"""
    UniformPigment(color)

It creates a uniform pigment, where color is a RGB type.
"""
struct UniformPigment <: Pigment
    color::RGB

    UniformPigment(color::RGB) = new(color)
end

Base.:≈(pig1::UniformPigment, pig2::UniformPigment) = pig1.color ≈ pig2.color

"""
    get_color(un_pig, Vec2D)
    get_color(im_pig, Vec2D)
    get_color(check_pig, Vec2D)

It returns the **COLOR** of the pixel with coordinates (u,v) of the Vec2D given. It works with every pigment implemented
(uniform, image, checkered).
"""
get_color(un_pig::UniformPigment, vec2d::Vec2D) = return un_pig.color

## Image Pigment

"""
    ImagePigment(img)

It creates a Image Pigment, where img is HdrImage type. Thus, the image is transferred to the pixels.
"""
struct ImagePigment <: Pigment
    image::HdrImage

    ImagePigment(img::HdrImage) = new(img)
end

function get_color(im_pig:: ImagePigment, uv::Vec2D)
    col = convert(Int64, floor(uv.u * im_pig.image.width +1))
    row = convert(Int64, floor(uv.v * im_pig.image.height +1))

    if col >= im_pig.image.width
        col = im_pig.image.width
    end

    if row >= im_pig.image.height
        row = im_pig.image.height
    end

    return im_pig.image.pixels[get_pixel(im_pig.image, col, row)]
end

## Checkered Pigment

"""
    CheckeredPigment(color1, color2, n_steps)

It creates a Checkered Pigment, useful for debugging.
If not defined n_steps = 10.
"""
struct CheckeredPigment <: Pigment
    color1::RGB
    color2::RGB
    n_steps::Number

    CheckeredPigment(color1, color2, n_steps=10) = new(color1, color2, n_steps)
end

Base.:≈(pig1::CheckeredPigment, pig2::CheckeredPigment) = pig1.color1 ≈ pig2.color1 && pig1.color2 ≈ pig2.color2 && pig1.n_steps ≈ pig2.n_steps

function get_color(check_pig::CheckeredPigment, uv::Vec2D)
    int_u = convert(Int64, floor(uv.u * check_pig.n_steps))
    int_v = convert(Int64, floor(uv.v * check_pig.n_steps))

    if int_u % 2 == int_v % 2
        return check_pig.color1
    else
        return check_pig.color2
    end
end

## Code for BRDF type and its sons #####################################################

abstract type BRDF end
eval(brdf::BRDF, n::Normal, in_dir::Vec, out_dir::Vec, uv::Vec2D) = return RGB(0.f0, 0.f0, 0.f0)
scatter_ray(brdf::BRDF, pcg::PCG, incoming_dir::Vec, interaction_point::Point, normal::Normal, depth) = return "Abstract method, not implemented"

"""
    DiffuseBRDF(pigment, r)

It creates a Diffuse BRDF, where _r_ is the reflectance of the surface.
If not defined:
- pigment = UniformPigment(RGB(1.,1.,1.));
- r = 1.
"""
struct DiffuseBRDF <: BRDF
    pigment::Pigment
    reflectance::Number

    DiffuseBRDF(pigment::Pigment = UniformPigment(RGB(1.,1.,1.)),
                reflectance::Number = 1. ) =
                new(pigment, reflectance)
end

Base.:≈(brdf1::DiffuseBRDF, brdf2::DiffuseBRDF) = brdf1.pigment ≈ brdf2.pigment && brdf1.reflectance ≈ brdf2.reflectance

eval(brdf::DiffuseBRDF, n::Normal, in_dir::Vec, out_dir::Vec, uv::Vec2D) = return brdf.pigment.get_color(uv) * (brdf.reflectance / π)

function scatter_ray(brdf::DiffuseBRDF, pcg::PCG, incoming_dir::Vec, interaction_point::Point, normal::Normal, depth)
    e1, e2, e3 = create_onb(normal)
    cos_theta_sq = pcg_randf(pcg)
    cosθ, sinθ = sqrt(cos_theta_sq), sqrt(1.0 - cos_theta_sq)
    ϕ = 2.0 * π * pcg_randf(pcg)
    return Ray( interaction_point, e1*cos(ϕ)*cosθ + e2*sin(ϕ)*cosθ + e3*sinθ, 1.0e-3, Inf, depth)
end

"""
    SpecularBRDF(pigment, r)

It creates a Diffuse BRDF, where _r_ is the reflectance of the surface.
If not defined:
- pigment = UniformPigment(RGB(1.f0,1.f0,1.f0));
- r = 1.
"""
struct SpecularBRDF <: BRDF
    pigment::Pigment

    SpecularBRDF(pigment::Pigment = UniformPigment(RGB(1.,1.,1.)) ) = new(pigment)
end

Base.:≈(brdf1::SpecularBRDF, brdf2::SpecularBRDF) = brdf1.pigment ≈ brdf2.pigment

#eval(brdf::SpecularBRDF, n::Normal, in_dir::Vec, out_dir::Vec, uv::Vec2D) = return brdf.pigment.get_color(uv) * (brdf.reflectance / π)
##in realtà non serve main

function scatter_ray(brdf::SpecularBRDF, pcg::PCG, incoming_dir::Vec, interaction_point::Point, normal::Normal, depth)
    ray_dir = normalize(Vec(incoming_dir.x, incoming_dir.y, incoming_dir.z))
    normal = normalize(toVec(normal))
    dot_prod = normal * ray_dir

    return Ray(interaction_point, ray_dir - normal * 2 * dot_prod, 1e-5, Inf, depth)
end

## Code for MATERIAL and its sons #####################################################

function _LightSource(x::RGB)
    if x == RGB()
        return false
    else
        return true
    end
end

"""
    Material(brdf, e_r)

It creates a Material, where brdf is a generic BRDF and e_r is the emitted radiance of the surface (Pigment type).
If not defined:
- brdf = DiffuseBRDF();
- e_r = UniformPigment(RGB(0.f0,0.f0,0.f0)).
"""
struct Material
    brdf::BRDF
    emitted_radiance::Pigment
    emission_intensity::Float32
    is_light_source:: Bool

    Material(brdf::BRDF = DiffuseBRDF(),
             emitted_radiance::Pigment = UniformPigment(RGB(0.f0,0.f0,0.f0)),
             emission_intensity::Float32 = 1.f0,
             is_light_source::Bool = _LightSource(emitted_radiance.color) ) =
             new(brdf, emitted_radiance, emission_intensity, is_light_source)
end

Base.:≈(M1::Material,M2::Material) = M1.brdf ≈ M2.brdf && M1.emitted_radiance ≈ M2.emitted_radiance


text:
function Frame(arg0::Frame, arg1::Transform, arg2::JString)
    return Frame((Frame, Transform, JString), arg0, arg1, arg2)
end

function Frame(arg0::Frame, arg1::Transform, arg2::JString, arg3::jboolean)
    return Frame((Frame, Transform, JString, jboolean), arg0, arg1, arg2, arg3)
end

function Frame(arg0::Frame, arg1::TransformProvider, arg2::JString)
    return Frame((Frame, TransformProvider, JString), arg0, arg1, arg2)
end

function Frame(arg0::Frame, arg1::TransformProvider, arg2::JString, arg3::jboolean)
    return Frame((Frame, TransformProvider, JString, jboolean), arg0, arg1, arg2, arg3)
end


text:
@testset "Smoke Tests" begin
    mts = UnionAll[RNODE, FFJORD, Planar]
    cmts = UnionAll[CondRNODE, CondFFJORD, CondPlanar]
    crs = AbstractResource[CPU1()]
    if has_cuda_gpu()
        push!(crs, CUDALibs())
    end
    tps = DataType[Float64, Float32, Float16]
    opt_apps = ICNF.OptApp[FluxOptApp(), OptimOptApp(), SciMLOptApp()]
    go_oa = SciMLOptApp()
    go_ads = SciMLBase.AbstractADType[Optimization.AutoZygote(), Optimization.AutoForwardDiff()]
    go_mds = Any[ICNF.default_optimizer[FluxOptApp], ICNF.default_optimizer[OptimOptApp]]
    pfm = typeof(ICNF.default_optimizer[OptimOptApp])
    nvars_ = (1:2)
    n_epochs = 2
    batch_size = 8
    n_batch = 2
    n = n_batch*batch_size

    @testset "$(typeof(cr).name.name) | $tp | $nvars Vars | $mt" for
            cr in crs,
            tp in tps,
            nvars in nvars_,
            mt in mts
        data_dist = Beta{tp}(convert(Tuple{tp, tp}, (2, 4))...)
        r = convert(Matrix{tp}, rand(data_dist, nvars, n))

        if mt <: Planar
            nn = PlanarNN(nvars, tanh)
        else
            nn = Chain(
                Dense(nvars, nvars, tanh),
            )
        end
        icnf = mt{tp}(nn, nvars; acceleration=cr)

        @test !isnothing(inference(icnf, TestMode(), r))
        @test !isnothing(inference(icnf, TrainMode(), r))
        @test !isnothing(generate(icnf, TestMode(), n))
        @test !isnothing(generate(icnf, TrainMode(), n))

        @test !isnothing(icnf(r))
        @test !isnothing(loss(icnf, r))
        @test !isnothing(loss_pn(icnf, r))
        @test !isnothing(loss_pln(icnf, r))
        @test !isnothing(loss_f(icnf, FluxOptApp())(r))
        @test !isnothing(loss_f(icnf, OptimOptApp(), [(r,), nothing])(icnf.p))
        @test !isnothing(loss_f(icnf, SciMLOptApp())(icnf.p, SciMLBase.NullParameters(), r))

        @test !isnothing(agg_loglikelihood(icnf, r))

        diff_loss = x -> loss(icnf, r, x)
        @test !isnothing(Zygote.gradient(diff_loss, icnf.p))
        @test !isnothing(Zygote.jacobian(diff_loss, icnf.p))
        @test !isnothing(Zygote.forwarddiff(diff_loss, icnf.p))
        # @test !isnothing(Zygote.diaghessian(diff_loss, icnf.p))
        # @test !isnothing(Zygote.hessian(diff_loss, icnf.p))
        # @test !isnothing(Zygote.hessian_reverse(diff_loss, icnf.p))

        @test !isnothing(ForwardDiff.gradient(diff_loss, icnf.p))
        @test_broken !isnothing(ForwardDiff.jacobian(diff_loss, icnf.p))
        # @test !isnothing(ForwardDiff.hessian(diff_loss, icnf.p))

        d = ICNFDist(icnf)

        @test !isnothing(logpdf(d, r))
        @test !isnothing(pdf(d, r))
        @test !isnothing(rand(d))
        @test !isnothing(rand(d, n))
    end
    @testset "Fitting | $(typeof(cr).name.name) | $tp | $nvars Vars | $mt" for
            cr in crs,
            tp in tps,
            nvars in nvars_,
            mt in mts
        data_dist = Beta{tp}(convert(Tuple{tp, tp}, (2, 4))...)
        r = convert(Matrix{tp}, rand(data_dist, nvars, n))
        df = DataFrame(r', :auto)

        @testset "Using $(typeof(opt_app).name.name)" for
                opt_app in opt_apps
            if mt <: Planar
                nn = PlanarNN(nvars, tanh)
            else
                nn = Chain(
                    Dense(nvars, nvars, tanh),
                )
            end
            icnf = mt{tp}(nn, nvars; acceleration=cr)
            model = ICNFModel(icnf; n_epochs, batch_size, opt_app)
            mach = machine(model, df)
            if !(mt <: Planar) || !(opt_app isa OptimOptApp)
                @test !isnothing(fit!(mach))
                @test !isnothing(MLJBase.transform(mach, df))
            end
        end
        @testset "$(typeof(go_oa).name.name) | Using $(typeof(go_ad).name.name) & $(typeof(go_md).name.name)" for
                go_ad in go_ads,
                go_md in go_mds
            if mt <: Planar
                nn = PlanarNN(nvars, tanh)
            else
                nn = Chain(
                    Dense(nvars, nvars, tanh),
                )
            end
            icnf = mt{tp}(nn, nvars; acceleration=cr)
            model = ICNFModel(icnf; n_epochs, batch_size, opt_app=go_oa, adtype=go_ad, optimizer=go_md)
            mach = machine(model, df)
            if !(mt <: Planar) || !(go_md isa pfm)
                @test !isnothing(fit!(mach))
                @test !isnothing(MLJBase.transform(mach, df))
            end
        end
    end
    @testset "$(typeof(cr).name.name) | $tp | $nvars Vars | $mt" for
            cr in crs,
            tp in tps,
            nvars in nvars_,
            mt in cmts
        data_dist = Beta{tp}(convert(Tuple{tp, tp}, (2, 4))...)
        data_dist2 = Beta{tp}(convert(Tuple{tp, tp}, (4, 2))...)
        r = convert(Matrix{tp}, rand(data_dist, nvars, n))
        r2 = convert(Matrix{tp}, rand(data_dist, nvars, n))

        if mt <: CondPlanar
            nn = PlanarNN(nvars, tanh; cond=true)
        else
            nn = Chain(
                Dense(nvars*2, nvars, tanh),
            )
        end
        icnf = mt{tp}(nn, nvars; acceleration=cr)

        @test !isnothing(inference(icnf, TestMode(), r, r2))
        @test !isnothing(inference(icnf, TrainMode(), r, r2))
        @test !isnothing(generate(icnf, TestMode(), r2, n))
        @test !isnothing(generate(icnf, TrainMode(), r2, n))

        @test !isnothing(icnf(r, r2))
        @test !isnothing(loss(icnf, r, r2))
        @test !isnothing(loss_pn(icnf, r, r2))
        @test !isnothing(loss_pln(icnf, r, r2))
        @test !isnothing(loss_f(icnf, FluxOptApp())(r, r2))
        @test !isnothing(loss_f(icnf, OptimOptApp(), [(r, r2), nothing])(icnf.p))
        @test !isnothing(loss_f(icnf, SciMLOptApp())(icnf.p, SciMLBase.NullParameters(), r, r2))

        @test !isnothing(agg_loglikelihood(icnf, r, r2))

        diff_loss = x -> loss(icnf, r, r2, x)
        @test !isnothing(Zygote.gradient(diff_loss, icnf.p))
        @test !isnothing(Zygote.jacobian(diff_loss, icnf.p))
        @test !isnothing(Zygote.forwarddiff(diff_loss, icnf.p))
        # @test !isnothing(Zygote.diaghessian(diff_loss, icnf.p))
        # @test !isnothing(Zygote.hessian(diff_loss, icnf.p))
        # @test !isnothing(Zygote.hessian_reverse(diff_loss, icnf.p))

        @test !isnothing(ForwardDiff.gradient(diff_loss, icnf.p))
        @test_broken !isnothing(ForwardDiff.jacobian(diff_loss, icnf.p))
        # @test !isnothing(ForwardDiff.hessian(diff_loss, icnf.p))

        d = CondICNFDist(icnf, r2)

        @test !isnothing(logpdf(d, r))
        @test !isnothing(pdf(d, r))
        @test !isnothing(rand(d))
        @test !isnothing(rand(d, n))
    end
    @testset "Fitting | $(typeof(cr).name.name) | $tp | $nvars Vars | $mt" for
            cr in crs,
            tp in tps,
            nvars in nvars_,
            mt in cmts
        data_dist = Beta{tp}(convert(Tuple{tp, tp}, (2, 4))...)
        data_dist2 = Beta{tp}(convert(Tuple{tp, tp}, (4, 2))...)
        r = convert(Matrix{tp}, rand(data_dist, nvars, n))
        r2 = convert(Matrix{tp}, rand(data_dist, nvars, n))
        df = DataFrame(r', :auto)
        df2 = DataFrame(r2', :auto)

        @testset "Using $(typeof(opt_app).name.name)" for
                opt_app in opt_apps
            if mt <: CondPlanar
                nn = PlanarNN(nvars, tanh; cond=true)
            else
                nn = Chain(
                    Dense(nvars*2, nvars, tanh),
                )
            end
            icnf = mt{tp}(nn, nvars; acceleration=cr)
            model = CondICNFModel(icnf; n_epochs, batch_size, opt_app)
            mach = machine(model, (df, df2))
            if !(mt <: CondPlanar) || !(opt_app isa OptimOptApp)
                @test !isnothing(fit!(mach))
                @test !isnothing(MLJBase.transform(mach, (df, df2)))
            end
        end
        @testset "$(typeof(go_oa).name.name) | Using $(typeof(go_ad).name.name) & $(typeof(go_md).name.name)" for
                go_ad in go_ads,
                go_md in go_mds
            if mt <: CondPlanar
                nn = PlanarNN(nvars, tanh; cond=true)
            else
                nn = Chain(
                    Dense(nvars*2, nvars, tanh),
                )
            end
            icnf = mt{tp}(nn, nvars; acceleration=cr)
            model = CondICNFModel(icnf; n_epochs, batch_size, opt_app=go_oa, adtype=go_ad, optimizer=go_md)
            mach = machine(model, (df, df2))
            if !(mt <: CondPlanar) || !(go_md isa pfm)
                @test !isnothing(fit!(mach))
                @test !isnothing(MLJBase.transform(mach, (df, df2)))
            end
        end
    end
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using ParallelAccelerator
using DocOpt
using Images

@acc function harrisCornerDetect(Iin)
  (w, h) = size(Iin)

  Ix  = Array(Float32, w, h)
  Iy  = Array(Float32, w, h)
  Sxx = Array(Float32, w, h)
  Syy = Array(Float32, w, h)
  Sxy = Array(Float32, w, h)

  runStencil(Ix, Iin, :oob_dst_zero) do b, a
    b[0,0] = ((a[-1,-1] * -1.0f0) + (a[-1,0] * -2.0f0) + (a[-1,1] * -1.0f0) + a[1,-1] + (a[1,0] * 2.0f0) + a[1,1]) / 12.0f0
  end
  runStencil(Iy, Iin, :oob_dst_zero) do b, a
    b[0,0] = ((a[-1,-1] * -1.0f0) + (a[0,-1] * -2.0f0) + (a[1,-1] * -1.0f0) + a[-1,1] + (a[0,1] * 2.0f0) + a[1,1]) / 12.0f0
  end

  Ixx = Ix .* Ix
  Iyy = Iy .* Iy
  Ixy = Ix .* Iy

  runStencil(Sxx, Ixx, :oob_dst_zero) do b, a
    b[0,0] = (a[-1,-1] + a[-1,0] + a[-1,1] + a[0,-1] + a[0,0] + a[0,1] + a[1,-1] + a[1,0] + a[1,1])
  end
  runStencil(Syy, Iyy, :oob_dst_zero) do b, a
    b[0,0] = (a[-1,-1] + a[-1,0] + a[-1,1] + a[0,-1] + a[0,0] + a[0,1] + a[1,-1] + a[1,0] + a[1,1])
  end
  runStencil(Sxy, Ixy, :oob_dst_zero) do b, a
    b[0,0] = (a[-1,-1] + a[-1,0] + a[-1,1] + a[0,-1] + a[0,0] + a[0,1] + a[1,-1] + a[1,0] + a[1,1])
  end

  det    = (Sxx .* Syy) .- (Sxy .* Sxy)
  trace  = Sxx .+ Syy
  harris = det .- (0.04f0 .* trace .* trace)

  return harris
end

function main()
    doc = """harris.jl

Harris corner detection.

Usage:
  harris.jl -h | --help
  harris.jl [--img-file=<img-file>]

Options:
  -h --help              Show this screen.
  --img-file=<img-file>  Specify a path to an input (grayscale) image file [default: ../example.jpg].
"""
    arguments = docopt(doc)

    dir = dirname(@__FILE__)
    img_file = joinpath(dir, arguments["--img-file"])

    (fname, ext) = splitext(img_file)
    out_file = string(fname, "-corners", ".jpg")

    println("input file = ", img_file)
    println("output file = ", out_file)

    function harris(input_fname, output_fname)
        local img = convert(Array{Float32, 2}, load(input_fname))
        tic()
        res = harrisCornerDetect(Matrix{Float32}(0, 0))
        println("SELFPRIMED ", toq())

        tic()
        res = harrisCornerDetect(img)
        selftimed = toq()
        res = map((x,y) -> x > 0.001f0 ? 1.0f0 : (y / 2.0f0), res, img)
        save(output_fname, res)
        println("checksum: ", sum(res))
        println("SELFTIMED ", selftimed)
    end

    harris(img_file, out_file)

end

main()


text:
module GilbertCurves

export gilbertindices

"""
    gilbertindices(dims::Tuple{Int,Int}; majdim=dims[1] >= dims[2] ? 1 : 2)

Constructs a vector of `CartesianIndex` objects, orderd by a generalized Hilbert
space-filling, starting at `CartesianIndex(1,1)`.  It will end at
`CartesianIndex(dims[1],1)` if `majdim==1`, or `CartesianIndex(1,dims[2])` if `majdim==2`
(or the closest feasible point).
"""
gilbertindices(dims::Tuple{Int,Int}; kwargs...) =
    gilbertorder(CartesianIndices(dims); kwargs...)

"""
    gilbertorder(mat::AbstractMatrix; majdim=size(mat,1) >= size(mat,2) ? 1 : 2)

Constructs a vector of the elements of `mat`, ordered by a generalized Hilbert
space-filling curve. The list will start at `mat[1,1]`, and end at `mat[end,1]` if
`majdim==1` or `mat[1,end]` if `majdim==2`  (or the closest feasible point).
"""
function gilbertorder(mat::AbstractMatrix{T}; majdim=size(mat,1) >= size(mat,2) ? 1 : 2) where {T}
    list = sizehint!(T[], length(mat))
    if majdim == 1
        append_gilbert!(list, mat)
    else
        append_gilbert!(list, permutedims(mat,(2,1)))
    end
end

function append_gilbert!(list, mat::AbstractMatrix)
    # 1 |*    |
    #   | )   |
    # a |v    |
    a,b = size(mat)
    if a == 1 || b == 1
        # single in one dimension
        append!(list, mat)
    elseif 2a > 3b
        # long case: split into two
        #   +-----+
        # 1 |*    |
        #   ||    |
        # a2|v    |
        #   +-----+
        #   |*    |
        #   ||    |
        # a |v    |
        #   +-----+
        a2 = div(a,2)
        if isodd(a2) && a > 2
            a2 += 1
        end
        append_gilbert!(list, mat[1:a2,:])
        append_gilbert!(list, mat[a2+1:a,:])
    else
        # standard case: split into three
        #      b2
        #   +---+---+
        # 1 |*->|*   |
        #   |   ||   |
        # a2|   ||   |
        #   +---+|   |
        #   |   ||   |
        # a |<-*|v   |
        #   +---+----+
        a2 = div(a,2)
        b2 = div(b,2)
        if isodd(b2) && b > 2
            b2 += 1
        end
        append_gilbert!(list, permutedims(mat[1:a2,1:b2],(2,1)))
        append_gilbert!(list, mat[:,b2+1:b])
        append_gilbert!(list, permutedims(mat[a:-1:a2+1,b2:-1:1],(2,1)))
    end
end

"""
    linearindices(list::Vector{CartesianIndex{2}})

Construct an integer matrix `M` containing the integers `1:length(list)` such that
`M[list[i]] == i`.
"""
function linearindices(list::Vector{CartesianIndex{2}})
    cmax = maximum(list)
    L = zeros(Int,cmax.I)
    for (i,c) in enumerate(list)
        L[c] = i
    end
    return L
end

end


text:
Base.@kwdef struct SurfaceFlow{T,R,L}
    β::T = 0.6                              # constant in Manning's equation
    sl::Vector{T}                           # Slope [m m⁻¹]
    n::Vector{T}                            # Manning's roughness [sl m⁻⅓]
    dl::Vector{T}                           # Drain length [m]
    q::Vector{T} = fill(0.0, length(sl))    # Discharge [m³ s⁻¹]
    q_av::Vector{T} = fill(0.0, length(sl)) # Average discharge [m³ s⁻¹]
    qlat::Vector{T} = fill(0.0, length(sl)) # Lateral discharge [m³ s⁻¹]
    h::Vector{T} = fill(0.0, length(sl))    # Water level [m]
    h_av::Vector{T} = fill(0.0, length(sl)) # Average water level [m]
    Δt::T                                   # Model time step [s]
    its::Int                                # Number of fixed iterations
    width::Vector{T}                        # Flow width [m]
    alpha_term::Vector{T} = pow.(n ./ sqrt.(sl), β)  # Constant part of α
    alpha_pow::T = (2.0 / 3.0) * β          # Used in the power part of α
    α::Vector{T} = alpha_term .* pow.(width .+ 2.0 .* h, alpha_pow) # Constant in momentum equation A = αQᵝ, based on Manning's equation
    eps::T = 1e-03                          # Maximum allowed change in α, if exceeded cross sectional area and h is recalculated
    cel::Vector{T} = fill(0.0, length(sl))  # Celerity of the kinematic wave
    to_river::Vector{T} = fill(0.0, length(sl)) # Part of overland flow [m³ s⁻¹] that flows to the river
    rivercells::Vector{Bool} = fill(false, length(sl)) # Location of river cells (0 or 1)
    wb_pit::Vector{Bool} = fill(false, length(sl)) # Boolean location (0 or 1) of a waterbody (wb, reservoir or lake).
    reservoir_index::Vector{Int} = fill(0, length(sl)) # map cell to 0 (no reservoir) or i (pick reservoir i in reservoir field)
    lake_index::Vector{Int} = fill(0, length(sl))      # map cell to 0 (no lake) or i (pick lake i in lake field)
    reservoir::R = nothing                  # Reservoir model struct of arrays
    lake::L = nothing                       # Lake model struct of arrays

    # TODO unclear why this causes a MethodError
    # function SurfaceFlow{T,R,L}(args...) where {T,R,L}
    #     equal_size_vectors(args)
    #     return new(args...)
    # end
end

statevars(::SurfaceFlow) = (:q,:h,:h_av,)

function update(
    sf::SurfaceFlow,
    network;
    frac_toriver = nothing,
    river = nothing,
    do_iter = false,
    doy = 0,
)
    @unpack graph, order = network

    n = length(order)
    # two options for iteration, fixed or based on courant number.
    if do_iter
        if sf.its > 0
            its = sf.its
        else
            # calculate celerity
            courant = zeros(n)
            for v in order
                if sf.q[v] > 0.0
                    sf.cel[v] = 1.0 / (sf.α[v] * sf.β * pow(sf.q[v], (sf.β - 1.0)))
                    courant[v] = (sf.Δt / sf.dl[v]) * sf.cel[v]
                end
            end
            filter!(x->x≠0.0,courant)
            its = isempty(courant) ? 1 : ceil(Int, (1.25 * quantile!(courant, 0.95)))
        end
    else
        its = 1
    end

    # sub time step
    adt = sf.Δt / its

    q_sum = zeros(n)
    h_sum = zeros(n)
    sf.to_river .= 0.0

    for _ = 1:its
        for v in order
            upstream_nodes = inneighbors(graph, v)
            # for overland flow frac_toriver and river cells need to be defined
            if (frac_toriver !== nothing) && (river !== nothing)
                # for a river cell without a reservoir or lake (wb_pit is false) part of the upstream surface flow
                # goes to the river (frac_toriver) and part goes to the surface flow reservoir (1.0 - frac_toriver)
                # upstream nodes with a reservoir or lake are excluded
                if river[v] && !sf.wb_pit[v] && sf.width[v] > 0.0
                    qin = sum(
                        sf.q[i] * (1.0 - frac_toriver[i])
                        for i in upstream_nodes if !sf.wb_pit[i]
                    )
                    # TODO: if in a loop should also take average!!!
                    sf.to_river[v] += sum(
                        sf.q[i] * frac_toriver[i] for i in upstream_nodes if !sf.wb_pit[i]
                    )
                    # for a river cell with a reservoir or lake (wb_pit is true) all upstream surface flow goes
                    # to the river.
                elseif river[v] && sf.wb_pit[v] && sf.width[v] == 0.0
                    sf.to_river[v]  += sum_at(sf.q, upstream_nodes)
                    qin = 0.0
                else
                    qin = sum_at(sf.q, upstream_nodes)
                end
                # for all the other cells all upstream surface flow goes to the surface flow reservoir.
            else
                qin = sum_at(sf.q, upstream_nodes)
            end
            if !isnothing(sf.reservoir) && sf.reservoir_index[v] != 0
                # run reservoir model and copy reservoir outflow to river cell
                i = sf.reservoir_index[v]
                update(sf.reservoir, i, qin, adt)
                sf.q[v] = sf.reservoir.outflow[i]

            elseif !isnothing(sf.lake) && sf.lake_index[v] != 0
                # run lake model and copy lake outflow to river cell
                i = sf.lake_index[v]
                update(sf.lake, i, qin, doy)
                sf.q[v] = sf.lake.outflow[i]
            else
                sf.q[v] =
                    kinematic_wave(qin, sf.q[v], sf.qlat[v], sf.α[v], sf.β, adt, sf.dl[v])
            end

            # update alpha
            crossarea = sf.α[v] * pow(sf.q[v], sf.β)
            sf.h[v] = crossarea / sf.width[v]
            wetper = sf.width[v] + (2.0 * sf.h[v]) # wetted perimeter
            α = sf.α[v]
            sf.α[v] = sf.alpha_term[v] * pow(wetper, sf.alpha_pow)

            if abs(α - sf.α[v]) > sf.eps
                crossarea = sf.α[v] * pow(sf.q[v], sf.β)
                sf.h[v] = crossarea / sf.width[v]
            end

            q_sum[v] += sf.q[v]
            h_sum[v] += sf.h[v]

        end

    end
    sf.q_av .= q_sum ./ its
    sf.h_av .= h_sum ./ its
    sf.to_river .= sf.to_river ./ its

end

Base.@kwdef struct LateralSSF{T}
    kh₀::Vector{T}                          # Horizontal hydraulic conductivity at soil surface [mm Δt⁻¹]
    f::Vector{T}                            # A scaling parameter [mm⁻¹] (controls exponential decline of kh₀)
    soilthickness::Vector{T}                # Soil thickness [mm]
    θₑ::Vector{T}                           # Effective porosity [-]
    Δt::T = 1.0                               # model time step
    βₗ::Vector{T}                           # Slope [m m⁻¹]
    dl::Vector{T}                           # Drain length [mm]
    dw::Vector{T}                           # Flow width [mm]
    zi::Vector{T} = fill(mv, length(f))     # Pseudo-water table depth [mm] (top of the saturated zone)
    exfiltwater::Vector{T} = fill(mv, length(f))  # Exfiltration [mm]  (groundwater above surface level, saturated excess conditions)
    recharge::Vector{T} = fill(mv, length(f))     # Net recharge to saturated store [mm]
    ssf::Vector{T} # Subsurface flow [mm³ Δt⁻¹]
    ssfin::Vector{T} = fill(mv, length(f))
    ssfmax::Vector{T} = ((kh₀ .* βₗ) ./ f) .* (1.0 .- exp.(-f .* soilthickness))     # Maximum subsurface flow [mm² Δt⁻¹]
    to_river::Vector{T} = zeros(length(f))  # Part of subsurface flow [mm³ Δt⁻¹] that flows to the river
    wb_pit::Vector{Bool} = zeros(Bool, length(f)) # Boolean location (0 or 1) of a waterbody (wb, reservoir or lake).

    function LateralSSF{T}(args...) where {T}
        equal_size_vectors(args)
        return new(args...)
    end
end

# depends on ini file settings (optional: glaciers, snow, irrigation)
statevars(::LateralSSF) = (:ssf,)

function update(ssf::LateralSSF, network, frac_toriver, river)
    @unpack graph, order = network
    for v in order
        upstream_nodes = inneighbors(graph, v)
        # for a river cell without a reservoir or lake (wb_pit is false) part of the upstream subsurface flow
        # goes to the river (frac_toriver) and part goes to the subsurface flow reservoir (1.0 - frac_toriver)
        # upstream nodes with a reservoir or lake are excluded
        if river[v] && !ssf.wb_pit[v]
            ssf.ssfin[v] = sum(
                ssf.ssf[i] * (1.0 - frac_toriver[i])
                for i in upstream_nodes if !ssf.wb_pit[i]
            )
            ssf.to_river[v] =
                sum(ssf.ssf[i] * frac_toriver[i] for i in upstream_nodes if !ssf.wb_pit[i])
            # for a river cell with a reservoir or lake (wb_pit is true) all upstream subsurface flow goes
            # to the river.
        elseif river[v] && ssf.wb_pit[v]
            ssf.to_river[v] = sum_at(ssf.ssf, upstream_nodes)
            ssf.ssfin[v] = 0.0
            # for all the other cells all upstream subsurface flow goes to the subsurface flow reservoir.
        else
            ssf.ssfin[v] = sum_at(ssf.ssf, upstream_nodes)
        end
        ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf(
            ssf.ssfin[v],
            ssf.ssf[v],
            ssf.zi[v],
            ssf.recharge[v],
            ssf.kh₀[v],
            ssf.βₗ[v],
            ssf.θₑ[v],
            ssf.f[v],
            ssf.soilthickness[v],
            ssf.Δt,
            ssf.dl[v],
            ssf.dw[v],
            ssf.ssfmax[v],
        )
    end
end


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

## Broadcast for Taylor1 and TaylorN

import .Broadcast: BroadcastStyle, Broadcasted, broadcasted

# BroadcastStyle definitions and basic precedence rules
struct Taylor1Style{T} <: Base.Broadcast.AbstractArrayStyle{0} end
Taylor1Style{T}(::Val{N}) where {T, N} = Base.Broadcast.DefaultArrayStyle{N}()
BroadcastStyle(::Type{<:Taylor1{T}}) where {T} = Taylor1Style{T}()
BroadcastStyle(::Taylor1Style{T}, ::Base.Broadcast.DefaultArrayStyle{0}) where {T} = Taylor1Style{T}()
BroadcastStyle(::Taylor1Style{T}, ::Base.Broadcast.DefaultArrayStyle{1}) where {T} = Base.Broadcast.DefaultArrayStyle{1}()
#
struct HomogeneousPolynomialStyle{T} <: Base.Broadcast.AbstractArrayStyle{0} end
HomogeneousPolynomialStyle{T}(::Val{N}) where {T, N}= Base.Broadcast.DefaultArrayStyle{N}()
BroadcastStyle(::Type{<:HomogeneousPolynomial{T}}) where {T} = HomogeneousPolynomialStyle{T}()
BroadcastStyle(::HomogeneousPolynomialStyle{T}, ::Base.Broadcast.DefaultArrayStyle{0}) where {T} = HomogeneousPolynomialStyle{T}()
BroadcastStyle(::HomogeneousPolynomialStyle{T}, ::Base.Broadcast.DefaultArrayStyle{1}) where {T} = Base.Broadcast.DefaultArrayStyle{1}()
#
struct TaylorNStyle{T} <: Base.Broadcast.AbstractArrayStyle{0} end
TaylorNStyle{T}(::Val{N}) where {T, N}= Base.Broadcast.DefaultArrayStyle{N}()
BroadcastStyle(::Type{<:TaylorN{T}}) where {T} = TaylorNStyle{T}()
BroadcastStyle(::TaylorNStyle{T}, ::Base.Broadcast.DefaultArrayStyle{0}) where {T} = TaylorNStyle{T}()
BroadcastStyle(::TaylorNStyle{T}, ::Base.Broadcast.DefaultArrayStyle{1}) where {T} = Base.Broadcast.DefaultArrayStyle{1}()

# Precedence rules for mixtures
BroadcastStyle(::TaylorNStyle{Taylor1{T}}, ::Taylor1Style{T}) where {T} = TaylorNStyle{Taylor1{T}}()
BroadcastStyle(::Taylor1Style{TaylorN{T}}, ::TaylorNStyle{T}) where {T} = Taylor1Style{TaylorN{T}}()

# Extend eltypes so things like [1.0] .+ t work
Base.Broadcast.eltypes(t::Tuple{Taylor1,AbstractArray}) =
    Tuple{Base.Broadcast._broadcast_getindex_eltype([t[1]]), Base.Broadcast._broadcast_getindex_eltype(t[2])}
Base.Broadcast.eltypes(t::Tuple{AbstractArray,Taylor1}) =
    Tuple{Base.Broadcast._broadcast_getindex_eltype(t[1]), Base.Broadcast._broadcast_getindex_eltype([t[2]])}
Base.Broadcast.eltypes(t::Tuple{HomogeneousPolynomial,AbstractArray}) =
    Tuple{Base.Broadcast._broadcast_getindex_eltype([t[1]]), Base.Broadcast._broadcast_getindex_eltype(t[2])}
Base.Broadcast.eltypes(t::Tuple{AbstractArray,HomogeneousPolynomial}) =
    Tuple{Base.Broadcast._broadcast_getindex_eltype(t[1]), Base.Broadcast._broadcast_getindex_eltype([t[2]])}
Base.Broadcast.eltypes(t::Tuple{TaylorN,AbstractArray}) =
    Tuple{Base.Broadcast._broadcast_getindex_eltype([t[1]]), Base.Broadcast._broadcast_getindex_eltype(t[2])}
Base.Broadcast.eltypes(t::Tuple{AbstractArray,TaylorN}) =
    Tuple{Base.Broadcast._broadcast_getindex_eltype(t[1]), Base.Broadcast._broadcast_getindex_eltype([t[2]])}

# # We follow https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration-1
# "`A = find_taylor(As)` returns the first Taylor1 among the arguments."
# find_taylor(bc::Broadcasted) = find_taylor(bc.args)
# find_taylor(args::Tuple) = find_taylor(find_taylor(args[1]), Base.tail(args))
# find_taylor(x) = x
# find_taylor(a::Taylor1, rest) = a
# find_taylor(a::HomogeneousPolynomial, rest) = a
# find_taylor(a::TaylorN, rest) = a
# find_taylor(::AbstractArray, rest) = find_taylor(rest)
#
# # Extend similar
# function similar(bc::Broadcasted{Taylor1Style{S}}, ::Type{T}) where {S, T}
#     # Proper promotion
#     R = Base.Broadcast.combine_eltypes(bc.f, bc.args)
#     # Scan the inputs for the Taylor1:
#     A = find_taylor(bc)
#     # Create the output
#     return Taylor1(similar(A.coeffs, R), A.order)
# end
#
# function similar(bc::Broadcasted{HomogeneousPolynomialStyle{S}}, ::Type{T}) where {S, T}
#     # Proper promotion
#     # combine_eltypes(f, args::Tuple) = Base._return_type(f, eltypes(args))
#     R = Base.Broadcast.combine_eltypes(bc.f, bc.args)
#     # Scan the inputs for the HomogeneousPolynomial:
#     A = find_taylor(bc)
#     # Create the output
#     return HomogeneousPolynomial(similar(A.coeffs, R), A.order)
# end
#
# function similar(bc::Broadcasted{TaylorNStyle{S}}, ::Type{T}) where {S, T}
#     # Proper promotion
#     R = Base.Broadcast.combine_eltypes(bc.f, bc.args)
#     # Scan the inputs for the TaylorN:
#     A = find_taylor(bc)
#     # Create the output
#     return TaylorN(similar(A.coeffs, R), A.order)
# end

# Adapted from Base.Broadcast.copyto!, base/broadcasting.jl, line 832
for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
    @eval begin
        @inline function copyto!(dest::$T{T}, bc::Broadcasted) where {T}
            axes(dest) == axes(bc) || Base.Broadcast.throwdm(axes(dest), axes(bc))
            # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
            if bc.f === identity && bc.args isa Tuple{$T{T}} # only a single input argument to broadcast!
                A = bc.args[1]
                if axes(dest) == axes(A)
                    return copyto!(dest, A)
                end
            end
            bc′ = Base.Broadcast.preprocess(dest, bc)
            copyto!(dest, bc′[1])
            return dest
        end
    end
end

# Broadcasted extensions
@inline broadcasted(::Taylor1Style{T}, ::Type{Float32}, a::Taylor1{T}) where {T} =
    Taylor1(Float32.(a.coeffs), a.order)
@inline broadcasted(::TaylorNStyle{T}, ::Type{Float32}, a::TaylorN{T}) where {T} =
    convert(TaylorN{Float32}, a)

# # This prevents broadcasting being applied to the Taylor1/TaylorN params
# # for the mutating functions, and to act only in `k`
# for (T, TS) in ((:Taylor1, :Taylor1Style), (:TaylorN, :TaylorNStyle))
#     for f in (add!, subst!, sqr!, sqrt!, exp!, log!, identity!, zero!,
#         one!, abs!, abs2!, deg2rad!, rad2deg!)
#         @eval begin
#             @inline function broadcasted(::$TS{T}, fn::typeof($f), r::$T{T}, a::$T{T}, k) where {T}
#                 @inbounds for i in eachindex(k)
#                     fn(r, a, k[i])
#                 end
#                 nothing
#             end
#         end
#     end
#     for f in (sincos!, tan!, asin!, acos!, atan!, sinhcosh!, tanh!)
#         @eval begin
#             @inline function broadcasted(::$TS{T}, fn::typeof($f), r::$T{T}, a::$T{T}, b::Taylor1{T}, k) where {T}
#                 @inbounds for i in eachindex(k)
#                     fn(r, a, b, k[i])
#                 end
#                 nothing
#             end
#         end
#     end
# end


text:
import Random.shuffle
import Random.shuffle!

function swap_degree!(Y::BinaryNetwork)
   iy = interactions(Y)
   i1, i2 = sample(iy, 2, replace=false)
   n1, n2 = (from=i1.from, to=i2.to), (from=i2.from, to=i1.to)

   while (n2 ∈ iy)|(n1 ∈ iy)
      i1, i2 = sample(iy, 2, replace=false)
      n1, n2 = (from=i1.from, to=i2.to), (from=i2.from, to=i1.to)
   end

   for old_i in [i1, i2, n1, n2]
      i = something(findfirst(isequal(old_i.from), species(Y; dims=1)), 0)
      j = something(findfirst(isequal(old_i.to), species(Y; dims=2)), 0)
      Y.A[i,j] = !Y.A[i,j]
   end
end

function swap_fill!(Y::BinaryNetwork)
   iy = interactions(Y)
   i1 = sample(iy)
   n1 = (from=sample(species(Y; dims=1)), to=sample(species(Y; dims=2)))

   while (n1 ∈ iy)
      n1 = (from=sample(species(Y; dims=1)), to=sample(species(Y; dims=2)))
   end

   for old_i in [i1, n1]
      i = something(findfirst(isequal(old_i.from), species(Y; dims=1)), 0)
      j = something(findfirst(isequal(old_i.to), species(Y; dims=2)), 0)
      Y.A[i,j] = !Y.A[i,j]
   end
end

function swap_vulnerability!(Y::BinaryNetwork)
   iy = interactions(Y)
   i1 = sample(iy)
   n1 = (from=sample(species(Y; dims=1)), to=i1.to)

   while (n1 ∈ iy)
      n1 = (from=sample(species(Y; dims=1)), to=i1.to)
   end

   for old_i in [i1, n1]
      i = something(findfirst(isequal(old_i.from), species(Y; dims=1)), 0)
      j = something(findfirst(isequal(old_i.to), species(Y; dims=2)), 0)
      Y.A[i,j] = !Y.A[i,j]
   end
end

function swap_generality!(Y::BinaryNetwork)
   iy = interactions(Y)
   i1 = sample(iy)
   n1 = (from=i1.from, to=sample(species(Y; dims=2)))

   while (n1 ∈ iy)
      n1 = (from=i1.from, to=sample(species(Y; dims=2)))
   end

   for old_i in [i1, n1]
      i = something(findfirst(isequal(old_i.from), species(Y; dims=1)), 0)
      j = something(findfirst(isequal(old_i.to), species(Y; dims=2)), 0)
      Y.A[i,j] = !Y.A[i,j]
   end
end

"""
    shuffle(N::BinaryNetwork; constraint::Symbol=:degree)

Return a shuffled copy of the network (the original network is not modified).
See `shuffle!` for a documentation of the keyword arguments.

#### References

Fortuna, M.A., Stouffer, D.B., Olesen, J.M., Jordano, P., Mouillot, D., Krasnov,
B.R., Poulin, R., Bascompte, J., 2010. Nestedness versus modularity in
ecological networks: two sides of the same coin? Journal of Animal Ecology 78,
811–817. https://doi.org/10.1111/j.1365-2656.2010.01688.x
"""
function shuffle(N::BinaryNetwork; constraint::Symbol=:degree)
   Y = copy(N)
   shuffle!(Y; constraint=constraint)
   return Y
end

"""
    shuffle!(N::BinaryNetwork; constraint::Symbol=:degree)

Shuffles interactions inside a network (the network is *modified*), under the
following `constraint`:

- :degree, which keeps the degree distribution intact
- :generality, which keeps the out-degree distribution intact
- :vulnerability, which keeps the in-degree distribution intact
- :fill, which moves interactions around freely

The function will take two interactions, and swap the species establishing them.
By repeating the process a large enough number of times, the resulting network
should be relatively random. Note that this function will conserve the degree
(when appropriate under the selected constraint) of *every* species. Calling the
function will perform **a single** shuffle. If you want to repeat the shuffling
a large enough number of times, you can use something like:

    [shuffle!(n) for i in 1:10_000]

If the keyword arguments are invalid, the function will throw an
`ArgumentError`.

#### References

Fortuna, M.A., Stouffer, D.B., Olesen, J.M., Jordano, P., Mouillot, D., Krasnov,
B.R., Poulin, R., Bascompte, J., 2010. Nestedness versus modularity in
ecological networks: two sides of the same coin? Journal of Animal Ecology 78,
811–817. https://doi.org/10.1111/j.1365-2656.2010.01688.x
"""
function shuffle!(N::BinaryNetwork; constraint::Symbol=:degree)
   constraint ∈ [:degree, :generality, :vulnerability, :fill] || throw(ArgumentError("The constraint argument you specificied ($(constraint)) is invalid -- see ?shuffle! for a list."))

   f = EcologicalNetworks.swap_degree!
   constraint == :generality && (f = EcologicalNetworks.swap_generality!)
   constraint == :vulnerability && (f = EcologicalNetworks.swap_vulnerability!)
   constraint == :fill && (f = EcologicalNetworks.swap_fill!)

   f(N)
end


text:
Base.:+(a::T, b::V) where {T<:Rule,V<:Rule} = Ruleset(a, b)
Base.:+(a::T, b::V) where {T<:Ruleset,V<:Rule} = Ruleset(rules(a)..., b)
Base.:+(a::T, b::V) where {T<:Rule,V<:Ruleset} = Ruleset(a, rules(b)...)
Base.:+(a::T, b::V) where {T<:Ruleset,V<:Ruleset} = Ruleset(rules(a)..., rules(b)...)


text:
module A
    export foo
    foo() = true
    Foo() = true
    macro foo(x) x end
    macro Foo() true end
    module B
        export bar
        bar() = true
        Bar() = true
        macro bar(x, y) x, y end
        module C
            export baz
            baz() = true
            Baz() = true
            macro baz(x, y, z) x, y, z end
        end
    end
end


text:
include("test_1.jl")
include("test_2.jl")
include("test_3.jl")
include("test_4.jl")
include("test_7.jl")


text:
module TestEnsembleMethods

using Test
using Random
using TSML

function getprediction(model::TSLearner,data::Dict)
  Random.seed!(126)
  trfeatures = data[:trfeatures]
  tstfeatures = data[:tstfeatures]
  troutput = data[:troutput]
  tstoutput = data[:tstoutput]
  fit!(model,trfeatures,troutput)
  trresults = transform!(model,tstfeatures)
  sum(trresults .== tstoutput)/length(tstoutput)
end

function test_ensembles()
  modeldirname = joinpath(dirname(pathof(TSML)),"../data/realdatatsclassification/model")
  trdirname = joinpath(dirname(pathof(TSML)),"../data/realdatatsclassification/training")
  tstdirname = joinpath(dirname(pathof(TSML)),"../data/realdatatsclassification/testing")

  frange = 5:20
  trdata = getstats(trdirname)
  trfeatures = trdata[:,frange]
  troutput = trdata[:,:dtype]

  tstdata = getstats(tstdirname)
  tstfeatures = tstdata[:,frange]
  tstoutput = tstdata[:,:dtype]

  models = [VoteEnsemble(),StackEnsemble(),BestLearner()]

  data = Dict(:trfeatures => trfeatures,
              :tstfeatures => tstfeatures,
              :troutput => troutput,
              :tstoutput => tstoutput
             )
  #acc = [4//6,3//6,5//6]
  ndx=1
  for model in models
    @test getprediction(model,data) > 0.10
    ndx += 1
  end
end
@testset "Ensemble learners" begin
  test_ensembles()
end

end # module


text:
type Image
    ptr::Ptr{Void}

    function Image(ptr::Ptr{Void})
        i = new(ptr)
        finalizer(i, destroy)
        i
    end
end

function Image(filename::AbstractString)
    Image(ccall((:sfImage_createFromFile, libcsfml_graphics), Ptr{Void}, (Ptr{Cchar},), filename))
end

function Image(width::Integer, height::Integer, color::Color = SFML.black)
    Image(ccall((:sfImage_createFromColor, libcsfml_graphics), Ptr{Void}, (UInt32, UInt32, Color,), width, height, color))
end

function copy(image::Image)
    Image(ccall((:sfImage_copy, libcsfml_graphics), Ptr{Void}, (Ptr{Void},), image.ptr))
end

function destroy(image::Image)
    ccall((:sfImage_destroy, libcsfml_graphics), Void, (Ptr{Void},), image.ptr)
end

function save_to_file(image::Image, filename::AbstractString)
    ccall((:sfImage_saveToFile, libcsfml_graphics), Bool, (Ptr{Void}, Ptr{Cchar},), image.ptr, filename)
end

function set_pixel(image::Image, x::Integer, y::Integer, color::Color)
    ccall((:sfImage_setPixel, libcsfml_graphics), Void, (Ptr{Void}, UInt32, UInt32, Color,), image.ptr, x, y, color)
end

function get_pixel(image::Image, x::Integer, y::Integer)
    ccall((:sfImage_getPixel, libcsfml_graphics), Color, (Ptr{Void}, UInt32, UInt32,), image.ptr, x, y)
end

function get_pixels(image::Image)
    imgsize = get_size(image)
    ptr = ccall((:sfImage_getPixelsPtr, libcsfml_graphics), Ptr{UInt8}, (Ptr{Void},), image.ptr)
    # @compat unsafe_wrap(Vector{UInt8}, ptr, imgsize.x * imgsize.y)
    convert(Array{UInt8}, unsafe_string(ptr, imgsize.x * imgsize.y))
end

function get_size(image::Image)
    ccall((:sfImage_getSize, libcsfml_graphics), Vector2u, (Ptr{Void},), image.ptr)
end

function flip_horizontally(image::Image)
    ccall((:sfImage_flipHorizontally, libcsfml_graphics), Void, (Ptr{Void},), image.ptr)
end

function flip_vertically(image::Image)
    ccall((:sfImage_flipVertically, libcsfml_graphics), Void, (Ptr{Void},), image.ptr)
end


text:
n = 20
T = 1000
K = 100

#x = zeros(K, T, n)
x = zeros(n, T, K)
y = randn(n)

function fill(x, y)
   #K = size(x, 1)
   #T = size(x, 2)
   T = size(x, 2)
   K = size(x, 3)

   for k=1:K
      for t=1:T
         x[:,t,k] = y
         #x[k,t,:] = y
      end
   end
end

@time fill(x,y)
Profile.clear_malloc_data()
@time fill(x,y)


text:
# copied from Clang.jl
# normalize literal with a size suffix
function literally(tok)
    # note: put multi-character first, or it will break out too soon for those!
    literalsuffixes = ["ULL", "Ull", "uLL", "ull", "LLU", "LLu", "llU", "llu",
                       "LL", "ll", "UL", "Ul", "uL", "ul", "LU", "Lu", "lU", "lu",
                       "U", "u", "L", "l", "F", "f"]

    function literal_totype(literal, txt)
      literal = lowercase(literal)

      # Floats following http://en.cppreference.com/w/cpp/language/floating_literal
      float64 = occursin(".", txt) && occursin("l", literal)
      float32 = occursin("f", literal)

      if float64 || float32
        float64 && return "Float64"
        float32 && return "Float32"
      end

      # Integers following http://en.cppreference.com/w/cpp/language/integer_literal
      unsigned = occursin("u", literal)
      nbits = count(x -> x == 'l', literal) == 2 ? 64 : 32
      return "$(unsigned ? "U" : "")Int$nbits"
    end

    token_kind = kind(tok)
    txt = tok.text |> strip
    if token_kind == CXToken_Identifier || token_kind == CXToken_Punctuation
        # pass
    elseif token_kind == CXToken_Literal
        for sfx in literalsuffixes
            if endswith(txt, sfx)
                type = literal_totype(sfx, txt)
                txt = txt[1:end-length(sfx)]
                txt = "$(type)($txt)"
                break
            end
        end
    end
    return txt
end

function translate(cursor::CLIntegerLiteral)
    tok = tokenize(cursor)[1]
    if tok.text != "0" &&
       startswith(tok.text, "0") &&
       !startswith(lowercase(tok.text), "0x")
        ex = Meta.parse("0o"*literally(tok))
        return MetaExpr(ex, cursor)
    else
        ex = Meta.parse(literally(tok))
        return MetaExpr(ex, cursor)
    end
end

function translate(cursor::CLFloatingLiteral)
    ex = tokenize(cursor)[1] |> literally |> Meta.parse
    return MetaExpr(ex, cursor)
end

function translate(cursor::CLCharacterLiteral)
    ex = tokenize(cursor)[1].text |> Meta.parse
    return MetaExpr(ex, cursor)
end

function translate(cursor::CLStringLiteral)
    ex = tokenize(cursor)[1].text |> Meta.parse
    return MetaExpr(ex, cursor)
end


text:
module TestDash

using ComoniconTypes
using ComoniconTargetExpr
using ComoniconTargetExpr: emit_expr, emit_body, emit_norm_body, emit_dash_body
using Test

const test_args = Ref{Vector{Any}}()
const test_kwargs = Ref{Vector{Any}}()

function foo(a, b=2, c...; kwargs...)
    test_args[] = [a, b, c...]
    test_kwargs[] = [kwargs...]
end

cmd = Entry(;
    version=v"1.1.0",
    root=LeafCommand(;
        fn=foo,
        name="leaf",
        args=[
            Argument(;name="a", type=Int),
            Argument(;name="b", type=Int, require=false),
        ],
        vararg=Argument(;name="c", type=Int, vararg=true),
        options=Dict(
            "option-a" => Option(;sym=:option_a, hint="int", type=Int, short=true),
            "option-b" => Option(;sym=:option_b, hint="float64", type=Float64),
        ),
        flags=Dict(
            "flag-a" => Flag(;sym=:flag_a, short=true),
            "flag-b" => Flag(;sym=:flag_b),
        ),
    )
)

eval(emit_expr(cmd))

@testset "test leaf optional argument" begin
    @test command_main(["3", "2", "5", "6", "7", "--option-a=2", "--option-b", "2.3"]) == 0
    @test test_args[] == Any[3, 2, 5, 6, 7]
    @test test_kwargs[] == [:option_a=>2, :option_b=>2.3]

    @test command_main(["--option-a=2", "--option-b=2.3", "--", "3", "2"]) == 0
    @test test_args[] == Any[3, 2]
    @test test_kwargs[] == [:option_a=>2, :option_b=>2.3]

    @test command_main(["3", "2"]) == 0
    @test test_args[] == Any[3, 2]

    @test command_main(["cmd", "--", "3", "2"]) == 1
end

end


text:
module Libtask

export CTask, consume, produce, TArray, get, tzeros, tfill, TRef

# Try to load the binary dependency
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("Libtask is not properly installed. Please run `import Pkg; Pkg.build(\"Libtask\")`")
end
include(depsjl_path)

function __init__()
    check_deps()
end

@static if VERSION < v"1.0.9999" # (-, v1.1)
    const libtask = libtask_v1_0
elseif VERSION < v"1.1.9999" # [v1.1, v1.2)
    const libtask = libtask_v1_1
elseif VERSION < v"1.2.9999" # [v1.2, v1.3)
    const libtask = libtask_v1_2
else # [v1.3, +)
    const libtask = libtask_v1_3
end

include("ctask.jl")
include("tarray.jl")

end


text:
#############################################################################
#Copyright (c) 2016-2017 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# This files contains functions to remove the signal of a cristal from the glass Raman signal.
#
#
#############################################################################

"""
	ctxremoval(liste,in_path,out_path,roi_all;input_properties=('\t',0),algorithm="FastICA",plot_intermediate_show = "no",plot_mixing_show = "yes",plot_final_show = "no",save_fig_switch = "yes", shutdown = 1300.,scaling=100.)

INPUTS

	liste: Array{Float64}, an array contaning the information for getting the spectra.

		Column 1: the name and relative path for the crystal spectra;

		Column 2: the name and relative path for the (mixed) glass spectra;

		Column 3: smo_water, the coefficient of smothing for the spline that fits the backgrous below the signal of water (if any) in the glass spectra

		Column 4: number of iteration for the generation of new mixed spectra

		Column 5: the K_Start parameter, set at 0.0

		Column 6: the K_increament parameter, for mixing the signals

	in_path: String, the relative location of the data, e.g. "./raw/"

	out_path: String, the relative location where you want to save the corrected spectra, e.g. "./treated/"

	roi_all: Tuple, contains 2 arrays and 2 Float64 numbers. The 2 arrays indicate the regions of interest where the background correction is applied, for the cristal and the glass. The  2 float numbers indicate the starting and ending frequency of the peak used to correct the spectra from any shift in frequency. For instance:

		Those are the roi for fitting the baseline on the crystal (roi_ctx) and glass (roi_glass) signals:

		roi_ctx = [1260. 2000.;2000. 4000.]

		roi_glass = [1260. 2000.;2000. 3000.;3750. 4000.]

		We have a strong peak from the crystal at ~650 cm-1 that we can use to correct the spectra from any shift in frequency. So we indicate the values here:

		roi_xshift_low = 655.

		roi_xshift_high = 670.

		Then we construct the final tuple as:

		roi_all = (roi_ctx,roi_glass,roi_xshift_low,roi_xshift_high)

OPTIONS

	input_properties: Tuple, this tuple contains the delimiter and the number of lines to skip in the raw data files. Default = ('\t',0);

	algorithm: String, This indicates if the FastICA or the Non-negative Matrix Factorisation (NMF) algorithms from SciKit Learn will be used. Default = "FastICA";

	plot_intermediate_show: String, This should be equal to "yes" or "no". It displays the intermediate figures. Default = "no";

	plot_mixing_show: String, This should be equal to "yes" or "no". It displays the figures showing the mixing step. Default = "yes";

	plot_final_show: String, This should be equal to "yes" or "no". It displays the final figures, showing the background subtraction and the retrieved signals. Default = "yes";

	save_fig_switch: String, This should be equal to "yes" or "no". It indicates if you want to save the final figures in the location indicated by out_path;

	shutdown: Float64, indicates where you consider the signals from silicate units to stop. Default = 1300.0;

	scaling: Float64, the retrieved spectra are scaled to the original spectra using the Boson peak, located ~ 60-80 cm-1. This parameters indicates where you consider the Boson peak to stop for the scaling procedure. No need to put a too high value, as you might get strong crystal signals at frequencies > 100-150 cm-1.

OUTPUTS

	All the corrected spectra and figures are saved in the location indicated in out_path. No direct outputs in Julia.
"""
function ctxremoval(liste,in_path,out_path,roi_all;input_properties=('\t',0),algorithm="FastICA",plot_intermediate_show = "no",plot_mixing_show = "yes",plot_final_show = "no",save_fig_switch = "yes", shutdown = 1300.,scaling=100.)

	# Model for final adjustement
	model_ctx_adj(x, p) = p[1]*x

	# GRABING THE REGIONS OF INTEREST
	roi_ctx = roi_all[1]
	roi_glass = roi_all[2]
	roi_xshift_low = roi_all[3]
	roi_xshift_high = roi_all[4]

	# NUMBER OF FILES?
	nb_exp = size(liste,1)

	# SARTING THE LOOP TO TREAT EACH COUPLE OF SPECTRA
	for i = 1:nb_exp

	    # IMPORTING DATA
	    ctx = readdlm(in_path*liste[i,1],input_properties[1],skipstart=input_properties[2])
	    glass = readdlm(in_path*liste[i,2],input_properties[1],skipstart=input_properties[2])

	    #CHECKING THAT X IS INCREASING
	    if glass[end,1] < glass[1,1]
	        glass = flipdim(glass,1)
	    end

	    if ctx[end,1] < ctx[1,1]
	        ctx = flipdim(ctx,1)
	    end

	    # THE COMMON X AXIS, TAKEN AS THE GLASS ONE
	    true_x = collect(minimum(glass[:,1]):0.8: maximum(glass[:,1]))

	    # RESAMPLING THE TWO SPECTRA WITH THE SAME X
	    data = zeros(size(true_x,1),2)

	    spl_ctx = Spline1D(ctx[:,1],ctx[:,2])
	    spl_gls = Spline1D(glass[:,1],glass[:,2])

	    data[:,1] = evaluate(spl_ctx, true_x)
	    data[:,2] = evaluate(spl_gls, true_x)

	    # A VARIABLE TO STORE THE BASELINE CORRECTED SPECTRA
	    bas_matrix = zeros(size(true_x,1),2)

	    # WE GET THE SMOOTHING SPLINE COEFFICIENT
	    smo_water = liste[i,3]

	    # BASELINE CTX
	    ycorr_ctx_1, baseline1_ctx = baseline(true_x[true_x[:].<roi_ctx[1,2]], data[true_x[:].<roi_ctx[1,2],1],reshape(roi_ctx[1,:],1,2),"poly",p=1.0)
	    ycorr_ctx_2, baseline2_ctx = baseline(true_x[true_x[:].>=roi_ctx[1,2]],data[true_x[:].>=roi_ctx[1,2],1],reshape(roi_ctx[2,:],1,2),"gcvspline",p=smo_water)
	    baseline_ctx = [baseline1_ctx;baseline2_ctx]
	    bas_matrix[:,1] = [ycorr_ctx_1;ycorr_ctx_2]

	    # BASELINE GLASS
	    ycorr_gls_1, baseline1_gls = baseline(true_x[true_x[:].<roi_glass[1,2]], data[true_x[:].<roi_glass[1,2],2],reshape(roi_glass[1,:],1,2),"poly",p=1.0)
	    ycorr_gls_2, baseline2_gls = baseline(true_x[true_x[:].>=roi_glass[1,2]],data[true_x[:].>=roi_glass[1,2],2],roi_glass[2:end,:],"gcvspline",p=smo_water)
	    baseline_gls = [baseline1_gls;baseline2_gls]
	    bas_matrix[:,2] = [ycorr_gls_1;ycorr_gls_2]

	    # DISPLAYING INTERMEDIATE PLOTS
	    if plot_intermediate_show == "yes"
	        figure()

	        title("Raw spectra and baseline")
	        xlabel(L"Raman shift, cm$^{-1}$")
	        ylabel("Intensity, counts")

	        plot(true_x,data[:,1],color="black",label="Ctx file "*liste[i,1])
	        plot(true_x,data[:,2],color="blue",label="Glass file "*liste[i,2])

	        plot(true_x,baseline_ctx,color="red")
	        plot(true_x,baseline_gls,color="orange")
	        legend(loc="best",fancybox=true)
			show()

	        figure()

	        title("Raw spectra and baseline")
	        xlabel(L"Raman shift, cm$^{-1}$")
	        ylabel("Intensity, counts")

	        plot(true_x,data[:,1],color="black",label="Ctx file "*liste[i,1])
	        plot(true_x,data[:,2],color="blue",label="Glass file "*liste[i,2])

	        plot(true_x,baseline_ctx,color="red")
	        plot(true_x,baseline_gls,color="orange")
	        legend(loc="best",fancybox=true)
			xlim(2600,4000)
			ylim(0,maximum(data[2600 .<true_x.<4000,:]))
	        show()
	    end

	    if plot_intermediate_show == "yes"
	        figure()
	        title("Frequency shift correction")
	        			plot(true_x[roi_xshift_low.<true_x.<roi_xshift_high],bas_matrix[roi_xshift_low.<true_x.<roi_xshift_high,1]./maximum(bas_matrix[roi_xshift_low.<true_x.<roi_xshift_high,1]),color="blue",label="shifted")
	        			plot(true_x[roi_xshift_low.<true_x.<roi_xshift_high],bas_matrix[roi_xshift_low.<true_x.<roi_xshift_high,2]./maximum(bas_matrix[roi_xshift_low.<true_x.<roi_xshift_high,2]),color="red",label="reference")
	    end

	    # AUTOMATIC ADJUSTEMENT OF ANY X SHIFT BETWEEN THE SPECTRA, the glass spectrum is the reference
	    ~, bas_matrix[:,1], p = xshift_correction(true_x,bas_matrix[:,1],true_x[roi_xshift_low .< true_x .<roi_xshift_high], data[roi_xshift_low .< true_x .<roi_xshift_high,2],data[roi_xshift_low .< true_x .<roi_xshift_high,1])

	    # DISPLAYING INTERMEDIATE PLOTS
	    if plot_intermediate_show == "yes"
	        			plot(true_x[roi_xshift_low.<true_x.<roi_xshift_high],bas_matrix[roi_xshift_low.<true_x.<roi_xshift_high,1]./maximum(bas_matrix[roi_xshift_low.<true_x.<roi_xshift_high,1]),color="green",label="corrected")
	        			plot(true_x[roi_xshift_low.<true_x.<roi_xshift_high],bas_matrix[roi_xshift_low.<true_x.<roi_xshift_high,2]./maximum(bas_matrix[roi_xshift_low.<true_x.<roi_xshift_high,2]),color="red",label="reference")
	        legend(loc="best",fancybox=true)
	        xlabel(L"Raman shift, cm$^{-1}$")
	        ylabel("Intensity, counts")
	    	show()

	        figure()
	        title("Spectra after the baseline and X shift corrections")
	        plot(true_x,bas_matrix[:,1],color="black",label="Cristal file "*liste[i,1])
	        plot(true_x,bas_matrix[:,2],color="blue",label="Glass file "*liste[i,2])
	        xlabel(L"Raman shift, cm$^{-1}$")
	        ylabel("Intensity, counts")
	        legend(loc="best",fancybox=true)
			show()
	    end

	    ######################################################################
	    # CRYSTAL SIGNAL REMOVAL ALGORITHM
	    # The following code creates new spectra with mixed glass-ctx signals
	    # Those spectra are used in FastICA to recovert the pure glass and ctx components

	    # READING THE GOOD VALUES FOR THE ITERATIVE PROCESS
	    nb_iter = liste[i,4] # the number of iterations
	    K_start = liste[i,5] # for the iteration, starting mixing point (fraction of ctx signal)
	    K_increment = liste[i,6] # for the iteration, increamentation mixing point (fraction of ctx signal)

	    # CREATING USEFULL FUNCTIONS
	    function f1(k,ctx,g_c)
	        g = (g_c - k.*ctx)./(1.0-k)
	        return g
	    end

	    function f2(k,ctx,g) # maybe we do not need this one ?
	        g_c = k.*ctx+g
	        return g_c
	    end

	    # WE ARE ONLY RECODING SIGNALS BELOW 1300 cm-1 (strictly)
	    x_for_ica = true_x[true_x.<shutdown]
	    y_for_ica = ones(size(x_for_ica,1),nb_iter)

	    # LAST FEW THINGS BEFORE ITERATIONS
	    results = zeros(size(bas_matrix)) # the final matrix with mixed signals
	    K = K_start # we start with K_start ctx signal

	    for iter = 1:nb_iter # iterative removal, should be good in a few iterations
	        if iter == 1
	            results[:,:] = bas_matrix[:,:]
	            results[:,1] = K*bas_matrix[:,1]
	        else
	            results[:,1] = K.*bas_matrix[:,1]
	            results[:,2] = f1(K,bas_matrix[:,1],bas_matrix[:,2])
				#results[:,2] = f2(K,bas_matrix[:,1],bas_matrix[:,2])
	            results[true_x.>shutdown,2] = bas_matrix[true_x.>shutdown,2]
	        end

	        if plot_mixing_show == "yes" # for showing intermediate plots
				if iter == 1
				    figure()
				    title("The generation of spectra for FastICA")
		            subplot(311)
		            title("Ctx, cristal signals",fontsize=18,fontname="Arial")
		            plot(true_x,results[:,1],color=[iter/(nb_iter+0.0),0.,0.])
		            xlim(0,shutdown)
		            subplot(312)
		            title("G, unmixed glass signals = (G_C - K*Ctx)/(1-K)",fontsize=18,fontname="Arial")
		            ylabel("Intensity, counts",fontsize=18,fontname="Arial")
		            xlim(0,shutdown)
		            plot(true_x,results[:,2],color=[iter/(nb_iter+0.0),0.,0.])
		            subplot(313)
		            title("G_C, mixed glass signals = G + K*CTX",fontsize=18,fontname="Arial")
		            plot(true_x,results[:,1]+(1-K).*results[:,2],color=[iter/(nb_iter+0.0),0.,0.])
		            xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
		            xlim(0,shutdown)
					tight_layout()
				else
					subplot(311)
					plot(true_x,results[:,1],color=[iter/(nb_iter+0.0),0.,0.])
					subplot(312)
					plot(true_x,results[:,2],color=[iter/(nb_iter+0.0),0.,0.])
					subplot(313)
					plot(true_x,results[:,1]+(1-K).*results[:,2],color=[iter/(nb_iter+0.0),0.,0.])
				end

	        end

			y_for_ica[:,iter] = results[true_x.<shutdown,2] # old way: actually works well
			#y_for_ica[:,iter] = K*bas_matrix[true_x.<shutdown,1] + bas_matrix[true_x.<shutdown,2]# New trial: only positive addition: may improve NMF

	        # INCREMENTING K
	        K = K + K_increment

	    end
		if plot_mixing_show == "yes"
			show()
		end

		# We test the algorithm wanted by the user
		if algorithm == "FastICA"

		    # PREPROCESSING (STANDARDIZATION)
		    X_scaler = preprocessing[:StandardScaler]()
		    X_scaler[:fit](y_for_ica)
		    y_for_ica_sc = X_scaler[:transform](y_for_ica)

	    	# USING FastICA ON THE G SIGNALS
	    	model = decomposition[:FastICA](n_components=2, fun = "cube")
	    	S = model[:fit_transform](y_for_ica_sc)

	    	# TREATING THE OUTPUTS, FastICA can return negative components so we treat the S array for that
	    	S_corr = S[:,:]
	    	for col = 1:2
	        	if maximum(S[990. .< x_for_ica .< 1050,col]) < maximum(S[1200 .< x_for_ica .<shutdown,col]) # checking if something is negative
	            	S_corr[:,col] = -S[:,col]
	        	end
	        	S_corr[:,col] = (S_corr[:,col]-minimum(S_corr[:,col])) / maximum(S_corr[:,col]-minimum(S_corr[:,col])) # normalisation
	    	end
		elseif algorithm == "NMF"

		    # PREPROCESSING (STANDARDIZATION)
		    X_scaler = preprocessing[:MinMaxScaler]()
		    X_scaler[:fit](y_for_ica)
		    y_for_ica_sc = X_scaler[:transform](y_for_ica)

			# with sklearn
	    	model = decomposition[:NMF](n_components=2,init="nndsvda")
	    	S_corr = model[:fit_transform](y_for_ica_sc)

			# with julia NMF implementation: nndsvda init + alspgrad algo
			# commented at 10/2018 because NMF is not compatible with Julia 0.7/1.0
			#r = nnmf(transpose(y_for_ica), 2; init=:nndsvda, alg=:alspgrad, maxiter=100, verbose=false)
			#H = r.H; W = r.W; S_corr = transpose(H./maximum(H,2))

		else
			error("Not implemented, choose between NMF and FastICA")
		end

	    # WHICH SIGNAL IS WAHT? Now we determine who is the glass? => the Boson peak is more intense!
	    if maximum(S_corr[40 .< x_for_ica .< 100,2]) < maximum(S_corr[40 .< x_for_ica .<100,1])
	        glass_ica = S_corr[:,1] - minimum(S_corr[:,1])
	        ctx_ica = S_corr[:,2] - minimum(S_corr[:,2])
	    else
	        glass_ica = S_corr[:,2] - minimum(S_corr[:,2])
	        ctx_ica = S_corr[:,1] - minimum(S_corr[:,1])
	    end

	    # FINAL RECORDING MATRIX
		fit = curve_fit(model_ctx_adj, glass_ica[x_for_ica.<scaling], bas_matrix[true_x .<scaling,2], [1.0])
		coef_scale = fit.param
	    spectra_final = [true_x [model_ctx_adj(glass_ica,coef_scale);bas_matrix[true_x.>=shutdown,2]]]

	    # FINAL FIGURE, SHOWED
	    figure(figsize=(15,15))

	    subplot(3,2,(1,2)) # UPPER PLOT
		title("Glass sp. no $(liste[i,2]), Baseline correction",fontsize=18,fontname="Arial")
	    plot(true_x,bas_matrix[:,1],color="blue",label="Cristal")
	    plot(true_x,bas_matrix[:,2],color="black",label="Glass")
	    ylabel("Intensity, counts",fontsize=18,fontname="Arial")
	    xlim(0,4000)
	    legend(loc="best",fancybox=true)

	    subplot(3,2,3)
	    title("results from the algorithm "*algorithm,fontsize=18,fontname="Arial")
	    plot(x_for_ica,glass_ica,color="black")
	    plot(x_for_ica,ctx_ica,color="blue")
	    ylabel("Intensity, counts",fontsize=18,fontname="Arial")
	    xlim(0,shutdown)

	    subplot(3,2,4)
	    title("Final fit, silicate region",fontsize=18,fontname="Arial")
	    plot(x_for_ica, model_ctx_adj(glass_ica,coef_scale),color="green",label="Glass, corrected")
	    plot(x_for_ica, bas_matrix[true_x.<shutdown,2],color="black",label="Glass, initial")
	    xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
	    ylabel("Intensity, counts",fontsize=18,fontname="Arial")
	    xlim(0,shutdown)

	    subplot(3,2,(5,6)) # LOWER PLOT
		title("Final fit, entire spectra",fontsize=18,fontname="Arial")
		plot(true_x,bas_matrix[:,2],color="black",label="Glass, initial")
	    plot(spectra_final[:,1],spectra_final[:,2],color="green",label="Glass, corrected")
	    legend(loc="best",fancybox=true)
	    xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
	    ylabel("Intensity, counts",fontsize=18,fontname="Arial")
	    xlim(0,4000)

		tight_layout()
	    # AND IF THE USER WANTS TO SAVE THE FIGURE, WE OUTPUT THEM ALSO IN OUT_PATH
	    if save_fig_switch == "yes"
	        savefig(out_path*liste[i,2]*".pdf")
	    end

		if plot_final_show == "yes"
			show()
		else
			close()
		end

	    # WE WRITE THE CORRECTED SPECTRA IN THE CORRECT OUTPUT PATH (OUT_PATH)
	    writedlm(out_path*liste[i,2],spectra_final)

	end
end


text:
module LibCImGui

import Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("CImGui was not build properly. Please run Pkg.build(\"CImGui\").")
end
include(depsjl_path)
# Module initialization function
function __init__()
    check_deps()
end

using CSyntax.CEnum

const Ctm = Base.Libc.TmStruct
const Ctime_t = UInt
const Cclock_t = UInt
export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "..", "gen", "libcimgui_common.jl"))
include(joinpath(@__DIR__, "..", "gen", "libcimgui_api.jl"))
include(joinpath(@__DIR__, "..", "gen", "libhelper_api.jl"))

foreach(names(@__MODULE__, all=true)) do s
    if startswith(string(s), "ig") || startswith(string(s), "Im")
        @eval export $s
    end
end

end # module


text:
import Pkg; Pkg.activate(".")
using Revise

import FiveThreeOne
const F = FiveThreeOne

const MAXES_FILE = "maxes.yaml"
const BBB_ASSISTANCE_FILE = "assistance.yaml"

main_lifts_5pro(maxes_data, day, week, order) = F.main_lifts_5pro(maxes_data["training"][day], week, order)
main_lifts(maxes_data, day, week, order) = F.main_lifts_5pro(maxes_data["training"][day], week, order, false, maxes["e1rm"])
main_lifts_pr_sets(maxes_data, day, week, order) = F.main_lifts_5pro(maxes_data["training"][day], week, order, true, maxes["e1rm"])
boringbutbig_light(maxes_data, day, week, order) = F.boringbutbig_light(maxes_data["training"][day], week, order)

bbb = F.routine_from_file(
    days=["Squat", "Bench", "Deadlift", "Press"],
    main_lifts=[
        main_lifts_5pro,
        boringbutbig_light,
    ],
    maxes_file=MAXES_FILE,
    assistance_file=BBB_ASSISTANCE_FILE,
    order=F.Order351,
)

bbb(F.Week2)


text:
using PlotUtils
using Test
using Statistics: mean
using Dates
using Random
using StableRNGs

rng = StableRNG(42)

# TODO: real tests

# ----------------------
# colors

const C = RGBA{Float64}
const C0 = RGBA{PlotUtils.Colors.N0f8}

@testset "colors" begin

    @test plot_color(nothing) == C(0, 0, 0, 0)
    @test plot_color(false) == C(0, 0, 0, 0)
    @test_throws ErrorException plot_color(true)

    @test plot_color(:red) == parse(C, :red)
    @test plot_color("red") == parse(C, "red")
    @test_throws ErrorException plot_color("notacolor")

    @test plot_color(colorant"red") == C(1, 0, 0, 1)

    grad = cgrad()
    @test typeof(grad) == PlotUtils.ContinuousColorGradient
    @test plot_color(grad) === grad

    grad = cgrad([:red, "blue"])
    @test color_list(grad) == C[colorant"red", colorant"blue"]
    @test grad.values == collect(range(0, stop = 1, length = 2))

    grad = cgrad([:red, "blue"], alpha = 0.5)
    @test C0.(color_list(grad)) == C0[C(1, 0, 0, 0.5), C(0, 0, 1, 0.5)]
    @test grad.values == collect(range(0, stop = 1, length = 2))

    grad = cgrad([:red,:blue], [0,0.1,1])
    @test length(color_list(grad)) == 3
    @test grad.values == [0,0.1,1]

    cs = plot_color(rand(rng, 10))
    @test typeof(cs) == Vector{C}
    @test length(cs) == 10

    cs = plot_color(rand(rng, 4, 4))
    @test typeof(cs) == Matrix{C}
    @test length(cs) == 16
    @test size(cs) == (4, 4)

    cs = plot_color(rand(rng, 10), 0.5)
    @test typeof(cs) == Vector{C}
    @test length(cs) == 10
    for c in cs
        @test alpha(c) == 0.5
    end

    cs = plot_color(rand(rng, 4, 4), 0.5)
    @test typeof(cs) == Matrix{C}
    @test length(cs) == 16
    @test size(cs) == (4, 4)
    for c in cs
        @test alpha(c) == 0.5
    end
end

# ----------------------
# gradients

@testset "gradients" begin
    grad = cgrad(:inferno)
    @test length(grad) == 256
    @test RGB(grad.colors[1]) == RGB(0.001462, 0.000466, 0.013866)
    @test RGB(grad.colors[end]) == RGB(0.988362, 0.998364, 0.644924)
end

# ----------------------
# ticks

# Copied from Plots.is_uniformly_spaced to avoid dependency on recent version
# on Plots which is not used on Travis.
function is_uniformly_spaced(v; tol = 1e-6)
    dv = diff(v)
    maximum(dv) - minimum(dv) < tol * mean(abs.(dv))
end

@testset "ticks" begin
    @test optimize_ticks(-1, 2) == ([-1.0,0.0,1.0,2.0], -1.0, 2.0)
    dt1, dt2 = Dates.value(DateTime(2000)), Dates.value(DateTime(2100))
    @test optimize_datetime_ticks(dt1, dt2) == (
        [63113990400000, 63902908800000, 64691827200000, 65480745600000],
        ["2001-01-01", "2026-01-01", "2051-01-01", "2076-01-01"])

    @testset "small range" begin
        @testset "small range $x, $(i)ϵ" for x in exp10.(-12:12), i in -5:5
            y = x + i * eps(x)
            x, y = minmax(x, y)
            ticks = PlotUtils.optimize_ticks(x, y)[1]
            @test issorted(ticks)
            @test all(x .<= ticks .<= y)
            # Fails:
            # @test allunique(ticks)
        end
    end

    function test_ticks(x, y, ticks)
        @test issorted(ticks)
        @test all(x .<= ticks .<= y)
        if x < y
            @test length(ticks) >= 2
            @test is_uniformly_spaced(ticks)
        end
    end

    @testset "fixed ranges" begin
        @testset "fixed range $x..$y" for (x, y) in [(2, 14),(14, 25),(16, 36),(57, 69)]
            test_ticks(x, y, optimize_ticks(x, y)[1])
            test_ticks(-y, -x, optimize_ticks(-y, -x)[1])
        end
    end

    @testset "random ranges" begin
        r = [minmax(rand(rng, -100:100, 2)...) .* 10.0^i for _ = 1:10, i = -5:5]
        @testset "random range $x..$y" for (x, y) in r
            test_ticks(x, y, optimize_ticks(x, y)[1])
        end
    end

    # issue 86
    let x = -1.0, y = 13.0
        test_ticks(x, y, optimize_ticks(x, y, k_min = 4, k_max = 8)[1])
    end

    @testset "digits" begin
        @testset "digits $((10^n) - 1)*10^$i" for n in 1:9, i in -9:9
            y0 = 10^n
            x0 = y0 - 1
            x, y = (x0, y0) .* 10.0^i
            ticks = optimize_ticks(x, y)[1]
            test_ticks(x, y, ticks)
        end
    end
end

# ----------------------
# adapted grid

@testset "adapted grid" begin
    f = sin
    int = (0, π)
    xs, fs = adapted_grid(f, int)
    l = length(xs) - 1
    for i in 1:l
        for λ in 0:0.1:1
            # test that `f` is well approximated by a line
            # in the interval `(xs[i], xs[i+1])`
            x = λ * xs[i] + (1 - λ) * xs[i + 1]
            y = λ * fs[i] + (1 - λ) * fs[i + 1]
            @test y ≈ f(x) atol = 1e-2
        end
    end

    int = (2, 2)
    xs, fs = adapted_grid(f, int)
    @test xs == [2]
    @test fs == [f(2)]

    int = (2, 1)
    @test_throws ArgumentError adapted_grid(f, int)
end

@testset "zscale" begin
    bkg = 30 .* randn(rng, 8192) .+ 1000
    data = bkg .+ 100 .* randn(rng, 8192) .+ 2500
    defects = rand(rng, CartesianIndices(bkg), 500)
    data[defects] .= rand(rng, [0, 1e7], 500)
    cmin, cmax = zscale(data)
    # values calculated using IRAF
    @test cmin ≈ 2841.586 atol = 1e-3
    @test cmax ≈ 4171.648 atol = 1e-3
    @test cmin > minimum(data)
    @test cmax < maximum(data)

    data = vcat(1:100)
    cmin, cmax = zscale(data)
    @test cmin == 1
    @test cmax == 100

    # Make sure output is finite
    data = vcat(0:999, NaN)
    cmin, cmax = zscale(data)
    @test cmin == 0
    @test cmax == 999
end


text:
mutable struct PlotlinesModel{T₁ <: AbstractVector} <: AbstractModel
    data::T₁
end

function get_data(model::PlotlinesModel)
    model.data
end

function set_data!(model::PlotlinesModel, data::AbstractVector)
    model.data = data
end


text:
"""
    x,y,z ∈ 1,2,3
    x,y,z ∈ 1:3
    x,y,z ∈ apples, oranges, grapes
"""
function operatoreval(command, logicset::LogicalCombo; verbose=true)
    counter!()
    logicsetcopy = deepcopy(logicset)

    #println("operatoreval($command)")
    (sum(logicset[:]) == 0) && return logicset
    occursin(r"\{\{.*\}\}", command) &&
      return operatorspawn(command, logicset, verbose=verbose)

    n = 1:sum(logicset[:])

    # convert a = b|c to a |= b,c
    if occursin("|", command) & occursin(r"(\b| )[!|]*=+[!|]*(\b| )", command)
      command = replace(command, "|"=>",")
      command = replace(command, r",*=+"=>"|=")
    end

    command = replace(command, r"^(.*?)[ ]*([!&])[ ]*$"=>s"\1 \2 \1")
    command = replace(command, r"^[ ]*([!&])[ ]*(.*?)$"=>s"\2 \1 \2")

    m = match(r"^(.*?)(([><=|!^&]{1,3}))(.*?)(\{([0-9]+),?([0-9]+)*\})?$",command)
    left, right, operator, nmin, nmax  = m.captures[[1,4,3,6,7]]

    (nmin === nothing)  && (nmax === nothing)  && (nrange = 1:999)
    !(nmin === nothing) && (nmax === nothing)  && (nrange = integer(nmin):integer(nmin))
    !(nmin === nothing) && !(nmax === nothing) && (nrange = integer(nmin):integer(nmax))

    leftarg  = strip.(split(replace(left, " "=> ""),  r"[,&|!]"))
    rightarg = strip.(split(replace(right, " "=> ""), r"[,&|!]"))

    (((length(collect(eachmatch(r"[+-]", left)))>0) &&
      (length(collect(eachmatch(r"[*\\÷]", left)))>0)) ||
    ((length(collect(eachmatch(r"[+-]", right)))>0) &&
      (length(collect(eachmatch(r"[*\\÷]", right)))>0))) && verbose &&
      println(
      Crayons.Box.GREEN_FG("\nWarning: Math operators are evaluated from left to right. 3 + 1 * 2 = 8"))

    leftvals  = hcat([grab(L, logicset, command=command) for L in leftarg]...)
    rightvals = hcat([grab(R, logicset, command=command) for R in rightarg]...)

    if operator == ""
        throw("No operator found!")

    elseif operator == "|"
        lcheck = [any(isodd.(leftvals[i,:])) for i in n]
        rcheck = [any(isodd.(rightvals[i,:])) for i in n]
        ℧Δ = [all(lcheck[i,:]) | all(rcheck[i,:]) for i in n]

    elseif operator == "&"
        lcheck = [all(isodd.(leftvals[i,:])) for i in n]
        rcheck = [all(isodd.(rightvals[i,:])) for i in n]
        ℧Δ = [all(lcheck[i,:]) & all(rcheck[i,:]) for i in n]

    elseif operator == "^"
        ℧Δ = (isodd.(leftvals[:,1]) .& iseven.(rightvals[:,1])) .|
             (iseven.(leftvals[:,1]) .& isodd.(rightvals[:,1]))

    elseif operator == "!"
        lcheck = [all(iseven.(leftvals[i,:])) for i in n]
        rcheck = [all(iseven.(rightvals[i,:])) for i in n]
        ℧Δ = [all(lcheck[i,:]) & all(rcheck[i,:]) for i in n]

    elseif operator == "^="
         lcheck = [sum(leftvals[i,j]  .== rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
         rcheck = [sum(rightvals[i,j] .== leftvals[i,:])  for i in n, j in 1:size(rightvals)[2]]
         ℧Δ = [sum(lcheck[i,:]) + sum(rcheck[i,:]) == 2 for i in n]

    elseif operator  ∈ ["==","="]
        lcheck = [all(leftvals[i,j] .== rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == "!="
        lcheck = [all(leftvals[i,j] .!= rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator ∈ ["|=", "=|"]
        lcheck = [any(leftvals[i,j] .== rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [sum(lcheck[i,:]) ∈ nrange for i in n]

    elseif operator ∈ ["|!=", "|!", "!|", "!|=", "=|!", "|=!"]
        lcheck = [any(leftvals[i,j] .!= rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [sum(lcheck[i,:]) ∈ nrange for i in n]

    elseif operator == "<="
        lcheck = [all(leftvals[i,j] .<= rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == "<"
        lcheck = [all(leftvals[i,j] .< rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == "<<"
        lcheck = [all(leftvals[i,j] .< rightvals[i,:].-1) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == ">="
        lcheck = [all(leftvals[i,j] .>= rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == ">"
        lcheck = [all(leftvals[i,j] .> rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == ">>"
        lcheck = [all(leftvals[i,j] .> rightvals[i,:].+1) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]
    end

    logicsetcopy.logical[logicsetcopy[:]] = ℧Δ

    logicsetcopy
end

function grab(argument::AbstractString, logicset::LogicalCombo; command = "")

  (argument[1:1]=="'") && return fill(replace(argument, "'"=>""), length(logicset[:]))

  # matcher = r"^([a-zA-z][a-zA-z0-9.]*)*([0-9]+)*([+\-*/])*([a-zA-z][a-zA-z0-9]*)*([0-9]+)*$"
  # matcher = r"^(?:([a-zA-z][a-zA-z0-9.]*)|([0-9]+))\s*([+\-*/÷])*\s*(.*?)$"
  matcher = r"^(.*?)\s*([+\-*/÷%]*)\s*(?:([a-zA-z][a-zA-z0-9.]*)|([0-9]+))$"

  m = match(matcher, argument)
  (m==nothing) && throw("Argument $argument could not be parsed in $command")

  v1, o1, v2, n2 = m.captures

  !(v2 === nothing) && (right  = logicset[:,Symbol(v2),true])
  !(n2 === nothing) && (right  = fill(integer(n2), nfeasible(logicset)))

  (length(o1)>1) && throw("Math Operators only accepts singletons (+-÷*\\) and no negative multiplication (a*-b)!")
  (o1 == "")   && return right

  !(v1 === nothing) && (left = grab(v1, logicset))

  (o1 == "+") && return left .+ right
  (o1 == "-") && return left .- right
  (o1 == "/") && return left ./ right
  (o1 == "÷") && return left .÷ right
  (o1 == "*") && return left .* right
  (o1 == "%") && return left .% right
end


text:
using DataDeps

register(DataDep(
  "japan-emmigration",
  "http://www.stat.ufl.edu/~winner/data/japanemg.txt",
  "http://www.stat.ufl.edu/~winner/data/japanemg.dat",
  "dc3f7d7cb5dec406251776b17208b1bee0db923203e1c30fc0547298f5559cb8",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=2,
        feature_indices=3:7,
        id_indices=1,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data