text:
using Revise

using Dates, DataFrames, BenchmarkTools, Random

import Base: ==

struct JoinableRange{Tuniformrange<:AbstractRange, Tval}
  uniformrange::Tuniformrange
  first::Tval
  last::Tval
end

JoinableRange(uniformrange::AbstractRange) = JoinableRange(
  uniformrange, minimum(uniformrange), maximum(uniformrange))

==(jr::JoinableRange, v) = ifelse((jr.first ≤ v) & (v ≤ jr.last), true, false)
==(jr::JoinableRange, ::Missing) = missing
==(v, jr::JoinableRange) = jr == v

function testfuzzydatejoin(N)
  df1 = DataFrame(grp=rand(collect(1:N), N))

  startdt = rand(collect(Date(1900,1,1):Day(1):Date(2019,12,31)), N)
  enddt = (startdt->rand(startdt:Day(1):Date(2020,12,31))).(startdt)

  rngs = ((startdt,enddt)->startdt:Day(1):enddt).(startdt,enddt)
  df1.jr = JoinableRange.(rngs)

  df2idx = shuffle!(collect(1:N))

  df2=DataFrame(grp=df1.grp[df2idx]|>deepcopy, dt = (jr->rand(jr.uniformrange)).(df1.jr[df2idx]))
  describe(df2) |> println

  df12 = innerjoin(df1,df2, on=[:grp, :jr=>:dt])

  #describe(df12) |> println

end

#@time testfuzzydatejoin(10^6)

using DataFrames, Dates

function mergecrspcompkernel(;crsp=error("crsp is required"), comp=error("comp is required"))
  ###Merge one permno at a time, linking each crsp row to a compid (comp row)
  scrsps = groupby(crsp, :permno)
  scomps = groupby(comp, :permno)
  comp.compid = 1:size(comp,1) |> collect
  crsp.compid = Vector{Union{Int, Missing}}(undef, size(crsp,1))

  #faster to work within the rows for a particular permno
  Threads.@threads for permno ∈ keys(scomps)
    scomp = scomps[permno]
    crsppermnokey = permno |> Tuple
    (!haskey(scrsps,crsppermnokey)) && continue
    #because the key originates from comp, cant use it directly on crsp- instead need its value
    scrsp = scrsps[crsppermnokey]

    #the idea here is to select all of the crsp rows between the effective date
    #and the end date of the comp row
    for r ∈ eachrow(scomp)
      targetrows = view(scrsp, (r.effdate .≤ scrsp.date) .& (scrsp.date .≤ r.enddate), :compid)
      if length(targetrows) > 0
        #under no circumstances should this be assigned, since date and permno should be unique
        (targetrows .=== missing) |> all || error(
          "multiple crsp rows found for a particular compid!!")
        targetrows .= r.compid
      end
    end
  end

  crsp = crsp[crsp.compid .!== missing, :]

  #now merge the linked rows
  rename!(comp, :permno=>:lpermno)
  panel = innerjoin(crsp, comp, on=:compid)

  #sanity checks
  @assert size(panel,1) == size(crsp,1)
  @assert (crsp.permno .== panel.permno) |> all
  @assert (crsp.date .== panel.date) |> all

  return panel
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module MyPkg0

using DataFrames

export sog, tmpsog, zo, zo1, zo3, newsog, newsog2, newsog3, newsog4, select_if
"""
    sog(x)

The function `sog()` takes a `Vector` of any type `x` as its input with any possible length, and will return a `Bool Vector` with the same length as input. Values
 `true` or `false` in the output depends on whether one element in the `Vector` is the start of a group or not. The sog() will return
 'nothing' if the length of input is 0. As any type can be applied in the input, be careful! To check whether two elements `a` and `b` are treated
as the same, please call `isequal(a, b)` to see the result.

When taking `Vector` of `Vectors`, the length of output should the same as the length of every `Vector`, and what `sog()` do is to compare
one element with the previous one on the same position for each `Vector`.

Examples:
≡≡≡≡≡≡≡≡≡≡≡≡≡≡

    sog([1,2,2,2,1,3,3,1,1]) returns: [true, true, false, false, true, true, false, true, false]
    sog() returns: nothing
    sog(["a", "a", "a", 'a']) returns: [true, false, false, true]
    sog([[1, 1, 1, 2, 2], [1.5, 1.5, 20.0, 3.0, 3.0]] returns: [true, false, true, true, false]

"""
function sog(x=[])
    if x == []
        return nothing
    end

    iVector = x
    len = length(x)
    #lenv = length(x[1])

    if typeof(iVector[1]) <: Vector{}           #For Vector of Vectors
        lenv = length(x[1])
        if lenv == 0
            return nothing
        end

        if lenv == 1
            return [true]
        end

        if lenv >= 2
            r0 = map(sog, iVector)
            rInt = Vector{Vector{Int64}}(r0)
            oVector = Vector{Bool}(replace(!iszero, sum(rInt)))

            return oVector
        end
    end

    if len == 1
        return [true]
    end

    r = iVector[1:(len-1)] .=== iVector[2:len]
    r1 = Vector{Bool}(replace(iszero, r))
    oVector = Vector{Bool}(append!([1], r1))

    return oVector
end

function zo(re, xi)
    for j in 2:length(re)
        re[j] = re[j]==1 ? 1 : !(xi[j]===xi[j-1])
    end
    return re
end

function tmpsog(x)
    re = zeros(Bool, length(x[1]))
    re[1] = 1
    for i in 1:length(x)
        re = zo(re, x[i])
    end
    return re
end

function zo1(re, xi)
    for j in 2:length(re)
        @inbounds re[j] = re[j]==1 ? 1 : !(xi[j]===xi[j-1])
    end
    return re
end

function newsog(x)
    re = zeros(Bool, length(x[1]))
    re[1] = 1
    for i in 1:length(x)
        re = zo1(re, x[i])
    end
    return re
end

function ini0(x1)
    len = length(x1)
    re = zeros(Bool, len)
    re[1] = 1
    for j in 2:len
        @inbounds re[j] = !(x1[j]===x1[j-1])
    end
    return re
end

function newsog2(x)
    re = ini0(x[1])
    for i in 2:length(x)
        re = zo1(re, x[i])
    end
    return re
end

function zo3(re, xi)
    Threads.@threads for j in 2:length(re)
        @inbounds re[j] = re[j]==1 ? 1 : !(xi[j]===xi[j-1])
    end
    return re
end

function newsog3(x)
    re = ini0(x[1])
    for i in 2:length(x)
        re = zo3(re, x[i])
    end
    return re
end

function newsog4(x, orders = eachindex(x))
    if length(eachindex(x)) != length(orders)
        println("Please specify the order of columns correctly! Number of columns is the length of order.")
        return nothing
    end
    re = ini0(x[orders[1]])
    for i in orders[2:end]
        re = zo3(re, x[i])
    end
    return re
end

"""
    select_if(df, predicate, elementwise_or_not, any_or_all)

This function is aim to select `DataFrame` columns based on a predicate applied to the columns or a logical vector.

Keyword arguments

≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡

-  `df` is the DataFrame.

-  `predicate` can be a predicate function for columns or a `Bool` `Vector` with the length of number of columns.

-  `elementwise_or_not` takes a `Bool` value that is `true` in default. Change `elementwise_or_not` to `false` will
   let the predicate apply to the DataFrame columnwise.

-  `any_or_all` also takes a `Bool` value with a default value `true` and it determinesto apply any or all predicate
   to the DataFrame when `elementwise_or_not` = `true`.

Examples

≡≡≡≡≡≡≡≡≡≡≡≡≡≡

```
julia> df1 = DataFrame([missing 1 2 3; missing 2 missing 4; 1 4 2 5], :auto)
3×4 DataFrame
 Row │ x1       x2      x3       x4
     │ Int64?   Int64?  Int64?   Int64?
─────┼──────────────────────────────────
   1 │ missing       1        2       3
   2 │ missing       2  missing       4
   3 │ 1             4        2       5

julia> select_if(df1, [false, true, true, true])
3×3 DataFrame
 Row │ x2      x3       x4
     │ Int64?  Int64?   Int64?
─────┼─────────────────────────
   1 │      1        2       3
   2 │      2  missing       4
   3 │      4        2       5

julia> select_if(df1, !ismissing, 1, 0)
3×2 DataFrame
 Row │ x2      x4
     │ Int64?  Int64?
─────┼────────────────
   1 │      1       3
   2 │      2       4
   3 │      4       5

julia> pre1 = (function(col) return mean(skipmissing(col)) >= 4 end)

julia> select_if(df1, pre1, 0)
3×1 DataFrame
 Row │ x4
     │ Int64?
─────┼────────
   1 │      3
   2 │      4
   3 │      5

```
"""
function select_if(df::DataFrame, predicate, elementwise_or_not = true, any_or_all = true)
    if typeof(predicate) == Vector{Bool}
        indices = predicate
    elseif Bool(elementwise_or_not) == true
        if any_or_all == true
            indices = map(x -> any(predicate, x), eachcol(df))
        else
            indices = map(x -> all(predicate, x), eachcol(df))
        end
    else
       indices = predicate.(eachcol(df))
    end
    return df[:, indices]
end

#Base.names(df::AbstractDataFrame, fun::Function) = typeof(fun) isa Regex ? filter!(fun, names(df)) : names(df, fun.(eachcol(df)))

end

ext: jl


text:
struct parametrization_tendencies{T<:AbstractFloat}

#Just a placeholder

end

"""
Compute physical parametrization tendencies
"""
function parametrization_tendencies{T}(utend, vtend, ttend, qtend) where {T<:AbstractFloat}

end


text:
# %% value of τ demonstration
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

d = 3
τ1 = 1
τ2 = 17
τ3 = 3τ2
θ = -30
e = 10
lo = Systems.lorenz([0, 10, 0.0])
A = trajectory(lo, 1000; Ttr=1000, reltol = 1e-12, abstol = 1e-12)
w = A[:, 1]
close("all")
N = 5000
fig = figure(figsize = (figx, 0.75figy))

for (i, τ) in enumerate((0, τ1, τ2, τ3))
    R = τ == 0 ? A : embed(w, d, τ)
    ax = subplot(1, 4, i, projection = "3d")
    ax.plot(R[1:N, 1], R[1:N, 2], R[1:N, 3],
    color = i == 1 ? "C0" : "C1", lw = 1,
    zorder = 1)
    ax.view_init(elev=e, azim=θ)
    ax.set_title(i == 1 ? "original" : "\$R:\\,\\tau=$(τ)\$", pad = -40)
    for a in (ax.xaxis, ax.yaxis, ax.zaxis); a.set_ticklabels([]); end
    dj = 167
    for j in (827, )
        χ, ψ, ζ = R[j:dj:j+dj, 1], R[j:dj:j+dj, 2], R[j:dj:j+dj, 3]
        ax.scatter3D(χ, ψ, ζ, depthshade = false, c = "C2", zorder = 99, s = 50)
    end
    ax.dist = 8
end

# ax = subplot(1, 4, 1, projection = "3d")
# ax.plot3D(x[1:N], y[1:N], z[1:N], lw = 1.0, color = "C1")
# ax.view_init(elev=e, azim=θ)
# for a in (ax.xaxis, ax.yaxis, ax.zaxis); a.set_ticklabels([]); end

fig.tight_layout(pad=0.1)
fig.subplots_adjust(wspace = 0.0001)#, bottom = -0.1, top = 1, left = 0.0, right = 1)
# wsave(plotsdir("6", "taudemonstration"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# const fmt = "[{date} | {level} | {name}]: {msg}"
const fmt = x -> Dates.format(x, "yyyy-mm-dd HH:MM:SS")
const ui_interface = Ref{Function}((x,y) -> nothing)
const logging = Dict{String,Any}()
logging["log_to_file"] = false
logging["file_logger"] = SimpleLogger()

function csinfo(msg, suppress_messages::Bool = false)
    log_to_file = logging["log_to_file"]
    file_logger = logging["file_logger"]
    msg = string(fmt(Dates.now())) * " : " * msg
    !suppress_messages && @info(msg)
    ui_interface[](msg, :info)
    if log_to_file
        with_logger(file_logger) do
            @info(msg)
        end
    end
end

function cswarn(msg)
    log_to_file = logging["log_to_file"]
    file_logger = logging["file_logger"]
    msg = string(fmt(Dates.now())) * " : " * msg
    @warn(msg)
    ui_interface[](msg, :warn)
    if log_to_file
        with_logger(file_logger) do
            @warn(msg)
        end
    end
end

function update_logging!(cfg)

    log_level = cfg["log_level"]
    log_file = cfg["log_file"]

    if log_level in DEBUG
        Logging.LogLevel(Logging.Debug)
    elseif log_level in WARNING
        Logging.LogLevel(Logging.Warn)
    end

    if log_file != "None"
        logging["file_logger"] = SimpleLogger(open(log_file, "w+"))
        logging["log_to_file"] = true
        csinfo("Logs will recorded to file: log_file", cfg["suppress_messages"] in TRUELIST)
    else
        logging["file_logger"] = SimpleLogger()
        logging["log_to_file"] = false
    end
end


text:
using LinearAlgebra
using Combinatorics

using Polyhedra

# Hypercube from `α * ones(n)` to `β * ones(n)`.
function vhypercubetest(lib::Polyhedra.Library, n, α, β)
    row = 0
    V = α * ones(Int, 2^n, n)
    for k in 0:n
        for p in combinations(1:n, k)
            row += 1
            V[row, p] = β * ones(Int, length(p))
        end
    end
    v = vrep(V)

    poly = polyhedron(v, lib)

    h = hrep([Matrix(1I, n, n); Matrix(-1I, n, n)], [β * ones(Int, n); -α * ones(Int, n)])

    inequality_fulltest(poly, h)
    generator_fulltest(poly, v)
end


text:
#
# This file is part of the GenServers.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#
"""
    GenServers

Implements a generic server `Actors` protocol.

It lets you write an implementation module with a
user API for a server and callback functions which
are called by the server if certain requests are
issued.

It lets you write purely sequential code while the
`:genserver` actor cares about the concurrency and
distributed part. This simplifies greatly the
development of servers.

The current stable, registered version is installed with
```julia
pkg> add GenServers
```

The development version is installed with:
```julia
pkg> add "https://github.com/JuliaActors/GenServers.jl"
```
"""
module GenServers

"Gives the package version."
const version = v"0.2.0"

using Reexport, Distributed

@reexport using Actors
import Actors: spawn, _ACT, Init, Call, Cast

include("genserver.jl")
include("protocol.jl")
include("callbacks.jl")

export genserver

end


text:
module TypeVars

export typevar, settypevar

const TYPEVARS = IdDict()

function typevar_key(::Type{T}) where T
    T
end

"""
    settypevar(::Type{T}, obj)

Associate the object with the type. The actual key used is
`typevar_key(T)`, which defaults to identity.
"""
function settypevar(::Type{T}, obj) where T
    get!(TYPEVARS, typevar_key(T), obj)
end

"""
    typevar(::Type)

Return the type varable associated with the type. See also `settypevar`.
"""
@generated function typevar(::Type{T}) where T
    t = TYPEVARS[typevar_key(T)]
    :( $t )
end

end


text:
@testset "TSP" begin

using STMO.TSP

X = [1 1; 1 3; 2 2; 2 1]

tsp = TravelingSalesmanProblem(X)

@test tsp isa TravelingSalesmanProblem{Int}

tour = [4, 1, 2, 3]

n = length(tour)

@test n == 4

@test Set(cities(tsp)) == Set(tour)

@test isvalid(tsp, tour)
@test !isvalid(tsp, [1, 2, 3])
@test !isvalid(tsp, [1, 2, 3, 4, 5])
@test !isvalid(tsp, [3, 1, 2, 3])

@test dist(tsp, 2, 3) ≈ sqrt(2)
@test computecost(tsp, tour) ≈ 2 + sqrt(2) + 1 + 1

@testset "swapping" begin
    cost = computecost(tsp, tour)

    for i in 1:n
        for j in 1:n
            Δc = deltaswapcost(tsp, tour, i, j)
            swap!(tour, i, j)
            @test computecost(tsp, tour) ≈ cost + Δc
            swap!(tour, i, j)
        end
    end
end

@testset "flipping" begin
    cost = computecost(tsp, tour)

    for i in 1:n
        for j in 1:n
            Δc = deltaflipcost(tsp, tour, i, j)
            flip!(tour, i, j)
            @test computecost(tsp, tour) ≈ cost + Δc
            flip!(tour, i, j)
        end
    end
end

end


text:
using MyInterpolations
using Base.Test

const _lin_interp = LinearInterpolation

@testset "Testing linear interporation" begin
    @testset "case with two grid points" begin
        grid = [1, 2]
        vals = [2, 0]
        itp = _lin_interp(grid, vals)

        @testset "scalar input" begin
            x, y = 1.25, 1.5
            @test isapprox(itp(x), y)
        end

        @testset "end points" begin
            for (x, y) in [(grid[1], grid[end]), (vals[1], vals[end])]
                @test itp(x) == y
            end
        end

        @testset "array inputs" begin
            xys = [
                ([4/3, 1.8], [4/3, 0.4]),
                (linspace(1, 2, 5), [2.0 - 0.5*i for i in 0:4])
            ]
            for (x, y) in xys
                @test isapprox(itp.(x), y)
            end
        end
    end

    @testset "general case" begin
        grid = 0:2
        vals = [1, 3, 2.5]
        itp = _lin_interp(grid, vals)

        @testset "scalar input" begin
            x, y = 1.2, 2.9
            @test isapprox(itp(x), y)
        end

        @testset "end points" begin
            for (x, y) in [(grid[1], vals[1]), (grid[end], vals[end])]
                @test itp(x) == y
            end
        end

        @testset "array inputs" begin
            xys = [
                ([0.5, 1.2], [2, 2.9]),  # Vector
                (grid, vals),  # UnitRange
                (0:2:2, vals[[1, 3]]),  # StepRange
                (0:0.5:1.5, [1, 2, 3, 2.75]),  # FloatRange
                (linspace(0, 1.5, 4), [1, 2, 3, 2.75])  # LinSpace
            ]
            for (x, y) in xys
                @test isapprox(itp.(x), y)
            end
        end
    end
end


text:
using BrainFlow

data = Array{Float64,3}(undef, 2, 4, 2)
data[:, :, 1] = [[6,3] [3,0] [1,5] [5,1]]
data[:, :, 2] = [[1,5] [5,1] [6,2] [2,2]]
data = permutedims(data, [3, 1, 2])
labels = [0.0, 1.0]

filters_true = [[-0.313406,-0.280803] [0.0792153, -0.480046]]
filters, eigvals = BrainFlow.get_csp(data, labels)

println("filters = ")
println(filters)

println("true filters =")
println(filters_true)


text:
using GAIO
using ForwardDiff
using StaticArrays

# example functions
# domain (-5,5)^2
function example_g1(x)
    fac1 = x[1]*x[1] + x[2] - 11
    fac2 = x[1] + x[2]*x[2] - 7

    return SVector(4*x[1]*fac1 + 2*fac2, 2*fac1 + 4*x[2]*fac2)
end

# domain (-5,5)^2
function example_g2(x)
    return SVector(x[1]*x[1]*x[1] - 3*x[1]*x[2]*x[2] - x[1] + 1/sqrt(2),
                    -x[2]*x[2]*x[2] + 3*x[1]*x[1]*x[2] - x[2])
end

# demanding example, newton could diverge at some point
# domain (-3,3)^2
function example_g3(x)
    para = [
        map(i -> 0.1 + (-10 + i) / 1000, StaticArrays.SUnitRange(1, 20));
        map(i -> 0.9 + (-30 + i) / 1000, StaticArrays.SUnitRange(21, 40))
    ]

    return SVector(sin(4*x[2]) * prod(x[1] .- para), sin(4*x[1]) * prod(x[2] .- para))
end

# domain (-40,40)^n, 3^n roots in domain,
function example_g4(x)
    return 100*x + x .^ 2 - x .^ 3 .- sum(x)
end

# domain (-0.3,0.8)^n
function example_g5(x)
    n = length(x)
    c = cos.(x)

    return StaticArrays.SUnitRange(1, n) .* (1 .- c) - sin.(x) .+ (n - sum(c))
end

dim = 3
g = example_g4
Dg = x -> ForwardDiff.jacobian(g, x)

center, radius = [0.0 for _ in 1:dim], [40.0 for _ in 1:dim]
P = BoxPartition(Box(center, radius))

R = cover_roots(g, Dg, P[:], depth=dim*8)

plot(R)


text:
module LL

using Tokenize

const IDENTIFIER = Tokenize.Tokens.IDENTIFIER
const WHITESPACE = Tokenize.Tokens.WHITESPACE
const STRING = Tokenize.Tokens.STRING
const INTEGER = Tokenize.Tokens.INTEGER
const FLOAT = Tokenize.Tokens.FLOAT
const DEF_ID = Symbol("def")
const DEFUN_ID = Symbol("defun")
const DEFUN_END_ID = Symbol("close")
kind(t) = Tokenize.Tokens.kind(t)
untokenize(t) = Tokenize.Tokens.untokenize(t)
null_lexpr(L) = length(L) == 0
isliteral(t) = kind(t) === STRING || kind(t) === INTEGER || kind(t) === FLOAT

parse_lexpr(source) = collect(tokenize(source))

function make_lexpr(source)
    expr = []
    for t in parse_lexpr(source)
        kind(t) === WHITESPACE && continue
        if kind(t) === IDENTIFIER
            push!(expr, Symbol(untokenize(t)))
        elseif isliteral(t)
            push!(expr, untokenize(t))
        end
    end
    expr
end

function eval_lexpr(L)
    null_lexpr(L) && return L
    atom = L[1]
    if atom === DEF_ID
        if length(L) < 3
            println("Invalid def expression at")
            println(t)
            return []
        end
        let target = L[2], value = L[3]
            @eval Primitives begin
                $(target) = $(value)
            end
        end
        L[4:end]
    elseif atom === DEFUN_ID
        close_i = findnext(x -> x == LL.DEFUN_END_ID, L, 1)
        let name = L[2], param = L[3], body = L[4:close_i-1]
            @eval Primitives begin
                function $(name)($(param))
                    @show $(body)
                    $(program)($(body), true)
                end
            end
        end
        return L[close_i+1:end]
    else
        value = Base.eval(Primitives, atom)
        if value isa String
            return [value; L[2:end]]
        else
            # @show value L[2:end]
            return value(L[2:end])
        end
    end
end

function program(L, quiet=False)
    while !null_lexpr(L)
        # @show L
        L = let new_L = eval_lexpr(L)
            # @show new_L
            if length(new_L) == length(L)
                println(new_L[1])
                new_L = new_L[2:end]
            end
            new_L
        end
    end
end

baremodule Primitives

import Base
import ..LL

function eval(L)
    LL.eval_lexpr(L)
end

function print(L)
    Base.println("My print")
    L = eval(L)
    Base.println(L[1])
    L[2:end]
end

end

end # module


text:
################################################################################
#    loggingOptimisation Problem
#
# You are the decision-maker of a production unit (here a logging company) and have a set of possible activities to undergo.
# How do you choose the activities that lead to a maximisation of the unit's profit?
# The objective of this problem is hence to find the optimal level of activities to maximise the company profit, given (a) the profitability (gross margin) of each activity, (b) the resources available to the company and (c) the matrix of technical coefficients that link each activity to the resources required (positive coefficient) or provided (negative coefficient) by that specific activity.

# The problem is the same as those in the SpreadSheet file "Optimal production mix.ods" (available in the `data` folder) and can be solved in LibreOffice by using the "Solver tool" as indicated in that file. You can use it to check your results.
#
#

# ### Environment set-up and data loading

# 1) Start by setting the working directory to the directory of this file and activate it. If you have the `Manifest.toml` file in the directory, just run `Pkg.instantiate()`, otherwise manually add the packages JuMP, GLPK, DataFrames, CSV, Pipe, HTTP.

# 2) Load the packages DelimitedFiles, JuMP, GLPK, DataFrames, CSV, Pipe, HTTP

# 3) Load from internet or from local files the following data:

urlActivities   = "https://raw.githubusercontent.com/sylvaticus/IntroSPMLJuliaCourse/main/lessonsMaterial/02_JULIA2/loggingOptimisation/data/activities.csv"
urlResources    = "https://raw.githubusercontent.com/sylvaticus/IntroSPMLJuliaCourse/main/lessonsMaterial/02_JULIA2/loggingOptimisation/data/resources.csv"
urlCoefficients = "https://raw.githubusercontent.com/sylvaticus/IntroSPMLJuliaCourse/main/lessonsMaterial/02_JULIA2/loggingOptimisation/data/coefficients.csv"

# The result must be:
# - `activities`, a dataframe with the columns `label`, `gm` and `integer`
# - `resources`, a dataframe with the columns `label`, `initial` and `initial2`
# - `coeff`, a 10x12 Matrix of Float64

# For example to download from internet and import the matrix you can use something like:
coef  = @pipe HTTP.get(urlCoefficients).body |> readdlm(_,';')

# 4) Determine `nA` and `nR` as the number of activities and the number of resources (use the `size` function)

# ### Optimisation model definition

# 5) Define `profitModel` as a model to be optimised using the `GLPK.Optimizer`
# 6) [OPTIONAL] set `GLPK.GLP_MSG_ALL` as the `msg_lev` of GLPK

# ### Model's endogenous variables definition

# 7) Define the non-negative model `x` variable, indexed by the positions between 1 and `nA` (i.e. `x[1:nA] >= 0`) (you could use the @variables macro). Don't set the `x` variable to be integer at this step, as some variables are continuous, just set them to be non-negative.

# 8) Set the variables for which the corresponding `integer` column in the `activity` dataframe is equal to 1 as a integer variable.
# To set the specific vaciable `x[i]` as integer use  `set_integer(x[i])`

# ### Model's constraint definition

# 9) Define the `resLimit[r in 1:nR]` family of contraints, such that when you sum `coef[r,a]*x[a]` for all the `1:nA` activities you must have a value not greater than `resources.initial[r]`

# ### Objective definition

# 10) Define the objective as the maximisation of the profit given by summing for each of the `1:nA` activities `activities.gm[a] * x[a]`

# ### Model resolution
# 11) [OPTIONAL] Print the model to check it
# 12) Optimize the model
# 13) Check with the function `status = termination_status(profitModel)` that the status is `OPTIMAL` (it should be!)

# ### Print optimal level of activities

# 14) Run the following code to print the results (optimal level of activities)

if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED || status == MOI.TIME_LIMIT) && has_values(profitModel)
    println("#################################################################")
    if (status == MOI.OPTIMAL)
        println("** Problem solved correctly **")
    else
        println("** Problem returned a (possibly suboptimal) solution **")
    end
    println("- Objective value (total costs): ", objective_value(profitModel))
    println("- Optimal Activities:\n")
    optValues = value.(x)
    for a in 1:nA
      println("* $(activities.label[a]):\t $(optValues[a])")
    end
    if JuMP.has_duals(profitModel)
        println("\n\n- Shadow prices of the resources:\n")
        for r in 1:nR
            println("* $(resources.label[r]):\t $(dual(resLimit[r]))")
        end
    end
else
    println("The model was not solved correctly.")
    println(status)
end

# 15) [OPTIONAL] Observe the emergence of scale effects
# Optionally re-run the model with `initial2` initial resources (instead of `initial`) and notice how this larger company can afford to perform different types of activities (logging high forest instead of coppices in this example) and obtain a better profitability per unit of resource employed.
# Can you guess which are the aspects of the optimisation model that allow for the emergence of these scale effects?

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advectionvelocity = (1.0, 1.0)
# advectionvelocity = (0.2, -0.3)
equations = LinearScalarAdvectionEquation2D(advectionvelocity)

initial_condition = initial_condition_convergence_test

# you can either use a single function to impose the BCs weakly in all
# 1*ndims == 2 directions or you can pass a tuple containing BCs for
# each direction
# Note: "boundary_condition_periodic" indicates that it is a periodic boundary and can be omitted on
#       fully periodic domains. Here, however, it is included to allow easy override during testing
boundary_conditions = boundary_condition_periodic

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

# The initial condition is 2-periodic
coordinates_min = (-1.5, 1.3) # minimum coordinates (min(x), min(y))
coordinates_max = ( 0.5, 5.3) # maximum coordinates (max(x), max(y))

cells_per_dimension = (19, 37)

# Create curved mesh with 19 x 37 elements
mesh = CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max)

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.0
tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan);

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(entropy, energy_total))

# The AliveCallback prints short status information in regular intervals
alive_callback = AliveCallback(analysis_interval=analysis_interval)

# The SaveRestartCallback allows to save a file from which a Trixi simulation can be restarted
save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
cfl=1.6
stepsize_callback = StepsizeCallback(cfl=cfl)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_restart, save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks, maxiters=1e5);

# Print the timer summary
summary_callback()


text:
function solve()
    n = parse(Int, readline())
    ret = 0
    for i=1:n
        a, b = [parse(Int, x) for x in split(readline())]
        ret += (b-a+1)*(a+b)÷2
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text: const cosmic_fl = "../data/CosmicCompleteExport.tsv"


text:
function triangulate_cell(cell_connectivity, cell_type::Int)
    if cell_type == 5 || cell_type == 6 || cell_type == 10 || cell_type == 24
        return NTuple{3, Int}.(decompose_cell(cell_connectivity, cell_type, target = "Faces")[1])
    elseif cell_type == 22
        return NTuple{3, Int}[(cell_connectivity[1], cell_connectivity[2], cell_connectivity[3])]
    elseif cell_type == 8
        return NTuple{3, Int}[(cell_connectivity[1], cell_connectivity[2], cell_connectivity[4]),
            (cell_connectivity[1], cell_connectivity[4], cell_connectivity[3])]
    elseif cell_type == 9 || cell_type == 24
        return NTuple{3, Int}[(cell_connectivity[1], cell_connectivity[2], cell_connectivity[3]),
            (cell_connectivity[1], cell_connectivity[3], cell_connectivity[4])]
    elseif cell_type == 11 || cell_type == 12 || cell_type == 13 || cell_type == 14
        tris = NTuple{3, Int}[]
        for (_cell, _cell_type) in zip(decompose_cell(cell_connectivity, cell_type, target = "Faces")...)
            append!(tris, triangulate_cell(_cell, _cell_type))
        end
        return tris
    else
        throw("Cell type $cell_type is not supported.")
    end
end

function triangulate_cell_glmesh(cell_connectivity, cell_type::Int)
    if cell_type == 5 || cell_type == 6 || cell_type == 10
        return map(GLFace{3}, decompose_cell(cell_connectivity, cell_type, target = "Faces")[1] .- 1)
    elseif cell_type == 8
        return GLFace{3}[GLFace{3}(cell_connectivity[1], cell_connectivity[2], cell_connectivity[4]), GLFace{3}(cell_connectivity[1], cell_connectivity[4], cell_connectivity[3])]
    elseif cell_type == 9
        return GLFace{3}[GLFace{3}(cell_connectivity[1], cell_connectivity[2], cell_connectivity[3]), GLFace{3}(cell_connectivity[1], cell_connectivity[3], cell_connectivity[4])]
    elseif cell_type == 11 || cell_type == 12 || cell_type == 13 || cell_type == 14
        tris = GLFace{3}[]
        for (_cell, _cell_type) in zip(decompose_cell(cell_connectivity, cell_type, target = "Faces")...)
            append!(tris, triangulate_cell_glmesh(_cell, _cell_type))
        end
        return tris
    else
        throw("Cell type $cell_type is not supported.")
    end
end


text:
N = 4
E = 5
adj = [0 1 1 1;
       1 0 1 0;
       1 1 0 1;
       1 0 1 0]
adj2 = [0 1 0 1;
        1 0 1 1;
        0 1 0 1;
        1 1 1 0]

ug = SimpleGraph(N)
add_edge!(ug, 1, 2); add_edge!(ug, 1, 3); add_edge!(ug, 1, 4)
add_edge!(ug, 2, 3); add_edge!(ug, 3, 4)

nf = rand(3, N)
ef = rand(5, E)
gf = rand(7)

@testset "featuredgraph" begin
    ng = NullGraph()
    @test FeaturedGraph() == ng
    @test has_graph(ng) == false
    @test has_node_feature(ng) == false
    @test has_edge_feature(ng) == false
    @test has_global_feature(ng) == false
    @test isnothing(graph(ng))
    @test isnothing(node_feature(ng))
    @test isnothing(edge_feature(ng))
    @test isnothing(global_feature(ng))
    @test isnothing(mask(ng))

    fg = FeaturedGraph(adj)
    @test has_graph(fg)
    @test has_node_feature(fg) == false
    @test has_edge_feature(fg) == false
    @test has_global_feature(fg) == false
    @test graph(fg) == adj
    @test node_feature(fg) == Fill(0, (0, N))
    @test edge_feature(fg) == Fill(0, (0, E))
    @test global_feature(fg) == Fill(0, 0)
    @test mask(fg) == zeros(N, N)

    fg = FeaturedGraph(adj, nf)
    @test has_graph(fg)
    @test has_node_feature(fg)
    @test has_edge_feature(fg) == false
    @test has_global_feature(fg) == false
    @test graph(fg) == adj
    @test node_feature(fg) == nf
    @test edge_feature(fg) == Fill(0., (0, E))
    @test global_feature(fg) == zeros(0)
    @test mask(fg) == zeros(N, N)

    fg = FeaturedGraph(ug, nf)
    @test has_graph(fg)
    @test has_node_feature(fg)
    @test has_edge_feature(fg) == false
    @test has_global_feature(fg) == false
    @test graph(fg) == ug
    @test node_feature(fg) == nf
    @test edge_feature(fg) == Fill(0., (0, E))
    @test global_feature(fg) == zeros(0)
    @test mask(fg) == zeros(N, N)

    fg = FeaturedGraph(ug, nf, ef ,gf)
    @test has_graph(fg)
    @test has_node_feature(fg)
    @test has_edge_feature(fg)
    @test has_global_feature(fg)
    @test graph(fg) == ug
    @test node_feature(fg) == nf
    @test edge_feature(fg) == ef
    @test global_feature(fg) == gf
    @test mask(fg) == zeros(N, N)

    fg = FeaturedGraph(adj, nf, ef ,gf)
    @test has_graph(fg)
    @test has_node_feature(fg)
    @test has_edge_feature(fg)
    @test has_global_feature(fg)
    @test graph(fg) == adj
    @test node_feature(fg) == nf
    @test edge_feature(fg) == ef
    @test global_feature(fg) == gf
    @test mask(fg) == zeros(N, N)

    T = Matrix{Float32}
    fg = FeaturedGraph{T,T,T,Vector{Float32}}(adj, nf, ef ,gf, zeros(N, N), :adjm, true)
    @test typeof(graph(fg)) == T
    @test typeof(node_feature(fg)) == T
    @test typeof(edge_feature(fg)) == T
    @test typeof(global_feature(fg)) == Vector{Float32}

    # fg = FeaturedGraph(adj, nf, ef ,gf)
    # fg.graph = adj2
    # @test fg.graph == adj2
    # @test_throws DimensionMismatch fg.graph = [0 1; 0 1]

    # nf_ = rand(10, N)
    # fg.nf = nf_
    # @test fg.nf == nf_
    # @test_throws DimensionMismatch fg.nf = rand(10, 11)

    # ef_ = rand(10, E)
    # fg.ef = ef_
    # @test fg.ef == ef_
    # @test_throws DimensionMismatch fg.ef = rand(10, 11)
end


text:
abstract type SpectrumFilter end

# From: https://gist.github.com/jiahao/b8b5ac328c18b7ae8a17

struct SavitzkyGolayFilter{M,N} <: SpectrumFilter end
@generated function (::SavitzkyGolayFilter{M,N})(data::AbstractVector{T}) where {M,N,T}
    # Create Jacobian matrix
    J = zeros(2M + 1, N + 1)
    for i = 1:2M+1, j = 1:N+1
        J[i, j] = (i - M - 1)^(j - 1)
    end
    e₁ = zeros(N + 1)
    e₁[1] = 1.0

    # Compute filter coefficients
    C = J' \ e₁
    # Evaluate filter on data matrix

    To = typeof(C[1] * one(T)) # Calculate type of output
    expr = quote
        n = size(data, 1)
        smoothed = zeros($To, n)
        @inbounds for i in eachindex(smoothed)
            smoothed[i] += $(C[M+1]) * data[i]
        end
        smoothed
    end

    for j = 1:M
        insert!(
            expr.args[6].args[3].args[2].args,
            1,
            :(
                if i - $j ≥ 1
                    smoothed[i] += $(C[M+1-j]) * data[i-$j]
                end
            ),
        )
        push!(
            expr.args[6].args[3].args[2].args,
            :(
                if i + $j ≤ n
                    smoothed[i] += $(C[M+1+j]) * data[i+$j]
                end
            ),
        )
    end
    return expr
end

"""
    apply(filt::SavitzkyGolayFilter, spec::Spectrum, applyLLD=false)

Applys a function to the channel data in `spec` (with/wo the low-level discriminator.)
The function can only be a function of the counts data and can not change the
energy scale or spectrum properties.  The result is a Spectrum.

Example:

    apply(SavitzkyGolayFilter{6,4}(),spec)
"""
function apply(filt::SavitzkyGolayFilter, spec::Spectrum, applyLLD = false)
    res = Spectrum(
        spec.energy,
        filt(counts(spec, eltype(spec), applyLLD)),
        copy(spec.properties),
    )
    res[:Name] = repr(filt)[1:end-1] * "$(spec[:Name]))"
    res[:Parent] = spec
    res
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
using SDWBA, Docile, Lexicon

Lexicon.save("API.md", SDWBA)
Lexicon.save("Models.md", SDWBA.Models)
cd("../")
run(`mkdocs build --clean`)
cd("docs")


text:
function geometrixmatrix(l::Float64)
  gm = zeros(12, 12)
  pt1 = 0.1
  opt2 = 1.2
  two = 2.0
  d15 = 15.0
  d30 = 30.0

  gm[1, 1]   = opt2 / l
  gm[1, 5]   = pt1
  gm[1, 7]   = -opt2 / l
  gm[1, 11]  = pt1

  gm[2, 2]   = opt2 / l
  gm[2, 4]   = -pt1
  gm[2, 8]   = -opt2 / l
  gm[2, 10]  = -pt1

  gm[4, 2]   = -pt1
  gm[4, 4]   = two * l / d15
  gm[4, 8]   = pt1
  gm[4, 10]  = -l / d30

  gm[5, 1]   = pt1
  gm[5, 5]   = two * l / d15
  gm[5, 7]   = -pt1
  gm[5, 11]  = -l / d30

  gm[7, 1]   = -opt2 / l
  gm[7, 5]   = -pt1
  gm[7, 7]   = opt2 / l
  gm[7, 11]  = -pt1

  gm[8, 2]   = -opt2 / l
  gm[8, 4]   = pt1
  gm[8, 8]   = opt2 / l
  gm[8, 10]  = pt1

  gm[10, 2]  = -pt1
  gm[10, 4]  = -l / d30
  gm[10, 8]  = pt1
  gm[10, 10] = two * l / d15

  gm[11, 1]  = pt1
  gm[11, 5]  = -l / d30
  gm[11, 7]  = -pt1
  gm[11, 11] = two * l / d15

  gm
end


text:
# # [Basic Epidemiology Models](@id epidemiology_example)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/covid/epidemiology.ipynb)

using AlgebraicPetri
using AlgebraicPetri.Epidemiology

using LabelledArrays
using OrdinaryDiffEq
using Plots

using Catlab
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra
using Catlab.Programs.RelationalPrograms

display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>".75"));

# #### SIR Model:

# define model
sir = @relation (s,i,r) begin
    infection(s,i)
    recovery(i,r)
end
display_uwd(sir)
#-
p_sir = apex(oapply_epi(sir))
Graph(p_sir)

# define initial states and transition rates, then
# create, solve, and visualize ODE problem

u0 = LVector(S=10, I=1, R=0);
p = LVector(inf=0.4, rec=0.4);

# The C-Set representation has direct support for generating a DiffEq vector field

prob = ODEProblem(vectorfield(p_sir),u0,(0.0,7.5),p);
sol = solve(prob,Tsit5())

plot(sol)

# #### SEIR Model:

# define model
seir = @relation (s,e,i,r) begin
    exposure(s,i,e)
    illness(e,i)
    recovery(i,r)
end
display_uwd(seir)
#-
p_seir = apex(oapply_epi(seir))
Graph(p_seir)

# define initial states and transition rates, then
# create, solve, and visualize ODE problem

u0 = LVector(S=10, E=1, I=0, R=0);
p = LVector(exp=.9, ill=.2, rec=.5);

prob = ODEProblem(vectorfield(p_seir),u0,(0.0,15.0),p);
sol = solve(prob,Tsit5())

plot(sol)

# #### SEIRD Model:

# define model
seird = @relation (s,e,i,r,d) begin
    exposure(s,i,e)
    illness(e,i)
    recovery(i,r)
    death(i,d)
end
display_uwd(seird)
#-
p_seird = apex(oapply_epi(seird))
Graph(p_seird)

# define initial states and transition rates, then
# create, solve, and visualize ODE problem

u0 = LVector(S=10, E=1, I=0, R=0, D=0);
p = LVector(exp=0.9, ill=0.2, rec=0.5, death=0.1);

prob = ODEProblem(vectorfield(p_seird),u0,(0.0,15.0),p);
sol = solve(prob,Tsit5())

plot(sol)

max_stars_repo_name: AlgebraicJulia/AlgebraicPetri

max_issues_repo_name: AlgebraicJulia/AlgebraicPetri

max_forks_repo_name: AlgebraicJulia/AlgebraicPetri


text:
using JuMP
using MathOptInterface
using GLPK
using Test
using MathOptSetDistances
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOD = MathOptSetDistances

variable_primal(model) = x -> MOI.get(model, MOI.VariablePrimal(), x)

model = Model(GLPK.Optimizer)
@variable(model, x>=1)
@variable(model, y>=0)
@constraint(model,c1, x + y >= 3)
@constraint(model,c2, x + y >= 2)
@objective(model, Min, x + y)
optimize!(model)

varval = variable_primal(backend(model))
f = JuMP.MOI.get(model, JuMP.MOI.ConstraintFunction(), c2)

delete(model,c2)
val  = MOIU.eval_variables(varval, f)

# m = HeuristicOptimizer.heuristicCreate(model)
# HeuristicOptimizer.removeConstrains!(m)
# optimize!(m.model)
# HeuristicOptimizer.updateFunctionValue!(m)
# HeuristicOptimizer.updateViolation!(m)
# @show m.violation

c = f, s


text:
# define an abstrat tree node type - concrete types are TreeSplit and TreeLeaf
abstract type Node{T<:AbstractFloat} end

abstract type ModelType end
abstract type GradientRegression <: ModelType end
abstract type L1Regression <: ModelType end
abstract type QuantileRegression <: ModelType end
abstract type MultiClassRegression <: ModelType end
abstract type GaussianRegression <: ModelType end
struct Linear <: GradientRegression end
struct Poisson <: GradientRegression end
struct Logistic <: GradientRegression end
struct L1 <: L1Regression end
struct Quantile <: QuantileRegression end
struct Softmax <: MultiClassRegression end
struct Gaussian <: GaussianRegression end

# store perf info of each variable
mutable struct SplitInfo{L, T<:AbstractFloat, S<:Int}
    gain::T
    ∑δL::SVector{L,T}
    ∑δ²L::SVector{L,T}
    ∑𝑤L::SVector{1,T}
    ∑δR::SVector{L,T}
    ∑δ²R::SVector{L,T}
    ∑𝑤R::SVector{1,T}
    gainL::T
    gainR::T
    𝑖::S
    feat::S
    cond::T
end

struct TreeNode{L, T<:AbstractFloat, S<:Int, B<:Bool}
    left::S
    right::S
    feat::S
    cond::T
    pred::SVector{L,T}
    split::B
end

TreeNode(left::S, right::S, feat::S, cond::T, L::S) where {T<:AbstractFloat, S<:Int} = TreeNode{L,T,S,Bool}(left, right, feat, cond, zeros(SVector{L,T}), true)
TreeNode(pred::SVector{L,T}) where {L,T} = TreeNode(0, 0, 0, 0.0, pred, false)

mutable struct EvoTreeRegressor{T<:AbstractFloat, U<:ModelType, S<:Int} #<: MLJBase.Deterministic
    loss::U
    nrounds::S
    λ::T
    γ::T
    η::T
    max_depth::S
    min_weight::T # real minimum number of observations, different from xgboost (but same for linear)
    rowsample::T # subsample
    colsample::T
    nbins::S
    α::T
    metric::Symbol
    seed::S
    K::S # length of predictions: 1 by default, > 1 for multiclassif or maxloglikelihood
end

function EvoTreeRegressor(;
    loss=:linear,
    nrounds=10,
    λ=0.0, #
    γ=0.0, # gamma: min gain to split
    η=0.1, # eta: learning rate
    max_depth=5,
    min_weight=1.0, # minimal weight, different from xgboost (but same for linear)
    rowsample=1.0,
    colsample=1.0,
    nbins=64,
    α=0.5,
    metric=:mse,
    seed=444,
    K=1)

    if loss == :linear model_type = Linear()
    elseif loss == :logistic model_type = Logistic()
    elseif loss == :poisson model_type = Poisson()
    elseif loss == :L1 model_type = L1()
    elseif loss == :quantile model_type = Quantile()
    elseif loss == :softmax model_type = Softmax()
    elseif loss == :gaussian model_type = Gaussian()
    end

    # override K for gaussian
    if loss == :gaussian
        K = 2
    end

    model = EvoTreeRegressor(model_type, nrounds, λ, γ, η, max_depth, min_weight, rowsample, colsample, nbins, α, metric, seed, K)
    # message = MLJBase.clean!(model)
    # isempty(message) || @warn message
    return model
end

# For R-package
function EvoTreeRegressorR(
    loss,
    nrounds,
    λ,
    γ,
    η,
    max_depth,
    min_weight,
    rowsample,
    colsample,
    nbins,
    α,
    metric,
    seed,
    K)

    if loss == :linear model_type = Linear()
    elseif loss == :logistic model_type = Logistic()
    elseif loss == :poisson model_type = Poisson()
    elseif loss == :L1 model_type = L1()
    elseif loss == :quantile model_type = Quantile()
    elseif loss == :softmax model_type = Softmax()
    elseif loss == :gaussian model_type = Gaussian()
    end

    if loss == :gaussian
        K = 2
    end

    model = EvoTreeRegressor(model_type, nrounds, λ, γ, η, max_depth, min_weight, rowsample, colsample, nbins, α, metric, seed, K)
    # message = MLJBase.clean!(model)
    # isempty(message) || @warn message
    return model
end

# single tree is made of a root node that containes nested nodes and leafs
struct TrainNode{L, T<:AbstractFloat, S<:Int}
    depth::S
    ∑δ::SVector{L,T}
    ∑δ²::SVector{L,T}
    ∑𝑤::SVector{1,T}
    gain::T
    𝑖::BitSet
    𝑗::Vector{S}
end

# single tree is made of a root node that containes nested nodes and leafs
struct Tree{L, T<:AbstractFloat, S<:Int}
    nodes::Vector{TreeNode{L,T,S,Bool}}
end

# eval metric tracking
struct Metric
    iter::Vector{Int}
    metric::Vector{Float64}
end
Metric() = Metric([0], [Inf])

# gradient-boosted tree is formed by a vector of trees
struct GBTree{L, T<:AbstractFloat, S<:Int}
    trees::Vector{Tree{L,T,S}}
    params::EvoTreeRegressor
    metric::Metric
end


text:
using PGFPlots
using Colors
using ColorSchemes
using LinearAlgebra
pasteljet = ColorMaps.RGBArrayMap(ColorSchemes.viridis, interpolation_levels=500, invert=false);
pastelRedBlue = ColorMaps.RGBArrayMap([RGB(246/255, 21/255, 92/255),
                                       RGB(1.0,1.0,1.0),
                                       RGB( 27/255,161/255,234/255)], interpolation_levels=500);
pushPGFPlotsPreamble("\\pgfplotscreateplotcyclelist{pastelcolors}{%")
pushPGFPlotsPreamble("  solid, pastelPurple, mark=none\\\\%")
pushPGFPlotsPreamble("  solid, pastelBlue, mark=none\\\\%")
pushPGFPlotsPreamble("  solid, pastelGreen, mark=none\\\\%")
pushPGFPlotsPreamble("  solid, pastelRed, mark=none\\\\%")
pushPGFPlotsPreamble("  solid, pastelMagenta, mark=none\\\\%")
pushPGFPlotsPreamble("  solid, pastelOrange, mark=none\\\\%")
pushPGFPlotsPreamble("  solid, pastelSeaGreen, mark=none\\\\%")
pushPGFPlotsPreamble("}")
pushPGFPlotsPreamble("\\usepackage{amsmath}")
pushPGFPlotsPreamble("\\usepgfplotslibrary{fillbetween}")
pushPGFPlotsPreamble("\\usepgfplotslibrary{statistics}")
pushPGFPlotsPreamble("\\definecolor{pastelMagenta}{HTML}{FF48CF}")
pushPGFPlotsPreamble("\\definecolor{pastelPurple}{HTML}{8770FE}")
pushPGFPlotsPreamble("\\definecolor{pastelBlue}{HTML}{1BA1EA}")
pushPGFPlotsPreamble("\\definecolor{pastelSeaGreen}{HTML}{14B57F}")
pushPGFPlotsPreamble("\\definecolor{pastelGreen}{HTML}{3EAA0D}")
pushPGFPlotsPreamble("\\definecolor{pastelOrange}{HTML}{C38D09}")
pushPGFPlotsPreamble("\\definecolor{pastelRed}{HTML}{F5615C}")
pushPGFPlotsPreamble("\\tikzset{myarrow/.style={line width = 0.05cm, ->, rounded corners=5mm}}")
pushPGFPlotsPreamble("\\newcommand{\\vect}[1]{\\boldsymbol{\\mathbf{#1}}}")
pushPGFPlotsPreamble("\\pgfplotsset{")
pushPGFPlotsPreamble("   	colormap={pasteljet}{")
pushPGFPlotsPreamble("		rgb=(0.26700,0.00487,0.32942)")
pushPGFPlotsPreamble("		rgb=(0.26851,0.00961,0.33543)")
pushPGFPlotsPreamble("		rgb=(0.26994,0.01463,0.34138)")
pushPGFPlotsPreamble("		rgb=(0.27131,0.01994,0.34727)")
pushPGFPlotsPreamble("		rgb=(0.27259,0.02556,0.35309)")
pushPGFPlotsPreamble("		rgb=(0.27381,0.03150,0.35885)")
pushPGFPlotsPreamble("		rgb=(0.27495,0.03775,0.36454)")
pushPGFPlotsPreamble("		rgb=(0.27602,0.04417,0.37016)")
pushPGFPlotsPreamble("		rgb=(0.27702,0.05034,0.37572)")
pushPGFPlotsPreamble("		rgb=(0.27794,0.05632,0.38119)")
pushPGFPlotsPreamble("		rgb=(0.27879,0.06214,0.38659)")
pushPGFPlotsPreamble("		rgb=(0.27957,0.06784,0.39192)")
pushPGFPlotsPreamble("		rgb=(0.28027,0.07342,0.39716)")
pushPGFPlotsPreamble("		rgb=(0.28089,0.07891,0.40233)")
pushPGFPlotsPreamble("		rgb=(0.28145,0.08432,0.40741)")
pushPGFPlotsPreamble("		rgb=(0.28192,0.08967,0.41241)")
pushPGFPlotsPreamble("		rgb=(0.28233,0.09495,0.41733)")
pushPGFPlotsPreamble("		rgb=(0.28266,0.10020,0.42216)")
pushPGFPlotsPreamble("		rgb=(0.28291,0.10539,0.42690)")
pushPGFPlotsPreamble("		rgb=(0.28309,0.11055,0.43155)")
pushPGFPlotsPreamble("		rgb=(0.28320,0.11568,0.43611)")
pushPGFPlotsPreamble("		rgb=(0.28323,0.12078,0.44058)")
pushPGFPlotsPreamble("		rgb=(0.28319,0.12585,0.44496)")
pushPGFPlotsPreamble("		rgb=(0.28307,0.13090,0.44924)")
pushPGFPlotsPreamble("		rgb=(0.28288,0.13592,0.45343)")
pushPGFPlotsPreamble("		rgb=(0.28262,0.14093,0.45752)")
pushPGFPlotsPreamble("		rgb=(0.28229,0.14591,0.46151)")
pushPGFPlotsPreamble("		rgb=(0.28189,0.15088,0.46541)")
pushPGFPlotsPreamble("		rgb=(0.28141,0.15583,0.46920)")
pushPGFPlotsPreamble("		rgb=(0.28087,0.16077,0.47290)")
pushPGFPlotsPreamble("		rgb=(0.28025,0.16569,0.47650)")
pushPGFPlotsPreamble("		rgb=(0.27957,0.17060,0.48000)")
pushPGFPlotsPreamble("		rgb=(0.27883,0.17549,0.48340)")
pushPGFPlotsPreamble("		rgb=(0.27801,0.18037,0.48670)")
pushPGFPlotsPreamble("		rgb=(0.27713,0.18523,0.48990)")
pushPGFPlotsPreamble("		rgb=(0.27619,0.19007,0.49300)")
pushPGFPlotsPreamble("		rgb=(0.27519,0.19490,0.49600)")
pushPGFPlotsPreamble("		rgb=(0.27413,0.19972,0.49891)")
pushPGFPlotsPreamble("		rgb=(0.27301,0.20452,0.50172)")
pushPGFPlotsPreamble("		rgb=(0.27183,0.20930,0.50443)")
pushPGFPlotsPreamble("		rgb=(0.27059,0.21407,0.50705)")
pushPGFPlotsPreamble("		rgb=(0.26931,0.21882,0.50958)")
pushPGFPlotsPreamble("		rgb=(0.26797,0.22355,0.51201)")
pushPGFPlotsPreamble("		rgb=(0.26658,0.22826,0.51435)")
pushPGFPlotsPreamble("		rgb=(0.26515,0.23296,0.51660)")
pushPGFPlotsPreamble("		rgb=(0.26366,0.23763,0.51876)")
pushPGFPlotsPreamble("		rgb=(0.26214,0.24229,0.52084)")
pushPGFPlotsPreamble("		rgb=(0.26057,0.24692,0.52283)")
pushPGFPlotsPreamble("		rgb=(0.25897,0.25154,0.52474)")
pushPGFPlotsPreamble("		rgb=(0.25732,0.25613,0.52656)")
pushPGFPlotsPreamble("		rgb=(0.25565,0.26070,0.52831)")
pushPGFPlotsPreamble("		rgb=(0.25394,0.26525,0.52998)")
pushPGFPlotsPreamble("		rgb=(0.25219,0.26978,0.53158)")
pushPGFPlotsPreamble("		rgb=(0.25043,0.27429,0.53310)")
pushPGFPlotsPreamble("		rgb=(0.24863,0.27877,0.53456)")
pushPGFPlotsPreamble("		rgb=(0.24681,0.28324,0.53594)")
pushPGFPlotsPreamble("		rgb=(0.24497,0.28768,0.53726)")
pushPGFPlotsPreamble("		rgb=(0.24311,0.29209,0.53852)")
pushPGFPlotsPreamble("		rgb=(0.24124,0.29648,0.53971)")
pushPGFPlotsPreamble("		rgb=(0.23935,0.30085,0.54084)")
pushPGFPlotsPreamble("		rgb=(0.23744,0.30520,0.54192)")
pushPGFPlotsPreamble("		rgb=(0.23553,0.30953,0.54294)")
pushPGFPlotsPreamble("		rgb=(0.23360,0.31383,0.54391)")
pushPGFPlotsPreamble("		rgb=(0.23167,0.31811,0.54483)")
pushPGFPlotsPreamble("		rgb=(0.22974,0.32236,0.54571)")
pushPGFPlotsPreamble("		rgb=(0.22780,0.32659,0.54653)")
pushPGFPlotsPreamble("		rgb=(0.22586,0.33081,0.54731)")
pushPGFPlotsPreamble("		rgb=(0.22393,0.33499,0.54805)")
pushPGFPlotsPreamble("		rgb=(0.22199,0.33916,0.54875)")
pushPGFPlotsPreamble("		rgb=(0.22006,0.34331,0.54941)")
pushPGFPlotsPreamble("		rgb=(0.21813,0.34743,0.55004)")
pushPGFPlotsPreamble("		rgb=(0.21621,0.35153,0.55063)")
pushPGFPlotsPreamble("		rgb=(0.21430,0.35562,0.55118)")
pushPGFPlotsPreamble("		rgb=(0.21240,0.35968,0.55171)")
pushPGFPlotsPreamble("		rgb=(0.21050,0.36373,0.55221)")
pushPGFPlotsPreamble("		rgb=(0.20862,0.36775,0.55268)")
pushPGFPlotsPreamble("		rgb=(0.20676,0.37176,0.55312)")
pushPGFPlotsPreamble("		rgb=(0.20490,0.37575,0.55353)")
pushPGFPlotsPreamble("		rgb=(0.20306,0.37972,0.55393)")
pushPGFPlotsPreamble("		rgb=(0.20124,0.38367,0.55429)")
pushPGFPlotsPreamble("		rgb=(0.19943,0.38761,0.55464)")
pushPGFPlotsPreamble("		rgb=(0.19764,0.39153,0.55497)")
pushPGFPlotsPreamble("		rgb=(0.19586,0.39543,0.55528)")
pushPGFPlotsPreamble("		rgb=(0.19410,0.39932,0.55556)")
pushPGFPlotsPreamble("		rgb=(0.19236,0.40320,0.55584)")
pushPGFPlotsPreamble("		rgb=(0.19063,0.40706,0.55609)")
pushPGFPlotsPreamble("		rgb=(0.18892,0.41091,0.55633)")
pushPGFPlotsPreamble("		rgb=(0.18723,0.41475,0.55655)")
pushPGFPlotsPreamble("		rgb=(0.18556,0.41857,0.55675)")
pushPGFPlotsPreamble("		rgb=(0.18390,0.42238,0.55694)")
pushPGFPlotsPreamble("		rgb=(0.18226,0.42618,0.55712)")
pushPGFPlotsPreamble("		rgb=(0.18063,0.42997,0.55728)")
pushPGFPlotsPreamble("		rgb=(0.17902,0.43376,0.55743)")
pushPGFPlotsPreamble("		rgb=(0.17742,0.43753,0.55756)")
pushPGFPlotsPreamble("		rgb=(0.17584,0.44129,0.55768)")
pushPGFPlotsPreamble("		rgb=(0.17427,0.44504,0.55779)")
pushPGFPlotsPreamble("		rgb=(0.17272,0.44879,0.55788)")
pushPGFPlotsPreamble("		rgb=(0.17118,0.45253,0.55797)")
pushPGFPlotsPreamble("		rgb=(0.16965,0.45626,0.55803)")
pushPGFPlotsPreamble("		rgb=(0.16813,0.45999,0.55808)")
pushPGFPlotsPreamble("		rgb=(0.16662,0.46371,0.55812)")
pushPGFPlotsPreamble("		rgb=(0.16512,0.46742,0.55814)")
pushPGFPlotsPreamble("		rgb=(0.16362,0.47113,0.55815)")
pushPGFPlotsPreamble("		rgb=(0.16214,0.47484,0.55814)")
pushPGFPlotsPreamble("		rgb=(0.16067,0.47854,0.55812)")
pushPGFPlotsPreamble("		rgb=(0.15919,0.48224,0.55807)")
pushPGFPlotsPreamble("		rgb=(0.15773,0.48593,0.55801)")
pushPGFPlotsPreamble("		rgb=(0.15627,0.48962,0.55794)")
pushPGFPlotsPreamble("		rgb=(0.15482,0.49331,0.55784)")
pushPGFPlotsPreamble("		rgb=(0.15336,0.49700,0.55772)")
pushPGFPlotsPreamble("		rgb=(0.15192,0.50069,0.55759)")
pushPGFPlotsPreamble("		rgb=(0.15048,0.50437,0.55743)")
pushPGFPlotsPreamble("		rgb=(0.14904,0.50805,0.55725)")
pushPGFPlotsPreamble("		rgb=(0.14761,0.51173,0.55705)")
pushPGFPlotsPreamble("		rgb=(0.14618,0.51541,0.55682)")
pushPGFPlotsPreamble("		rgb=(0.14476,0.51909,0.55657)")
pushPGFPlotsPreamble("		rgb=(0.14334,0.52277,0.55629)")
pushPGFPlotsPreamble("		rgb=(0.14194,0.52645,0.55599)")
pushPGFPlotsPreamble("		rgb=(0.14054,0.53013,0.55566)")
pushPGFPlotsPreamble("		rgb=(0.13915,0.53381,0.55530)")
pushPGFPlotsPreamble("		rgb=(0.13777,0.53749,0.55491)")
pushPGFPlotsPreamble("		rgb=(0.13641,0.54117,0.55448)")
pushPGFPlotsPreamble("		rgb=(0.13507,0.54485,0.55403)")
pushPGFPlotsPreamble("		rgb=(0.13374,0.54853,0.55354)")
pushPGFPlotsPreamble("		rgb=(0.13244,0.55222,0.55302)")
pushPGFPlotsPreamble("		rgb=(0.13117,0.55590,0.55246)")
pushPGFPlotsPreamble("		rgb=(0.12993,0.55958,0.55186)")
pushPGFPlotsPreamble("		rgb=(0.12873,0.56327,0.55123)")
pushPGFPlotsPreamble("		rgb=(0.12757,0.56695,0.55056)")
pushPGFPlotsPreamble("		rgb=(0.12645,0.57063,0.54984)")
pushPGFPlotsPreamble("		rgb=(0.12539,0.57432,0.54909)")
pushPGFPlotsPreamble("		rgb=(0.12440,0.57800,0.54829)")
pushPGFPlotsPreamble("		rgb=(0.12346,0.58169,0.54744)")
pushPGFPlotsPreamble("		rgb=(0.12261,0.58537,0.54656)")
pushPGFPlotsPreamble("		rgb=(0.12183,0.58905,0.54562)")
pushPGFPlotsPreamble("		rgb=(0.12115,0.59274,0.54464)")
pushPGFPlotsPreamble("		rgb=(0.12057,0.59642,0.54361)")
pushPGFPlotsPreamble("		rgb=(0.12009,0.60010,0.54253)")
pushPGFPlotsPreamble("		rgb=(0.11974,0.60379,0.54140)")
pushPGFPlotsPreamble("		rgb=(0.11951,0.60746,0.54022)")
pushPGFPlotsPreamble("		rgb=(0.11942,0.61114,0.53898)")
pushPGFPlotsPreamble("		rgb=(0.11948,0.61482,0.53769)")
pushPGFPlotsPreamble("		rgb=(0.11970,0.61849,0.53635)")
pushPGFPlotsPreamble("		rgb=(0.12008,0.62216,0.53495)")
pushPGFPlotsPreamble("		rgb=(0.12064,0.62583,0.53349)")
pushPGFPlotsPreamble("		rgb=(0.12138,0.62949,0.53197)")
pushPGFPlotsPreamble("		rgb=(0.12231,0.63315,0.53040)")
pushPGFPlotsPreamble("		rgb=(0.12344,0.63681,0.52876)")
pushPGFPlotsPreamble("		rgb=(0.12478,0.64046,0.52707)")
pushPGFPlotsPreamble("		rgb=(0.12633,0.64411,0.52531)")
pushPGFPlotsPreamble("		rgb=(0.12809,0.64775,0.52349)")
pushPGFPlotsPreamble("		rgb=(0.13007,0.65138,0.52161)")
pushPGFPlotsPreamble("		rgb=(0.13227,0.65501,0.51966)")
pushPGFPlotsPreamble("		rgb=(0.13469,0.65864,0.51765)")
pushPGFPlotsPreamble("		rgb=(0.13734,0.66225,0.51557)")
pushPGFPlotsPreamble("		rgb=(0.14021,0.66586,0.51343)")
pushPGFPlotsPreamble("		rgb=(0.14330,0.66946,0.51121)")
pushPGFPlotsPreamble("		rgb=(0.14662,0.67305,0.50894)")
pushPGFPlotsPreamble("		rgb=(0.15015,0.67663,0.50659)")
pushPGFPlotsPreamble("		rgb=(0.15389,0.68020,0.50417)")
pushPGFPlotsPreamble("		rgb=(0.15785,0.68376,0.50169)")
pushPGFPlotsPreamble("		rgb=(0.16202,0.68732,0.49913)")
pushPGFPlotsPreamble("		rgb=(0.16638,0.69086,0.49650)")
pushPGFPlotsPreamble("		rgb=(0.17095,0.69438,0.49380)")
pushPGFPlotsPreamble("		rgb=(0.17571,0.69790,0.49103)")
pushPGFPlotsPreamble("		rgb=(0.18065,0.70140,0.48819)")
pushPGFPlotsPreamble("		rgb=(0.18578,0.70489,0.48527)")
pushPGFPlotsPreamble("		rgb=(0.19109,0.70837,0.48228)")
pushPGFPlotsPreamble("		rgb=(0.19657,0.71183,0.47922)")
pushPGFPlotsPreamble("		rgb=(0.20222,0.71527,0.47608)")
pushPGFPlotsPreamble("		rgb=(0.20803,0.71870,0.47287)")
pushPGFPlotsPreamble("		rgb=(0.21400,0.72211,0.46959)")
pushPGFPlotsPreamble("		rgb=(0.22012,0.72551,0.46623)")
pushPGFPlotsPreamble("		rgb=(0.22640,0.72889,0.46279)")
pushPGFPlotsPreamble("		rgb=(0.23281,0.73225,0.45928)")
pushPGFPlotsPreamble("		rgb=(0.23937,0.73559,0.45569)")
pushPGFPlotsPreamble("		rgb=(0.24607,0.73891,0.45202)")
pushPGFPlotsPreamble("		rgb=(0.25290,0.74221,0.44828)")
pushPGFPlotsPreamble("		rgb=(0.25986,0.74549,0.44447)")
pushPGFPlotsPreamble("		rgb=(0.26694,0.74875,0.44057)")
pushPGFPlotsPreamble("		rgb=(0.27415,0.75199,0.43660)")
pushPGFPlotsPreamble("		rgb=(0.28148,0.75520,0.43255)")
pushPGFPlotsPreamble("		rgb=(0.28892,0.75839,0.42843)")
pushPGFPlotsPreamble("		rgb=(0.29648,0.76156,0.42422)")
pushPGFPlotsPreamble("		rgb=(0.30415,0.76470,0.41994)")
pushPGFPlotsPreamble("		rgb=(0.31193,0.76782,0.41559)")
pushPGFPlotsPreamble("		rgb=(0.31981,0.77091,0.41115)")
pushPGFPlotsPreamble("		rgb=(0.32780,0.77398,0.40664)")
pushPGFPlotsPreamble("		rgb=(0.33588,0.77702,0.40205)")
pushPGFPlotsPreamble("		rgb=(0.34407,0.78003,0.39738)")
pushPGFPlotsPreamble("		rgb=(0.35236,0.78301,0.39264)")
pushPGFPlotsPreamble("		rgb=(0.36074,0.78596,0.38781)")
pushPGFPlotsPreamble("		rgb=(0.36921,0.78889,0.38291)")
pushPGFPlotsPreamble("		rgb=(0.37778,0.79178,0.37794)")
pushPGFPlotsPreamble("		rgb=(0.38643,0.79464,0.37289)")
pushPGFPlotsPreamble("		rgb=(0.39517,0.79748,0.36776)")
pushPGFPlotsPreamble("		rgb=(0.40400,0.80027,0.36255)")
pushPGFPlotsPreamble("		rgb=(0.41291,0.80304,0.35727)")
pushPGFPlotsPreamble("		rgb=(0.42191,0.80577,0.35191)")
pushPGFPlotsPreamble("		rgb=(0.43098,0.80847,0.34648)")
pushPGFPlotsPreamble("		rgb=(0.44014,0.81114,0.34097)")
pushPGFPlotsPreamble("		rgb=(0.44937,0.81377,0.33538)")
pushPGFPlotsPreamble("		rgb=(0.45867,0.81636,0.32973)")
pushPGFPlotsPreamble("		rgb=(0.46805,0.81892,0.32400)")
pushPGFPlotsPreamble("		rgb=(0.47750,0.82144,0.31820)")
pushPGFPlotsPreamble("		rgb=(0.48703,0.82393,0.31232)")
pushPGFPlotsPreamble("		rgb=(0.49661,0.82638,0.30638)")
pushPGFPlotsPreamble("		rgb=(0.50627,0.82879,0.30036)")
pushPGFPlotsPreamble("		rgb=(0.51599,0.83116,0.29428)")
pushPGFPlotsPreamble("		rgb=(0.52578,0.83349,0.28813)")
pushPGFPlotsPreamble("		rgb=(0.53562,0.83579,0.28191)")
pushPGFPlotsPreamble("		rgb=(0.54552,0.83804,0.27563)")
pushPGFPlotsPreamble("		rgb=(0.55548,0.84025,0.26928)")
pushPGFPlotsPreamble("		rgb=(0.56550,0.84243,0.26288)")
pushPGFPlotsPreamble("		rgb=(0.57556,0.84457,0.25642)")
pushPGFPlotsPreamble("		rgb=(0.58568,0.84666,0.24990)")
pushPGFPlotsPreamble("		rgb=(0.59584,0.84872,0.24333)")
pushPGFPlotsPreamble("		rgb=(0.60604,0.85073,0.23671)")
pushPGFPlotsPreamble("		rgb=(0.61629,0.85271,0.23005)")
pushPGFPlotsPreamble("		rgb=(0.62658,0.85464,0.22335)")
pushPGFPlotsPreamble("		rgb=(0.63690,0.85654,0.21662)")
pushPGFPlotsPreamble("		rgb=(0.64726,0.85840,0.20986)")
pushPGFPlotsPreamble("		rgb=(0.65764,0.86022,0.20308)")
pushPGFPlotsPreamble("		rgb=(0.66805,0.86200,0.19629)")
pushPGFPlotsPreamble("		rgb=(0.67849,0.86374,0.18950)")
pushPGFPlotsPreamble("		rgb=(0.68894,0.86545,0.18272)")
pushPGFPlotsPreamble("		rgb=(0.69942,0.86712,0.17597)")
pushPGFPlotsPreamble("		rgb=(0.70990,0.86875,0.16926)")
pushPGFPlotsPreamble("		rgb=(0.72039,0.87035,0.16260)")
pushPGFPlotsPreamble("		rgb=(0.73089,0.87192,0.15603)")
pushPGFPlotsPreamble("		rgb=(0.74139,0.87345,0.14956)")
pushPGFPlotsPreamble("		rgb=(0.75188,0.87495,0.14323)")
pushPGFPlotsPreamble("		rgb=(0.76237,0.87642,0.13706)")
pushPGFPlotsPreamble("		rgb=(0.77285,0.87787,0.13111)")
pushPGFPlotsPreamble("		rgb=(0.78331,0.87928,0.12540)")
pushPGFPlotsPreamble("		rgb=(0.79376,0.88068,0.12001)")
pushPGFPlotsPreamble("		rgb=(0.80418,0.88205,0.11496)")
pushPGFPlotsPreamble("		rgb=(0.81458,0.88339,0.11035)")
pushPGFPlotsPreamble("		rgb=(0.82494,0.88472,0.10622)")
pushPGFPlotsPreamble("		rgb=(0.83527,0.88603,0.10265)")
pushPGFPlotsPreamble("		rgb=(0.84556,0.88732,0.09970)")
pushPGFPlotsPreamble("		rgb=(0.85581,0.88860,0.09745)")
pushPGFPlotsPreamble("		rgb=(0.86601,0.88987,0.09595)")
pushPGFPlotsPreamble("		rgb=(0.87617,0.89112,0.09525)")
pushPGFPlotsPreamble("		rgb=(0.88627,0.89237,0.09537)")
pushPGFPlotsPreamble("		rgb=(0.89632,0.89362,0.09634)")
pushPGFPlotsPreamble("		rgb=(0.90631,0.89485,0.09813)")
pushPGFPlotsPreamble("		rgb=(0.91624,0.89609,0.10072)")
pushPGFPlotsPreamble("		rgb=(0.92611,0.89733,0.10407)")
pushPGFPlotsPreamble("		rgb=(0.93590,0.89857,0.10813)")
pushPGFPlotsPreamble("		rgb=(0.94564,0.89982,0.11284)")
pushPGFPlotsPreamble("		rgb=(0.95530,0.90107,0.11813)")
pushPGFPlotsPreamble("		rgb=(0.96489,0.90232,0.12394)")
pushPGFPlotsPreamble("		rgb=(0.97442,0.90359,0.13021)")
pushPGFPlotsPreamble("		rgb=(0.98387,0.90487,0.13690)")
pushPGFPlotsPreamble("		rgb=(0.99325,0.90616,0.14394)")
pushPGFPlotsPreamble("	  }")
pushPGFPlotsPreamble("}");

max_stars_repo_name: sisl/VerifyGAN

max_issues_repo_name: sisl/VerifyGAN

max_forks_repo_name: sisl/VerifyGAN


text:
using IsoRank
using Base.Test

for k = 1:2:100
    G = sprand(k,k,0.4); G = G'+G; G = Int.(G - Diagonal(G) .> 0)
    R, res, L = isorank(G,G,maxiter=25,tol=1e-9,details=true,verbose=false)
    @test norm(L*vec(R) - res[1].*vec(R)) < 0.1
    b = rand(k,k);
    R, res, L = isorank(G,G,0.5,b,maxiter=25,tol=1e-9,details=true,verbose=false)
    @test norm(L*vec(R) - res[1].*vec(R)) < 0.1
end

function test1()
    G = sparse([1,2,3,1,3,2],[2,1,1,3,2,3],1,3,3)
    R, res, L = isorank(G,G,maxiter=25,tol=1e-9,details=true,verbose=false)
    norm(L*vec(R) .- res[1].*vec(R)) < 1e-2
end
@test test1()

function test2()
    G = sparse([1,2],[2,1],1,2,2)
    R,res,L = isorank(G, G, 0.85, ones(2,2), maxiter=25, tol=1e-9, details=true,verbose=false)
    vec(R) == [0.25,0.25,0.25,0.25] && norm(L*vec(R) - res[1].*vec(R)) == 0.0
end
@test test2()

function test3()
    G = sparse([1,2,3,1,3,2],[2,1,1,3,2,3],1,3,3)
    R,res,L = isorank(G, G, 0.85, ones(size(G,1),size(G,1)), maxiter=25, tol=1e-9, details=true,verbose=false)
    sum(vec(R) .== ones(9)./9)==0 && norm(L*vec(R) - res[1].*vec(R)) < 1e-2
end
@test test3()

function test4()
    diamond = sparse(vec([1 1 2 2 3]), vec([2 3 3 4 4]), 1.0,4,4)
    x = IsoRank.pagerank(diamond,verbose=false)
    isapprox(x, vec([0.128414 0.182991 0.260762 0.427833]), atol=0.01)
end
@test test4()

function test5()
    G1 = sprand(10,10,0.4); G1 = G1'+G1; G1 = Int.(G1 - Diagonal(G1) .> 0);
    G2 = sprand(15,15,0.4); G2 = G2'+G2; G2 = Int.(G2 - Diagonal(G2) .> 0);
    R, res, L = isorank(G1,G2,1.0,maxiter=25,tol=1e-9,details=true,verbose=false);
    A = Float64.(kron(G2,G1)); S = 1.0 ./ (A' * ones(Float64,size(A,2))); scale!(A,S);
    @test norm(A*vec(R) - res[1].*vec(R)) < 0.001
    els,evs = eigs(A);
    @test norm(A*evs[:,1] - L*evs[:,1]) < 0.001
end
test5()


text:
# Julian program to Find the n'th Catalan number
# To know more about Catalan Numbers Visit Here : https://en.wikipedia.org/wiki/Catalan_number

function find_catalan(n)
    if (n == 0 || n == 1)
        return 1
    end

    # DP Table
    cat = [0 for i in 1:(n + 1)]

    # Initialize the first to values
    cat[1] = 1
    cat[2] = 1
    for i in 3:(n + 1)
        for j in 1:(i - 1)
            cat[i] = cat[i] +  (cat[j] * cat[i - j])
        end
    end
    return cat[n + 1]
end

print("Enter the number: ")
num = readline()
num = parse(Int, num)
catalan_number = find_catalan(abs(num))
println("The $num'th catalan number is $catalan_number.")

"""
Time Complexity: O(num^2), where 'num' is the given number
Space Complexity: O(num)

SAMPLE INPUT AND OUTPUT

SAMPLE 1
Enter the number: 12
The 12'th catalan number is 208012.

SAMPLE 2
Enter the number: 2
The 2'th catalan number is 2.

"""


text:
"""
An `OperatedDict` represents a set that is acted on by an operator, for example
the differentiation operator. The `OperatedDict` has the dimension of the source
set of the operator, but each basis function is acted on by the operator.

In practice the operator `D` acts on the coefficients. Thus, an expansion
in an `OperatedDict` with coefficients `c` corresponds to an expansion in the
underlying dictionary with coefficients `Dc`. If `D` represents differentiation,
then the `OperatedDict` effectively represents the dictionary of derivatives of
the underlying dictionary elements.
"""
struct OperatedDict{S,T} <: HighlySimilarDerivedDict{S,T}
    "The operator that acts on the set"
    op          ::  DictionaryOperator{T}

    scratch_src
    scratch_dest

    function OperatedDict{S,T}(op::DictionaryOperator) where {S,T}
        scratch_src = zeros(src(op))
        scratch_dest = zeros(dest(op))
        new(op, scratch_src, scratch_dest)
    end
end

superdict(dict::OperatedDict) = src(dict)

function OperatedDict(op::DictionaryOperator{T}) where {T}
    S = domaintype(src(op))
    OperatedDict{S,T}(op)
end

OperatedDict(op::IdentityOperator) = src(op)

## Printing

stencil_parentheses(dict::OperatedDict) = true
object_parentheses(dict::OperatedDict) = true
Display.displaystencil(d::OperatedDict) = _stencil(d, superdict(d), operator(d))
_stencil(d::OperatedDict, dict, op) = [op, " * ", dict]

src(s::OperatedDict) = src(s.op)

dest(s::OperatedDict) = dest(s.op)

domaintype(s::OperatedDict) = domaintype(src(s))

operator(set::OperatedDict) = set.op

function similar(d::OperatedDict, ::Type{T}, dims::Int...) where {T}
    # We enforce an equal length since we may be able to resize dictionaries,
    # but we can not in general resize the operator (we typically do not store what
    # the operator means in its type).
    @assert length(d) == prod(dims)
    # not sure what to do with dest(d) here - in the meantime invoke similar
    OperatedDict(similar_operator(operator(d), similar(src(d), T), similar(dest(d),T) ))
end

function similardictionary(d::OperatedDict, dict::Dictionary)
    if dict == src(d)
        d
    elseif src(d) == dest(d)
        OperatedDict(wrap_operator(dict, dict, operator(d)))
    end
end

==(op1::OperatedDict, op2::OperatedDict) =
	operator(op1) == operator(op2) && src(op1) == src(op2)

≈(op1::OperatedDict, op2::OperatedDict) =
	operator(op1) ≈ operator(op2) && src(op1) == src(op2)

for op in (:support, :length)
    @eval $op(s::OperatedDict) = $op(src(s))
end

# We don't know in general what the support of a specific basis functions is.
# The safe option is to return the support of the set itself for each element.
for op in (:support,)
    @eval $op(s::OperatedDict, idx) = $op(src(s))
end

dict_in_support(set::OperatedDict, i, x) = in_support(superdict(set), x)

zeros(::Type{T}, s::OperatedDict) where {T} = zeros(T, src(s))
tocoefficientformat(a, d::OperatedDict) = tocoefficientformat(a, src(d))

##########################
# Indexing and evaluation
##########################

ordering(dict::OperatedDict) = ordering(src(dict))

unsafe_eval_element(dict::OperatedDict, i, x) =
    _unsafe_eval_element(dict, i, x, operator(dict), dict.scratch_src, dict.scratch_dest)

function _unsafe_eval_element(dict::OperatedDict, idxn, x, op, scratch_src, scratch_dest)
    idx = linear_index(dict, idxn)

    scratch_src[idx] = 1
    apply!(op, scratch_dest, scratch_src)
    scratch_src[idx] = 0
    eval_expansion(dest(dict), scratch_dest, x)
end

function _unsafe_eval_element(dict::OperatedDict, idxn, x, op::ScalingOperator, scratch_src, scratch_dest)
    idx = linear_index(dict, idxn)
    diag(op, idx) * unsafe_eval_element(src(dict), idxn, x)
end

evaluation(::Type{T}, dict::OperatedDict, gb::GridBasis, grid::AbstractGrid; options...) where {T} =
    wrap_operator(dict, gb, evaluation(T, dest(dict), gb, grid; options...) * operator(dict))

## Properties

# Perhaps we can extend the underlying dictionary, but in general we can not extend the operator
hasextension(dict::OperatedDict) = false

isreal(dict::OperatedDict) = isreal(operator(dict))

for op in (:isbasis, :isframe)
    @eval $op(set::OperatedDict) = isdiag(operator(set)) && $op(src(set))
end

#################
# Transform
#################

hastransform(dict::OperatedDict) = isdiag(operator(dict)) && hastransform(src(dict))
hastransform(dict::OperatedDict, dgs::GridBasis) = isdiag(operator(dict)) && hastransform(src(dict), dgs)

transform_dict(dict::OperatedDict; options...) = transform_dict(superdict(dict); options...)

function transform_to_grid(T, src::OperatedDict, dest::GridBasis, grid; options...)
	op = transform_to_grid(T, superdict(src), dest, grid; options...) * operator(src)
	wrap_operator(src, dest, op)
end

function transform_from_grid(T, src::GridBasis, dest::OperatedDict, grid; options...)
	 op = inv(operator(dest)) * transform_from_grid(T, src, superdict(dest), grid; options...)
	 wrap_operator(src, dest, op)
 end

#################
# Differentiation
#################

diff(dict::OperatedDict, args...; options...) = differentiation(dest(dict), args...; options...) * dict

hasderivative(dict::OperatedDict) = hasderivative(superdict(dict))
hasderivative(dict::OperatedDict, order) = hasderivative(superdict(dict), order)

# TODO: allow a diagonal operator
hasantiderivative(dict::OperatedDict) = false

function derivative_dict(dict::OperatedDict, order; options...)
	if orderiszero(order)
		dict
	else
		derivative_dict(dest(dict), order; options...)
	end
end

differentiation(::Type{T}, dsrc::OperatedDict, ddest::Dictionary, order; options...) where {T} =
	differentiation(T, dest(dsrc), ddest, order; options...) * operator(dsrc)

function unsafe_eval_element_derivative(dict::OperatedDict, i, x, order)
	@assert order == 1
    _unsafe_eval_element_derivative(dict, i, x, order, operator(dict), dict.scratch_src, dict.scratch_dest)
end

function _unsafe_eval_element_derivative(dict::OperatedDict, idxn, x, order, op, scratch_src, scratch_dest)
    idx = linear_index(dict, idxn)
    if isdiag(op)
        diag(op, idx) * unsafe_eval_element_derivative(src(dict), idxn, x, order)
    else
        scratch_src[idx] = 1
        apply!(op, scratch_dest, scratch_src)
        scratch_src[idx] = 0
		D = differentiation(superdict(dict))
        eval_expansion(dest(D), D*scratch_dest, x)
    end
end

#################
# Projections
#################

hasmeasure(dict::OperatedDict) = hasmeasure(superdict(dict))

measure(dict::OperatedDict) = measure(superdict(dict))

for f in (:isorthogonal, :isorthonormal, :isbiorthogonal)
    @eval $f(dict::OperatedDict, measure::Measure) =
        $f(superdict(dict), measure) && $f(operator(dict))
end

gram1(T, dict::OperatedDict, measure; options...) =
	adjoint(operator(dict)) * gram(T, dest(operator(dict)), measure; options...) * operator(dict)

function innerproduct1(d1::OperatedDict, i, d2, j, measure; options...)
	op = operator(d1)
	if isdiag(op)
		idx = linear_index(d1, i)
		conj(diag(op, idx)) * innerproduct(superdict(d1), i, d2, j, measure; options...)
	else
		innerproduct2(d1, i, d2, j, measure; options...)
	end
end

function innerproduct2(d1, i, d2::OperatedDict, j, measure; options...)
	op = operator(d2)
	if isdiag(op)
		idx = linear_index(d2, j)
		diag(op, idx) * innerproduct(d1, i, superdict(d2), j, measure; options...)
	else
		default_dict_innerproduct(d1, i, d2, j, measure; options...)
	end
end

function mixedgram1(T, dict1::OperatedDict, dict2, measure; options...)
	G = mixedgram(T, superdict(dict1), dict2, measure; options...)
	wrap_operator(dict2, dict1, adjoint(operator(dict1)) * G)
end

function mixedgram2(T, dict1, dict2::OperatedDict, measure; options...)
	G = mixedgram(T, dict1, superdict(dict2), measure; options...)
	wrap_operator(dict2, dict1, G *  operator(dict2))
end

#################
# Special cases
#################
(*)(a::Number, s::Dictionary) = ScalingOperator(s, a) * s
(*)(a::Number, s::OperatedDict) = (ScalingOperator(dest(s), a) * operator(s)) * superdict(s)

function (*)(op::DictionaryOperator, dict::Dictionary)
    @assert src(op) == dict
    OperatedDict(op)
end

function (*)(op::DictionaryOperator, dict::OperatedDict)
    @assert src(op) == dest(dict)
    OperatedDict(op*operator(dict))
end


text:
using Compat # for Nullable
using Colors
using Lazy

#### Model ####

@defonce immutable Board{lost}
    uncovered::AbstractMatrix
    mines::AbstractMatrix
end

newboard(m, n, minefraction=0.2) =
    Board{false}(fill(-1, (m, n)), rand(m, n) .< minefraction)

function mines_around(board, i, j)
    m, n = size(board.mines)

    a = max(1, i-1)
    b = min(i+1, m)
    c = max(1, j-1)
    d = min(j+1, n)

    sum(board.mines[a:b, c:d])
end

### Update ###

update(board::Board{true}, move) = board
function update(board, move)
    i, j = move
    if board.mines[i, j]
        return Board{true}(board.uncovered, board.mines) # Game over
    else
        uncovered = copy(board.uncovered)
        if uncovered[i, j] == -1
            uncovered[i, j] = mines_around(board, i, j)
        end
        return Board{false}(uncovered, board.mines)
    end
end

moves_signal = Signal((0, 0))
initial_board_signal = Signal(Board, newboard(10, 10))
board_signal = flatten(
    map(initial_board_signal) do b
        foldp(update, b, moves_signal; typ=Board)
    end
)

### View ###

colors = ["#fff"; colormap("reds", 7)]

box(content, color) =
    inset(Escher.middle,
        fillcolor(color, size(4em, 4em, empty)),
        Escher.fontsize(2em, content)) |> paper(1) |> Escher.pad(0.2em)

number(x) = box(x == -1 ? "" : string(x) |> fontweight(800), colors[x+2])
mine = box(icon("report"), "#e58")
block(board::Board{true}, i, j) =
    board.mines[i, j] ? mine :
        number(board.uncovered[i, j])

block(board, i, j) =
     intent(constant((i, j)), clickable(number(board.uncovered[i, j]))) >>> moves_signal

gameover = vbox(
        title(2, "Game Over!") |> Escher.pad(1em),
        intent(_ -> newboard(10, 10), button("Start again")) >>> initial_board_signal
    ) |> Escher.pad(1em) |> fillcolor("white")

function showboard{lost}(board::Board{lost})
    m, n = size(board.mines)
    b = hbox([vbox([block(board, i, j) for j in 1:m]) for i in 1:n])
    lost ? inset(Escher.middle, b, gameover) : b
end

function main(window)
    push!(window.assets, "widgets")

    vbox(
       vskip(2em),
       title(3, "minesweeper"),
       vskip(2em),
       map(showboard, board_signal, typ=Tile),
    ) |> packacross(center)
end