text:
function xxz_hamiltonian(lx::Int64, delta::Float64=1.0;
                         boundary::Symbol=:OBC, mode::Symbol=:FULL)
    @assert lx > 1
    if mode == :FULL
        return xxz_explicit(lx, delta, boundary)
    elseif mode == :U1
        return xxz_u1sym(lx, delta, boundary=boundary)
    else
        error("unrecognized mode :", mode)
    end
end

function xxz_explicit(lx::Int64, delta::Float64, boundary::Symbol)

    sz, sp, sm = sparse.(spinoperators(1/2))
    heis_term = 0.5 * (kron(sp, sm) + kron(sm, sp)) +
        delta * kron(sz,sz)

    Hmat = heis_term
    if lx > 2
        # recursive generation of open chain
        for i=3:lx
            Hmat = kron(Hmat, I(2)) + kron(I(2^(i-2)), heis_term)
        end
    end

    if boundary == :OBC
        return Hmat
    elseif boundary == :PBC
        Hmat = Hmat +
            0.5 * kron(kron(sp, I(2^(lx-2))), sm) +
            0.5 * kron(kron(sm, I(2^(lx-2))), sp) +
            delta * kron(kron(sz, I(2^(lx-2))), sz)
        return Hmat
    else
        error("unrecognized boundary condition :", boundary)
    end
end

"""
    xxz_u1sym(Lx, delta; boundary, zsector)
"""

# # Generates the XXZ Hamiltonian of size `lx` in the given `zsector`
# # block. This is possible because ``\left[ S^z_{\text{tot}},
# # H\right]=0`` and its eigenvalues can be used to label separate
# # magnetization blocks in the Hamiltonian matrix and each block can be
# # diagonalized separately to significantly reduce time complexity. Each
# # `zsector` block consists of Ising vectors with exactly `M = (zsector
# # + lx)/2` spins up, so the size of the block is ``\binom{lx}{M}``. Using
# # Sterling formula the size of the largest block, ``M=lx/2``, is ``√{π
# # lx/2}`` times smaller than the full Hilbert space, which For typical
# # sizes of ED calculation ,``lx∼ 30``, is roughly about ``∼\!7`` times.

# """
function xxz_u1sym(lx::Int, delta::Float64=1.0;
                   boundary::Symbol=:OBC, zsector::Int=(lx % 2))

    @assert lx > 1
    @assert (lx-zsector) % 2 == 0

    if boundary == :OBC
        bond_range = 0:lx-2
    elseif boundary == :PBC
        bond_range = 0:lx-1
    else
        error("unrecognized boundary condition :", boundary)
    end

    M = div(lx + zsector, 2)
    block_size = binomial(lx, M)

    I = Int32[]
    J = Int32[]
    V = Float64[]

    szblock_states = Vector{Int}(undef, block_size)
    index = 0
    for state=0:2^lx-1
        if count_ones(state) == M
            index += 1
            szblock_states[index] = state
        end
    end

    # For all states (of magnetization M):
    for state_index=1:block_size
        state = szblock_states[state_index]
        # For all bond terms of the Hamiltonian
        for site=bond_range
            site_nn = (site + 1) % lx
            # if spin are parallel:
            if xor((state >> site) & 1, (state >> site_nn) & 1) == 0
                push!(I, state_index)
                push!(J, state_index)
                push!(V, delta * 0.25)
                # if spin are not parallel:
            else
                push!(I, state_index)
                push!(J, state_index)
                push!(V, - delta * 0.25)

                flipped_state = xor(state, 1 << site | 1 << site_nn)
                flipped_state_index = searchsortedfirst(szblock_states, flipped_state)
                push!(I, state_index)
                push!(J, flipped_state_index)
                push!(V, 0.5)
            end
        end
    end

    return sparse(I, J, V, block_size, block_size, +), (szblock_states .+ 1)
 end

"""
    j1j2_explicit(Lx, j1, j2)

explicit construction of j1j2 sparse full Hamiltonian.

"""
function j1j2_explicit(Lx::Int64,
                       j1::Float64,
                       j2::Float64;
                       boundary::Symbol=:open)
    @assert Lx > 2
    Sz, Sp, Sm = sparse.(spinoperators(1/2))
    I2 = I(2)

    heis_term1 = j1 * 0.5 * (kron(Sp, Sm) + kron(Sm, Sp)) + j1 * kron(Sz,Sz)
    heis_term2 = j2 * 0.5 * (kron(kron(Sp, I2), Sm) + kron(kron(Sm, I2), Sp)) +
        j2 * kron(kron(Sz, I2), Sz)

    Hmat = heis_term1
    for i=3:Lx
        Hmat = kron(Hmat, I2) + kron(I(2^(i-2)), heis_term1) +
            kron(I(2^(i-3)), heis_term2)
    end

    if boundary == :open
        return Hmat
    elseif boundary == :periodic
        Hmat = Hmat +
            j1 * 0.5 * kron(kron(Sp,I(2^(Lx-2))), Sm) +
            j1 * 0.5 * kron(kron(Sm,I(2^(Lx-2))), Sp) +
            j1 * kron(kron(Sz,I(2^(Lx-2))), Sz) +
            j2 * 0.5 * kron(I2, kron(kron(Sp,I(2^(Lx-3))), Sm)) +
            j2 * 0.5 * kron(I2, kron(kron(Sm,I(2^(Lx-3))), Sp)) +
            j2 * kron(I2, kron(kron(Sz,I(2^(Lx-3))), Sz)) +
            j2 * 0.5 * kron(kron(kron(Sp,I(2^(Lx-3))), Sm), I2) +
            j2 * 0.5 * kron(kron(kron(Sm,I(2^(Lx-3))), Sp), I2) +
            j2 * kron(kron(kron(Sz,I(2^(Lx-3))), Sz), I2)
        return Hmat
    else
        error("unrecognized boundary conditions :", boundary)
    end
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
using Proquint
using Test

@testset "Proquint.jl" begin
    for i in 1:100
        r = rand(UInt32)
        q = uint2quint(r)
        @test quint2uint(q) == r
    end
    for i in 1:100
        r = rand(UInt32)
        q = uint2quint(r, short=true)
        @test quint2uint(q) == r
    end
end

ext: jl


text:
"""
    tolatex(iter; order = true)

Generate a TikzPicture object visualizing the order on a quadratic neighborhood
"""
function tolatex(iter; order = true)
    preprocess = preprocessor()
    Φ = preprocess(iter)
    sb = IOBuffer()
    el_dim = length(first(iter))
    max_val = 0
    for i in 1:el_dim
        max_val = max(max_val, maximum(map(x->x[i], iter)))
    end
    neigh_header(sb) # Add header containing appropriate tikzstyles
    quad_neigh(sb, max_val) # Draw parameter space
    order && neigh_order(sb, Φ) # Add elements and arrows
    latex_string = String(take!(sb))
    close(sb)
    return TikzPicture(latex_string)
end

as_tex(tp, fname) = save(TEX(fname, limit_to=:picture), tp)
as_svg(tp, fname) = save(SVG(fname), tp)

function neigh_header(sb)
    println(sb, "\\tikzstyle{arrow} = [thick, ->, >=stealth]")
    println(sb, "\\tikzstyle{op} = [circle, draw=black, fill=black, minimum width = .05cm, inner sep = 1];")
end

function quad_neigh(sb, n)
    grid =
    """
    \\foreach \\x in {1,...,$n}{
        \\foreach \\y in {1,...,$n}{
            \\node[op] (oper \\x_\\y) at (\\x,  -\\y) {};
        }
    }
    \\foreach \\x in {1,...,$n}{
	\\draw (0, -\\x) node[text=gray] {\\x};
	\\draw (\\x, 0) node[text=gray] {\\x};
}
    """
    println(sb, grid)
end

function neigh_order(sb, Φ)
    edges =
    """
    \\foreach \\s / \\d in {$(Φ)}{
        \\draw[arrow] (oper \\s) -- (oper \\d);
    }
    """
    println(sb, edges)
end

"""
    preprocessor()

Returns a function which -- upon executing on an order -- returns a tikz-compatible string.
"""
function preprocessor()
    zipper = iter -> begin v = collect(iter); Iterators.zip(v[1:end-1], v[2:end]) end
    parser = param -> join(map(x->join(x, "_"), param), "/")
    return Φ -> join(map(Φ -> parser(Φ), zipper(Φ)), ", ")
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz


text:
using Test
using MiniTB

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

# @testset begin
 p1 = Point2D( 1.011,2.011)
 p2 = Point2D( 1.011, 2.011 )
 p3 = Point2D( 1.011, 3.011 )
 p4 = Point2D( -1.011, -3.011 )
 p5 = Point2D( 10.011, 13.011 )

 bc = BudgetConstraint( [p1,p2,p3,p4,p5])

sort!( bc )

@test size( bc )[1] == 5

ps = PointsSet([p1,p2,p3,p4,p5] )
push!( ps, p1 )
push!( ps, p2 )
push!( ps, p3 )
push!( ps, p4 )
push!( ps, p5 )
println( "ps = $ps ")
push!( ps, Point2D( 1.0011, 1.999999 ))

bc2 = censor( ps )

println( "bc2 = $bc2 ")
@test size( bc2 )[1] == 5
pers = DEFAULT_PERSON
pars = DEFAULT_PARAMS
res = calculate( pers, pars )
println( "res=$res" )

 # @test size( ps ) == 4

function makebc( pers :: Person, params :: Parameters ) :: BudgetConstraint

    function getnet( gross :: Float64 ) :: Float64
        persedit = modifiedcopy( pers, wage=gross )
        # println( "getnet; made person $persedit")
        rc = calculate( persedit, params )
        return rc[:netincome]
    end

    bc = makebc( getnet )

    return bc
end

bc = makebc( DEFAULT_PERSON, DEFAULT_PARAMS )
println(pointstoarray( bc))
bc = makebc( DEFAULT_PERSON, ZERO_PARAMS )
println( pointstoarray( bc))
pars = DEFAULT_PARAMS
pars.it_allow = 999.0
person = DEFAULT_PERSON
#person.wage -= 1
rc = calculate( DEFAULT_PERSON, pars )
println(rc)
bc = makebc( person, pars )
println( pointstoarray( bc ))


text:
using Bokeh

plotfile("functions.html")

# Bokeh.jl tries hard to follow TextPlots.jl's API particularly for functions
# (in future there might be "use textplots" switch in Bokeh.jl)
# the following examples are all equivilent
plot(sin)
plot(sin, -10:10)
plot(sin, -10, 10)

plot([sin, cos, tan, sinh, cosh, tanh], -1:1)
showplot()


text:
using CSV, DataFrames

const names_url = "https://files.osf.io/v1/resources/jum2f/providers/osfstorage/5d4c2bc722a88600165004f6?action=download&direct&version=2"
const bottom_up_url = "https://files.osf.io/v1/resources/jum2f/providers/osfstorage/5d4ebd7d0f488d0016909bfd?action=download&direct&version=1"
const top_down_url = "https://files.osf.io/v1/resources/jum2f/providers/osfstorage/5d4ebd8422a886001750ee0b?action=download&direct&version=1"
const manual_url = "https://files.osf.io/v1/resources/jum2f/providers/osfstorage/5d4ebd7e22a8860019504e43?action=download&direct&version=1"

load_tsv(file) = DataFrame(CSV.File(file, header = true, delim = '\t'))
download_tsv(url) = load_tsv(download(url))

function map_categories(df::DataFrame, object_names::Vector{String})
    mapped_categories = Dict()

    for (i, row) in enumerate(eachrow(df))
        positive_match = findfirst(x -> x == 1, row)
        category = isnothing(positive_match) ? "" : positive_match

        push!(mapped_categories, object_names[i] => string(category))
    end

    return mapped_categories
end

function get_files(dir::String)
    return collect(Base.Iterators.flatten([files for (_, _, files) in walkdir(dir)]))
end

function main(distractor_dir::String, out_file::String)
    object_names = download_tsv(names_url).uniqueID
    df_bottom_up = download_tsv(bottom_up_url)
    df_top_down = download_tsv(top_down_url)
    df_manual = download_tsv(manual_url)

    mapped_categories = Dict(
        "bottom up" => map_categories(df_bottom_up, object_names),
        "top_down" => map_categories(df_top_down, object_names),
        "manual" => map_categories(df_manual, object_names)
    )

    distractor_categories = DataFrame(file_names = [], bottom_up = [], top_down = [], manual = [])

    for distractor in get_files(distractor_dir)
        distractor == "." && continue
        distractor == ".DS_Store" && continue

        distractor_name = distractor[1:(end - 8)] # _[[:digit:]][[:digit:]].\.jpg

        temp_data = [distractor]
        for (map_type, mapped_dict) in mapped_categories
            distractor_category = mapped_dict[distractor_name]
            push!(temp_data, distractor_category)
        end

        push!(distractor_categories, temp_data)
    end

    CSV.write(out_file, distractor_categories)

    return distractor_categories
end

main("/Users/jakeireland/projects/GorillaSCExperiment/data/Original/", "/Users/jakeireland/projects/GorillaSCExperiment/data/distractor_categories.csv")

max_stars_repo_name: jakewilliami/GorillaSCExperiment

max_issues_repo_name: jakewilliami/GorillaSCExperiment

max_forks_repo_name: jakewilliami/GorillaSCExperiment


text:
"""
    univariate(k::Int64, l::Int64, low::N, high::N)::Vector{N} where {N<:AbstractFloat}

Compute the Bernstein coefficients of a univariate monomial.

### Input

- `k`    -- degree of the given monomial
- `l`    -- degree of the Bernstein polynomial
- `low`  -- the lower bound of the interval where the Bernstein coefficients are computed
- `high` -- the upper bound of the interval the Bernstein coefficients are computed

### Output

A vector with floating point entries containing the Bernstein coefficients.
"""
function univariate(k::Int64, l::Int64, low::N, high::N)::Vector{N} where {N<:AbstractFloat}
    m = l - k
    b = zeros(N, l+1)
    @inbounds @fastmath for i in 0:l
        if k < l
            for j in max(0, i-m):min(k, i)
                aux =  binomial(m, i-j) * binomial(k, j) / binomial(k+m, i)
                b[i+1] += aux * low^(k-j) * high^j
            end
        else
            b[i+1] = low^(k-i) * high^i
        end
    end
    return b
end

"""
    univariate(k::Int64, l::Int64, low::Rational, high::Rational)::Vector{Rational}

Compute *exactly* the Bernstein coefficients of a univariate monomial.

### Input

- `k`    -- degree of the given monomial
- `l`    -- degree of the Bernstein polynomial
- `low`  -- the lower bound of the interval where the Bernstein coefficients are computed
- `high` -- the upper bound of the interval the Bernstein coefficients are computed

### Output

A vector with rational entries containing the Bernstein coefficients.
"""
function univariate(k::Int64, l::Int64, low::Rational, high::Rational)::Vector{Rational}
    m = l - k
    b = zeros(Rational, l+1)
    @inbounds for i in 0:l
        if k < l
            for j in max(0, i-m):min(k, i)
                aux =  binomial(m, i-j) * binomial(k, j) // binomial(k+m, i)
                b[i+1] += aux * low^(k-j) * high^j
            end
        else
            b[i+1] = low^(k-i) * high^i
        end
    end
    return b
end

"""
    multivariate(k::Vector{Int64}, l::Vector{Int64},
                          low::Vector{N}, high::Vector{N})::ImplicitForm{N} where {N<:Number}

Compute the Bernstein coefficients of a multivariate monomial.

### Input

- `k`    -- vector of degrees for each monomial
- `l`    -- vector of degrees of the Bernstein polynomial for each monomial
- `low`  -- the lower bounds of the interval where the Bernstein coefficients are computed
- `high` -- the upper bounds of the interval the Bernstein coefficients are computed

### Output

A Bernstein implicit form holding the Bernstein coefficients.

### Examples

```julia
julia> m = multivariate([3,2],[3,2],[1.0,2],[2,4.0]);
julia> m.array
2-element Array{Array{Float64,1},1}:
 [1.0, 2.0, 4.0, 8.0]
 [4.0, 8.0, 16.0]
```
"""
function multivariate(k::Vector{Int64}, l::Vector{Int64},
                      low::Vector{N}, high::Vector{N})::ImplicitForm{N} where {N<:Number}
    n = length(low)
    ncoeffs = 1
    B = Vector{Vector{N}}(n)
    @inbounds for i in 1:n
        @views B[i] = univariate(k[i], l[i], low[i], high[i])
        ncoeffs = ncoeffs * length(B[i])
    end
    return ImplicitForm(B, n, ncoeffs)
end


text:
## Using basic pyqtgraph installed as follows:
## download source from http://www.pyqtgraph.org/
## installed globally with sudo python setup.py install

## Mostly From ScatterPlot.py example

using PySide
reload(Pkg.dir("PySide", "src", "pyqtgraph.jl"))
using PyQtGraph

w = Widget()
lyt = VBoxLayout(w)
setLayout(w, lyt)

## Graphics Layout Widget make grid of plot devices
win = GraphicsLayoutWidget(w)
addWidget(lyt, win)
set_size(w, 800, 600)

raise(w)

## Here we create plot objects and lay them out first:
w0 = addPlot(win)
nextRow(win)

w1 = addPlot(win)

## add an API for this?
w2 = addViewBox(win)
w2[:setAspectLocked](true)

nextRow(win)

w3 = addPlot(win)
w4 = addPlot(win)

raise(win)

## This follows the `pyqtgraph` example that this example is taken from:

## There are a few different ways we can draw scatter plots; each is optimized for different types of data:
## 0) The easiest way, using addPoints addtion to GraphicsPlot objects
x = 1:10; y = rand(10)
w0.addPoints(x, y)
w0.addLine(x=0, y=1)

## 1) All spots identical and transform-invariant (top-left plot).
## In this case we can get a huge performance boost by pre-rendering the spot
## image and just drawing that image repeatedly.

n = 300
x = randn(n); y = randn(n)
s1 = ScatterPlotItem(x, y, size=10, pen=nothing, brush=(255, 255,255,120))
w1.addItem(s1)

## Make all plots clickable
lastClicked = nothing
function clicked(plot, points)
    global lastClicked
    println(typeof(lastClicked))
    if !isa(lastClicked, Nothing)
        map(u -> u[:resetPen](), lastClicked)
    end
    map(u -> u[:setPen]("b", width=2), points)
    lastClicked = points
end

function get_lastclicked()
    x = map(i -> qinvoke(i, [:pos, :x]), lastClicked) |> float
    y = map(i -> qinvoke(i, [:pos, :y]), lastClicked) |> float
    [x y]
end

qconnect(s1, :sigClicked, clicked)

## 2) Spots are transform-invariant, but not identical (top-right plot).
## In this case, drawing is as fast as 1), but there is more startup overhead
## and memory usage since each spot generates its own pre-rendered image.

scale = 1e-5
x = rand(n)*scale; y = rand(n)*scale
s2 = ScatterPlotItem(x, y, size=[5 + (1:n)/10], pen="w", pxMode=true)
w2.addItem(s2)
qconnect(s2, :sigClicked, clicked)
#s2.sigClicked.connect(clicked)

## 3) Spots are not transform-invariant, not identical (bottom-left).
## This is the slowest case, since all spots must be completely re-drawn
## every time because their apparent transformation may have changed.

## Set pxMode=False to allow spots to transform with the win

x = [1e-6*i for  i in 1:10, j in 1:10] |> u -> reshape(u, 100)
y = [1e-6*j for  i in 1:10, j  in 1:10] |> u -> reshape(u, 100)
brush = [pyqtgraph.intColor(i*10+j, 100) for i in 1:10, j in 1:10] |> u -> reshape(u, 100)

s3 = ScatterPlotItem(x, y, pxMode=false, size=1e-6, pen={:color => "w", :width=>2}, brush=brush)

w3.addItem(s3)
qconnect(s3, :sigClicked,clicked)

## Test performance of large scatterplots
## This gets slow. Could not call as in exmaple
     ## addPoints(x=pos[0], y=pos[1])
n = 10^3                         # 1e4 is slow to respond
s4 = ScatterPlotItem(randn(n), randn(n), data=string(1:n), size=10, pen=nothing, brush=(255, 255, 255, 20))
w4.addItem(s4)
qconnect(s4, :sigClicked, clicked)


text:
@testset "generate" begin
    root = "assets"

    for theme in ("grid", )
        templates = DemoCards.cardtheme(Val(Symbol(theme)))[1]

        sec = DemoSection(joinpath(root, "section", "default"))
        @test_reference joinpath("references", "generate_section_$theme.txt") generate(sec, templates) by=ignore_CR

        page = DemoPage(joinpath(root, "page", "default"))
        @test_reference joinpath("references", "generate_page_$theme.txt") generate(page, templates) by=ignore_CR
    end

    abs_root = joinpath(pwd(), root, "page", "default")
    mktempdir() do root
        tmp_root = joinpath(root, "page", "default")
        mkpath(tmp_root)
        cp(abs_root, tmp_root, force=true)
        templates, theme = cardtheme(root=root)
        path, post_process_cb = @suppress_err makedemos(tmp_root, templates, root=root)
        @test @suppress_err path == joinpath("default", "index.md")
        @test theme == joinpath("democards", "gridtheme.css")
    end

    @testset "hidden" begin
        page_dir = @suppress_err preview_demos(joinpath(root, "page", "hidden"); theme="grid", require_html=false)
        md_index = joinpath(page_dir, "index.md")
        # test only one card is shown in index.md
        @test_reference joinpath("references", "hidden_index.txt") read(md_index, String)

        # check that assets are still normally generated even if they are hidden from index.md
        @test readdir(page_dir) == ["covers", "index.md", "sec"]
        @test readdir(joinpath(page_dir, "covers")) == ["democards_logo.svg"]
        @test readdir(joinpath(page_dir, "sec")) == ["hidden1.ipynb", "hidden1.jl", "hidden1.md", "hidden2.md", "normal.md"]
    end

    @testset "throw_error" begin
        @suppress_err @test_throws LoadError preview_demos(joinpath(root, "broken_demo.jl"); throw_error=true)
    end

    @testset "badges" begin
        @testset "author" begin
            for (dir, file) in [
                ("julia", joinpath(root, "card", "julia", "author_with_md_url.jl")),
                ("markdown", joinpath(root, "card", "markdown", "author_with_md_url.md"))
            ]
                page_dir = @suppress_err preview_demos(file; theme="grid", require_html=false)
                md_file = joinpath(page_dir, dir, "author_with_md_url.md")
                @test isfile(md_file)
                line = read(md_file, String)
                @test occursin("[![Author](https://img.shields.io/badge/Author-Johnny%20Chen-blue)](https://github.com/johnnychen94)", line)
                @test occursin("![Author](https://img.shields.io/badge/Author-Jane%20Doe-blue)", line)
            end
        end
    end
end


text:
# tests for kernel matrix methods

include("simple.jl")


text:
@testset "Constructors" begin
    @testset "Inner constructors" begin
        @test (@inferred Tensor{Tuple{2,2}, Int, 2, 4}((1,2,3,4)))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
        @test (@inferred Tensor{Tuple{2,2}, Float64, 2, 4}((1,2,3,4)))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,3.0,4.0)
        @test (@inferred Tensor{Tuple{2,2}, Int, 2, 4}((1.0,2,3,4)))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
        @test (@inferred Tensor{Tuple{2,2}, Float64, 2, 4}((1.0,2,3,4)))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,3.0,4.0)
        @test (@inferred Tensor{Tuple{2,Symmetry{Tuple{2,2}}}, Int, 3, 6}((1,2,3,4,5,6)))::Tensor{Tuple{2,Symmetry{Tuple{2,2}}}, Int, 3, 6} |> Tuple == (1,2,3,4,5,6)
        @test (@inferred Tensor{Tuple{2,Symmetry{Tuple{2,2}}}, Float64, 3, 6}((1,2,3,4,5,6)))::Tensor{Tuple{2,Symmetry{Tuple{2,2}}}, Float64, 3, 6} |> Tuple == (1.0,2.0,3.0,4.0,5.0,6.0)

        # bad input
        @test_throws Exception Tensor{Tuple{2,2}, Int, 2, 4}((1,2,3))
        @test_throws Exception Tensor{Tuple{2,2}, Int, 2, 4}(())
        @test_throws Exception Tensor{Tuple{2}, Vector{Int}, 1, 2}(([1,2],[3,4]))

        # bad parameters
        @test_throws Exception Tensor{Tuple{Int,2}, Int, 2, 4}((1,2,3,4)) # bad size
        @test_throws Exception Tensor{Tuple{2,2}, Int, 2, 5}((1,2,3,4,5)) # bad ncomponents
        @test_throws Exception Tensor{Tuple{2,2}, Int, 3, 4}((1,2,3,4))   # bad ndims
        @test_throws Exception Tensor{Tuple{2,Symmetry{Int,2}}, Int, 3, 6}((1,2,3,4,5,6)) # bad Symmetry size
        @test_throws Exception Tensor{Tuple{2,Symmetry{2,2}}, Int, 3, 5}((1,2,3,4,5))     # bad ncomponents
        @test_throws Exception Tensor{Tuple{2,Symmetry{2,2}}, Int, 2, 6}((1,2,3,4,5,6))   # bad ndims
    end
    @testset "Outer constructors" begin
        data = (1,2,3,4)
        @testset "From Tuple" begin
            # same type
            @test (@inferred Tensor{Tuple{2,2}, Int, 2}(data))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
            @test (@inferred Tensor{Tuple{2,2}, Int}(data))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
            @test (@inferred Tensor{Tuple{2,2}}(data))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
            # different type
            @test (@inferred Tensor{Tuple{2,2}, Float64, 2}(data))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,3.0,4.0)
            @test (@inferred Tensor{Tuple{2,2}, Float64}(data))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,3.0,4.0)
            @test (@inferred Tensor{Tuple{2,2}}((1.0,2,3,4)))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,3.0,4.0)
        end
        @testset "From Vararg" begin
            # same type
            @test (@inferred Tensor{Tuple{2,2}, Int, 2}(data...))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
            @test (@inferred Tensor{Tuple{2,2}, Int}(data...))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
            @test (@inferred Tensor{Tuple{2,2}}(data...))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
            # different type
            @test (@inferred Tensor{Tuple{2,2}, Float64, 2}(data...))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,3.0,4.0)
            @test (@inferred Tensor{Tuple{2,2}, Float64}(data...))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,3.0,4.0)
            @test (@inferred Tensor{Tuple{2,2}}(1.0,2,3,4))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,3.0,4.0)
        end
        @testset "From Function" begin
            # same type
            @test (@inferred Tensor{Tuple{2,2}, Int, 2}((i,j) -> i))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,1,2)
            @test (@inferred Tensor{Tuple{2,2}, Int}((i,j) -> i))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,1,2)
            @test (@inferred Tensor{Tuple{2,2}}((i,j) -> i))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,1,2)
            # different type
            @test (@inferred Tensor{Tuple{2,2}, Float64, 2}((i,j) -> i))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,1.0,2.0)
            @test (@inferred Tensor{Tuple{2,2}, Float64}((i,j) -> i))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,1.0,2.0)
        end
        @testset "From AbstractArray" begin
            A = [1 3; 2 4]
            # same type
            @test (@inferred Tensor{Tuple{2,2}, Int, 2}(A))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
            @test (@inferred Tensor{Tuple{2,2}, Int}(A))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
            @test (@inferred Tensor{Tuple{2,2}}(A))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,2,3,4)
            # different type
            @test (@inferred Tensor{Tuple{2,2}, Float64, 2}(A))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,3.0,4.0)
            @test (@inferred Tensor{Tuple{2,2}, Float64}(A))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,2.0,3.0,4.0)
        end
        @testset "Macro" begin
            a = 2 # check escaping
            @test (@Vec [1,a,3])::Vec{3,Int} == [1,a,3]
            @test (@Vec ones(a))::Vec{a,Float64} == ones(a)
            @test (@Mat [1 a; 3 4])::Mat{2,2,Int} == [1 a; 3 4]
            @test (@Mat ones(a,2))::Mat{a,2,Float64} == ones(a,2)
            @test (@Tensor ones(a,3,3))::Tensor{Tuple{a,3,3},Float64} == ones(a,3,3)
        end
        @testset "zero" begin
            @test (@inferred zero(Tensor{Tuple{2,2}, Int, 2}))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (0,0,0,0)
            @test (@inferred zero(Tensor{Tuple{2,2}, Int}))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (0,0,0,0)
            @test (@inferred zero(Tensor{Tuple{2,2}}))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (0.0,0.0,0.0,0.0)
            for TensorType in (SecondOrderTensor, FourthOrderTensor,
                               SymmetricSecondOrderTensor, SymmetricFourthOrderTensor,
                               Vec)
                @test all(==(0), (@inferred zero(TensorType{2, Int}))::TensorType{2, Int})
                @test all(==(0), (@inferred zero(TensorType{2}))::TensorType{2, Float64})
            end
            @test all(==(0), (@inferred zero(Mat{2, 3, Int}))::Mat{2, 3, Int})
            @test all(==(0), (@inferred zero(Mat{2, 3}))::Mat{2, 3, Float64})
        end
        @testset "ones" begin
            @test (@inferred ones(Tensor{Tuple{2,2}, Int, 2}))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,1,1,1)
            @test (@inferred ones(Tensor{Tuple{2,2}, Int}))::Tensor{Tuple{2,2}, Int, 2, 4} |> Tuple == (1,1,1,1)
            @test (@inferred ones(Tensor{Tuple{2,2}}))::Tensor{Tuple{2,2}, Float64, 2, 4} |> Tuple == (1.0,1.0,1.0,1.0)
            for TensorType in (SecondOrderTensor, FourthOrderTensor,
                               SymmetricSecondOrderTensor, SymmetricFourthOrderTensor,
                               Vec)
                @test all(==(1), (@inferred ones(TensorType{2, Int}))::TensorType{2, Int})
                @test all(==(1), (@inferred ones(TensorType{2}))::TensorType{2, Float64})
            end
            @test all(==(1), (@inferred ones(Mat{2, 3, Int}))::Mat{2, 3, Int})
            @test all(==(1), (@inferred ones(Mat{2, 3}))::Mat{2, 3, Float64})
        end
        @testset "rand/randn" begin
            for op in (rand, randn)
                @test (@inferred op(Tensor{Tuple{2,2}, Float32, 2}))::Tensor{Tuple{2,2}, Float32, 2, 4} |> unique |> length != 1
                @test (@inferred op(Tensor{Tuple{2,2}, Float32}))::Tensor{Tuple{2,2}, Float32, 2, 4} |> unique |> length != 1
                @test (@inferred op(Tensor{Tuple{2,2}}))::Tensor{Tuple{2,2}, Float64, 2, 4} |> unique |> length != 1
                for TensorType in (SecondOrderTensor, FourthOrderTensor,
                                   SymmetricSecondOrderTensor, SymmetricFourthOrderTensor,
                                   Vec)
                    @test (@inferred op(TensorType{2, Float32}))::TensorType{2, Float32} |> unique |> length != 1
                    @test (@inferred op(TensorType{2}))::TensorType{2, Float64} |> unique |> length != 1
                end
                @test (@inferred op(Mat{2, 3, Float32}))::Mat{2, 3, Float32} |> unique |> length != 1
                @test (@inferred op(Mat{2, 3}))::Mat{2, 3, Float64} |> unique |> length != 1
            end
        end
        @testset "aliases" begin
            for T in (Float32, Float64), dim in 1:3
                for TT in (SecondOrderTensor{dim, T},
                           FourthOrderTensor{dim, T},
                           SymmetricSecondOrderTensor{dim, T},
                           SymmetricFourthOrderTensor{dim, T},
                           Mat{dim, dim, T},
                           Vec{dim, T})
                    data = ntuple(i -> T(1), Tensorial.ncomponents(Space(TT)))
                    @test (@inferred TT(data))::TT |> Tuple == data
                    @test (@inferred TT(data...))::TT |> Tuple == data
                end
            end
        end
        @testset "Identity tensor" begin
            # second order tensor
            for TensorType in (SecondOrderTensor, SymmetricSecondOrderTensor)
                for dim in 1:3
                    v = rand(Vec{dim})
                    @test (@inferred one(TensorType{dim, Int}))::TensorType{dim, Int} ⋅ v ≈ v
                    @test (@inferred one(TensorType{dim}))::TensorType{dim, Float64} ⋅ v ≈ v
                end
            end
            # fourth order tensor
            for dim in 1:3
                I = one(SecondOrderTensor{dim})
                A = rand(SecondOrderTensor{dim})
                As = rand(SymmetricSecondOrderTensor{dim})
                II = (@inferred one(FourthOrderTensor{dim, Float32}))::FourthOrderTensor{dim, Float32}
                II = (@inferred one(FourthOrderTensor{dim}))::FourthOrderTensor{dim, Float64}
                @test (II ⊡ A)::SecondOrderTensor{dim} ≈ A
                @test (II ⊡ As)::SecondOrderTensor{dim} ≈ As
                IIs = (@inferred one(SymmetricFourthOrderTensor{dim, Float32}))::SymmetricFourthOrderTensor{dim, Float32}
                IIs = (@inferred one(SymmetricFourthOrderTensor{dim}))::SymmetricFourthOrderTensor{dim, Float64}
                @test (IIs ⊡ A)::SymmetricSecondOrderTensor{dim} ≈ (A+A')/2
                @test (IIs ⊡ As)::SymmetricSecondOrderTensor{dim} ≈ As
            end
        end
        @testset "UniformScaling" begin
            for T in (Float32, Float64)
                @test (@inferred Mat{2, 2}(I))::Mat{2, 2, Bool} == Matrix(I, 2, 2)
                @test (@inferred Mat{2, 2}(2I))::Mat{2, 2, Int} == Matrix(2I, 2, 2)
                @test (@inferred Mat{2, 2, T}(2I))::Mat{2, 2, T} == Matrix(2I, 2, 2)
                @test (@inferred Mat{2, 3, T}(2I))::Mat{2, 3, T} == Matrix(2I, 2, 3)
            end
        end
    end
end

@testset "Symmetric tensors" begin
    x = rand(Tensor{Tuple{2, @Symmetry{2,2}, 3}})
    @test Tensorial.ncomponents(x) == 18
    for i in axes(x, 1), l in axes(x, 4)
        for j in axes(x, 2), k in axes(x, 3)
            @test x[i,j,k,l] == x[i,k,j,l]
        end
    end
    x = rand(Tensor{Tuple{2, @Symmetry{3,3,3}}})
    @test Tensorial.ncomponents(x) == 2 * Tensorial.ncomponents(@Symmetry{3,3,3}())
    for i in axes(x, 1)
        for j in axes(x, 2), k in axes(x, 3), l in axes(x, 4)
            @test x[i,j,k,l] == x[i,k,j,l]
            @test x[i,j,k,l] == x[i,j,l,k]
            @test x[i,j,k,l] == x[i,l,k,j]
        end
    end
end

@testset "Indices" begin
    x = rand(Tensor{Tuple{2, @Symmetry{2,2}, 3, @Symmetry{2,2}}})
    n = Tensorial.ncomponents(x)
    @test (@inferred Tensorial.independent_indices(x))::SArray{Tuple{size(x)...}, Int} |> unique == 1:n
    inds = (@inferred Tensorial.indices(x))::SVector{n, Int}
    @test x[inds] == unique(x[inds])
    dups = (@inferred Tensorial.duplicates(x))::SVector{n, Int}
    for i in eachindex(inds)
        v = x[inds[i]]
        @test count(==(v), x) == dups[i]
    end
end

@testset "Conversion" begin
    @testset "Tensor -> Tensor" begin
        A = Tensor{Tuple{3,3}}(1:9...)
        S = Tensor{Tuple{@Symmetry{3,3}}}(1:6...)
        v = Vec(1:3...)
        for T in (Float32, Float64)
            # convert eltype
            Adata = map(T, (1:9...,))
            @test (@inferred convert(Tensor{Tuple{3,3}, T}, A))::Tensor{Tuple{3,3}, T} |> Tuple == Adata
            @test (@inferred convert(SecondOrderTensor{3, T}, A))::Tensor{Tuple{3,3}, T} |> Tuple == Adata
            @test (@inferred convert(Mat{3, 3, T}, A))::Tensor{Tuple{3,3}, T} |> Tuple == Adata
            Sdata = map(T, (1:6...,))
            @test (@inferred convert(Tensor{Tuple{@Symmetry{3,3}}, T}, S))::Tensor{Tuple{@Symmetry{3,3}}, T} |> Tuple == Sdata
            @test (@inferred convert(SymmetricSecondOrderTensor{3, T}, S))::Tensor{Tuple{@Symmetry{3,3}}, T} |> Tuple == Sdata
            vdata = map(T, (1:3...,))
            @test (@inferred convert(Vec{3, T}, v))::Vec{3, T} |> Tuple == vdata
            # convert symmetric tensor to tensor
            (@inferred convert(Tensor{Tuple{3, 3}, T}, S))::Tensor{Tuple{3, 3}, T} == Array(S)
            (@inferred convert(Mat{3, 3, T}, S))::Tensor{Tuple{3, 3}, T} == Array(S)
            # convert tensor to symmetric tensor
            @test_throws Exception convert(Tensor{Tuple{@Symmetry{3, 3}}, T}, A)
        end
    end
    @testset "AbstractArray -> Tensor" begin
        A = [1 3; 2 4]
        v = [1, 2]
        for T in (Float32, Float64)
            Adata = map(T, (1,2,3,4))
            @test (@inferred convert(Tensor{Tuple{2,2}, T}, A))::Tensor{Tuple{2,2}, T} |> Tuple == Adata
            @test (@inferred convert(SecondOrderTensor{2, T}, A))::Tensor{Tuple{2,2}, T} |> Tuple == Adata
            @test (@inferred convert(Mat{2, 2, T}, A))::Tensor{Tuple{2,2}, T} |> Tuple == Adata
            Sdata = map(T, (1,2,4))
            @test (@inferred convert(Tensor{Tuple{@Symmetry{2,2}}, T}, A))::Tensor{Tuple{@Symmetry{2,2}}, T} |> Tuple == Sdata
            @test (@inferred convert(SymmetricSecondOrderTensor{2, T}, A))::Tensor{Tuple{@Symmetry{2,2}}, T} |> Tuple == Sdata
            vdata = map(T, (1,2))
            @test (@inferred convert(Vec{2, T}, v))::Vec{2, T} |> Tuple == vdata
        end
    end
    @testset "Tuple -> Tensor" begin
        for T in (Float32, Float64)
            Adata = map(T, (1,2,3,4))
            @test (@inferred convert(Vec{4, T}, Adata))::Vec{4, T} |> Tuple == Adata
            @test (@inferred convert(Vec{4}, Adata))::Vec{4, T} |> Tuple == Adata
            @test (@inferred convert(Mat{2, 2, T}, Adata))::Mat{2, 2, T} |> Tuple == Adata
            @test (@inferred convert(Mat{2, 2}, Adata))::Mat{2, 2, T} |> Tuple == Adata
        end
    end
end

@testset "Promotion" begin
    @test (@inferred promote_rule(Vec{2, Float64}, Vec{2, Float32})) === Vec{2, Float64}
    @test (@inferred promote_rule(Mat{2, 2, Int, 4}, Mat{2, 2, Float32, 4})) === Mat{2, 2, Float32, 4}
end

@testset "Tensor misc" begin
    # size
    TT = Tensor{Tuple{2, 3}}
    @test (@inferred size(TT))::Tuple{Int, Int} == (2,3)
    TT = Tensor{Tuple{2, 3, @Symmetry{3, 3}}}
    @test (@inferred size(TT))::Tuple{Int, Int, Int, Int} == (2,3,3,3)
end

@testset "Promotion" begin
    x = Vec(1,2,3)
    for T in (Float32, Float64)
        res = (@inferred Tensorial.promote_elements(x, 3, rand(Mat{3,3,T})))
        @test promote_type(map(eltype, res)...) == T
    end
end


text:
# Synthetic data, not from a particular puzzle
@testset "entirely state abbreviations" begin
    f = best_feature(["wahine", "almond", "invade", "mandarin", "inland", "mainland", "ganymede"])
    @test f.description == "can be completely broken down into US state abbreviations"
end

@testset "state abbreviations" begin
    f = best_feature(["wahinex", "yalmond", "invadez", "wmandar", "inpland", "mainqla"])
    @test f.description == "contains 3 US state abbreviations"
end


text:
#=
New DA structure that allows minimal manipulation
=#

struct EdgeTreeRates2
    tree::Tree                          # the grafted tree
    stem_rate::Array{Float64,1}         # λ_i : every other rates are given divided by λ_i
    tip_rate::Float64                   # the tip rate
    tip_id::Int64                       # tip_id
    tip_rates::Array{Float64,1}         # rates at t_i (for internal branches)
    effective_bl::Float64               # t_j * λ_j / λ_i
    tip_number::Int64                   # extant tips
    scale::Float64
    #relative_rates::Array{Float64,1}
    depth::Float64
    stem_depth::Float64
    parent_edge::Int64
end

function EdgeTreeRates2(tree::Tree, stem_rate::Float64, tip_rate::Float64,
    tip_id::Int64, tip_rates::Array{Float64,1}, tip_number::Int64, depth::Float64,stem_depth::Float64, parent_edge::Int64)
    effective_bl = sum(extract_branch_lengths(tree) .* extract_rates(tree)) / stem_rate
    return EdgeTreeRates2(tree, [stem_rate], tip_rate, tip_id, tip_rates,
        effective_bl, tip_number, depth, stem_depth,parent_edge)
end

function EdgeTreeRates2(tree::Tree, stem_rate::Float64, tip_rate::Float64,
    tip_id::Int64, tip_rates::Array{Float64,1},
    effective_bl::Float64, tip_number::Int64,depth::Float64,stem_depth::Float64,
    parent_edge::Int64)
    return EdgeTreeRates2(tree, [stem_rate], tip_rate, tip_id, tip_rates, effective_bl, tip_number, depth,stem_depth, parent_edge)
end

function EdgeTreeRates2(tree::Tree, stem_rate::Array{Float64,1}, tip_rate::Float64,
    tip_id::Int64, tip_rates::Array{Float64,1},
    effective_bl::Float64, tip_number::Int64,
    depth::Float64,stem_depth::Float64,parent_edge::Int64)

    return EdgeTreeRates2(tree, stem_rate, tip_rate, tip_id, tip_rates, effective_bl,
        tip_number, stem_rate[1], depth,stem_depth, parent_edge)
end

#=
A few functions to manipulate it
=#

function change_bl(tree::Tree, scale::Float64)
    function aux(sub_tree)
        if sub_tree.n_nodes < 2
            return Tree(sub_tree.offsprings,
                sub_tree.branch_length / scale,
                sub_tree.attributes.*scale,
                sub_tree.n_nodes,
                sub_tree.extant)
        else
            left = aux(sub_tree.offsprings[1])
            right = aux(sub_tree.offsprings[2])
            return Tree([left,right],
                sub_tree.branch_length / scale,
                sub_tree.attributes.*scale,
                sub_tree.n_nodes,
                sub_tree.extant)
        end
    end

    aux(tree)
end

function change_bl(tree::Tree, scale::Float64, stem_rate::Float64)
    function aux(sub_tree)
        if sub_tree.n_nodes < 2
            return Tree(sub_tree.offsprings,
                sub_tree.branch_length / scale,
                sub_tree.attributes.*stem_rate,
                sub_tree.n_nodes,
                sub_tree.extant)
        else
            left = aux(sub_tree.offsprings[1])
            right = aux(sub_tree.offsprings[2])
            return Tree([left,right],
                sub_tree.branch_length / scale,
                sub_tree.attributes.*stem_rate,
                sub_tree.n_nodes,
                sub_tree.extant)
        end
    end

    aux(tree)
end

function graft_edge_trees(tree, edge_trees::Array{EdgeTreeRates2,1})

    function aux(subtree, sub_edge_trees)
        if subtree.n_nodes < 2
            return change_bl(sub_edge_trees[1].tree, sub_edge_trees[1].scale, sub_edge_trees[1].stem_rate[1])
        else
            n_edges_left = subtree.offsprings[1].n_nodes + 1
            sub_left = deepcopy(sub_edge_trees[2:n_edges_left])
            sub_right = deepcopy(sub_edge_trees[(n_edges_left+1):subtree.n_nodes])
            subtree_left = aux(subtree.offsprings[1], sub_left)
            subtree_right = aux(subtree.offsprings[2], sub_right)
            return graft_to_edge(Tree([subtree_left, subtree_right], 0., subtree.attributes),
                change_bl(sub_edge_trees[1].tree, sub_edge_trees[1].scale, sub_edge_trees[1].stem_rate[1]), 0, sub_edge_trees[1].tip_id)
        end
    end

    n_edges_left = tree.offsprings[1].n_nodes
    left = deepcopy(edge_trees[1:n_edges_left])
    right = deepcopy(edge_trees[(n_edges_left+1):(tree.n_nodes - 1)])
    tree_left = aux(tree.offsprings[1], left)
    tree_right = aux(tree.offsprings[2], right)
    return Tree([tree_left, tree_right], tree.branch_length, tree.attributes)
end

function node_depths_base(tree::Tree, edge_trees::Array{EdgeTreeRates2,1})
    function aux(subtree, node_times, current_depth, att)
        if subtree.n_nodes < 2
            pushfirst!(node_times, current_depth)
        else
            aux(subtree.offsprings[2],node_times, current_depth + subtree.branch_length/att, att)
            aux(subtree.offsprings[1],node_times, current_depth + subtree.branch_length/att, att)
            pushfirst!(node_times,current_depth)
        end
    end

    node_times = node_depths_base(tree)[2:end]
    node_times_full = Array{Float64,1}(undef,0)
    for i in length(edge_trees):-1:1
        aux(edge_trees[i].tree, node_times_full, node_times[i], edge_trees[i].scale)
    end
    #pushfirst!(node_times,0.)
    return node_times_full
end

function node_depths(tree::Tree, edge_trees::Array{EdgeTreeRates2,1})
    function aux(subtree, node_times, current_depth, sa)
        if subtree.n_nodes < 2
            pushfirst!(node_times, current_depth + subtree.branch_length/sa)
        else
            aux(subtree.offsprings[2],node_times, current_depth + subtree.branch_length/sa, sa)
            aux(subtree.offsprings[1],node_times, current_depth + subtree.branch_length/sa, sa)
            pushfirst!(node_times,current_depth+ subtree.branch_length/sa)
        end
    end

    node_times = node_depths_base(tree)[2:end]
    node_times_full = Array{Float64,1}(undef,0)
    for i in length(edge_trees):-1:1
        aux(edge_trees[i].tree, node_times_full, node_times[i], edge_trees[i].scale)
    end
    #pushfirst!(node_times,0.)
    return node_times_full
end

function node_depths_base(tree::Tree)
    function aux(subtree, node_times, current_depth)
        if subtree.n_nodes < 2
            pushfirst!(node_times, current_depth)
        else
            aux(subtree.offsprings[2],node_times, current_depth + subtree.branch_length)
            aux(subtree.offsprings[1],node_times, current_depth + subtree.branch_length)
            pushfirst!(node_times,current_depth)
        end
    end

    node_times = Array{Float64,1}(undef,0)
    aux(tree, node_times, 0.)
    #pushfirst!(node_times,0.)
    return node_times
end

function node_depths(tree::Tree)
    function aux(subtree, node_times, current_depth)
        if subtree.n_nodes < 2
            pushfirst!(node_times, current_depth + subtree.branch_length)
        else
            aux(subtree.offsprings[2],node_times, current_depth + subtree.branch_length)
            aux(subtree.offsprings[1],node_times, current_depth + subtree.branch_length)
            pushfirst!(node_times,current_depth+ subtree.branch_length)
        end
    end

    node_times = Array{Float64,1}(undef,0)
    aux(tree, node_times, 0.)
    #pushfirst!(node_times,0.)
    return node_times
end

function graft_to_tip(tree,subtree, tip_id)

    function aux(t1, t2, tip)
        if t1.n_nodes == 1
            return Tree(t2.offsprings, t2.branch_length + t1.branch_length, t1.attributes, t2.n_nodes, t2.extant)
        else
            t11 = t1.offsprings[1]
            t12 = t1.offsprings[2]
            t11_tips = (t11.n_nodes + 1) / 2
            if t11_tips < tip
                t12 = aux(t12, t2, tip - t11_tips)
                return Tree([t11, t12], t1.branch_length, t1.attributes)
            else
                t11 = aux(t11, t2, tip)
                return Tree([t11, t12], t1.branch_length, t1.attributes)
            end
        end
    end

    aux(tree,subtree, tip_id)
end

function graft_to_edge(tree, sub_tree, edge, tip)

    if sub_tree.n_nodes < 2
        return Tree(tree.offsprings, tree.branch_length + sub_tree.branch_length, tree.attributes, tree.n_nodes, sub_tree.extant)
    end

    function aux(t1, t2, edge_id)
        if t1.n_nodes == 1
            return t2
        elseif edge_id == 1
            new_t1 = Tree(t1.offsprings, 0., t1.attributes, t1.n_nodes, t1.extant)
            return graft_to_tip(t2, new_t1, tip)
        else
            t11 = t1.offsprings[1]
            t12 = t1.offsprings[2]
            if t11.n_nodes < (edge_id-1)
                t12 = aux(t12, t2, edge_id-1-t11.n_nodes)
                return Tree([t11, t12], t1.branch_length, t1.attributes)
            else
                t11 = aux(t11, t2, edge_id-1)
                return Tree([t11, t12], t1.branch_length, t1.attributes)
            end
        end
    end

    return aux(tree, sub_tree, edge+1)
end

function graft_to_tips(tree, subtree_list, tips)
    if length(tips) == 0
        return tree
    end

    ntips = Int64((tree.n_nodes + 1) / 2)
    subtrees = fill(Tree(),ntips)
    #println("$tips, $ntips")

    for i in 1:ntips
        for j in 1:length(tips)
            if i==tips[j]
                subtrees[i] = subtree_list[j]
                break
            end
        end
    end

    function aux(t1, t2)
        if t1.n_nodes <= 1
            #=if ! t2[1].extant
                println("is it ok ? $(t2[1]) ; $(t1.extant)")
            end=#
            return Tree(t2[1].offsprings, t2[1].branch_length + t1.branch_length, t1.attributes, t2[1].extant)
        else
            t11 = t1.offsprings[1]
            t12 = t1.offsprings[2]
            t11_tips = Int64((t11.n_nodes + 1) / 2)

            t11 = aux(t11, t2[1:t11_tips])
            t12 = aux(t12, t2[(t11_tips+1):end])

            return Tree([t11, t12], t1.branch_length, t1.attributes)
        end
    end

    new_tree = aux(tree,subtrees)
    return new_tree
end

function init_edge_tree_rates2(tree, rates)
    bl = extract_branch_lengths(tree)
    edge_trees = Array{EdgeTreeRates2,1}(undef,0)
    lefts = n_left(tree)
    for i in 2:tree.n_nodes
        nd = get_node_depth(tree, i-1)
        parent_edge = get_parent_edge(tree, i-1)

        ne = 0
        if lefts[i]==0
            ne = 1
        end

        push!(edge_trees,
            EdgeTreeRates2(Tree(Array{Tree,1}(undef,0), bl[i] * rates[i], 1.),       # tree
                rates[i],                                                           # stem_rate
                1.,                                                                 # tip_rate
                1,                                                                  # tip_id
                [1.],                                                               # tip_rates
                bl[i],                                                              # effective_bl
                ne,                                                                  # tip_number
                nd,
                nd + bl[i],
                parent_edge))
    end
    return edge_trees
end

function sim_ClaDS2_time(root_age,σ,α,ε,λ0,u,sf,lf ; return_if_extinct = true, max_node_number = Inf, make_tree = true, prune_extinct = false, return_if_max = true)
    # accesory functions that will be called latter
    bernouilli=Bernoulli(ε/(1+ε))
    function dead()
        return (rand(bernouilli)[1] == 1)
    end

    lambda_law = LogNormal(log(α),σ)
    function new_rates(λ)
        λ * rand(lambda_law, 2)
    end
    #=function new_rates(λ::Float64)
        X = (rand(lambda_law, 2))
        λ * X#* sort2!(X)
    end=#
    #println()
    function aux_ext(time, rate, n_max, s)
        #println(n_max)
        if n_max == 0
            #println("max")
            return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0, -Inf
        end
        node_time = randexp()/(rate*(1+ε))
        if node_time <= 0
            return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0, -Inf
        end
        if node_time > time
            s += lf
            #print(" $s $u ;")
            if u<s
                return Tree(Array{Tree,1}(undef,0), time, [rate], 1, true), n_max + 1, s
            else
                #print(" !")
                return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0, s
            end
        else
                is_dead = dead()
                left_time = time - node_time
                if time == left_time
                    return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0
                end
                if is_dead
                    #print("dead")
                    return Tree(Array{Tree,1}(undef,0), node_time, [rate], 1, false), n_max + 1, s
                else
                    offspring_rates = new_rates(rate)
                    if offspring_rates[1] > offspring_rates[2]
                        offspring_rates = [offspring_rates[2],offspring_rates[1]]
                    end
                    #print("$depth ;")
                    left_tree = aux_ext(left_time, offspring_rates[1], n_max - 1, s)
                    #print(left_tree[2])
                    if left_tree[1].n_nodes == -1
                        return left_tree
                    end
                    right_tree = aux_ext(left_time, offspring_rates[2], left_tree[2], left_tree[3])
                    if right_tree[1].n_nodes == -1
                        return right_tree
                    end
                    return Tree([left_tree[1], right_tree[1]], node_time, [rate]), right_tree[2], right_tree[3]
                end
        end
    end
    tree = aux_ext(root_age, λ0, max_node_number, sf)

    if prune_extinct
        if (tree[1].n_nodes == -1) & return_if_max
            return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), -Inf
        else
            return prune_extinct_lineages(tree[1]), tree[3]
        end
    else
        if (tree[1].n_nodes == -1) & return_if_max
            return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), -Inf
        else
            return tree[1], tree[3]
        end
    end
end

function sim_ClaDS2_time_rates(root_age,σ,α,ε,λ0 ; return_if_extinct = true, max_node_number = 1_000,
    make_tree = true, prune_extinct = false, return_if_max = true, return_na = false)

    # accesory functions that will be called latter
    bernouilli=Bernoulli(ε/(1+ε))
    function dead()
        return (rand(bernouilli)[1] == 1)
    end

    lambda_law = LogNormal(log(α),σ)
    function new_rates(λ)
        λ * rand(lambda_law, 2)
    end
    #=function new_rates(λ::Float64)
        X = (rand(lambda_law, 2))
        λ * X#* sort2!(X)
    end=#
    if return_na
        function aux_na(time, rate, n_max, rates, live)
            #println(n_max)
            if n_max == 0
                return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0
            end
            node_time = randexp()/(rate*(1+ε))
            if node_time <= 0
                return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0
            end
            if node_time > time
                #offsprings, branch_length, attributes, n_nodes
                push!(rates, rate)
                push!(live, true)
                return Tree(Array{Tree,1}(undef,0), time, [rate], 1, true), n_max
            else
                is_dead = dead()
                if is_dead
                    #print("dead")
                    push!(live, false)
                    return Tree(Array{Tree,1}(undef,0), node_time, [rate], 1, false), n_max + 1
                else
                    offspring_rates = new_rates(rate)
                    left_time = time - node_time
                    if time == left_time
                        return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0
                    end
                    #print("$depth ;")
                    left_tree = aux_na(left_time, offspring_rates[1], n_max - 1, rates, live)
                    #print(left_tree[2])
                    if left_tree[1].n_nodes == -1
                        return left_tree
                    end
                    right_tree = aux_na(left_time, offspring_rates[2], left_tree[2], rates, live)
                    if right_tree[1].n_nodes == -1
                        return right_tree
                    end
                    return Tree([left_tree[1], right_tree[1]], node_time, [rate]), right_tree[2]
                end
            end
        end

        rates = Array{Float64,1}(undef,0)
        alive = Array{Bool,1}(undef,0)

        tree = aux_na(root_age, λ0, max_node_number, rates, alive)[1]
        if (tree.n_nodes == -1) & return_if_max
            return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), Array{Float64,1}(undef,0), Array{Bool,1}(undef,0)
        else
            return tree, rates, alive
        end
    else
        function aux_ext(time, rate, n_max, rates)
            #println(n_max)
            if n_max == 0
                return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0
            end
            node_time = randexp()/(rate*(1+ε))
            if node_time > time
                #offsprings, branch_length, attributes, n_nodes
                push!(rates, rate)
                return Tree(Array{Tree,1}(undef,0), time, [rate], 1, true), n_max
            else
                is_dead = dead()
                if is_dead
                    #print("dead")
                    return Tree(Array{Tree,1}(undef,0), node_time, [rate], 1, false), n_max + 1
                else
                    offspring_rates = new_rates(rate)
                    left_time = time - node_time
                    #print("$depth ;")
                    left_tree = aux_ext(left_time, offspring_rates[1], n_max - 1, rates)
                    #print(left_tree[2])
                    if left_tree[1].n_nodes == -1
                        return left_tree
                    end
                    right_tree = aux_ext(left_time, offspring_rates[2], left_tree[2], rates)
                    if right_tree[1].n_nodes == -1
                        return right_tree
                    end
                    return Tree([left_tree[1], right_tree[1]], node_time, [rate]), right_tree[2]
                end
            end
        end

        rates = Array{Float64,1}(undef,0)

        tree = aux_ext(root_age, λ0, max_node_number, rates)[1]
        if (tree.n_nodes == -1) & return_if_max
            return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), Array{Float64,1}(undef,0)
        else
            return tree, rates
        end
    end

end

function sim_ClaDS2_time_unsampled_rates(root_age,σ,α,ε,λ0, u,fl ;
    sampling_proba = 1., return_if_extinct = true, max_node_number = 1_000,
    make_tree = false, prune_extinct = false)

    while true
        tree, rates, s = sim_ClaDS2_time_rates_tip(root_age,σ,α,ε,λ0,u,sampling_proba,fl,
            return_if_extinct = true, max_node_number = max_node_number,
            prune_extinct = false)
        if tree.n_nodes == -1
            return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), false, Inf, rates
        else
            n_alive = length(rates)
            if s < u
                return tree, false, n_alive, rates
            else
                return tree, true, n_alive, rates
            end
        end
    end
end

function sim_ClaDS2_time_rates_tip(root_age,σ,α,ε,λ0,u,sf,fl ; return_if_extinct = true, max_node_number = 1_000,
    make_tree = true, prune_extinct = false, return_if_max = true, return_na = false)

    # accesory functions that will be called latter
    #print(" $(ε/(1+ε))")
    bernouilli=Bernoulli(ε/(1+ε))
    function dead()
        return (rand(bernouilli)[1] == 1)
    end

    lambda_law = LogNormal(log(α),σ)
    function new_rates(λ)
        λ * rand(lambda_law, 2)
    end
    #=function new_rates(λ::Float64)
        X = (rand(lambda_law, 2))
        λ * X#* sort2!(X)
    end=#
    if return_na
    else
        function aux_ext(time, rate, n_max, rates, s ,n)
            if n_max == 0
                return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0,0,0
            end
            node_time = randexp()/(rate*(1+ε))
            #println(node_time)
            if node_time <= 0
                return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0,0,0
            end#print(" $s $n $u;")# $time $node_time $rate;")
            if node_time > time
                #print("zzzz")
                new_s = s
                #offsprings, branch_length, attributes, n_nodes
                push!(rates, rate)
                if n>0
                    new_s = s + log(1-sf)+log(n+1)-log(n)
                else
                    new_s = s + log(sf)
                end
                if new_s >= s
                    #print("a")
                    return Tree(Array{Tree,1}(undef,0), time, [rate], 1, true), n_max, new_s , n+1
                elseif u < new_s
                    #print("b")
                    return Tree(Array{Tree,1}(undef,0), time, [rate], 1, true), n_max, new_s , n+1
                else
                    #print("c")
                    return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0,new_s ,n+1
                end
            else
                is_dead = dead()
                if is_dead
                    return Tree(Array{Tree,1}(undef,0), node_time, [rate], 1, false), n_max + 1, s, n
                else
                    offspring_rates = new_rates(rate)
                    left_time = time - node_time
                    if time == left_time
                        return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), 0,0,0
                    end
                    #print("$depth ;")
                    left_tree = aux_ext(left_time, offspring_rates[1], n_max - 1, rates, s, n)
                    #print(left_tree[2])
                    if left_tree[1].n_nodes == -1
                        return left_tree
                    end
                    right_tree = aux_ext(left_time, offspring_rates[2], left_tree[2], rates, left_tree[3], left_tree[4])
                    if right_tree[1].n_nodes == -1
                        return right_tree
                    end
                    return Tree([left_tree[1], right_tree[1]], node_time, [rate]), right_tree[2], right_tree[3], right_tree[4]
                end
            end
        end

        rates = Array{Float64,1}(undef,0)

        tree = aux_ext(root_age, λ0, max_node_number, rates,fl,0)
        #if ((tree[1].n_nodes == -1) & return_if_max) | (tree[4] == 0) | (u > tree[3])
        #    return Tree(Array{Tree,1}(undef,0), -1., [-1], -1, false), Array{Float64,1}(undef,0)
        #else
        return tree[1], rates, tree[3]
        #end
    end

end

function update_edges_ETR2!(tree::Tree, edge_trees::Array{EdgeTreeRates2,1}, σ::Float64, α::Float64, ε::Float64, rates::Array{Float64,1},
    lefts::Array{Int64,1}; it_rates = 1, prior_ε = "lognormal", logε0 = 0., sd = 0.5)

    if prior_ε == "uniform"
        ε = draw_ε_crown_priorUnif(tree, edge_trees, lefts)
    elseif prior_ε == "uniformInf"
        ε = draw_ε_crown(tree, edge_trees, lefts)
    elseif prior_ε == "ClaDS0"
        ε = 0.
    else
        ε = draw_ε_crown_priorln(tree, edge_trees, lefts, logε0 = logε0, sd = sd)
    end

    for j in 1:it_rates

        draw_λi_quad!(rates, edge_trees, σ, α, ε, tree, lefts)
        relative_rates = extract_relative_rates(tree, edge_trees, rates)

        σ = draw_σ(relative_rates, α, β0 = 0.05, α0 = 0.5)
        α = draw_α(relative_rates, σ,  α_0 = -0.05, σ_0 = 0.1)

        if prior_ε == "uniform"
            ε = draw_ε_crown_priorUnif(tree, edge_trees, lefts)
        elseif prior_ε == "uniformInf"
            ε = draw_ε_crown(tree, edge_trees, lefts)
        elseif prior_ε == "ClaDS0"
            ε = 0.
        else
            ε = draw_ε_crown_priorln(tree, edge_trees, lefts, logε0 = logε0, sd = sd)
        end
    end

    return ε, σ, α
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
function constraint_converter_losses(pm::_PM.AbstractLPACModel, n::Int, i::Int, a, b, c, plmax)
    pconv_ac = _PM.var(pm, n, :pconv_ac, i)
    pconv_dc = _PM.var(pm, n, :pconv_dc, i)
    iconv = _PM.var(pm, n, :iconv_ac, i)

    JuMP.@constraint(pm.model, pconv_ac + pconv_dc == a + b*iconv)
end

function constraint_conv_transformer(pm::_PM.AbstractLPACModel, n::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    ptf_fr = _PM.var(pm, n, :pconv_tf_fr, i)
    qtf_fr = _PM.var(pm, n, :qconv_tf_fr, i)
    ptf_to = _PM.var(pm, n, :pconv_tf_to, i)
    qtf_to = _PM.var(pm, n, :qconv_tf_to, i)

    phi = _PM.var(pm, n, :phi, acbus)
    phi_vmf = _PM.var(pm, n, :phi_vmf, i)
    va = _PM.var(pm, n, :va, acbus)
    vaf = _PM.var(pm, n, :vaf, i)
    cs = _PM.var(pm, n, :cs_vaf,i)

    ztf = rtf + im*xtf
    if transformer
        ytf = 1/(rtf + im*xtf)
        gtf = real(ytf)
        btf = imag(ytf)
        c1, c2, c3, c4 = lpac_power_flow_constraints(pm.model, gtf, btf, phi, phi_vmf, va, vaf, ptf_fr, ptf_to, qtf_fr, qtf_to, tm, cs)
        c5 = constraint_cos_angle_diff_PWL(pm, n, cs, va, vaf)
    else
        JuMP.@constraint(pm.model, ptf_fr + ptf_to == 0)
        JuMP.@constraint(pm.model, qtf_fr + qtf_to == 0)
        JuMP.@constraint(pm.model, va == vaf)
        JuMP.@constraint(pm.model, (1+phi) == (1+phi_vmf))
    end
end

"constraints for a voltage magnitude transformer + series impedance"

function lpac_power_flow_constraints(model, g, b, phi_fr, phi_to, va_fr, va_to, p_fr, p_to, q_fr, q_to, tm, cs)

    c1 = JuMP.@constraint(model, p_fr ==  g/(tm^2)*(1.0 + 2*phi_fr) + (-g/tm)*(cs + phi_fr + phi_to) + (-b/tm)*(va_fr-va_to))
    c2 = JuMP.@constraint(model, q_fr == -b/(tm^2)*(1.0 + 2*phi_fr) - (-b/tm)*(cs + phi_fr + phi_to) + (-g/tm)*(va_fr-va_to))
    c3 = JuMP.@constraint(model, p_to ==  g*(1.0 + 2*phi_to) + (-g/tm)*(cs + phi_fr + phi_to) + (-b/tm)*-(va_fr-va_to))
    c4 = JuMP.@constraint(model, q_to == -b*(1.0 + 2*phi_to) - (-b/tm)*(cs + phi_fr + phi_to) + (-g/tm)*-(va_fr-va_to))
    return c1, c2, c3, c4
end

function constraint_conv_reactor(pm::_PM.AbstractLPACModel, n::Int, i::Int, rc, xc, reactor)
    pconv_ac  = _PM.var(pm, n, :pconv_ac, i)
    qconv_ac  = _PM.var(pm, n, :qconv_ac, i)
    ppr_to = - pconv_ac
    qpr_to = - qconv_ac
    ppr_fr = _PM.var(pm, n, :pconv_pr_fr, i)
    qpr_fr = _PM.var(pm, n, :qconv_pr_fr, i)

    phi_vmc = _PM.var(pm, n, :phi_vmc, i)
    phi_vmf = _PM.var(pm, n, :phi_vmf, i)
    vac = _PM.var(pm, n, :vac, i)
    vaf = _PM.var(pm, n, :vaf, i)
    cs = _PM.var(pm, n, :cs_vac,i)

    phi_vmc_ub = JuMP.upper_bound(phi_vmc)
    ppr_to_ub = JuMP.upper_bound(_PM.var(pm, n)[:pconv_ac][i])
    qpr_to_ub = JuMP.upper_bound(_PM.var(pm, n)[:qconv_ac][i])
    Smax = sqrt(ppr_to_ub^2 + qpr_to_ub^2)
    zc = rc + im*xc
    if reactor
        yc = 1/(zc)
        gc = real(yc)
        bc = imag(yc)
        c1, c2, c3, c4 = lpac_power_flow_constraints(pm.model, gc, bc, phi_vmf, phi_vmc, vaf, vac, ppr_fr, ppr_to, qpr_fr, qpr_to, 1, cs)
        c5 = constraint_cos_angle_diff_PWL(pm, n, cs, vaf, vac)
        c6 = constraint_conv_capacity_PWL(pm, n, ppr_to, qpr_to, ppr_to_ub, qpr_to_ub, Smax)
   else
        JuMP.@constraint(pm.model, ppr_fr + ppr_to == 0)
        JuMP.@constraint(pm.model, qpr_fr + qpr_to == 0)
        JuMP.@constraint(pm.model, vac == vaf)
        JuMP.@constraint(pm.model, (1+phi_vmf) == (1+phi_vmc))

    end
end

function constraint_conv_filter(pm::_PM.AbstractLPACModel, n::Int, i::Int, bv, filter)
    ppr_fr = _PM.var(pm, n, :pconv_pr_fr, i)
    qpr_fr = _PM.var(pm, n, :qconv_pr_fr, i)
    ptf_to = _PM.var(pm, n, :pconv_tf_to, i)
    qtf_to = _PM.var(pm, n, :qconv_tf_to, i)
    phi_vmf = _PM.var(pm, n, :phi_vmf, i)

    JuMP.@constraint(pm.model, ppr_fr + ptf_to == 0 )
    JuMP.@constraint(pm.model, qpr_fr + qtf_to + -bv*filter*(1+2*phi_vmf) == 0)
end

function constraint_converter_current(pm::_PM.AbstractLPACModel, n::Int, i::Int, Umax, Imax)
    phi_vmc = _PM.var(pm, n, :phi_vmc, i)
    pconv_ac = _PM.var(pm, n, :pconv_ac, i)
    qconv_ac = _PM.var(pm, n, :qconv_ac, i)
    iconv = _PM.var(pm, n, :iconv_ac, i)

    JuMP.@constraint(pm.model, iconv <= Imax)
end

function variable_converter_filter_voltage(pm::_PM.AbstractLPACModel; kwargs...)
    variable_converter_filter_voltage_magnitude(pm; kwargs...)
    variable_converter_filter_voltage_angle_cs(pm; kwargs...)
    variable_converter_filter_voltage_angle(pm; kwargs...)
end

function variable_converter_internal_voltage(pm::_PM.AbstractLPACModel; kwargs...)
    variable_converter_internal_voltage_magnitude(pm; kwargs...)
    variable_converter_internal_voltage_angle_cs(pm; kwargs...)
    variable_converter_internal_voltage_angle(pm; kwargs...)
end

function variable_converter_filter_voltage_angle_cs(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded = true, report = true)
    csvaf = _PM.var(pm, nw)[:cs_vaf] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_cs_vaf",
    start = 0
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(csvaf[c],  0)
            JuMP.set_upper_bound(csvaf[c],  1)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :cs_vaf, _PM.ids(pm, nw, :convdc), csvaf)
end

function variable_converter_internal_voltage_angle_cs(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded = true, report = true)
    csvac = _PM.var(pm, nw)[:cs_vac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_cs_vac",
    start = 0
    )

    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(csvac[c],  0)
            JuMP.set_upper_bound(csvac[c],  1)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :cs_vac, _PM.ids(pm, nw, :convdc), csvac)
end

function variable_converter_filter_voltage_magnitude(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded = true, report = true)
    phivmf = _PM.var(pm, nw)[:phi_vmf] = JuMP.@variable(pm.model,
            [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_phi_vmf",
            start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "phi_start")
        )

        if bounded
            for (c, convdc) in _PM.ref(pm, nw, :convdc)
                JuMP.set_lower_bound(phivmf[c],  convdc["Vmmin"] - 1.0)
                JuMP.set_upper_bound(phivmf[c],  convdc["Vmmax"] - 1.0)
            end
        end

        report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :phi_vmf, _PM.ids(pm, nw, :convdc), phivmf)
end

function variable_converter_internal_voltage_magnitude(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded = true, report = true)
    phivmc = _PM.var(pm, nw)[:phi_vmc] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_phi_vmc",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "phi_start")
        )

        if bounded
            for (c, convdc) in _PM.ref(pm, nw, :convdc)
                JuMP.set_lower_bound(phivmc[c],  convdc["Vmmin"] - 1.0)
                JuMP.set_upper_bound(phivmc[c],  convdc["Vmmax"] - 1.0)
            end
        end

        report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :phi_vmc, _PM.ids(pm, nw, :convdc), phivmc)
end

function constraint_conv_capacity_PWL(pm::_PM.AbstractLPACModel, n::Int, ppr_to, qpr_to, Umax, Imax, Smax)
    np = 20 #no. of segments, can be passed as an argument later
    l = 0
    for i = 1:np
        a= Smax*sin(l)
        b = Smax*cos(l)
        c6 = JuMP.@constraint(pm.model, a*ppr_to + b*qpr_to <= Smax^2) #current and voltage bounds to be proper to use Umax*Imax because Umax*Imax == Smax
        l = l + 2*pi/np
    end
end

function constraint_cos_angle_diff_PWL(pm::_PM.AbstractLPACModel, n::Int, cs, va_fr, va_to)
    nb = 20 #no. of segments, can be passed as an argument later
    l = -pi/6
    h = pi/6
    inc = (h-l)/(nb+1)
    a = l + inc
    diff = va_fr - va_to
    for i = 1:nb
        c5 = JuMP.@constraint(pm.model, cs <= -sin(a)*(diff-a) + cos(a))
        a = a + inc
    end
end

function add_dcconverter_voltage_setpoint(sol, pm::_PM.AbstractLPACModel)
    _PM.add_setpoint!(sol, pm, "convdc", "vmconv", :phi_vmc, status_name="islcc", inactive_status_value = 4, scale = (x,item,cnd) -> 1.0+x)
    _PM.add_setpoint!(sol, pm, "convdc", "vmfilt", :phi_vmf, status_name="islcc", inactive_status_value = 4, scale = (x,item,cnd) -> 1.0+x)
    _PM.add_setpoint!(sol, pm, "convdc", "vaconv", :vac, status_name="islcc", inactive_status_value = 4)
    _PM.add_setpoint!(sol, pm, "convdc", "vafilt", :vaf, status_name="islcc", inactive_status_value = 4)
end

############ TNEP constraints ###############
function constraint_converter_losses_ne(pm::_PM.AbstractLPACModel, n::Int, i::Int, a, b, c, plmax)
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne, i)
    pconv_dc = _PM.var(pm, n, :pconv_dc_ne, i)
    iconv = _PM.var(pm, n, :iconv_ac_ne, i)
    z = _PM.var(pm, n, :conv_ne, i)

    JuMP.@constraint(pm.model, pconv_ac + pconv_dc == a*z + b*iconv)
end

function constraint_conv_transformer_ne(pm::_PM.AbstractLPACModel, n::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    ptf_fr = _PM.var(pm, n, :pconv_tf_fr_ne, i)
    qtf_fr = _PM.var(pm, n, :qconv_tf_fr_ne, i)
    ptf_to = _PM.var(pm, n, :pconv_tf_to_ne, i)
    qtf_to = _PM.var(pm, n, :qconv_tf_to_ne, i)

    phi = _PM.var(pm, n, :phi, acbus)
    phi_du = _PM.var(pm, n, :phi_du, i)
    phi_vmf = _PM.var(pm, n, :phi_vmf_ne, i)
    va = _PM.var(pm, n, :va, acbus)
    va_du = _PM.var(pm, n, :va_du, i)
    JuMP.set_upper_bound(va, pi/2)
    JuMP.set_lower_bound(va, -pi/2)
    vaf = _PM.var(pm, n, :vaf_ne, i)
    cs = _PM.var(pm, n, :cs_vaf_ne,i)
    z = _PM.var(pm, n, :conv_ne)[i]

    ztf = rtf + im*xtf
    if transformer
        ytf = 1/(rtf + im*xtf)
        gtf = real(ytf)
        btf = imag(ytf)
        c1, c2, c3, c4 = lpac_power_flow_constraints(pm.model, gtf, btf, phi_du, phi_vmf, va_du, vaf, ptf_fr, ptf_to, qtf_fr, qtf_to, tm, cs, z)
        c5 = constraint_cos_angle_diff_PWL(pm, n, cs, va, vaf)
        relaxation_semicont_variable_on_off(pm.model, cs, z) # cos not yet semicont. but once put ub and lb acc. to angle difference, it can be semi-cont
    else
        JuMP.@constraint(pm.model, ptf_fr + ptf_to == 0)
        JuMP.@constraint(pm.model, qtf_fr + qtf_to == 0)
        JuMP.@constraint(pm.model, va_du == vaf)
        JuMP.@constraint(pm.model, (1*z+phi_du) == (1*z+phi_vmf))
    end

    # relaxation_variable_on_off(pm.model, phi, phi_du, z)
    _IM.relaxation_equality_on_off(pm.model, phi, phi_du, z)

    # relaxation_variable_on_off(pm.model, va, va_du, z)
    _IM.relaxation_equality_on_off(pm.model, va, va_du, z)

    relaxation_semicont_variable_on_off(pm.model, phi_vmf, z)
    relaxation_semicont_variable_on_off(pm.model, vaf, z)
end

"constraints for a voltage magnitude transformer + series impedance"

function lpac_power_flow_constraints(model, g, b, phi_fr, phi_to, va_fr, va_to, p_fr, p_to, q_fr, q_to, tm, cs, z)
    c1 = JuMP.@constraint(model, p_fr ==  g/(tm^2)*(1.0*z + 2*phi_fr) + (-g/tm)*(cs + phi_fr + phi_to) + (-b/tm)*(va_fr-va_to))
    c2 = JuMP.@constraint(model, q_fr == -b/(tm^2)*(1.0*z + 2*phi_fr) - (-b/tm)*(cs + phi_fr + phi_to) + (-g/tm)*(va_fr-va_to))
    c3 = JuMP.@constraint(model, p_to ==  g*(1.0*z + 2*phi_to) + (-g/tm)*(cs + phi_fr + phi_to) + (-b/tm)*-(va_fr-va_to))
    c4 = JuMP.@constraint(model, q_to == -b*(1.0*z + 2*phi_to) - (-b/tm)*(cs + phi_fr + phi_to) + (-g/tm)*-(va_fr-va_to))
    return c1, c2, c3, c4
end

function constraint_conv_reactor_ne(pm::_PM.AbstractLPACModel, n::Int, i::Int, rc, xc, reactor)
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne, i)
    qconv_ac = _PM.var(pm, n, :qconv_ac_ne, i)
    ppr_to = - pconv_ac
    qpr_to = - qconv_ac
    ppr_fr = _PM.var(pm, n, :pconv_pr_fr_ne, i)
    qpr_fr = _PM.var(pm, n, :qconv_pr_fr_ne, i)

    phi_vmc = _PM.var(pm, n, :phi_vmc_ne, i)
    phi_vmf = _PM.var(pm, n, :phi_vmf_ne, i)
    vac = _PM.var(pm, n, :vac_ne, i)
    vaf = _PM.var(pm, n, :vaf_ne, i)
    cs = _PM.var(pm, n, :cs_vac_ne,i)
    z = _PM.var(pm, n, :conv_ne)[i]

    phi_vmc_ub = JuMP.upper_bound(phi_vmc)
    ppr_to_ub = JuMP.upper_bound(_PM.var(pm, n)[:pconv_ac_ne][i])
    qpr_to_ub = JuMP.upper_bound(_PM.var(pm, n)[:qconv_ac_ne][i])
    Smax = sqrt(ppr_to_ub^2 + qpr_to_ub^2)
    zc = rc + im*xc
    if reactor
        yc = 1/(zc)
        gc = real(yc)
        bc = imag(yc)
        c1, c2, c3, c4 = lpac_power_flow_constraints(pm.model, gc, bc, phi_vmf, phi_vmc, vaf, vac, ppr_fr, ppr_to, qpr_fr, qpr_to, 1, cs, z)
        c5 = constraint_cos_angle_diff_PWL(pm, n, cs, vaf, vac)
        relaxation_semicont_variable_on_off(pm.model, cs, z) # cos not yet semicont. but once put ub and lb acc. to angle difference, it can be semi-cont
    else
        JuMP.@constraint(pm.model, ppr_fr + ppr_to == 0)
        JuMP.@constraint(pm.model, qpr_fr + qpr_to == 0)
        JuMP.@constraint(pm.model, vac == vaf)
        JuMP.@constraint(pm.model, (1*z+phi_vmf) == (1*z+phi_vmc))
   end
    c6 = constraint_conv_capacity_PWL(pm, n, ppr_to, qpr_to, ppr_to_ub, qpr_to_ub, Smax)
    relaxation_semicont_variable_on_off(pm.model, phi_vmc, z)
    relaxation_semicont_variable_on_off(pm.model, vac, z)
end

function constraint_conv_filter_ne(pm::_PM.AbstractLPACModel, n::Int, i::Int, bv, filter)
    ppr_fr = _PM.var(pm, n, :pconv_pr_fr_ne, i)
    qpr_fr = _PM.var(pm, n, :qconv_pr_fr_ne, i)
    ptf_to = _PM.var(pm, n, :pconv_tf_to_ne, i)
    qtf_to = _PM.var(pm, n, :qconv_tf_to_ne, i)
    phi_vmf = _PM.var(pm, n, :phi_vmf_ne, i)
    z = _PM.var(pm, n, :conv_ne)[i]

    JuMP.@constraint(pm.model, ppr_fr + ptf_to == 0 )
    JuMP.@constraint(pm.model, qpr_fr + qtf_to + -bv*filter*(1*z+2*phi_vmf) == 0)
end

function constraint_converter_current_ne(pm::_PM.AbstractLPACModel, n::Int, i::Int, Umax, Imax)
    phi_vmc = _PM.var(pm, n, :phi_vmc_ne, i)
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne, i)
    qconv_ac = _PM.var(pm, n, :qconv_ac_ne, i)
    iconv = _PM.var(pm, n, :iconv_ac_ne, i)
    # @NLconstraint(pm.model, pconv_ac^2 + qconv_ac^2 <= (1.0+phi_vmc)^2 * (iconv)^2)
    JuMP.@constraint(pm.model, iconv <= Imax)
end

function variable_converter_filter_voltage_ne(pm::_PM.AbstractLPACModel; kwargs...)
    variable_converter_filter_voltage_magnitude_ne(pm; kwargs...)
    variable_converter_filter_voltage_angle_cs_ne(pm; kwargs...)
    variable_converter_filter_voltage_angle_ne(pm; kwargs...)
end

function variable_converter_internal_voltage_ne(pm::_PM.AbstractLPACModel; kwargs...)
    variable_converter_internal_voltage_magnitude_ne(pm; kwargs...)
    variable_converter_internal_voltage_angle_cs_ne(pm; kwargs...)
    variable_converter_internal_voltage_angle_ne(pm; kwargs...)
end

function variable_converter_filter_voltage_angle_cs_ne(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded = true)
    _PM.var(pm, nw)[:cs_vaf_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_cs_vaf_ne",
    lower_bound =  0, # using the maximum bound allowed in cosine linearization. can be changed to bus angle +/- 10 degree?
    upper_bound =  1,
    start = 0
    )
end

function variable_converter_internal_voltage_angle_cs_ne(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded = true)
    _PM.var(pm, nw)[:cs_vac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_cs_vac_ne",
    lower_bound =  0, # using the maximum bound allowed in cosine linearization. can be changed to bus angle +/- 10 degree?
    upper_bound =  1,
    start = 0
    )
end

function variable_converter_filter_voltage_magnitude_ne(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded = true)
    _PM.var(pm, nw)[:phi_vmf_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_phi_vmf_ne",
    lower_bound = _PM.ref(pm, nw, :convdc_ne, i, "Vmmin") - 1.0,
    upper_bound = _PM.ref(pm, nw, :convdc_ne, i, "Vmmax") - 1.0,
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "phi_start")
    )
end

function variable_converter_internal_voltage_magnitude_ne(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded = true)
    _PM.var(pm, nw)[:phi_vmc_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_phi_vmc_ne",
    lower_bound = _PM.ref(pm, nw, :convdc_ne, i, "Vmmin") - 1.0,
    upper_bound = _PM.ref(pm, nw, :convdc_ne, i, "Vmmax") - 1.0,
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "phi_start")
    )
end

function variable_voltage_slack(pm::_PM.AbstractLPACModel; nw::Int=_PM.nw_id_default, bounded::Bool = true)
    _PM.var(pm, nw)[:phi_du] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_phi_du",
    lower_bound = -0.2,
    upper_bound = 0.2,
    start = 0,
    )
    _PM.var(pm, nw)[:va_du] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_va_du",
    lower_bound = -2*pi,
    upper_bound = 2*pi,
    start = 0,
    )
end

function add_dcconverter_voltage_setpoint_ne(sol, pm::_PM.AbstractLPACModel)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "vmconv", :phi_vmc_ne, status_name="islcc", inactive_status_value = 4, scale = (x,item,cnd) -> 1.0+x)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "vmfilt", :phi_vmf_ne, status_name="islcc", inactive_status_value = 4, scale = (x,item,cnd) -> 1.0+x)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "vaconv", :vac_ne, status_name="islcc", inactive_status_value = 4)
    _PM.add_setpoint!(sol, pm, "convdc_ne", "vafilt", :vaf_ne, status_name="islcc", inactive_status_value = 4)

end


text:
### Day 6 - Problem 1 & 2
## Signals and Noise
## Author: Thanasis Georgiou

f = open("input.txt")
input = readlines(f)
close(f)

# Find length
packetLen = length(chomp(input[1]))

# Create an array of hashmaps
packet = Array(Dict{Char, Int}, packetLen)
map!(d -> Dict{Char, Int}(), packet)

# Read every line
for line in input
    arr = collect(chomp(line))
    for i = 1:length(arr)
        lastValue = get(packet[i], arr[i], -1)
        if lastValue == -1
            packet[i][arr[i]] = 1
        else
            packet[i][arr[i]] += 1
        end
    end
end

# Find max for each hashmap
sorter(a, b) = last(a) < last(b)

resultFalse = map(d -> sort(collect(d), lt = sorter, rev = true)[1][1], packet)
resultTrue = map(d -> sort(collect(d), lt = sorter, rev = false)[1][1], packet)

println("False packet: $(resultFalse)")
println("True packet: $(resultTrue)")


text:
"""
    ComputeWeightFunction(z::Float64, i::Int64, ConvolvedDensity::AbstractConvolvedDensity,
    AnalitycalDensity::AnalitycalDensity, InstrumentResponse::InstrumentResponse,
    w0waCDMCosmology::w0waCDMCosmology, GCWeightFunction::GCWeightFunction)

This function returns the Galaxy Clustering Weight function for a given redshift
``z`` and tomographic bin ``i``.
"""
function ComputeWeightFunction(z::Float64, i::Int64,
    ConvolvedDensity::AbstractConvolvedDensity, AnalitycalDensity::AnalitycalDensity,
    InstrumentResponse::InstrumentResponse, Cosmology::AbstractCosmology,
    GCWeightFunction::GCWeightFunction)
    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                       #physical constants involved in calculations
    return ComputeConvolvedDensity(z, i, ConvolvedDensity,
    AnalitycalDensity, InstrumentResponse) *
    ComputeBias(z, GCWeightFunction.BiasKind, ConvolvedDensity) *
    ComputeHubbleFactor(z, Cosmology) / c_0
end

"""
    ComputeWeightFunctionGrid!(GCWeightFunction::GCWeightFunction,
    ConvolvedDensity::AbstractConvolvedDensity, CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities, Cosmology::AbstractCosmology)

This function evaluates the Galaxy Clustering Weight Function over the ``z``
grid and for all tomographic bins ``i``.
"""
function ComputeWeightFunctionGrid!(GCWeightFunction::GCWeightFunction,
    ConvolvedDensity::AbstractConvolvedDensity, CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities, Cosmology::AbstractCosmology)
    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                       #physical constants involved in calculations
    #TODO use avx here to accelerate
    for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
        for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
            GCWeightFunction.WeightFunctionArray[idx_ZBinArray, idx_ZArray] =
            GCWeightFunction.BiasArray[idx_ZBinArray, idx_ZArray]*
            ConvolvedDensity.DensityGridArray[idx_ZBinArray,
            idx_ZArray] *
            BackgroundQuantities.HZArray[idx_ZArray] / c_0
        end
    end
end

"""
    ComputeLensingEfficiency(z::Float64, i::Int64,
    ConvolvedDensity::AbstractConvolvedDensity, AnalitycalDensity::AnalitycalDensity,
    InstrumentResponse::InstrumentResponse, w0waCDMCosmology::w0waCDMCosmology,
    CosmologicalGrid::CosmologicalGrid, ::WLWeightFunction)

This function returns the Lensing efficiency, for a given redshift
``z`` and tomographic bin ``i``.
"""
function ComputeLensingEfficiency(z::Float64, i::Int64,
    ConvolvedDensity::AbstractConvolvedDensity,
    AnalitycalDensity::AnalitycalDensity,
    InstrumentResponse::InstrumentResponse,
    Cosmology::AbstractCosmology,
    CosmologicalGrid::CosmologicalGrid,
    ::WLWeightFunction)
    int, err = QuadGK.quadgk(x -> CosmoCentral.ComputeConvolvedDensity(
    x, i, ConvolvedDensity, AnalitycalDensity, InstrumentResponse)*
    ((Computeχ(x, Cosmology) -
    Computeχ(z, Cosmology))/
    (Computeχ(x, Cosmology) *
    Computeχ(z, Cosmology)  ) ), z,
    last(CosmologicalGrid.ZArray) , rtol=1e-12)
    return int
end

function ComputeLensingEfficiency(z::Float64, i::Int64,
    ConvolvedDensity::AbstractConvolvedDensity,
    AnalitycalDensity::AnalitycalDensity,
    InstrumentResponse::InstrumentResponse,
    Cosmology::AbstractCosmology,
    CosmologicalGrid::CosmologicalGrid,
    LensingSourceFunction::LensingSourceFunction)
    int, err = QuadGK.quadgk(x -> CosmoCentral.ComputeConvolvedDensity(
    x, i, ConvolvedDensity, AnalitycalDensity, InstrumentResponse)*
    ((Computeχ(x, Cosmology) -
    Computeχ(z, Cosmology))/
    (Computeχ(x, Cosmology) *
    Computeχ(z, Cosmology)  ) ), z,
    last(CosmologicalGrid.ZArray) , rtol=1e-12)
    return int
end

"""
    ComputeLensingEfficiencyGrid!(LensingFunction::WLWeightFunction,
    AnalitycalDensity::AnalitycalDensity, InstrumentResponse::InstrumentResponse,
    ConvolvedDensity::AbstractConvolvedDensity, CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities, Cosmology::AbstractCosmology,
    ::StandardLensingEfficiency)

This function evaluates the Lensing Efficiency over the ``z``
grid and for all tomographic bins ``i``.
"""
function ComputeLensingEfficiencyGrid!(LensingFunction::WLWeightFunction,
    AnalitycalDensity::AnalitycalDensity, InstrumentResponse::InstrumentResponse,
    ConvolvedDensity::AbstractConvolvedDensity, CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities, Cosmology::AbstractCosmology,
    ::StandardLensingEfficiency)
    for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
        for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
            LensingFunction.LensingEfficiencyArray[idx_ZBinArray, idx_ZArray] =
            ComputeLensingEfficiency(CosmologicalGrid.ZArray[idx_ZArray],
            idx_ZBinArray, ConvolvedDensity, AnalitycalDensity,
            InstrumentResponse, Cosmology, CosmologicalGrid,
            LensingFunction)
        end
    end
end

function ComputeLensingEfficiencyGrid!(
    LensingSourceFunction::LensingSourceFunction,
    AnalitycalDensity::AnalitycalDensity,
    InstrumentResponse::InstrumentResponse,
    ConvolvedDensity::AbstractConvolvedDensity,
    CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities,
    Cosmology::AbstractCosmology,
    ::StandardLensingEfficiency)
    for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
        for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
            LensingSourceFunction.LensingEfficiencyArray[idx_ZBinArray, idx_ZArray] =
            ComputeLensingEfficiency(CosmologicalGrid.ZArray[idx_ZArray],
            idx_ZBinArray, ConvolvedDensity, AnalitycalDensity,
            InstrumentResponse, Cosmology, CosmologicalGrid,
            LensingSourceFunction)
        end
    end
end

"""
    ComputeLensingEfficiencyGrid!(
    LensingFunction::WLWeightFunction, ConvolvedDensity::AbstractConvolvedDensity,
    CosmologicalGrid::CosmologicalGrid, BackgroundQuantities::BackgroundQuantities,
    w0waCDMCosmology::w0waCDMCosmology, ::CustomLensingEfficiency)

This function evaluates the Lensing Efficiency over the ``z`` grid and for all tomographic
bins ``i``.
"""
function ComputeLensingEfficiencyGrid!(
    LensingFunction::WLWeightFunction,
    ConvolvedDensity::AbstractConvolvedDensity,
    CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities,
    Cosmology::AbstractCosmology,
    ::CustomLensingEfficiency)
    LensingFunction.LensingEfficiencyArray .*= 0
    Weight_Matrix = SimpsonWeightMatrix(length(CosmologicalGrid.ZArray))
    @avx for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
        for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
            for idx_ZArrayInt in 1:length(CosmologicalGrid.ZArray)
                LensingFunction.LensingEfficiencyArray[idx_ZBinArray,
                idx_ZArray] += ConvolvedDensity.DensityGridArray[idx_ZBinArray,
                idx_ZArrayInt] * (BackgroundQuantities.rZArray[idx_ZArrayInt] -
                BackgroundQuantities.rZArray[idx_ZArray]) /
                BackgroundQuantities.rZArray[idx_ZArrayInt] /
                BackgroundQuantities.rZArray[idx_ZArray] *
                Weight_Matrix[idx_ZArray, idx_ZArrayInt]
            end
        end
    end
    LensingFunction.LensingEfficiencyArray .*=
    (last(CosmologicalGrid.ZArray)-first(CosmologicalGrid.ZArray))/
    (length(CosmologicalGrid.ZArray)-1)
end

"""
    ComputeWeightFunction(z::Float64, i::Int64,ConvolvedDensity::AbstractConvolvedDensity,
    AnalitycalDensity::AnalitycalDensity, InstrumentResponse::InstrumentResponse,
    Cosmology::AbstractCosmology, CosmologicalGrid::CosmologicalGrid,
    LensingFunction::WLWeightFunction)

This function returns the Weak Lensing Weight Function, for a given redshift
``z`` and tomographic bin ``i``.
"""
function ComputeWeightFunction(z::Float64, i::Int64,
    ConvolvedDensity::AbstractConvolvedDensity, AnalitycalDensity::AnalitycalDensity,
    InstrumentResponse::InstrumentResponse, Cosmology::AbstractCosmology,
    CosmologicalGrid::CosmologicalGrid, LensingFunction::WLWeightFunction)
    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                       #physical constants involved in calculations
    return 1.5 * ComputeLensingEfficiency(z, i, ConvolvedDensity,
    AnalitycalDensity,InstrumentResponse, Cosmology, CosmologicalGrid,
    LensingFunction) * (Cosmology.H0/c_0)^2 * Cosmology.ΩM *
    (1. + z) * Computeχ(z, Cosmology)^2
end

"""
    ComputeWeightFunctionGrid!(LensingFunction::WLWeightFunction,
    ConvolvedDensity::AbstractConvolvedDensity, CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities, Cosmology::AbstractCosmology)

This function evaluates the Weak Lensing Weight Function over the ``z``
grid and for all tomographic bins ``i``.
"""
function ComputeWeightFunctionGrid!(
    LensingFunction::WLWeightFunction,
    ConvolvedDensity::AbstractConvolvedDensity,
    CosmologicalGrid::CosmologicalGrid,
    BackgroundQuantities::BackgroundQuantities,
    Cosmology::AbstractCosmology)
    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                       #physical constants involved in calculations
    for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
        for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
            LensingFunction.WeightFunctionArray[idx_ZBinArray, idx_ZArray] =
            1.5 * (Cosmology.H0/c_0)^2 * Cosmology.ΩM *
            (1. + CosmologicalGrid.ZArray[idx_ZArray]) *
            BackgroundQuantities.rZArray[idx_ZArray]^2 *
            LensingFunction.LensingEfficiencyArray[idx_ZBinArray, idx_ZArray] +
            LensingFunction.IntrinsicAlignmentArray[idx_ZBinArray, idx_ZArray]
        end
    end
end