text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# 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.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# 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.
# *********************************************************************************
"""
	add_generator_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")

Adds the Generator results to the dictionary passed back from `run_reopt` using the solved model `m` and the `REoptInputs` for node `_n`.
Note: the node number is an empty string if evaluating a single `Site`.

Generator results:
- `size_kw` Optimal generator capacity
- `lifecycle_fixed_om_cost_after_tax` Lifecycle fixed operations and maintenance cost in present value, after tax
- `year_one_fixed_om_cost__before_tax` fixed operations and maintenance cost over the first year, before considering tax benefits
- `lifecycle_variable_om_cost_after_tax` Lifecycle variable operations and maintenance cost in present value, after tax
- `year_one_variable_om_cost_before_tax` variable operations and maintenance cost over the first year, before considering tax benefits
- `lifecycle_fuel_cost_after_tax` Lifecycle fuel cost in present value, after tax
- `year_one_fuel_cost_before_tax` Fuel cost over the first year, before considering tax benefits
- `average_annual_fuel_used_gal` Gallons of fuel used in each year
- `year_one_to_battery_series_kw` Vector of power sent to battery in year one
- `year_one_to_grid_series_kw` Vector of power sent to grid in year one
- `year_one_to_load_series_kw` Vector of power sent to load in year one
- `year_one_energy_produced_kwh` Total energy produced in year one
- `average_annual_energy_produced_kwh` Average annual energy produced over analysis period
"""
function add_generator_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
    r = Dict{String, Any}()

	GenPerUnitSizeOMCosts = @expression(m, p.third_party_factor * p.pwf_om * sum(m[:dvSize][t] * p.om_cost_per_kw[t] for t in p.techs.gen))

	GenPerUnitProdOMCosts = @expression(m, p.third_party_factor * p.pwf_om * p.hours_per_time_step *
		sum(m[:dvRatedProduction][t, ts] * p.production_factor[t, ts] * p.s.generator.om_cost_per_kwh
			for t in p.techs.gen, ts in p.time_steps)
	)
	r["size_kw"] = round(value(sum(m[:dvSize][t] for t in p.techs.gen)), digits=2)
	r["lifecycle_fixed_om_cost_after_tax"] = round(value(GenPerUnitSizeOMCosts) * (1 - p.s.financial.owner_tax_pct), digits=0)
	r["lifecycle_variable_om_cost_after_tax"] = round(value(m[:TotalPerUnitProdOMCosts]) * (1 - p.s.financial.owner_tax_pct), digits=0)
	r["lifecycle_fuel_cost_after_tax"] = round(value(m[:TotalGenFuelCosts]) * (1 - p.s.financial.offtaker_tax_pct), digits=2)
	r["year_one_fuel_cost_before_tax"] = round(value(m[:TotalGenFuelCosts]) / p.pwf_fuel["Generator"], digits=2)
	r["year_one_variable_om_cost_before_tax"] = round(value(GenPerUnitProdOMCosts) / (p.pwf_om * p.third_party_factor), digits=0)
	r["year_one_fixed_om_cost_before_tax"] = round(value(GenPerUnitSizeOMCosts) / (p.pwf_om * p.third_party_factor), digits=0)

	if !isempty(p.s.storage.types.elec)
	generatorToBatt = @expression(m, [ts in p.time_steps],
		sum(m[:dvProductionToStorage][b, t, ts] for b in p.s.storage.types.elec, t in p.techs.gen))
	else
		generatorToBatt = zeros(length(p.time_steps))
	end
	r["year_one_to_battery_series_kw"] = round.(value.(generatorToBatt), digits=3)

	generatorToGrid = @expression(m, [ts in p.time_steps],
		sum(m[:dvProductionToGrid][t, u, ts] for t in p.techs.gen, u in p.export_bins_by_tech[t])
	)
	r["year_one_to_grid_series_kw"] = round.(value.(generatorToGrid), digits=3)

	generatorToLoad = @expression(m, [ts in p.time_steps],
		sum(m[:dvRatedProduction][t, ts] * p.production_factor[t, ts] * p.levelization_factor[t]
			for t in p.techs.gen) -
			generatorToBatt[ts] - generatorToGrid[ts]
	)
	r["year_one_to_load_series_kw"] = round.(value.(generatorToLoad), digits=3)

    GeneratorFuelUsed = @expression(m, sum(m[:dvFuelUsage][t, ts] for t in p.techs.gen, ts in p.time_steps) / KWH_PER_GAL_DIESEL)
	r["average_annual_fuel_used_gal"] = round(value(GeneratorFuelUsed), digits=2)

	Year1GenProd = @expression(m,
		p.hours_per_time_step * sum(m[:dvRatedProduction][t,ts] * p.production_factor[t, ts]
			for t in p.techs.gen, ts in p.time_steps)
	)
	r["year_one_energy_produced_kwh"] = round(value(Year1GenProd), digits=0)
	AverageGenProd = @expression(m,
		p.hours_per_time_step * sum(m[:dvRatedProduction][t,ts] * p.production_factor[t, ts] *
		p.levelization_factor[t]
			for t in p.techs.gen, ts in p.time_steps)
	)
	r["average_annual_energy_produced_kwh"] = round(value(AverageGenProd), digits=0)

	d["Generator"] = r
    nothing
end

function add_generator_results(m::JuMP.AbstractModel, p::MPCInputs, d::Dict; _n="")
    r = Dict{String, Any}()

	r["variable_om_cost"] = round(value(m[:TotalPerUnitProdOMCosts]), digits=0)
	r["fuel_cost"] = round(value(m[:TotalGenFuelCosts]), digits=2)

    if p.s.storage.attr["ElectricStorage"].size_kw > 0
        generatorToBatt = @expression(m, [ts in p.time_steps],
            sum(m[:dvProductionToStorage][b, t, ts] for b in p.s.storage.types.elec, t in p.techs.gen))
        r["to_battery_series_kw"] = round.(value.(generatorToBatt), digits=3).data
    else
        generatorToBatt = zeros(length(p.time_steps))
    end

	generatorToGrid = @expression(m, [ts in p.time_steps],
		sum(m[:dvProductionToGrid][t, u, ts] for t in p.techs.gen, u in p.export_bins_by_tech[t])
	)
	r["to_grid_series_kw"] = round.(value.(generatorToGrid), digits=3).data

	generatorToLoad = @expression(m, [ts in p.time_steps],
		sum(m[:dvRatedProduction][t, ts] * p.production_factor[t, ts] * p.levelization_factor[t]
			for t in p.techs.gen) -
			generatorToBatt[ts] - generatorToGrid[ts]
	)
	r["to_load_series_kw"] = round.(value.(generatorToLoad), digits=3).data

    GeneratorFuelUsed = @expression(m, sum(m[:dvFuelUsage][t, ts] for t in p.techs.gen, ts in p.time_steps) / KWH_PER_GAL_DIESEL)
	r["average_annual_fuel_used_gal"] = round(value(GeneratorFuelUsed), digits=2)

	Year1GenProd = @expression(m,
		p.hours_per_time_step * sum(m[:dvRatedProduction][t,ts] * p.production_factor[t, ts]
			for t in p.techs.gen, ts in p.time_steps)
	)
	r["energy_produced_kwh"] = round(value(Year1GenProd), digits=0)

	d["Generator"] = r
    nothing
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
#This module works as a example of solving linear programming problems for Ju4ja
module RereDmlLpSolver
using JuMP
using GLPK
export solveDmlLp

  #solve linear programming
  function solveDmlLp(c::Vector,A::Matrix,b::Vector,ifReg::Bool=true)
    model = JuMP.Model(with_optimizer(GLPK.Optimizer))
    #println(length(c))
    @variable(model, x[i=1:length(c)])
    @constraint(model, con, A * x .>= b)
    @constraint(model, con0, x .>= 0)
    @objective(model, Min, c'*x )
    status = optimize!(model)
    result=JuMP.value.(x)
    obj=JuMP.objective_value(model)
    println("x = ", result)
    println("Objective value: ", obj)
    return result
  end

end

ext: jl


text:
@testset "gradient.jl" begin

end


text:
#' ## Run this cell to load the graphics packages
using Plots
gr()

#' ## Run this cell to use a couple of definitions from the previous
#' Tools notebook

sigma(x) = 1 / (1 + exp(-x))
f(x, w) = sigma(w * x)

#' ## Multiple function parameters
#'
#' In notebook 6, we saw how we could adjust a parameter to make a curve fit
#' a single data point. What if there is more data to fit?
#'
#' We'll see that as we acquire more data to fit, it can sometimes be useful
#' to add complexity to our model via additional parameters.

#' ## Multiple function parameters
#'
#' In notebook 6, we saw how we could adjust a parameter to make a curve fit
#' a single data point. What if there is more data to fit?
#'
#' We'll see that as we acquire more data to fit, it can sometimes be useful
#' to add complexity to our model via additional parameters.

#' ## Adding more data
#' Suppose there are now two data points to fit, the previous $(x_0, y_0)$
#' used in notebook 6 and also $(x_1, y_1) = (-3, 0.3)$.

#' #### Exercise 1
#
# Make an interactive plot of the function $f_w$ together with the two data
# points. Can you make the graph of $f_w$ pass through *both* data points
# at the *same* time?

x0, y0 = 2, 0.8
x1, y1 = -3, 0.3
w = 0.4
plot(x->f(x, w), -5, 5, ylims=(0, 1), label="f")
scatter!([x0,x1], [y0,y1], label="data")

#' You should have found that it's actually *impossible* to fit both data points
#' at the same time! The best we could do is to *minimise* how far away the
#' function is from the data. To do so, we need to somehow balance the distance
#' from each of the two data points.

#' #### Exercise 2
#'
#' Play with the slider to find the value $w^*$ of $w$ that you think has the "least error" in an intuitive sense.

w=0.45
plot(x->f(x, w), -5, 5, ylim=(0,1))
scatter!([x0, x1], [y0, y1])

#' ## Defining a loss function
#' How can we *quantify* some kind of overall measure of the distance from *all* of the data? Well, we just need to define a new loss function! One way to do so would be to sum up the loss functions for each data point, i.e. the sum of the squared vertical distances from the graph to each data point:
#'
#' $$L(w) = [y_0 - f_w(x_0)]^2 + [y_1 - f_w(x_1)]^2.$$
#'
#' Since the two pieces that we are adding have the same mathematical "shape" or structure, we can abbreviate this by writing it as a sum:
#'
#' $$L(w) = \sum_{i=0}^1 [y_i - f_w(x_i)]^2.$$

#' So now we want to find the value $w^*$ of $w$ that minimizes this new
#' function $L$!

#' #### Exercise 3
#'
#' Make a visualization to show the function $f_w$ and to visualize the
#' distance from each of the data points.
xs = [2, -3]
ys = [0.8, 0.3]

L(w) = sum( (ys .- f.(xs, w)) .^ 2 )

begin
    w = 0.3 # Try manipulating w between -2 and 2

    plot(x->f(x, w), -5, 5, ylims=(0, 1), lw = 3, label="f_w")

    scatter!(xs, ys, label="data")

    for i in 1:2
        plot!([xs[i], xs[i]], [ys[i], f(xs[i], w)], c=:green)
    end

    title!("L(w) =  $(round(L(w); sigdigits = 5))")
end

#' #### Exercise 4
#'
#' After playing with this for a while, it is intuitively obvious that we
# cannot make the function pass through both data points for any value of
# $w$. In other words, our loss function, `L(w)`, is never zero.
#'
#' What is the minimum value of `L(w)` that you can find by altering `w`?
#' What is the corresponding value of `w`?
w = 0.42
plot(x->f(x, w), -5, 5, ylims=(0, 1), lw = 3, label="f_w", legend=(0.8,0.1))

scatter!(xs, ys, label="data")

for i in 1:2
    plot!([xs[i], xs[i]], [ys[i], f(xs[i], w)], c=:green)
end

title!("L($w) =  $(round(L(w); sigdigits = 5))")

#' ### Sums in Julia
#' To generate the above plot we used the `sum` function.
#' `sum` can add together all the elements of a collection or range,
#' or it can add together the outputs of applying a function to all the
#' elements of a collection or range.
#'
#' Look up the docs for `sum` via
#'
#' ```julia
#' ?sum
#' ```
#' if you need more information.

#' #### Exercise 6
#'
#' What is the sum of the absolute values of all integers between -3 and 3? Use `sum` and the `abs` function.

#' #### Solution

sum(abs, -3:3)

#' ## What does the loss function $L$ look like?
#' In our last attempt to minimize `L(w)` by varying `w`, we saw that `L(w)`
#' always seemed to be greater than 0.  Is that true? Let's plot it to find out!

#' #### Exercise 7
#'
#' Plot the new loss function $L(w)$ as a function of $w$.
#' #### Solution
plot(L, -2, 2, xlabel="w", ylabel="L(w)", ylims=(0, 1.2))

#' ### Features of the loss function
#' The first thing to notice is that $L$ is always positive.
#' Since it is the sum of squares, and squares cannot be negative,
#' the sum cannot be negative either!
#'
#' However, we also see that although $L$ dips close to $0$ for a single,
#' special value $w^* \simeq 0.4$, it never actually *reaches* 0!
#' Again we could zoom in on that region of the graph to estimate it more
#' precisely.

#' We might start suspecting that there should be a better way of using the
#' computer to minimize $L$ to find the location $w^*$ of the minimum,
#' rather than still doing everything by eye. Indeed there is, as we will
#' see in the next two notebooks!

#' ## Adding more parameters to the model

#' If we add more parameters to a function, we may be able to improve how it
#' fits to data. For example, we could define a new function $g$ with another
#' parameter, a shift or **bias**:
#'
#' $$g(x; w, b) := \sigma(w \, x + b).$$
#+ results="hidden"
g(x, w, b) = sigma(w*x + b)
#+

xs = [2, -3]
ys = [0.8, 0.3]

L2D(w, b) = sum( (ys .- g.(xs, w, b)) .^ 2)
w = 0.45 # Try manipulating w between -2 and 2
b = 0.5 # Try manipulating b between -2 and 2
plot(x -> g(x, w, b), -5, 5, ylims=(0, 1), lw=3)

scatter!(xs, ys)

for i in 1:2
    plot!([xs[i], xs[i]], [ys[i], g(xs[i], w, b)])
end

title!("L2D(w, b) = $(L2D(w, b))")

#' You should be able to convince yourself that we can now make the curve pass
#' through both points simultaneously.

#' #### Exercise 9
#'
#' For what values of `w` and `b` does the line pass through both points?

w = 0.45
b = 0.48
plot(x -> g(x, w, b), -5, 5, ylims=(0, 1), lw=3)

scatter!(xs, ys)

for i in 1:2
    plot!([xs[i], xs[i]], [ys[i], g(xs[i], w, b)])
end

title!("L2D(w, b) = $(L2D(w, b))")
#' ## Fitting both data points: a loss function
#' Following the procedure that we used when we had a single parameter,
#' we can think of the fitting procedure once again as minimizing a suitable
#' *loss function*. The expression for the loss function is almost the same,
#' except that now $f$ has two parameters $w$ and $b$, so the loss function
#' $L_2$ is itself a function of $w$ *and* $b$:
#'
#' $$L_2(w, b) = \sum_i [y_i - f(x_i; w, b)]^2.$$
#'
#' So we want to minimize this loss function over *both* of the parameters
#' $w$ and $b$! Let's plot it.

#' #### Exercise 10
#'
#' Define the function `L2` in Julia.

#' #### Solution

L2(w, b) = sum( (ys .- g.(xs, w, b)) .^ 2 )

#' To plot the loss as a function of *two* variables (the weights *and* the bias),
#' we will make use of the `surface` function.  To get a nice interactive
# 3D plot, we will use the PlotlyJS "backend" (plotting engine):

gr()  # load the backend

#-

ws = -2:0.05:2
bs = -2:0.05:2

contour(ws, bs, L2, levels=0.001:0.025:1.25, xlabel="value of w", ylabel="value of b")

#' We can see that indeed there is a unique point $(w^*, b^*)$ where the
#' function $L_2$ attains its minimum. You can try different ranges for the
#' $x$, $y$ and $z$ coordinates to get a better view.

#' ## More data

#' If we add more data, however, we will again not be able to fit all of the
#' data; we will only be able to attain a "best fit".
#'
#' Let's create `xs` and `ys` with some more data:

xs = [2, -3, -1, 1]
ys = [0.8, 0.3, 0.4, 0.4]

#' #### Exercise 11
#'
#' a) Make an interactive visualization of the function $f$ and the data. Try
#' to find the values of $w$ and $b$ that give the best fit.
#'
#' b) Define the loss function and plot it.
#'
#' ```julia
#' sigma(x) = 1 / (1 + exp(-x))
#' f(x, w) = sigma(w * x)
#' g(x, w, b) = sigma(w*x + b)
#' ```
L2(w, b) = sum( (ys .- g.(xs, w, b)) .^ 2 )

w = 0.4 # Try manipulating w between -2 and 2
b = -0.1 # Try manipulating b between -2 and 2
plot(x->g(x, w, b), -5, 5, ylims=(-0.2, 1))

scatter!(xs, ys)

for i in 1:length(xs)
    plot!([xs[i], xs[i]], [ys[i], g(xs[i], w, b)])
end

title!("L2(w, b) = $(round(L2(w, b); sigdigits = 5))")

#' Let's define the loss function that we're using above so that we can plot
#' it as a function of the parameters `w` and `b`.
#' #### Exercise 12
#'
#' We've seen the loss function, $$L_{2D}(w, b) = \sum_i[\mathrm{ys}_i - g(\mathrm{xs}_i, w, b)]^2,$$ written in Julia as
#
# ```julia
# L2(w, b) = sum( (ys .- g.(xs, w, b)) .^ 2 )
# ```
#
# a few times now. To ensure you understand what this function is doing, implement your own loss function using the commented code below. Do this without using `sum`, `abs2`, or broadcasting dot syntax (for example, `.-`). Hint: you'll want to use a `for` loop to do this.

function myL2(w, b)
    loss = 0.0
    for (x, y) in zip(xs, ys)
        loss += (y - g(x, w, b))^2
    end
    return loss
end

#' Now that you've defined `L2`, we can plot it using either `surface` or
#' `contour` to see how the loss changes as `w` and `b` change!

#' ## plotlyjs() # Use Plotly for 3d surface plots (if possible)

ws = -2:0.05:2
bs = -2:0.05:2

#' ## surface(ws, bs, L2, alpha=0.8, zlims=(0,3))  # alpha gives transparency
contour(ws, bs, L2, levels=0.055:0.05:1.5, xlabel="value of w", ylabel="value of b")


text:
export Timestamp

using Printf

"""
    struct Timestamp

A simple time-keeping struct with space for minutes and (fractional) seconds.
"""
mutable struct Timestamp
    minutes::Int64
    seconds::Real

    function Timestamp(minutes::Int64, seconds::Real)
        new(minutes, seconds)
    end

    function Timestamp(seconds::Real)
        seconds_over = mod(seconds, 60)
        minutes = floor(Int64, (seconds - seconds_over) / 60)
        new(minutes, seconds_over)
    end
end

function Base.show(io::IO, ts::Timestamp)
    seconds_fmt = @sprintf "%05.2f" ts.seconds
    print(io, "$(ts.minutes) mins $(seconds_fmt) secs")
end

function Base.:+(x::Timestamp, y::Timestamp)
    total_minutes = x.minutes + y.minutes
    total_seconds = x.seconds + y.seconds
    if total_seconds ≥ 60
        total_minutes += 1
        total_seconds -= 60
    end
    return Timestamp(total_minutes, total_seconds)
end

function Base.:+(x::Timestamp, y::Real)
    ts2 = Timestamp(y)
    return x + ts2
end
# TODO: Figure out why we don't automatically get commutation here...
Base.:+(x::Real, y::Timestamp) = y + x

function Base.:/(x::Timestamp, y::Real)
    total_seconds = x.minutes*60 + x.seconds
    divided_seconds = total_seconds / y
    return Timestamp(divided_seconds)
end

function Base.:%(x::Timestamp, y::Real)
    total_seconds = x.minutes*60 + x.seconds
    modulo = total_seconds % y
    return Timestamp(modulo)
end

function Base.:(==)(x::Timestamp, y::Timestamp)
    return (x.minutes == y.minutes) && (x.seconds == y.seconds)
end

function Base.:(<)(x::Timestamp, y::Timestamp)
    x_tot_s = x.minutes*60 + x.seconds
    y_tot_s = y.minutes*60 + y.seconds
    return x_tot_s < y_tot_s
end

function Base.:(<=)(x::Timestamp, y::Timestamp)
    return (x == y) || (x < y)
end


text:
module FundamentalsNumericalComputation

export FNC
FNC = FundamentalsNumericalComputation

using Reexport

@reexport using LinearAlgebra
@reexport using SparseArrays
@reexport using Polynomials
@reexport using NLsolve
@reexport using DifferentialEquations
@reexport using Plots
@reexport using PrettyTables
@reexport using QuadGK
@reexport using Interpolations

@info "Exporting: LinearAlgebra,SparseArrays,Polynomials,NLsolve,Interpolations,DifferentialEquations,DataFrames,Plots"

include("chapter01.jl")
include("chapter02.jl")
include("chapter03.jl")
include("chapter04.jl")
include("chapter05.jl")
include("chapter06.jl")
include("chapter08.jl")
include("chapter09.jl")
include("chapter10.jl")
include("chapter11.jl")
include("chapter13.jl")

function magic(n)
	# Adapted from an implementation by Phillip Berndt
	if isodd(n)
		p = 1:n
		A = p .- div(n+3, 2) .+ p'
		B = 2p .- 2 .+ p'
		return n*mod.(A,n) + mod.(B,n) .+ 1
	elseif iseven(n) && iseven(div(n,2))
		J = @. div((1:n)%4, 2)
		K = J' .== J
		M = (0:n-1) .+ (1:n:n^2)'
		M[K] = n^2 + 1 .- M[K]
		return M
	else
		p = div(n,2)
		M = magic(p)
		M = [M M.+2p^2; M.+3p^2 M.+p^2]
		(n == 2) && return M
		i = 1:p
		k = div(n-2,4)
		j = [(1:k); ((n-k+2):n)]
		M[[i; i.+p],j] = M[[i.+p; i],j]
		i = k+1
		j = [1; i]
		M[[i; i+p],j] = M[[i+p; i],j]
		return M
	end
end

end


text:
@testset "RSSItem" begin
    rss = F.RSSItem(
        "title", "www.link.com", "description", "", "author@author.com", "category",
        "www.comments.com", "enclosure", Date(2012,12,12))
    @test rss.title == "title"
    @test rss.link == "www.link.com"
    @test rss.description == "description"
    @test rss.author == "author@author.com"
    @test rss.category == "category"
    @test rss.comments == "www.comments.com"
    @test rss.enclosure == "enclosure"
    @test rss.pubDate == Date(2012,12,12)
end

@testset "RSSbasics" begin
    fs()
    empty!(F.RSS_DICT)
    F.def_GLOBAL_VARS!()
    F.set_var!(F.GLOBAL_VARS, "website_title", "Website title")
    F.set_var!(F.GLOBAL_VARS, "website_descr", "Website descr")
    F.set_var!(F.GLOBAL_VARS, "website_url", "https://github.com/tlienart/Franklin.jl/")

    # Page 1 with tags = ["foo", "t.a g"]
    F.def_LOCAL_VARS!()
    set_curpath("hey/hello.md")
    F.set_var!(F.LOCAL_VARS, "rss_title", "title 1")
    F.set_var!(F.LOCAL_VARS, "rss", "Page with tag foo.")
    F.set_var!(F.LOCAL_VARS, "rss_pubdate", Date(2020, 10, 27))
    F.set_var!(F.LOCAL_VARS, "tags", ["foo", "t.a g"])

    item, tags = F.add_rss_item()
    @test item.title == "title 1"
    @test item.description == "Page with tag foo.\n"
    @test item.author == ""
    @test item.pubDate == Date(2020, 10, 27)
    @test tags == ["foo", "t.a g"]
    @test F.RSS_DICT["/hey/hello/index.html"][1].description == item.description
    @test F.RSS_DICT["/hey/hello/index.html"][2] == ["foo", "t.a g"]

    # Page 2 with tags = ["foo", "bar"]
    F.def_LOCAL_VARS!()
    set_curpath("hey/ho.md")
    F.set_var!(F.LOCAL_VARS, "rss_title", "title 2")
    F.set_var!(F.LOCAL_VARS, "rss", "A **description** done. Page with tags foo and bar.")
    F.set_var!(F.LOCAL_VARS, "rss_pubdate", Date(2020, 10, 30))
    F.set_var!(F.LOCAL_VARS, "rss_author", "chuck@norris.com")
    F.set_var!(F.LOCAL_VARS, "tags", ["foo", "bar"])

    item, tags = F.add_rss_item()
    @test item.title == "title 2"
    @test item.description == "A description done. Page with tags foo and bar.\n"
    @test item.author == "chuck@norris.com"
    @test item.pubDate == Date(2020, 10, 30)
    @test tags == ["foo", "bar"]
    @test F.RSS_DICT["/hey/ho/index.html"][1].description == item.description
    @test F.RSS_DICT["/hey/ho/index.html"][2] == ["foo", "bar"]

    # Generation
    F.PATHS[:folder] = td
    F.rss_generator()
    ## Global feed
    feed = joinpath(F.PATHS[:site], "feed.xml")
    @test isfile(feed)
    fc = read(feed, String)
    @test occursin("<title><![CDATA[title 2]]></title>", fc)
    @test occursin("<title><![CDATA[title 1]]></title>", fc)
    @test occursin("<description><![CDATA[A description done. Page with tags foo and bar.", fc)
    @test occursin("<description><![CDATA[Page with tag foo.", fc)
    @test occursin("<link>https://github.com/tlienart/Franklin.jl/hey/ho/</link>", fc)
    @test occursin("<link>https://github.com/tlienart/Franklin.jl/hey/hello/</link>", fc)
    @test occursin("<guid>https://github.com/tlienart/Franklin.jl/hey/ho/</guid>", fc)
    @test occursin("<guid>https://github.com/tlienart/Franklin.jl/hey/hello/</guid>", fc)
    @test occursin("
Read more", fc)
    @test occursin("
Read more", fc)
    @test occursin("<pubDate>Fri, 30 Oct 2020 00:00:00 UT</pubDate>", fc)
    @test occursin("<pubDate>Tue, 27 Oct 2020 00:00:00 UT</pubDate>", fc)
    @test occursin("<atom:link href=\"https://github.com/tlienart/Franklin.jl/feed.xml\" rel=\"self\" type=\"application/rss+xml\" />", fc)
    @test findfirst("Fri, 30 Oct 2020", fc) < findfirst("27 Oct 2020", fc) # ordered by pubDate

    ## Tag filtered feeds
    ### foo tag
    feed = joinpath(F.PATHS[:tag], "foo", "feed.xml")
    @test isfile(feed)
    foo_feed = read(feed, String)
    @test occursin("<title><![CDATA[title 2]]></title>", foo_feed)
    @test occursin("<title><![CDATA[title 1]]></title>", foo_feed)
    @test occursin("<description><![CDATA[A description done. Page with tags foo and bar.", foo_feed)
    @test occursin("<description><![CDATA[Page with tag foo.", foo_feed)
    @test occursin("<pubDate>Fri, 30 Oct 2020 00:00:00 UT</pubDate>", foo_feed)
    @test occursin("<pubDate>Tue, 27 Oct 2020 00:00:00 UT</pubDate>", foo_feed)
    @test occursin("<atom:link href=\"https://github.com/tlienart/Franklin.jl/tag/foo/feed.xml\" rel=\"self\" type=\"application/rss+xml\" />", foo_feed)
    @test findfirst("Fri, 30 Oct 2020", foo_feed) < findfirst("27 Oct 2020", foo_feed) # ordered by pubDate
    ### bar tag
    feed = joinpath(F.PATHS[:tag], "bar", "feed.xml")
    @test isfile(feed)
    bar_feed = read(feed, String)
    @test occursin("<title><![CDATA[title 2]]></title>", bar_feed)
    @test !occursin("<title><![CDATA[title 1]]></title>", bar_feed)
    @test occursin("<description><![CDATA[A description done. Page with tags foo and bar.", bar_feed)
    @test !occursin("<description><![CDATA[Page with tag foo.", bar_feed)
    @test occursin("<pubDate>Fri, 30 Oct 2020 00:00:00 UT</pubDate>", bar_feed)
    @test !occursin("<pubDate>Tue, 27 Oct 2020 00:00:00 UT</pubDate>", bar_feed)
    @test occursin("<atom:link href=\"https://github.com/tlienart/Franklin.jl/tag/bar/feed.xml\" rel=\"self\" type=\"application/rss+xml\" />", bar_feed)

    ### t.a g tag (check handling of special characters)
    feed = joinpath(F.PATHS[:tag], "ta_g", "feed.xml")
    @test isfile(feed)
    tag_feed = read(feed, String)
    @test occursin("<atom:link href=\"https://github.com/tlienart/Franklin.jl/tag/ta_g/feed.xml\" rel=\"self\" type=\"application/rss+xml\" />", tag_feed)
end


text:
function baselinetransform(Fₕ, Fₗ=zeros(size(Fₕ)))
    return F -> F |> yeojohnson(Fₕ)
end
export baselinetransform


text:
translate = Dict(
    "black"  => 0,
    "brown"  => 1,
    "red"    => 2,
    "orange" => 3,
    "yellow" => 4,
    "green"  => 5,
    "blue"   => 6,
    "violet" => 7,
    "grey"   => 8,
    "white"  => 9,
)

function label(colors)
    n = [translate[c] for c in colors]
    str = "$(n[1])$(n[2])$(repeat('0', n[3])) "
    replace(str, r"000 " => " kilo") * "ohms"
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
@testset "utils" begin
    #=
    @testset "dist_matrix" begin
        input = [2, 2, 2, 2]
        dist = dist_matrix(input, input; dist_met = (x,y)->abs(x-y))
        @test dist == [[0.0, 0.0, 0.0, 0.0],
                      [0.0, 0.0, 0.0, 0.0],
                      [0.0, 0.0, 0.0, 0.0],
                      [0.0, 0.0, 0.0, 0.0]]
    end
    =#

    @testset "degap" begin
        input = "GGG-G-A-CT"
        degapped = degap(input)
        @test degapped == "GGGGACT"
    end

    @testset "dash_count" begin
        input = "GGG-G-A-CT"
        count = dash_count(input)
        @test count == 3
    end

    @testset "single_gap" begin
        input = "GGG-G-A-CT"
        dash = single_gap(input, 4)
        @test dash == true
    end

    @testset "logsum" begin
        out = logsum(1, 2)
        @test out == 2.313261687518223
    end
end


text:
# This is reference to some of the ideas in Applied Numerical Linear Algebra,
# We are investigating a Lanczos the self-correct the lost of orthogonality
# using ritz vectors projections.
# Here is what we need to do:
#   * Testing for the loss of orthgonality during the runtime of the Lanczos algorithm
#   * Once it's detected, start a Eigensolve on SymTridiagonal and then look for well converged Ritz Vectors.
#   * Store all of those Ritz Vectors as a matrix for reorthogonalizations.

mutable struct IterativeLanczosModified{T<:Number}
    A::Function                         # Linear Operator on Tensor as a function.
    T::Dict{Tuple{Int64, Int64}, T}     # Mapping coordinates [i, j] to elements in the Tridiaognal Matrix.
    k::Int64                            # when k=1, the algorithm haven't started yet.
    Aq::Vector{T}                       # For temp storage to reduce GC time.
    Q::Dict{Int64, Vector{T}}           # A Map for storing all the Lanczos Vectors.
    Q_blocks::Vector{Matrix{T}}         # Blocks of Qs from each restarts
    T_blocks::Vector{Dict}              # Blocks of Ts from each restarts

    tensor_shape::Tuple{Int64}          # Shape of the input for the linear operator.
    V::Vector{Union{Matrix{T}, Vector{T}}}
                                        # Eigen system of T.
    Y::Vector{Vector{T}}                # Well converged ritz vectors from last restart.
    ritz_values::Vector{T}              # The list of ritz values.

    function IterativeLanczosModified(A::Function, q1)
        q1 /= norm(q1)
        Aq = A(q1)
        α = dot(q1, reshape(Aq, :))
        D = typeof(Aq[1])        # Generic Type info
        this = new{D}()          # Referenc to instance created here.
        this.tensor_shape = size(Aq)
        this.Aq = reshape(Aq, :)        # flatten it.
        this.A = A
        this.Q = Dict{Int64, Vector{D}}()
        this.T = Dict{Tuple{Int64}, D}()
        this.T[1, 1] = α
        this.T[0, 1] = this.T[1, 0] = 0
        this.T_blocks = Vector{AbstractArray{D}}()
        this.Q[1] = q1
        this.Q[0] = zeros(D, length(this.Aq))
        this.Q_blocks = Vector{Matrix{D}}()
        this.k = 1
        this.V = Vector{Union{Matrix{D}, Vector{D}}}()
        this.Y = Vector{Vector{D}}()
        this.ritz_values = Vector{D}()
        return this
    end

    function IterativeLanczosModified(A::AbstractMatrix, q1)
        IterativeLanczosModified((x) -> A*x, q1)
    end

end

# Override the () operator.
"""
    It performs one iteration of the Lanczos Orthogonalization.
"""
function (this::IterativeLanczosModified)()
    Aq = this.Aq
    k = this.k
    Q = this.Q
    T = this.T

    q̃ = Aq - T[k - 1, k]*Q[k - 1]
    q̃ -= T[k, k]*Q[k]

    for y in this.Y
        q̃ .-= dot(y, q̃)*y
    end

    β = norm(q̃)
    q = q̃/β

    this.Aq = this.A(q)
    α = dot(q, this.Aq)
    T[k + 1, k + 1] = α
    T[k, k + 1] = T[k + 1, k] = β
    Q[k + 1] = q
    this.k += 1
    return q
end

"""
    Restart and figure out the Ritz Vectors from the SymTridiaognal Matrix,
    get the converged ones and put them into the V vectors of matrices.
"""
function Restart(this::IterativeLanczosModified)
    if this.k == 1                      #  TODO: Investigate Edge Case.
        return
    end

    T = GetPreviousT(this, true)        # T is partial
    Q = GetPreviousQ(this)              # Q is full
    D = typeof(this.Aq[1])
    k = this.k

    this.V = Vector{Union{Matrix{D}, Vector{D}}}()
    this.Y = Vector{Vector{D}}()
    this.ritz_values = Vector{D}()
    # Θ, V = eigen(T, max(k - 8, 1):k - 1)
    Θ, V = eigen(T)
    converged = 0
    for (θ, j) in zip(Θ, 1:length(Θ))
        v = V[:, j]
        error = abs(this.T[k, k - 1]*v[k - 1])
        if error < sqrt(eps(Float32))*norm(T)    # store the ritz values and the ritz vector.
            y = Q[:, 1:end - 1]*v
            ŷ = y/norm(y)
            projVal = 0
            for y in this.Y
                projVal += dot(y, ŷ)
            end

            push!(this.Y, ŷ)
            push!(this.ritz_values, θ)
            converged += 1
        end
    end
    sort!(this.ritz_values)
    return converged
end

"""
    Get the orthogonality error between the newest lanczos vectors
    and the most recent lanczos vector after the last restart.
"""
function GetOrthogonalError(this::IterativeLanczosModified)
    return abs(dot(this.Q[1], this.Q[this.k]))
end

"""
    Get the Q matrix generated from the Iterative Lanczos Process.
"""
function GetQMatrix(this::IterativeLanczosModified)
    return hcat([this.Q_blocks..., GetPreviousQ(this)]...)
end

"""
    Get the Q Matrix since the last restart of the Lanczos Algorithm.
"""
function GetPreviousQ(this::IterativeLanczosModified)
    return hcat([this.Q[idx] for idx in 1: this.k]...)
end

"""
    Get the Symmetric Tridiagonal matrix from the Iterative Lanczos Process.
"""
function GetTMatrix(this::IterativeLanczosModified)
    αs = Vector{Number}()
    βs = Vector{Number}()
    for (T, j) in zip([this.T_blocks..., this.T], 1:length(this.T_blocks) + 1)
        blockSize = maximum(tuple[1] for tuple in keys(T))
        α = [T[idx, idx] for idx in 1:blockSize]
        middle = j > 1 && j <= length(this.T_blocks) + 1
        middleStart = middle ? 0 : 1
        β = [T[idx + 1, idx] for idx in middleStart: blockSize - 1]
        append!(αs, α)
        append!(βs, β)
    end
    return SymTridiagonal(αs, βs)
end

"""
    Get the T matrix since last restart of Lanczos Algorithm.
    parameter:
        `ignore_last`: ignore the last row of matrix T_k.
"""
function GetPreviousT(this::IterativeLanczosModified, ignore_last=false)
    endAt = ignore_last ? this.k - 1 : this.k
    α = [this.T[idx, idx] for idx in 1: endAt]
    β = [this.T[idx + 1, idx] for idx in 1: endAt]
    return SymTridiagonal(α, β)
end

# Basic Testing here -----------------------------------------------------------
using LinearAlgebra, Plots
n = 512
eigenValues = collect(LinRange(1e-3, 1, n)).^2
A = Diagonal(eigenValues);
b = rand(n)
ilm = IterativeLanczosModified(A, b)
converged = 0
for itr in 1: n - 1
    ilm()
    # if itr % 1 == 0
    if GetOrthogonalError(ilm) > 1e-10
        converged = Restart(ilm)
        println("Restarts at itr=$itr, converged = $converged")
        # if converged >= n
        #     break
        # end
    end
end

Q = GetQMatrix(ilm)
T = GetTMatrix(ilm)
println()
heatmap(Q'*Q .|> abs .|> log2)

# LANSO,


text:
module IndependentRandomSequences

using Distributions

import Base: length
import Distributions: _rand!,_logpdf,
                      insupport,
                      mean,var,cov,entropy,
                      mgf,cf

export IIDRandomSequence,INIDRandomSequence,AbstractIDSequence,distributions

include("abstract.jl")
include("IID.jl")
include("INID.jl")

end


text:
@counted function simple1(x::Vector)
    return -x[1] * x[2]
end

@counted function simple1_gradient(x::Vector)
    diff = ForwardDiff.gradient(simple1, x)
    if haskey(COUNTERS, "simple1")
        COUNTERS["simple1"] -= 1
    end
    return diff
end

@counted function simple1_constraints(x::Vector)
    return [x[1] + x[2]^2 - 1,
            -x[1] - x[2]]
end

function simple1_init()
    return rand(2) * 2.0
end

@counted function simple2(x::Vector)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

@counted function simple2_gradient(x::Vector)
    storage = zeros(2)
    storage[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
    storage[2] = 200.0 * (x[2] - x[1]^2)
    return storage
end

@counted function simple2_constraints(x::Vector)
    return [(x[1]-1)^3 - x[2] + 1,
            x[1] + x[2] - 2]
end

function simple2_init()
    return rand(2) .* 2.0 .- 1.0
end

@counted function simple3(x::Vector)
    return x[1] - 2*x[2] + x[3]
end

@counted function simple3_gradient(x::Vector)
    diff = ForwardDiff.gradient(simple3, x)
    if haskey(COUNTERS, "simple3")
        COUNTERS["simple3"] -= 1
    end
    return diff
end

@counted function simple3_constraints(x::Vector)
    return [x[1]^2 + x[2]^2 + x[3]^2 - 1]
end

function simple3_init()
    b = 2.0 .* [1.0, -1.0, 0.0]
    a = -2.0 .* [1.0, -1.0, 0.0]
    return rand(3) .* (b-a) + a
end

const simple_problems = Dict(1=>(f=simple1, g=simple1_gradient, c=simple1_constraints, x0=simple1_init, n=1000),
                             2=>(f=simple2, g=simple2_gradient, c=simple2_constraints, x0=simple2_init, n=1000),
                             3=>(f=simple3, g=simple3_gradient, c=simple3_constraints, x0=simple3_init, n=1000))


text:
using NonlinearEigenproblemsTest
using Test
using NonlinearEigenproblems.RKHelper

@testset "discretizepolygon" begin
@bench @testset "concave polygon" begin
    poly = [0.0; 0+10im; 5+5im; 10+10im; 10+0im]

    expected_boundary_points = [
        0 + 0im, 0 + 2.2071068im, 0 + 4.4142136im,
        0 + 6.6213203im, 0 + 8.8284271im, 0.73223305 + 9.267767im,
        2.2928932 + 7.7071068im, 3.8535534 + 6.1464466im, 5.4142136 + 5.4142136im,
        6.9748737 + 6.9748737im, 8.5355339 + 8.5355339im, 10 + 9.863961im,
        10 + 7.6568542im, 10 + 5.4497475im, 10 + 3.2426407im,
        10 + 1.0355339im, 8.8284271 + 0im, 6.6213203 + 0im,
        4.4142136 + 0im, 2.2071068 + 0im]

    nr_boundary = 20
    nr_interior = 100
    boundary, interior = discretizepolygon(poly, true, nr_boundary, nr_interior)

    # returned result contains boundary points + polygon vertices + first vertex duplicated
    @test length(boundary) == nr_boundary + length(poly) + 1
    @test boundary[1:nr_boundary] ≈ expected_boundary_points
    @test boundary[nr_boundary+1:nr_boundary+length(poly)] == poly
    @test boundary[nr_boundary+length(poly)+1] == poly[1]

    # points on polygon's boundary should count as inside; expand polygon a tiny bit to counteract rounding errors
    t = eps()*5
    expanded_poly = poly + [-t-t*im, -t+t*im, 0+t*im, t+t*im, t-t*im]
    @test all(p -> inpolygon(real(p), imag(p), real.(expanded_poly), imag.(expanded_poly)), boundary)

    @test length(interior) >= nr_interior
    @test all(p -> inpolygon(real(p), imag(p), real.(poly), imag.(poly)), interior)
end

@bench @testset "narrow Σ" begin
    boundary, interior = discretizepolygon([-10.0-2im, 10-2im, 10+2im, -10+2im], true, 100, 5)
    @test length(boundary) == 100 + 5
    @test length(interior) >= 5
end

@bench @testset "too narrow Σ" begin
    @test_throws ErrorException discretizepolygon([-10.0-0.2im, 10-0.2im, 10+0.2im, -10+0.2im], true, 100, 5)
end

@bench @testset "unit disk" begin
    boundary, interior = discretizepolygon(Vector{ComplexF64}(), true, 100, 100)

    @test length(boundary) == 101
    @test all(isapprox.(abs.(boundary[1:100]), 1.0, rtol = 100*eps()))

    @test length(interior) >= 100
    @test all(abs.(interior) .< 1)
end

@bench @testset "Chebyshev points" begin
    p1 = -2.0 - 1im
    p2 = 2.0 + 1im
    boundary, interior = discretizepolygon([p1, p2], true, 100, 100)

    @test length(boundary) == 102
    @test all(imag((boundary .- p1) / (p2-p1)) .== 0)

    @test length(interior) >= 100
    @test all(imag((interior .- p1) / (p2-p1)) .== 0)
end
end

# To visualize, use something like this:
#using Plots
#import GR
#scatter(real(zz), imag(zz))
#scatter(real(Z), imag(Z))


text:
import Base: +, length, size, ==, isempty, show

Lar = LinearAlgebraicRepresentation

#-------------------------------------------------------------------------------
#   MODEL STRUCT DEFINITION
#-------------------------------------------------------------------------------

mutable struct Model
    G::Lar.Points
    T::Array{Lar.ChainOp, 1}

    """
        Model(V::Lar.Points, T::Array{Lar.ChainOp, 1})

    Generic constructor for CAGD.Model.
    Coherency checks are performed between Topology and Geometry
    """
    function Model(V::Lar.Points, T::Array{Lar.ChainOp, 1})
        dim, npts = size(V)
        dim > 0 ||
            throw(ArgumentError("At least one point is needed."))
        length(T) == dim ||
            throw(ArgumentError("Topology is not coherent with Geometry."))

        size(T[1], 2) == npts ||
            throw(ArgumentError("Topology not coherent with Geometry."))
        for i = 2 : dim
            isempty(T[i-1]) || size(T[i], 2) == size(T[i-1],1) ||
                throw(ArgumentError("Topology not coherent with Topology."))
        end

        new(V, T)
    end

    function Model(V::Lar.Points, T::Array{Lar.Cells, 1})
        m = Model(V)
        addModelCells!(m, 1, T[1], signed=true)
        cFE = convert(Lar.ChainOp, Lar.coboundary_1(m.G, T[2], T[1]))
        addModelCells!(m, 2, cFE)
        return m
    end

    """
        Model(V::Lar.Points, EV::Lar.Cells)

    Builds a signed CAGD.Model with vertices and edges only.
    """
    function Model(V::Lar.Points, EV::Lar.Cells)

        I = Array{Int,1}()
        J = Array{Int,1}()
        K = Array{Int8,1}()
        for i = 1 : length(EV)
            let sign = -1;  for j = 1 : length(EV[i])
                push!(I, i)
                push!(J, EV[i][j])
                push!(K, sign)
                sign = 1
            end  end
        end

        T = [SparseArrays.sparse(I, J, K, length(EV), size(V, 2))]
        if size(V, 1) > 1  push!(T, SparseArrays.spzeros(Int8, 0, length(EV)))  end
        Ts = SparseArrays.spzeros(Int8, 0, 0)
        for i = 3 : size(V, 1)  push!(T, Ts)  end
        CAGD.Model(V, T)
    end

    """
        Model(V::Lar.Points)

    Constructor for CAGD.Model with geometry only.
    Topology is set to void by default.
    """
    function Model(V::Lar.Points)
        T = convert(Array{Lar.ChainOp,1},
            [SparseArrays.spzeros(Int8, 0, 0) for i = 1 : size(V, 1)]
        )
        T[1] = convert(Lar.ChainOp, SparseArrays.spzeros(Int8, 0, size(V,2)))
        CAGD.Model(V, T)
    end

    """
        Model()

    Void constructor for CAGD.Model.
    Returns `nothing`
    """
    function Model()
        nothing
    end
end

#-------------------------------------------------------------------------------
#   BASIC PROPERTIES
#-------------------------------------------------------------------------------

length(m::CAGD.Model)               = size(m.G, 1)
size(m::CAGD.Model)		            = size(m.G)
size(m::CAGD.Model, i::Int)         = i == 0 ? size(m.G   ) : size(m.T[i]   )
size(m::CAGD.Model, i::Int, j::Int) = i == 0 ? size(m.G, j) : size(m.T[i], j)
==(m1::CAGD.Model, m2::CAGD.Model)  = m1.G == m2.G && m1.T == m2.T
isempty(m::CAGD.Model, d::Int)      = isempty(m.T[d])
function show(io::IO, m::CAGD.Model)
    println(io, "$(length(m))D model with:")
    println(io, " - $(size(m, 1, 2)) points")
    for d = 1 : length(m)  println(io, " - $(size(m, d, 1)) $d-cells")  end
    display(m.G)
    for d = 1 : length(m)  display(Matrix(m.T[d]))  end
end

#-------------------------------------------------------------------------------
#   OTHER BASE REIMPLEMENTATIONS
#-------------------------------------------------------------------------------

Base.copy(m::CAGD.Model)     = CAGD.Model(m.G, m.T)
Base.deepcopy(m::CAGD.Model) = CAGD.Model(Base.deepcopy(m.G), Base.deepcopy(m.T))

#-------------------------------------------------------------------------------
#   GEOMETRY MANIPULATION
#-------------------------------------------------------------------------------

"""
    addModelVertices!(m::CAGD.Model, V::Lar.Points)::Nothing

Adds `V` vertices to model `m` in tail of `m.G`.
Coherently updates `m.T[1]`.
Points dimension must match.

---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0;0.0 1.0]);

julia> CAGD.addModelVertices!(m, [1.0 1.0; 0.0 1.0])

julia> m
2D model with:
 - 4 points
 - 0 1-cells
 - 0 2-cells
2×4 Array{Float64,2}:
 0.0  0.0  1.0  1.0
 0.0  1.0  0.0  1.0
0×4 Array{Int8,2}
0×0 Array{Int8,2}
```
"""
function addModelVertices!(m::CAGD.Model, V::Lar.Points)::Nothing
    length(m) == size(V, 1) ||
        throw(ArgumentError("Point dimension mismatch."))
    m.G = [m.G V]

    m.T[1] = SparseArrays.sparse(
        findnz(m.T[1])...,
        size(m, 1, 1),
        size(m, 1, 2) + size(V, 2)
    )
    return
end

"""
    addModelVertex!(m::CAGD.Model, v::Array{Float64,1})::Nothing

Adds `v` vertex to model `m` in tail of `m.G`.
Coherently updates `m.T[1]`.
Point dimension must match.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0;0.0 1.0]);

julia> CAGD.addModelVertex!(m, [1.0; 1.0])

julia> m
2D model with:
 - 3 points
 - 0 1-cells
 - 0 2-cells
2×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  1.0  1.0
0×3 Array{Int8,2}
0×0 Array{Int8,2}
```
"""
function addModelVertex!(m::CAGD.Model, v::Array{Float64, 1})::Nothing
    return addModelVertices!(m, v[:, :])
end

"""
    deleteModelVertices!(m::CAGD.Model, vs::Array{Int,1})::Nothing

Delete `vs` indexed vertices from model `m`.
Coherently deletes cells containing such points.
---
# Examples
```jldoctest
julia>  m = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
           SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
           SparseArrays.spzeros(Int8, 0, 3)
       ]);

julia> CAGD.deleteModelVertices!(m, [1; 2])

julia> m
2D model with:
 - 1 points
 - 0 1-cells
 - 0 2-cells
2×1 Array{Float64,2}:
 1.0
 0.0
0×1 Array{Int8,2}
0×3 Array{Int8,2}
```
"""
function deleteModelVertices!(m::CAGD.Model, vs::Array{Int,1})::Nothing
    min(vs...) > 0 || throw(ArgumentError("Vertex indices are positive integers"))
    max(vs...) ≤ size(m, 0, 2) || throw(ArgumentError("There are not $(max(vs...)) vertices"))

    tokeepV = setdiff(collect(1 : size(m, 0, 2)), vs)
    m.G     = m.G[:, tokeepV]

    if !isempty(m.T[1])
        # From EV must be deleted all the columns linked to a deleted vertex
        #  and all rows related to a dangling edge
        todelE = abs.(m.T[1][:, vs])
        todelE = [l for l = 1 : todelE.m if sum(todelE[l, :]) != 0]
        m.T[1] = m.T[1][:, tokeepV]
        if !isempty(todelE) CAGD.deleteModelCells!(m, 1, todelE) end
    end

    return
end

"""
    deleteModelVertex!(m::CAGD.Model, v::Int)::Nothing

Delete `v` indexed vertex from model `m`.
Coherently deletes cells containing such point.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
    SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
    SparseArrays.spzeros(Int8, 0, 3)
]);

julia> CAGD.deleteModelVertex!(m, 1)

julia> m
2D model with:
 - 2 points
 - 1 1-cells
 - 0 2-cells
2×2 Array{Float64,2}:
 0.0  1.0
 1.0  0.0
1×2 Array{Int8,2}:
 -1  1
0×3 Array{Int8,2}
```
"""
function deleteModelVertex!(m::CAGD.Model, v::Int)::Nothing
    deleteModelVertices!(m, [v])
end

#-------------------------------------------------------------------------------
#   TOPOLOGY MANIPULATION
#-------------------------------------------------------------------------------

"""
    addModelCells!(m::CAGD.Model, deg::Int, cs::Lar.ChainOp)::Nothing

Adds to model `m` the `deg` cells specified by `cs` given as sparse matrix.
Coherently updates `deg+1` cells.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0]);

julia> CAGD.addModelCells!(m, 1, SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]))

julia> m
2D model with:
 - 3 points
 - 3 1-cells
 - 0 2-cells
2×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  1.0  0.0
3×3 Array{Int8,2}:
 -1   1  0
 -1   0  1
  0  -1  1
0×3 Array{Int8,2}
```
"""
function addModelCells!(m::CAGD.Model, deg::Int, cs::Lar.ChainOp)::Nothing
    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    (cs.n == size(m, deg - 1, 1) && deg > 1) ||
        (cs.n == size(m, deg - 1, 2) && deg == 1) ||
        throw(ArgumentError("Incoherent Chain Operator size"))

    # Adding input cells to the topological structure
    m.T[deg] = [m.T[deg]; cs]

    # Adding slot for the new cells to the higher order cells
    if deg < length(m)
        m.T[deg + 1] = [
            m.T[deg + 1] SparseArrays.spzeros(Int8, size(m, deg + 1, 1), cs.m)
        ]
    end

    return
end

"""
    addModelCell!(m::CAGD.Model, deg::Int, c::Lar.Cell)::Nothing

Adds to model `m` the `deg` cell specified by `c` as sparse vector.
Coherently updates `deg+1` cells.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0]);

julia> CAGD.addModelCell!(m, 1, SparseArrays.sparse(Int8[-1; 1; 0]))

julia> m
2D model with:
 - 3 points
 - 1 1-cells
 - 0 2-cells
2×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  1.0  0.0
1×3 Array{Int8,2}:
 -1  1  0
0×1 Array{Int8,2}
```
"""
function addModelCell!(m::CAGD.Model, deg::Int, c::Lar.Cell)::Nothing
    CAGD.addModelCells!(m, deg, convert(Lar.ChainOp, c'))
end

"""
    addModelCells!(m::CAGD.Model, deg::Int, cs::Lar.Cells; signed=false)::Nothing

Adds to model `m` the `deg` cells specified by `cs` as array of arrays.
Coherently updates `deg+1` cells.
If `signed` is set to `true`, cells are evaluated with sign.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0]);

julia> CAGD.addModelCells!(m, 1, [[1, 2], [1, 3], [2, 3]])

julia> CAGD.addModelCells!(m, 1, [[1, 2], [1, 3], [2, 3]], signed = true)

julia> m
2D model with:
 - 3 points
 - 6 1-cells
 - 0 2-cells
2×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  1.0  0.0
6×3 Array{Int8,2}:
  1   1  0
  1   0  1
  0   1  1
 -1   1  0
 -1   0  1
  0  -1  1
0×6 Array{Int8,2}
```
"""
function addModelCells!(
        m::CAGD.Model, deg::Int, cs::Lar.Cells; signed=false
    )::Nothing

    if signed
        if deg == 1
            I = Array{Int,1}()
            J = Array{Int,1}()
            K = Array{Int8,1}()
            for i = 1 : length(cs)
                let sign = -1;  for j = 1 : length(cs[i])
                    push!(I, i)
                    push!(J, cs[i][j])
                    push!(K, sign)
                    sign = 1
                end  end
            end

            scs = SparseArrays.sparse(I, J, K, length(cs), size(m, deg, 2))
        elseif deg == 2
            EV = Lar.cop2lar(m.T[1]);
            FV = [collect(Set(vcat([EV[e] for e in face]...))) for face in cs]
            scs = convert(Lar.ChainOp, Lar.coboundary_1(m.G, FV, EV))
        else
            println("No Methods for $deg-signed-coboundary.")
            println("Using unsigned coboundary instead.")
            signed = false
        end
    end

    if !signed
        I = Array{Int,1}()
        J = Array{Int,1}()
        K = Array{Int8,1}()
        for i = 1 : length(cs)
            for j = 1 : length(cs[i])
                push!(I, i)
                push!(J, cs[i][j])
                push!(K, 1)
            end
        end

        scs = SparseArrays.sparse(I, J, K, length(cs), size(m, deg, 2))
    end

    return CAGD.addModelCells!(m, deg, scs)
end

"""
    addModelCell!(m::CAGD.Model, deg::Int, c::Array{Int64,1}; signed=false)::Nothing

Adds to model `m` the `deg` cell specified by `c` as index of `deg-1` cells.
Coherently updates `deg+1` cells.
If `signed` is set to `true`, the cell is computed with sign.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0]);

julia> CAGD.addModelCell!(m, 1, [1, 2])

julia> CAGD.addModelCell!(m, 1, [1, 3], signed = true)

julia> m
2D model with:
 - 3 points
 - 2 1-cells
 - 0 2-cells
2×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  1.0  0.0
2×3 Array{Int8,2}:
  1  1  0
 -1  0  1
0×2 Array{Int8,2}
```
"""
function addModelCell!(m::CAGD.Model, deg::Int, c::Array{Int64,1}; signed=false)::Nothing
    CAGD.addModelCells!(m, deg, [c], signed = signed)
end

"""
    deleteModelCells!(m::CAGD.Model, deg::Int, cs::Array{Int, 1})::Nothing

Delete `cs` indexed `deg`-cells from model `m`.
Coherently deletes higher order cells containing such cells.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
    SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
    SparseArrays.sparse(Int8[1 -1 1])
]);

julia> CAGD.deleteModelCells!(m, 1, [1; 2])

julia> m
2D model with:
 - 3 points
 - 1 1-cells
 - 0 2-cells
2×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  1.0  0.0
1×3 Array{Int8,2}:
 0  -1  1
0×1 Array{Int8,2}
```
"""
function deleteModelCells!(m::CAGD.Model, deg::Int, cs::Array{Int, 1})::Nothing
    !isempty(cs) || return
    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    !isempty(m.T[deg]) ||
        throw(ArgumentError("There are no cells of degree $deg"))
    max(cs...) ≤ size(m, deg, 1) ||
        throw(ArgumentError("There are not $(max(cs...)) cells of degree $deg"))

    # Removing `cs` rows from `m.T[deg]`
    tokeep     = setdiff(collect(1 : size(m, deg, 1)), cs)
    m.T[deg]   = m.T[deg][tokeep, :]

    if  deg == length(m) return  end

    # Removing `cs` cols from `m.T[deg+1]` by checking if some `deg+1` cell
    #  has to be removed too.
    # An higher ord cell must be removed if it contains a removed lower ord cell
    todelHo    = m.T[deg + 1][:, cs]
    m.T[deg+1] = m.T[deg + 1][:, tokeep]
    if isempty(m.T[deg+1]) return end
    todelHo, _, _ = SparseArrays.findnz(todelHo)
    if !isempty(todelHo) CAGD.deleteModelCells!(m, deg+1, unique(todelHo)) end
    return
end

"""
    deleteModelCells!(m::CAGD.Model, deg::Int, c::Int)::Nothing

Delete `c` indexed `deg`-cell from model `m`.
Coherently deletes higher order cells containing such cell.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
    SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
    SparseArrays.sparse(Int8[1 -1 1])
]);

julia> CAGD.deleteModelCell!(m, 1, 1)

julia> m
2D model with:
 - 3 points
 - 2 1-cells
 - 0 2-cells
2×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  1.0  0.0
2×3 Array{Int8,2}:
 -1   0  1
  0  -1  1
0×2 Array{Int8,2}
```
"""
function deleteModelCell!(m::CAGD.Model, deg::Int, c::Int)::Nothing
    CAGD.deleteModelCells!(m, deg, [c])
end

"""
    getModelLoCell(m::CAGD.Model, deg::Int, c::Int)::Array{Int, 1}
    getModelLoCell(m::CAGD.Model, deg::Int, cs::Array{Int,1})::Array{Int, 1}

Determines the indices of `deg-1` cells forming the boundary of `c`/`cs` indexed `deg` cell(s) of `m`
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0], [
    SparseArrays.sparse(Int8[-1 1 0 0; -1 0 1 0; 0 -1 1 0; 0 -1 0 1; 0 0 -1 1]),
    SparseArrays.sparse(Int8[1 -1 1 0 0; 0 0 1 -1 1])
])
2D model with:
 - 4 points
 - 5 1-cells
 - 2 2-cells
2×4 Array{Float64,2}:
 0.0  0.0  1.0  1.0
 0.0  1.0  0.0  1.0
5×4 Array{Int8,2}:
 -1   1   0  0
 -1   0   1  0
  0  -1   1  0
  0  -1   0  1
  0   0  -1  1
2×5 Array{Int8,2}:
 1  -1  1   0  0
 0   0  1  -1  1

julia> CAGD.getModelLoCell(m, 2, 2)
3-element Array{Int64,1}:
 3
 4
 5

julia> CAGD.getModelLoCell(m, 1, [3, 4, 5])
3-element Array{Int64,1}:
 2
 3
 4
```
"""
function getModelLoCell(m::CAGD.Model, deg::Int, c::Int)::Array{Int, 1}
    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    c > 0 || throw(ArgumentError("Cell indices are positive integers"))
    c ≤ size(m, deg, 1) || throw(ArgumentError("There are not $c $deg-cells"))
    return m.T[deg][c, :].nzind
end

function getModelLoCell(m::CAGD.Model, deg::Int, cs::Array{Int, 1})::Array{Int, 1}
    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    min(cs...) > 0 || throw(ArgumentError("Cell indices are positive integers"))
    max(cs...) ≤ size(m, deg, 1) || throw(ArgumentError("There are not $(max(cs...)) $deg-cells"))
    return ∪([m.T[deg][c, :].nzind for c in cs]...)
end

"""
    getModelCellVertices(m::CAGD.Model, deg::Int, c::Int, ret_idx=false)
    getModelCellVertices(m::CAGD.Model, deg::Int, cs::Array{Int,1}, ret_idx=false)
        ::Union{Array{Array{Float64,1},1}, Tuple{Array{Array{Float64,1},1}, Array{Int,1}}}

Determines the list of vertices of `c`/`cs` indexed `deg` cell(s) of `m`.
If `ret_idx` is set true, it also returns the vertices list of index.

See also [`CAGD.getModelCellGeometry`](@ref) to get vertices as `Lar.Points`.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0], [
    SparseArrays.sparse(Int8[-1 1 0 0; -1 0 1 0; 0 -1 1 0; 0 -1 0 1; 0 0 -1 1]),
    SparseArrays.sparse(Int8[1 -1 1 0 0; 0 0 1 -1 1])
])
2D model with:
 - 4 points
 - 5 1-cells
 - 2 2-cells
2×4 Array{Float64,2}:
 0.0  0.0  1.0  1.0
 0.0  1.0  0.0  1.0
5×4 Array{Int8,2}:
 -1   1   0  0
 -1   0   1  0
  0  -1   1  0
  0  -1   0  1
  0   0  -1  1
2×5 Array{Int8,2}:
 1  -1  1   0  0
 0   0  1  -1  1

julia> CAGD.getModelCellVertices(m, 2, 2)
3-element Array{Array{Float64,1},1}:
 [0.0, 1.0]
 [1.0, 0.0]
 [1.0, 1.0]

julia> CAGD.getModelCellVertices(m, 1, [3, 4, 5], true)
(Array{Float64,1}[[0.0, 1.0], [1.0, 0.0], [1.0, 1.0]], [2, 3, 4])
```
"""
function getModelCellVertices(
        m::CAGD.Model, deg::Int, cs::Array{Int,1}, ret_idx=false
    )::Union{Array{Array{Float64,1},1}, Tuple{Array{Array{Float64,1},1}, Array{Int,1}}}

    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    min(cs...) > 0 || throw(ArgumentError("Cell indices are positive integers"))
    max(cs...) ≤ size(m, deg, 1) || throw(ArgumentError("There are not $(max(cs...)) $deg-cells"))

    for d = deg : -1 : 1
        cs = ∪([m.T[d][el, :].nzind for el in cs]...)
    end

    if ret_idx == true
        return (map(i -> m.G[:, cs[i]], 1:length(cs)), cs)
    end
    return map(i -> m.G[:, cs[i]], 1:length(cs))
end

function getModelCellVertices(
        m::CAGD.Model, deg::Int, c::Int, ret_idx=false
    )::Union{Array{Array{Float64,1},1}, Tuple{Array{Array{Float64,1},1}, Array{Int,1}}}
    return getModelCellVertices(m, deg, [c], ret_idx)
end

"""
    getModelCellGeometry(m::CAGD.Model, deg::Int, c::Int, ret_idx=false)
    getModelCellGeometry(m::CAGD.Model, deg::Int, cs::Array{Int,1}, ret_idx=false)
        ::Union{Lar.Points, Tuple{Lar.Points, Array{Int,1}}}

Determines the geometry of `c`/`cs` indexed `deg` cell(s) of `m`.
If `ret_idx` is set true, it also returns the vertices list of index.

See also [`CAGD.getModelCellVertices`](@ref) to get vertices as a list.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0], [
    SparseArrays.sparse(Int8[-1 1 0 0; -1 0 1 0; 0 -1 1 0; 0 -1 0 1; 0 0 -1 1]),
    SparseArrays.sparse(Int8[1 -1 1 0 0; 0 0 1 -1 1])
])
2D model with:
 - 4 points
 - 5 1-cells
 - 2 2-cells
2×4 Array{Float64,2}:
 0.0  0.0  1.0  1.0
 0.0  1.0  0.0  1.0
5×4 Array{Int8,2}:
 -1   1   0  0
 -1   0   1  0
  0  -1   1  0
  0  -1   0  1
  0   0  -1  1
2×5 Array{Int8,2}:
 1  -1  1   0  0
 0   0  1  -1  1

julia> CAGD.getModelCellGeometry(m, 2, 2)
2×3 Array{Float64,2}:
 0.0  1.0  1.0
 1.0  0.0  1.0

julia> CAGD.getModelCellGeometry(m, 1, [3, 4, 5], true)
([0.0 1.0 1.0; 1.0 0.0 1.0], [2, 3, 4])
```
"""
function getModelCellGeometry(
        m::CAGD.Model, deg::Int, cs::Array{Int,1}, ret_idx=false
    )::Union{Lar.Points, Tuple{Lar.Points, Array{Int,1}}}

    deg > 0 || throw(ArgumentError("Degree must be a non negative value"))
    deg ≤ length(m) || throw(ArgumentError("The model do not have degree $deg"))
    min(cs...) > 0 || throw(ArgumentError("Cell indices are positive integers"))
    max(cs...) ≤ size(m, deg, 1) || throw(ArgumentError("There are not $(max(cs...)) $deg-cells"))
    for d = deg : -1 : 1
        cs = ∪([m.T[d][el, :].nzind for el in cs]...)
    end

    if ret_idx == true
        return (m.G[:, cs], cs)
    end
    return m.G[:, cs]
end

function getModelCellGeometry(
        m::CAGD.Model, deg::Int, c::Int, ret_idx=false
    )::Union{Lar.Points, Tuple{Lar.Points, Array{Int,1}}}

    return getModelCellGeometry(m, deg, [c], ret_idx)
end

#-------------------------------------------------------------------------------
#   MODEL MANIPULATION
#-------------------------------------------------------------------------------

"""
    modelPurge!(m::CAGD.Model[; depth::Int = 0])::Nothing

Purges the model Geometry from unused elements.

Purges each Geometry element that is not related to a topological Cell.
If the optional argument `depth` is specified than it also purges all the
topological cells of degree lower or equal to `depth` that are not related
to an higher order topological cell.
---
# Examples
```jldoctest
julia> m = CAGD.Model([0.0 0.0 1.0 0.0 2.0; 0.0 1.0 0.0 2.0 0.0], [
    SparseArrays.sparse(Int8[-1 1 0 0 0; -1 0 1 0 0; 0 -1 1 0 0; 0 -1 0 1 0]),
    SparseArrays.sparse(Int8[1 -1 1 0 ])
])
2D model with:
 - 5 points
 - 4 1-cells
 - 1 2-cells
2×5 Array{Float64,2}:
 0.0  0.0  1.0  0.0  2.0
 0.0  1.0  0.0  2.0  0.0
4×5 Array{Int8,2}:
 -1   1  0  0  0
 -1   0  1  0  0
  0  -1  1  0  0
  0  -1  0  1  0
1×4 Array{Int8,2}:
 1  -1  1  0

julia> CAGD.modelPurge!(m)

julia> m
2D model with:
 - 4 points
 - 4 1-cells
 - 1 2-cells
2×4 Array{Float64,2}:
 0.0  0.0  1.0  0.0
 0.0  1.0  0.0  2.0
4×4 Array{Int8,2}:
 -1   1  0  0
 -1   0  1  0
  0  -1  1  0
  0  -1  0  1
1×4 Array{Int8,2}:
 1  -1  1  0

julia> CAGD.modelPurge!(m, depth = 1)

julia>m
2D model with:
 - 3 points
 - 3 1-cells
 - 1 2-cells
2×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  1.0  0.0
3×3 Array{Int8,2}:
 -1   1  0
 -1   0  1
  0  -1  1
1×3 Array{Int8,2}:
 1  -1  1
```
"""
function modelPurge!(m::CAGD.Model; depth::Int = 0)::Nothing
    depth < length(m) ||
        throw(ArgumentError("The model do not have degree $(depth+1)"))

    for d = depth : -1 : 1
        todel = [i for i = 1 : size(m, d+1, 2) if sum(abs.(m.T[d+1][:,i])) == 0]
        if !isempty(todel) CAGD.deleteModelCells!(m, d, todel) end
    end
    todel = [i for i = 1 : size(m, 1, 2) if sum(abs.(m.T[1][:, i])) == 0]
    if !isempty(todel) CAGD.deleteModelVertices!(m, todel) end
    return
end

"""
    mergeModelVertices(m::CAGD.Model[; err=1e-6[, signed_merge = false]])::CAGD.Model

Evaluates Geometry congruence with atol `err` and coherently updates Topology via CCE algorithm of `m`.

The method regroups vertices closer than `err` in a single vertex via KD-tree.
Then it coherently updates topology cells, eventually mashing them up if needed.
A topological gift wrapping (see [`CAGD.tgw`](@ref)) might be needed after the application,
if cells merged loosing in resolution.

If `signed_merge` is set to true, than `m` is expected signed and it returns a signed model.

See also:
 - [`CAGD.vcongruence`](@ref)
 - [`CAGD.cellcongruence`](@ref)
 - [`CAGD.signedCellCongruence`](@ref)
---
# Examples
```jldoctest

```
"""
function mergeModelVertices(m::CAGD.Model; err=1e-6, signed_merge = false)
    V, cls = CAGD.vcongruence(m.G, epsilon = err)
    if signed_merge  lo_sign = [ones(Int8, length(cl)) for cl in cls]  end
    congModel = CAGD.Model(V)
    for d = 1 : length(m)
        if isempty(m.T[d])  break  end
        if signed_merge
            cop, cls, lo_sign = CAGD.signedCellCongruence(m.T[d], cls, lo_sign, dim=d)
            if  d > 1 cop = -cop  end
        else
            cop, cls = CAGD.cellcongruence(m.T[d], cls, dim=d)
            cop = convert(Lar.Cells, cop)
        end
        CAGD.addModelCells!(congModel, d, cop)
    end

    return congModel
end

#-------------------------------------------------------------------------------
#   MODEL INTERACTION
#-------------------------------------------------------------------------------

"""
    uniteModels!(model::CAGD.Model, m2::CAGD.Model)::Nothing

Fuses `m2` in `model` by adding coherently geometry and topology.
---
# Examples
```jldoctest
julia> m1 = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
    SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
    SparseArrays.sparse(Int8[1 -1 1])
]);

julia> m2 = CAGD.Model([0.0 0.0 -1.0; 0.0 -1.0 0.0], [
    SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
    SparseArrays.sparse(Int8[1 -1 1])
]);

julia> CAGD.uniteModels!(m1, m2)

julia> m1
2D model with:
 - 6 points
 - 6 1-cells
 - 2 2-cells
2×6 Array{Float64,2}:
 0.0  0.0  1.0  0.0   0.0  -1.0
 0.0  1.0  0.0  0.0  -1.0   0.0
6×6 Array{Int8,2}:
 -1   1  0   0   0  0
 -1   0  1   0   0  0
  0  -1  1   0   0  0
  0   0  0  -1   1  0
  0   0  0  -1   0  1
  0   0  0   0  -1  1
2×6 Array{Int8,2}:
 1  -1  1  0   0  0
 0   0  0  1  -1  1
```
"""
function uniteModels!(model::CAGD.Model, m2::CAGD.Model)::Nothing

    length(model) == length(m2) ||
        throw(ArgumentError("ERROR: Inconsistent models dimension!"))

    CAGD.addModelVertices!(model, m2.G)
    for d = 1 : length(model)
        shift = size(model, d, 2) - size(m2, d, 2)
        I, J, K = SparseArrays.findnz(m2.T[d])
        J = J .+ shift
        CAGD.addModelCells!(model, d, SparseArrays.sparse(
            I, J, K, size(m2, d, 1), size(m2, d, 2) + shift
        ))
    end

    return
end

"""
    uniteModels(m1::CAGD.Model, m2::CAGD.Model)::CAGD.Model

Fuses `m1` and `m2` in a new model by adding coherently geometry and topology.
---
# Examples
```jldoctest
julia> m1 = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
    SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
    SparseArrays.sparse(Int8[1 -1 1])
]);

julia> m2 = CAGD.Model([0.0 0.0 -1.0; 0.0 -1.0 0.0], [
    SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
    SparseArrays.sparse(Int8[1 -1 1])
]);

julia> CAGD.uniteModels(m1, m2)
2D model with:
 - 6 points
 - 6 1-cells
 - 2 2-cells
2×6 Array{Float64,2}:
 0.0  0.0  1.0  0.0   0.0  -1.0
 0.0  1.0  0.0  0.0  -1.0   0.0
6×6 Array{Int8,2}:
 -1   1  0   0   0  0
 -1   0  1   0   0  0
  0  -1  1   0   0  0
  0   0  0  -1   1  0
  0   0  0  -1   0  1
  0   0  0   0  -1  1
2×6 Array{Int8,2}:
 1  -1  1  0   0  0
 0   0  0  1  -1  1
```
"""
function uniteModels(m1::CAGD.Model, m2::CAGD.Model)::CAGD.Model
    model = CAGD.deepcopy(m1)
    CAGD.uniteModels!(model, m2)
    return model
end

"""
    uniteMultipleModels(models::Array{CAGD.Model,1})::CAGD.Model

Fuses all the `models` in a new model by adding coherently Geometry and Topology
---
# Examples
```jldoctest
julia> models = [
    CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
        SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
        SparseArrays.sparse(Int8[1 -1 1])
    ]);

    CAGD.Model([0.0 0.0 -1.0; 0.0 -1.0 0.0], [
        SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
        SparseArrays.sparse(Int8[1 -1 1])
    ]);

    CAGD.Model([0.0 -1.0 -1.0; 1.0 0.0 1.0], [
        SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
        SparseArrays.sparse(Int8[1 -1 1])
    ])
];

julia> CAGD.uniteMultipleModels(models)
2D model with:
 - 9 points
 - 9 1-cells
 - 3 2-cells
2×9 Array{Float64,2}:
 0.0  0.0  1.0  0.0   0.0  -1.0  0.0  -1.0  -1.0
 0.0  1.0  0.0  0.0  -1.0   0.0  1.0   0.0   1.0
9×9 Array{Int8,2}:
 -1   1  0   0   0  0   0   0  0
 -1   0  1   0   0  0   0   0  0
  0  -1  1   0   0  0   0   0  0
  0   0  0  -1   1  0   0   0  0
  0   0  0  -1   0  1   0   0  0
  0   0  0   0  -1  1   0   0  0
  0   0  0   0   0  0  -1   1  0
  0   0  0   0   0  0  -1   0  1
  0   0  0   0   0  0   0  -1  1
3×9 Array{Int8,2}:
 1  -1  1  0   0  0  0   0  0
 0   0  0  1  -1  1  0   0  0
 0   0  0  0   0  0  1  -1  1
```
"""
function uniteMultipleModels(models::Array{CAGD.Model,1})::CAGD.Model

    model = deepcopy(models[1])

    for i = 2 : length(models)
        CAGD.uniteModels!(model, models[i])
    end

    return model
end

"""
    mergeMultipleModels(models::Array{CAGD.Model,1}[; err=1e-6[, signed_merge=false]])::CAGD.Model

Unite `models` and performs CCE algorithm of the result with atol `err`.

If `signed_merge` is set true, `models` are expected signed and it returns a signed model.

See also:
 - [`CAGD.uniteMultipleModels`](@ref) for model unions;
 - [`CAGD.mergeModelVertices`](@ref) for CEE algorithm.
---
# Examples
```jldoctest
julia> models = [
    CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
        SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
        SparseArrays.sparse(Int8[1 -1 1])
    ]);

    CAGD.Model([0.0 0.0 -1.0; 0.0 -1.0 0.0], [
        SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
        SparseArrays.sparse(Int8[1 -1 1])
    ]);

    CAGD.Model([0.0 -1.0 -1.0; 1.0 0.0 1.0], [
        SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
        SparseArrays.sparse(Int8[1 -1 1])
    ])
];

julia> m = CAGD.mergeMultipleModels(models, signed_merge = true)
2D model with:
 - 6 points
 - 9 1-cells
 - 3 2-cells
2×6 Array{Float64,2}:
 0.0  0.0  1.0   0.0  -1.0  -1.0
 0.0  1.0  0.0  -1.0   0.0   1.0
9×6 Array{Int8,2}:
 -1   0  0   0   1  0
 -1   0  0   1   0  0
 -1   0  1   0   0  0
 -1   1  0   0   0  0
  0  -1  0   0   0  1
  0  -1  0   0   1  0
  0  -1  1   0   0  0
  0   0  0  -1   1  0
  0   0  0   0  -1  1
3×9 Array{Int8,2}:
 1  -1  0   0  0   0   0  -1   0
 0   0  1  -1  0   0  -1   0   0
 0   0  0   0  1  -1   0   0  -1

julia> Matrix(tgw(m, 2)[1])
4×9 Array{Int8,2}:
  1  -1   0   0   0   0   0  -1   0
  0   0   1  -1   0   0  -1   0   0
  0   0   0   0   1  -1   0   0  -1
 -1   0   0   1   0   1   0   0   0
  0   1  -1   0  -1   0   1   1   1
```
"""
function mergeMultipleModels(models::Array{CAGD.Model,1}; err=1e-6, signed_merge=false)::CAGD.Model
    return CAGD.mergeModelVertices(CAGD.uniteMultipleModels(models), err=err, signed_merge=signed_merge)
end

#-------------------------------------------------------------------------------
#   MODEL EXPORT
#-------------------------------------------------------------------------------

"""
    jlexportModel(model::CAGD.Model, filename::String[, modelname::String = "model"])::Nothing

Exports `model` by creating a file named `filename` that, if executed,
reloads the model with identifier `modelname` (wich is model, by default).
---
# Examples
```jldoctest
julia> modelOut = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
    SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
    SparseArrays.sparse(Int8[1 -1 1])
]);

julia> CAGD.jlexportModel(modelOut, "test.jl", "modelIn")

julia> include("./test.jl");

julia> modelIn
2D model with:
 - 3 points
 - 3 1-cells
 - 1 2-cells
2×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  1.0  0.0
3×3 Array{Int8,2}:
 -1   1  0
 -1   0  1
  0  -1  1
1×3 Array{Int8,2}:
 1  -1  1

julia> rm("./testModelExport.jl")
```
"""
function jlexportModel(
        model::CAGD.Model,
        filename::String,
        modelname::String = "model"
    )::Nothing

    open(filename, "w") do f
        # Write imports
        write(f, "# using LinearAlgebraicRepresentation\n")
        write(f, "using SparseArrays\n")
        write(f, "# using CAGD\n\n")

        # Write Geometry
        write(f, "#=== GEOMETRY ===#\n")
        write(f, "$modelname = CAGD.Model([")
        for d = 1 : length(model)
            write(f, "\n\t")
            for p in model.G[d, :]  write(f, "$p ")  end
        end
        write(f, "\n]);\n\n")

        # Write Topology
        for d = 1 : length(model)
            write(f, "#=== TOPOLOGY : $d-CELLS ===#\n")
            I, J, K = findnz(model.T[d])
            write(f, "I = Array{Int64}([")
            for i in I  write(f, "$i, ")  end
            write(f, "]);\nJ = Array{Int64}([")
            for j in J  write(f, "$j, ")  end
            write(f, "]);\nK = Array{Int8}([")
            for k in K  write(f, "$k, ")  end
            write(f, "]);\n\n")
            write(f, "CAGD.addModelCells!($modelname, $d, SparseArrays.sparse(")
            write(f, "I, J, K, $(size(model, d, 1)), $(size(model, d, 2))")
            write(f, "))\n\n\n")
        end
    end

    return
end


text:
## Great

using ApproxManifoldProducts
using Test
using LinearAlgebra

using KernelDensityEstimate

# pts = [-pi+0.15; -pi+0.2; +pi-0.2; pi-0.149]
#
# pcm = kde!(pts, [0.05], (addtheta,), (difftheta,) )
#
#
#
# pcm.bt.centers

## but wait

# TODO ensure this manifold product is still valid
# @testset "basic tests..." begin
#
#
# mbe1 = ManifoldBelief(EuclideanManifold, 0.0)
# mbe2 = ManifoldBelief(EuclideanManifold, 0.0)
#
# *([mbe1;mbe2])
#
# end

##

@testset "test circular difftheta..." begin

th1 = -0.1
th2 = +0.1

# should result in +0.2
@test -1e-10 < difftheta(th2, th1) - 0.2 < 1e-10
# should result in -0.2
@test -1e-10 < difftheta(th1, th2) + 0.2 < 1e-10

th1 = +0.1
th2 = -0.1

# should result in -0.2
@test -1e-10 < difftheta(th2, th1) + 0.2 < 1e-10
# should result in +0.2
@test -1e-10 < difftheta(th1, th2) - 0.2 < 1e-10

th1 = pi-0.1
th2 = -pi+0.1

# should result in +0.2
@test -1e-10 < difftheta(th2, th1) - 0.2 < 1e-10
# should result in -0.2
@test -1e-10 < difftheta(th1, th2) + 0.2 < 1e-10

th1 = pi-0.1 + 2pi
th2 = -pi+0.1 +2pi

# should result in +0.2
@test -1e-10 < difftheta(th2, th1) - 0.2 < 1e-10
# should result in -0.2
@test -1e-10 < difftheta(th1, th2) + 0.2 < 1e-10

th1 = pi-0.1  -2pi
th2 = -pi+0.1 -2pi

# should result in +0.2
@test -1e-10 < difftheta(th2, th1) - 0.2 < 1e-10
# should result in -0.2
@test -1e-10 < difftheta(th1, th2) + 0.2 < 1e-10

th1 = pi-0.1  +2pi
th2 = -pi+0.1 -2pi

# should result in +0.2
@test -1e-10 < difftheta(th2, th1) - 0.2 < 1e-10
# should result in -0.2
@test -1e-10 < difftheta(th1, th2) + 0.2 < 1e-10

th1 = pi-0.1  -2pi
th2 = -pi+0.1 +2pi

# should result in +0.2
@test -1e-10 < difftheta(th2, th1) - 0.2 < 1e-10
# should result in -0.2
@test -1e-10 < difftheta(th1, th2) + 0.2 < 1e-10

th1 = -0.1 + pi/2.0
th2 = +0.1 + pi/2.0

# should result in +0.2
@test -1e-10 < difftheta(th2, th1) - 0.2 < 1e-10
# should result in -0.2
@test -1e-10 < difftheta(th1, th2) + 0.2 < 1e-10

eps_l = 1e-5
th1 = 0.0
th2 = -pi + eps_l

@test 0.0 < difftheta(th2, th1) + pi < 2*eps_l
@test -2*eps_l < difftheta(th1, th2) - pi < 0.0

end

##

@testset "test circular manifold BallTree construction small positive..." begin

pts = [0.15; 0.2; 0.35; 0.4]

p = kde!(pts, [0.05])

pcp = kde!(pts, [0.05], (addtheta,) )

pcm = kde!(pts, [0.05], (addtheta,), (difftheta,) )

@test norm(p.bt.centers - pcp.bt.centers) < 1e-10
@test norm(p.bt.centers - pcm.bt.centers) < 1e-10

@test norm(p.bt.ranges - pcp.bt.ranges) < 1e-10
@test norm(p.bt.ranges - pcm.bt.ranges) < 1e-10

@test norm(p.bt.weights - pcp.bt.weights) < 1e-10
@test norm(p.bt.weights - pcm.bt.weights) < 1e-10

@test norm(p.bt.left_child - pcp.bt.left_child) < 1e-10
@test norm(p.bt.left_child - pcm.bt.left_child) < 1e-10

@test norm(p.bt.right_child - pcp.bt.right_child) < 1e-10
@test norm(p.bt.right_child - pcm.bt.right_child) < 1e-10

@test norm(p.bt.lowest_leaf - pcp.bt.lowest_leaf) < 1e-10
@test norm(p.bt.lowest_leaf - pcm.bt.lowest_leaf) < 1e-10

@test norm(p.bt.highest_leaf - pcp.bt.highest_leaf) < 1e-10
@test norm(p.bt.highest_leaf - pcm.bt.highest_leaf) < 1e-10

@test norm(p.bt.permutation - pcp.bt.permutation) < 1e-10
@test norm(p.bt.permutation - pcm.bt.permutation) < 1e-10

end

##

@testset "test circular manifold BallTree construction small negative..." begin

pts = -[0.15; 0.2; 0.35; 0.4]

p = kde!(pts, [0.05])

pcp = kde!(pts, [0.05], (addtheta,) )

pcm = kde!(pts, [0.05], (addtheta,), (difftheta,) )

# p.bt.centers
@test norm(p.bt.centers - pcp.bt.centers) < 1e-10
@test norm(p.bt.centers - pcm.bt.centers) < 1e-10

@test norm(p.bt.ranges - pcp.bt.ranges) < 1e-10
@test norm(p.bt.ranges - pcm.bt.ranges) < 1e-10

@test norm(p.bt.weights - pcp.bt.weights) < 1e-10
@test norm(p.bt.weights - pcm.bt.weights) < 1e-10

@test norm(p.bt.left_child - pcp.bt.left_child) < 1e-10
@test norm(p.bt.left_child - pcm.bt.left_child) < 1e-10

@test norm(p.bt.right_child - pcp.bt.right_child) < 1e-10
@test norm(p.bt.right_child - pcm.bt.right_child) < 1e-10

@test norm(p.bt.lowest_leaf - pcp.bt.lowest_leaf) < 1e-10
@test norm(p.bt.lowest_leaf - pcm.bt.lowest_leaf) < 1e-10

@test norm(p.bt.highest_leaf - pcp.bt.highest_leaf) < 1e-10
@test norm(p.bt.highest_leaf - pcm.bt.highest_leaf) < 1e-10

@test norm(p.bt.permutation - pcp.bt.permutation) < 1e-10
@test norm(p.bt.permutation - pcm.bt.permutation) < 1e-10

end

##

@testset "test circular manifold BallTree construction around 0rad..." begin

pts = [-0.2; 0.15; 0.15; 0.2]

p = kde!(pts, [0.05])

pcp = kde!(pts, [0.05], (addtheta,) )

pcm = kde!(pts, [0.05], (addtheta,), (difftheta,) )

# p.bt.centers
@test norm(p.bt.centers - pcp.bt.centers) < 1e-10
@test norm(p.bt.centers - pcm.bt.centers) < 1e-10

@test norm(p.bt.ranges - pcp.bt.ranges) < 1e-10
@test norm(p.bt.ranges - pcm.bt.ranges) < 1e-10

@test norm(p.bt.weights - pcp.bt.weights) < 1e-10
@test norm(p.bt.weights - pcm.bt.weights) < 1e-10

@test norm(p.bt.left_child - pcp.bt.left_child) < 1e-10
@test norm(p.bt.left_child - pcm.bt.left_child) < 1e-10

@test norm(p.bt.right_child - pcp.bt.right_child) < 1e-10
@test norm(p.bt.right_child - pcm.bt.right_child) < 1e-10

@test norm(p.bt.lowest_leaf - pcp.bt.lowest_leaf) < 1e-10
@test norm(p.bt.lowest_leaf - pcm.bt.lowest_leaf) < 1e-10

@test norm(p.bt.highest_leaf - pcp.bt.highest_leaf) < 1e-10
@test norm(p.bt.highest_leaf - pcm.bt.highest_leaf) < 1e-10

@test norm(p.bt.permutation - pcp.bt.permutation) < 1e-10
@test norm(p.bt.permutation - pcm.bt.permutation) < 1e-10

end

##

@testset "test circular manifold BallTree construction around π wrap point..." begin

# pts = [pi-0.2; pi-0.15; -pi+0.149; -pi+0.2]
pts = [-pi+0.15; -pi+0.2; +pi-0.2; pi-0.149]

# p = kde!(pts, [0.05])

# pcp = kde!(pts, [0.05], (addtheta,) )

pcm = kde!(pts, [0.05], (addtheta,), (difftheta,) )

@test abs(pcm.bt.centers[1] + pi) < 1e-3

@test abs(pcm.bt.centers[2] - pi + 0.175) < 1e-3
@test abs(pcm.bt.centers[3] + pi - 0.175) < 1e-3

@test abs(pcm.bt.centers[5] - pts[3]) < 1e-3
@test abs(pcm.bt.centers[6] - pts[4]) < 1e-3
@test abs(pcm.bt.centers[7] - pts[1]) < 1e-3
@test abs(pcm.bt.centers[8] - pts[2]) < 1e-3

end


text:
function orig_eval(exp::Expr)
    xs = 0.0:0.1:5
    ys = 0.0:0.1:5
    ls = []
    fun = @eval (xg, yg) -> $exp
    # fun = eval(exp)
    for x in xs
        for y in ys
            if Base.invokelatest(fun, x, y)
                push!(ls, (x, y))
            end
        end
    end
    return ls
end

function get_pose_from_fxn!(exp::Expr, ls::Vector{Any})
    xs = 0.0:0.1:5
    ys = 0.0:0.1:5
    fun = @eval (xg, yg) -> $exp
    for x in xs
        for y in ys
            if Base.invokelatest(fun, x, y)
                push!(ls, (x, y))
            end
        end
    end
end

function get_poses_from_fxn(exps::Vector{Expr})
    ls = []
    for exp in exps
        get_pose_from_fxn!(exp, ls)
    end
    return collect(Set(ls))
end

function and_get_pose_from_fxn(exps::Vector{Expr})
    xs = 0.0:0.1:5
    ys = 0.0:0.1:5
    ls = []
    for xg in xs
        for yg in ys
            sat = true
            for exp in exps
                @eval f(xg, yg)=$exp
                sat = sat && f(xg, yg)
            end
            if sat push!(ls, (xg, yg)) end
        end
    end
    return ls
end

function or_get_pose_from_fxn(exps::Vector{Expr})
    xs = 0.0:0.1:5
    ys = 0.0:0.1:5
    ls = []
    for exp in exps
        for xg in xs
            for yg in ys
                @eval f(xg, yg)=$exp
                if f(xg, yg)
                    push!(ls, (xg, yg))
                end
            end
        end
        println("length ls ",length(ls))
    end
    return ls
end
# xg=0.; yg=0.
get_poses_from_fxn([:(0.0<=xg<=0.05&&0.0<=yg<=1.5) , :(2xg+3yg>=0.7&&2xg+3yg<=1.7)] )

# orig_eval(:(0.0<=xg<=0.05&&0.0<=yg<=1.5))


text:
function solve_node(τ1::Array{R, N},
                    α::Vector{R},
                    β::Vector{R},
                    τ0::Array{R, N},
                    ∇τ0::Array{T},
                    tags::Array{UInt8, N},
                    κ2::Array{R, N},
                    x::S,
                    dx::Vector{<:AbstractFloat},
                    cs::Vector{S},
                    I1::S,
                    Iend::S) where {N, R <: AbstractFloat, T <: Array{R, N}, S <: CartesianIndex{N}}
    #loop over dimensions to set \alpha and \beta ...
    for (i, s) in enumerate(cs)
        fwdok = check1fwd(τ1, τ0, tags, x, s, Iend)
        bwdok = check1bwd(τ1, τ0, tags, x, s, I1)
        if fwdok && bwdok
            if @inbounds τ1[x+s]*τ0[x+s] < τ1[x-s]*τ0[x-s]
                bwkok = false
            else
                fwdok = false
            end
        end
        # Note Here: working out the maths gives an additional 1/\alpha for the \beta terms
        # to reduce total operations we redefine \beta = \beta' / \alpha where \beta' is the \beta in the paper
        if fwdok
            if check2fwd(τ1, τ0, tags, x, s, Iend)
                @inbounds α[i] = 1.5*τ0[x]/dx[i]-∇τ0[i][x]
                @inbounds β[i] = (2*τ1[x+s]-0.5*τ1[x+2*s])*τ0[x]/dx[i]
            else
                @inbounds α[i] = τ0[x]/dx[i]-∇τ0[i][x]
                @inbounds β[i] = τ0[x]*τ1[x+s]/dx[i]
            end
        elseif bwdok
            if check2bwd(τ1, τ0, tags, x, s, I1)
                @inbounds α[i] = 1.5*τ0[x]/dx[i]+∇τ0[i][x]
                @inbounds β[i] = (2*τ1[x-s]-0.5*τ1[x-2*s])*τ0[x]/dx[i]
            else
                @inbounds α[i] = τ0[x]/dx[i]+∇τ0[i][x]
                @inbounds β[i] = τ0[x]*τ1[x-s]/dx[i]
            end
        else
            @inbounds α[i] = zero(R)
            @inbounds β[i] = zero(R)
        end
    end
    #All of the \alpha and \beta should be set, so now we will proceed with solving the solve_quadratic
    #@inbounds solve_piecewise_quadratic(α..., β..., κ2[x])
    if N == 1
        @inbounds τ1y = solve_piecewise_quadratic(α[1], β[1], κ2[x])
    elseif N == 2
        @inbounds τ1y = solve_piecewise_quadratic(α[1], α[2], β[1], β[2], κ2[x])
    elseif N == 3
        @inbounds τ1y = solve_piecewise_quadratic(α[1], α[2], α[3], β[1], β[2], β[3], κ2[x])
    end
    τ1y
end

function fefmm_loop!(τ1::Array{R, N},
                     ocount::Integer,
                     α::Vector{R},
                     β::Vector{R},
                     τ0::Array{R, N},
                     ∇τ0::Array{T},
                     tags::Array{UInt8, N},
                     front::BinaryMinHeap{Node{R, N}},
                     κ2::Array{R, N},
                     dx::Vector{<:AbstractFloat},
                     cs::Vector{S},
                     xn::Vector{S},
                     I1::S,
                     Iend::S,
                     LI::LinearIndices{N}) where {N, R <: AbstractFloat, T <: Array{R, N}, S <: CartesianIndex{N}}
    while isempty(front) == false
        x = pop!(front).ind
        @inbounds if tags[x] < 0x3# as we are using a non-mutable minheap, we might have already set this node to known on another pass (i.e. this could be an old worthless version of the node). This check makes sure we don't accidentally override something we already have computed
            ocount += 1
            @inbounds tags[x] = 0x3
            set_neighbours!(xn, x, tags, cs, I1, Iend)
            for y in xn
                if y != x # to avoid allocating neighbours repeatedly we reuse xn but have y = x as the "neighbour not acceptable" flag value
                    @inbounds tags[y] = 0x2 #move neighbours to front
                    τ1y = solve_node(τ1, α, β, τ0, ∇τ0, tags, κ2, y, dx, cs, I1, Iend)
                    @inbounds if τ1y < τ1[y]
                            #only update and add new node to heap if the new solution is smaller
                        @inbounds τ1[y] = τ1y
                        @inbounds push!(front, Node(y, τ1[y]*τ0[y]))
                    end
                end
            end
        end
    end
end

"""
Inputs to fefmm:
κ - array of slowness
dx - array of grid spacings
x0 - linear index of source location

Outputs of fefmm:
τ1τ0 = travel time array from source

Tags:
1 = unknown
2 = front
3 = known
"""

function fefmm(κ2::Array{R, N},
               dx::Vector{<:AbstractFloat},
               xs::CartesianIndex{N}) where {R <: AbstractFloat, N}
    #initialization
    cs = cartstrides(κ2)
    inds = CartesianIndices(κ2)
    I1 = first(inds)
    Iend = last(inds)
    τ0 = ones(R, size(κ2))
    mul_analytic!(τ0, dx..., size(κ2)..., Tuple(xs)...)
    ∇τ0 = grad_analytic(τ0, dx..., size(κ2)..., Tuple(xs)...)
    τ1 = R(Inf).*ones(R, size(κ2))
    τ1[xs] = sqrt(κ2[xs])
    tags = ones(UInt8,size(κ2))
    tags[xs] = 0x2
    front = BinaryMinHeap{Node{R,N}}()
    push!(front, Node(inds[xs], τ1[xs]*τ0[xs]))
    α = Array{R}(undef,length(size(τ1)))
    β = similar(α)
    LI = LinearIndices(τ1)
    xn = Array{typeof(I1)}(undef, N*2)
    #main loop
    fefmm_loop!(τ1, 1, α, β, τ0, ∇τ0, tags, front, κ2, dx, cs, xn, I1, Iend, LI)
    τ1.*τ0
end


text:
using Documenter, ProjectedGPE

makedocs(
    modules = [ProjectedGPE],
    doctest=true
)

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/AshtonSBradley/ProjectedGPE.jl.git",
    julia  = "0.5.2",
    osname = "macos")

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
using Test, StaticArrays
using LogisticRegression: sigmoid

@test isapprox(sigmoid(0.0), 0.5);
@test isapprox(sigmoid(2.0), 0.88079708);
@test isapprox(sigmoid.([0.0, 2.0]), [0.5, 0.88079708])
@test isapprox(sigmoid.(@SVector [0.0, 2.0]), [0.5, 0.88079708])


text:
# Interface

abstract type AbstractAccumulator{T}; end

Base.getindex(s::T) where {T<:AbstractAccumulator} = error("getindex not defined for $T")
Base.setindex!(s::T, val) where {T<:AbstractAccumulator} = error("setindex! not defined for $T")

function update!(acc::T, rw, sd) where T<:AbstractAccumulator
    error("update! not defined for $T")
end

################################################################################

to_sign(b::Bool) = ifelse(b, 1, -1)

mutable struct ValueAccumulator{T} <: AbstractAccumulator{T}
    γ::T
    value::T
end

ValueAccumulator(γ::T) where {T} = ValueAccumulator(γ, zero(T))

Base.getindex(s::ValueAccumulator) = s.value
Base.setindex!(s::ValueAccumulator, val) = (s.value = val; s.value)

function update!(acc::ValueAccumulator{T}, rw, sd)::T where T
    setindex!(acc, (one(T)-acc.γ) * acc[] + acc.γ * rw * to_sign(sd))
end

################################################################################

function to_diff(s::T)::T where T
    x = 1-1/(1+s)
    x - (1 - x)
end

function to_ratio(s::T)::T where T
    x = (s+1)/2
    1 / (1-x) - 1
end

mutable struct InferenceAccumulator{T} <: AbstractAccumulator{T}
    prob::TaskStats{T}
    value::T
    side::Union{Bool, Nothing}
    reward_evidence::T
    failure_evidence::T
    ratio::T
    function InferenceAccumulator(prob::TaskStats{T}, value::T = zero(T), side::Union{Bool, Nothing} = nothing;
                                  reward_evidence::T = zero(T),
                                  failure_evidence::T = one(T) / (one(T) - prob.rwd)) where {T}
        new{T}(prob, value, side, reward_evidence, failure_evidence, to_ratio(value/prob.rwd))
    end
end

Base.getindex(s::InferenceAccumulator) = s.value
Base.setindex!(s::InferenceAccumulator, val) = (s.value = val; s.value)

function update!(acc::InferenceAccumulator{T}, rwd, sd)::T where T
    prob = acc.prob
    ratio = acc.ratio
    prior = if acc.side === nothing
        ratio
    else
        prevlow = acc.side ? one(T) / ratio : ratio
        lowprior = (prevlow + prob.dpl) / (one(T) - prob.dpl)
        (acc.side == sd) ? lowprior : one(T) / lowprior
    end
    evidence = ifelse(rwd, acc.reward_evidence, acc.failure_evidence)
    posterior = prior * evidence
    acc.side = sd
    acc.ratio = sd ? one(T) / posterior : posterior
    setindex!(acc, to_diff(acc.ratio)*prob.rwd)
end

###############################################################################

struct AccumulateFromData{T, I, W<:AbstractAccumulator{T}}
    acc::W
    itr::I
    AccumulateFromData(acc::AbstractAccumulator{T}, itr::I) where {T, I} =
        new{T, I, typeof(acc)}(acc, itr)
end

getiterator(s::AccumulateFromData) = s.itr
getaccumulator(s::AccumulateFromData) = s.acc

function Base.iterate(acc::AccumulateFromData, args...)
    next = iterate(getiterator(acc), args...)
    next === nothing && return nothing
    (val, status) = next
    rwd, sd = val
    return update!(getaccumulator(acc), rwd, sd), status
end

Base.IteratorSize(::Type{AccumulateFromData{T, I}}) where {T, I} = Base.IteratorSize(I)
Base.IteratorEltype(::Type{AccumulateFromData{T, I}}) where {T, I} = Base.HasEltype()

Base.length(acc::AccumulateFromData) = length(getiterator(acc))
Base.size(acc::AccumulateFromData) = size(getiterator(acc))
Base.eltype(::Type{<:AccumulateFromData{T, I}}) where {T, I} = T

Base.getindex(s::AccumulateFromData) = getindex(getaccumulator(s))
Base.setindex!(s::AccumulateFromData, val) = setindex!(getaccumulator(s), val)

accumulatefromdata(acc::AbstractAccumulator, rewards, sides) =
    collect(AccumulateFromData(acc, zip(rewards, sides)))

function valuediff(γ::T, rewards::AbstractVector{Bool}, sides::AbstractVector{Bool}) where T
    accumulatefromdata(ValueAccumulator(γ), rewards, sides)
end

function valuediff(prob::TaskStats{T}, rewards::AbstractVector{Bool}, sides::AbstractVector{Bool}) where T
    ratios = probabilityratio(prob, rewards, sides)
    to_diff.(ratios) .* prob.rwd
end

function probabilityratio(prob::TaskStats{T}, rewards::AbstractVector{Bool}, sides::AbstractVector{Bool}) where T
    accumulatefromdata(InferenceAccumulator(prob), rewards, sides)
end


text:
using DataFrames
using Statistics # for using mean
using Pipe: @pipe # for @pipe macro
using RDatasets # for iris dataset

iris = dataset("datasets", "iris");

import FloatingTableView
FloatingTableView.browse(iris)

# R"""
# iris %>%
#   group_by(Species) %>%
#   summarise(across(starts_with("Sepal"), mean))
# """

filter(startswith("Sepal"), names(iris))

@pipe iris |>
      groupby(_, :Species) |>
      combine(_, filter(startswith("Sepal"), names(_)) .=> mean)

# R"""
# iris %>%
# as_tibble() %>%
# mutate(across(where(is.factor), as.character))
# """

typeof.(eachcol(iris)) .<: CategoricalArray
n = names(iris)[typeof.(eachcol(iris)).<:CategoricalArray]

catarr = categorical(["abc", "def", "abc"])

Vector{String}(catarr)

# you can use the `string` function to convert any type to string
@pipe iris |> transform(_, n .=> Vector{String} .=> n)

# Will come back to this case at the end to show you how to simplify this!

# A purrr-style formula
# iris %>%
#   group_by(Species) %>%
#   summarise(across(starts_with("Sepal"), ~mean(.x, na.rm = TRUE)))
@pipe iris |>
      groupby(_, :Species) |>
      combine(_, filter(startswith("Sepal"), names(_)) .=> x -> mean(x |> skipmissing))

# A named list of functions
# iris %>%
#   group_by(Species) %>%
#   summarise(across(starts_with("Sepal"), list(mean = mean, sd = sd)))

# this doesn't work; because it is not a cartesian product
@pipe iris |>
      groupby(_, :Species) |>
      combine(_, filter(startswith("Sepal"), names(iris)) .=> (mean, std))

@pipe iris |>
      groupby(_, :Species) |>
      combine(
          _,
          [name => fn for
          name in filter(startswith("Sepal"), names(_)), fn in (mean, std)]...,
      )

function dplyr_across(names, fns)
    (name => fn for name in names, fn in fns)
end

@pipe iris |>
      groupby(_, :Species) |>
      combine(_, dplyr_across(filter(startswith("Sepal"), names(_)), (mean, std))...)

struct DplyrAcross
    innames::Any
    fn::Any
    outnames::Any
end

DplyrAcross(innames, fn) = DplyrAcross(innames, fn, nothing)

import DataFrames.combine

function DataFrames.combine(df, da::DplyrAcross)
    input_cols = filter(da.innames, names(df))

    if isnothing(da.outnames)
        pairs = [col => fn for col in input_cols, fn in da.fn]
        return combine(df, pairs...)
    else
        error("not implemented yet")
    end
end

@pipe iris |>
      groupby(_, :Species) |>
      combine(_, DplyrAcross(startswith("Sepal"), (mean, std)))

# Use the .names argument to control the output names
# iris %>%
#   group_by(Species) %>%
#   summarise(across(starts_with("Sepal"), mean, .names = "mean_{col}"))

DplyrAcross(innames, fn; names = nothing) = DplyrAcross(innames, fn, names)

function DataFrames.combine(df, da::DplyrAcross)
    input_cols = filter(da.innames, names(df))

    if isnothing(da.outnames)
        pairs = [col => fn for col in input_cols, fn in da.fn]
        return combine(df, pairs...)
    else
        outcols = [
            replace(replace(da.outnames, "{col}" => col), "{fn}" => string(fn))
            for col in input_cols, fn in keys(da.fn)
        ]

        in_fn_out_iterator = zip(Iterators.product(input_cols, da.fn), outcols)

        pairs = [col => fn => outcol for ((col, fn), outcol) in in_fn_out_iterator]

        return combine(df, pairs...)
    end
end

@pipe iris |>
      groupby(_, :Species) |>
      combine(_, DplyrAcross(startswith("Sepal"), [mean]; names = "mean_{col}"))

# iris %>%
#   group_by(Species) %>%
#   summarise(across(starts_with("Sepal"), list(mean = mean, sd = sd), .names = "{col}.{fn}"))

@pipe iris |>
      groupby(_, :Species) |>
      combine(
          _,
          DplyrAcross(startswith("Sepal"), (mean = mean, std = std); names = "{col}.{fn}"),
      )

@pipe iris |>
      groupby(_, :Species) |>
      combine(
          _,
          DplyrAcross(startswith("Sepal"), (mean = mean, std = std); names = "{col}_{fn}"),
      )

# iris %>%
#   group_by(Species) %>%
#   summarise(across(starts_with("Sepal"), list(mean, sd), .names = "{col}.fn{fn}"))

@pipe iris |>
      groupby(_, :Species) |>
      combine(_, DplyrAcross(startswith("Sepal"), (mean, std); names = "{col}_fn{fn}"))

using Dply

@pipe iris |>
      groupby(_, :Species) |>
      combine(_, Dply.Across(startswith("Sepal"), (mean, std); names = "{col}_fn{fn}"))

@pipe iris |>
      groupby(_, :Species) |>
      combine(
          _,
          Dply.across(
              startswith("Sepal"),
              (mean = mean, std = std);
              names = "{col}_fn_{fn}",
          ),
      )

#######################################
## Do you recall this case we skipped?
#######################################

# R"""
# iris %>%
# as_tibble() %>%
# mutate(across(where(is.factor), as.character))
# """

typeof.(eachcol(iris)) .<: CategoricalArray
n = names(iris)[typeof.(eachcol(iris)).<:CategoricalArray]

@pipe iris |> transform(_, n .=> Vector{String} .=> n)

@pipe iris |> transform(_, n .=> (col -> string.(col)) .=> n)

using Dply

@pipe iris |> transform(
    _,
    Dply.across(Dply.Where(col -> typeof(col) <: CategoricalArray), Vector{String}),
)

# define a helper function
iscatarray(col) = typeof(col) <: CategoricalArray

@pipe iris |> transform(_, Dply.across(Dply.Where(iscatarray), Vector{String}))

using Dply: across, where

@pipe iris |> transform(_, across(where(iscatarray), Vector{String}))

@pipe iris |> transform(_, across(where(iscatarray), col -> string.(col)))


text:
export PDBGeom

abstract type PDBGeom end

"""
    function create(geom::PDBGeom; resolution=0.1, rand_=0.05, type=1)

This create a PDB geometry.
"""
function PDMesh.create(geom::PDBGeom; resolution=0.1, rand_=0.05, type=1)
    error("Not implemented for type **$(typeof(geom))** yet.")
end

include("./bar.jl")

max_stars_repo_name: ravinderbhattoo/PDBenchmark

max_issues_repo_name: ravinderbhattoo/PDBenchmark

max_forks_repo_name: ravinderbhattoo/PDBenchmark


text:
# ------------------------------------------------------------------------------
#                           BLAS SERVICE FUNCTION
# ------------------------------------------------------------------------------
#=
enum CBLAS_LAYOUT {
   CblasRowMajor=101,  /* row-major arrays */
   CblasColMajor=102};   /* column-major arrays */
enum CBLAS_TRANSPOSE {
   CblasNoTrans=111,     /* trans='N' */
   CblasTrans=112,       /* trans='T' */
   CblasConjTrans=113};  /* trans='C' */
enum CBLAS_UPLO {
   CblasUpper=121,        /* uplo ='U' */
   CblasLower=122};       /* uplo ='L' */
enum CBLAS_DIAG {
   CblasNonUnit=131,      /* diag ='N' */
   CblasUnit=132};        /* diag ='U' */
enum CBLAS_SIDE {
   CblasLeft=141,         /* side ='L' */
   CblasRight=142};       /* side ='R' */

   #ifndef __MKL_CBLAS_H__
   #define __MKL_CBLAS_H__
   #include <stddef.h>

   #include "mkl_types.h"

   #ifdef __cplusplus
   extern "C" {            /* Assume C declarations for C++ */
   #endif /* __cplusplus */

   /*
    * Enumerated and derived types
    */
   #define CBLAS_INDEX size_t  /* this may vary between platforms */

   typedef enum {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER;
   typedef enum {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113} CBLAS_TRANSPOSE;
   typedef enum {CblasUpper=121, CblasLower=122} CBLAS_UPLO;
   typedef enum {CblasNonUnit=131, CblasUnit=132} CBLAS_DIAG;
   typedef enum {CblasLeft=141, CblasRight=142} CBLAS_SIDE;
=#

immutable BLAS

# credit : Simon Kornblith (VML.jl)
immutable VMLAccuracy
    mode::UInt
end

const VML_LA = VMLAccuracy(0x00000001)
const VML_HA = VMLAccuracy(0x00000002)
const VML_EP = VMLAccuracy(0x00000003)

Base.show(io::IO, m::VMLAccuracy) = print(io, m == VML_LA ? "VML_LA" :
                                              m == VML_HA ? "VML_HA" : "VML_EP")

vml_get_mode() = ccall((:_vmlGetMode, libvml), Cuint, ())
vml_set_mode(mode::Integer) = (ccall((:_vmlSetMode, libvml), Cuint, (UInt,), mode); nothing)

vml_set_accuracy(m::VMLAccuracy) = vml_set_mode((vml_get_mode() & ~0x03) | m.mode)
vml_get_accuracy() = VMLAccuracy(vml_get_mode() & 0x3)

vml_set_mode((vml_get_mode() & ~0x0000FF00))

function vml_check_error()
    vml_error = ccall((:_vmlClearErrStatus, libvml), Cint, ())
    if vml_error != 0
        if vml_error == 1
            error(DomainError())
        elseif vml_error == 2 || vml_error == 3 || vml_error == 4
            # Singularity, overflow, or underflow
            # I don't think Base throws on these
        elseif vml_error == 1000
            warn("VML does not support $(vml_get_accuracy); lower accuracy used instead")
        else
            error("an unexpected error occurred in VML ($vml_error)")
        end
    end
end