text:
## Iteration utilities

## some algorithms are easier to write using this neat
## package. However, for complicated usage, it proved
## to have issues with recursive usage and precompilation
using ResumableFunctions

# iterator over tuples of length n that sum to tot all non-negative
# used in non-commutative matching algorithm
@resumable function fixed_sum(n, tot, ps=())
    if n == 0
        @yield ()
        return
#    n > tot && return ()
    ## for all (k1, k2, ..., kseq) st sum(ks) = nfree
    elseif n == 1
        @yield (tot, ps...)
    else
        for i = 0:tot
            for v in collect(fixed_sum(n-1, tot-i, tuple(i, ps...)))
                @yield v
            end
        end
    end
end

# This comes from https://github.com/HPAC/matchpy/blob/master/matchpy/utils.py
# MIT licensed
# solve ax + by = c
# XXX: return (-1,-1) for nothing?
@resumable function base_solution_linear(a,b,c)

    @assert a > 0 && b > 0 && c >= 0

    NO_SOL = nothing

    d = gcd(a,b,c)
    a,b,c = div.((a,b,c), d)

    if c == 0
        @yield (0,0)
    else
        d, x0, y0 = gcdx(a, b)  # d = x0*a + y0*b

        !iszero(rem(c,d)) && return NO_SOL

        x, y = c*x0, c*y0

        if x <= 0
            while y >= 0
                if x >= 0
                    @yield (x,y)
                end
                x += b
                y -= a
            end
        else
            while x >=0
                if y >= 0
                    @yield (x,y)
                end
                x -= b
                y += a
            end
        end
    end
    NO_SOL
end

## solve x1*a + x2*b + ... = c for non-negative solutions
@resumable function solve_linear_diophantine(c, xs)
    NO_SOL = nothing

    n = length(xs)
    n == 0 && return NO_SOL
    if n == 1
        d,r = divrem(c, xs[1])
        r == 0 && @yield (d,)
        return NO_SOL
    end
    if n == 2
        u = base_solution_linear(xs[1], xs[2], c)
        for v in u
            @yield v
        end
        return NO_SOL
    end

    d = gcd(xs...)

    # solve xs[1] * x + d * y = c
    for val in base_solution_linear(xs[1], d, c)
        nxs = div.(xs[2:end], d)
        x, y = val
        for val in solve_linear_diophantine(y, nxs)
            @yield tuple(x, val...)
        end
    end

    NO_SOL
end

# trim iterator to yield only unique values
@resumable function just_unique(itr)
    theta = Any[]
    for i in itr
        if !(i in theta)
            @yield i
            push!(theta, i)
        end
    end
end

# Product generator
# for i in g1, j in g2, k in g3 ...
# do something with (i,j,k,...)
# end
# This is Iterators.product, but that was failing or some
# resumable functions
function gen_product(gens...)
    if length(gens) == 1
        [(i,) for i in collect(gens[1])] # not just collect(gens[1])
    elseif length(gens) == 2
        [tuple(g,h) for g in collect(gens[1]), h in gens[2]]
    else
        return [tuple(g, h...) for h in gen_product(gens[2:end]...), g in collect(gens[1])]
    end
end

## Create generator to cover case of
## for i in f1()
##    for j in f2(i)
##       for k in f3(j) ...
##
## Uses a backtracking algorithm
mutable struct GeneratorChain
fs
init
gs
sts
vals
end

function generator_chain(init, fs...)
    GeneratorChain(fs, init, Any[], Any[], Any[])
end

Base.iteratorsize(itertype::GeneratorChain) = Base.SizeUnknown()
function Base.start(o::GeneratorChain)

    g1 = o.fs[1](o.init...)
    st =start(g1)

    empty!(o.gs)
    empty!(o.sts)
    empty!(o.vals)
    push!(o.gs, g1)
    push!(o.sts, st)
    st
end

# advance to next tip
function Base.done(o::GeneratorChain, st)
    # is vals empty? If not, pop!
    !isempty(o.vals) && return false

    # okay, work back down
    while done(o.gs[end], o.sts[end])
        pop!(o.gs)
        pop!(o.sts)
        isempty(o.gs) && return true
    end

    # now back up
    i, n = length(o.gs), length(o.fs)

    val, o.sts[end] = next(o.gs[end], o.sts[end])
    for fn in o.fs[(i+1):n]
        gi = fn(val)
        st = start(gi)
        done(gi, st) && return done(o, st)
        val, st = next(gi, st)
        push!(o.gs, gi)
        push!(o.sts, st)
    end

    push!(o.vals, val)
    return false
end

function Base.next(o::GeneratorChain, st)
    pop!(o.vals), st
end


text:
include("noalloc_scalar.jl")
include("return.jl")

ext: jl


text:
# Solução de Sistemas de Equações Lineares: Métodos iterativos
# Gauss-Jacobi e Gauss-Seidel
using LinearAlgebra, Plots, Printf

# Função para atualização do vetor x
###################################################################################################
cn_update_x(A, b, x, i) = (b[i] - A[i, :] ⋅ x + A[i, i] * x[i]) / A[i, i]
cn_gaussJacobi(A, b, x, x0, i) = cn_update_x(A, b, x0, i)
cn_gaussSeidel(A, b, x, x0, i) = cn_update_x(A, b, x, i)

# Função de erro
###################################################################################################
function error(x, x_)
    nx, nΔx = norm(x), norm(x - x_)
    if nx < eps(); return nΔx else; return nΔx / nx; end
end

# Função para solução de sistema de equações lineares por métodos iterativos
###################################################################################################
function cn_slIter(A, b, x0=zeros(length(b)); met=cn_gaussSeidel, tol=1e-8, maxit=50)
    # Matriz A deve ser quadrada
    n = length(b)
    @assert size(A) == (n, n)
    # Inicialização
    err, nite, x = tol+1, 0, copy(x0)
    x_hist, err_hist = zeros(n, maxit+1), zeros(maxit)
    x_hist[:, 1] = x
    # Procedimento iterativo
    while err > tol && nite < maxit
        nite += 1
        # Atualização dos valores das variáveis, x
        for i = 1:n; x[i] = met(A, b, x, x_hist[:, nite], i); end
        # Atualização do erro
        err = error(x, x_hist[:, nite])
        # Atualização dos históricos
        x_hist[:, nite+1], err_hist[nite] = x, err
    end
    # Retorna históricos com valores das variáveis e erros por iteração
    return x_hist[:, 1:nite+1], err_hist[1:nite]
end

# Função com os critérios de convergência
###################################################################################################
function cn_conv(A)
    abs_A, n = abs.(A), size(A, 1)
    # Critério das linhas e colunas
    row = vec(sum(abs_A, dims=2)) ./ diag(abs_A) .- 1
    col = vec(sum(abs_A, dims=1)) ./ diag(abs_A) .- 1
    # Critério de Sassenfeld
    sas = zeros(n); sas[1] = row[1]
    for i = 2:n
        sas[i] = (sum(abs_A[i, i+1:end]) + sas[1:i-1] ⋅ abs_A[i, 1:i-1]) / abs_A[i, i]
    end
    return maximum(row), maximum(col), maximum(sas)
end

# Exemplos
###################################################################################################
exemplo = 3
if exemplo == 1
    A = [3 -0.1 -0.2; 0.1 7 -0.3; 0.3 -0.2 10];
    b = [7.85, -19.3, 71.4]
    x0 = [0., 0, 0]; tol = 1e-4; maxit = 50
elseif exemplo == 2
    A = [2. 1; 1 -2]; b = [2., -2]
    x0 = [0., 0]; tol = 1e-2; maxit = 50
elseif exemplo == 3
    A = [10. 2 1; 1 5 1; 2 3 10]; b = [14., 11, 8]
    x0 = [0., 0, 0]; tol = 1e-2; maxit = 50
elseif exemplo == 4
    A = [10. 5 3 1; 0 10 1 -2; 4 8 2 -1; 2 1 1 2];
    b = [9., -5, 3, 12]
    x0 = [0., 0, 0, 0]; tol = 1e-5; maxit = 200
end

# Cálculo da solução e apresentação dos resultados
###################################################################################################
fmt(output) = round.(output, sigdigits=4)
println("\nA: ", fmt(A)); println("b: ", fmt(b))

row, col, sas = cn_conv(A)
println("\nCritérios de convergência:")
@printf("linha: %.2f, coluna: %.2f, Sassenfeld: %.2f\n", row, col, sas)

println("\n### Gauss-Jacobi ###")
xGJ_hist, errGJ_hist = cn_slIter(A, b, x0, met=cn_gaussJacobi, tol=tol, maxit=maxit)
println("$(length(errGJ_hist)) iterações")
println("x_GJ: ", fmt(xGJ_hist))
println("err_GJ: ", fmt(errGJ_hist))

println("\n### Gauss-Seidel ###")
xGS_hist, errGS_hist = cn_slIter(A, b, x0, met=cn_gaussSeidel, tol=tol, maxit=maxit)
println("$(length(errGS_hist)) iterações")
println("x_GS: ", fmt(xGS_hist))
println("err_GS: ", fmt(errGS_hist))

gr(title="Histórico de erros", xlabel="Iteração", ylabel="Erro", yaxis=:log)
plot(errGJ_hist, label="GJ");
plot!(errGS_hist, label="GS")

max_stars_repo_name: ricardoaf/calcnum

max_issues_repo_name: ricardoaf/calcnum

max_forks_repo_name: ricardoaf/calcnum


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

# XXX: trailsign keyword not implemented, not sure it's worth the effort.
# Possible strategy: hold the sign of number in a variable (you can use
# copysign(1, number)), set result[1] equal to dd and multiply the first
# non-zero element (use findfirst for that) by the sign of number.
function sixty(number::T) where {T<:AbstractFloat}
    dd = trunc(abs(number))
    mm = trunc(abs(60 * number) - 60 * dd)
    ss = abs(3600 * number)
    return SVector{3,T}(trunc(number), mm, ss - 3600 * dd - 60 * mm)
end

"""
    sixty(number) -> [deg, min, sec]

### Purpose ###

Converts a decimal number to sexagesimal.

### Explanation ###

The reverse of `ten` function.

### Argument ###

* `number`: decimal number to be converted to sexagesimal.

### Output ###

An array of three `AbstractFloat`, that are the sexagesimal counterpart
(degrees, minutes, seconds) of `number`.

### Example ###

```jldoctest
julia> using AstroLib

julia> sixty(-0.615)
3-element StaticArrays.SArray{Tuple{3},Float64,1,3}:
 -0.0
 36.0
 54.0
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
sixty(n::Real) = sixty(float(n))


text:
include("dg/dg.jl")
include("dg_curved/dg.jl")
include("dg_common.jl")


text:
# SPDX-License-Identifier: MIT

push!(LOAD_PATH, "../src/")

using Documenter
using ExpFamilyDistributions

DocMeta.setdocmeta!(ExpFamilyDistributions, :DocTestSetup,
                    :(using ExpFamilyDistributions), recursive = true)

makedocs(
    sitename="ExpFamilyDistributions.jl",
    modules = [ExpFamilyDistributions],
    pages = [
        "Home" => "index.md",
        "Exponential Family Distributions" => "expfamily.md",
        "Distributions" => "dists.md",
    ],
)
deploydocs(
    repo = "github.com/lucasondel/ExpFamilyDistributions.jl.git",
)

max_stars_repo_name: lucasondel/ExpFamilyDistributions-jl

max_issues_repo_name: BUTSpeechFIT/ExpFamilyDistributions

max_forks_repo_name: lucasondel/ExpFamilyDistributions


text:
struct Pulse{T,F}
    Tp :: T
    f  :: F
end

(p::Pulse)(t) = p.f(t)

"""cosinepulse(T ; normalize = true)

A cosine pulse of duration `T`. This pulse is also known as a "time-domain
raised cosine" pulse."""

function cosinepulse(T ; normalize = true)
    p1(t) = begin
        if 0 ≤ t < T
            return 1.0-cospi(2*t/T)
        else
            return 0.0
        end
    end
    N = 1.0
    if normalize
        area, _ = quadgk(t -> p1(t)^2, 0.0, T)
        N = 1.0/sqrt(area)
    end
    return Pulse(T, t -> N*p1(t))
end

"""halfinepulse(T ; normalize = true)

A half-sine pulse of duration `T`."""
function  halfsinepulse(T ; normalize = true)
    p1(t) = begin
        if 0 ≤ t < T
            return sinpi(t/T)
        else
            return 0.0
        end
    end
    N = 1.0
    if normalize
        area, _ = quadgk(t -> p1(t)^2, 0.0, T)
        N = 1.0/sqrt(area)
    end
    return Pulse(T, t -> N*p1(t))
end

"""
    rcpulse(T, β, duration = 10T ; normalize = true)

A causal raised-cosine pulse with main lobe duration `T`, roll-off
factor (excess bandwidth) `β`, and duration `d`.
"""
function rcpulse(T, β, duration = 10T ; normalize = true)
    if β > 1 || β < 0
        @error "Excess bandwidth β must be between 0 and 1, inclusive." β
    end
    p1(t) = begin
        t ≥ duration/2 && return 0.0
        t < -duration/2 && return 0.0
        β == 0 && return sinc(t)
        if t ≈ T/(2β) || t ≈ -T/(2β)
            return (π/(4T))*sinc(1/(2β))
        end
        return (1/T)*sinc(t/T)*cospi(β*t/T)/(1-(2β*t/T)^2)
    end
    N = 1.0
    if normalize
        area, _ = quadgk(t -> p1(t)^2, -duration/2, duration/2)
        N = 1.0 / sqrt(area)
    end
    return Pulse(T, t -> N * p1(t - duration/2))
end

"""
    srrcpulse(T, β, duration = 10T ; normalize = true)

A causal square-root raised cosine pulse with main lobe duration `T`, roll-off
factor (excess bandwidth) `β``, and duration `d`.
"""
function srrcpulse(T, β, duration = 10T ; normalize = true)
    if β > 1 || β < 0
        @error "Excess bandwidth β must be between 0 and 1, inclusive." β
    end
    p1(t) = begin
        t ≥ duration/2 && return 0.0
        t < -duration/2 && return 0.0
        β == 0 && return sinc(t)
        t ≈ 0 && return (1 - β + 4β / π) / sqrt(T)
        if t ≈ T / (4β) || t ≈ -T / (4β)
            x = (β * (π + 2)) / (π * sqrt(2T)) * sinpi(1 / (4β))
            y = (β * (π - 2)) / (π * sqrt(2T)) * cospi(1 / (4β))
            return x + y
        else
            x = sinpi((1 - β) * t / T)
            y = (4β * t / T) * cospi((1 + β) * t / T)
            z = sqrt(T) * (π * t / T) * (1 - (4β * t / T)^2)
            return (x + y) / z
        end
    end
    N = 1.0
    if normalize
        area, _ = quadgk(t -> p1(t)^2, -duration/2, duration/2)
        N = 1.0 / sqrt(area)
    end
    return Pulse(T, t -> N * p1(t - duration/2))
end

"""gaussianpulse(B, T, duration=10 ; normalize = true)

A causal gaussian pulse with 3-dB bandwidth `B`, symbol time
`T`, and the specified `duration`.

See Molisch, "Wireless Communications", 11.2.2 and 11.3.9.
"""
function gaussianpulse(B, T, duration = 4T; normalize = true)
    p1(t) = begin
        t ≥ duration/2 && return 0.0
        t < -duration/2 && return 0.0
        σ = sqrt(log(2))/(2π*B*T)
        return 1/(sqrt(2π)*σ*T)*exp(-t^2/(2σ^2*T^2))
    end
    N = 1.0
    if normalize
        area, _ = quadgk(t -> p1(t)^2, -duration / 2, duration / 2)
        N = 1.0 / sqrt(area)
    end
    return Pulse(T, t -> N * p1(t - duration / 2))
end


text:
using Base.Test
include("../src/transducers.jl")
using Transducers

@test transduce(ttake(4), tpush!, Any[], 1:100) == Any[1, 2, 3, 4]

@test transduce(treplace(Dict{Any,Any}([(1, "A"),(4, "D")])), tpush!, Any[], 1:5) == Any["A", 2, 3, "D", 5]

@test transduce(dedupe, tpush!, Any[], [1, 1, 2, 2, 3, 3]) == Any[1, 2, 3]

@test transduce(partition_all(4), tpush!, Any[], 1:8) == Any[[1,2,3,4],[5,6,7,8]]
@test transduce(partition_all(4), tpush!, Any[], 1:5) == Any[[1,2,3,4],[5]]

# this is a dorky test, but failures should be rare
@test size(transduce(random_sample(0.99999), tpush!, Any[], 1:10))[1] == 10


text:
function byvolume{T,N}(f::TimeArray{T,N}, n::Int)

    x = f["Volume"].values
    increments = Int[]
    volumes = Int[]
    fm = 1; m = 1

    while  m < length(f)
        next_increment = findfirst(x->x > n,cumsum(x[fm:end]))
        push!(increments, next_increment)
        #push!(volumes, cumsum(x[fm:end]))
        push!(volumes, fm)
        fm += next_increment
        m  += 1
    end

    res = cumsum(increments[increments .!= 0])
    f[res] #, volumes[1:length(res)]
end


text:
module CyclicData

    using Dates
    using NCDatasets
    using CFTime
    using Formatting

    export CyclicDataManager, interpData!, makeDataContainer

    function toSec(dt)
        return Dates.Millisecond(dt).value / 1000.0
    end

    mutable struct CyclicDataManager

        timetype  :: DataType

        var_file_map:: Dict

        beg_time    :: AbstractCFDateTime
        end_time    :: AbstractCFDateTime
        align_time  :: AbstractCFDateTime
        period      :: Float64  # end_time - beg_time in seconds

        # The first and last index with in beg_time and end_time
        time_idx_beg   :: Integer
        time_idx_end   :: Integer

        t_vec_raw       :: AbstractArray{AbstractCFDateTime, 1}  # Contains the untrimmed time vector from input file
        t_vec           :: AbstractArray{Float64, 1}  # Actual time axis in seconds, has length
                                                      # (time_idx_end - time_idx_beg + 1).

        # Same as t_vec but appended 0 and period at
        # the beginning and end for interpolating usage.
        phantom_t_vec  :: AbstractArray{Float64, 1}

        # Precomputed arrays for time bound detection
        # - If t lies within [0, t_vec[1]] then interpolation
        #   will use data of idx_l_arr[1] and idx_r_arr[2].
        # - If t lies within [t_vec[i], t_vec[i+1]] then interpolation
        #   will use data of idx_l_arr[i+1] and idx_r_arr[i+1].
        # - If t lies within [t_vec[N], period] then interpolation
        #   will use data of idx_l_arr[N+1] and idx_r_arr[N+1]
        #   where N=length(t_vec)
        #
        # I use phantom_t_vec to simplify the algorithm above
        # - If t lies within [phantom_t_vec[i], phantom_t_vec[i+1]]
        #   then use data of idx_l_arr[i] and idx_r_arr[i] for interpolation.
        idx_l_arr      :: AbstractArray{Integer, 1}
        idx_r_arr      :: AbstractArray{Integer, 1}

        t_l_arr      :: AbstractArray{Float64, 1}
        t_r_arr      :: AbstractArray{Float64, 1}

        # Time pointer. Indicating the last data read. If it is zero then no previous data is read
        t_ptr_l   :: Integer
        t_ptr_r   :: Integer

        sub_yrng  :: Union{Colon, UnitRange}

        data_l   :: Union{Dict, Nothing}     # Data on the lower point. Used for interpolation
        data_r   :: Union{Dict, Nothing}     # Data on the upper point. Used for interpolation

        function CyclicDataManager(;
            timetype         :: DataType,
            var_file_map :: Dict,
            beg_time         :: AbstractCFDateTime,
            end_time         :: AbstractCFDateTime,
            align_time       :: AbstractCFDateTime,
            sub_yrng         :: Union{UnitRange, Colon} = Colon(),
            varname_time     :: String = "time",
        )

            data = Dict()
            t_vec = nothing

            if typeof(beg_time) != timetype
                throw(ErrorException("User specifies time type of {:s}, but beg_time has time type of {:s}", string(timetype), string(typeof(beg_time))))
            end

            if typeof(end_time) != timetype
                throw(ErrorException("User specifies time type of {:s}, but end_time has time type of {:s}", string(timetype), string(typeof(end_time))))
            end

            if typeof(align_time) != timetype
                throw(ErrorException("User specifies time type of {:s}, but align_time has time type of {:s}", string(timetype), string(typeof(align_time))))
            end

            period = toSec(end_time - beg_time)
            if period <= 0
                throw(ErrorException("End time cannot be earlier than beg time"))
            end

            local t_vec_raw

            local _compare_filename = nothing
            local _compare_t_vec = nothing
            local t_vec_raw

            for (varname, filename) in var_file_map

                Dataset(filename, "r") do ds
                    t_vec_raw = nomissing(ds[varname_time][:])
                    t_attrib = ds[varname_time].attrib

                    if typeof(t_vec_raw[1]) != timetype
                        throw(ErrorException("User specifies time type of $(string(timetype)), but data file has time type of $(string(typeof(t_vec_raw[1])))."))
                    end

                    if ! haskey(ds, varname)
                        throw(ErrorException("File $(filename) does not have variable $(varname)."))
                    end

                    if _compare_t_vec == nothing
                        _compare_filename = "$(filename)"
                        _compare_t_vec = copy(t_vec_raw)
                    end

                    if length(t_vec_raw) != length(_compare_t_vec)
                        throw(ErrorException("File $(filename) does not agree with other time dimension length in $(_compare_filename)."))
                    end

                    for i=1:length(t_vec_raw)
                        if t_vec_raw[i] != _compare_t_vec[i]
                            throw(ErrorException("`time` variable in file $(filename) is not the same as $(_compare_filename) at index $(i)"))
                        end
                    end
                end

            end

            if any(t_vec_raw[2:end] .<= t_vec_raw[1:end-1])
                throw(ErrorException("Time dimension has to be monotonically increasing."))
            end

            # Search for time_idx_beg and time_idx_end

            test = ( beg_time .<= t_vec_raw .<= end_time)
            time_idx_beg = findfirst(test)
            time_idx_end = findlast(test)

            #println("$(time_idx_beg) : $(time_idx_end)")

            if time_idx_end == nothing || time_idx_beg == nothing
                println(beg_time)
                println(end_time)
                println(t_vec_raw)
                throw(ErrorException("Time range is wrong that no time is within scope."))
            end

#            println("($(time_idx_beg), $time_idx_end)")

            t_vec = [
                mod(toSec(t_vec_raw[i] - beg_time), period) for i=time_idx_beg:time_idx_end
            ]

            if t_vec[1] == t_vec[end]
                throw(ErrorException("The first and last data of the obtained time series overlap. Please check."))
            end

#            for i in 1:length(t_vec)
#                println("$(t_vec_raw[i]) - $(align_time) = $(t_vec_raw[i] - align_time)")
#                println("t_vec[$(i)] = $(t_vec[i]/86400)")
#            end
#            println("Period = $(period / 86400)")

            #
            # - General idea of interpolation
            #
            # phantom
            #     1       2        3                N        N+1      N+2
            #     0 --- t[1] --- t[2] --- ... --- t[N-1] --- t[N] --- p(=0) --- t[1]   (time)
            #           d[1]     d[2]             d[N-1]     d[N]               d[1]   (data)
            # box    1        2       3       N-1        N        N+1
            # idx_l  N        1       2       N-2       N-1        N
            # idx_r  1        2       3       N-1        N         1
            #
            # t_l   t[N]-p   t[1]    t[2]                        t[N]
            # t_r   t[1]     t[2]    t[3]                        t[1]+p

            N = length(t_vec)

            phantom_t_vec = zeros(Float64, N+2)
            phantom_t_vec[1]       = 0.0
            phantom_t_vec[2:end-1] = t_vec
            phantom_t_vec[end]     = period

#            for (i, phantom_t) in enumerate(phantom_t_vec)
#                println("[$(i)] = $(phantom_t/86400)")
#            end

            # idx_r_arr[i] stores the index of right-point of box i
            idx_r_arr = collect(1:N+1)
            idx_r_arr[end] = 1

            # idx_l_arr[i] stores the index of left-point of box i
            idx_l_arr = circshift(idx_r_arr, (1, ))
            idx_l_arr[1] = N

            # Important: Need to offset the idx because we select the data
            # range according to beg_time and end_time
            idx_r_arr .+= time_idx_beg - 1
            idx_l_arr .+= time_idx_beg - 1

            # t_l_arr[i] stores the time used for left-point of box i
            t_l_arr = zeros(Float64, N+1)
            t_l_arr[1] = t_vec[end] - period
            t_l_arr[2:end] = t_vec

            # t_r_arr[i] stores the time used for right-point of box i
            t_r_arr = zeros(Float64, N+1)
            t_r_arr[1:end-1] = t_vec
            t_r_arr[end] = t_vec[1] + period

            obj = new(
                timetype,
                var_file_map,
                beg_time,
                end_time,
                align_time,
                period,
                time_idx_beg,
                time_idx_end,
                t_vec_raw,
                t_vec,
                phantom_t_vec,
                idx_l_arr,
                idx_r_arr,
                t_l_arr,
                t_r_arr,
                0,
                0,
                sub_yrng,
                nothing,
                nothing,
            )

            obj.data_l = makeDataContainer(obj)
            obj.data_r = makeDataContainer(obj)

            return obj
        end

    end

    function detectTimeBoundary(
        cdm :: CyclicDataManager,
        t   :: AbstractCFDateTime,
    )

        # lcr = left, center, right
        # Determine interpolation position
        local t_l, t_r, t_c, idx_l, idx_r

        t_c = mod( Second(t - cdm.align_time).value, cdm.period)
        for i = 1:length(cdm.phantom_t_vec)-1
#            println(format("t_c={}, left={}, right={} ", t_c,  cdm.phantom_t_vec[i], cdm.phantom_t_vec[i+1]))

            if cdm.phantom_t_vec[i] <= t_c <= cdm.phantom_t_vec[i+1]
                idx_l = cdm.idx_l_arr[i]
                idx_r = cdm.idx_r_arr[i]
                t_l = cdm.t_l_arr[i]
                t_r = cdm.t_r_arr[i]

#                println(format("Found $(i): t_c={}, left={}, right={} ", t_c,  t_l, t_r))
                break
            end
        end

        return idx_l, idx_r, t_l, t_c, t_r

    end

    function makeDataContainer(
        cdm      :: CyclicDataManager,
    )

        local data = Dict()

        for (varname, filename) in cdm.var_file_map
            Dataset(filename, "r") do ds

                s = size(ds[varname])

                Ny = s[2]

                if cdm.sub_yrng != Colon()
                    Ny = length(cdm.sub_yrng)
                end

                if length(s) == 4  # 3D case
                    data[varname] = zeros(Float64, s[3], s[1], Ny)
                elseif length(s) == 3  # 2D case
                    data[varname] = zeros(Float64, 1, s[1], Ny)
                else
                    throw(ErrorException("Unknown dimension: " * string(s)))
                end

            end
        end

        return data
    end

    function loadData!(
        cdm      :: CyclicDataManager,
        t_idx    :: Int64,
        data     :: Dict
    )

        missing_data = 0.0

        for (varname, filename) in cdm.var_file_map
            Dataset(filename, "r") do ds
                var = ds[varname]
                s = size(var)

                if length(s) == 4  # 3D case

                    data[varname][:, :, :] = permutedims(nomissing(ds[varname][:, cdm.sub_yrng, :, t_idx], missing_data), [3,1,2])
                elseif length(s) == 3  # 2D case
                    #println(varname, "; cdm.sub_yrng: ", cdm.sub_yrng)
                    #println("size of data container: ", size(data[varname]))
                    #println("size of ds : ", size(ds[varname][:, cdm.sub_yrng, t_idx]))
                    _tmp = nomissing(ds[varname][:, cdm.sub_yrng, t_idx], missing_data)
                    data[varname][:, :, :] = reshape(_tmp, 1, size(_tmp)...)
                else
                    throw(ErrorException("Unknown dimension: " * string(s)))
                end

            end
        end
    end

    function interpData!(
        cdm      :: CyclicDataManager,
        t        :: AbstractCFDateTime,
        data     :: Dict;
        create   :: Bool = false,
    )

        idx_l, idx_r, t_l, t_c, t_r = detectTimeBoundary(cdm, t)

        #println(format("{:d}, {:d}, {:.1f}", idx_l, idx_r, t_c))

        if cdm.t_ptr_l == 0 && cdm.t_ptr_r == 0 # Initialization

            # load idx_l into data_l
            # load idx_r into data_r

            cdm.t_ptr_l, cdm.t_ptr_r = idx_l, idx_r

            loadData!(cdm, cdm.t_ptr_l, cdm.data_l)
            loadData!(cdm, cdm.t_ptr_r, cdm.data_r)

        elseif cdm.t_ptr_r == idx_l

            # move data_r to data_l
            cdm.t_ptr_l, cdm.t_ptr_r = cdm.t_ptr_r, idx_r

            # swap data_l and data_r
            cdm.data_l, cdm.data_r = cdm.data_r, cdm.data_l

            # load idx_r into data_r
            loadData!(cdm, cdm.t_ptr_r, cdm.data_r)

        elseif cdm.t_ptr_l == idx_l && cdm.t_ptr_r == idx_r

            # do nothing.

        else
            throw(ErrorException(format("Unknown situation. cdm.t_ptr_l = {:d}, cdm.t_ptr_r = {:d}", cdm.t_ptr_l, cdm.t_ptr_r)))
        end

        Δt_lr = t_r - t_l
        Δt_lc = t_c - t_l
        Δt_cr = t_r - t_c

        coe_r = Δt_lc / Δt_lr
        coe_l = Δt_cr / Δt_lr

        if any([Δt_lr, Δt_lc, Δt_cr] .< 0)
            println("[Δt_lr, Δt_lc, Δt_cr] = ", [Δt_lr, Δt_lc, Δt_cr])
            throw(ErrorException("Δt sign error."))
        end

        for varname in keys(cdm.var_file_map)

            tmp = view(data[varname], :)

            # interpolation happens here
            data_l = view(cdm.data_l[varname], :)
            data_r = view(cdm.data_r[varname], :)
            @. tmp = data_l * coe_l + data_r * coe_r

        end
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
################################################################################
#                           TEST: Compose Rotations
################################################################################

import Base: isapprox

# Define the function `isapprox` for `EulerAngleAxis` to make the comparison
# easier.
function isapprox(x::EulerAngleAxis, y::EulerAngleAxis; keys...)
    a = isapprox(x.a, y.a; keys...)
    v = isapprox(x.v, y.v; keys...)

    a && v
end
#
# Define the function `isapprox` for `EulerAngles` to make the comparison
# easier.
function isapprox(x::EulerAngles, y::EulerAngles; keys...)
    a1 = isapprox(x.a1, y.a1; keys...)
    a2 = isapprox(x.a2, y.a2; keys...)
    a3 = isapprox(x.a3, y.a3; keys...)
    r  = x.rot_seq == y.rot_seq

    a1 && a2 && a3 && r
end

for i = 1:samples
    # DCMs
    # ==========================================================================

    # Sample 8 DCMs.
    for sym in [:D1, :D2, :D3, :D4, :D5, :D6, :D7, :D8]
        @eval $sym = angle_to_dcm(EulerAngles(-pi + 2*pi*rand(),
                                              -pi + 2*pi*rand(),
                                              -pi + 2*pi*rand(),
                                              rand($rot_seq_array)))
    end

    # Test the function `compose_rotation`.
    @test                   D2*D1 ≈ compose_rotation(D1,D2)                   atol=1e-14
    @test                D3*D2*D1 ≈ compose_rotation(D1,D2,D3)                atol=1e-14
    @test             D4*D3*D2*D1 ≈ compose_rotation(D1,D2,D3,D4)             atol=1e-14
    @test          D5*D4*D3*D2*D1 ≈ compose_rotation(D1,D2,D3,D4,D5)          atol=1e-14
    @test       D6*D5*D4*D3*D2*D1 ≈ compose_rotation(D1,D2,D3,D4,D5,D6)       atol=1e-14
    @test    D7*D6*D5*D4*D3*D2*D1 ≈ compose_rotation(D1,D2,D3,D4,D5,D6,D7)    atol=1e-14
    @test D8*D7*D6*D5*D4*D3*D2*D1 ≈ compose_rotation(D1,D2,D3,D4,D5,D6,D7,D8) atol=1e-14

    # Euler angle and axis
    # ==========================================================================

    # Sample 8 Euler angle and axis.
    for sym in [:ea1, :ea2, :ea3, :ea4, :ea5, :ea6, :ea7, :ea8 ]
        α   = -π + 2π*rand()
        v   = randn(3)
        v   = v/norm(v)

        @eval $sym = EulerAngleAxis($α,$v)
    end

    ear1 = ea2*ea1
    eac1 = compose_rotation(ea1,ea2)

    ear2 = ea3*ea2*ea1
    eac2 = compose_rotation(ea1,ea2,ea3)

    ear3 = ea4*ea3*ea2*ea1
    eac3 = compose_rotation(ea1,ea2,ea3,ea4)

    ear4 = ea5*ea4*ea3*ea2*ea1
    eac4 = compose_rotation(ea1,ea2,ea3,ea4,ea5)

    ear5 = ea6*ea5*ea4*ea3*ea2*ea1
    eac5 = compose_rotation(ea1,ea2,ea3,ea4,ea5,ea6)

    ear6 = ea7*ea6*ea5*ea4*ea3*ea2*ea1
    eac6 = compose_rotation(ea1,ea2,ea3,ea4,ea5,ea6,ea7)

    ear7 = ea8*ea7*ea6*ea5*ea4*ea3*ea2*ea1
    eac7 = compose_rotation(ea1,ea2,ea3,ea4,ea5,ea6,ea7,ea8)

    # Test the function `compose_rotation`.
    @test ear1 ≈ eac1 atol = 1e-14
    @test ear2 ≈ eac2 atol = 1e-14
    @test ear3 ≈ eac3 atol = 1e-14
    @test ear4 ≈ eac4 atol = 1e-14
    @test ear5 ≈ eac5 atol = 1e-14
    @test ear6 ≈ eac6 atol = 1e-14
    @test ear7 ≈ eac7 atol = 1e-14

    # Euler angles
    # ==========================================================================

    # Sample 8 Euler angles.
    for sym in [:Θ1, :Θ2, :Θ3, :Θ4, :Θ5, :Θ6, :Θ7, :Θ8 ]
        θx      = -π + 2π*rand()
        θy      = -π + 2π*rand()
        θz      = -π + 2π*rand()
        rot_seq = Meta.quot(rand([:XYZ,:ZYX]))

        @eval $sym = EulerAngles($θx, $θy, $θz, $rot_seq)
    end

    Θr1 = Θ2*Θ1
    Θc1 = compose_rotation(Θ1,Θ2)

    Θr2 = Θ3*Θ2*Θ1
    Θc2 = compose_rotation(Θ1,Θ2,Θ3)

    Θr3 = Θ4*Θ3*Θ2*Θ1
    Θc3 = compose_rotation(Θ1,Θ2,Θ3,Θ4)

    Θr4 = Θ5*Θ4*Θ3*Θ2*Θ1
    Θc4 = compose_rotation(Θ1,Θ2,Θ3,Θ4,Θ5)

    Θr5 = Θ6*Θ5*Θ4*Θ3*Θ2*Θ1
    Θc5 = compose_rotation(Θ1,Θ2,Θ3,Θ4,Θ5,Θ6)

    Θr6 = Θ7*Θ6*Θ5*Θ4*Θ3*Θ2*Θ1
    Θc6 = compose_rotation(Θ1,Θ2,Θ3,Θ4,Θ5,Θ6,Θ7)

    Θr7 = Θ8*Θ7*Θ6*Θ5*Θ4*Θ3*Θ2*Θ1
    Θc7 = compose_rotation(Θ1,Θ2,Θ3,Θ4,Θ5,Θ6,Θ7,Θ8)

    # Test the function `compose_rotation`.
    @test Θr1 ≈ Θc1 atol = 1e-14
    @test Θr2 ≈ Θc2 atol = 1e-14
    @test Θr3 ≈ Θc3 atol = 1e-14
    @test Θr4 ≈ Θc4 atol = 1e-14
    @test Θr5 ≈ Θc5 atol = 1e-14
    @test Θr6 ≈ Θc6 atol = 1e-14
    @test Θr7 ≈ Θc7 atol = 1e-14

    # Quaternions
    # ==========================================================================

    # Sample 8 quaternions.
    for sym in [:q1, :q2, :q3, :q4, :q5, :q6, :q7, :q8]
        @eval $sym = angle_to_quat(EulerAngles(-pi + 2*pi*rand(),
                                               -pi + 2*pi*rand(),
                                               -pi + 2*pi*rand(),
                                               rand($rot_seq_array)))
    end

    # Test the function `compose_rotation`.
    @test                   (q1*q2)[:] ≈ compose_rotation(q1,q2)[:]                   atol=1e-14
    @test                (q1*q2*q3)[:] ≈ compose_rotation(q1,q2,q3)[:]                atol=1e-14
    @test             (q1*q2*q3*q4)[:] ≈ compose_rotation(q1,q2,q3,q4)[:]             atol=1e-14
    @test          (q1*q2*q3*q4*q5)[:] ≈ compose_rotation(q1,q2,q3,q4,q5)[:]          atol=1e-14
    @test       (q1*q2*q3*q4*q5*q6)[:] ≈ compose_rotation(q1,q2,q3,q4,q5,q6)[:]       atol=1e-14
    @test    (q1*q2*q3*q4*q5*q6*q7)[:] ≈ compose_rotation(q1,q2,q3,q4,q5,q6,q7)[:]    atol=1e-14
    @test (q1*q2*q3*q4*q5*q6*q7*q8)[:] ≈ compose_rotation(q1,q2,q3,q4,q5,q6,q7,q8)[:] atol=1e-14
end


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# Breadth-first search / traversal

#################################################
#
#  Breadth-first visit
#
#################################################
"""
**Conventions in Breadth First Search and Depth First Search**
VertexColorMap :
- color == 0    => unseen
- color < 0     => examined but not closed
- color > 0     => examined and closed

EdgeColorMap :
- color == 0    => unseen
- color == 1     => examined
"""

type BreadthFirst <: SimpleGraphVisitAlgorithm
end

function breadth_first_visit_impl!(
    graph::SimpleGraph,                 # the graph
    queue::Vector{Int},                 # an (initialized) queue that stores the active vertices
    vertexcolormap::AbstractVertexMap,   # an (initialized) color-map to indicate status of vertices (-1=unseen, otherwise distance from root)
    edgecolormap::AbstractEdgeMap,        # an (initialized) color-map to indicate status of edges
    visitor::SimpleGraphVisitor,            # the visitor
    dir::Symbol)                        # direction [:in,:out]

    fneig = dir == :out ? out_neighbors : in_neighbors
    while !isempty(queue)
        u = shift!(queue)
        open_vertex!(visitor, u)
        u_color = vertexcolormap[u]

        for v in fneig(graph, u)
            v_color = get(vertexcolormap, v, 0)
            v_edge = Edge(u,v)
            e_color = get(edgecolormap, v_edge, 0)
            examine_neighbor!(visitor, u, v, u_color, v_color, e_color) || return
            edgecolormap[v_edge] = 1
            if v_color == 0
                vertexcolormap[v] = u_color - 1
                discover_vertex!(visitor, v) || return
                push!(queue, v)
            end
        end
        close_vertex!(visitor, u)
        vertexcolormap[u] *= -1
    end
end

function traverse_graph!(
    graph::SimpleGraph,
    alg::BreadthFirst,
    source,
    visitor::SimpleGraphVisitor;
    vertexcolormap::AbstractVertexMap = Dict{Int, Int}(),
    edgecolormap::AbstractEdgeMap = DummyEdgeMap(),
    queue = Vector{Int}(),
    dir = :out)

    for s in source
        vertexcolormap[s] = -1
        discover_vertex!(visitor, s) || return
        push!(queue, s)
    end

    breadth_first_visit_impl!(graph, queue, vertexcolormap, edgecolormap
            , visitor, dir)
end

#################################################
#
#  Useful applications
#
#################################################

###########################################
# Get the map of the (geodesic) distances from vertices to source by BFS                  #
###########################################

immutable GDistanceVisitor <: SimpleGraphVisitor
end

"""
    gdistances!(g, source, dists) -> dists

Fills `dists` with the geodesic distances of vertices in  `g` from vertex/vertices `source`.
`dists` can be either a vector or a dictionary.
"""
function gdistances!(g::SimpleGraph, source, dists)
    visitor = GDistanceVisitor()
    traverse_graph!(g, BreadthFirst(), source, visitor, vertexcolormap=dists)
    for i in eachindex(dists)
        dists[i] -= 1
    end
    return dists
end

"""
    gdistances(g, source) -> dists

Returns a vector filled with the geodesic distances of vertices in  `g` from vertex/vertices `source`.
For vertices in disconnected components the default distance is -1.
"""
gdistances(g::SimpleGraph, source) = gdistances!(g, source, fill(0,nv(g)))

###########################################
# Constructing BFS trees                  #
###########################################

"""TreeBFSVisitorVector is a type for representing a BFS traversal
of the graph as a parents array. This type allows for a more performant implementation.
"""
type TreeBFSVisitorVector <: SimpleGraphVisitor
    tree::Vector{Int}
end

function TreeBFSVisitorVector(n::Int)
    return TreeBFSVisitorVector(fill(0, n))
end

"""tree converts a parents array into a DiGraph"""
function tree(parents::AbstractVector)
    n = length(parents)
    t = DiGraph(n)
    for i in 1:n
        parent = parents[i]
        if parent > 0  && parent != i
            add_edge!(t, parent, i)
        end
    end
    return t
end

tree(parents::TreeBFSVisitorVector) = tree(parents.tree)

function examine_neighbor!(visitor::TreeBFSVisitorVector, u::Int, v::Int,
                            ucolor::Int, vcolor::Int, ecolor::Int)
    if u != v && vcolor == 0
        visitor.tree[v] = u
    end
    return true
end

function bfs_tree!(visitor::TreeBFSVisitorVector,
        g::SimpleGraph,
        s::Int;
        vertexcolormap = Dict{Int,Int}(),
        queue = Vector{Int}())
    # this version of bfs_tree! allows one to reuse the memory necessary to compute the tree
    # the output is stored in the visitor.tree array whose entries are the vertex id of the
    # parent of the index. This function checks if the scratch space is too small for the graph.
    # and throws an error if it is too small.
    # the source is represented in the output by a fixed point v[root] == root.
    # this function is considered a performant version of bfs_tree for useful when the parent
    # array is more helpful than a DiGraph struct, or when performance is critical.
    nvg = nv(g)
    length(visitor.tree) >= nvg || error("visitor.tree too small for graph")
    visitor.tree[s] = s
    traverse_graph!(g, BreadthFirst(), s, visitor; vertexcolormap=vertexcolormap, queue=queue)
end

"""Provides a breadth-first traversal of the graph `g` starting with source vertex `s`,
and returns a directed acyclic graph of vertices in the order they were discovered.

This function is a high level wrapper around bfs_tree!, use that function for more performance.
"""
function bfs_tree(g::SimpleGraph, s::Int)
    nvg = nv(g)
    visitor = TreeBFSVisitorVector(nvg)
    bfs_tree!(visitor, g, s)
    return tree(visitor)
end

############################################
# Connected Components with BFS            #
############################################
"""Performing connected components with BFS starting from seed"""
type ComponentVisitorVector <: SimpleGraphVisitor
    labels::Vector{Int}
    seed::Int
end

function examine_neighbor!(visitor::ComponentVisitorVector, u::Int, v::Int,
                            ucolor::Int, vcolor::Int, ecolor::Int)
    if u != v && vcolor == 0
        visitor.labels[v] = visitor.seed
    end
    return true
end

############################################
# Test graph for bipartiteness             #
############################################
type BipartiteVisitor <: SimpleGraphVisitor
    bipartitemap::Vector{UInt8}
    is_bipartite::Bool
end

BipartiteVisitor(n::Int) = BipartiteVisitor(zeros(UInt8,n), true)

function examine_neighbor!(visitor::BipartiteVisitor, u::Int, v::Int,
        ucolor::Int, vcolor::Int, ecolor::Int)
    if vcolor == 0
        visitor.bipartitemap[v] = (visitor.bipartitemap[u] == 1) ? 2 : 1
    else
        if visitor.bipartitemap[v] == visitor.bipartitemap[u]
            visitor.is_bipartite = false
        end
    end
    return visitor.is_bipartite
end

"""
    is_bipartite(g)
    is_bipartite(g, v)

Will return `true` if graph `g` is [bipartite](https://en.wikipedia.org/wiki/Bipartite_graph).
If a node `v` is specified, only the connected component to which it belongs is considered.
"""
function is_bipartite(g::SimpleGraph)
    cc = filter(x->length(x)>2, connected_components(g))
    vmap = Dict{Int,Int}()
    for c in cc
        _is_bipartite(g,c[1], vmap=vmap) || return false
    end
    return true
end

is_bipartite(g::SimpleGraph, v::Int) = _is_bipartite(g, v)

_is_bipartite(g::SimpleGraph, v::Int; vmap = Dict{Int,Int}()) = _bipartite_visitor(g, v, vmap=vmap).is_bipartite

function _bipartite_visitor(g::SimpleGraph, s::Int; vmap=Dict{Int,Int}())
    nvg = nv(g)
    visitor = BipartiteVisitor(nvg)
    for v in keys(vmap) #have to reset vmap, otherway problems with digraphs
        vmap[v] = 0
    end
    traverse_graph!(g, BreadthFirst(), s, visitor, vertexcolormap=vmap)
    return visitor
end

"""
If the graph is bipartite returns a vector `c`  of size `nv(g)` containing
the assignment of each vertex to one of the two sets (`c[i] == 1` or `c[i]==2`).
If `g` is not bipartite returns an empty vector.
"""
function bipartite_map(g::SimpleGraph)
    cc = connected_components(g)
    visitors = [_bipartite_visitor(g, x[1]) for x in cc]
    !all([v.is_bipartite for v in visitors]) && return zeros(Int, 0)
    m = zeros(Int, nv(g))
    for i=1:nv(g)
        m[i] = any(v->v.bipartitemap[i] == 1, visitors) ? 2 : 1
    end
    m
end


text:
#!/usr/bin/env julia

using Plots; plotlyjs()

h1=histogram(randn(1000))

gui(h1)

readline()

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
export DeMatosCutPruner, DeMatosPruningAlgo

"""
    DeMatosPruningAlgo <: AbstractCutPruningAlgo

Removes the cuts with lower trust where the trust is the number
of points `x` associated to the cuts.
The more points are associated, the higher is the trust.

We refer to [1] for further details.

[1] De Matos, Vitor L., Andy B. Philpott, and Erlon C. Finardi.
"Improving the performance of stochastic dual dynamic programming."
Journal of Computational and Applied Mathematics 290 (2015): 196-208.
"""
struct DeMatosPruningAlgo <: AbstractCutPruningAlgo
    # maximum number of cuts
    maxncuts::Int
    function DeMatosPruningAlgo(maxncuts::Int)
        new(maxncuts)
    end
end

mutable struct DeMatosCutPruner{N, T} <: AbstractCutPruner{N, T}
    # used to generate cuts
    isfun::Bool
    islb::Bool
    lazy_minus::Bool
    A::AbstractMatrix{T}
    b::AbstractVector{T}

    maxncuts::Int

    trust::Vector{Float64}
    ids::Vector{Int} # small id means old
    id::Int # current id

    #set of states where cut k is active
    territories::Vector{Vector{Tuple{Int, T}}}
    nstates::Int
    states::Matrix{T}

    # peform exhaustive check before adding cuts
    excheck::Bool

    # tolerance to check redundancy between two cuts
    TOL_EPS::Float64

    function (::Type{DeMatosCutPruner{N, T}})(sense::Symbol, maxncuts::Int, lazy_minus::Bool=false, tol=1e-6, excheck::Bool=false) where {N, T}
        isfun, islb = gettype(sense)
        new{N, T}(isfun, islb, lazy_minus, spzeros(T, 0, N), T[], maxncuts, Tuple{Int, T}[], Int[], 0, [], 0, zeros(T, 0, N), excheck, tol)
    end
end

(::Type{CutPruner{N, T}})(algo::DeMatosPruningAlgo, sense::Symbol, lazy_minus::Bool=false) where {N, T} = DeMatosCutPruner{N, T}(sense, algo.maxncuts, lazy_minus)

getnreplaced(man::DeMatosCutPruner, R, ncur, nnew, mycut) = nnew, length(R)

hasterritories(man::DeMatosCutPruner) = true

"""
    addposition!(man::DeMatosCutPruner, position::Matrix)

Update territories with cuts previously computed during backward pass.

# Arguments
* `man::DeMatosCutPruner`
    Pruner to update.
* `position::Array{T, 2}`
    New visited positions, corresponding to a collection of points.
"""
function addposition!(man::DeMatosCutPruner, position::Matrix)
    # get number of new positions to analyse:
    nx = size(position, 1)

    for i in 1:nx
        addstate!(man, position[i, :])
    end

    updatetrust!(man)
end

function addposition!(man::DeMatosCutPruner, position::Vector)
    addstate!(man, position)
    updatetrust!(man)
end

"""
    addstate!(man::DeMatosCutPruner, x::Vector)

Add a new state to test and accordingly update territories of each cut.
"""
function addstate!(man::DeMatosCutPruner, x::Vector)
    # update number of states
    man.nstates += 1
    # Add `x` to the list of visited state:
    man.states = vcat(man.states, x')

    giveterritory!(man, man.nstates, x)
end

function giveterritory!(man::DeMatosCutPruner, ix::Int, x::Vector=man.states[ix, :])
    if ncuts(man) > 0
        # Get cut which is active at point `x`:
        bcost, bcuts = optimalcut(man, x)
        # Add `x` with index nstates  to the territory of cut with index `bcuts`:
        push!(man.territories[bcuts], (ix, bcost))
    end
end

"""
    optimalcut{T}(man::DeMatosCutPruner, xf::Vector{T})

Find active cut at point `xf`.

# Arguments
* `man::DeMatosCutPruner`:
    CutPruner
* `xf::Vector{T}`:

# Return
* `bestcost::T`
    Value of supporting cut at point `xf`
* `bestcut::Int`
    Index of supporting cut at point `xf`
"""
function optimalcut(man::DeMatosCutPruner,
                       xf::Vector{T}) where {T}
    bestcost = -Inf::Float64
    bestcut = -1
    dimstates = length(xf)
    nc = ncuts(man)

    @inbounds for i in 1:nc
        cost = cutvalue(man, i, xf)
        if cost > bestcost
            bestcost = cost
            bestcut = i
        end
    end
    return bestcost, bestcut
end

"""
    updateterritory!(man::DeMatosCutPruner, indcut::Int)

Update territories (i.e. the set of tested states where
    a given cut is active) considering new cut given by index `indcut`.

# Arguments
* `man::DeMatosCutPruner`:
* `indcut::Int`:
    new cut index
"""
function updateterritory!(man::DeMatosCutPruner, indcut::Int)
    @assert length(man.territories) == ncuts(man)
    for k in 1:ncuts(man)
        if k == indcut
            continue
        end
        todelete = []
        for (num, (ix, cost)) in enumerate(man.territories[k])
            x = man.states[ix, :]

            costnewcut = cutvalue(man, indcut, x)

            if costnewcut > cost
                push!(todelete, num)
                push!(man.territories[indcut], (ix, costnewcut))
            end
        end
        deleteat!(man.territories[k], todelete)
    end
end

"""
    cutvalue{T}(man::DeMatosCutPruner, indc::Int, x::Vector{T})

Get value of cut with index `indc` at point `x`.

# Arguments
- `man::DeMatosCutPruner`
    Approximation of the value function as linear cuts
- `indc::Int`
    Index of cut to consider
- `x::Vector{T}`
    Coordinates of state

# Return
`cost::T`
    Value of cut `indc` at point `x`.
    If `man` is a polyhedral function, then it is the value of the cut at `x`,
    otherwise, it is the distance between `x` and the cut.
    As a rule of thumb, the higher the `cutvalue` is, the less it is redundant.
"""
function cutvalue(man::DeMatosCutPruner, indc::Int, x::Vector{T}) where {T}
    β = man.b[indc]
    a = @view man.A[indc, :]
    ax = dot(a, x)
    if man.lazy_minus
        ax = -ax
    end
    cost = isfun(man) ? ax + β : (β - ax) / norm(a, 2)
    islb(man) ? cost : -cost
end

flength(a)::Float64 = length(a)

function updatetrust!(man)
    @assert length(man.territories) == ncuts(man)
    if ncuts(man) == length(man.trust)
        # Avoid new allocation. Avoiding this allocation is the whole point of
        # doint updatetrust! instead of doing trust = nothing
        for i in 1:ncuts(man)
            man.trust[i] = length(man.territories[i])
        end
    else
        man.trust = flength.(man.territories)
    end
end

function replacecuts!(man::DeMatosCutPruner{N, T}, K::AbstractVector{Int}, A, b, mycut::AbstractVector{Bool}) where {N, T}
    @assert length(man.territories) == ncuts(man)
    # FIXME If K is 1:ncuts, then checkconsistency will be true and trust will not be recomputed by gettrust
    _replacecuts!(man, K, A, b)
    # Do not do view here since will will modify the entries
    freeterritories = man.territories[K]
    for k in K
        man.territories[k] = Tuple{Int, T}[]
    end
    for k in K
        updateterritory!(man, k)
    end
    for freet in freeterritories
        for (ik, _) in freet
            giveterritory!(man, ik)
        end
    end
    updatetrust!(man)
    @assert length(man.territories) == ncuts(man)
end

"""Push new cut in CutPruner `man`."""
function appendcuts!(man::DeMatosCutPruner{N, T}, A, b, mycut::AbstractVector{Bool}) where {N, T}
    @assert length(man.territories) == ncuts(man)
    oldncuts = ncuts(man)
    _appendcuts!(man, A, b)
    nnew = length(b)
    newterritories = [Tuple{Int, T}[] for _ in 1:nnew]
    if isempty(man.territories)
        man.territories = newterritories
        # No cut yet, we share territories among new cuts
        for i in 1:size(man.states, 1)
            giveterritory!(man, i)
        end
    else
        man.territories = vcat(man.territories, newterritories)
        for k in oldncuts .+ (1:nnew)
            updateterritory!(man, k)
        end
    end
    updatetrust!(man)
    @assert length(man.territories) == ncuts(man)
end


text:
using SAIMFitKit

#-------------- OPTICAL PARAMETERS----------------------

opt = SAIMOptics()
opt.nB = 1.33;	                    #The refractive index of the ambient buffer / cytoplasm
opt.dOx = 1924.9;	                #The thickness of the SiO2 layer in units of nm
opt.λ_Ex_1 = 642.0;		        	#The wavelength of the microscope excitation laser #1 in units of nm
opt.nOx_1 = 1.4719;		            #The refractive index of SiO2
opt.nSi_1 = 3.8660 + 0.017933im;    #The complex refractive index of Si

#-------------- FIT PARAMETERS----------------------
p0 = [0.5, 1.0, 75.0]		    #Initial guesses for parameters [A, B, H]
lb = [0.0, 0.0, 0.0]			#Lower bounds for [A, B, H]
ub = [10000., 10000., 200.]		#Upper bounds for [A, B, H]

#-------------- IMAGE STACKS ---------------------
angles = range(5.0, length=32, stop=43.75)	    #Incidece angle in degrees for each image frame
files = ["21TR_5c"]                             #Names of image stacks; multiple experiments can be listed to run in batch
path = "C:\\Users\\matth\\.julia\\dev\\SAIMFitKit\\example\\TestImages" #Path to image files

#-------------- CONDUCT FITS ---------------------
for f in files
    fit_SAIM(path, f, opt, angles, p0, lb, ub, glb=true, step=40.0, show=true)
end


text:
module Spaces

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, template, register_normal_element

export space

register_normal_element("q__space", context = @__MODULE__)

function space(args...; kwargs...)
  q__space(args...; kw(kwargs)...)
end

end


text:
function createnodes(p::Array{Float64,2},t::Array{Int,2},porder::Int)
# %CREATEDGNODES Computes the Coordinates of the DG nodes.
# %   DGNODES=CREATENODES(MESH,FD,FPARAMS)
# %
# %      MESH:      Mesh Data Structure
# %      FD:        Distance Function d(x,y)
# %      FPARAMS:   Additional parameters passed to FD
# %      DGNODES:   Triangle indices (NPL,2,NT). The nodes on
# %                 the curved boundaries are projected to the
# %                 true boundary using the distance function FD
# %

# % npv : number of nodes per volume element
# % nfv : number of faces per volume element
# % npf : number of nodes per face element

nve,ne = size(t);
nd=size(p,1);

elemtype = 0;
if (nd==2) && (nve==4)
    elemtype=1;
end
if (nd==3) && (nve==8)
    elemtype=1;
end

plocal,~,~,~,perm = masternodes(porder,nd,elemtype);
npl = size(plocal,1);

philocal = zeros(Float64, npl, nve);

if nd==1
    xi  = plocal[:,1];
    philocal[:,1] = 1.0 .- xi;
    philocal[:,2] = xi;
elseif nd==2 && nve==3 # tri
    xi  = plocal[:,1];
    eta = plocal[:,2];
    philocal[:,1] = 1.0 .- xi .- eta;
    philocal[:,2] = xi;
    philocal[:,3] = eta;
elseif nd==2 && nve==4 # quad
    xi  = plocal[:,1];
    eta = plocal[:,2];
    philocal[:,1] = (1.0 .- xi).*(1.0 .- eta);
    philocal[:,2] = xi.*(1.0 .- eta);
    philocal[:,3] = xi.*eta;
    philocal[:,4] = (1.0 .- xi).*eta;
elseif nd==3 && npv==4 # tet
    xi   = plocal[:,1];
    eta  = plocal[:,2];
    zeta = plocal[:,3];
    philocal[:,1] = 1.0 - xi - eta - zeta;
    philocal[:,2] = xi;
    philocal[:,3] = eta;
    philocal[:,4] = zeta;
elseif nd==3 && npv==8 # hex
    xi   = plocal[:,1];
    eta  = plocal[:,2];
    zeta = plocal[:,3];
    philocal[:,1] = (1.0 .-xi).*(1.0 .-eta).*(1.0 .-zeta);
    philocal[:,2] = xi.*(1.0 .-eta).*(1.0 .-zeta);
    philocal[:,3] = xi.*eta.*(1.0 .-zeta);
    philocal[:,4] = (1.0 .-xi).*eta.*(1.0 .-zeta);
    philocal[:,5] = (1.0 .-xi).*(1.0 .-eta).*(zeta);
    philocal[:,6] = xi.*(1.0 .-eta).*(zeta);
    philocal[:,7] = xi.*eta.*(zeta);
    philocal[:,8] = (1.0 .-xi).*eta.*(zeta);
end

# Allocate nodes
dgnodes = zeros(Float64,npl,ne,nd);
for dim=1:nd
  for node=1:nve
    dp = reshape(philocal[:,node],npl,1)*reshape(p[dim,t[node,:]],1,ne);
    dgnodes[:,:,dim] = dgnodes[:,:,dim] + reshape(dp,npl,ne);
  end
end

dgnodes = permutedims(dgnodes,[1, 3, 2]);

return dgnodes, elemtype, perm

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
mutable struct ComplementarityType
    lb::AbstractFloat
    var::JuMP.VariableRef
    ub::AbstractFloat
    F::JuMP.NonlinearExpression
    raw_idx::Integer
    var_name::String
    F_name::String
    result_value::AbstractFloat
end

function MCPModel()
    m = JuMP.Model()
    m.ext[:MCP] = ComplementarityType[]
    return m
end

function get_MCP_data(m::JuMP.Model)
    if haskey(m.ext, :MCP)
        return m.ext[:MCP]::Array
    else
        error("The 'get_MCP_data' function is only for MCP models as in ComplementarityType.jl")
    end
end

raw_index(v::JuMP.VariableRef) = JuMP.index(v).value

function get_bounds_in_raw_index(mcp_data)
    n = length(mcp_data)
    lb = zeros(n)
    ub = ones(n)
    for i in 1:n
        lb[raw_index(mcp_data[i].var)] = mcp_data[i].lb
        ub[raw_index(mcp_data[i].var)] = mcp_data[i].ub
    end
    return lb, ub
end

function get_names_in_raw_index(mcp_data)
    n = length(mcp_data)
    var_name = Array{String}(undef, n)
    F_name = Array{String}(undef, n)
    for i in 1:n
        var_name[raw_index(mcp_data[i].var)] = mcp_data[i].var_name
        F_name[raw_index(mcp_data[i].var)] = mcp_data[i].F_name
    end
    return var_name, F_name
end

function get_initial_values_in_raw_index(m, mcp_data)
    n = length(mcp_data)
    initial_values = Array{Float64}(undef, n)
    for i in 1:n
        init_val = JuMP.start_value(mcp_data[i].var)
        if init_val == nothing
            init_val = mcp_data[i].lb
        end
        initial_values[raw_index(mcp_data[i].var)] = init_val
    end
    return initial_values
end

function solveMCP(m::JuMP.Model; solver=:PATH, method=:trust_region, linear=false)
    if solver == :PATH
        return _solve_path(m, linear=linear)
    elseif solver == :NLsolve
        return _solve_nlsolve(m, method=method)
    end
end

function solveLCP(m::JuMP.Model; solver=:PATH, method=:trust_region)
    solveMCP(m, solver=solver, method=method, linear=true)
end

function sortperm_MCP_data(obj::Array{ComplementarityType,1})
    n = length(obj)
    ref = Array{Int}(undef, n)
    for i in 1:n
        ref[i] = obj[i].raw_idx
    end

    return sortperm(ref)
end

# Using PATHSolver
function _solve_path(m::JuMP.Model; linear=false)

    function myfunc(z)
        # z is in RawIndex, passed from PATHSolver
        d = JuMP.NLPEvaluator(m)
        MOI.initialize(d, [:Grad])
        F_val = zeros(n)
        MOI.eval_constraint(d, F_val, z)

        # F_val also should be in RawIndex
        # since it is the order in which constraints are added

        return F_val
    end

    function myjac(z)
        # z is in RawIndex, passed from PATHSolver
        d = JuMP.NLPEvaluator(m)
        MOI.initialize(d, [:Grad])
        J_struct = MOI.jacobian_structure(d)

        I = first.(J_struct)
        J = last.(J_struct)

        jac_val = zeros(length(J))
        MOI.eval_constraint_jacobian(d, jac_val, z)

        # return matrix also should be in RawIndex
        # since it is the order in which constraints are added

        return sparse(I, J, jac_val)
    end

    mcp_data = get_MCP_data(m)
    n = length(mcp_data)

    # Two Indices
    # MCP_Index: the order stored in MCPModel = array index of Array{ComplementarityType}
    # RawIndex: the order used in JuMP / MathOptInteface

    # Declaring MCP mapping F as constraints
    # in order to query Jacobian using AutoDiff thru MathProgBase
    # i = RawIndex
    # Add constraint in the order of RawIndex
    p = sortperm_MCP_data(mcp_data)
    JuMP.@NLconstraint(m, [i=1:n], mcp_data[p[i]].F == 0)

    # lb and ub in RawIndex
    lb, ub = get_bounds_in_raw_index(mcp_data)

    # var_name, F_name in RawIndex
    var_name, F_name = get_names_in_raw_index(mcp_data)

    # initial values
    initial_values = get_initial_values_in_raw_index(m, mcp_data)

    # Solve the MCP using PATHSolver
    # ALL inputs to PATHSolver must be in RawIndex
    if linear==true
        J0 = myjac(zeros(size(lb)))
        Jr = myjac(100*rand(Float64, size(ub)))

        if !isapprox(J0, Jr)
            error("The mappings do not seem linear. Rerun 'solveMCP()' after removing 'linear=true'.")
        end
        status, z, f = PATHSolver.solveLCP(myfunc, J0, lb, ub, initial_values, var_name, F_name)
    else
        status, z, f = PATHSolver.solveMCP(myfunc, myjac, lb, ub, initial_values, var_name, F_name)
    end
    # z, f are in RawIndex

    # After solving set the values in m::JuMP.Model to the solution obtained.
    for i in 1:n
        # setvalue(mcp_data[i].var, z[mcp_data[i].raw_idx])
        # MOI.set(m, MOI.VariablePrimal(), mcp_data[i].var, z[mcp_data[i].raw_idx])
        # mcp_data[i].result_value = z[mcp_data[i].raw_idx]
        set_result_value(mcp_data[i], z[mcp_data[i].raw_idx])
    end

    # Cleanup. Remove all dummy @NLconstraints added,
    # so that the model can be re-used for multiple runs
    # Array{JuMP.NonlinearConstraint,1}(undef, 0)
    # m.nlp_data.nlconstr = JuMP.NonlinearConstraint[]
    m.nlp_data.nlconstr = []

    # This function has changed the content of m already.
    return status
end

function _solve_nlsolve(m::JuMP.Model; method=:trust_region)

    function myfunc!(fvec, z)
        # z is in RawIndex, passed from PATHSolver
        d = JuMP.NLPEvaluator(m)
        MOI.initialize(d, [:Grad])
        F_val = zeros(n)
        MOI.eval_constraint(d, F_val, z)

        copyto!(fvec, F_val)
    end

    function myjac!(fjac, z)
        # z is in RawIndex, passed from PATHSolver
        d = JuMP.NLPEvaluator(m)
        MOI.initialize(d, [:Grad])
        J_struct = MOI.jacobian_structure(d)

        I = first.(J_struct)
        J = last.(J_struct)

        jac_val = zeros(length(J))
        MOI.eval_constraint_jacobian(d, jac_val, z)

        sparse_fjac = sparse(I, J, jac_val)
        copyto!(fjac, full(sparse_fjac))
    end

    mcp_data = get_MCP_data(m)
    n = length(mcp_data)

    # Two Indices
    # MCP_Index: the order stored in MCPModel = array index of Array{ComplementarityType}
    # RawIndex: the order used in JuMP / MathProgBase

    # Declaring MCP mapping F as constraints
    # in order to query Jacobian using AutoDiff thru MathProgBase
    # i = RawIndex
    # Add constraint in the order of RawIndex
    p = sortperm_MCP_data(mcp_data)
    JuMP.@NLconstraint(m, [i=1:n], mcp_data[p[i]].F == 0)

    # lb and ub in RawIndex
    lb, ub = get_bounds_in_raw_index(mcp_data)

    # initial values
    initial_values = get_initial_values_in_raw_index(m, mcp_data)

    # Solve the MCP using NLsolve
    # ALL inputs to NLsolve must be in RawIndex

    r = NLsolve.mcpsolve(myfunc!, myjac!, lb, ub, initial_values, method = method,
        iterations = 10_000)
    # function mcpsolve{T}(f,
    #                   j,
    #                   lower::Vector,
    #                   upper::Vector,
    #                   initial_x::AbstractArray{T};
    #                   method::Symbol = :trust_region,
    #                   reformulation::Symbol = :smooth,
    #                   xtol::Real = zero(T),
    #                   ftol::Real = convert(T,1e-8),
    #                   iterations::Integer = 1_000,
    #                   store_trace::Bool = false,
    #                   show_trace::Bool = false,
    #                   extended_trace::Bool = false,
    #                   linesearch = LineSearches.BackTracking(),
    #                   factor::Real = one(T),
    #                   autoscale = true,
    #                   inplace = true)

    # julia> fieldnames(typeof(r))
    # 12-element Array{Symbol,1}:
    #  :method
    #  :initial_x
    #  :zero
    #  :residual_norm
    #  :iterations
    #  :x_converged
    #  :xtol
    #  :f_converged
    #  :ftol
    #  :trace
    #  :f_calls
    #  :g_calls

    # r.zero is in RawIndex
    for i in 1:n
        set_result_value(mcp_data[i], r.zero[mcp_data[i].raw_idx])
    end

    # Cleanup. Remove all dummy @NLconstraints added,
    # so that the model can be re-used for multiple runs
    # m.nlp_data.nlconstr =  Array{JuMP.NonlinearConstraint,1}(undef, 0)
    m.nlp_data.nlconstr =  []

    # This function has changed the content of m already.
    return r
end

function set_result_value(mcp_data::ComplementarityType, value::Float64)
    mcp_data.result_value = value
end

function set_result_value(v::JuMP.VariableRef, value::Float64)
    mcp_data = get_MCP_data(v.model)
    for i in 1:length(mcp_data)
        if mcp_data[i].var == v
            mcp_data.result_value = value
            break
        end
    end
end

function result_value(v::JuMP.VariableRef)
    mcp_data = get_MCP_data(v.model)
    result_value = NaN
    for i in 1:length(mcp_data)
        if mcp_data[i].var == v
            result_value = mcp_data[i].result_value
            break
        end
    end
    return result_value
end

# https://stackoverflow.com/questions/50084877/how-to-alias-a-macro-in-julia#comment87277306_50085297
@eval const $(Symbol("@mapping")) = $(Symbol("@NLexpression"))

function add_complementarity(m::JuMP.Model, var::JuMP.VariableRef, F::JuMP.NonlinearExpression, F_name::String)

  lb = JuMP.has_lower_bound(var) ? JuMP.lower_bound(var) : -Inf
  ub = JuMP.has_upper_bound(var) ? JuMP.upper_bound(var) : Inf
  var_name = JuMP.name(var)
  new_dimension = ComplementarityType(lb, var, ub, F, raw_index(var), var_name, F_name, NaN)
  mcp_data = get_MCP_data(m)
  push!(mcp_data, new_dimension)
end

macro complementarity(m, F, var)
  F_base_name = string(F)
  F_sym = string(F)
  var_sym = string(var)

  m = esc(m)
  F = esc(F)
  var = esc(var)

  quote
    # if isa($F, JuMP.JuMPArray) || isa($F, Array)
    #   @assert length($F) == length($var)
    # end

    # when var is a single JuMP variable
    if isa($var, JuMP.VariableRef)
      ex_var = Meta.parse(JuMP.name($var))
      ex_F = Meta.parse($F_base_name)

      if typeof(ex_var) == Symbol || typeof(ex_F) == Symbol
        add_complementarity($m, $var, $F, $F_base_name)

      elseif typeof(ex_var) == Expr && typeof(ex_F) == Expr
        @assert ex_var.head == :ref
        @assert ex_F.head == :ref

        ex_F.args[2] = ex_var.args[2]
        F_name = string(ex_F)

        add_complementarity($m, $var, $F, F_name)

      else
        error("Error in @complementarity. Please file an issue with an example.")
      end

    # when var is a single dimensional Array of JuMP.Variable
    elseif isa($var, Array{JuMP.VariableRef})
      for idx in CartesianIndices(size($var))
        idx_name = idx.I
        F_name = string($F_base_name, "[", idx_name, "]")
        var_idx = $var[idx]
        F_idx = $F[idx]

        add_complementarity($m, var_idx, F_idx, F_name)
      end

    else # isa($var, JuMP.JuMPArray) && length(($var).indexsets) == 1 or > 1
      # when var is a multi-dimensional JuMP variable array
      ex = :(Base.product())
      for i in 1:length($var.axes)
        push!(ex.args, $var.axes[i])
      end
      idx_list = collect(eval(ex))

      for idx in idx_list
        F_name = string($F_base_name, "[", join(idx,","), "]")

        var_idx = $var[idx...]
        F_idx = $F[idx...]

        add_complementarity($m, var_idx, F_idx, F_name)
      end

    end # end of if

  end # end of quote

end # end of @complementarity

#


text:
using Test, CharFuncPricing
include("hestonnemo.jl")
include("hestoncum.jl")
include("schobelzhu.jl")
include("hestonfloat.jl")
#include("hestonbench.jl")
#Testing https://financepress.com/2019/02/15/heston-model-reference-prices/


text:
export iauLtp
"""
Long-term precession matrix.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   epj     double         Julian epoch (TT)

Returned:
   rp      double[3][3]   precession matrix, J2000.0 to date

Notes:

   1. The matrix is in the sense

      P_date = rp x P_J2000,

   where P_J2000 is a vector with respect to the J2000.0 mean
   equator and equinox and P_date is the same vector with respect to
   the equator and equinox of epoch epj.

   2. The Vondrak et al. (2011, 2012) 400 millennia precession model
   agrees with the IAU 2006 precession at J2000.0 and stays within
   100 microarcseconds during the 20th and 21st centuries.  It is
   accurate to a few arcseconds throughout the historical period,
   worsening to a few tenths of a degree at the end of the
   +/- 200,000 year time span.

Called:
   iauLtpequ    equator pole, long term
   iauLtpecl    ecliptic pole, long term
   iauPxp       vector product
   iauPn        normalize vector

References:

   Vondrak, J., Capitaine, N. and Wallace, P., 2011, New precession
   expressions, valid for long time intervals, Astron.Astrophys. 534,
   A22

   Vondrak, J., Capitaine, N. and Wallace, P., 2012, New precession
   expressions, valid for long time intervals (Corrigendum),
   Astron.Astrophys. 541, C1

This revision:  2015 December 6

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauLtp(double epj, double rp[3][3])
function iauLtp(epj::Real)
   rp = zeros(Float64, 3, 3)
   ccall((:iauLtp, libsofa_c), Cvoid,
         (Cdouble, Ref{Cdouble}),
          convert(Float64, epj),
          rp)

   return SMatrix{3,3}(rp')
end


text:
"""

Julia の Clustering.hclust() の分析結果を水平方向のデンドログラムとして描く。

StatsPlots パッケージの plot() は垂直方向のデンドログラムを描く。

"""

using Plots
function plot_hclust_horizontal(hc)
    function get(i, x)
        x[i] < 0 && return(ord[abs(x[i])])
        get(x[i], x)
    end
    n = length(hc.order)
    apy = collect(n:-1:1) .+ 0.5
    apx = zeros(n)
    ord = sortperm(hc.order)
    plot(yshowaxis=false, yticks=false, tick_direction=:out, grid=false,
         xlims=(-0.05, maximum(hc.height)), ylims=(0, apy[1]),
         xlabel="height", label="")
    for i in 1:n
        annotate!(0, apy[i], text(string(hc.labels[hc.order[i]]) * " ", 8, :right))
    end
    for i in 1:n-1
        c1 = get(i, hc.merge[:,1])
        c2 = get(i, hc.merge[:,2])
        plot!([apx[c1], hc.height[i], hc.height[i], apx[c2]],
              [apy[c1], apy[c1], apy[c2], apy[c2]], color=:black, label="")
        apx[c1] = apx[c2] = hc.height[i]
        apy[c1] = apy[c2] = (apy[c1] + apy[c2]) / 2
    end
    plot!()
end

"""
使用法

julia> using RDatasets

julia> iris = RDatasets.dataset("datasets", "iris");

julia> function distancematrix(X)
           nr,  nc = size(X)
           d = zeros(nr, nr)
           for i in 1:nr-1
               for j in i+1:nr
                   d[i, j] = d[j, i] = sqrt(sum((X[i, :] .- X[j, :]).^2))
               end
           end
           d
       end;

julia> # using Random; Random.seed!(123) # 毎回同じ結果にするためには乱数の種を設定する
       using Clustering

julia> x = Matrix(iris[1:20, 1:4]);

julia> D = distancematrix(x);

julia> hc = hclust(D, linkage=:ward);

julia> plot_hclust_horizontal(hc)

"""

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
using LinearAlgebra

function template!(g, H, x)
end

function trs_hard!(g, cH::Cholesky, Δ, v)
    n = length(g)
    v′ = v .* sqrt(norm(cH)) ./ n
    lowrankupdate!(cH, v′)
    q = ldiv!(cH, g)
    return v .* sqrt(Δ ^ 2 - dot(q, q)) .- q
end

function trs_hard!(g, H::AbstractMatrix, Δ, v)
    n = length(g)
    H′ = mul!(H, v, v', norm(H) / n ^ 2, 1)
    q = ldiv!(H′, g)
    return v .* sqrt(Δ ^ 2 - dot(q, q)) .- q
end

function trs_eig!(trA, g, H; Δ, gaptol=eps()*10^8)
    n = length(g)

    Λ, V = eigen(trA)
    λ1, λ2 = float(Λ[1]), float(Λ[2])
    y2 = float.(V[n+1:end, 1])
    y1 = float.(V[1:n, 1])

    if norm(y1) <= gaptol / sqrt(λ1 - λ2)
        v = y2 / norm(y2)
        return trs_hard!(g, H, Δ, v)
    else
        return -sign(dot(g, y2)) .* Δ .* y1 ./ norm(y1)
    end
end

function trust_region_sub!(g, H; Δ, gaptol=eps()*10^8)
    n = length(g)

    w = g ./ Δ
    trA = [H Matrix(1.0 * I, n, n); w * w' H]

    cH = cholesky(H; check=false)
    if issuccess(cH)
        p = cH \ g
        if norm(p) <= Δ
            return p, false
        else
            return trs_eig!(trA, g, cH; Δ=Δ, gaptol=gaptol), true
        end
    end
    return trs_eig!(trA, g, H; Δ=Δ, gaptol=gaptol), true
end

function _tr_interp!()
end

function trust_region_newton(
    x0, feval!;
    maxiter=100, Δmax=Inf, monitor=(x, rnorm, Δ) -> nothing
)
    n = length(x0)
    x = copy(x0)
    g = similar(x)
    H = similar(x, n, n)

    f = feval!((), g, H, x)

    cH = cholesky(H; check=false)
    nH = opnorm(H)
    p = -H \ g

    Δ = 1.2 * dot(p, p)
    active_constraint = false
    return nothing
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
type NaiveBayes
    n_class::Integer
    class_mean::Matrix
    class_var::Matrix
    class_priors::Vector
    class_::Vector
end

function NaiveBayes(;
                    class_mean=zeros(2,2),
                    class_var=zeros(2,2),
                    class_priors=zeros(2),
                    n_class = 2,
                    class_ = zeros(2))
    return NaiveBayes(n_class, class_mean, class_var, class_priors,class_)
end

function train!(model::NaiveBayes, X::Matrix, y::Vector)
    n_feature = size(X,2)
    n_sample = size(X,1)
    model.class_ = sort(unique(y))
    model.n_class = length(unique(y))
    model.class_mean = zeros(n_feature, model.n_class)
    model.class_var = zeros(n_feature, model.n_class)
    model.class_priors = zeros(model.n_class)

    for (i,class) in enumerate(model.class_)
        X_temp = X[vec(y.==class),:]
        model.class_mean[:,i] = mean(X_temp,1)
        model.class_var[:,i] = var(X_temp,1)
        model.class_priors[i] = size(X_temp,1)/n_sample
    end

end

function predict(model::NaiveBayes,
                 X::Matrix)
    n_sample = size(X,1)
    res = zeros(n_sample)
    for i = 1:n_sample
        res[i] = predict(model, vec(X[i,:]))
    end
    return res
end

function predict(model::NaiveBayes,
                 X::Vector)
    temp = zeros(model.n_class)
    for i = 1:model.n_class
        prior = log(model.class_priors[i])
        cond = sum(log_pdf(model, X, i))
        temp[i] = prior + cond
    end
    res = softmax(temp)
    class = model.class_[res]
    return class
end

function log_pdf(model::NaiveBayes,
                 X::Vector,n::Integer)

    mean_ = vec(model.class_mean[:,n])
    var_ = vec(model.class_var[:, n])
    d = sqrt(2*pi*var_)
    n = exp(-((X-mean_).^2./(2*var_)))
    return log(n./d)
end

function test_naive()
    X_train, X_test, y_train, y_test = make_cla()
    model = NaiveBayes()
    train!(model,X_train, y_train)
    predictions = predict(model,X_test)
    print("classification accuracy", accuracy(y_test, predictions))

    #PCA

    pca_model = PCA()
    train!(pca_model, X_test)
    plot_in_2d(pca_model, X_test, predictions, "Naive Bayes")
end


text:
import ..CausalityTests: JointDistanceDistributionTTest

################################################################
# Integration with UncertainData.jl
################################################################
function causality(source::Vector{<:UVT1}, target::Vector{<:UVT2}, p::JointDistanceDistributionTTest) where {
        UVT1<:AbstractUncertainValue,
        UVT2<:AbstractUncertainValue}
    causality(resample.(source), resample.(target), p)
end

function causality(source, target::Vector{<:UVT}, p::JointDistanceDistributionTTest) where {
        UVT<:AbstractUncertainValue}
    causality(source, resample.(target), p)
end

function causality(source::Vector{<:UVT}, target, p::JointDistanceDistributionTTest) where {
        UVT<:AbstractUncertainValue}
    causality(resample.(source), target, p)
end

function causality(source::UVT1, target::UVT2, p::JointDistanceDistributionTTest) where {
        UVT1<:AbstractUncertainValueDataset,
        UVT2<:AbstractUncertainValueDataset}
    causality(resample.(source), resample.(target), p)
end

function causality(source, target::UVT, p::JointDistanceDistributionTTest) where {
        UVT<:AbstractUncertainValueDataset}
    causality(source, resample.(target), p)
end

function causality(source::UVT, target, p::JointDistanceDistributionTTest) where {
        UVT<:AbstractUncertainValueDataset}
    causality(resample.(source), target, p)
end


text:
# Test Rand index

using Test
using Clustering

@testset "randindex() (Rand index)" begin

a1 = [1, 1, 1, 2, 2, 2, 3, 3, 3, 3]
a2 = [1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
a3 = [3, 3, 3, 2, 2, 2, 1, 1, 1, 1]

(ARI, RI, MI, HI) = randindex(a1, a1)

@test ARI ≈ 1.0 atol=1.0e-12
@test RI  ≈ 1.0 atol=1.0e-12
@test MI  ≈ 0.0 atol=1.0e-12
@test HI  ≈ 1.0 atol=1.0e-12

(ARI, RI, MI, HI) = randindex(a1, a3)

@test ARI ≈ 1.0 atol=1.0e-12
@test RI  ≈ 1.0 atol=1.0e-12
@test MI  ≈ 0.0 atol=1.0e-12
@test HI  ≈ 1.0 atol=1.0e-12

(ARI, RI, MI, HI) = randindex(a1, a2)

@test ARI ≈ 0.403669 atol=1.0e-5
@test RI  ≈ 0.711111 atol=1.0e-5
@test MI  ≈ 0.288888 atol=1.0e-5
@test HI  ≈ 0.422222 atol=1.0e-5

@test randindex(a1, a2) == randindex(a2, a1)

end


text:
using Plots
using ControlSystems, RobustAndOptimalControl
"""
This is a simple SISO example with a pole in the origin, corresponding to the
DC servos used in the Lund laboratories. It serves to exeplify how the syntheis
can be done for simple SISO systems, and also demonstrates how we chan verify
if the problem is feasible to solve using the ARE method.

The example can be set to visualize plots using the variables
  makeplots - true/false (true if plots are to be generated, false for testing)
"""
makeplots = true

# Define the process
Gtrue   = tf([11.2], [1, 0.12,0])

# Sensitivity weight function
M, wB, A = 1.5, 20.0, 1e-8
WS = tf([1/M, wB],[1, wB*A])

# Output sensitivity weight function
WU = ss(1.0)

# Complementary sensitivity weight function
WT = []

# Form the P in the LFT F_l(P,C) as a partitioned state-space object
P = hinfpartition(Gtrue, WS, WU, WT)

# Check if the system is feasible for synythesis
flag = hinfassumptions(P, verbose=false)

# Since it is not, modify the plant desciption
epsilon = 1e-5
G = tf([11.2], [1, 0.12]) * tf([1], [1, epsilon])

# Form the P in the LFT Fl(P,C) as a partitioned state-space object
P = hinfpartition(G, WS, WU, WT)

# Check if the problem is feasible
flag = hinfassumptions(P)

# Synthesize the H-infinity optimal controller
flag, C, γ = hinfsynthesize(P, γrel=1)

# Extract the transfer functions defining some signals of interest
Pcl, S, CS, T = hinfsignals(P, G, C)

## Plot the specifications
if makeplots
  specificationplot([S, CS, T], [WS, WU, WT], γ) |> display
  ## Plot the closed loop gain from w to z
  specificationplot(Pcl, γ; s_labels=["\$\\sigma(P_{cl}(j\\omega))\$"], w_labels=["\$\\gamma\$"])
  ylims!((0.1, 10))

end


text:
function load_xyz_data(file_path :: String; atom_names = nothing, radii = [3.5], weights = [1.0], masses = [1.0], boundary_type = ["p", "p", "p"], units = "lj")
    # frames = read_frames(file_path)
    # n = length(frames)
    # r = Vector{Configuration}(undef, n)

    # if atom_names == nothing
    #     atom_names = unique(frames[1]["arrays"]["species"])
    # end
    # num_atom_types = length(atom_names)

    # if length(radii) != num_atom_types
    #     radii = radii[1] .* ones(num_atom_types)
    # end
    # if length(weights) != num_atom_types
    #     weights = weights .* ones(num_atom_types)
    # end
    # if length(masses) !=  num_atom_types
    #     masses = masses .* ones(num_atom_types)
    # end

    # for i = 1:n
    #     frame = frames[i]
    #     num_atoms = frame["N_atoms"]
    #     Atoms = Vector{Atom}(undef, num_atoms)
    #     for j = 1:num_atoms
    #         arrays = frame["arrays"]
    #         name = arrays["species"][j]
    #         id = findall(name .== atom_names)[1]
    #         pos = arrays["pos"][:, j]
    #         if "vel" in keys(arrays)
    #             vel = arrays[vel][:, j]
    #         else
    #             vel = zeros(3)
    #         end
    #         a = Atom(masses[id], pos, vel, name)
    #         Atoms[j] = a

    #     c = Configuration(Atoms, num_atom_types,
    #                       Angles, num_angle_types,
    #                       Bonds, num_bond_types,
    #                       Impropers, num_improper_types,
    #                       Dihedrals, num_dihedral_types,
    #                       atom_names, masses, radii, weights,
    #                       domain, units)
    #     r[i] = c
    #     end
    # end

    println("This is not yet implemented")
end

function save_xyz_data(c::Configuration, file_path::String)
    println("This is not yet implemented")
end


text:
using InteractiveChaos, Makie, OrdinaryDiffEq
diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)

hh = Systems.henonheiles()

potential(x, y) = 0.5(x^2 + y^2) + (x^2*y - (y^3)/3)
energy(x,y,px,py) = 0.5(px^2 + py^2) + potential(x,y)
const E = energy(get_state(hh)...)

function complete(y, py, x)
    V = potential(x, y)
    Ky = 0.5*(py^2)
    Ky + V ≥ E && error("Point has more energy!")
    px = sqrt(2(E - V - Ky))
    ic = [x, y, px, py]
    return ic
end

plane = (1, 0.0) # first variable crossing 0

# %%

state, scene = interactive_poincaresos(hh, plane, (2, 4), complete;
markersizes = (-5, -1), labels = ("q₂" , "p₂"), diffeq...);

# lines!(scene.children[2], [Point2f0(0.06465084105730057, -0.4 + 0.2*i) for i in 0:4]);
# lines!(scene.children[2], [Point2f0(-0.3 + 0.15*i, -0.18481573462486267) for i in 0:6]);
# display(scene);

# %% Coloring points using a custom function
# Here I use the first momentum
momentum1(u) = RGBf0((0.5*u[3]^2)/E, 0, 0)

state, scene = interactive_poincaresos(hh, plane, (2, 4), complete;
markersizes = (-5, -1), color = momentum1, diffeq...)

# %% Coloring points using the Lyapunov exponent
function λcolor(u)
    λ = lyapunovs(hh, 4000; u0 = u)[1]
    λmax = 0.1
    return RGBf0(0, 0, clamp(λ/λmax, 0, 1))
end

state, scene = interactive_poincaresos(hh, plane, (2, 4), complete;
markersizes = (-5, -1), labels = ("q₂" , "p₂"), color = λcolor, diffeq...)


text:
function test_crmr()
  crmr_tol = 1.0e-4;  # We're tolerant just so random tests don't fail.

  function test_crmr(A, b; λ=0.0, M=opEye(size(A,1)))
    (nrow, ncol) = size(A);
    (x, stats) = crmr(A, b, λ=λ, M=M);
    r = b - A * x;
    if λ > 0
      s = r / sqrt(λ);
      r = r - sqrt(λ) * s;
    end
    resid = norm(r) / norm(b);
    @printf("CRMR: residual: %7.1e\n", resid);
    return (x, stats, resid)
  end

  # Underdetermined consistent.
  A = rand(10, 25); b = A * ones(25);
  (x, stats, resid) = test_crmr(A, b);
  @test(resid <= crmr_tol);
  @test(stats.solved);
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x);
  @test(norm(xI - xmin) <= cond(A) * crmr_tol * xmin_norm);

  # Underdetermined inconsistent.
  A = ones(10, 25); b = rand(10); b[1] = -1.0;
  (x, stats, resid) = test_crmr(A, b);
  @test(stats.inconsistent);
  @test(stats.Aresiduals[end] <= crmr_tol);

  # Square consistent.
  A = rand(10, 10); b = A * ones(10);
  (x, stats, resid) = test_crmr(A, b);
  @test(resid <= crmr_tol);
  @test(stats.solved);
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x);
  @test(norm(xI - xmin) <= cond(A) * crmr_tol * xmin_norm);

  # Square inconsistent.
  A = ones(10, 10); b = rand(10); b[1] = -1.0;
  (x, stats, resid) = test_crmr(A, b);
  @test(stats.inconsistent);
  @test(stats.Aresiduals[end] <= crmr_tol);

  # Overdetermined consistent.
  A = rand(25, 10); b = A * ones(10);
  (x, stats, resid) = test_crmr(A, b);
  @test(resid <= crmr_tol);
  @test(stats.solved);
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x);
  @test(norm(xI - xmin) <= cond(A) * crmr_tol * xmin_norm);

  # Overdetermined inconsistent.
  A = ones(5, 3); b = rand(5); b[1] = -1.0;
  (x, stats, resid) = test_crmr(A, b);
  @test(stats.inconsistent);
  @test(stats.Aresiduals[end] <= crmr_tol);

  # With regularization, all systems are underdetermined and consistent.
  (x, stats, resid) = test_crmr(A, b, λ=1.0e-3);
  @test(resid <= crmr_tol);
  @test(stats.solved);
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x, λ=1.0e-3);
  @test(norm(xI - xmin) <= cond(A) * crmr_tol * xmin_norm);

  # Code coverage.
  (x, stats, resid) = test_crmr(sparse(A), b, λ=1.0e-3);
  show(stats);

  # Test b == 0
  (x, stats) = crmr(A, zeros(size(A,1)))
  @test x == zeros(size(A,2))
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = [I; rand(1:10, 2, 3)]
  b = A * ones(Int, 3)
  (x, stats) = crmr(A, b)
  @test stats.solved

  # Test preconditioner with an under-determined problem:
  # Find the least norm force that transfers mass unit distance with zero final velocity
  A = 0.5 * [19.0 17.0 15.0 13.0 11.0 9.0 7.0 5.0 3.0 1.0;
              2.0  2.0  2.0  2.0  2.0 2.0 2.0 2.0 2.0 2.0]
  b = [1.0; 0.0]
  M = LinearOperator(Diagonal(1 ./ (A * A')))
  (x, stats, resid) = test_crmr(A, b, M=M)
  @test(resid ≤ crmr_tol)
  @test(stats.solved)
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x)
  @test(norm(xI - xmin) ≤ cond(A) * crmr_tol * xmin_norm)
end

test_crmr()


text:
# indicator of a simplex

"""
  IndSimplex(a::Union{Real, Integer}=1.0)

Returns the function `g = ind{x : x ⩾ 0, sum(x) = a}`.
"""

immutable IndSimplex{T <: Union{Real, Integer}} <: IndicatorConvex
  a::T
  function IndSimplex(a::T)
    if a <= 0
      error("parameter a must be positive")
    else
      new(a)
    end
  end
end

IndSimplex{T <: Union{Real, Integer}}(a::T=1.0) = IndSimplex{T}(a)

function (f::IndSimplex){T <: Real}(x::AbstractArray{T,1})
  if all(x .>= 0) && abs(sum(x)-f.a) <= 1e-14
    return 0.0
  end
  return +Inf
end

function prox!{T <: Real}(f::IndSimplex, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
# Implements Algorithm 1 in Condat, "Fast projection onto the simplex and the l1 ball", Mathematical Programming, 158:575–585, 2016.
# We should consider implementing the other algorithms reviewed there, and the one proposed in the paper.
  n = length(x)
  p = []
  if ndims(x) == 1
    p = sort(x, rev=true)
  else
    p = sort(x[:], rev=true)
  end
  s = 0
  for i = 1:n-1
    s = s + p[i]
    tmax = (s - f.a)/i
    if tmax >= p[i+1]
      @inbounds for j in eachindex(y)
        y[j] = x[j] < tmax ? 0.0 : x[j] - tmax
      end
      return 0.0
    end
  end
  tmax = (s + p[n] - f.a)/n
  @inbounds for j in eachindex(y)
    y[j] = x[j] < tmax ? 0.0 : x[j] - tmax
  end
  return 0.0
end

fun_name(f::IndSimplex) = "indicator of the probability simplex"
fun_dom(f::IndSimplex) = "AbstractArray{Real}"
fun_expr(f::IndSimplex) = "x ↦ 0 if x ⩾ 0 and sum(x) = a, +∞ otherwise"
fun_params(f::IndSimplex) = "a = $(f.a)"

function prox_naive{T <: Real}(f::IndSimplex, x::AbstractArray{T}, gamma::Real=1.0)
  low = minimum(x)
  upp = maximum(x)
  v = x
  s = Inf
  for i = 1:100
    if abs(s)/f.a <= 1e-15
      break
    end
    alpha = (low+upp)/2
    v = max.(x - alpha, 0.0)
    s = sum(v) - f.a
    if s <= 0
      upp = alpha
    else
      low = alpha
    end
  end
  return v, 0.0
end


text:
using LightGraphs
using Random
using GraphPlot
import Cairo, Fontconfig
using Compose
using Colors

#exercice 1
function plot_graphs(nb_nodes=400, nb_clusters=4, q=0.1)
    p = [0.1, 0.3, 0.5, 0.9]
    for p_i in p
        g, y = generate_graph(p_i, q, nb_nodes, nb_clusters)
        draw(PNG(string("../res/", p_i, "_clusters.png"), 16cm, 16cm), gplot(g))
    end
end

function generate_graph(p::Float64, q::Float64, nb_nodes::Int64, nb_clusters::Int64)
    clusters = Vector{Int64}()
    for j in 1:nb_clusters
        for i in 1:Int(floor(nb_nodes//nb_clusters))
            push!(clusters, j)
        end
    end

    edges = Vector{Tuple{Int64, Int64}}()
    for i1 in 1:nb_nodes
        for i2 in i1:nb_nodes
            r = rand()
            if clusters[i1] == clusters[i2]
                if r < p
                    push!(edges, (i1, i2))
                end
            elseif r < q
                push!(edges, (i1, i2))
            end
        end
    end

    g = SimpleGraph(nb_nodes)
    for edge in edges
        add_edge!(g, edge[1], edge[2])
    end
    return g, clusters
end

plot_graphs()

#exercice 2

function fy_shuffle!(tab)
    #Fisher-Yates shuffle
    list_range = [1:1:length(tab);]
    for i in list_range
        j = rand(list_range[1]: list_range[end])
        tab[i], tab[j] = tab[j], tab[i]
    end
    return(tab)
end

function neighbour_label(graph, node, labels)
    #find label occurring with the highest frequency among neighbours of node in graph
    unique_labels = unique(labels)
    values = zeros(length(unique_labels))
    scores = Dict(zip(unique_labels, values))

    for neighb in  neighbors(graph, node)
        neighb_label = labels[neighb]
        scores[neighb_label] += 1
    end
    a = argmax(scores)

    return a
end

function label_propagation_step_three(graph,labels)

    new_labels =zeros(length(labels))
    graph_nodes = [1:1:nv(graph);]
    shuffled_nodes = fy_shuffle!(graph_nodes)

    #for each node in the network (in this random order)
    # set its label to a label occurring with the highest frequency among its neighbours
    for node in shuffled_nodes
        new_labels[node] = neighbour_label(graph,node,labels)
    end

    return new_labels
end

function label_propagation!(graph, init_labels, max_iter=1000)

    test = true #bool that indicates if the we keep iterating the algorithm
    n_iter = 0
    labels = init_labels

    while test
        #step three: label propagation
        labels = label_propagation_step_three(graph,labels)

        #stoping criteria
        test = false
        for node in vertices(graph)
            # Check if one node doesn't belong to the same communities than the majority of its neighbors
            if labels[node] != neighbour_label(graph,node,labels)
                test = true# If one node doesn't satisfy the condition, we need to continue
            end
        end

        if n_iter > max_iter
            test = false
        end

        n_iter+=1
    end

    return labels
end

function main_label_propagation(nb_nodes=400, nb_clusters=4, q=0.1)
    p = [0.1, 0.3, 0.5, 0.9]
#     p = [0.1]
    for p_i in p
        graph, true_y = generate_graph(p_i, q, nb_nodes, nb_clusters)

        ##
        init_y = [1:1:nv(graph);]
        init_y = fy_shuffle!(init_y)

        label_prop_y = label_propagation!(graph, init_y)
        membership = unique(label_prop_y)
        chosen=[]
        for l in label_prop_y
            f = findall(membership .== l)
            push!(chosen, f[1])
        end

        nodecolor = [colorant"lightseagreen", colorant"orange", colorant"blue", colorant"purple"]
        nodefillc = nodecolor[chosen]
        draw(PNG(string("../res/", p_i, "_propagation.png"), 16cm, 16cm),gplot(graph, nodefillc=nodefillc))

	nodefillc2 = nodecolor[true_y]
	draw(PNG(string("../res/", p_i, "_true.png"), 16cm, 16cm),gplot(graph, nodefillc=nodefillc2))
    end
end

main_label_propagation()

max_stars_repo_name: quentin-burthier/GrandsGraphes

max_issues_repo_name: quentin-burthier/GrandsGraphes

max_forks_repo_name: quentin-burthier/GrandsGraphes


text:
"""
Fisher information for magnetometry with continuously monitored spin systems, with independent Markovian noise acting on each spin
"""
module ContinuousMeasurementFI

    using SparseArrays
    using LinearAlgebra

    export Eff_QFI, Molmer_QFI_GHZ, Molmer_qfi_transverse, uncond_QFI_transverse

    include("NoiseOperators.jl")
    include("States.jl")
    include("Fisher.jl")
    include("Molmer_QFI.jl")
    include("Uncond_qfi_transverse.jl")

    include("Eff_QFI_PD.jl")
    include("Eff_QFI_PD_pure.jl")
    include("Eff_QFI_HD.jl")
    include("Eff_QFI_HD_pure.jl")

    """
        (t, FI, QFI) = Eff_QFI(kwargs...)

    Evaluate the continuous-time FI and QFI of a final strong measurement for the
    estimation of the frequency ω with continuous monitoring of each half-spin
    particle affected by noise at an angle θ, with efficiency η using SME
    (stochastic master equation).

    The function returns a tuple `(t, FI, QFI)` containing the time vector and the
    vectors containing the FI and average QFI

    # Arguments

    * `Nj`: number of spins
    * `Ntraj`: number of trajectories for the SSE
    * `Tfinal`: final time of evolution
    * `measurement = :pd` measurement (either `:pd` or `:hd`)
    * `dt`: timestep of the evolution
    * `κ = 1`: the noise coupling
    * `θ = 0`: noise angle (0 parallel, π/2 transverse)
    * `ω = 0`: local value of the frequency
    * `η = 1`: measurement efficiency
    """
    function Eff_QFI(args::Dict{String, Any})
        kwargs = map(x -> (Symbol(x[1]), x[2]), collect(args))
        @show kwargs
        tic()
        (t, fi, qfi) = Eff_QFI(; kwargs...)
        eval_time = toc()
        return merge(args, Dict("t" => t, "fi" => fi,"qfi" => qfi, "eval_time" => eval_time))
    end

    function Eff_QFI(; Nj=1,        # Number of spins
                    Ntraj=1,        # Number of trajectories
                    Tfinal=1.,      # Final time
                    dt=0.01,        # Time step
                    measurement=:pd,# Measurement type
                    η = 1.,         # Measurement efficiency
                    κ = 1.,         # Field rate
                    θ = pi/2,       # Angle of the noise
                    ω = 1.,         # Magnetic field
                    _watherver...)

            if measurement == :pd
                if η == 1.
                    f = Eff_QFI_PD_pure
                else
                    f = Eff_QFI_PD
                end
            elseif measurement == :hd
                if η == 1.
                    f = Eff_QFI_HD_pure
                else
                    f = Eff_QFI_HD
                end
            else
                throw(ArgumentError("Measurement must be either :pd or :hd"))
            end

            if η == 1.
                f(Nj, Ntraj, Tfinal, dt;
                    κ = κ,
                    θ = θ,
                    ω = ω)
            else
                f(Nj, Ntraj, Tfinal, dt;
                κ = κ,
                θ = θ,
                ω = ω,
                η = η)
            end

    end
end


text:
function poissonpvec(μ::Float64, n::Int)
    # Poisson probabilities, from 0 to n
    pv = Vector{Float64}(undef, n+1)
    @inbounds pv[1] = p = exp(-μ)
    for i = 1:n
        @inbounds pv[i+1] = (p *= (μ / i))
    end
    return pv
end

# Naive sampler by counting exp variables
#
# Suitable for small μ
#
struct PoissonCountSampler{T<:Real} <: Sampleable{Univariate,Discrete}
    μ::T
end

PoissonCountSampler(d::Poisson) = PoissonCountSampler(rate(d))

function rand(rng::AbstractRNG, s::PoissonCountSampler)
    μ = s.μ
    T = typeof(μ)
    n = 0
    c = randexp(rng, T)
    while c < μ
        n += 1
        c += randexp(rng, T)
    end
    return n
end

# Algorithm from:
#
#   J.H. Ahrens, U. Dieter (1982)
#   "Computer Generation of Poisson Deviates from Modified Normal Distributions"
#   ACM Transactions on Mathematical Software, 8(2):163-179
#
#   For μ sufficiently large, (i.e. >= 10.0)
#
struct PoissonADSampler{T<:Real} <: Sampleable{Univariate,Discrete}
    μ::T
    s::T
    d::T
    L::Int
end

PoissonADSampler(d::Poisson) = PoissonADSampler(rate(d))

function PoissonADSampler(μ::Real)
    s = sqrt(μ)
    d = 6 * μ^2
    L = floor(Int, μ - 1.1484)

    PoissonADSampler(promote(μ, s, d)..., L)
end

function rand(rng::AbstractRNG, sampler::PoissonADSampler)
    μ = sampler.μ
    s = sampler.s
    d = sampler.d
    L = sampler.L
    μType = typeof(μ)

    # Step N
    G = μ + s * randn(rng, μType)

    if G >= zero(G)
        K = floor(Int, G)
        # Step I
        if K >= L
            return K
        end

        # Step S
        U = rand(rng, μType)
        if d * U >= (μ - K)^3
            return K
        end

        # Step P
        px, py, fx, fy = procf(μ, K, s)

        # Step Q
        if fy * (1 - U) <= py * exp(px - fx)
            return K
        end
    end

    while true
        # Step E
        E = randexp(rng, μType)
        U = 2 * rand(rng, μType) - one(μType)
        T = 1.8 + copysign(E, U)
        if T <= -0.6744
            continue
        end

        K = floor(Int, μ + s * T)
        px, py, fx, fy = procf(μ, K, s)
        c = 0.1069 / μ

        # Step H
        if c*abs(U) <= py*exp(px + E) - fy*exp(fx + E)
            return K
        end
    end
end

# Procedure F
function procf(μ, K::Int, s)
    # can be pre-computed, but does not seem to affect performance
    ω = 0.3989422804014327/s
    b1 = 0.041666666666666664/μ
    b2 = 0.3*b1*b1
    c3 = 0.14285714285714285*b1*b2
    c2 = b2 - 15 * c3
    c1 = b1 - 6 * b2 + 45 * c3
    c0 = 1 - b1 + 3 * b2 - 15 * c3

    if K < 10
        px = -μ
        py = μ^K/factorial(K)
    else
        δ = 0.08333333333333333/K
        δ -= 4.8*δ^3
        V = (μ-K)/K
        px = K*log1pmx(V) - δ # avoids need for table
        py = 0.3989422804014327/sqrt(K)

    end
    X = (K-μ+0.5)/s
    X2 = X^2
    fx = -X2 / 2 # missing negation in pseudo-algorithm, but appears in fortran code.
    fy = ω*(((c3*X2+c2)*X2+c1)*X2+c0)
    return px,py,fx,fy
end