text:
# evaluation of objective function

"""
# TermFunction{F}  with F<: Function

This structure stores a function that compute a term of serie representing the objective function given
the position and the term of the serie to evaluate.

This **eval** function can be a closure caching observations data (see regressions examples)

## Fields

- eval is a function of signature   **Fn(Array{Float64, N}, Int64)::Float64**
    taking as arguments (in this order):
    - position
    - a term rank

    This function returns the value of the component of rank term in the sum of objective function
        at position position given in args

- nbterms in the sum

- dims : characterize the dimensions on variable for which we do a minimization
    It is not the same as observations which is always a one dimensional array.

    For example on a logistic regression with nbclass we used (Cf applis) dims = (d , nbclass-1)
    with *d* the length(observations) and nbclass -1 beccause of solvability constraints.

    It is also the dimension of gradients we compute.

Note:
    When instantiating a TermFunction you just write:
    *tf = TermFunction{typeof(yourfunction)}(yourfunction, nbterms, dims)*

"""
mutable struct TermFunction{F <: Function}
    eval :: F
    nbterms :: Int64
    dim :: Dims
    function TermFunction{F}(evalf :: F, nbterms:: Int64, d::Dims) where{F}
        # check signature
        z = zeros(d)
        evalf(z, 1)
        # we can initialize object
        new(evalf, nbterms, d)
    end
end

"""
#  function compute_value(tf :: TermFunction{F}, position :: Vector{Float64})

This function computes value of function at a given position summing over all terms

"""
function compute_value(tf :: TermFunction{F}, position :: Array{Float64,N}) where {F,N}
    nbterm = tf.nbterms
    value = compute_value(tf, position, Vector{Int64}((1:nbterm)))
    value
end

"""
#  function compute_value(tf :: TermFunction{F}, position :: Array{Float64,N}, terms::Vector{Int64})

This function computes value of function at a given position summing over all terms passed as arg

"""
function compute_value(tf :: TermFunction{F}, position :: Array{Float64,N}, terms::Vector{Int64}) where {F,N}
    nbterms = length(terms)
    # split in blocks
    batchsize = 2000
    nbblocks = floor(Int64, nbterms / batchsize)
    nbblocks = nbterms % batchsize > 0  ? nbblocks+1 : nbblocks
    blockvalue = zeros(Float64, nbblocks)
    #
    Threads.@threads for i in 1:nbblocks
        first = (i-1) * batchsize +1
        last = min(i*batchsize, nbterms)
        blockvalue[i] = mapreduce(i->tf.eval(position, i), + , terms[first:last])
    end
    # pamp with batch size = 1000. check speed versus a mapreduce
    value = sum(blockvalue)/nbterms
    value
end

"""
# struct TermGradient{G}  with G <: Function

This structure is the building block for all gradient computations.
It stores a **eval** function that computes the gradient of term of serie representing the objective function given
the position and the term of the serie to evaluate.

This function **eval** can be a closure caching observations data (see regressions examples)

## Fields

- eval is a function of signature **Fn(Array{Float64,N}, Int64, Array{Float64,N})**
    taking as arguments (in this order):
    - position
    - a term rank
    - a vector for gradient to be returned in. This avoid reallocations as we loop on this function summing
            the resulting Array on required terms.
            Array for gradient has the same dimension as array for position (...)

- nbterms : the number of terms in the sum defining the objective function
- dims : characterize the dimensions on variable for which we do a minimization

"""
struct TermGradient{G<:Function}
    eval :: G
    nbterms :: Int64
    dim :: Dims
    function TermGradient{G}(evalg :: G, nbterms :: Int64, d::Dims) where {G}
        # we check that signature of function defined by the user is correct
        z = zeros(d)
        grad = zeros(d)
        evalg(z, 1, grad)
        # we can allocate our object
        new(evalg, nbterms, d)
    end
end

"""
# function compute_gradient!(termg::TermGradient{G}, position::Array{Float64,N} , term:: Int64, gradient:: Array{Float64,N}) where {G,N}

the function  dispatches to TermGradient  the actual gradient computation for a term.

"""
function compute_gradient!(termg::TermGradient{G}, position::Array{Float64,N} , term:: Int64, gradient:: Array{Float64,N}) where {G,N}
    termg.eval(position, term, gradient)
end

"""
# function compute_gradient!(termg::TermGradient{G}, position :: Array{Float64,N}, terms::Vector{Int64}, gradient::Array{Float64,N}) where {G,N}

This function computes a gradient Array at a given position summing over all terms passed as arg.

NOTA: This function computes the mean of gradients returned by the termg function returned on each term.
So that the gradient computed on a batch is an estimator of the gradient computed on the whole objective function.

It is multithreaded and computes gradient by blocks of 1500 terms.

"""
function compute_gradient!(termg::TermGradient{G}, position :: Array{Float64,N}, terms::Vector{Int64}, gradient::Array{Float64,N}) where {G,N}
        #
        batchsize=1500
        nbterms = length(terms)
         # split in blocks
        nbblocks = floor(Int64, nbterms / batchsize)
        nbblocks = nbterms % batchsize > 0  ? nbblocks+1 : nbblocks
        # we must allocate an array of nbblocks arrays of dimensions size(gradient)
        dimg = ndims(gradient)
        gradblocks = Vector{Array{Float64,dimg}}(undef, nbblocks)
        # CAVEAT to be threaded
        Threads.@threads for i in 1:nbblocks
            first = (i-1) * batchsize +1
            last = min(i*batchsize, nbterms)
            # must respect dimensions of gradient
            gradtmp = zeros(Float64, size(gradient))
            gradblocks[i] = zeros(Float64, size(gradient))
            for j in first:last
                termg.eval(position, terms[j], gradtmp)
                gradblocks[i] = gradblocks[i] + gradtmp
            end
        end
        # recall that in julia is column oriented so summing along rows is sum(,dims=2)
        copy!(gradient, sum(gradblocks)/nbterms)
 #       @debug "gradient sum blocks" gradient
end

"""
# function compute_gradient!(termg::TermGradient{G}, position :: Array{Float64,N}, gradient::Array{Float64,N}) where {G,N}

This function computes a gradient Array at a given position summing over all terms
"""
function compute_gradient!(termg::TermGradient{G}, position :: Array{Float64,N}, gradient::Array{Float64,N}) where {G,N}
    nbterms = termg.nbterms
    compute_gradient!(termg, position, collect(1:nbterms), gradient)
end

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

"""
#  Evaluator{F,G}

   This structure contains all that is necessary to compute function value at any term
   and gradient total or any partial sum of terms.
   It can be passed as argument in any algorithm for our stochastic gradients.

## Fields

- compute\\_term\\_value : a TermFunction{F}
- compute\\_term\\_gradient: a TermGradient{G}

It is associated to various functions dispatching computations to it.

"""
struct Evaluator{F,G}
    #  A vector of observation , associated value
    compute_term_value :: TermFunction{F}
    compute_term_gradient :: TermGradient{G}
    function Evaluator{F,G}(compute_term_value :: TermFunction{F}, compute_term_gradient :: TermGradient{G}) where {F,G}
        new(compute_term_value, compute_term_gradient)
    end
end

"""
# function compute_gradient!(evaluator::Evaluator{F,G}, position :: Array{Float64,N}, term::Int64 , gradient :: Array{Float64,N}) where {F,G,N}

the function computes gradient given an evaluator, a position, and a term

"""
function compute_gradient!(evaluator::Evaluator{F,G}, position :: Array{Float64,N}, term::Int64 , gradient :: Array{Float64,N}) where {F,G,N}
    compute_gradient!(evaluator.compute_term_gradient, position, term, gradient)
end

"""
# function compute_gradient!(evaluator::Evaluator{F,G}, position :: Array{Float64,N}, terms::Vector{Int64}, gradient :: Array{Float64,N}) where {F,G,N}

the function computes a gradient given an evaluator, a position, and a vector of rank term

"""
function compute_gradient!(evaluator::Evaluator{F,G}, position :: Array{Float64,N}, terms::Vector{Int64}, gradient :: Array{Float64,N}) where {F,G,N}
    compute_gradient!(evaluator.compute_term_gradient, position, terms, gradient)
end

"""
# function compute_gradient!(evaluator::Evaluator{F,G}, position :: Array{Float64,N}, gradient::Array{Float64,N}) where {F, G,N}

This function computes a gradient Array at a given position summing over all terms
"""
function compute_gradient!(evaluator::Evaluator{F,G} , position :: Array{Float64,N}, gradient::Array{Float64,N}) where {F, G,N}
    compute_gradient!(evaluator.compute_term_gradient, position, gradient)
end

"""
# function compute_value(evaluator::Evaluator{F,G}, position :: Array{Float64,N}) where {F,G,N}

this function computes a value given an evaluator, and a position using all terms
"""
function compute_value(evaluator::Evaluator{F,G}, position :: Array{Float64,N}) where {F,G,N}
    compute_value(evaluator.compute_term_value, position)
end

"""
# function compute_value(evaluator::Evaluator{F,G}, position :: Array{Float64,N}, terms::Vector{Int64})  where {F,G,N}

this function computes a value given an evaluator, and a position and a list of terms

"""
function compute_value(evaluator::Evaluator{F,G}, position :: Array{Float64,N}, terms::Vector{Int64})  where {F,G,N}
    compute_value(evaluator.compute_term_value, position, terms)
end

"""
# function get_nbterms(evaluator::Evaluator{F,G}) where {F,G}

retrieves number of terms in the sum defining objective function
"""
function get_nbterms(evaluator::Evaluator{F,G}) where {F,G}
    evaluator.compute_term_gradient.nbterms
end


text:
# To run this script, `cd` to the `./test/fixtures` directory and then, from the Julia terminal, `include("./runner.jl")`.

import JSON

function gen( x, name )
	y = Array( Float64, length( x ) );
	for i in eachindex(x)
		y[i] = erfinv( x[i] );
	end

	data = Dict([
		("x", x),
		("expected", y)
	]);

	outfile = open( name, "w" );
	write( outfile, JSON.json(data) );
	close( outfile );
end

# |x| <= 0.5
x = linspace( -0.5, 0.5, 3000 );
gen( x, "x_-0.5_0.5.json" );

# 1-|x| > 0.25
x = linspace( 0.5, 0.75, 500 );
gen( x, "x_0.5_0.75.json" );

x = linspace( -0.5, -0.75, 500 );
gen( x, "x_-0.5_-0.75.json" );

# 0.75 < |x| < 0.9998765901959134
x = linspace( 0.75, 0.9998, 500 );
gen( x, "x_0.75_0.9998.json" );

# 0.9998765901959134 < |x| < 0.9999999999999998
x = linspace( 0.9998, 0.9999999999999998, 500 );
gen( x, "x_0.9998_0.9999..8.json" );

# 0.9999999999999998 < |x| < 1
x = linspace( 0.9999999999999998, 1, 500 );
gen( x, "x_0.9999..8_1.json" );

ext: jl

max_stars_repo_name: math-io/erfinv

max_issues_repo_name: math-io/erfinv

max_forks_repo_name: math-io/erfinv


text:
using LinearAlgebra, MixedModels, StatsModels, Random, Test

xtx(X) = X'X  # creat the symmetric matrix X'X from X

const rng = Random.MersenneTwister(4321234)

const simdat = (
    G = repeat('A':'T', inner=10),
    H = repeat('a':'e', inner=2, outer=20),
    U = repeat(0.:9, outer=20),
    V = repeat(-4.5:4.5, outer=20),
    Y = 0.1 * randn(rng, 200),
    Z = rand(rng, 200)
)

@testset "fullranknumeric" begin
    XtX = xtx(modelmatrix(@formula(Y ~ 1 + U), simdat))
    ch = statscholesky(Symmetric(XtX, :U))
    @test ch.rank == 2
    @test ch.piv == 1:2
    @test iszero(ch.info)
    @test isapprox(xtx(ch.U), XtX[ch.piv, ch.piv])
end

@testset "fullrankcategorical" begin
    XtX = xtx(modelmatrix(@formula(Y ~ 1 + G*H), simdat))
    ch = statscholesky(Symmetric(XtX, :U))
    @test ch.rank == 100
    @test ch.piv == 1:100
    @test iszero(ch.info)
    @test isapprox(xtx(ch.U), XtX)
end

@testset "dependentcolumn" begin
    XtX = xtx(modelmatrix(@formula(Y ~ 1 + U + V + Z), simdat))
    ch = statscholesky(Symmetric(XtX, :U))
    perm = [1,2,4,3]
    @test ch.rank == 3
    @test ch.piv == perm
    @test isapprox(xtx(ch.U), XtX[perm, perm])
end

@testset "missingcells" begin
    XtX = xtx(modelmatrix(@formula(Y ~ 1 + G*H), simdat)[5:end,:])
    ch = statscholesky(Symmetric(XtX, :U))
    perm = [1:42; 44:100; 43]
    @test ch.rank == 98
    @test ch.piv == perm
    @test isapprox(xtx(ch.U), XtX[perm, perm], atol=0.00001)
end


text:
module MvNormalMeanPrecisionNodeTest

using Test
using ReactiveMP
using Random

import ReactiveMP: make_node

@testset "MvNormalMeanPrecisionNode" begin

    @testset "Creation" begin
        node = make_node(MvNormalMeanPrecision)

        @test functionalform(node)          === MvNormalMeanPrecision
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :Λ)
        @test factorisation(node)           === ((1, 2, 3), )
        @test localmarginalnames(node)      === (:out_μ_Λ, )
        @test metadata(node)                === nothing

        node = make_node(MvNormalMeanPrecision, meta = 1)

        @test functionalform(node)          === MvNormalMeanPrecision
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :Λ)
        @test factorisation(node)           === ((1, 2, 3), )
        @test localmarginalnames(node)      === (:out_μ_Λ, )
        @test metadata(node)                === 1

        node = make_node(MvNormalMeanPrecision, factorisation = ((1, ), (2, 3)))

        @test functionalform(node)          === MvNormalMeanPrecision
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :Λ)
        @test factorisation(node)           === ((1, ), (2, 3, ))
        @test localmarginalnames(node)      === (:out, :μ_Λ, )
        @test metadata(node)                === nothing
    end

    @testset "AverageEnergy" begin

        begin
            q_out = PointMass([ 1.0, 1.0 ])
            q_μ   = MvNormalMeanPrecision([ 1.0, 1.0 ], [ 1.0 0.0; 0.0 1.0 ])
            q_Λ   = Wishart(3, [ 2.0 0.0; 0.0 2.0 ])

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), G in (Wishart, )
                marginals = (Marginal(q_out, false, false), Marginal(convert(N, q_μ), false, false), Marginal(convert(G, q_Λ), false, false))
                @test score(AverageEnergy(), MvNormalMeanPrecision, Val{(:out, :μ, :Λ)}, marginals, nothing) ≈ 6.721945550750932
            end
        end

        begin
            q_out = PointMass([0.3873049736301686, 0.6250550344669505, 0.9865681758036855])
            q_μ   = MvNormalMeanPrecision([0.21393068431529905, 0.5271266845217295, 0.0032162989237940476], [1.1178940049080035 1.004643185289189 0.5726053202206195; 1.004643185289189 0.9113807582278761 0.515956570765438; 0.5726053202206195 0.515956570765438 0.3402071079911484])
            q_Λ   = Wishart(4, [0.4720138588203935 0.49760446060666863 0.4405032929933653; 0.49760446060666863 0.6757755009955487 0.636955076077564; 0.4405032929933653 0.636955076077564 0.6660144475917664])

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), G in (Wishart, )
                marginals = (Marginal(q_out, false, false), Marginal(convert(N, q_μ), false, false), Marginal(convert(G, q_Λ), false, false))
                @test score(AverageEnergy(), MvNormalMeanPrecision, Val{(:out, :μ, :Λ)}, marginals, nothing) ≈ 57.362404928342
            end
        end

        begin
            q_out = MvNormalMeanPrecision([0.8625589157256603, 0.6694783342639599], [1.0014322413749484 0.7989099036521625; 0.7989099036521625 1.0976639268696966])
            q_μ   = MvNormalMeanPrecision([0.9334416739853251, 0.38318522701093105], [0.21867945696266933 0.5704895781120056; 0.5704895781120056 1.5321190185800933])
            q_Λ   = Wishart(3, [1.2509510680086597 1.3662166757006484; 1.3662166757006484 1.4925451500802907])

            for N1 in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), N2 in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), G in (Wishart, )
                marginals = (Marginal(convert(N1, q_out), false, false), Marginal(convert(N2, q_μ), false, false), Marginal(convert(G, q_Λ), false, false))
                @test score(AverageEnergy(), MvNormalMeanPrecision, Val{(:out, :μ, :Λ)}, marginals, nothing) ≈ 114.57678973411974
            end
        end

        begin
            q_out = MvNormalMeanPrecision([0.9580573744420284, 0.26086767384401943, 0.28770712080127914], [0.9182022104547379 0.6106281139047683 1.0779505969017997; 0.6106281139047683 1.032830918722213 0.9771172696798327; 1.0779505969017997 0.9771172696798327 1.5109842607623276])
            q_μ   = MvNormalMeanPrecision([0.7234472062526713, 0.2824687847180598, 0.2746840635757333], [0.3590359808189558 0.3331759830932437 0.7133660103700628; 0.3331759830932437 0.9384711630986695 0.972435803354717; 0.7133660103700628 0.972435803354717 1.8219531599263405])
            q_Λ   = Wishart(4, [0.32860287160871043 0.11075216999194112 0.19071629144893898; 0.11075216999194112 0.15595676288444096 0.1700073984759777; 0.19071629144893898 0.1700073984759777 0.2755688590223264])

            for N1 in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), N2 in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), G in (Wishart, )
                marginals = (Marginal(convert(N1, q_out), false, false), Marginal(convert(N2, q_μ), false, false), Marginal(convert(G, q_Λ), false, false))
                @test score(AverageEnergy(), MvNormalMeanPrecision, Val{(:out, :μ, :Λ)}, marginals, nothing) ≈ 14.92006071729396
            end
        end

        begin
            q_out = MvNormalMeanPrecision([0.8436742420779029, 0.10538074240979411, 0.6165907670265702], [1.9016529880577202 0.10379920228081287 1.8123623785462684; 0.10379920228081287 0.608931528715827 0.31789777698874117; 1.8123623785462684 0.31789777698874117 1.812858220779861])
            q_μ   = MvNormalMeanPrecision([0.2069541691083061, 0.5589468611740172, 0.4820571901630184], [0.4530148040808207 0.3039331260702957 0.35941114218805115; 0.3039331260702957 0.4389148552241666 0.23999683261779386; 0.35941114218805115 0.23999683261779386 0.3727021595284692])
            q_Λ   = PointMass([0.705864402777104 1.0323593788454744 0.9684301882026219; 1.0323593788454744 1.7349617392942005 1.7059080462235048; 0.9684301882026219 1.7059080462235048 1.816579920226959])

            for N1 in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), N2 in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision)
                marginals = (Marginal(convert(N1, q_out), false, false), Marginal(convert(N2, q_μ), false, false), Marginal(q_Λ, false, false))
                @test score(AverageEnergy(), MvNormalMeanPrecision, Val{(:out, :μ, :Λ)}, marginals, nothing) ≈ 38.31309509463591
            end
        end

        begin
            q_out_μ = MvNormalMeanPrecision([0.2932046487282065, 0.7716085147100042, 0.03978072440454361, 0.2814883836121471], [1.0684808331872628 0.6721958601342372 0.7164104160110533 0.2869444930570181; 0.6721958601342372 1.4472786772965438 1.3516272546828385 0.5533932426057602; 0.7164104160110533 1.3516272546828385 1.2958919150214063 0.5171784879755998; 0.2869444930570181 0.5533932426057602 0.5171784879755998 0.30898006058959576])
            q_Λ   = PointMass([0.6531636197654624 0.8023009003361893; 0.8023009003361893 1.0001231468452463])

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision)
                marginals = (Marginal(convert(N, q_out_μ), false, false), Marginal(q_Λ, false, false))
                @test score(AverageEnergy(), MvNormalMeanPrecision, Val{(:out_μ, :Λ)}, marginals, nothing) ≈ 87.02279699224562
            end
        end

        begin
            q_out_μ = MvNormalMeanPrecision([0.35156676223859784, 0.6798203100143094, 0.954485919235333, 0.9236981452828203], [1.3182839156957349 0.9159049032047119 1.170482409249098 1.132202025059748; 0.9159049032047119 1.4737964254194567 1.4024254322343757 0.7350293025705011; 1.170482409249098 1.4024254322343757 2.0577570913647105 1.3137472032115916; 1.132202025059748 0.7350293025705011 1.3137472032115916 1.2083880803032556])
            q_Λ   = PointMass([0.6202203324986401 0.6037236520027125; 0.6037236520027125 1.140352058799863])

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision)
                marginals = (Marginal(convert(N, q_out_μ), false, false), Marginal(q_Λ, false, false))
                @test score(AverageEnergy(), MvNormalMeanPrecision, Val{(:out_μ, :Λ)}, marginals, nothing) ≈ 11.240486862933556
            end
        end

        begin
            q_out_μ = MvNormalMeanPrecision([0.35156676223859784, 0.6798203100143094, 0.954485919235333, 0.9236981452828203], [1.3182839156957349 0.9159049032047119 1.170482409249098 1.132202025059748; 0.9159049032047119 1.4737964254194567 1.4024254322343757 0.7350293025705011; 1.170482409249098 1.4024254322343757 2.0577570913647105 1.3137472032115916; 1.132202025059748 0.7350293025705011 1.3137472032115916 1.2083880803032556])
            q_Λ   = Wishart(3, [0.6202203324986401 0.6037236520027125; 0.6037236520027125 1.140352058799863])

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), G in (Wishart, )
                marginals = (Marginal(convert(N, q_out_μ), false, false), Marginal(convert(G, q_Λ), false, false))
                @test score(AverageEnergy(), MvNormalMeanPrecision, Val{(:out_μ, :Λ)}, marginals, nothing) ≈ 28.552276936591902
            end
        end

        begin
            q_out_μ = MvNormalMeanPrecision([0.5500318675144917, 0.817183383468487, 0.5347013723720735, 0.5114047207869656], [1.1343108996797893 0.47265813998866796 0.9272779005094741 1.0428711245585718; 0.47265813998866796 0.715974520938224 0.7635970006454071 1.193450872579602; 0.9272779005094741 0.7635970006454071 1.392689362184066 1.5818118006093156; 1.0428711245585718 1.193450872579602 1.5818118006093156 2.158417891827849])
            q_Λ   = Wishart(3, [1.545789934380036 1.1984283177361577; 1.1984283177361577 1.1061038557496674])

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), G in (Wishart, )
                marginals = (Marginal(convert(N, q_out_μ), false, false), Marginal(convert(G, q_Λ), false, false))
                @test score(AverageEnergy(), MvNormalMeanPrecision, Val{(:out_μ, :Λ)}, marginals, nothing) ≈ 1938.6631975673586
            end
        end

    end

end
end


text:
    ######### getting all together in Housedorff
    # julia -g 2 C:\GitHub\GitHub\NuclearMedEval\test\distanceMetr\HousdorffTest\integrationHousedorff.jl
    using Revise, Parameters, Logging, Test
    using CUDA
    includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
    using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
    using ..MainLoopKernel,..PrepareArrtoBool,..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils,..ResultListUtils, ..Housdorff

    CUDA.allowscalar(true)

    mainArrDims= (120,120,120);
    mainArrCPU= zeros(UInt8,mainArrDims);
    refArrCPU = zeros(UInt8,mainArrDims);
    ##### we will create two planes 20 units apart from each
    mainArrCPU[10:50,10:50,10].= 1;
    refArrCPU[10:50,10:50,80].= 1;
    # mainArrCPU[10:12,10:12,10].= 1;
    # refArrCPU[10:12,10:12,30].= 1;

    goldGPU = CuArray(mainArrCPU);
    segmGPU= CuArray(mainArrCPU);

    robustnessPercent= 0.95
    numberToLooFor=1
    boolKernelArgs, mainKernelArgs,threadsBoolKern,blocksBoolKern ,threadsMainKern,blocksMainKern ,shmemSizeBool,shmemSizeMain=    preparehousedorfKernel(goldGPU,segmGPU,robustnessPercent,numberToLooFor)
    mainArrDims,dataBdim,metaData,metaDataDims,reducedGoldA,reducedSegmA,loopXinPlane,loopYinPlane,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp ,numberToLooFor,inBlockLoopXZIterWithPadding,shmemblockDataLoop,shmemblockDataLenght,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength=boolKernelArgs

        gold3d = CuArray(mainArrCPU)
        segm3d= CuArray(refArrCPU)
        numberToLooFor =UInt8(1)
        robustnessPercent = 0.9

    # globalIterationNumb= Housdorff.getHousedorffDistance(goldGPU,segmGPU,boolKernelArgs,mainKernelArgs,threadsBoolKern,blocksBoolKern ,threadsMainKern,blocksMainKern,shmemSizeBool,shmemSizeMain)
    # globalIterationNumb[1]

        boolKernelArgs, mainKernelArgs,threadsBoolKern,blocksBoolKern ,threadsMainKern,blocksMainKern ,shmemSizeBool,shmemSizeMain=    preparehousedorfKernel(gold3d,segm3d,robustnessPercent,numberToLooFor)
        mainArrDims,dataBdim,metaData,metaDataDims,reducedGoldA,reducedSegmA,loopXinPlane,loopYinPlane,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp ,numberToLooFor,inBlockLoopXZIterWithPadding,shmemblockDataLoop,shmemblockDataLenght,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength=boolKernelArgs

        golddd = CuArray(mainArrCPU);
        segmmm= CuArray(refArrCPU);

    @cuda threads=threadsBoolKern blocks=blocksBoolKern shmem=shmemSizeBool  cooperative=true boolKernelLoad(golddd,segmmm,boolKernelArgs...)

    Int64(sum(mainArrCPU))
    Int64(fn[1])
    Int64(fp[1])
    Int64(minxRes[1])
    Int64(maxxRes[1])
    Int64( minyRes[1])
    Int64(maxyRes[1])
    Int64( minzRes[1])
    Int64(maxzRes[1])
size(reducedGoldA)
size(reducedSegmA)
size(metaData)
    #now time to get data structures dependent on bool kernel like for example loading subsections of meta data, creating work queue ...
    #some arrays needs to be instantiated only after we know the number of the false and true positives
 referenceArrs,dilatationArrs, mainArrDims,dataBdim,numberToLooFor,metaDataDims,metaDataNull,iterThrougWarNumb,robustnessPercent,shmemSumLengthMaxDiv4,globalFpResOffsetCounter,globalFnResOffsetCounter    ,globalIterationNumber,globalCurrentFnCount,globalCurrentFpCount    ,globalIterationNumb,workQueue,workQueueCounter,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed    ,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop    ,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength, fn,fp,resList,inBlockLoopXZIterWithPadding,paddingStore,shmemblockDataLenght,shmemblockDataLoop= mainKernelArgs
        size(metaData)

    metaData,reducedGoldA  ,reducedSegmA ,paddingStore,resList,workQueue,workQueueCounter= getBigGPUForHousedorffAfterBoolKernel(metaData,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp,reducedGoldA,reducedSegmA,dataBdim)

    reducedSegmA[:,:,3]

    metaDataDims = size(metaData)
    metaDataLength= metaDataDims[1]*metaDataDims[2]*metaDataDims[3]
    loopMeta= fld(metaDataLength,blocksMainKern )
    loopWarpMeta= cld(metaDataLength,(blocksMainKern*threadsMainKern[1] ))

    dilatationArrs= (reducedGoldA,reducedSegmA)
    #reverse order as when dilatating gold we want to establish do we cover segm voxels
    referenceArrs=(golddd,segmmm )
    shmemblockDataLenght = threadsMainKern[1]*threadsMainKern[2]*2
    shmemblockDataLoop = fld(shmemblockDataLenght,threadsMainKern[1]*threadsMainKern[2])

using  BenchmarkTools
    BenchmarkTools.DEFAULT_PARAMETERS.samples = 2
    BenchmarkTools.DEFAULT_PARAMETERS.seconds =600
    BenchmarkTools.DEFAULT_PARAMETERS.gcsample = true

    #main calculations

    # @benchmark CUDA.@sync @cuda threads=threadsMainKern blocks=blocksMainKern shmem=shmemSizeMain cooperative=true mainKernelLoad( referenceArrs,dilatationArrs, mainArrDims,dataBdim
    @benchmark CUDA.@sync @cuda threads=threadsMainKern blocks=blocksMainKern shmem=shmemSizeMain cooperative=true mainKernelLoad( referenceArrs,dilatationArrs, mainArrDims,dataBdim
    ,numberToLooFor,metaDataDims,metaData,iterThrougWarNumb,robustnessPercent
    ,shmemSumLengthMaxDiv4,globalFpResOffsetCounter,globalFnResOffsetCounter
    ,globalIterationNumber,globalCurrentFnCount,globalCurrentFpCount
    ,globalIterationNumb,workQueue,workQueueCounter
    ,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed
    ,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY
    ,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop
    ,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength, fn,fp,resList,inBlockLoopXZIterWithPadding,paddingStore,shmemblockDataLenght,shmemblockDataLoop)

    globalIterationNumb[1]

    Int64(workQueueCounter[1])
    Int64.(workQueue)

#     (threadIdxX())+(innerWarpNumb+21)*33
#     32+(16+21)*33
# dataBdim[1]* dataBdim[2]*2

resList

for i in 1:16
    # println(" i $(i)  $((Int64(metaData[2,2,2, (getOldCountersBeg()+1)])))")
    println(" i $(i)  $((Int64(metaData[2,2,2, (getNewCountersBeg()+1)])))")
end
mainArrDims


text:
export Op

struct Op; f; args; kwargs; end

Op(f,x...;o...) = Op(f, x, o)

(o::Op)(x) = o.f(x, o.args...; o.kwargs...)


text:
@testset "Create types" begin

    GB = 10
    PER = 20
    baseindex = 100.0

    v = rand(PER, GB) .- 0.25
    ipc = capitalize(v, baseindex)
    w = rand(GB)

    basedate = Date(2010,12)
    enddate = basedate + Month(PER-1)
    dates = basedate:Month(1):enddate

    ## Create individual types with same arrays: FullCPIBase, VarCPIBase, IndexCPIBase
    fullcpi = FullCPIBase(ipc, v, w, dates, baseindex)
    vcpi = VarCPIBase(v, w, dates, baseindex)
    indexcpi = IndexCPIBase(ipc, w, dates, baseindex)

    # Test they share the same arrays
    @test fullcpi.ipc === indexcpi.ipc
    @test fullcpi.v === vcpi.v
    @test fullcpi.w === vcpi.w === indexcpi.w

    # Create copies with constructors from FullCPIBase and test they hold different arrays
    vcpi2 = VarCPIBase(fullcpi)
    indexcpi2 = IndexCPIBase(fullcpi)

    @test vcpi.v == vcpi2.v     # same values in matrix
    @test !(vcpi.v === vcpi2.v) # different reference in memory

    @test indexcpi.ipc == indexcpi2.ipc
    @test !(indexcpi.ipc === indexcpi2.ipc)

    @test vcpi.dates == vcpi2.dates
    @test indexcpi.dates == indexcpi2.dates

    ## Test convert methods with different precision floats
    fullcpi32 = convert(Float32, deepcopy(fullcpi))
    vcpi32 = convert(Float32, deepcopy(vcpi))
    indexcpi32 = convert(Float32, deepcopy(indexcpi))

    # as matrixs are different type, == does not apply, check with isapprox
    @test fullcpi32.ipc ≈ fullcpi.ipc
    @test fullcpi32.v == vcpi32.v

    @test indexcpi32.ipc ≈ indexcpi.ipc
    @test fullcpi32.ipc == indexcpi32.ipc

    @test vcpi32.v ≈ vcpi.v
    @test fullcpi32.baseindex == vcpi32.baseindex == indexcpi32.baseindex

    ## Create types with different base indexes
    baseindex = rand(100:0.5:110, GB)

    fullcpi_b = FullCPIBase(ipc, v, w, dates, baseindex)
    vcpi_b = VarCPIBase(fullcpi_b)
    indexcpi_b = IndexCPIBase(vcpi_b)

    # Test they hold the same base indexes, but not same arrays
    @test vcpi_b.baseindex == indexcpi_b.baseindex
    @test !(vcpi_b.baseindex === indexcpi_b.baseindex)
    @test length(fullcpi_b.baseindex) > 1
    @test length(vcpi_b.baseindex) > 1
    @test length(indexcpi_b.baseindex) > 1

    ## Create CountryStructure

    # These hold the same v array values (not same array objects) with different base indexes
    cs = MixedCountryStructure((vcpi, vcpi_b))

    @test cs[1].v == cs[2].v
    @test !(cs[1].v === cs[2].v)
    @test typeof(cs[1].baseindex) != typeof(cs[2].baseindex) # different containers
    @test eltype(cs[1].baseindex) == eltype(cs[2].baseindex) # same types
    @test length(cs[1].baseindex) != length(cs[2].baseindex) # different length

    # Check conversion of precision types
    cs32 = convert(Float32, cs)

    @test cs32[1].v ≈ cs[1].v
    @test eltype(cs32[1].baseindex) == Float32

end


text:
### Unary broadcasting

function broadcasted(::DefaultArrayStyle{N}, op, r::AbstractFill{T,N}) where {T,N}
    return Fill(op(getindex_value(r)), size(r))
end

### Binary broadcasting

function broadcasted(::DefaultArrayStyle, op, a::AbstractFill, b::AbstractFill)
    val = op(getindex_value(a), getindex_value(b))
    return Fill(val, broadcast_shape(size(a), size(b)))
end

function _broadcasted_zeros(a, b)
    return Zeros{promote_type(eltype(a), eltype(b))}(broadcast_shape(size(a), size(b)))
end
function _broadcasted_ones(a, b)
    return Ones{promote_type(eltype(a), eltype(b))}(broadcast_shape(size(a), size(b)))
end

broadcasted(::DefaultArrayStyle, ::typeof(+), a::Zeros, b::Zeros) = _broadcasted_zeros(a, b)
broadcasted(::DefaultArrayStyle, ::typeof(+), a::Ones, b::Zeros) = _broadcasted_ones(a, b)
broadcasted(::DefaultArrayStyle, ::typeof(+), a::Zeros, b::Ones) = _broadcasted_ones(a, b)

broadcasted(::DefaultArrayStyle, ::typeof(*), a::Zeros, b::Zeros) = _broadcasted_zeros(a, b)

for op in (:*, :/)
    @eval begin
        broadcasted(::DefaultArrayStyle, ::typeof($op), a::Zeros, b::Ones) = _broadcasted_zeros(a, b)
        broadcasted(::DefaultArrayStyle, ::typeof($op), a::Zeros, b::Fill) = _broadcasted_zeros(a, b)
        broadcasted(::DefaultArrayStyle, ::typeof($op), a::Zeros, b::Number) = _broadcasted_zeros(a, b)
        broadcasted(::DefaultArrayStyle, ::typeof($op), a::Zeros, b::AbstractRange) = _broadcasted_zeros(a, b)
        broadcasted(::DefaultArrayStyle, ::typeof($op), a::Zeros, b::AbstractArray) = _broadcasted_zeros(a, b)
        broadcasted(::DefaultArrayStyle{1}, ::typeof($op), a::Zeros, b::AbstractRange) = _broadcasted_zeros(a, b)
    end
end

for op in (:*, :\)
    @eval begin
        broadcasted(::DefaultArrayStyle, ::typeof($op), a::Ones, b::Zeros) = _broadcasted_zeros(a, b)
        broadcasted(::DefaultArrayStyle, ::typeof($op), a::Fill, b::Zeros) = _broadcasted_zeros(a, b)
        broadcasted(::DefaultArrayStyle, ::typeof($op), a::Number, b::Zeros) = _broadcasted_zeros(a, b)
        broadcasted(::DefaultArrayStyle, ::typeof($op), a::AbstractRange, b::Zeros) = _broadcasted_zeros(a, b)
        broadcasted(::DefaultArrayStyle, ::typeof($op), a::AbstractArray, b::Zeros) = _broadcasted_zeros(a, b)
        broadcasted(::DefaultArrayStyle{1}, ::typeof($op), a::AbstractRange, b::Zeros) = _broadcasted_zeros(a, b)
    end
end

broadcasted(::DefaultArrayStyle, ::typeof(*), a::Ones, b::Ones) = _broadcasted_ones(a, b)
broadcasted(::DefaultArrayStyle, ::typeof(/), a::Ones, b::Ones) = _broadcasted_ones(a, b)
broadcasted(::DefaultArrayStyle, ::typeof(\), a::Ones, b::Ones) = _broadcasted_ones(a, b)

# special case due to missing converts for ranges
_range_convert(::Type{AbstractVector{T}}, a::AbstractRange{T}) where T = a
_range_convert(::Type{AbstractVector{T}}, a::AbstractUnitRange) where T = convert(T,first(a)):convert(T,last(a))
_range_convert(::Type{AbstractVector{T}}, a::AbstractRange) where T = convert(T,first(a)):step(a):convert(T,last(a))

function broadcasted(::DefaultArrayStyle{1}, ::typeof(*), a::Ones{T}, b::AbstractRange{V}) where {T,V}
    broadcast_shape(size(a), size(b)) # Check sizes are compatible.
    return _range_convert(AbstractVector{promote_type(T,V)}, b)
end

function broadcasted(::DefaultArrayStyle{1}, ::typeof(*), a::AbstractRange{V}, b::Ones{T}) where {T,V}
    broadcast_shape(size(a), size(b)) # Check sizes are compatible.
    return _range_convert(AbstractVector{promote_type(T,V)}, a)
end

function broadcasted(::DefaultArrayStyle{1}, ::typeof(*), a::AbstractFill, b::AbstractRange)
    broadcast_shape(size(a), size(b)) # Check sizes are compatible.
    return broadcasted(*, getindex_value(a), b)
end

function broadcasted(::DefaultArrayStyle{1}, ::typeof(*), a::AbstractRange, b::AbstractFill)
    broadcast_shape(size(a), size(b)) # Check sizes are compatible.
    return broadcasted(*, a, getindex_value(b))
end

broadcasted(::DefaultArrayStyle{N}, op, r::AbstractFill{T,N}, x::Number) where {T,N} = Fill(op(getindex_value(r),x), size(r))
broadcasted(::DefaultArrayStyle{N}, op, x::Number, r::AbstractFill{T,N}) where {T,N} = Fill(op(x, getindex_value(r)), size(r))
broadcasted(::DefaultArrayStyle{N}, op, r::AbstractFill{T,N}, x::Ref) where {T,N} = Fill(op(getindex_value(r),x[]), size(r))
broadcasted(::DefaultArrayStyle{N}, op, x::Ref, r::AbstractFill{T,N}) where {T,N} = Fill(op(x[], getindex_value(r)), size(r))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""PodStatus represents information about the status of a pod. Status may trail the actual state of a system, especially if the node that hosts the pod cannot contact the control plane.

    IoK8sApiCoreV1PodStatus(;
        conditions=nothing,
        containerStatuses=nothing,
        ephemeralContainerStatuses=nothing,
        hostIP=nothing,
        initContainerStatuses=nothing,
        message=nothing,
        nominatedNodeName=nothing,
        phase=nothing,
        podIP=nothing,
        podIPs=nothing,
        qosClass=nothing,
        reason=nothing,
        startTime=nothing,
    )

    - conditions::Vector{IoK8sApiCoreV1PodCondition} : Current service state of pod. More info: https://kubernetes.io/docs/concepts/workloads/pods/pod-lifecycle#pod-conditions
    - containerStatuses::Vector{IoK8sApiCoreV1ContainerStatus} : The list has one entry per container in the manifest. Each entry is currently the output of `docker inspect`. More info: https://kubernetes.io/docs/concepts/workloads/pods/pod-lifecycle#pod-and-container-status
    - ephemeralContainerStatuses::Vector{IoK8sApiCoreV1ContainerStatus} : Status for any ephemeral containers that have run in this pod. This field is alpha-level and is only populated by servers that enable the EphemeralContainers feature.
    - hostIP::String : IP address of the host to which the pod is assigned. Empty if not yet scheduled.
    - initContainerStatuses::Vector{IoK8sApiCoreV1ContainerStatus} : The list has one entry per init container in the manifest. The most recent successful init container will have ready = true, the most recently started container will have startTime set. More info: https://kubernetes.io/docs/concepts/workloads/pods/pod-lifecycle#pod-and-container-status
    - message::String : A human readable message indicating details about why the pod is in this condition.
    - nominatedNodeName::String : nominatedNodeName is set only when this pod preempts other pods on the node, but it cannot be scheduled right away as preemption victims receive their graceful termination periods. This field does not guarantee that the pod will be scheduled on this node. Scheduler may decide to place the pod elsewhere if other nodes become available sooner. Scheduler may also decide to give the resources on this node to a higher priority pod that is created after preemption. As a result, this field may be different than PodSpec.nodeName when the pod is scheduled.
    - phase::String : The phase of a Pod is a simple, high-level summary of where the Pod is in its lifecycle. The conditions array, the reason and message fields, and the individual container status arrays contain more detail about the pod's status. There are five possible phase values:  Pending: The pod has been accepted by the Kubernetes system, but one or more of the container images has not been created. This includes time before being scheduled as well as time spent downloading images over the network, which could take a while. Running: The pod has been bound to a node, and all of the containers have been created. At least one container is still running, or is in the process of starting or restarting. Succeeded: All containers in the pod have terminated in success, and will not be restarted. Failed: All containers in the pod have terminated, and at least one container has terminated in failure. The container either exited with non-zero status or was terminated by the system. Unknown: For some reason the state of the pod could not be obtained, typically due to an error in communicating with the host of the pod.  More info: https://kubernetes.io/docs/concepts/workloads/pods/pod-lifecycle#pod-phase
    - podIP::String : IP address allocated to the pod. Routable at least within the cluster. Empty if not yet allocated.
    - podIPs::Vector{IoK8sApiCoreV1PodIP} : podIPs holds the IP addresses allocated to the pod. If this field is specified, the 0th entry must match the podIP field. Pods may be allocated at most 1 value for each of IPv4 and IPv6. This list is empty if no IPs have been allocated yet.
    - qosClass::String : The Quality of Service (QOS) classification assigned to the pod based on resource requirements See PodQOSClass type for available QOS classes More info: https://git.k8s.io/community/contributors/design-proposals/node/resource-qos.md
    - reason::String : A brief CamelCase message indicating details about why the pod is in this state. e.g. 'Evicted'
    - startTime::IoK8sApimachineryPkgApisMetaV1Time : RFC 3339 date and time at which the object was acknowledged by the Kubelet. This is before the Kubelet pulled the container image(s) for the pod.
"""
mutable struct IoK8sApiCoreV1PodStatus <: SwaggerModel
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PodCondition} } # spec name: conditions
    containerStatuses::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ContainerStatus} } # spec name: containerStatuses
    ephemeralContainerStatuses::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ContainerStatus} } # spec name: ephemeralContainerStatuses
    hostIP::Any # spec type: Union{ Nothing, String } # spec name: hostIP
    initContainerStatuses::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ContainerStatus} } # spec name: initContainerStatuses
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    nominatedNodeName::Any # spec type: Union{ Nothing, String } # spec name: nominatedNodeName
    phase::Any # spec type: Union{ Nothing, String } # spec name: phase
    podIP::Any # spec type: Union{ Nothing, String } # spec name: podIP
    podIPs::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PodIP} } # spec name: podIPs
    qosClass::Any # spec type: Union{ Nothing, String } # spec name: qosClass
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    startTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: startTime

    function IoK8sApiCoreV1PodStatus(;conditions=nothing, containerStatuses=nothing, ephemeralContainerStatuses=nothing, hostIP=nothing, initContainerStatuses=nothing, message=nothing, nominatedNodeName=nothing, phase=nothing, podIP=nothing, podIPs=nothing, qosClass=nothing, reason=nothing, startTime=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("containerStatuses"), containerStatuses)
        setfield!(o, Symbol("containerStatuses"), containerStatuses)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("ephemeralContainerStatuses"), ephemeralContainerStatuses)
        setfield!(o, Symbol("ephemeralContainerStatuses"), ephemeralContainerStatuses)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("hostIP"), hostIP)
        setfield!(o, Symbol("hostIP"), hostIP)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("initContainerStatuses"), initContainerStatuses)
        setfield!(o, Symbol("initContainerStatuses"), initContainerStatuses)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("nominatedNodeName"), nominatedNodeName)
        setfield!(o, Symbol("nominatedNodeName"), nominatedNodeName)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("phase"), phase)
        setfield!(o, Symbol("phase"), phase)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("podIP"), podIP)
        setfield!(o, Symbol("podIP"), podIP)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("podIPs"), podIPs)
        setfield!(o, Symbol("podIPs"), podIPs)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("qosClass"), qosClass)
        setfield!(o, Symbol("qosClass"), qosClass)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("startTime"), startTime)
        setfield!(o, Symbol("startTime"), startTime)
        o
    end
end # type IoK8sApiCoreV1PodStatus

const _property_map_IoK8sApiCoreV1PodStatus = Dict{Symbol,Symbol}(Symbol("conditions")=>Symbol("conditions"), Symbol("containerStatuses")=>Symbol("containerStatuses"), Symbol("ephemeralContainerStatuses")=>Symbol("ephemeralContainerStatuses"), Symbol("hostIP")=>Symbol("hostIP"), Symbol("initContainerStatuses")=>Symbol("initContainerStatuses"), Symbol("message")=>Symbol("message"), Symbol("nominatedNodeName")=>Symbol("nominatedNodeName"), Symbol("phase")=>Symbol("phase"), Symbol("podIP")=>Symbol("podIP"), Symbol("podIPs")=>Symbol("podIPs"), Symbol("qosClass")=>Symbol("qosClass"), Symbol("reason")=>Symbol("reason"), Symbol("startTime")=>Symbol("startTime"))
const _property_types_IoK8sApiCoreV1PodStatus = Dict{Symbol,String}(Symbol("conditions")=>"Vector{IoK8sApiCoreV1PodCondition}", Symbol("containerStatuses")=>"Vector{IoK8sApiCoreV1ContainerStatus}", Symbol("ephemeralContainerStatuses")=>"Vector{IoK8sApiCoreV1ContainerStatus}", Symbol("hostIP")=>"String", Symbol("initContainerStatuses")=>"Vector{IoK8sApiCoreV1ContainerStatus}", Symbol("message")=>"String", Symbol("nominatedNodeName")=>"String", Symbol("phase")=>"String", Symbol("podIP")=>"String", Symbol("podIPs")=>"Vector{IoK8sApiCoreV1PodIP}", Symbol("qosClass")=>"String", Symbol("reason")=>"String", Symbol("startTime")=>"IoK8sApimachineryPkgApisMetaV1Time")
Base.propertynames(::Type{ IoK8sApiCoreV1PodStatus }) = collect(keys(_property_map_IoK8sApiCoreV1PodStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodStatus[property_name]

function check_required(o::IoK8sApiCoreV1PodStatus)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PodStatus }, name::Symbol, val)
end


text:
export linear_interpolation!, face_average!, upwind_values!, Neumann_conditions_impose!, volume_gradient!, face_gradient!, flux_integral!

"""
Linearly interpolate values from volumes to faces

* `msh`: the mesh at hand
* `at_volumes`: array of values at volume centers (input)
* `at_faces`: array of values at face centers (output)
"""
function linear_interpolation!(msh::FinVolMesh{N}, at_volumes::Vector{Float64}, at_faces::Vector{Float64}) where N
    v1=0.0
    v2=0.0

    for (i, (fc_adj, fc_dists)) in enumerate(zip(msh.face_connectivity, msh.face_dists))
        if fc_adj[1]==fc_adj[2] && fc_adj[1]==0
            throw(
                error(
                    "linear_interpolation!:ERROR:unconnected face detected"
                )
            )
        end

        v1=(fc_adj[1]==0 ? 0.0 : at_volumes[fc_adj[1]])
        v2=(fc_adj[2]==0 ? 0.0 : at_volumes[fc_adj[2]])

        at_faces[i]=(v1*fc_dists[2]+v2*fc_dists[1])/(fc_dists[1]+fc_dists[2])
    end
end

"""
Average volume values at face centers.
When a face is connected to only one volume, it incorporates the value of that volume

* `msh`: the mesh at hand
* `at_volumes`: array of values at volume centers (input)
* `at_faces`: array of values at face centers (output)
"""
function face_average!(msh::FinVolMesh{N}, at_volumes::Vector{Float64}, at_faces::Vector{Float64}) where N
    for (fc, cons) in enumerate(msh.face_connectivity)
        if (cons[1]==0 || cons[2]==0)
            if (cons[1]==cons[2])
                throw(
                    error(
                        "face_average!:ERROR:face $(fc) not connected to any volumes"
                    )
                )
            end

            if cons[1]==0
                at_faces[fc]=at_volumes[cons[2]]
            else
                at_faces[fc]=at_volumes[cons[1]]
            end
        else
            at_faces[fc]=(at_volumes[cons[1]]+at_volumes[cons[2]])/2
        end
    end
end

"""
Fetch values at faces from upstream volumes

* `msh`: the mesh at hand
* `at_volumes`: values at volume centers (input)
* `at_faces`: values at face centers (output)
* `u`, `v`, ... and any other set of velocities **at faces** with which to determine
    flow direction
"""
function upwind_values!(msh::FinVolMesh{N}, at_volumes::Vector{Float64}, at_faces::Vector{Float64}, vels...) where N
    nv=0.0

    for (fc, fccons) in enumerate(msh.face_connectivity)
        if fccons[1]==0
            if fccons[2]==0
                throw(
                    error(
                        "upwind_values!:ERROR:unconnected face detected"
                    )
                )
            end

            at_faces[fc]=at_volumes[fccons[2]]
        elseif fccons[2]==0
            at_faces[fc]=at_volumes[fccons[1]]
        else
            nv=0.0

            for (n, v) in zip(msh.face_normals[fc], vels)
                nv+=n*v[fc]
            end

            if nv>=0.0
                at_faces[fc]=at_volumes[fccons[1]]
            else
                at_faces[fc]=at_volumes[fccons[2]]
            end
        end
    end
end

"""
Impose Neumann conditions for certain faces

* `msh`: the mesh at hand
* `at_volumes`: values at volume centers, from which to apply variations
* `at_faces`: vector with values at all face centers
* `Neumann_conditions`: tuple with an array of face indexes
    and an array of float arrays with the imposed gradients at those faces

    Example:

    ```
    Neumann_conditions=(
        [1, 3, 4], # face indexes
        [ # values of gradients at faces
            [1.0, 0.0], # ∇f=(1, 0)
            [0.0, 1.0], # ∇f=(0, 1)
            [0.5, 0.5] # ∇f=(1/2, 1/2)
        ]
    )
    ```
"""
function Neumann_conditions_impose!(msh::FinVolMesh{N},
    at_volumes::Vector{Float64}, at_faces::Vector{Float64},
    Neumann_conditions::Tuple{Vector{Int64}, Vector{Vector{Float64}}}) where N
    ccell=0

    for (fc, nc) in zip(Neumann_conditions...)
        ccell=msh.face_connectivity[fc][1]

        if ccell==0
            ccell=msh.face_connectivity[fc][2]
        end

        if ccell==0
            throw(
                error(
                    "Neumann_conditions_impose!:ERROR:face disconnected from any cells"
                )
            )
        end

        at_faces[fc]=at_volumes[ccell]

        for (nd, d) in enumerate(nc)
            at_faces[fc]+=d*(msh.face_centers[fc][nd]-msh.volume_centers[ccell][nd])
        end
    end
end

"""
Green-Gauss gradient at cell centroids

* `msh`: the mesh at hand
* `at_volumes`: gradients at cell centroids (output)
* `at_faces`: values at faces (input)
* `dim`: the dimension along which to obtain the gradient
"""
function volume_gradient!(msh::FinVolMesh{N}, at_volumes::Vector{Float64}, at_faces::Vector{Float64}, dim::Int64) where N
    at_volumes.=0.0

    for (i, (cons, shat)) in enumerate(zip(msh.face_connectivity, msh.face_normals))
        if cons[1]!=0
            at_volumes[cons[1]]+=shat[dim]*at_faces[i]
        end

        if cons[2]!=0
            at_volumes[cons[2]]-=shat[dim]*at_faces[i]
        end
    end

    for (i, v) in enumerate(msh.volume_volumes)
        at_volumes[i]/=v
    end
end

"""
Gradients at faces (averaged from finite differences with adjacent volumes)

* `msh`: the mesh at hand
* `at_volumes`: array of values at volume centers (input)
* `at_faces`: array of values at face centers (input and output)
* `dim`: spatial dimension for gradient
"""
function face_gradient!(msh::FinVolMesh{N}, at_volumes::Vector{Float64}, at_faces::Vector{Float64}, dim::Int64) where N
    for (fc, (cons, shat, dists)) in enumerate(zip(msh.face_connectivity, msh.face_normal_unit, msh.face_dists))
        if (cons[1]==0 || cons[2]==0)
            if (cons[1]==cons[2])
                throw(
                    error(
                        "face_average!:ERROR:face $(fc) not connected to any volumes"
                    )
                )
            end

            if cons[1]==0
                at_faces[fc]=(at_volumes[cons[2]]-at_faces[fc])*shat[dim]/dists[2]
            else
                at_faces[fc]=(at_faces[fc]-at_volumes[cons[1]])*shat[dim]/dists[1]
            end
        else
            at_faces[fc]=((at_volumes[cons[2]]-at_faces[fc])*shat[dim]/dists[2]+
                (at_faces[fc]-at_volumes[cons[1]])*shat[dim]/dists[1])/2
        end
    end
end

"""
Accumulate fluxes from faces into variations at cell centroids

Example:

Given `f1, f2, f3` in a three-dimensional mesh,
`du/dt=-1/v∫∇⋅fdΩ` is calculated at cell centroids and stored at the given vectors

* `msh`: the mesh at hand
* `at_volumes`: flux integral according to Gauss's theorem (output)
* `u, v, w...` and remaining arguments: fluxes at each spatial dimension, as vectors for face values (input)
"""
function flux_integral!(msh::FinVolMesh{N}, at_volumes::Vector{Float64}, fs...) where N
    at_volumes.=0.0

    for (i, (cons, n)) in enumerate(zip(msh.face_connectivity, msh.face_normals))
        if cons[1]!=0
            for (j, f) in enumerate(fs)
                at_volumes[cons[1]]-=n[j]*f[i]
            end
        end

        if cons[2]!=0
            for (j, f) in enumerate(fs)
                at_volumes[cons[2]]+=n[j]*f[i]
            end
        end
    end

    for (i, v) in enumerate(msh.volume_volumes)
        at_volumes[i]/=v
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    GitRemote(repo::GitRepo, rmt_name::AbstractString, rmt_url::AbstractString) -> GitRemote

Look up a remote git repository using its name and URL. Uses the default fetch refspec.

# Examples
```julia
repo = LibGit2.init(repo_path)
remote = LibGit2.GitRemote(repo, "upstream", repo_url)
```
"""
function GitRemote(repo::GitRepo, rmt_name::AbstractString, rmt_url::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_create, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Cstring),
                rmt_ptr_ptr, repo.ptr, rmt_name, rmt_url)
    return GitRemote(repo, rmt_ptr_ptr[])
end

"""
    GitRemote(repo::GitRepo, rmt_name::AbstractString, rmt_url::AbstractString, fetch_spec::AbstractString) -> GitRemote

Look up a remote git repository using the repository's name and URL,
as well as specifications for how to fetch from the remote
(e.g. which remote branch to fetch from).

# Examples
```julia
repo = LibGit2.init(repo_path)
refspec = "+refs/heads/mybranch:refs/remotes/origin/mybranch"
remote = LibGit2.GitRemote(repo, "upstream", repo_url, refspec)
```
"""
function GitRemote(repo::GitRepo, rmt_name::AbstractString, rmt_url::AbstractString, fetch_spec::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_create_with_fetchspec, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Cstring, Cstring),
                rmt_ptr_ptr, repo.ptr, rmt_name, rmt_url, fetch_spec)
    return GitRemote(repo, rmt_ptr_ptr[])
end

"""
    GitRemoteAnon(repo::GitRepo, url::AbstractString) -> GitRemote

Look up a remote git repository using only its URL, not its name.

# Examples
```julia
repo = LibGit2.init(repo_path)
remote = LibGit2.GitRemoteAnon(repo, repo_url)
```
"""
function GitRemoteAnon(repo::GitRepo, url::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_create_anonymous, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring),
                rmt_ptr_ptr, repo.ptr, url)
    return GitRemote(repo, rmt_ptr_ptr[])
end

"""
    lookup_remote(repo::GitRepo, remote_name::AbstractString) -> Nullable{GitRemote}

Determine if the `remote_name` specified exists within the `repo`. Returns a
[`Nullable`](@ref), which will be null if the requested remote does not exist. If the remote
does exist, the `Nullable` contains a [`GitRemote`](@ref) to the remote name.

# Examples
```julia
repo = LibGit2.GitRepo(path)
remote_name = "test"
isnull(LibGit2.lookup_remote(repo, remote_name)) # will return true
```
"""
function lookup_remote(repo::GitRepo, remote_name::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    err = ccall((:git_remote_lookup, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring),
                rmt_ptr_ptr, repo.ptr, remote_name)
    if err == Int(Error.GIT_OK)
        return Nullable{GitRemote}(GitRemote(repo, rmt_ptr_ptr[]))
    elseif err == Int(Error.ENOTFOUND)
        return Nullable{GitRemote}()
    else
        throw(Error.GitError(err))
    end
end

function get(::Type{GitRemote}, repo::GitRepo, rmt_name::AbstractString)
    rmt_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_remote_lookup, :libgit2), Cint,
                (Ptr{Ptr{Void}}, Ptr{Void}, Cstring),
                rmt_ptr_ptr, repo.ptr, rmt_name)
    return GitRemote(repo, rmt_ptr_ptr[])
end

"""
    url(rmt::GitRemote)

Get the fetch URL of a remote git repository.

# Examples
```julia-repl
julia> repo_url = "https://github.com/JuliaLang/Example.jl";

julia> repo = LibGit2.init(mktempdir());

julia> remote = LibGit2.GitRemote(repo, "origin", repo_url);

julia> LibGit2.url(remote)
"https://github.com/JuliaLang/Example.jl"
```
"""
function url(rmt::GitRemote)
    url_ptr = ccall((:git_remote_url, :libgit2), Cstring, (Ptr{Void},), rmt.ptr)
    url_ptr == C_NULL && return ""
    return unsafe_string(url_ptr)
end

"""
    push_url(rmt::GitRemote)

Get the push URL of a remote git repository.

# Examples
```julia-repl
julia> repo_url = "https://github.com/JuliaLang/Example.jl";

julia> repo = LibGit2.init(mktempdir());

julia> LibGit2.set_remote_push_url(repo, "origin", repo_url);

julia> LibGit2.push_url(LibGit2.get(LibGit2.GitRemote, repo, "origin"))
"https://github.com/JuliaLang/Example.jl"
```
"""
function push_url(rmt::GitRemote)
    url_ptr = ccall((:git_remote_pushurl, :libgit2), Cstring, (Ptr{Void},), rmt.ptr)
    url_ptr == C_NULL && return ""
    return unsafe_string(url_ptr)
end

"""
    name(rmt::GitRemote)

Get the name of a remote repository, for instance `"origin"`.
If the remote is anonymous (see [`GitRemoteAnon`](@ref))
the name will be an empty string `""`.

# Examples
```julia-repl
julia> repo_url = "https://github.com/JuliaLang/Example.jl";

julia> repo = LibGit2.clone(cache_repo, "test_directory");

julia> remote = LibGit2.GitRemote(repo, "origin", repo_url);

julia> name(remote)
"origin"
```
"""
function name(rmt::GitRemote)
    name_ptr = ccall((:git_remote_name, :libgit2), Cstring, (Ptr{Void},), rmt.ptr)
    name_ptr == C_NULL && return ""
    return unsafe_string(name_ptr)
end

"""
    fetch_refspecs(rmt::GitRemote) -> Vector{String}

Get the *fetch* refspecs for the specified `rmt`. These refspecs contain
information about which branch(es) to fetch from.

# Examples
```julia-repl
julia> remote = LibGit2.get(LibGit2.GitRemote, repo, "upstream");

julia> LibGit2.add_fetch!(repo, remote, "upstream");

julia> LibGit2.fetch_refspecs(remote)
String["+refs/heads/*:refs/remotes/upstream/*"]
```
"""
function fetch_refspecs(rmt::GitRemote)
    sa_ref = Ref(StrArrayStruct())
    @check ccall((:git_remote_get_fetch_refspecs, :libgit2), Cint,
                 (Ptr{StrArrayStruct}, Ptr{Void}), sa_ref, rmt.ptr)
    res = convert(Vector{String}, sa_ref[])
    free(sa_ref)
    res
end

"""
    push_refspecs(rmt::GitRemote) -> Vector{String}

Get the *push* refspecs for the specified `rmt`. These refspecs contain
information about which branch(es) to push to.

# Examples
```julia-repl
julia> remote = LibGit2.get(LibGit2.GitRemote, repo, "upstream");

julia> LibGit2.add_push!(repo, remote, "refs/heads/master");

julia> close(remote);

julia> remote = LibGit2.get(LibGit2.GitRemote, repo, "upstream");

julia> LibGit2.push_refspecs(remote)
String["refs/heads/master"]
```
"""
function push_refspecs(rmt::GitRemote)
    sa_ref = Ref(StrArrayStruct())
    @check ccall((:git_remote_get_push_refspecs, :libgit2), Cint,
                 (Ptr{StrArrayStruct}, Ptr{Void}), sa_ref, rmt.ptr)
    res = convert(Vector{String}, sa_ref[])
    free(sa_ref)
    res
end

"""
    add_fetch!(repo::GitRepo, rmt::GitRemote, fetch_spec::String)

Add a *fetch* refspec for the specified `rmt`. This refspec will contain
information about which branch(es) to fetch from.

# Examples
```julia-repl
julia> LibGit2.add_fetch!(repo, remote, "upstream");

julia> LibGit2.fetch_refspecs(remote)
String["+refs/heads/*:refs/remotes/upstream/*"]
```
"""
function add_fetch!(repo::GitRepo, rmt::GitRemote, fetch_spec::String)
    @check ccall((:git_remote_add_fetch, :libgit2), Cint,
                 (Ptr{Void}, Cstring, Cstring), repo.ptr,
                 name(rmt), fetch_spec)
end

"""
    add_push!(repo::GitRepo, rmt::GitRemote, push_spec::String)

Add a *push* refspec for the specified `rmt`. This refspec will contain
information about which branch(es) to push to.

# Examples
```julia-repl
julia> LibGit2.add_push!(repo, remote, "refs/heads/master");

julia> remote = LibGit2.get(LibGit2.GitRemote, repo, branch);

julia> LibGit2.push_refspecs(remote)
String["refs/heads/master"]
```

!!! note
    You may need to [`close`](@ref) and reopen the `GitRemote`
    in question after updating its push refspecs in order for
    the change to take effect and for calls to [`push`](@ref)
    to work.
"""
function add_push!(repo::GitRepo, rmt::GitRemote, push_spec::String)
    @check ccall((:git_remote_add_push, :libgit2), Cint,
                 (Ptr{Void}, Cstring, Cstring), repo.ptr,
                 name(rmt), push_spec)
end

"""
    fetch(rmt::GitRemote, refspecs; options::FetchOptions=FetchOptions(), msg="")

Fetch from the specified `rmt` remote git repository, using `refspecs` to
determine which remote branch(es) to fetch.
The keyword arguments are:
  * `options`: determines the options for the fetch, e.g. whether to prune afterwards.
    See [`FetchOptions`](@ref) for more information.
  * `msg`: a message to insert into the reflogs.
"""
function fetch(rmt::GitRemote, refspecs::Vector{<:AbstractString};
               options::FetchOptions = FetchOptions(),
               msg::AbstractString="")
    msg = "libgit2.fetch: $msg"
    @check ccall((:git_remote_fetch, :libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{FetchOptions}, Cstring),
                 rmt.ptr, isempty(refspecs) ? C_NULL : refspecs, Ref(options), msg)
end

"""
    push(rmt::GitRemote, refspecs; force::Bool=false, options::PushOptions=PushOptions())

Push to the specified `rmt` remote git repository, using `refspecs` to
determine which remote branch(es) to push to.
The keyword arguments are:
  * `force`: if `true`, a force-push will occur, disregarding conflicts.
  * `options`: determines the options for the push, e.g. which proxy headers to use.
    See [`PushOptions`](@ref) for more information.

!!! note
    You can add information about the push refspecs in two other ways: by setting
    an option in the repository's `GitConfig` (with `push.default` as the key) or
    by calling [`add_push!`](@ref). Otherwise you will need to explicitly specify
    a push refspec in the call to `push` for it to have any effect, like so:
    `LibGit2.push(repo, refspecs=["refs/heads/master"])`.
"""
function push(rmt::GitRemote, refspecs::Vector{<:AbstractString};
              force::Bool = false, options::PushOptions = PushOptions())
    @check ccall((:git_remote_push, :libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{PushOptions}),
                 rmt.ptr, isempty(refspecs) ? C_NULL : refspecs, Ref(options))
end

"""
    remote_delete(repo::GitRepo, remote_name::AbstractString) -> Void

Delete the `remote_name` from the git `repo`.
"""
function remote_delete(repo::GitRepo, remote_name::AbstractString)
    @check ccall((:git_remote_delete, :libgit2), Cint,
                 (Ptr{Void}, Cstring),
                 repo.ptr, remote_name)
end

Base.show(io::IO, rmt::GitRemote) = print(io, "GitRemote:\nRemote name: ", name(rmt), " url: ", url(rmt))

"""
    set_remote_fetch_url(repo::GitRepo, remote_name, url)
    set_remote_fetch_url(path::String, remote_name, url)

Set the fetch `url` for the specified `remote_name` for the [`GitRepo`](@ref) or the git repository
located at `path`. Typically git repos use `"origin"` as the remote name.
"""
function set_remote_fetch_url end

function set_remote_fetch_url(repo::GitRepo, remote_name::AbstractString, url::AbstractString)
    @check ccall((:git_remote_set_url, :libgit2), Cint,
                 (Ptr{Void}, Cstring, Cstring),
                 repo.ptr, remote_name, url)
end

function set_remote_fetch_url(path::AbstractString, remote_name::AbstractString, url::AbstractString)
    with(GitRepo, path) do repo
        set_remote_fetch_url(repo, remote_name, url)
    end
end

"""
    set_remote_push_url(repo::GitRepo, remote_name, url)
    set_remote_push_url(path::String, remote_name, url)

Set the push `url` for the specified `remote_name` for the [`GitRepo`](@ref) or the git repository
located at `path`. Typically git repos use `"origin"` as the remote name.
"""
function set_remote_push_url end

function set_remote_push_url(repo::GitRepo, remote_name::AbstractString, url::AbstractString)
    @check ccall((:git_remote_set_pushurl, :libgit2), Cint,
                 (Ptr{Void}, Cstring, Cstring),
                 repo.ptr, remote_name, url)
end

function set_remote_push_url(path::AbstractString, remote_name::AbstractString, url::AbstractString)
    with(GitRepo, path) do repo
        set_remote_push_url(repo, remote_name, url)
    end
end

"""
    set_remote_url(repo::GitRepo, remote_name, url)
    set_remote_url(repo::String, remote_name, url)

Set both the fetch and push `url` for `remote_name` for the [`GitRepo`](@ref) or the git repository
located at `path`. Typically git repos use `"origin"` as the remote name.

# Examples
```julia
repo_path = joinpath(tempdir(), "Example")
repo = LibGit2.init(repo_path)
LibGit2.set_remote_url(repo, "upstream", "https://github.com/JuliaLang/Example.jl")
LibGit2.set_remote_url(repo_path, "upstream2", "https://github.com/JuliaLang/Example2.jl")
```
"""
function set_remote_url end

function set_remote_url(repo::GitRepo, remote_name::AbstractString, url::AbstractString)
    set_remote_fetch_url(repo, remote_name, url)
    set_remote_push_url(repo, remote_name, url)
end

function set_remote_url(path::AbstractString, remote_name::AbstractString, url::AbstractString)
    with(GitRepo, path) do repo
        set_remote_url(repo, remote_name, url)
    end
end

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
using SAG
using Test

@testset "SAG.jl" begin
    # Test constructors
    dims = (3, 2)
    n = 3
    ∇s = [randn(dims) for _ in 1:n]
    sg = StochasticGradient(∇s, zeros(dims), overwrite=false)
    @test sg ≈ zeros(dims)
    @test sg ≈ sg.∇
    sg = StochasticGradient(∇s, zeros(dims))
    @test sg ≈ sum(∇s)
    sg = StochasticGradient(∇s)
    @test sg ≈ sum(∇s)
    ∇ = randn(dims)
    sg = StochasticGradient(copy(∇), n)
    @test sg ≈ ∇
    @test sum(sg.∇s) ≈ ∇
    sg = StochasticGradient(n, dims)
    @test sg ≈ zeros(dims)
    @test sum(sg.∇s) ≈ zeros(dims)
    sg = StochasticGradient(Int, n, dims)
    @test sg ≈ zeros(dims)
    @test sum(sg.∇s) ≈ zeros(dims)
    @test eltype(sg) == Int
    @test eltype(sg.∇) == Int
    @test eltype(sg.∇s[1]) == Int

    # Test updates
    dims = (3, 2)
    n = 3
    i = 1
    ∇i = randn(dims)
    sg = StochasticGradient(n, dims)
    update!(sg, i, ∇i)
    @test sg.∇s[i] ≈ ∇i
    @test sg ≈ ∇i
    @test initialized_fraction(sg) ≈ 1/3
    sg = StochasticGradient(n, dims)
    update!(sg, i=>∇i)
    @test sg.∇s[i] ≈ ∇i
    @test sg ≈ ∇i
    @test initialized_fraction(sg) ≈ 1/3
    sg = StochasticGradient(n, dims)
    update!(sg, [1=>∇i, 2=>∇i])
    @test sg.∇s[1] ≈ ∇i
    @test sg.∇s[2] ≈ ∇i
    @test sg ≈ 2∇i
    @test initialized_fraction(sg) ≈ 2/3
    sg = StochasticGradient(n, dims)
    update!(sg, zip([1,2], [∇i, ∇i]))
    @test sg.∇s[1] ≈ ∇i
    @test sg.∇s[2] ≈ ∇i
    @test sg ≈ 2∇i
    @test initialized_fraction(sg) ≈ 2/3
end


text:
using DynamicalSystems
using CausalityTools
using Plots
plotlyjs()

s = DynamicalSystems.Systems.coupledstandardmaps(2)
ts_length = 500
orbit = trajectory(s, ts_length - 1, Ttr = 1000)
x1, x2, y1, y2 = ([orbit[:, i] for i = 1:4]...,)

px = plot(x1, label = "x1", lc = :blue)
plot!(x2, label = "x2", lc = :blue, ls = :dash)
ylabel!("Value")

py = plot(y1, label = "y1", lc = :green)
plot!(y2, label = "y2", lc = :green, ls = :dash)
xlabel!("Time step")
ylabel!("Value")

plot(px, py, layout = (2, 1))

E = StateSpaceReconstruction.embed([x1, x2, y2], [1, 2, 3], [0, -10, -5])
E1, E2, E3 = ([E.points[i, :] for i = 1:3]...,)
scatter(E1, E2, E3, legend = false)

using RecipesBase
@recipe function f(E::StateSpaceReconstruction.AbstractEmbedding; vars = [1, 2, 3],
        mc = :black, ms = 2, lc = :black, lw = 1)
    if size(E.points, 1) < 3
        error("Cannot plot reconstruction with < 3 variables.")
    end
    legend --> false
    xlabel --> "x1"
    ylabel --> "x2"
    zlabel --> "x3"
    markercolor --> mc
    markersize --> ms
    ([E.points[i, :] for i in vars]...,)
end

scatter(E, ms = 0.5)

#########################################################
# Helper functions to plot 3D rectangles.
#########################################################
"""
    rectangle3dpts(x, y, z, ϵx, ϵy, ϵz)

Compute the coordinates of the vertices of a 3D rectangle, given the
coordinates of the origin (`x`, `y`, and `z`) and the corresponding edge lengths
(`ϵx`, `ϵy`, `ϵz`).
"""
function rectangle3dpts(x, y, z, ϵx, ϵy, ϵz)
    v1b = [x,    y,    z]
    v2b = [x+ϵx, y,    z]
    v3b = [x+ϵx, y+ϵy, z]
    v4b = [x,    y+ϵy, z]
    v1t = [x,    y,    z+ϵz]
    v2t = [x+ϵx, y,    z+ϵz]
    v3t = [x+ϵx, y+ϵy, z+ϵz]
    v4t = [x,    y+ϵy, z+ϵz]
    (v1b, v2b, v3b, v4b, v1t, v2t, v3t, v4t)
end

"""
    rectangle3dpts(o, ϵ)

Given the origin `o` (3-element vector) and edge lengths `ϵ` (also a
3-element vector) of a 3D rectangle, return the coordinates of its
vertices.
"""
function rectangle3dpts(o, ϵ)
    x, y, z = (o...,)
    ϵx, ϵy, ϵz = (ϵ...,)
    rectangle3dpts(x, y, z, ϵx, ϵy, ϵz)
end

"""
    connectvertices(o, ϵ)

Given the origin `o` (3-element vector) and edge lengths `ϵ` (also a
3-element vector) of a 3D rectangle, return a vector of line segments
(each a dim-by-n_vertices array) that when plotted together yields
the entire rectangle.
"""
function connectvertices(o, ϵ)
    # Get the vertices
    v1b, v2b, v3b, v4b, v1t, v2t, v3t, v4t = rectangle3dpts(o, ϵ)

    # Connect vertices in top and bottom planes
    connect_top = zeros(3, 5)
    connect_bottom = zeros(3, 5)

    # Connect corners of the top and bottom planes
    corner1 = zeros(3, 2)
    corner2 = zeros(3, 2)
    corner3 = zeros(3, 2)
    corner4 = zeros(3, 2)
    connect_bottom[:, 1:5] = [v1b v2b v3b v4b v1b]
    connect_top[:, 1:5]    = [v1t v2t v3t v4t v1t]
    corner1[:, 1:2] = [v1b v1t]
    corner2[:, 1:2] = [v2b v2t]
    corner3[:, 1:2] = [v3b v3t]
    corner4[:, 1:2] = [v4b v4t]

    return [connect_top, connect_bottom,
            corner1, corner2, corner3, corner4]
end

"""
    splitaxes(x)

Return a vector of the individual components of an array of points
provided as an array where each point is a column.
"""
splitaxes(x) = ([x[k, :] for k = 1:size(x, 1)]...,)

"""
    plot_3D_rect!(p, origin, edgelengths;
             lc = :black, lw = 1.0, ls = :solid)

Append a a 3D rectangle to an existing plot `p`, given the
origin `o` (3-element vector) and edge lengths
`ϵ` (3-element vector) of the rectangle.
"""
function plot_3D_rect!(p, o, ϵ;
             lc = :black, lw = 1.0, ls = :solid, lα = 0.7)

    linesegments = connectvertices(o, ϵ)
    for segment in linesegments
        plot!(p, splitaxes(segment),
            lc = lc, lw = lw, ls = ls, lα = lα)
    end
end

"""
Given the `origin` of a hyperrectangle and the `stepsizes` along each
coordinate axis, divide each axis into `N` stepsizes starting
from `origin[i] + (stepsizes[i]/N)/2`
"""
function fill_hyperrectangle_3D(origin, stepsizes, N = 10)
    xs, ys, zs = ([LinRange(origin[i] + stepsizes[i]/N/2, origin[i] + stepsizes[i] - stepsizes[i]/N/2, N) for i = 1:length(origin)]...,)
    pts = zeros(Float64, 3, length(xs)^3)
    i = 0
    for x in xs
        for y in ys
            for z in zs
                i += 1
                pts[:, i] = [x, y, z]
            end
        end
    end
    pts
end

@recipe function f(N, riv::RectangularInvariantMeasure;
        lw = 0.8, lc = :black, lα = 0.5, ls = :do,
        ms = 0.2, mc = :black, mα = 0.3)

    pts = riv.E.points
    ϵ = riv.ϵ
    axisminima, stepsizes = minima_and_stepsizes(pts, ϵ)

    # All bins are assigned a bin origin.
    v = unique(riv.visited_bins_coordinates, dims = 2)

    # Bins may be visited by several times, so we'll get rid
    # of the repetitions.
    n_visited_boxes = size(v, 2)
    D = size(v, 1)

    # Plot the partition grid over the points of the reconstructed
    # orbit.
    legend --> false

    for i = 1:n_visited_boxes
        origin = riv.visited_bins_coordinates[:, i]
        # Get the measure of the box, and fill the box with a
        # number of points scaled to the measure.
        μ = riv.measure.dist[i] * 50
        fillpts = fill_hyperrectangle_3D(origin, stepsizes, N)
        @series begin
            seriestype := :scatter3d
            markersize --> μ/2
            markeralpha --> μ
            markercolor --> mc
            linecolor --> mc
            fillcolor --> :transparent
            splitaxes(fillpts)
        end

        linesegments = connectvertices(origin, stepsizes)

        for segment in linesegments
            @series begin
                seriestype := :path
                linealpha --> lα
                linewidth --> lw
                linecolor --> lc
                splitaxes(segment)
            end
        end
    end
end

@recipe function f(r::RectangularBinningTransferOperator)
    seriestype  :=  :heatmap
    xlabel --> "state #i"
    ylabel --> "state #j"
    r.transfermatrix
end

@recipe function f(id::InvariantDistribution)
    seriestype  :=  :bar
    xlabel --> "state # (i)"
    ylabel --> "invariant density"
    linecolor --> :black
    fillcolor --> :black
    fillalpha --> 0.5
    legend --> :none
    id.dist
end

@recipe function plot_grid(pts::AbstractArray{T, 2}, riv::RectangularInvariantMeasure; N = 5,
        lw = 0.8, lc = :black, lα = 0.5, ls = :do,
        ms = 0.2, mc = :black, mα = 0.3) where {T}

    axisminima, stepsizes = minima_and_stepsizes(pts, riv.ϵ)

    # All bins are assigned a bin origin.
    v = unique(riv.visited_bins_coordinates, dims = 2)

    # Bins may be visited by several times, so we'll get rid
    # of the repetitions.
    n_visited_boxes = size(v, 2)
    D = size(v, 1)

    # Plot the partition grid over the points of the reconstructed
    # orbit.
    legend --> false

    for i = 1:n_visited_boxes
        origin = riv.visited_bins_coordinates[:, i]
        # Get the measure of the box, and fill the box with a
        # number of points scaled to the measure.
        μ = riv.measure.dist[i]

        println(ceil(Int, μ*100))
        fillpts = fill_hyperrectangle_3D(origin, stepsizes, ceil(Int, μ*100))
        @series begin
            seriestype := :scatter3d
            markersize --> μ
            markeralpha --> μ
            markercolor --> mc
            linecolor --> mc
            fillcolor --> mc
            splitaxes(fillpts)
        end

        linesegments = connectvertices(origin, stepsizes)

        for segment in linesegments
            @series begin
                seriestype := :path
                linealpha --> lα
                linewidth --> lw
                linecolor --> lc
                splitaxes(segment)
            end
        end
    end
end

invm = rectangularinvariantmeasure(E, 4)

gr()
p1 = plot(invm.transfermatrix);
p2 = plot(invm.measure);
p3 = scatter(E, ms = 1.0);
p4 = plot(5, invm, lα = 0.2, ls = :dot,
                    ms = 0.5, mα = 0.5, mc = :blue);

plot(plot(p1, p2, link = :x, layout = (2, 1)),
    plot(p3, p4, layout = (2, 1)), layout = (1, 2))
savefig("docs/src/scripts/test.pdf")


text:
@testset "Lake State" begin
    model = initialise(experiment = Experiment(nutrient_series = Constant()))
    @test LimnoSES.nutrient_load!(model, model.nutrient_series) == nothing

    model = initialise(experiment = Experiment(nutrient_series = Dynamic()))
    model[1].affectors = 80
    @test model.init_nutrients < LimnoSES.nutrient_load!(model, model.nutrient_series)

    model = initialise(
        experiment = Experiment(
            target_nutrients = 3.0,
            nutrient_series = TransientUp(start_year = 5, post_target_series = Constant()),
        ),
    )
    LimnoSES.nutrient_load!(model, model.nutrient_series)
    @test model.lake.p.nutrients == model.init_nutrients
    model.year = 6
    @test LimnoSES.nutrient_load!(model, model.nutrient_series) ≈ 0.8

    model = initialise(
        experiment = Experiment(
            target_nutrients = 0.1,
            nutrient_series = TransientDown(
                start_year = 5,
                post_target_series = Constant(),
            ),
        ),
    )
    LimnoSES.nutrient_load!(model, model.nutrient_series)
    @test model.lake.p.nutrients == model.init_nutrients
    model.year = 6
    @test LimnoSES.nutrient_load!(model, model.nutrient_series) ≈ 0.6
end


text:
using Flux, FluxJS, DataFlow, MacroTools
using FluxJS: traceλ
using JSExpr, JSON, WebIO
using Blink
using Test
using Flux.Tracker: data

atomshell = Blink.AtomShell.isinstalled()

if !atomshell
    Blink.AtomShell.install()
end

include("./helpers.jl")

@testset "FluxJS" begin

m = Dense(10,5)
v = traceλ(m, rand(10))
ex = prettify(DataFlow.syntax(v))

@test @capture ex _ -> (+).(matVecMul(_,_),_)

w = setupWindow()

x = rand(Float32, 2)
@testset "identity" begin
    testjs(w, identity, x)
end

@testset "*" begin
    m = (x) -> ones(2, 2) * x
    testjs(w, m, x)
    testjs(w, m, rand(2, 1))
end

@testset "+" begin
    m = (x) -> x + ones(2)
    testjs(w, m, x)
end

@testset "vcat" begin
    m = (x) -> vcat(x, ones(2))
    testjs(w, m, x)

    m = (x) -> vcat(x, ones(2, 2))
    testjs(w, m, rand(3, 2))
end

@testset "softmax" begin
    m = x -> softmax(x)
    testjs(w, m, x)
end

@testset "Dense Layer" begin
    m = Chain(Dense(2, 2))
    testjs(w, m, x)
end

@testset "Chain" begin
    m = Chain(identity, identity)
    testjs(w, m, x)
end

@testset "broadcast" begin
    b = ones(2)

    m = x -> x .+ b
    testjs(w, m, x)

    m = x -> x .* b
    testjs(w, m, x)

    m = x -> x ./ b
    testjs(w, m, x)

    m = x -> x .- b
    testjs(w, m, x)

    m = x -> exp.(x)
    testjs(w, m, x)

    m = x -> log.(x)
    testjs(w, m, x)

    m = x -> x .^ [2]
    testjs(w, m, x)

    m = x -> σ.(x)
    testjs(w, m, x)

    m = x -> tanh.(x)
    testjs(w, m, x)

    m = x -> relu.(x)
    testjs(w, m, x)

    m = x -> leakyrelu.(x)
    testjs(w, m, x)

    m = x -> copy.(x)
    testjs(w, m, x)
end

@testset "Conv" begin
    m = Chain(Conv((2, 2), 2=>2))
    x = rand(Float32, 2,2,2,2)
    testjs(w, m, x)
end

@testset "maxpool" begin
    m = Chain(Conv((2, 2), 2=>2), MaxPool((2, 2)))
    x = rand(Float32, 4,4,2,2)
    testjs(w, m, x)
end

@testset "BatchNorm" begin
    m = Chain(BatchNorm(10))
    Flux.testmode!(m)
    x = rand(10, 1)
    testjs(w, m, x)
end

@testset "reshape" begin
    m = x -> reshape(x, :, size(x, 3))
    x = rand(1, 2, 3)
    testjs(w, m, x)
end

@testset "RNN" begin
    m = Chain(RNN(10, 10))
    x = rand(10)
    testjs(w, m, x)
end

@testset "LSTM" begin
    m = x -> Flux.gate(x, 4, 1)
    x = rand(10)
    testjs(w, m, x)

    m = Chain(LSTM(10, 10))
    x = rand(10)
    testjs(w, m, x)
end

end

if !atomshell
    Blink.AtomShell.uninstall()
end


text:
using SortingLab
using Base.Test
using SortingAlgorithms
import SortingAlgorithms: load_bits

# write your own tests here
# @test 1 == 2
N = 1_000_000
K = 100
samplespace = "id".*dec.(1:N÷K,10);
srand(1)
svec = rand(samplespace, N);

using BenchmarkTools

# overall sorting is faster too
function ghi(svec)
    csvec = copy(svec)
    gc()
    @time a = @elapsed SortingLab.sort_spointer!(csvec, 1, length(csvec), StringRadixSort, Base.Forward)

    csvec = copy(svec)
    gc()
    @time b = @elapsed SortingAlgorithms.sort!(csvec, 1, length(csvec), StringRadixSort, Base.Forward)
    a,b
end
@time [ghi(svec) for i=repeat([4], outer=5)]

# sorting vectors shows sorttwo! is faster
function def(svec, skipbytes)
    vs = load_bits.(UInt, svec, skipbytes);
    csvec = copy(svec)
    gc()
    @time a = @elapsed SortingLab.sorttwo!(vs, csvec);

    vs = load_bits.(UInt, svec, skipbytes);
    csvec = copy(svec)
    gc()
    @time b = @elapsed SortingAlgorithms.sorttwo!(vs, csvec);
    a,b
end
@time [def(svec, i) for i=repeat([4], outer=5)]

# sorting pointers shows that sorttwo! is fasters
function abc(svec, skipbytes)
    vs = load_bits.(UInt, svec, skipbytes);
    # index = FastGroupBy.fcollect(length(vs));
    # csvec = copy(svec)
    csvec = pointer.(svec)
    gc()
    @time a = @elapsed SortingLab.sorttwo!(vs, csvec);

    vs = load_bits.(UInt, svec, skipbytes);
    # index = FastGroupBy.fcollect(length(vs));
    # csvec = copy(svec)
    csvec = pointer.(svec)
    gc()
    # b = @elapsed SortingLab.sorttwo_old!(vs, index);
    @time b = @elapsed SortingAlgorithms.sorttwo!(vs, csvec);
    a,b

    # csvec = copy(svec)
    # @time a = @elapsed SortingLab.sort_spointer!(csvec, 1, length(csvec), StringRadixSort, Base.Forward)

    # csvec = copy(svec)
    # @time b = @elapsed SortingAlgorithms.sort!(csvec, 1, length(csvec), StringRadixSort, Base.Forward)
    # a,b
end

using Base.Threads
println(nthreads());

@time [abc(svec, i) for i=repeat([4], outer=2)]

csvec = [randstring(8) for i =1:100];
@time SortingLab.sort_spointer!(csvec, 1, length(csvec), StringRadixSort, Base.Forward);

csvec = [randstring(8) for i =1:100];
@time SortingAlgorithms.sort!(csvec, 1, length(csvec), StringRadixSort, Base.Forward);

lo = 1
hi = length(svec)

# [abc(svec) for i  = 1:20]
# 20-element Array{Tuple{Float64,Float64},1}:
#  (3.30447, 2.27538)
#  (1.56587, 2.38168)
#  (1.72477, 2.29592)
#  (1.51554, 2.19927)
#  (1.58114, 2.37205)
#  (1.67712, 2.27833)
#  (1.54081, 2.25543)
#  (1.6349, 2.38872)
#  (1.62107, 2.27291)
#  (1.60027, 2.393)
#  (1.5623, 2.34241)
#  (1.64421, 2.36433)
#  (1.65756, 2.29755)
#  (1.58756, 2.27567)
#  (1.6318, 2.37143)
#  (1.6035, 2.31926)
#  (1.58667, 2.42479)
#  (1.6024, 2.33549)
#  (1.58388, 2.39197)
#  (1.64123, 2.29582)

function count_missing(A::Vector{Union{Bool,Missing}})
    t, m = 0, 0
    @inbounds for v in A
        ismissing(v) && (m += 1; continue)
        v && (t += 1)
    end
    return t, length(A) - t - m, m
end

a = [missing, true, false]
x = rand(a, 2^31-1)

function basic_bench(x)
    [@elapsed(count_missing(x)) for i = ]
end

using BenchmarkTools
@benchmark count_missing(x);


text:
# Fizz BUzz with Julia

for num in 1:100
    if num % 15 == 0
        println("FizzBuzz")
    elseif num % 3 == 0
        println("Fizz")
    elseif num % 5 == 0
        println("Buzz")
    else
        println(num)
    end
end

max_stars_repo_name: capturemathan/Fizz-Buzz

max_issues_repo_name: capturemathan/Fizz-Buzz

max_forks_repo_name: capturemathan/Fizz-Buzz


text:
using HorizonSideRobots
include("../basics.jl")

md"""
            task_3 c перегородками
"""

include("../task_3/solution.jl")

 """Решение задачи 16"""
 function solve_16!(r)
    solve_problem_3!(r)
end


text:
"""
Constructs a block diagonal matrix from blocks

    MakeBlockDiagonalMatrixR(
        model::SFFM.Model,
        mesh::DGMesh,
        Blocks,
        Factors::Array,
    )

# Aguments
- `model`: A Model object
- `mesh`: A DGMesh object from `DGMesh()`
- `Blocks(x::Array{Float64}, i::Int)`: a function wich returns a
    `NBases(mesh)×NBases(mesh)` array to put along the diagonal. The input
    argument `x` is a column vector of corresponding to the nodes in each cell,
    i.e. `SFFM.CellNodes(mesh)[:,n]`. `i` denotes the phase.
- `Factors::Array{<:Real,1}`: a `NIntervals(mesh)×1` vector of factors which multiply blocks

# Output
- `BlockMatrix::Array{Float64,2}`: `TotalNBases(mesh)×TotalNBases(mesh)` the
        block matrix
"""
function MakeBlockDiagonalMatrixR(
    model::SFFM.Model,
    mesh::DGMesh,
    Blocks,
    Factors::Array,
)
    BlockMatrix = Array{SparseArrays.SparseMatrixCSC{Float64,Int64},1}(undef,NPhases(model))
    for i in 1:NPhases(model)
        BlockMatrix[i] = SparseArrays.spzeros(Float64, TotalNBases(mesh), TotalNBases(mesh))
    end
    for i = 1:NIntervals(mesh), j = 1:NPhases(model)
        idx = (1:NBases(mesh)) .+ (i - 1) * NBases(mesh)
        BlockMatrix[j][idx, idx] = Blocks(SFFM.CellNodes(mesh)[:, i], j) * Factors[i]
    end
    return (BlockMatrix = BlockMatrix)
end

"""
Constructs the flux matrices for DG

    MakeFluxMatrixR(
        mesh::DGMesh,
        model::SFFM.Model,
        Phi,
    )

# Arguments
- `mesh`: a Mesh object from MakeMesh
- `model`: A Model object
- `Phi::Array{Float64,2}`: where `Phi[1,:]` and `Phi[1,:]` are the basis
    function evaluated at the left-hand and right-hand edge of a cell,
    respectively

# Output
- `F::Array{SparseArrays.SparseMatrixCSC{Float64,Int64},1}(undef,NPhases(model))`:
    an array with index `i ∈ 1:NPhases(model)`, of sparse arrays which are
    `TotalNBases×TotalNBases` flux matrices for phase `i`.
"""
function MakeFluxMatrixR(
    mesh::DGMesh,
    model::SFFM.Model,
    Phi,
)
    ## Create the blocks
    PosDiagBlock = -Phi[end, :] * Phi[end, :]'
    NegDiagBlock = Phi[1, :] * Phi[1, :]'
    UpDiagBlock = Phi[end, :] * Phi[1, :]'
    LowDiagBlock = -Phi[1, :] * Phi[end, :]'

    ## Construct global block diagonal matrix
    F = Array{SparseArrays.SparseMatrixCSC{Float64,Int64},1}(undef,NPhases(model))
    for i = 1:NPhases(model)
        F[i] = SparseArrays.spzeros(Float64, TotalNBases(mesh), TotalNBases(mesh))
        for k = 1:NIntervals(mesh)
            idx = (1:NBases(mesh)) .+ (k - 1) * NBases(mesh)
            if model.C[i] > 0
                xright = SFFM.CellNodes(mesh)[end, k]
                R = 1.0 ./ model.r.a(xright)[i]
                F[i][idx, idx] = PosDiagBlock * R
            elseif model.C[i] < 0
                xleft = SFFM.CellNodes(mesh)[1, k]
                R = 1.0 ./ model.r.a(xleft)[i]
                F[i][idx, idx] = NegDiagBlock * R
            end # end if C[i]
            if k > 1
                idxup = (1:NBases(mesh)) .+ (k - 2) * NBases(mesh)
                if model.C[i] > 0
                    xright = SFFM.CellNodes(mesh)[end, k-1]
                    R = 1.0 ./ model.r.a(xright)[i]
                    η = (Δ(mesh)[k] / NBases(mesh)) / (Δ(mesh)[k-1] / NBases(mesh))
                    F[i][idxup, idx] = UpDiagBlock * η * R
                elseif model.C[i] < 0
                    xleft = SFFM.CellNodes(mesh)[1, k]
                    R = 1.0 ./ model.r.a(xleft)[i]
                    η = (Δ(mesh)[k-1] / NBases(mesh)) / (Δ(mesh)[k] / NBases(mesh))
                    F[i][idx, idxup] = LowDiagBlock * η * R
                end # end if C[i]
            end # end if k>1
        end # for k in ...
    end # end for i in NPhases

    return (F = F)
end

"""
Creates the Local and global mass, stiffness and flux matrices to compute `D(s)`
directly.

    MakeMatricesR(
        model::SFFM.Model,
        mesh::DGMesh,
    )

# Arguments
- `model`: A Model object
- `mesh`: A mesh object

# Output
- A tuple of tuples
    - `:Global`: a tuple with fields
      - `:G::Array{SparseArrays.SparseMatrixCSC{Float64,Int64},1}`:
            a `NPhases` length array containing `TotalNBases×TotalNBases`
            dimensional global stiffness matrices
      - `:M::Array{SparseArrays.SparseMatrixCSC{Float64,Int64},1}`:
            a `NPhases` length array containing `TotalNBases×TotalNBases`
            dimensional global mass matrices
      - `:MInv::Array{SparseArrays.SparseMatrixCSC{Float64,Int64},1}`:
            a `NPhases` length array containing `TotalNBases×TotalNBases`
            dimensional matrices, the inverse of `Global.M`
      - `:F::Array{SparseArrays.SparseMatrixCSC{Float64,Int64},1}`:
            a `NPhases` length array containing `TotalNBases×TotalNBases`
            dimensional global flux matrices
      - `:Q::Array{SparseArrays.SparseMatrixCSC{Float64,Int64},1}`:
            a `NPhases` length array containing `TotalNBases×TotalNBases`
            dimensional global DG drift operator
    - `:Local`: a tuple with fields
      - `:G::Array{Float64,2}`: `NBases×NBases` Local stiffness matrix
      - `:M::Array{Float64,2}`: `NBases×NBases` Local mass matrix
      - `:MInv::Array{Float64,2}`: the inverse of `Local.M`
      - `:V::NamedTuple`: as output from SFFM.vandermonde
 """
function MakeMatricesR(
    model::SFFM.Model,
    mesh::DGMesh,
)
    ## Construct blocks
    V = vandermonde(NBases(mesh))
    if Basis(mesh) == "legendre"
        MLocal = function (x::Array{Float64}, i::Int)
            # Numerical integration of ϕᵢ(x)|r(x)|ϕⱼ(x) over Dk with Gauss-Lobatto
            # quadrature
            # Inputs:
            #   - x a vector of Gauss-Lobatto points on Dk
            #   - i a phase
            V.V' * LinearAlgebra.diagm(V.w ./ model.r.a(x)[:, i]) * V.V
        end
        GLocal = function (x::Array{Float64}, i::Int)
            # Numerical integration of ϕᵢ(x)|r(x)|ϕⱼ'(x) over Dk with Gauss-Lobatto
            # quadrature
            # Inputs:
            #   - x a vector of Gauss-Lobatto points on Dk
            #   - i a phase
            V.V' * LinearAlgebra.diagm(V.w ./ model.r.a(x)[:, i]) * V.D
        end
        MInvLocal = function (x::Array{Float64}, i::Int)
            MLocal(x, i)^-1
        end
        Phi = V.V[[1; end], :]

    elseif Basis(mesh) == "lagrange"
        MLocal = function (x::Array{Float64}, i::Int)
            # Numerical integration of ϕᵢ(x)|r(x)|ϕⱼ(x) over Dk with Gauss-Lobatto
            # quadrature
            # Inputs:
            #   - x a vector of Gauss-Lobatto points on Dk
            #   - i a phase
            LinearAlgebra.diagm(V.w ./ model.r.a(x)[:, i])
        end
        GLocal = function (x::Array{Float64}, i::Int)
            # Numerical integration of ϕᵢ(x)|r(x)|ϕⱼ'(x) over Dk with Gauss-Lobatto
            # quadrature
            # Inputs:
            #   - x a vector of Gauss-Lobatto points on Dk
            #   - i a phase
            MLocal(x, i) * V.D * V.inv
        end
        MInvLocal = function (x::Array{Float64}, i::Int)
            MLocal(x, i)^-1
        end
        Phi = (V.inv*V.V)[[1; end], :]
    end

    ## Assemble into block diagonal matrices
    G = SFFM.MakeBlockDiagonalMatrixR(
        model,
        mesh,
        GLocal,
        ones(NIntervals(mesh)),
    )
    M = SFFM.MakeBlockDiagonalMatrixR(
        model,
        mesh,
        MLocal,
        Δ(mesh) * 0.5,
    )
    MInv = SFFM.MakeBlockDiagonalMatrixR(
        model,
        mesh,
        MInvLocal,
        2.0 ./ Δ(mesh),
    )

    F = SFFM.MakeFluxMatrixR(mesh, model, Phi)

    ## Assemble the DG drift operator
    Q = Array{SparseArrays.SparseMatrixCSC{Float64,Int64},1}(undef,NPhases(model))
    for i = 1:NPhases(model)
        Q[i] = model.C[i] * (G[i] + F[i])
    end

    Local = (G = GLocal, M = MLocal, MInv = MInvLocal, V = V, Phi = Phi)
    Global = (G = G, M = M, MInv = MInv, F = F, Q = Q)
    out = (Local = Local, Global = Global)
    println("UPDATE: MatricesR object created with keys ", keys(out))
    println("UPDATE:    MatricesR[:",keys(out)[1],"] object created with keys ", keys(out[1]))
    println("UPDATE:    MatricesR[:",keys(out)[2],"] object created with keys ", keys(out[2]))
    return out
end

"""
Construct the operator `D(s)` directly.

    MakeDR(
        Matrices,
        MatricesR,
        model::SFFM.Model,
        mesh::DGMesh,
        B,
    )

# Arguments
- `Matrices`: a tuple as output from
    `Matrices = SFFM.MakeMatrices(model=model,mesh=mesh,probTransform=false)`
    note, must have `probTransform=false`.
- `MatricesR`: a tuple as output from
    `MatricesR = SFFM.MakeMatricesR(model=approxModel,mesh=mesh)`
- `model`: a Model object
- `mesh`: a Mesh object as constructed by `MakeMesh`

# Output
- a tuple with keys `:DDict` and `:DR`
    - `DDict::Dict{String,Function(s::Real)}`: a dictionary of functions. Keys are
      of the for `"ℓm"` where `ℓ,m∈{+,-}`. Values are functions with one argument.
      Usage is along the lines of `D["+-"](s=1)`.
    - `DR`: a function which evaluates `D(s)`,
      Usage is along the lines of `DR(s=1)`. Returns an array structured as
      `D = [D["++"](s=1) D["+-"](s=1); D["-+"](s=1) D["--"](s=1)]`.
"""
function MakeDR(
    Matrices,
    MatricesR,
    model::SFFM.Model,
    mesh::DGMesh,
    B,
)
    N₊ = sum(model.C .>= 0)
    N₋ = sum(model.C .<= 0)

    BigN = NPhases(model) * TotalNBases(mesh) + N₊ + N₋
    MR = SparseArrays.spzeros(Float64, NPhases(model) * TotalNBases(mesh), NPhases(model) * TotalNBases(mesh))
    Minv =
        SparseArrays.spzeros(Float64, NPhases(model) * TotalNBases(mesh), NPhases(model) * TotalNBases(mesh))
    FGR = SparseArrays.spzeros(Float64, NPhases(model) * TotalNBases(mesh), NPhases(model) * TotalNBases(mesh))
    for i = 1:NPhases(model)
        idx = ((i-1)*TotalNBases(mesh)+1:i*TotalNBases(mesh))
        MR[idx, idx] = MatricesR.Global.M[i]
        Minv[idx, idx] = Matrices.Global.MInv
        FGR[idx, idx] = MatricesR.Global.Q[i] * Matrices.Global.MInv
    end

    # Interior behaviour
    T = SparseArrays.kron(model.T, SparseArrays.sparse(LinearAlgebra.I,TotalNBases(mesh),TotalNBases(mesh)))
    BR = SparseArrays.spzeros(Float64, BigN, BigN)
    BR[(N₋+1):(end-N₊), (N₋+1):(end-N₊)] = MR * T * Minv + FGR

    # Boundary behaviour
    # Lower boundary
    # At boundary
    BR[1:N₋, 1:N₋] = (1.0./ model.r.a(mesh.Nodes[1])'.*model.T)[model.C.<=0, model.C.<=0]
    # Out of boundary
    idxup = ((1:NBases(mesh)).+TotalNBases(mesh)*(findall(model.C .> 0) .- 1)')[:] .+ N₋
    BR[1:N₋, idxup] = kron(
        (1.0./model.r.a(mesh.Nodes[1])'.*model.T)[model.C.<=0, model.C.>0],
        Matrices.Local.Phi[1, :]' * Matrices.Local.MInv * 2 ./ Δ(mesh)[1],
    )
    # Into boundary
    idxdown = ((1:NBases(mesh)).+TotalNBases(mesh)*(findall(model.C .<= 0) .- 1)')[:] .+ N₋
    BR[idxdown, 1:N₋] = LinearAlgebra.kron(
        LinearAlgebra.diagm(
            model.C[model.C.<=0] ./ model.r.a(mesh.Nodes[1].+sqrt(eps()))[model.C.<=0],
        ),
        -Matrices.Local.Phi[1, :],
    )

    # Upper boundary
    # At boundary
    BR[(end-N₊+1):end, (end-N₊+1):end] =
        (1.0./model.r.a(mesh.Nodes[end])'.*model.T)[model.C.>=0, model.C.>=0]
    # Out of boundary
    idxdown =
        ((1:NBases(mesh)).+TotalNBases(mesh)*(findall(model.C .< 0) .- 1)')[:] .+
        (N₋ + TotalNBases(mesh) - NBases(mesh))
    BR[(end-N₊+1):end, idxdown] = kron(
        (1.0./model.r.a(mesh.Nodes[end])'.*model.T)[model.C.>=0, model.C.<0],
        Matrices.Local.Phi[end, :]' * Matrices.Local.MInv * 2 ./ Δ(mesh)[end],
    )
    # Into boundary
    idxup =
        ((1:NBases(mesh)).+TotalNBases(mesh)*(findall(model.C .>= 0) .- 1)')[:] .+
        (N₋ + TotalNBases(mesh) - NBases(mesh))
    BR[idxup, (end-N₊+1):end] = LinearAlgebra.kron(
        LinearAlgebra.diagm(
            model.C[model.C.>=0] ./ model.r.a(mesh.Nodes[end].-sqrt(eps()))[model.C.>=0],
        ),
        Matrices.Local.Phi[end, :],
    )

    idx0 = [mesh.Fil["p0"]; repeat(mesh.Fil["0"]', NBases(mesh), 1)[:]; mesh.Fil["q0"]]
    bullet = [
        (mesh.Fil["p+"] .| mesh.Fil["p-"])
        repeat((mesh.Fil["+"] .| mesh.Fil["-"])', NBases(mesh), 1)[:]
        (mesh.Fil["q+"] .| mesh.Fil["q-"])
    ]
    MR = [
        LinearAlgebra.I(N₋) SparseArrays.spzeros(Float64, N₋, BigN - N₋)
        SparseArrays.spzeros(Float64, BigN - N₊ - N₋, N₋) MR SparseArrays.spzeros(Float64, BigN - N₊ - N₋, N₊)
        SparseArrays.spzeros(Float64, N₊, BigN - N₊) LinearAlgebra.I(N₊)
    ]
    Minv = [
        LinearAlgebra.I(N₋) SparseArrays.spzeros(Float64, N₋, BigN - N₋)
        SparseArrays.spzeros(Float64, BigN - N₊ - N₋, N₋) Minv SparseArrays.spzeros(Float64, BigN - N₊ - N₋, N₊)
        SparseArrays.spzeros(Float64, N₊, BigN - N₊) LinearAlgebra.I(N₊)
    ]

    # BR[idx0, :] = B.B[idx0, :]
    DR = function (;s=0)
        BR[bullet, bullet] -
        MR[bullet, bullet] * s * SparseArrays.sparse(LinearAlgebra.I,sum(bullet),sum(bullet)) * Minv[bullet, bullet] +
        BR[bullet, idx0] *
        Matrix(s * SparseArrays.sparse(LinearAlgebra.I,sum(idx0),sum(idx0)) - B.B[idx0, idx0])^-1 *
        B.B[idx0, bullet]
    end
    # D = function (s)
    #     MR[bullet, bullet] *
    #     (T[bullet, bullet] - s * LinearAlgebra.I(sum(bullet))) *
    #     Minv[bullet, bullet] +
    #     FGR[bullet, bullet] +
    #     MR[bullet, bullet] *
    #     T[bullet, idx0] *
    #     (B.BDict["00"] - s * LinearAlgebra.I(sum(idx0)))^-1 *
    #     T[idx0, bullet]
    # end

    DDict = Dict{String,Any}()
    for ℓ in ["+", "-"], m in ["+", "-"]
        FlBases = mesh.Fil[ℓ][.!mesh.Fil["0"]]
        FmBases = mesh.Fil[m][.!mesh.Fil["0"]]
        FlBases = [
            mesh.Fil["p"*ℓ][.!mesh.Fil["p0"]]
            repeat(FlBases', NBases(mesh), 1)[:]
            mesh.Fil["q"*ℓ][.!mesh.Fil["q0"]]
        ]
        FmBases = [
            mesh.Fil["p"*m][.!mesh.Fil["p0"]]
            repeat(FmBases', NBases(mesh), 1)[:]
            mesh.Fil["q"*m][.!mesh.Fil["q0"]]
        ]
        DDict[ℓ*m] = function (; s = 0)#::Array{Float64}
            DR(s=s)[FlBases, FmBases]
        end # end function
    end # end for ℓ, m ...
    out = (DDict = DDict, DR = DR)
    println("UPDATE: D(s) operator created with keys ", keys(out))
    println("UPDATE:    DDict[key](s) operator created with keys ", keys(DDict))
    println("UPDATE:    DR(s) operator created (a matrix)")
    return out
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM