text:
###### This file is part of the MendelImpute.jl package.
###### It contains relevant code to perform phasing by dynamic programming,
###### which finds the global solution to the minimum number of phase switches.

"""
Helper function to calculate the squared Hamming difference between 2
unordered pair of integers.

# Inputs
- `pair1`: tuple (a, b) where a, b are integers
- `pair2`: tuple (c, d) where c, d are integers

# Optional Inputs
- `λ`: Error each switch contributes. Defaults to 1.0

# Examples
- `pair_error((1, 2), (2, 3)) = pair_error((2, 1), (2, 3)) = λ`
- `pair_error((2, 5), (5, 2)  = 0`
- `pair_error((1, 2), (3, 4)) = 4λ`
"""
function pair_error(
    pair1::T,
    pair2::T;
    λ::Real = 1.0
    ) where T <: Tuple{Int32, Int32}
    # parallel connections
    # a b
    # | |
    # c d
    parallel_diff = (pair1[1] != pair2[1]) + (pair1[2] != pair2[2])
    # a b
    #  X
    # c d
    crossover_diff = (pair1[1] != pair2[2]) + (pair1[2] != pair2[1])
    return λ * abs2(min(parallel_diff, crossover_diff))
end

"""
Finds the optimal sequence of haplotype pairs across all windows
such that number of switch points is minimized. Windows with
too little typed SNPs will inherit haplotype pairs from the closest
window containing feasible haplotype pairs.

# Inputs
- `haplotype_set`: A person's possible haplotype pairs in each window.

# Optional input:
- `λ`: Error each switch contributes. Defaults to 1.0

# Output:
- `sol_path`: Optimal sequence of haplotype pairs indices (in the complete
    haplotype pool) in each window
- `memory`: Vector of dictionary storing the optimal error for each haplotype
    pair in each window
- `path_err`: Error for each window induced by `sol_path`
- `best_err`: Osverall error induced by `sol_path`. Equals λ times number of
    switch points.
"""
function connect_happairs!(
    haplotype_set::Vector{Vector{T}};
    λ::Float64 = 1.0
    ) where T <: Tuple{Int32, Int32}

    # allocate working arrays
    windows  = length(haplotype_set)
    sol_path = Vector{T}(undef, windows)
    next_pair = [Int32[] for i in 1:windows]
    subtree_err = [Float64[] for i in 1:windows]

    # computational routine
    best_err = connect_happairs!(sol_path, next_pair, subtree_err,
        haplotype_set, λ = λ)

    return sol_path, next_pair, subtree_err, best_err
end

"""
Finds the optimal sequence of haplotype pairs across all windows
such that number of switch points is minimized. Windows with
too little typed SNPs will inherit haplotype pairs from the closest
window containing feasible haplotype pairs.

# Inputs
- `sol_path`: Optimal sequence of haplotype pairs indices (in the complete
    haplotype pool) in each window
- `memory`: Vector of dictionary storing the optimal error for each haplotype
    pair in each window
- `path_err`: Error for each window induced by `sol_path`
- `haplotype_set`: A vector of vectors. `haplotype_set[1]` stores all pairs of
    haplotypes in window 1 in a vector, and so on.
- `λ`: Error each switch contributes. Defaults to 1.0

# Optional input:
- `λ`: Error each switch contributes. Defaults to 1.0
- `tolerance`: A heuristic to reduce search space. At each window, pairs that
    are more than `tolerance` worse than the optimal pair gets deleted.

# Output
- `best_err`: Overall error induced by `sol_path`. Equals λ times number of
    switch points.
"""
function connect_happairs!(
    sol_path::Vector{T},
    next_pair::Vector{Vector{Int32}},
    subtree_err::Vector{Vector{Float64}},
    haplotype_set::Vector{Vector{T}};
    λ::Float64 = 1.0,
    tolerance::Float64 = Inf,
    ) where T <: Tuple{Int32, Int32}

    windows = length(haplotype_set)

    # reset storage
    empty!.(next_pair)
    empty!.(subtree_err)

    # base case: last window induces no error and connects to nothing
    @inbounds for pair in haplotype_set[windows]
        push!(next_pair[windows], 0)
        push!(subtree_err[windows], 0.0)
    end

    # search for best haplotype pair in each window bottom-up
    @inbounds for w in Iterators.reverse(1:(windows - 1))
        win_best_err = Inf

        # first pass to compute each pair's optimal error
        for (j, happair) in enumerate(haplotype_set[w])
            # search all pairs in next window
            best_err = Inf
            best_next_pair = 0
            for (i, nextpair) in enumerate(haplotype_set[w + 1])
                err = pair_error(happair, nextpair) + subtree_err[w + 1][i]
                if err < best_err
                    best_err = err
                    best_next_pair = i
                end
            end
            if best_err < win_best_err
                win_best_err = best_err
            end
            push!(subtree_err[w], best_err)
            push!(next_pair[w], best_next_pair)
        end

        # Heuristic to reduce dynamic programming search space
        tol = win_best_err + tolerance
        for (i, err) in enumerate(subtree_err[w])
            if err > tol
                deleteat!(subtree_err[w], i)
                deleteat!(next_pair[w], i)
                deleteat!(haplotype_set[w], i)
            end
        end
    end

    # find best solution path by forward-tracing
    best_err, cur_idx = findmin(subtree_err[1])
    @inbounds for w in 1:windows
        sol_path[w] = haplotype_set[w][cur_idx]
        cur_idx = next_pair[w][cur_idx]
    end

    return best_err
end

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
using StaticArrays
using DocStringExtensions
using Optim

#--- Gaussian Processes

"""
An abstract Gaussian Process type. Every child type should implement the
following functions.

    `?`
"""
abstract type AbstractGP{T <: Number} end

export AbstractGP

#--- Generic functions used by multiple types

"""
Set the hyperparameters of an object.
"""
function sethyperparameters!
end

"""
Return the hyperparameters of an object.
"""
function gethyperparameters
end

"""
Return the number of hyperparameters an object has.
"""
function hyperparametercount
end

"""
Optimize the hyperparameters of a GP using Optim.jl.
"""
function optimizehyperparameters!
end

"""
`loglikelihood(gp)`

Return the negative log likelihood for a GP. Used in optimizing the
hyperparameters.

`loglikelihood(gp, x, y)`

Return the negative log likelihood of an observation given a GP.

"""
function loglikelihood
end

"""
Return the derivative of the negative log likelihood with respect to the
hyperparameters.
"""
function loglikelihood_dθ
end

"""
Return the derivative of the negative log likelihood with respect to the
hyperparameters inplace.
"""
function loglikelihood_dθ!
end

"""
Return the derivative of the negative log likelihood with respect to the
sample positions.
"""
function loglikelihood_dx
end

"""
Return the derivative of the negative log likelihood with respect to the
sample positions inplace.
"""
function loglikelihood_dx!
end

export sethyperparameters!, gethyperparameters, hyperparametercount,
    optimizehyperparameters!, loglikelihood, loglikelihood_dθ,
    loglikelihood_dθ!, loglikelihood_dx, loglikelihood_dx!

#--- Data conversion routines

"""
$(SIGNATURES)

Convert the input to a vector of SVectors.
"""
toSVector(x::Vector{T}) where {T <: Number} = reinterpret(SVector{1, T}, x, (length(x),))
toSVector(x::Vector{SVector{N, T}}) where {N, T <: Number} = x
toSVector(x::Matrix{T}) where {T <: Number} = reinterpret(SVector{size(x, 1), T}, x, (size(x, 2),))
toSVector(x::AbstractVector{T}) where {T <: Number} = [SVector{1}(el) for el in x] # catch-all for generators, etc
toSVector(x::AbstractVector) = toSVector(hcat(x...)) # going to be horrible!

export toSVector

"""
$(SIGNATURES)

Convert the input to a Matrix.
"""
toMatrix(x::Vector{T}) where {T <: Number} = reshape(x, (1, length(x)))
toMatrix(x::Vector{SVector{N, T}}) where {N, T <: Number} = reinterpret(T, x, (N, length(x)))
toMatrix(x::Matrix{T}) where {T <: Number} = x
toMatrix(x::AbstractVector) = hcat(x...)

export toMatrix

#--- Optimization

function optimizehyperparameters!(gp, solver = LBFGS(), options = Optim.Options())
    # Define objective functions with exponential scalings on the
    # hyperparameters as that seems to work best
    f = ϕ -> (θ = exp.(ϕ) ; θ != gethyperparameters(gp) ? sethyperparameters!(gp, θ) : nothing ; loglikelihood(gp))
    g! = (dϕ, ϕ) -> (θ = exp.(ϕ) ; println("ϕ: $(ϕ), θ: $(θ)") ; θ != gethyperparameters(gp) ? sethyperparameters!(gp, θ) : nothing ; loglikelihood_dθ!(dϕ, gp) ; dϕ .*= θ)

    θ₀ = Vector(gethyperparameters(gp))
    res = optimize(f, g!, log.(θ₀), solver, options)
    if Optim.converged(res)
        ϕ = Optim.minimizer(res)
        θ₀ .= exp.(ϕ)
        θ₀ != gethyperparameters(gp) ? sethyperparameters!(gp, θ₀) : nothing
        return true
    else
        sethyperparameters!(gp, θ₀)
        return false
    end
end

ext: jl


text:
using FastqSubsample
using Base.Test

ifastq = joinpath(dirname(@__FILE__), "input.fastq.gz")
ofastq = joinpath(dirname(@__FILE__), "output.fastq.gz")
ofastq2 = joinpath(dirname(@__FILE__), "output2.fastq.gz")
nSubsample = 40

fastqsubsample(ifastq, ofastq, nSubsample, seed = 1234)
fastqsubsample(ifastq, ofastq2, nSubsample, seed = 1234)

# Test for number of subsampled reads
@test FastqSubsample.countFastq(ofastq) == nSubsample
@test FastqSubsample.countFastq(ofastq2) == nSubsample

# Test for reproducibility
mdfive1 = split(chomp(readstring(`md5sum-lite $ofastq`)), r"\s+")[1]
mdfive2 = split(chomp(readstring(`md5sum-lite $ofastq2`)), r"\s+")[1]
@test mdfive1 == mdfive2

max_stars_repo_name: yuifu/FastqSubsample

max_issues_repo_name: yuifu/FastqSubsample

max_forks_repo_name: yuifu/FastqSubsample


text:
################################################# {BITS,ES} #### RSH HEX VALUE ########## NEG ## INV # EXP ## RSH FRAC VALUE #### UBIT ##
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,0}(0x4000)) == (false, false, 0, 0x0000_0000_0000_0000, false)
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x4000)) == (false, false, 0, 0x0000_0000_0000_0000, false)
#0b0_10_1_000000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x5000)) == (false, false, 1, 0x0000_0000_0000_0000, false)
#0b0_10_1_100000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x5800)) == (false, false, 1, 0x8000_0000_0000_0000, false)
#0b0_110_0_00000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x6000)) == (false, false, 2, 0x0000_0000_0000_0000, false)
#0b0_110_0_10000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x6400)) == (false, false, 2, 0x8000_0000_0000_0000, false)
#0b0_110_1_00000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x6800)) == (false, false, 3, 0x0000_0000_0000_0000, false)
#0b0_110_1_10000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x6C00)) == (false, false, 3, 0x8000_0000_0000_0000, false)
#0b0_1110_0_0000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x7000)) == (false, false, 4, 0x0000_0000_0000_0000, false)
#0b0_111111111111111
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x7FFF)) == (false, false, 28, 0x0000_0000_0000_0000, false)
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x7FFF)) == (false, false, 56, 0x0000_0000_0000_0000, false)
#0b0_111111111111110
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x7FFE)) == (false, false, 52, 0x0000_0000_0000_0000, false)
#0b0_111111111111101[0]
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x7FFD)) == (false, false, 50, 0x0000_0000_0000_0000, false)
#0b0_1111111111110_00
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x7FFD)) == (false, false, 50, 0x0000_0000_0000_0000, false)

#repeat some of these tests with valid numbers
#0b0_10_1_00000000000_1
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,1}(0x5001)) == (false, false, 1, 0x0000_0000_0000_0000, true)
#0b0_110_1_0000000000_1
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,1}(0x6801)) == (false, false, 3, 0x0000_0000_0000_0000, true)
#0b0_11111111111111_1
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,1}(0x7FFF)) == (false, false, 26, 0x0000_0000_0000_0000, true)

#repeat these tests, except with inverted numbers
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,0}(0x2000)) == (false, true, -1, 0x0000_0000_0000_0000, false)
#0b0_01_0000000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,0}(0x3000)) == (false, true, -1, 0x8000_0000_0000_0000, false)
#0b0_01_1000000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,0}(0x1000)) == (false, true, -2, 0x0000_0000_0000_0000, false)
#0b0_000000000000001
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,0}(0x0001)) == (false, true, -14, 0x0000_0000_0000_0000, false)
#0b0_01_0_00000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x2000)) == (false, true, -2, 0x0000_0000_0000_0000, false)
#0b0_01_1_00000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x3000)) == (false, true, -1, 0x0000_0000_0000_0000, false)
#0b0_01_00_0000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x2000)) == (false, true, -4, 0x0000_0000_0000_0000, false)
#0b0_01_01_0000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x2800)) == (false, true, -3, 0x0000_0000_0000_0000, false)
#0b0_01_10_0000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x3000)) == (false, true, -2, 0x0000_0000_0000_0000, false)
#0b0_01_11_0000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x3800)) == (false, true, -1, 0x0000_0000_0000_0000, false)
#0b0_0000000000001_0
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x0002)) == (false, true, -52, 0x0000_0000_0000_0000, false)
#0b0_0000000000001_1
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x0003)) == (false, true, -50, 0x0000_0000_0000_0000, false)

#and now with valids
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,0}(0x2000)) == (false, true, -1, 0x0000_0000_0000_0000, false)
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,0}(0x2001)) == (false, true, -1, 0x0000_0000_0000_0000, true)
#and then go really extreme.
#0b0_0000000000001_0
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,1}(0x0002)) == (false, true, -26, 0x0000_0000_0000_0000, false)
#0b0_0000000000001_1
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,1}(0x0003)) == (false, true, -26, 0x0000_0000_0000_0000, true)
#semantically different due to the a higher ES
#0b0_0000000000001_0
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,2}(0x0002)) == (false, true, -52, 0x0000_0000_0000_0000, false)
#0b0_0000000000001_1
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,2}(0x0003)) == (false, true, -52, 0x0000_0000_0000_0000, true)

## all these tests, but with negative numbers.
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,0}(0xC000)) == (true, false, 0, 0x0000_0000_0000_0000, false)
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0xC000)) == (true, false, 0, 0x0000_0000_0000_0000, false)
#0b0_10_1_000000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0xB000)) == (true, false, 1, 0x0000_0000_0000_0000, false)
#0b0_10_1_100000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0xA800)) == (true, false, 1, 0x8000_0000_0000_0000, false)
#0b0_110_0_00000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0xA000)) == (true, false, 2, 0x0000_0000_0000_0000, false)
#0b0_110_0_10000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x9C00)) == (true, false, 2, 0x8000_0000_0000_0000, false)
#0b0_110_1_00000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x9800)) == (true, false, 3, 0x0000_0000_0000_0000, false)
#0b0_110_1_10000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x9400)) == (true, false, 3, 0x8000_0000_0000_0000, false)
#0b0_1110_0_0000000000
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x9000)) == (true, false, 4, 0x0000_0000_0000_0000, false)
#0b0_111111111111111
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,1}(0x8001)) == (true, false, 28, 0x0000_0000_0000_0000, false)
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x8001)) == (true, false, 56, 0x0000_0000_0000_0000, false)
#0b0_111111111111110
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x8002)) == (true, false, 52, 0x0000_0000_0000_0000, false)
#0b0_111111111111101[0]
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x8003)) == (true, false, 50, 0x0000_0000_0000_0000, false)
#0b0_1111111111110_00
@test SigmoidNumbers.__numeric_breakdown_wrapper(Posit{16,2}(0x8003)) == (true, false, 50, 0x0000_0000_0000_0000, false)

#repeat some of these tests with valid numbers
#0b0_10_1_00000000000_1
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,1}(0xAFFF)) == (true, false, 1, 0x0000_0000_0000_0000, true)
#0b0_110_1_0000000000_1
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,1}(0x97FF)) == (true, false, 3, 0x0000_0000_0000_0000, true)
#0b0_11111111111111_1
@test SigmoidNumbers.__numeric_breakdown_wrapper(Vnum{16,1}(0x8001)) == (true, false, 26, 0x0000_0000_0000_0000, true)

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
using Images, TestImages, ImageFiltering, MosaicViews, BenchmarkTools
using FFTW

function my_diffy(img)
    d = similar(img)
    d[:, 1:end-1, :] .= @views img[:, 2:end, :] .- img[:, 1:end-1, :]
    d[:, end:end, :] .= @views img[:, 1:1, :] .- img[:, end:end, :]
    d
end

function my_diffyt(img)
    d = similar(img)
    d[:, 2:end, :] .= @views img[:, 1:end-1, :] .- img[:, 2:end, :]
    d[:, 1:1, :] .= @views img[:, end:end, :] .- img[:, 1:1, :]
    d
end

function my_diffx(img)
    d = similar(img)
    d[:, :, 1:end-1] .= @views img[:, :, 2:end] .- img[:, :, 1:end-1]
    d[:, :, end:end] .= @views img[:, :, 1:1] .- img[:, :, end:end]
    d
end

function my_diffxt(img)
    d = similar(img)
    d[:, :, 2:end] .= @views img[:, :, 1:end-1] .- img[:, :, 2:end]
    d[:, :, 1:1] .= @views img[:, :, end:end] .- img[:, :, 1:1]
    d
end

expanded_channelview(img::AbstractArray{T}) where T<:Colorant = channelview(img)
expanded_channelview(img::AbstractArray{T}) where T<:Gray = reshape(channelview(img), 1, size(img)...)

function l0smoothing(img,  λ=2e-2, κ=2.0)
    S = float.(expanded_channelview(img))
    βmax = 1e5
    fx = [1 -1]
    fy = [1, -1]
    D, N, M = size(S)
    sizeI2D = (N, M)
    sizeI2D_t = (M, N)
    otfFx = freqkernel(centered(fx), sizeI2D)
    otfFy = transpose(freqkernel(centered(transpose(fy)), sizeI2D_t))
    Normin1 = fft(S, (2, 3))
    Denormin2 = abs.(otfFx).^2 + abs.(otfFy ).^2
    if D > 1
        Denormin2 = repeat(reshape(Denormin2, 1, size(Denormin2)...), inner=(1, 1, 1), outer=(3, 1, 1))
    else
        Denormin2 = reshape(Denormin2, 1, size(Denormin2)...)
    end
    β = 2*λ
    while β < βmax
        Denormin = 1 .+ β*Denormin2

        h = my_diffx(S)
        v = my_diffy(S)

        if D == 1
            t = (h.^2+v.^2) .< λ/β
        else
            t = sum((h.^2+v.^2), dims=1) .< λ/β
            t = repeat(t, inner=(1, 1, 1), outer=(D, 1, 1))
        end

        h[t] .= 0
        v[t] .= 0

        Normin2 = my_diffxt(h)
        Normin2 = Normin2 + my_diffyt(v)
        FS = (Normin1 + β*fft(Normin2, (2, 3))) ./ Denormin
        S = real(ifft(FS, (2, 3)))
        β = β*κ
    end
    if D == 1
        return colorview(Gray, S[1,:,:])
    else
        return colorview(RGB, S)
    end
end

img1 = load("pflower.jpg")
img2 = testimage("cameraman")

S1 = @btime l0smoothing(img1)
S2 = @btime l0smoothing(img2)

mosaicview(img1,img2, S1, S2;nrow=2, ncol=2, rowmajor=true)


text:
module UnstructuredGridTopologiesTests

using Test
using Gridap.Helpers
using Gridap.Arrays
using Gridap.Fields
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Geometry: GridTopologyMock

m = GridTopologyMock()

g = UnstructuredGridTopology(
  get_vertex_coordinates(m),
  get_cell_vertices(m),
  get_cell_type(m),
  get_polytopes(m))

test_grid_topology(g)
test_grid_topology(g)

@test num_faces(g,0) == num_faces(m,0)
@test num_faces(g,1) == num_faces(m,1)
@test num_faces(g,2) == num_faces(m,2)
@test get_reffaces(Polytope{1},g) == get_reffaces(Polytope{1},m)
@test get_reffaces(Polytope{0},g) == get_reffaces(Polytope{0},m)
@test get_reffaces(Polytope{2},g) == get_reffaces(Polytope{2},m)
@test get_face_type(g,1) == get_face_type(m,1)
@test get_isboundary_face(g,0) == get_isboundary_face(m,0)
@test get_isboundary_face(g,1) == get_isboundary_face(m,1)
@test get_isboundary_face(g,2) == get_isboundary_face(m,2)
@test get_faces(g,2,0) == get_faces(m,2,0)
@test get_faces(g,2,1) == get_faces(m,2,1)
@test get_faces(g,1,2) == get_faces(m,1,2)
@test get_faces(g,2,1) == get_faces(m,2,1)
@test get_faces(g,0,1) == get_faces(m,0,1)
@test get_faces(g,1,0) == get_faces(m,1,0)
@test get_faces(g,1,1) == get_faces(m,1,1)
@test get_vertex_coordinates(g) == get_vertex_coordinates(m)
@test is_oriented(g) == false

g = UnstructuredGridTopology(
  get_vertex_coordinates(m),
  get_cell_vertices(m),
  get_cell_type(m),
  get_polytopes(m),
  Val{true}())

@test is_oriented(g) == true

g = UnstructuredGridTopology(
  get_vertex_coordinates(m),
  [get_face_vertices(m,d) for d in 0:num_cell_dims(m)],
  get_cell_type(m),
  get_polytopes(m))

@test get_faces(g,2,0) == get_faces(m,2,0)
@test get_faces(g,2,1) == get_faces(m,2,1)
@test get_faces(g,1,2) == get_faces(m,1,2)
@test get_faces(g,2,1) == get_faces(m,2,1)
@test get_faces(g,0,1) == get_faces(m,0,1)
@test get_faces(g,1,0) == get_faces(m,1,0)
@test get_faces(g,1,1) == get_faces(m,1,1)
@test is_oriented(g) == false

test_grid_topology(g)
test_grid_topology(g)

g = UnstructuredGridTopology(
  get_vertex_coordinates(m),
  [get_face_vertices(m,d) for d in 0:num_cell_dims(m)],
  get_cell_type(m),
  get_polytopes(m),
  Val{true}())

@test is_oriented(g) == true

m = GridTopologyMock()
topo = UnstructuredGridTopology(m)
test_grid_topology(topo)
@test topo === UnstructuredGridTopology(topo)

end # module


text:
export EmptyAccelerator
# ---------------------------
# EmptyAccelerator
# ---------------------------

"""
    EmptyAccelerator <: AbstractAccelerator

An accelerator that does nothing.
"""
struct EmptyAccelerator <: AbstractAccelerator
  function EmptyAccelerator()
      return new()
  end
end
EmptyAccelerator(args...; kwargs...) = EmptyAccelerator()

function update!(ea::EmptyAccelerator, args...)
  return nothing
end

function accelerate!(g::AbstractVector{T}, x::AbstractVector{T}, ea::EmptyAccelerator, args... ) where {T <: AbstractFloat}
  return false
end
restart!(ea::EmptyAccelerator) = nothing
log!(ea, args...; kwargs...) = nothing
was_successful(ea::EmptyAccelerator) = false

max_issues_repo_name: oxfordcontrol/COSMOAccelerators

max_forks_repo_name: oxfordcontrol/COSMOAccelerators


text:
using Unitful
using DFTK

lattice = load_lattice(joinpath(@__DIR__, "WFe.json"))
atoms   = load_atoms(  joinpath(@__DIR__, "WFe.json"))
atoms   = attach_pseudos(atoms, Fe="hgh/pbe/fe-q16.hgh", W="hgh/pbe/w-q14.hgh")
magnetic_moments = build_magnetic_moments(atoms; Fe=3.0, W=2.0)

cases = standard_cases()
cases["standard_06"] = (damping=DFTK.FixedDamping(0.6), tol=1e-12, )

config[@__DIR__] = (
    directory        =  @__DIR__,
    lattice          = lattice,
    atoms            = atoms,
    magnetic_moments = magnetic_moments,
    #
    kwargs_model     = (smearing=Smearing.Gaussian(), temperature=0.01),
    Ecut             = 25,
    kgrid            = (7, 7, 1),
    supersampling    = 1.5,
    potential        = nothing,
    kwargs_scf       = (; ),
    mixing           = KerkerMixing(),
    cases            = cases,
)

max_stars_repo_name: mfherbst/supporting-adaptive-damping

max_issues_repo_name: mfherbst/supporting-adaptive-damping

max_forks_repo_name: mfherbst/supporting-adaptive-damping


text:
using Documenter, IterableTables

makedocs(
	modules = [IterableTables],
	format = :html,
	sitename = "IterableTables.jl",
	pages = [
		"Introduction" => "index.md",
        "User Guide" => "userguide.md",
        "Integration Guide" => "integrationguide.md",
        "Developer Guide" => "developerguide.md"]
)

deploydocs(
    deps = nothing,
    make = nothing,
    target = "build",
    repo = "github.com/queryverse/IterableTables.jl.git",
    julia = "0.7"
)


text:
using DeepDiffs
using Test

if isdefined(Base, :have_color)
    # Capture the original state of the global flag
    orig_color = Base.have_color
end

@testset "DeepDiff Tests" begin
    include("arrays.jl")
    include("dicts.jl")
    include("display.jl")
    include("simplediff.jl")
    include("strings.jl")
end

max_stars_repo_name: ssfrr/StructureDiffs

max_issues_repo_name: ssfrr/StructureDiffs

max_forks_repo_name: ssfrr/StructureDiffs


text:
using SequentialFit, LsqFit, Plots
default(; label=nothing)

function benchmark(; noise_level=0.01, k_0=5, K=200, lims=(-5, 5))
    p_true = randn(3 * 1)#randn(3 * rand(1:10))
    p_0 = p_true .+ 0.1 * randn(size(p_true))
    x_0 = collect(range(lims...; length=k_0))

    f(x) = model(x, p_true) .+ noise_level * randn(size(x))

    s = SeqFit(f, model, x_0, p_0, lims)

    p = Array{Float64}(undef, length(p_true), 2)
    error = Array{Float64}(undef, K - k_0, 2)
    for k in (k_0 + 1):K
        update!(s)
        error[k - k_0, 1] = errorcalc(s.fit.param, p_true) #sqrt(sum(abs2, s.fit.param .- p_true))
        x = collect(range(lims...; length=k))
        fit = curve_fit(model, x, f.(x), p_0)
        error[k - k_0, 2] = errorcalc(fit.param, p_true) #sqrt(sum(abs2, fit.param .- p_true))
        p .= [s.fit.param fit.param]
    end

    x = range(lims...; length=100)
    pl1 = plot((k_0 + 1):K, error; color=[1 2], label=["Sequential fit" "Linear space"], yscale=:log10, xguide="k", yguide="error")

    pl2 = plot(x, model(x, p_true); color=:black, label="Ground truth", xguide="x", yguide="f")
    plot!(pl2, x, f(x); color=:gray, label="Noisy signal")
    plot!(pl2, s.X, s.Y; color=1, seriestype=:scatter)
    plot!(pl2, x, model(x, p[:, 1]); color=1, label="Sequential fit")
    plot!(pl2, x, model(x, p[:, 2]); color=2, label="Linear sample fit")
    plot!(pl2; ylims=extrema(model(x, p_true)), yexpand=true)
    return plot(pl1, pl2; layout=(2, 1))
end

function sumofgaussians(x, p)
    return sum(
               exp.(-(x .- m) .^ 2 ./ (0.2w .^ 2)) * a for
               (a, m, w) in zip(p[1:3:end], p[2:3:end], p[3:3:end])
              )
end

function model(x, p)
    return sumofgaussians.(x, Ref(p))
end

function errorcalc(p, p_0)
    (a_0, m_0, w_0) = (p_0[1:3:end], p_0[2:3:end], p_0[3:3:end])
    (a, m, w) = (p[1:3:end], p[2:3:end], p[3:3:end])
    i_0 = sortperm(m_0)
    i = sortperm(m)
    return sqrt(
                sum(abs2, vcat(a[i] .- a_0[i_0], m[i] .- m_0[i_0], w[i] .- w_0[i_0]))
               )
end


text:
mutable struct DataXMLElement <: MyXMLElement
    el::XMLOrNothing
    data_mats::Vector{<:AbstractMatrix{Float64}}
    trait_names::Vector{S} where S <: AbstractString
    taxa::Vector{<:AbstractString}
    node_times::Vector{Float64}
    use_times::Bool

    function DataXMLElement(mat::AbstractMatrix{Float64},
                taxa::Vector{S};
                trait_name::AbstractString = bn.DEFAULT_TRAIT_NAME
                ) where S <: AbstractString

        node_times = Float64[]
        return new(nothing, [mat], [trait_name], taxa, node_times, false)
    end

    function DataXMLElement(mats::Vector{<:AbstractMatrix{Float64}},
                trait_names::Vector{T} where T <: AbstractString,
                taxa::Vector{S}
                ) where S <: AbstractString

        node_times = Float64[]

        return new(nothing, mats, trait_names, taxa, node_times, false)
    end

end

function DataXMLElement(dm::DataModel)
    n_partitions = length(dm.partitions)
    n = size(dm.data, 1)
    data = [@view dm.data[1:n, dm.partitions[i]] for i = 1:n_partitions]
    return DataXMLElement(data, dm.trait_names, dm.taxa)
end

function make_xml(dl::DataXMLElement)
    dl.el = make_taxa_el(dl)
end

function make_taxa_el(dxl::DataXMLElement)

    data_mats = dxl.data_mats
    el = new_element(bn.TAXA)
    n, p = size(data_mats[1])
    for i = 1:length(data_mats)
        @assert size(data_mats[i], 1) == n
    end
    for i = 1:n
        add_taxon(dxl, el, i)
    end
    return el
end

function add_taxon(dxl::DataXMLElement, pel::XMLElement, ind::Int)

    taxa = dxl.taxa
    data_mats = dxl.data_mats
    trait_names = dxl.trait_names
    node_times = dxl.node_times

    el = new_element(bn.TAXON)
    set_attribute(el, bn.ID, taxa[ind])
    if dxl.use_times
        date_el = new_child(el, bn.DATE)
        set_attributes(date_el, [(bn.VALUE, string(node_times[ind])), (bn.DIRECTION, bn.FORWARDS)])
    end
    for i = 1:length(data_mats)
        attr_el = new_element(bn.ATTR)
        set_attribute(attr_el, bn.NAME, trait_names[i])
        s_data = format_data_line(data_mats[i], ind)
        set_content(attr_el, s_data)
        add_child(el, attr_el)
    end
    add_child(pel, el)
    return el
end

function format_data_line(data::AbstractMatrix{Float64}, ind::Int)
    p = size(data, 2)
    s_data = Vector{String}(undef, p)
    for j = 1:p
        val = data[ind, j]
        if isnan(val)
            s_data[j] = bn.MISSING_VAL
        else
            s_data[j] = string(val)
        end
    end
    return join(s_data, ' ')
end

# function get_node_times(taxa::Vector{String}, newick::String)
#     tree = PhyloNetworks.readTopology(newick)
#     root_dists = TreeUtils.leaf_distances(tree)
#     tree_taxa = [tree.leaf[i].name for i = 1:tree.numTaxa]
#     node_times = zeros(tree.numTaxa)
#     for i = 1:length(taxa)
#         taxon = taxa[i]
#         ind = findfirst(x -> x == taxon, tree_taxa)
#         node_times[i] = root_dists[ind]
#     end
#     return node_times
# end

function add_trait!(dxl::DataXMLElement, data::Matrix{Float64}, trait_name::String)
    push!(dxl.data_mats, data)
    push!(dxl.trait_names, trait_name)
end

function set_dates(xml::DataXMLElement, dates::AbstractVector{Float64})
    if length(xml.node_times) == 0
        xml.node_times = copy(dates)
    elseif length(dates) != length(xml.node_times)
        throw(DimensionMismatch("The number of dates must match the number of taxa."))
    else
        xml.node_times .= dates
    end
end


text:
################################################################################
#
#            GrpAb/Lattice.jl : Lattice of abelian groups
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2017 Tommy Hofmann
#
################################################################################

################################################################################
#
#  Graph with data at the edges
#
################################################################################

# A graph G consists of nested dictionaries to store edges with data attached.
# Additionally, we keep track of the degree of each vertex in the G.degrees
# field. Whenever we modify G by appending or deleting vertices, we check
# wether the degree of a vertex becomes < 2. If so, we store the vertex in
# G.new_low_degrees.

function Base.show(io::IO, G::Graph{T, M}) where {T, M}
  print("Graph with $(length(G.edges)) vertices of type $T and edge data of type $M\n")
end

# Append an element to a graph
function Base.append!(G::Graph{T, M}, a::T) where {T, M}
  if haskey(G.edges, a)
    error("Element already a vertex of the graph")
  end
  G.edges[a] = Dict{T, M}()
  G.degrees[a] = 0
  G.new_low_degrees[a] = nothing
end

# For debugging purposes
function _compute_degrees(G::Hecke.Graph{T, N}) where {T, N}
  l = Dict{T, Int}()
  for e in keys(G.edges)
    l[e] = length(G.edges[e])
    for a in keys(G.edges)
      if haskey(G.edges[a], e)
        l[e] += 1
      end
    end
  end
  return l
end

# Append an edge to a graph
function Base.append!(G::Graph{T, M}, e::Tuple{T, T}, data::M) where {T, M}
  if !haskey(G.edges, e[1])
    append!(G, e[1])
  end
  if !haskey(G.edges, e[2])
    append!(G, e[2])
  end
  G.edges[e[1]][e[2]] = data
  G.degrees[e[1]] += 1
  G.degrees[e[2]] += 1

  if G.degrees[e[1]] > 1
    delete!(G.new_low_degrees, e[1])
  else
    G.new_low_degrees[e[1]] = nothing
  end

  if G.degrees[e[2]] > 1
    delete!(G.new_low_degrees, e[2])
  else
    G.new_low_degrees[e[2]] = nothing
  end
end

# Delete a vertex from a graph
function Base.delete!(G::Graph{T, M}, a::T) where {T, M}
  Base.delete!(G.new_low_degrees, a)
  Base.delete!(G.degrees, a)

  if haskey(G.edges, a)
    for e in keys(G.edges[a])
      G.degrees[e] -= 1
      if G.degrees[e] < 2
        G.new_low_degrees[e] = nothing
      end
    end
    Base.delete!(G.edges, a)
  end

  #@show "iterating over $(length(keys(G.edges)))"
  for e in keys(G.edges)
    if haskey(G.edges[e], a)
      G.degrees[e] -= 1
      Base.delete!(G.edges[e], a)
      if G.degrees[e] < 2
        G.new_low_degrees[e] = nothing
      end
    end
  end
end

function Base.deepcopy_internal(G::Graph{T, M}, dict::IdDict) where {T, M}
  GG = Graph{T, M}()
  for g in keys(G.edges)
    GG.edges[g] = Base.deepcopy_internal(G.edges[g], dict)
  end
  GG.degrees = Base.deepcopy_internal(G.degrees, dict)
  GG.new_low_degrees = Base.deepcopy_internal(G.new_low_degrees, dict)
  return GG
end

function Base.delete!(G::Graph{T, M}, to_delete::Vector{T}) where {T, M}
  if length(to_delete) == 0
    return
  end
  #@show "have to delete $(length(to_delete)) many things"
  #@show "there will be $(length(G.edges) - length(to_delete)) left"
  if length(G.edges) - length(to_delete) < 1000
    #@show "Going through the fast path"
    #GG = deepcopy(G)
    _delete_fastpath(G, to_delete)
  else
    while length(to_delete) > 0
      a = pop!(to_delete)
      Base.delete!(G, a)
    end
  end
end

function _delete_fastpath(G::Graph{T, M}, to_delete::Vector{T}) where {T, M}
  _to_delete = T[]

  while length(to_delete) > 0
    a = pop!(to_delete)
    Base.delete!(G.new_low_degrees, a)
    Base.delete!(G.degrees, a)
    Base.delete!(G.edges, a)
    push!(_to_delete, a)
  end

  for e in keys(G.edges)
    for a in _to_delete
      delete!(G.edges[e], a)
    end
  end
  # Now compute the degrees new
  for e in keys(G.edges)
    G.degrees[e] = length(G.edges[e])
    for a in keys(G.edges)
      if haskey(G.edges[a], e)
        G.degrees[e] += 1
      end
    end
    if G.degrees[e] < 2
      G.new_low_degrees[e] = nothing
    end
  end
end

# Find the shortest path between root and target (if it exists).
# The first return value indicates on whether there exists a path.
# The second return value is an array of type Vector{T},
# with T[1] == target and T[end] == root (note that it is reversed).
function find_shortest(G::Graph{T, M}, root::T, target::T) where {T, M}
  if !(root in keys(G.edges) && target in keys(G.edges))
    error("Elements are not vertices of the graph")
  end
  S = T[]
  Swithparent = Tuple{T, T}[]
  Q = T[]
  push!(S, root)
  push!(Q, root)

  found = false

  while length(Q) > 0 && found == false
    current = pop!(Q)
    if current == target
      found = true
      break
    end
    for n in keys(G.edges[current])
      if !(n in S)
        push!(S, n)
        push!(Swithparent, (current, n))
        if n == target
          found = true
          break
        end
        prepend!(Q, n)
      end
    end
  end

  if !found
    return false, T[]
  end

  current = target
  path = [ target ]

  while current != root
    for s in Swithparent
      if s[2] == current
        current = s[1]
        push!(path, current)
        break
      end
    end
  end

  return true, path
end

# Given two vertices root1 and root2, check if the "connected components"
# intersect. If so, find two paths from root1 and root2 which end in a vertex
# of the intersect. The second and third return values are arrays T1
# and T2 with T1[1] == T2[1], T1[end] == root1 and T2[end] == root2.
function find_common(G::Graph{T, M}, root1::T, root2::T) where {T, M}
  S1 = T[]
  S1withparent = Tuple{T, T}[]
  Q1 = T[]
  push!(S1, root1)
  push!(Q1, root1)

  S2 = T[]
  S2withparent = Tuple{T, T}[]
  Q2 = T[]
  push!(S2, root2)
  push!(Q2, root2)

  found = false

  current = root1

  while length(Q1) > 0 && found == false
    current = pop!(Q1)

    if current in S2
      found = true
    end

    for n in keys(G.edges[current])
      if !(n in S1)
        push!(S1, n)
        push!(S1withparent, (current, n))
        prepend!(Q1, n)
      end
    end
  end

  while length(Q2) > 0  && found == false
    current = pop!(Q2)
    if current in S1
      found = true
      break
    end
    for n in keys(G.edges[current])
      if !(n in S2)
        push!(S2, n)
        push!(S2withparent, (current, n))
        prepend!(Q2, n)
      end
    end
  end

  if found
    target = current

    path1 = [ target ]

    while current != root1
      for s in S1withparent
        if s[2] == current
          current = s[1]
          push!(path1, current)
          break
        end
      end
    end

    path2 = [ target ]

    current = target

    while current != root2
      for s in S2withparent
        if s[2] == current
          current = s[1]
          push!(path2, current)
          break
        end
      end
    end
    return true, path1, path2
  end

  false, T[], T[]
end

################################################################################
#
#  Lattice of finitely generated abelian groups
#
################################################################################

# A lattice of groups L has basically two parts. The first is a weak dictionary
# L.weak_vertices, which hold references to all groups, which are currently
# in the lattice. The second part is an actual graph, whose vertices are the
# objectid's of the groups and whose edges represented maps (the additional
#fmpz_mat(0, 0) data at an edge is the fmpz_mat describing the map).
#
# Now things get complicated (interesting) due to the presence of the gc.
# Here is the most important rule:
#
#   A vertex in the lattice of groups is allowed to be removed only if
#   the degree is < 2.
#
# We achieve this by keeping an additional dictionary L.block_gc, which
# contains the groups which must not be gc'ed.
#
# More details:
# - L.weak_vertices_rev is the "inverse" of L.weak_vertices. For each
#   objectid it contains the corresponding group as a WeakRef.
# - Whenever a group (or a map) is added to a lattice, we check if the
#   degree of a vertex is > 1. If so, we add to L.block_gc.
# - Whenever a group is added to a lattice, a finalizer is added. If then
#   the group is gc'ed, its objectid is added to L.to_delete.
# - Before a group (or a map) is added to a lattice, we update the lattice:
#   - Remove all vertices from the graph whose vertices are in L.to_delete
#   - If the degree of a vertex drops below 2, remove the group from L.block_gc
#
# Minor but important note: Over time, it may happen that two groups (created
# at different times) have the same objectid. But this can happen only, if
# one of them was already gc'ed (so its objectid is contained in L.to_delete).
# Thus, it is important that we update the lattice before we add a new group.
# This way the correspondence G -> objectid(G) is always faithful.

# T is the type of the objects i the lattice
# D is the type of the data attached to the edges, usually some matrix
#  objects of type D need to support products
#  if would be useful if hom(::T, ::T, ::D) would exist

function Base.show(io::IO, L::RelLattice{T, D}) where {T, D}
  print("Relation lattice for $T with underlying graph \n$(L.graph)\n")
  print("In weak dict: $(length(L.weak_vertices))\n")
  print("In dict: $(length(L.block_gc))\n")
end

# The finalizer, which must be attached to a every group in the lattice.
# Note that it may happen to a lot of groups are gc'ed so that the
# L.to_delete queue get's too big. We cut it off at 100.
function finalizer_lattice(L::RelLattice{T, D}, A::T) where {T, D}
  # TODO: Find out why the following does not work
  #if length(L.to_delete) > 1000
  #  Hecke.update!(L)
  #end
  push!(L.to_delete, objectid(A))
end

# Add an object to a lattice of groups.
function Base.append!(L::RelLattice{T, D}, A::T) where {T, D}
  update!(L)
  L.weak_vertices[A] = nothing
  obid = objectid(A)
  L.weak_vertices_rev[obid] = WeakRef(A)
  append!(L.graph, obid)

  if L.graph.degrees[obid] > 1
    L.block_gc[A] = nothing
  end

  if !A.isfinalized
    finalizer(x -> finalizer_lattice(L, x), A)
    A.isfinalized = true
  end
end

# Add a map to a lattice of groups
# .. sugar for ablian groups ...
function Base.append!(L::GrpAbLattice, f::Map)
  return Base.append!(L, domain(f), codomain(f), f.map)
end

function Base.append!(L::RelLattice{T, D}, dom::T, co::T, f::D) where {T, D}
  update!(L)

  if !haskey(L.weak_vertices, dom)
    append!(L, dom)
  end

  if !haskey(L.weak_vertices, co)
    append!(L, co)
  end

  append!(L.graph, (objectid(dom), objectid(co)), f)

  if L.graph.degrees[objectid(dom)] > 1
    L.block_gc[dom] = nothing
  end

  if L.graph.degrees[objectid(co)] > 1
    L.block_gc[co] = nothing
  end
end

# Delete the group with objectid u from a lattice of groups.
function delete_from_lattice!(L::RelLattice, u::UInt)
  Base.delete!(L.graph, u)
  Base.delete!(L.weak_vertices_rev, u)
end

function delete_from_lattice!(L::RelLattice, us::Vector{UInt})
  Base.delete!(L.graph, us)
  for u in us
    Base.delete!(L.weak_vertices_rev, u)
  end
end

# Update a lattice of groups.
function update!(L::RelLattice)
  #while length(L.to_delete) > 0
  #  u = pop!(L.to_delete)
  #  delete_from_lattice!(L, u)
  #end
  delete_from_lattice!(L, L.to_delete)
  #L.to_delete = UInt[]

  for k in keys(L.graph.new_low_degrees)
    @assert haskey(L.graph.degrees, k)
    # TODO: Why does it crash without the following?
    if !haskey(L.weak_vertices_rev, k)
      delete!(L.graph.new_low_degrees, k)
      continue
    end
    @assert L.weak_vertices_rev[k].value != nothing
    a = L.weak_vertices_rev[k].value
    @assert k == objectid(a)
    @assert L.graph.degrees[k] < 2
    L.weak_vertices[a] = nothing
    Base.delete!(L.block_gc, a)
    delete!(L.graph.new_low_degrees, k)
  end

  return nothing
end

# For a given lattice L and groups G and H, check whether there exists a path
# from G to H. The second return value is a fmpz_mat describing the map (in
# case it exists).

function eval_path(L::RelLattice{T, D}, M::T, pG::Vector{UInt}) where {T, D}
  lG = length(pG)
  if lG == 1
    return L.make_id(M)::D
  end
  mG = L.graph.edges[pG[lG]][pG[lG - 1]]
  for i in lG-1:-1:2
    mG = L.mult(mG, (L.graph.edges[pG[i]][pG[i - 1]]))::D
  end
  return mG
end
# Todo: Reduce the entries of the final matrix.
function can_map_into(L::RelLattice{T, D}, G::T, H::T) where {T, D}
  if !(G in keys(L.weak_vertices) && H in keys(L.weak_vertices))
    return false, L.zero
  end

  b, p = find_shortest(L.graph, objectid(G), objectid(H))
  if b
    m = eval_path(L, G, p)
    return true, m
  else
    return false, L.zero
  end
end

# For a given lattice L and groups G and H, check whether there exists a common
# "overgroup" M. If so, the second return value is M and the third and fourth
# return values describe the map from G to M and H to M respectively.
function can_map_into_overstructure(L::RelLattice{T, D}, G::T, H::T) where {T, D}
  if !(G in keys(L.weak_vertices) && H in keys(L.weak_vertices))
    return false, G, L.zero, L.zero
  end
  b, pG, pH = find_common(L.graph, objectid(G), objectid(H))
  if b
    @assert pG[1] == pH[1]
    M = L.weak_vertices_rev[pG[1]].value::T
    @assert M != nothing

    mG = eval_path(L, M, pG)
    mH = eval_path(L, M, pH)

    return true, M, mG, mH
  else
    return false, G, L.zero, L.zero
  end
end


text:
# Implement experiments to measure the success of subgoal evolution as a function of:
#  1.  Funcs
#  2.  numinteriors
#  3.  numinputs #  4.  numoutputs
#  5.  Length of goallist
#  5.  max_steps
# Uses randomly generated goallist
# Keep track of number of "worse" and "same" updates of Chromosome since its last improvement
export run_mut_evolution, run_result, run_result_to_tuple
using DataFrames
using CSV
using Statistics
using Distributions
using Distributed
indiv_result_type = Main.CGP.indiv_result_type
#=
iterations = 4
numinputs = 2:2
numoutputs = 2:2
nodearity = 2
numinteriors = 18:18
numlevelsback = 6:6
ngoals = 4:4
goallistlength=8:8
levelsback=6:6
maxsteps = 400:400
run_mut_evolution( iterations, numinputs, numoutputs, numinteriors, goallistlength, maxsteps, levelsback, "testdata.csv" )
=#
#context = construct_contexts(numinputs)[numinputs]
#p = Parameters(numinputs, numoutputs, nodearity, numinteriors, numlevelsback)
#iterations = 20
#print_parameters(p)

function run_mut_evolution( numiterations::Int64, numinputs::IntRange, numoutputs::IntRange,
    numinteriors::IntRange, goallistlength::IntRange, maxsteps::IntRange,
    levelsback::IntRange, hamming_rng::IntRange, csvfile::String;
    base::Float64=2.0, active_only::Bool=false, avgfit_rng::IntRange=false )
  maxints_for_degen=20
  max_numinteriors = collect(numinteriors)[end]
  test_MyInt(max_numinteriors)
  ftf_param = 0.95
  nodearity = 2
  run_result_list = indiv_result_type[]
  df = DataFrame()
  df.numinputs=Int64[]
  df.numoutputs=Int64[]
  df.numints=Int64[]
  df.levelsback=Int64[]
  df.ngoals=Int64[]
  #df.hamming_sel=Bool[]
  #df.robust_sel=Bool[]
  #df.avgfitness=Bool[]
  #df.fault_tol=Bool[]
  #df.active_only=Bool[]
  df.maxsteps=Int64[]
  #df.gl_reps=Int64[]
  #df.fit_limit=Float64[]
  df.hdist=Float64[]
  df.steps=Int64[]
  df.same=Int64[]
  df.worse=Int64[]
  df.better=Int64[]
  df.nactive=Int64[]
  df.redundancy=Float64[]
  df.complexity=Float64[]
  df.degeneracy=Float64[]
  df.sdegeneracy=Float64[]
  println("size(df): ",size(df))
  for num_inputs = numinputs
    for num_outputs = numoutputs
      fit_limit = Float64(num_outputs)
      funcs = default_funcs(num_inputs)
      for num_interiors = numinteriors
        for num_goals = goallistlength
          println("numinputs: ",num_inputs,"  numoutputs: ",num_outputs,"  numints: ",num_interiors,"  numgoals: ",num_goals)
          for levsback = levelsback
            for hamming_sel = hamming_rng
              for avgfitness = avgfit_rng
                for active_only = [false]
                  #println("hamming_sel: ",hamming_sel,"  active_only: ",active_only)
                  for max_steps = maxsteps
                    for _ = 1:numiterations
                      p = Parameters( num_inputs, num_outputs, nodearity, num_interiors, levsback )
                      rr = run_result( p, num_goals, hamming_sel, active_only, max_steps )
                      push!(run_result_list,rr)
                      #Base.push!( df, new_row )
                    end
                  end
                end
              end
            end
          end
        end
      end
    end
  end
  new_run_result_list = pmap(r->run_mut_evolve!(r,maxints_for_degen=maxints_for_degen,base=base),run_result_list)
  #new_run_result_list = map(r->run_mut_evolve!(r,maxints_for_degen=maxints_for_degen,base=base),run_result_list)
  for r = new_run_result_list
    new_row = run_result_to_tuple(r)
    println("len row: ",length(new_row))
    Base.push!( df, new_row )
  end
  println(default_funcs(2))
  open( csvfile, "w" ) do f
    println(f,"# funcs: ", Main.CGP.default_funcs(numinputs[end]))
    println(f,"# nodearity: ",nodearity)
    #println(f,"# active_only: ",active_only)
    println(f,"# max_steps",maxsteps)
    CSV.write( f, df, append=true, writeheader=true )
  end
  #println(df)
  df
end

function run_mut_evolve!( rr::indiv_result_type; maxints_for_degen::Int64, gl_repetitions::Int64=1,
      ftf_param::Float64=0.95, base::Float64=2.0 )
  #println("run_mut_evolve! fault_tol: ",rr.fault_tol)
  nodearity = 2   # built-in default
  p = Parameters( numinputs=rr.numinputs, numoutputs=rr.numoutputs, numinteriors=rr.numints, numlevelsback=rr.levelsback )
  #print_parameters( p )
  gl = randgoallist(rr.ngoals,rr.numinputs,rr.numoutputs,repetitions=gl_repetitions)
  #println("gl: ",gl)
  funcs = default_funcs(rr.numinputs)
  c = random_chromosome( p, funcs )
  sav_c = deepcopy(c)
  (new_c,rr.steps,rr.worse,rr.same,rr.better,output,matched_goals,matched_goals_list) =
      mut_evolve(c,gl,funcs,rr.maxsteps,hamming_sel=rr.hamming_sel )
  rr.hdist = hamming_distance( output_values(c), output_values(new_c), rr.numinputs )
  rr.nactive = number_active( new_c )
  rr.redundancy = redundancy( new_c, base=base )
  rr.complexity = rr.numints <= maxints_for_degen ? complexity5( new_c, base=base ) : 0.0
  rr.degeneracy = rr.numints <= maxints_for_degen ? degeneracy( new_c, base=base ) : 0.0
  rr.sdegeneracy = rr.numints <= maxints_for_degen ? degeneracy( new_c, base=base, mutinf=mutinf2 ) : 0.0
  rr
end

function run_result( p::Parameters, num_goals::Int64, hamming_sel::Bool, active_only::Bool, max_steps::Int64 )
  indiv_result_type(
    p.numinputs,
    p.numoutputs,
    p.numinteriors,
    p.numlevelsback,
    num_goals,
    hamming_sel,
    #avgfitness,
    active_only,
    max_steps,
    #gl_reps,
    #fault_tol,
    #fit_limit,
    0.0,    # hdist
    0,      # steps
    0,      # same
    0,      # worse
    0,      # better
    0,      # nactive
    0.0,    # redundancy
    0.0,    # complexity
    0.0,    # degeneracy
    0.0     # sdegeneracy
  )
end

function run_result_to_tuple( rr::indiv_result_type )
 (
  rr.numinputs,
  rr.numoutputs,
  rr.numints,
  rr.levelsback,
  rr.ngoals,
  #rr.hamming_sel,
  #rr.avgfitness,
  #rr.fault_tol,
  #rr.active_only,
  rr.maxsteps,
  #rr.gl_reps,
  #rr.fit_limit,
  rr.hdist,
  rr.steps,
  rr.same,
  rr.worse,
  rr.better,
  rr.nactive,
  rr.redundancy,
  rr.complexity,
  rr.degeneracy,
  rr.sdegeneracy
 )
end


text:
α = 0.025

@testset "helpers" begin
  @test rand2cop(0.5, 2., "clayton", 0.5) ≈ 0.5463906428428872
  @test rand2cop(0.5, 2., "frank", 0.5) ≈ 0.5
  @test rand2cop(0.5, -.5, "amh", 0.5) ≈ 0.5061130556252271
end

@testset "exceptions" begin
  @test_throws DomainError testbivθ(-2., "clayton")
  @test_throws DomainError usebivρ(-.9, "amh", SpearmanCorrelation)
  @test_throws DomainError usebivρ(-.25, "amh", KendallCorrelation)
  @test_throws AssertionError Chain_of_Archimedeans([2., 3.], ["frank", "gumbel"])
  @test_throws BoundsError Chain_of_Archimedeans([2., 3., 4.], ["frank", "frank"])
  @test_throws DomainError Chain_of_Archimedeans([2., -3.], ["frank", "clayton"])
  @test_throws AssertionError Chain_of_Archimedeans([2., 3., 4.], "gumbel")
  @test_throws DomainError Chain_of_Archimedeans([2., -3.], "clayton")
  @test_throws DomainError Chain_of_Archimedeans([2., -.3], "clayton", KendallCorrelation)
  @test_throws DomainError Chain_of_Archimedeans([.2, -1.3], ["frank", "clayton"], KendallCorrelation)
  @test_throws AssertionError Chain_of_Archimedeans([0.2, 0.2], "gumbel", KendallCorrelation)
  @test_throws AssertionError Chain_of_Archimedeans([.2, .3], ["gumbel", "clayton"], KendallCorrelation)

  u = zeros(2,3)
  c = Chain_of_Archimedeans([2.], ["clayton"])
  @test_throws AssertionError simulate_copula!(u, c)
end

@testset "chain of Archimedean copulas" begin
  @testset "small example" begin
    Random.seed!(43)
    @test simulate_copula(1, Chain_of_Archimedeans([4., 11.], "frank")) ≈ [0.180975  0.492923  0.679345] atol=1.0e-5
    Random.seed!(43)
    @test simulate_copula(1, Chain_of_Archimedeans([4., 11.], ["frank", "clayton"])) ≈ [0.180975  0.492923  0.600322] atol=1.0e-5

    Random.seed!(43)
    u = zeros(1,3)
    simulate_copula!(u, Chain_of_Archimedeans([4., 11.], "frank"))
    @test u ≈ [0.180975  0.492923  0.679345] atol=1.0e-5

    Random.seed!(43)
    u = zeros(1,3)
    simulate_copula!(u, Chain_of_Archimedeans([4., 11.], ["frank", "clayton"]))
    @test u ≈ [0.180975  0.492923  0.600322] atol=1.0e-5
  end

  @testset "larger example" begin
    cops = ["clayton", "clayton", "clayton", "frank", "amh", "amh"]

    Random.seed!(43)
    x = simulate_copula(300000, Chain_of_Archimedeans([-0.9, 3., 2, 4., -0.3, 1.], cops))
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,4], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,5], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,6], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,7], Uniform(0,1))) > α
    @test tail(x[:,1], x[:,2], "l", 0.0001) ≈ 0
    @test tail(x[:,6], x[:,7], "l") ≈ 0.5 atol=1.0e-1
    @test tail(x[:,5], x[:,6], "r", 0.0001) ≈ 0
    @test tail(x[:,3], x[:,4], "l") ≈ 1/(2^(1/2)) atol=1.0e-1
    @test tail(x[:,3], x[:,4], "r", 0.0001) ≈ 0
    @test corkendall(x)[1,2] ≈ -0.9/(2-0.9) atol=1.0e-3
    @test corkendall(x)[2,3] ≈ 3/(2+3) atol=1.0e-3

    # further testing
    Random.seed!(43)
    x = simulate_copula(500000, Chain_of_Archimedeans([-0.9, 2.], "clayton"))
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test tail(x[:,2], x[:,3], "l") ≈ 1/(2^(1/2)) atol=1.0e-1
    @test tail(x[:,2], x[:,3], "r", 0.0001) ≈ 0

  end
end
@testset "correlations" begin
  Random.seed!(43)
  c = Chain_of_Archimedeans([0.6, -0.2], "clayton", SpearmanCorrelation)
  x = simulate_copula(500000, c)
  @test corspearman(x[:,1], x[:,2]) ≈ 0.6 atol=1.0e-2
  @test corspearman(x[:,2], x[:,3]) ≈ -0.2 atol=1.0e-2

  Random.seed!(43)
  c = Chain_of_Archimedeans([0.6, -0.2], ["clayton", "clayton"], KendallCorrelation)
  x = simulate_copula(500000, c)
  @test corkendall(x[:,1], x[:,2]) ≈ 0.6 atol=1.0e-3
  @test corkendall(x[:,2], x[:,3]) ≈ -0.2 atol=1.0e-3
end

@testset "chain of Frechet copulas" begin
  @testset "exceptions" begin
      @test_throws DomainError Chain_of_Frechet([1.1, 0.8])
      @test_throws DomainError Chain_of_Frechet([-0.1, 0.8])
      @test_throws DomainError Chain_of_Frechet([0.1, 0.1], [-0.1, 0.8])
      @test_throws DomainError Chain_of_Frechet([0.1, 0.3], [0.1, 0.8])
      @test_throws DomainError Chain_of_Frechet([0.1, -0.3], [0.1, 1.2])
      @test_throws AssertionError Chain_of_Frechet([0.1, 0.3], [0.1, 0.2, 0.3])
  end
  @testset "test on small example" begin
    @test fncopulagen([0.2, 0.4], [0.1, 0.1], [0.2 0.4 0.6; 0.3 0.5 0.7]) == [0.6 0.4 0.2; 0.7 0.5 0.3]
    Random.seed!(43)
    @test simulate_copula(1, Chain_of_Frechet([0.6, 0.4], [0.3, 0.5])) ≈ [0.888934  0.775377  0.180975] atol=1.0e-5
  end
  @testset "test on large example" begin

    # one parameter case
    Random.seed!(43)
    x = simulate_copula(500000, Chain_of_Frechet([0.9, 0.6, 0.2]))
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test corspearman(x) ≈ [1. 0.9 0.6 0.2; 0.9 1. 0.6 0.2; 0.6 0.6 1. 0.2; 0.2 0.2 0.2 1.] atol=1.0e-2
    @test tail(x[:,1], x[:,2], "r") ≈ 0.9 atol=1.0e-1
    @test tail(x[:,1], x[:,2], "l") ≈ 0.9 atol=1.0e-1
    @test tail(x[:,1], x[:,4], "r") ≈ 0.2 atol=1.0e-1

    # two parameters case
    Random.seed!(43)
    x = simulate_copula(500000, Chain_of_Frechet([0.8, 0.5], [0.2, 0.3]));
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test corspearman(x) ≈ [1. 0.6 0.2; 0.6 1. 0.2; 0.2 0.2 1.] atol=1.0e-3
    @test tail(x[:,1], x[:,2], "r") ≈ 0.8 atol=1.0e-1
    @test tail(x[:,2], x[:,3], "r") ≈ 0.5 atol=1.0e-1
  end
end

@testset "Big Float implementation" begin

  cops = ["clayton", "clayton", "clayton", "frank", "amh", "amh"]

  Random.seed!(43)
  θs = BigFloat.([-0.9, 3., 2, 4., -0.3, 1.])
  x = simulate_copula(1000, Chain_of_Archimedeans(θs, cops))
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,4], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,5], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,6], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,7], Uniform(0,1))) > α
  @test typeof(x) == Array{BigFloat,2}

  va = BigFloat.([0.9, 0.6, 0.2])
  Random.seed!(43)
  x = simulate_copula(1000, Chain_of_Frechet(va))
  @test typeof(x) == Array{BigFloat,2}
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α

  va = BigFloat.([0.79999, 0.5])
  vb = BigFloat.([0.2, 0.3])
  Random.seed!(43)
  x = simulate_copula(1000, Chain_of_Frechet(va, vb))
  @test typeof(x) == Array{BigFloat,2}
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α

end


text:
export generate_noise_mpi

import CorrNoise
using Random
using FITSIO

try
    import MPI
catch
end

"""
This structure holds informations about a certain chunk of data produced by a certain MPI process to send to another MPI process.

Field              | Type           | Meaning
:----------------- |:-------------- |:------------------------------------------------------------
`pol_number`       | Int            | ID number of the polarimeter
`first_idx`        | Int            | Index of the first element ofthe chunk
`last_idx`         | Int            | Index of the last element of the chunk
`dest_rank`        | Int            | Number of destination rank

# Example
Suppose that the following array of data: [1, 10, 20, 3, 4, 5, 7]
is measured by the polarimeter number 2.
If the chunk [20, 3, 4] needs to be sent to rank number 3, the corresponding structure will be:

chunk_to_send(2, 3, 5, 3)
"""
struct chunk_to_send
    pol_number ::Int
    first_idx::Int
    last_idx::Int
    dest_rank::Int
end

function generate_noisechunks_to_send(chunks, baseline_length_s, fsamp_hz, total_number_of_polarimeters)

    commsize = length(chunks)
    noisechunks_to_send = []

    for cur_rank in 0:(commsize-1) #loop on ranks
        pol_cur_rank  = cur_rank+1 #polarimeter assigned to cur_rank for noise production
        chunk_cur_rank = []

        #Only the first "total_number_of_polarimeters" ranks will generate noise, so they will have something to send.
        if  0 <= cur_rank < total_number_of_polarimeters
            for i in 1:length(chunks) #loop on chunks
                rank_it = i-1
                for j in 1:length(chunks[i]) #loop on detectors per rank
                    if chunks[i][j].pol_number == pol_cur_rank #if cur_rank is going to produce noise for the "correct" polarimeter
                        first_idx = (chunks[i][j].first_idx-1)*baseline_length_s*fsamp_hz+1
                        last_idx = first_idx + chunks[i][j].num_of_elements*baseline_length_s*fsamp_hz -1
                        dest_rank = rank_it
                        chunk = chunk_to_send(pol_cur_rank, first_idx, last_idx, dest_rank)
                        chunk_cur_rank = append!(chunk_cur_rank, [chunk])
                    end
                end
            end
            noisechunks_to_send= append!(noisechunks_to_send, [chunk_cur_rank])
        end
    end

    return noisechunks_to_send

end

"""
    function generate_noise_mpi(chunks, baseline_length_s, baselines_per_process, fsamp_hz, σ_k, fknee_hz, rank, comm)

    This MPI based function can be used to generate white noise and 1/f noise when the simulation is splitted in various MPI processes.

    It requires in input:
    - the data chunks (which can be obtained by using the function `split_tod_mpi`)
    - the array containing the number of 1/f baselines to simulate for each process.
            It can be obtained by using the function `split_into_n` in the following way:

            split_into_n(num_of_polarimeters*baselines_per_pol, num_of_MPI_proc)

            where baselines_per_pol = Int64(total_time/baseline_length_s) is the number of baselines of each polarimeter
    - the length (in seconds) of each 1/f noise baseline
    - the total duration (in seconds) of data acquisition
    - the sampling frequency (in Hz)
    - an array containing the RMS of the temperature (in K) for each polarimeter, σ_k[i]= tsys_k[i] / sqrt(β_hz[i] * τ_s)
    - an array containing the knee frequency (in Hz) for each polarimeter
    - an array containing the oof slope for each polarimeter
    - the MPI communicator
    - the random generator seed, as optional argument.

    It returns the noise TOD for the current rank.

    The function behaves differently according to the number of processes used:

        1. If the number of processes is greater than the number of polarimeters:
           the first "number of polarimeters" processes will simultaneously produce the noise, one process for one polarimeter
           (rank 0 for pol.1, rank 1 for pol.2 and so on).
           After that, the noise is redistributed between the ranks according to the input chunk structure.

        2. If the number of processes is lesser than the number of polarimeters,
           Each partial noise TOD is generated by rank 0 and then sent to the correspondent rank.

    N.B. if you want to use this function without MPI, remember to put rank = 0 and comm = missing

"""
function  generate_noise_mpi(chunks, baselines_per_process, baseline_length_s, total_time, fsamp_hz, σ_k, fknee_hz, slope, rank, comm, input_seed = missing)

    total_number_of_polarimeters = length(σ_k)
    commsize = length(chunks)

    if ismissing(input_seed) father_seed = rand(1:10000) else father_seed = input_seed end
    Random.seed!(father_seed)
    seeds = Int64[rand(1:10000) for i in 1:total_number_of_polarimeters]

    ################ IF MORE PROCESSES THAN POLARIMETERS ################
    if(commsize >= total_number_of_polarimeters)

        noisechunks_to_send = generate_noisechunks_to_send(chunks, baseline_length_s, fsamp_hz, total_number_of_polarimeters)

        ##### NOISE PRODUCTION #####
        pol_number = rank+1
        pol_noise = Float64[]
        noise = Array{Float64}(undef, baselines_per_process[rank+1]*baseline_length_s*fsamp_hz)

        if  0 <= rank < total_number_of_polarimeters  #if first "total_number_of_polarimeters" ranks
            rng = CorrNoise.OofRNG(Random.MersenneTwister(seeds[pol_number]), -abs(slope[pol_number]), 1.15e-5, fknee_hz[pol_number], fsamp_hz)
            samples_per_pol = total_time*fsamp_hz
            pol_noise = Float64[CorrNoise.randoof(rng) * σ_k[pol_number] for i in 1:(samples_per_pol)]
        end

        if(!ismissing(comm))
            MPI.Barrier(comm)
        end

        ##### SEND AND RECEIVE #####

        start_idx = 1
        end_idx = 0

        for i in 1:length(noisechunks_to_send)  #loop on ranks
            cur_rank = i-1

            for j in 1:length(noisechunks_to_send[i]) #loop on the noise chunks produced by cur_rank
                start_idx_noise = noisechunks_to_send[i][j].first_idx
                end_idx_noise = noisechunks_to_send[i][j].last_idx
                dest_rank = noisechunks_to_send[i][j].dest_rank

                cur_num_noise_samples = end_idx_noise-start_idx_noise+1

                #if rank has to send noise (since it produced it)
                if cur_rank == rank

                    #if the destination rank is itself, keep noise chunk
                    if dest_rank == rank
                        end_idx = start_idx + cur_num_noise_samples - 1
                        noise[start_idx:end_idx] = pol_noise[start_idx_noise:end_idx_noise]
                        start_idx += cur_num_noise_samples

                    #otherwise send to the destination rank
                    else
                        cur_noise_send = pol_noise[start_idx_noise:end_idx_noise]
                        MPI.Send(cur_noise_send, dest_rank , 0, comm)
                    end

                #if rank is a destination rank for someone
                elseif dest_rank == rank

                    cur_noise_rec = Array{Float64}(undef, cur_num_noise_samples)
                    end_idx = start_idx + cur_num_noise_samples -1
                    MPI.Recv!(cur_noise_rec, cur_rank, 0, comm)
                    noise[start_idx:end_idx] = cur_noise_rec
                    start_idx += cur_num_noise_samples
                end

            end

        end

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

    ################IF LESS PROCESSES THAN POLARIMETERS################
    else
        if rank == 0
            previous_detector = chunks[1][1].pol_number
            noise = Float64[]  #noise rank 0

            rng = CorrNoise.OofRNG(Random.MersenneTwister(seeds[1]), -abs(slope[1]), 1.15e-5, fknee_hz[1], fsamp_hz)

            for i in 1:length(chunks)   #loop on ranks
                num_noise_samples = baselines_per_process[i]*baseline_length_s*fsamp_hz

                for j in 1:length(chunks[i]) #loop on detectors per rank
                    cur_detector = chunks[i][j].pol_number

                    if cur_detector != previous_detector  #if new detector generate new noise
                        rng = CorrNoise.OofRNG(Random.MersenneTwister(seeds[cur_detector]), -abs(slope[cur_detector]), 1.15e-5, fknee_hz[cur_detector], fsamp_hz)
                    end
                    cur_num_noise_samples = chunks[i][j].num_of_elements*baseline_length_s*fsamp_hz
                    cur_noise = Float64[CorrNoise.randoof(rng) * σ_k[cur_detector] for i in 1:(cur_num_noise_samples)]
                    previous_detector =  cur_detector

                    if (i > 1)  #send to other ranks, not rank 0
                        MPI.Send(cur_noise, i-1, 0, comm)
                    else
                        noise = append!(noise, cur_noise)
                    end
                end
            end
        else
            noise = Float64[]
            for j in 1:length(chunks[rank+1])
                cur_noise = Array{Float64}(undef, chunks[rank+1][j].num_of_elements*baseline_length_s*fsamp_hz)
                MPI.Recv!(cur_noise, 0, 0, comm)
                noise = append!(noise, cur_noise)
            end
        end

    end

    return noise

end


text:
struct PinnedFoot <: Constraints
    n
    ineq
    q_pinned
end

function pinned_foot_constraint(model, q_pinned, T)
    n = 2 * (T + 1)
    ineq = (1:0)
    PinnedFoot(n, ineq, q_pinned)
end

function constraints!(c, Z, con::PinnedFoot, model, idx, h, T)
    k(y) = kinematics_2(model, y, body = :leg_1, mode = :ee)

    for t = 1:T
        if t == 1
            q⁻ = view(Z, idx.x[t][1:model.nq])
            c[1:2] = k(q⁻) - k(con.q_pinned)
        end
        q = view(Z, idx.x[t][model.nq .+ (1:model.nq)])
        c[2 + (t - 1) * 2 .+ (1:2)] = k(q) - k(con.q_pinned)
    end
    nothing
end

function constraints_jacobian!(∇c, Z, con::PinnedFoot, model, idx, h, T)
    k(y) = kinematics_2(model, y, body = :leg_1, mode = :ee)
    shift = 0
    for t = 1:T
        if t == 1
            q⁻ = view(Z, idx.x[t][1:model.nq])
            r_idx = (1:2)
            c_idx = idx.x[t][1:model.nq]
            len = length(r_idx) * length(c_idx)
            ∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(k, q⁻))
            shift += len
        end
        q = view(Z, idx.x[t][model.nq .+ (1:model.nq)])
        r_idx = 2 + (t - 1) * 2 .+ (1:2)
        c_idx = idx.x[t][model.nq .+ (1:model.nq)]
        len = length(r_idx) * length(c_idx)
        ∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(k, q))
        shift += len
    end

    return nothing
end

function constraints_sparsity(con::PinnedFoot, model, idx, T;
    shift_row = 0, shift_col = 0)
    row = []
    col = []

    for t = 1:T
        if t == 1
            r_idx = shift_row .+ (1:2)
            c_idx = shift_col .+ idx.x[t][1:model.nq]
            row_col!(row, col, r_idx, c_idx)
        end
        r_idx = shift_row + 2 + (t - 1) * 2 .+ (1:2)
        c_idx = shift_col .+ idx.x[t][model.nq .+ (1:model.nq)]
        row_col!(row, col, r_idx, c_idx)
    end

    return collect(zip(row, col))
end

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
# unit vectors and identity matrix
const e1 = SVector(1, 0, 0)
const e2 = SVector(0, 1, 0)
const e3 = SVector(0, 0, 1)
const I3 = @SMatrix [1 0 0; 0 1 0; 0 0 1]
const Z3 = @SMatrix zeros(3,3)

"""
    tilde(x)

Construct the cross product operator matrix
"""
@inline tilde(x) = @SMatrix [0 -x[3] x[2]; x[3] 0 -x[1]; -x[2] x[1] 0]

"""
    wiener_milenkovic(c)

Construct a Wiener-Milenkovic transformation matrix, given the three Wiener-Milenkovic
parameters in `c`.

Note that the corresponding rotation matrix is the transpose of this
transformation matrix.
"""
@inline function wiener_milenkovic(c)

    c0 = 2 - c'*c/8

    C = 1/(4-c0)^2*(@SMatrix [
        c0^2 + c[1]^2 - c[2]^2 - c[3]^2  2*(c[1]*c[2] + c0*c[3])          2*(c[1]*c[3]-c0*c[2]);
        2*(c[1]*c[2] - c0*c[3])          c0^2 - c[1]^2 + c[2]^2 - c[3]^2  2*(c[2]*c[3] + c0*c[1]);
        2*(c[1]*c[3] + c0*c[2])          2*(c[2]*c[3] - c0*c[1])          c0^2 - c[1]^2 - c[2]^2 + c[3]^2
        ]
    )

    return C
end

"""
    transform_properties(K, T)

Applies the transformation `T` to the stiffness or mass matrix `K`
"""
transform_properties(X, T) = [T' Z3; Z3 T'] * X * [T Z3; Z3 T]

"""
    rotation_parameter_scaling(θ)

Extract a scaling parameter which may be multiplied by the angular parameters
to yield the Wiener-Milenkovic rotation parameters.  Use of this scaling
parameter allows deflections greater than 360 degrees.
"""
@inline function rotation_parameter_scaling(θ)

    θ = real(θ)

    θnorm = sqrt(θ'*θ)

    if θnorm <= 4
        scaling = one(θnorm)
    else
        half_rotations = θnorm/4
        irot, m = floor(half_rotations/2), half_rotations%2
        if m <= 1
            scaling = 1 .- 8*irot/θnorm
        else
            scaling = 1 .- 8*(irot+1)/θnorm
        end
    end

    return scaling
end

"""
    get_C(θ)

Returns the transformation matrix `C` given the three angular parameters in `θ`.
"""
@inline function get_C(θ)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ

    return wiener_milenkovic(c)
end

"""
    get_C_t([C, ] θ, θ_t)

Calculate the derivative of the Wiener-Milenkovic transformation matrix `C` with
respect to the scalar parameter `t`. `θ_t` is the derivative of the angular parameter
`θ` with respect to `t`.
"""
@inline get_C_t(θ, θ_t) = get_C_t(get_C(θ), θ, θ_t)

@inline function get_C_t(C, θ, θ_t)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ
    c_t = scaling*θ_t

    c0 = 2 - c'*c/8
    c0dot = -c'*c_t/4
    tmp = 1/(4-c0)^2
    tmpdot = 2/(4-c0)^3*c0dot

    Cdot = tmpdot*C/tmp + tmp*(@SMatrix [
        2*c0*c0dot + 2*c[1]*c_t[1] - 2*c[2]*c_t[2] - 2*c[3]*c_t[3] (
        2*(c_t[1]*c[2] + c[1]*c_t[2] + c0dot*c[3] + c0*c_t[3])) (
        2*(c_t[1]*c[3] + c[1]*c_t[3] - c0dot*c[2] - c0*c_t[2]));
        2*(c_t[1]*c[2] + c[1]*c_t[2] - c0dot*c[3] - c0*c_t[3]) (
        2*c0*c0dot - 2*c[1]*c_t[1] + 2*c[2]*c_t[2] - 2*c[3]*c_t[3]) (
        2*(c_t[2]*c[3] + c[2]*c_t[3] + c0dot*c[1] + c0*c_t[1]));
        2*(c_t[1]*c[3] + c[1]*c_t[3] + c0dot*c[2] + c0*c_t[2]) (
        2*(c_t[2]*c[3] + c[2]*c_t[3] - c0dot*c[1] - c0*c_t[1])) (
        2*c0*c0dot - 2*c[1]*c_t[1] - 2*c[2]*c_t[2] + 2*c[3]*c_t[3])]
    )

    return Cdot
end

"""
    get_C_θ([C, ] θ)

Calculate the derivative of the Wiener-Milenkovic transformation matrix `C` with
respect to each of the rotation parameters in `θ`.
"""
@inline get_C_θ(c) = get_C_θ(get_C(c), c)

@inline function get_C_θ(C, θ)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ

    c0 = 2 - c'*c/8
    c0_θ = -c/4

    tmp = 1/(4-c0)^2
    tmp_θ = 2/(4-c0)^3*c0_θ

    C_θ1 = tmp_θ[1]*C/tmp + tmp*(@SMatrix [
        2*c0*c0_θ[1] + 2*c[1]    2*(c[2] + c0_θ[1]*c[3]) 2*(c[3] - c0_θ[1]*c[2]);
         2*(c[2] - c0_θ[1]*c[3])  2*c0*c0_θ[1] - 2*c[1]   2*(c0_θ[1]*c[1] + c0);
         2*(c[3] + c0_θ[1]*c[2]) -2*(c0_θ[1]*c[1] + c0)   2*c0*c0_θ[1] - 2*c[1]
        ]
    )

    C_θ2 = tmp_θ[2]*C/tmp + tmp*(@SMatrix [
        2*c0*c0_θ[2] - 2*c[2]   2*(c[1] + c0_θ[2]*c[3]) -2*(c0_θ[2]*c[2] + c0);
        2*(c[1] - c0_θ[2]*c[3]) 2*c0*c0_θ[2] + 2*c[2]    2*(c[3] + c0_θ[2]*c[1]);
         2*(c0_θ[2]*c[2] + c0)   2*(c[3] - c0_θ[2]*c[1])  2*c0*c0_θ[2] - 2*c[2]
        ]
    )

    C_θ3 = tmp_θ[3]*C/tmp + tmp*(@SMatrix [
         2*c0*c0_θ[3] - 2*c[3]   2*(c0_θ[3]*c[3] + c0)    2*(c[1] - c0_θ[3]*c[2]);
        -2*(c0_θ[3]*c[3] + c0)   2*c0*c0_θ[3] - 2*c[3]    2*(c[2] + c0_θ[3]*c[1]);
         2*(c[1] + c0_θ[3]*c[2]) 2*(c[2] - c0_θ[3]*c[1])  2*c0*c0_θ[3] + 2*c[3]
        ]
    )

    return scaling*C_θ1, scaling*C_θ2, scaling*C_θ3
end

"""
    get_C_t_θ(θ, θ_t)

Calculate the derivative of the time derivative of the Wiener-Milenkovic
transformation matrix `C` with respect to `θ`.
"""
@inline function get_C_t_θ(θ, θ_t)

    θ1, θ2, θ3 = θ
    C_t_θ1 = ForwardDiff.derivative(θ1 -> get_C_t(SVector(θ1, θ2, θ3), θ_t), θ1)
    C_t_θ2 = ForwardDiff.derivative(θ2 -> get_C_t(SVector(θ1, θ2, θ3), θ_t), θ2)
    C_t_θ3 = ForwardDiff.derivative(θ3 -> get_C_t(SVector(θ1, θ2, θ3), θ_t), θ3)

    return C_t_θ1, C_t_θ2, C_t_θ3
end

"""
    get_C_t_θdot([C, ] θ)

Calculate the derivative of the time derivative of the Wiener-Milenkovic
transformation matrix `C` with respect to each of the time derivatives of `θ`.
"""
get_C_t_θdot

@inline get_C_t_θdot(θ) = get_C_t_θdot(get_C(θ), θ)

@inline function get_C_t_θdot(C, θ)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ

    c0 = 2 - c'*c/8
    c0dot_θdot = -c/4
    tmp = 1/(4-c0)^2
    tmpdot_θdot = 2/(4-c0)^3*c0dot_θdot

    C_θdot1 = tmpdot_θdot[1]*C/tmp + tmp*(@SMatrix [
        2*c0*c0dot_θdot[1] + 2*c[1]   2*(c[2] + c[3]*c0dot_θdot[1]) 2*(c[3] - c[2]*c0dot_θdot[1]);
        2*(c[2] - c[3]*c0dot_θdot[1]) 2*c0*c0dot_θdot[1] - 2*c[1]   2*(c[1]*c0dot_θdot[1] + c0);
        2*(c[3] + c[2]*c0dot_θdot[1]) 2*(-c[1]*c0dot_θdot[1] - c0)  2*c0*c0dot_θdot[1] - 2*c[1]]
    )

    C_θdot2 = tmpdot_θdot[2]*C/tmp + tmp*(@SMatrix [
        2*c0*c0dot_θdot[2] - 2*c[2]   2*(c[1] + c[3]*c0dot_θdot[2]) 2*(-c[2]*c0dot_θdot[2] - c0);
        2*(c[1] - c[3]*c0dot_θdot[2]) 2*c0*c0dot_θdot[2] + 2*c[2]   2*(c[3] + c[1]*c0dot_θdot[2]);
        2*(c[2]*c0dot_θdot[2] + c0)   2*(c[3] - c[1]*c0dot_θdot[2]) 2*c0*c0dot_θdot[2]- 2*c[2]]
    )

    C_θdot3 = tmpdot_θdot[3]*C/tmp + tmp*(@SMatrix [
        2*c0*c0dot_θdot[3] - 2*c[3] 2*(c[3]*c0dot_θdot[3] + c0) 2*(c[1] - c[2]*c0dot_θdot[3]);
        2*(-c[3]*c0dot_θdot[3] - c0) 2*c0*c0dot_θdot[3] - 2*c[3] 2*(c[2] + c[1]*c0dot_θdot[3]);
        2*(c[1] + c[2]*c0dot_θdot[3]) 2*(c[2] - c[1]*c0dot_θdot[3]) 2*c0*c0dot_θdot[3] + 2*c[3]]
    )

    return scaling*C_θdot1, scaling*C_θdot2, scaling*C_θdot3
end

"""
    get_Q(θ)

Calculate the matrix Q as defined in the paper "Geometrically nonlinear analysis
of composite beams using Wiener-Milenković parameters" by Qi Wang and Wenbin Yu
given the rotational parameters in `θ`.
"""
@inline function get_Q(θ)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ

    c0 = 2 - c'*c/8

    return 1/(4-c0)^2*((4 - 1/4*c'*c)*I - 2*tilde(c) + 1/2*c*c')
end

"""
    get_Q_θ(θ)
    get_Q_θ(Q, θ)

Calculate the derivative of the matrix `Q` with respect to each of the rotation
parameters in `c`.
"""
@inline get_Q_θ(θ) = get_Q_θ(get_Q(θ), θ)

@inline function get_Q_θ(Q, θ)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ

    c0 = 2 - c'*c/8
    c0_θ = -c/4

    tmp = 1/(4-c0)^2
    tmp_θ = 2/(4-c0)^3*c0_θ

    Q_θ1 = tmp_θ[1]*Q/tmp + tmp*(-c[1]/2*I - 2*tilde(e1) + (e1*c' + c*e1')/2)

    Q_θ2 = tmp_θ[2]*Q/tmp + tmp*(-c[2]/2*I - 2*tilde(e2) + (e2*c' + c*e2')/2)

    Q_θ3 = tmp_θ[3]*Q/tmp + tmp*(-c[3]/2*I - 2*tilde(e3) + (e3*c' + c*e3')/2)

    return scaling*Q_θ1, scaling*Q_θ2, scaling*Q_θ3
end

"""
    get_Qinv(θ)

Calculate the matrix inverse `Qinv` as defined in the paper "Geometrically
nonlinear analysis of composite beams using Wiener-Milenković parameters" by
Qi Wang and Wenbin Yu given the rotational parameters in `θ`.
"""
@inline function get_Qinv(θ)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ

    Qinv = (1 - 1/16*c'*c)*I3 + 1/2*tilde(c) + 1/8*c*c'

    return Qinv
end

"""
    get_Qinv_θ(θ)

Calculate the derivative of the matrix inverse `Qinv` with respect to each of the
rotation parameters in `θ`.
"""
@inline function get_Qinv_θ(θ)

    scaling = rotation_parameter_scaling(θ)

    c = scaling*θ

    Qinv_θ1 = -c[1]/8*I3 + 1/2*tilde(e1) + 1/8*(e1*c' + c*e1')
    Qinv_θ2 = -c[2]/8*I3 + 1/2*tilde(e2) + 1/8*(e2*c' + c*e2')
    Qinv_θ3 = -c[3]/8*I3 + 1/2*tilde(e3) + 1/8*(e3*c' + c*e3')

    return scaling*Qinv_θ1, scaling*Qinv_θ2, scaling*Qinv_θ3
end

"""
    mul3(A_1, A_2, A_3, b)

Return the product of a 3x3x3 tensor represented by `A_1`, `A_2`, and `A_3` with
the vector `b`.
"""
@inline mul3(A_1, A_2, A_3, b) = hcat(A_1*b, A_2*b, A_3*b)

"""
    gauss_quadrature(f, a, b)

Default gauss-quadrature function used for integrating distributed loads.
"""
@inline function gauss_quadrature(f, a, b)
    h = b - a
    c = (a + b)/2
    x = h/2*GAUSS_NODES .+ c
    return h/2*GAUSS_WEIGHTS'*f.(x)
end

# this function is not used
# @inline function linf_norm_scaling(A)
#
#     droptol!(A, eps(eltype(A)))
#
#     m, n = size(A)
#     rows = rowvals(A)
#     vals = nonzeros(A)
#
#     # initialize outputs
#     rowmax = fill(typemin(eltype(A)), n)
#     rowmin = fill(typemax(eltype(A)), n)
#     @inbounds for j = 1:n
#         @inbounds for i in nzrange(A, j)
#            row = rows[i]
#            val = abs(vals[i])
#            # perform sparse wizardry...
#            if val > rowmax[row]
#                rowmax[row] = val
#            end
#            if val < rowmin[row]
#                rowmin[row] = val
#            end
#         end
#     end
#
#     # use storage already allocated for rowmax for r
#     r = rowmax
#     @inbounds for i = 1:m
#         if i in rows
#             r[i] = nextpow(2, 1/sqrt(rowmax[i]*rowmin[i]))
#         else
#             r[i] = 1
#         end
#     end
#
#     # use storage already allocated for rowmin for s
#     s = rowmin
#     @inbounds for j = 1:n
#         colrange = nzrange(A, j)
#         if isempty(colrange)
#             s[j] = 1
#         else
#             colmax = maximum(i -> r[rows[i]]*abs(vals[i]), colrange)
#             colmin = minimum(i -> r[rows[i]]*abs(vals[i]), colrange)
#             s[j] = nextpow(2, 1/sqrt(colmax*colmin))
#         end
#     end
#
#     return r, s
# end


text:
# 1D linear diffusion Julia MPI solver
# run: ~/.julia/bin/mpiexecjl -n 4 julia --project diffusion_1D_mpi.jl
using Plots, Printf, MAT
import MPI

# enable plotting by default
if !@isdefined do_save; do_save = true end

# MPI functions
@views function update_halo(A, neighbors_x, comm)
    # Send to / receive from neighbor 1 ("left neighbor")
    if neighbors_x[1] != MPI.MPI_PROC_NULL
        sendbuf = A[2]
        recvbuf = zeros(1)
        MPI.send(sendbuf, neighbors_x[1],0,comm)
        MPI.Recv!(recvbuf, neighbors_x[1],1,comm)
        A[1] = recvbuf
    end
    # Send to / receive from neighbor 2 ("right neighbor")
    if neighbors_x[2] != MPI.MPI_PROC_NULL
        sendbuf = A[end-1]
        recvbuf = zeros(1)
        MPI.send(sendbuf, neighbors_x[2],1,comm)
        MPI.Recv!(recvbuf, neighbors_x[2],0,comm)
        A[end] = recvbuf
    end
    return
end

@views function diffusion_1D_mpi(; do_save=false)
    # MPI
    MPI.Init()
    dims        = [0]
    comm        = MPI.COMM_WORLD
    nprocs      = MPI.Comm_size(comm)
    MPI.Dims_create!(nprocs, dims)
    comm_cart   = MPI.Cart_create(comm, dims, [0], 1)
    me          = MPI.Comm_rank(comm_cart)
    coords      = MPI.Cart_coords(comm_cart)
    neighbors_x = MPI.Cart_shift(comm_cart, 0, 1)
    if (me==0) println("nprocs=$(nprocs), dims[1]=$(dims[1])") end
    # Physics
    lx    = 10.0
    D     = 1.0
    nt    = 100
    # Numerics
    nx    = 32                 # local number of grid points
    nx_g  = dims[1]*(nx-2) + 2 # global number of grid points
    # Derived numerics
    dx    = lx/nx_g            # global
    dt    = dx^2/D/2.1
    # Array allocation
    qx    = zeros(nx-1)
    # Initial condition
    x0    = coords[1]*(nx-2)*dx
    xc    = [x0 + ix*dx - dx/2 - 0.5*lx  for ix=1:nx]
    C     = exp.(.-xc.^2)
    t_tic = 0.0
    # Time loop
    for it = 1:nt
        if (it==11) t_tic = Base.time() end
        qx         .= .-D*diff(C)/dx
        C[2:end-1] .= C[2:end-1] .- dt*diff(qx)/dx
        update_halo(C, neighbors_x, comm_cart)
    end
    t_toc = Base.time()-t_tic
    if (me==0) @printf("Time = %1.4e s, T_eff = %1.2f GB/s \n", t_toc, round((2/1e9*nx*sizeof(lx))/(t_toc/(nt-10)), sigdigits=2)) end
    if do_save file = matopen("$(@__DIR__)/mpi1D_out_C_$(me).mat", "w"); write(file, "C", Array(C)); close(file) end
    MPI.Finalize()
    return
end

diffusion_1D_mpi(; do_save=do_save)

max_stars_repo_name: Wilboludriver/Julia-practice

max_issues_repo_name: Wilboludriver/Julia-practice

max_forks_repo_name: Wilboludriver/Julia-practice


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "CUDDBuilder"
version = v"3.0.0"

# Collection of sources required to build CUDDBuilder
sources = [
    "http://davidkebo.com/source/cudd_versions/cudd-3.0.0.tar.gz" =>
    "b8e966b4562c96a03e7fbea239729587d7b395d53cadcc39a7203b49cf7eeb69",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd cudd-3.0.0

if [[ "${target}" == *-mingw* ]]; then
    sed -i 's/cross_compiling=no/cross_compiling=yes/' configure
    ./configure --prefix=$prefix --target=${target} --build=x86_64-w64-mingw32 --enable-silent-rules --enable-shared --enable-obj
elif [[ "${target}" == *-apple-* ]]; then
    sed -i 's/cross_compiling=no/cross_compiling=yes/' configure
    export CC=/opt/${target}/bin/${target}-gcc
    export CXX=/opt/${target}/bin/${target}-g++
    ./configure --prefix=$prefix --target=${target} --enable-silent-rules --enable-shared --enable-obj
else
    autoreconf -fi
    sed -i 's/cross_compiling=no/cross_compiling=yes/' configure
    export CC=/opt/${target}/bin/${target}-gcc
    export CXX=/opt/${target}/bin/${target}-g++
    export LD=/opt/${target}/bin/${target}-ld
    export LDFLAGS=-L/opt/${target}/${target}/lib64

    ./configure --prefix=$prefix --target=${target} --enable-silent-rules --enable-shared --enable-obj
fi

make -j${nproc} check
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc),
    MacOS(:x86_64),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libcudd", Symbol(""))
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
export summary_statistics, basic_summary_statistics

function data_size_cutoff(simplices::Vector{Int64}, nverts::Vector{Int64},
                          times::Vector{Int64}, min_simplex_size::Int64,
                          max_simplex_size::Int64)
    new_simplices = Int64[]
    new_nverts = Int64[]
    new_times = Int64[]
    curr_ind = 1
    for (nv, time) in zip(nverts, times)
        if min_simplex_size <= nv <= max_simplex_size
            append!(new_simplices, simplices[curr_ind:(curr_ind + nv - 1)])
            push!(new_nverts, nv)
            push!(new_times, time)
        end
        curr_ind += nv
    end
    return (new_simplices, new_nverts, new_times)
end

"""
summary_statistics
------------------------

Computes several statistics about the dataset.

summary_statistics(dataset::HONData)

Input parameter:
- dataset::HONData: The dataset.
"""
function summary_statistics(dataset::HONData)
    function _summary_statistics(simplices::Vector{Int64}, nverts::Vector{Int64}, times::Vector{Int64})
        if length(nverts) == 0
            return ("0,0,0,0,0,0.0,0.0,0,0,0,0,0.0,0.0,0.0,0.0,0,0")
        end

        A, At, B = basic_matrices(simplices, nverts)
        no, nc = num_open_closed_triangles(A, At, B)
        num_nodes = sum(sum(At, dims=1) .> 0)  # note: includes 1-node simplices
        density = nnz(B) / (num_nodes^2 - num_nodes)
        str1 = @sprintf("%d,%d,%d,%d,%d,%f,%e,%d,%d",
                        num_nodes, length(nverts), nnz(A), mean(nverts),
                        maximum(nverts), mean(nonzeros(B)), density, nc, no)

        # Backbone dataset
        bb_simplices, bb_nverts, bb_times = backbone(simplices, nverts, times)
        (A, At, C) = basic_matrices(bb_simplices, bb_nverts)
        str2 = @sprintf("%d,%d,%f,%f", length(bb_nverts), nnz(A),
                        mean(bb_nverts), mean(nonzeros(C)))

        # Random configuration
        config = configuration_sizes_preserved(bb_simplices, bb_nverts)
        (A, At, B) = basic_matrices(config, bb_nverts)
        num_nodes = sum(sum(At, dims=1) .> 0)
        density = nnz(B) / (num_nodes^2 - num_nodes)
        no, nc = num_open_closed_triangles(A, At, B)
        str3 = @sprintf("%f,%e,%d,%d",
                        mean(nonzeros(B)), density, nc, no)
        return "$str1,$str2,$str3"
    end

    dataset_name = dataset.name
    open("$(dataset_name)-statistics.csv", "w") do f
        write(f, "dataset,")
        write(f, "nnodes,nsimps,nconns,meansimpsize,maxsimpsize,meanprojweight,projdensity,nclosedtri,nopentri,")
        write(f, "nbbnodes,nbbconns,meanbbsimpsize,meanbbprojweight,")
        write(f, "meanbbconfigweight,bbconfigdensity,nbbconfigclosedtri,nbbconfigopentri\n")
        simplices, nverts, times = dataset.simplices, dataset.nverts, dataset.times
        stats = _summary_statistics(simplices, nverts, times)
        write(f, "$(dataset_name),$stats\n")

        # Get the same statistics but for the dataset restricted to 3-node
        # simplices.
        csimplices, cnverts, ctimes = data_size_cutoff(simplices, nverts, times, 3, 3)
        stats = _summary_statistics(csimplices, cnverts, ctimes)
        write(f, "$(dataset_name)-3-3,$stats\n")
    end
end

"""
basic_summary_statistics
------------------------

Prints some basic summary statistics of the dataset.

basic_summary_statistics(dataset::String)

Input parameter:
- dataset::HONData: The dataset.
"""
function basic_summary_statistics(dataset::HONData)
    simplices, nverts, times = dataset.simplices, dataset.nverts, dataset.times
    bb_simplices, bb_nverts, bb_times = backbone(simplices, nverts, times)
    A, At, B = basic_matrices(simplices, nverts)
    num_nodes = sum(sum(At, dims=1) .> 0)
    num_edges = nnz(B) / 2
    num_simps = length(nverts)
    num_bb_simps = length(bb_nverts)
    println("dataset & # nodes & # edges in proj. graph & # simplices & # unique simplices")
    @printf("%s & %d & %d & %d & %d\n", dataset.name, num_nodes, num_edges, num_simps, num_bb_simps)
end


text:
"""
    Sequential linear programming with trust region

This is based on the paper: https://doi.org/10.1016/j.epsr.2018.09.002

# Note
- The predicted reduction ϕ_pre implemented in the original paper failed to solve many test instances. So, we have modified the predicted reduction ϕ_pre to use the directional derivative.
- We have implemented the feasibility restoration, which is not available in the original paper.
"""
mutable struct SlpTR{T,Tv,Tt} <: AbstractSlpOptimizer
    problem::Model{T,Tv,Tt} # problem data

    x::Tv # primal solution
    p::Tv # search direction
    p_slack::Dict{Int,Tv} # search direction at feasibility restoration phase
    lambda::Tv # Lagrangian dual multiplier
    mult_x_L::Tv # reduced cost for lower bound
    mult_x_U::Tv # reduced cost for upper bound

    # Evaluations at `x`
    f::T # objective function
    df::Tv # gradient
    E::Tv # constraint
    dE::Tv # Jacobian

    phi::T # merit function value
    ν::Tv # penalty parameters for the merit function

    Δ::T # current trust region size
    Δ_max::T # maximum trust region size allowed
    alpha1::T # parameter to adjust trust region size
    alpha2::T # parameter to adjust trust region size

    prim_infeas::T # primal infeasibility at `x`
    dual_infeas::T # dual (approximate?) infeasibility
    compl::T # complementary slackness

    optimizer::Union{Nothing,AbstractSubOptimizer} # Subproblem optimizer

    options::Parameters

    feasibility_restoration::Bool # indicator for feasibility restoration
    iter::Int # iteration counter
    ret::Int # solution status
    start_time::Float64 # solution start time

    function SlpTR(problem::Model{T,Tv,Tt}) where {T,Tv<:AbstractArray{T},Tt}
        slp = new{T,Tv,Tt}()
        slp.problem = problem
        slp.x = Tv(undef, problem.n)
        slp.p = zeros(problem.n)
        slp.p_slack = Dict()
        slp.lambda = zeros(problem.m)
        slp.mult_x_L = zeros(problem.n)
        slp.mult_x_U = zeros(problem.n)
        slp.df = Tv(undef, problem.n)
        slp.E = Tv(undef, problem.m)
        slp.dE = Tv(undef, length(problem.j_str))
        slp.phi = Inf

        slp.ν = Tv(undef, problem.m)
        slp.Δ = problem.parameters.tr_size
        slp.Δ_max = 2.0
        slp.alpha1 = 0.1
        slp.alpha2 = 0.25

        slp.prim_infeas = Inf
        slp.dual_infeas = Inf
        slp.compl = Inf

        slp.options = problem.parameters
        slp.optimizer = nothing

        slp.feasibility_restoration = false
        slp.iter = 1
        slp.ret = -5
        slp.start_time = 0.0
        return slp
    end
end

"""
    run!

Run the trust-region SLP algorithm
"""
function run!(slp::SlpTR)

    slp.start_time = time()

    if slp.options.OutputFlag == 1
        sparsity_val = ifelse(
            slp.problem.m > 0,
            length(slp.problem.j_str) / (slp.problem.m * slp.problem.n),
            0.0,
        )
        @printf("LP subproblem sparsity: %e\n", sparsity_val)
        add_statistics(slp.problem, "sparsity", sparsity_val)
    else
        Logging.disable_logging(Logging.Info)
    end

    # Set initial point from MOI
    @assert length(slp.x) == length(slp.problem.x)
    slp.x .= slp.problem.x
    # Adjust the initial point to satisfy the column bounds
    for i = 1:slp.problem.n
        if slp.problem.x_L[i] > -Inf
            slp.x[i] = max(slp.x[i], slp.problem.x_L[i])
        end
        if slp.problem.x_U[i] > -Inf
            slp.x[i] = min(slp.x[i], slp.problem.x_U[i])
        end
    end

    slp.iter = 1
    while true

        # evaluate function, constraints, gradient, Jacobian
        eval_functions!(slp)

        LP_time_start = time()
        # solve LP subproblem (to initialize dual multipliers)
        slp.p, slp.lambda, slp.mult_x_U, slp.mult_x_L, slp.p_slack, status =
            sub_optimize!(slp, slp.Δ)
        # @show slp.x, slp.p, slp.lambda, slp.p_slack

        add_statistics(slp.problem, "LP_time", time() - LP_time_start)

        if status ∉ [MOI.OPTIMAL, MOI.INFEASIBLE]
            @warn("Unexpected LP subproblem solution status")
            slp.ret == -3
            if norm_violations(slp, slp.x) <= slp.options.tol_infeas
                slp.ret = 6
            end
            break
        elseif status == MOI.INFEASIBLE
            if slp.feasibility_restoration == true
                @info "Failed to find a feasible direction"
                if slp.prim_infeas <= slp.options.tol_infeas
                    slp.ret = 6
                else
                    slp.ret = 2
                end
                break
            else
                slp.feasibility_restoration = true
                continue
            end
        end

        compute_nu!(slp)

        slp.prim_infeas = norm_violations(slp, Inf)
        slp.dual_infeas = KT_residuals(slp)
        slp.compl = norm_complementarity(slp)

        print_header(slp)
        print(slp)
        collect_statistics(slp)

        # Check the first-order optimality condition
        if slp.prim_infeas <= slp.options.tol_infeas &&
           slp.compl <= slp.options.tol_residual &&
           norm(slp.p, Inf) <= slp.options.tol_direction

            if slp.feasibility_restoration
                slp.feasibility_restoration = false
                slp.iter += 1
                continue
            elseif slp.dual_infeas <= slp.options.tol_residual
                slp.ret = 0
                break
            end
        end

        # Iteration counter limit
        if slp.iter >= slp.options.max_iter
            slp.ret = -1
            if slp.prim_infeas <= slp.options.tol_infeas
                slp.ret = 6
            end
            break
        end

        # step size computation
        rho = step_quality(slp)
        if slp.ret ∈ [0, 2, 6]
            break
        end
        if rho >= 0
            # update primal points
            slp.x += slp.p
        end

        slp.iter += 1
    end
    slp.problem.obj_val = slp.problem.eval_f(slp.x)
    slp.problem.status = Int(slp.ret)
    slp.problem.x .= slp.x
    slp.problem.g .= slp.E
    slp.problem.mult_g .= slp.lambda
    slp.problem.mult_x_U .= slp.mult_x_U
    slp.problem.mult_x_L .= slp.mult_x_L
    add_statistic(slp.problem, "iter", slp.iter)
end

"""
    step_quality

Evaluate the next step and adjust the trust region size
"""
function step_quality(slp::SlpTR)
    slp.phi = compute_phi(slp, slp.x, slp.p) - compute_phi(slp, slp.x)
    # @show compute_phi(slp, slp.x, slp.p), compute_phi(slp, slp.x)

    ϕ_pre = compute_derivative(slp)

    ρ = 0.0
    # @show slp.phi, ϕ_pre
    if abs(ϕ_pre) > 0.0
        ρ = slp.phi / ϕ_pre
        if ρ <= 0
            slp.Δ *= slp.alpha1
        elseif ρ <= 0.25
            slp.Δ *= slp.alpha2
        elseif ρ > 0.75
            slp.Δ = min(2 * slp.Δ, slp.Δ_max)
        end
    else
        ρ = -slp.phi
        if abs(slp.phi) < 1.e-8
            if slp.feasibility_restoration
                slp.feasibility_restoration = false
            else
                if slp.prim_infeas <= slp.options.tol_infeas
                    if slp.dual_infeas <= slp.options.tol_residual &&
                       slp.compl <= slp.options.tol_residual
                        slp.ret = 0
                    else
                        slp.ret = 6
                    end
                else
                    slp.ret = 2
                end
            end
        end
    end

    return ρ
end

"""
    compute_phi

Evaluate and return the merit function value for the current point x.
"""
compute_phi(slp::SlpTR, x::Tv) where {T,Tv<:AbstractArray{T}} =
    compute_phi(slp, x, 0.0, slp.p)

"""
    compute_phi

Evaluate and return the merit function value for the next point x + p.
"""
compute_phi(slp::SlpTR, x::Tv, p::Tv) where {T,Tv<:AbstractArray{T}} =
    compute_phi(slp, x, 1.0, p)

"""
    print_header

Print the header of iteration information.
"""
function print_header(slp::SlpTR)
    if slp.options.OutputFlag == 0
        return
    end
    if (slp.iter - 1) % 25 == 0
        @printf("  %6s", "iter")
        @printf("  %15s", "f(x_k)")
        @printf("  %15s", "Δ")
        @printf("  %14s", "|p|")
        @printf("  %14s", "inf_pr")
        @printf("  %14s", "inf_du")
        @printf("  %14s", "compl")
        @printf("  %10s", "time")
        @printf("\n")
    end
end

"""
    print

Print iteration information.
"""
function print(slp::SlpTR)
    if slp.options.OutputFlag == 0
        return
    end
    st = ifelse(slp.feasibility_restoration, "FR", "  ")
    @printf("%2s%6d", st, slp.iter)
    @printf("  %+6.8e", slp.f)
    @printf("  %+6.8e", slp.Δ)
    @printf("  %6.8e", norm(slp.p, Inf))
    @printf("  %6.8e", slp.prim_infeas)
    @printf("  %6.8e", slp.dual_infeas)
    @printf("  %6.8e", slp.compl)
    @printf("  %10.2f", time() - slp.start_time)
    @printf("\n")
end

"""
    collect_statistics

Collect iteration information.
"""
function collect_statistics(slp::SlpTR)
    if slp.options.StatisticsFlag == 0
        return
    end
    add_statistics(slp.problem, "f(x)", slp.f)
    add_statistics(slp.problem, "|p|", norm(slp.p, Inf))
    add_statistics(slp.problem, "|J|2", norm(slp.dE, 2))
    add_statistics(slp.problem, "|J|inf", norm(slp.dE, Inf))
    add_statistics(slp.problem, "inf_pr", slp.prim_infeas)
    # add_statistics(slp.problem, "inf_du", dual_infeas)
    add_statistics(slp.problem, "compl", slp.compl)
    add_statistics(slp.problem, "time_elapsed", time() - slp.start_time)
end

max_issues_repo_name: exanauts/ActiveSetMethods