text:
import Splittings: Geometry
using LinearAlgebra
import LinearAlgebra.LAPACK: pttrf!, pttrs!

mutable struct SplineNN

    geom       :: Geometry
    a1x        :: Float64
    a2x        :: Float64
    a3x        :: Float64
    a4x        :: Float64
    a1y        :: Float64
    a2y        :: Float64
    a3y        :: Float64
    a4y        :: Float64
    axd        :: Vector{Float64}
    ayd        :: Vector{Float64}
    axod       :: Vector{Float64}
    ayod       :: Vector{Float64}
    aym1gamma1 :: Vector{Float64}
    aym1gamma2 :: Vector{Float64}
    axm1gamma1 :: Vector{Float64}
    axm1gamma2 :: Vector{Float64}
    coef       :: Array{Float64,2}
    bcoef      :: Array{Float64,2}

    function SplineNN(geom::Geometry)

        n1   = geom.n1
        n2   = geom.n2
        n1p1 = n1+1
        n1p2 = n1+2
        n2p1 = n2+1
        n2p2 = n2+2

	delta1 = geom.delta1
	delta2 = geom.delta2

        axm1gamma1 = zeros(Float64, n1)
        axm1gamma2 = zeros(Float64, n1)
        aym1gamma1 = zeros(Float64, n2)
        aym1gamma2 = zeros(Float64, n2)
        axd        = 4ones(Float64, n1)
        axod       =  ones(Float64, n1-1)
        ayd        = 4ones(Float64, n2)
        ayod       =  ones(Float64, n2-1)

        coef       = zeros(Float64,(n1p2,n2p2))
        bcoef      = zeros(Float64,(n2p2,n1p2))

	pttrf!(axd, axod)
	pttrf!(ayd, ayod)

        axm1gamma2[1]  = 1
        axm1gamma1[n1] = 1

        pttrs!(axd, axod, axm1gamma1)
        pttrs!(axd, axod, axm1gamma2)

        # compute Ay-1.gamma
        aym1gamma2[1]  = 1
        aym1gamma1[n2] = 1

        pttrs!(ayd, ayod, aym1gamma1)
        pttrs!(ayd, ayod, aym1gamma2)

        aa1x = 3 /  delta1
        aa1y = 3 /  delta2
        aa2x = 6 / (delta1*delta1)
        aa2y = 6 / (delta2*delta2)

        # assemblage de la matrice 2x2 pour la spline dans la direction Ox

        a1x = aa2x*(1. -   axm1gamma1[n1-1] + 2*axm1gamma1[n1])
        a2x = aa2x*(   -   axm1gamma2[n1-1] + 2*axm1gamma2[n1])
        a3x = aa2x*(     2*axm1gamma1[1]    -   axm1gamma1[2])
        a4x = aa2x*(1. + 2*axm1gamma2[1]    -   axm1gamma2[2])

        # assemblage de la matrice 2x2 pour spline naturels (direction Oy)

        a1y = aa2y*(1. -   aym1gamma1[n2-1] + 2*aym1gamma1[n2])
        a2y = aa2y*(   -   aym1gamma2[n2-1] + 2*aym1gamma2[n2])
        a3y = aa2y*(     2*aym1gamma1[1]    -   aym1gamma1[2])
        a4y = aa2y*(1. + 2*aym1gamma2[1]    -   aym1gamma2[2])

        new( geom, a1x, a2x, a3x, a4x, a1y, a2y, a3y, a4y,
	     axd , ayd , axod , ayod,
             aym1gamma1, aym1gamma2, axm1gamma1, axm1gamma2,
             coef, bcoef )

    end

end

function interpolate!(spline :: SplineNN,
                      f      :: Array{Float64,2},
		      x1     :: Array{Float64,2},
		      x2     :: Array{Float64,2})

    nat_x!(spline, f)
    nat_y!(spline)
    evaltab!(spline, x1, x2, f)

end

"""
interpolation par spline periodique dans les deux directions.
Les points d'interpolation sont definis grace a depx et depy
qui definissent le deplacement par rapport au maillage.
 - f contient les valeurs de la fonction de distribution
 - depx et depy : deplacements par rapport au maillage
     des points dans les quels on veut evaluer la spline.
"""
function interpolate!(spline :: SplineNN,
                      f      :: Array{Float64,2},
                      depx   :: Float64,
                      depy   :: Float64)

    nat_x!(spline, f)
    nat_y!(spline)
    evaldep!(spline, depx, depy, f)

end

"""
    natural splines

"""
function nat_x!( spline :: SplineNN,
                 gtau   :: Array{Float64,2} )

    n1     = spline.geom.n1
    n2     = spline.geom.n2
    delta1 = spline.geom.delta1
    delta2 = spline.geom.delta2

    axm1f = similar(gtau)

    n1p2 = n1+2
    n1p1 = n1+1
    det  = spline.a1x*spline.a4x - spline.a2x*spline.a3x

    for j=1:n2
       for i=1:n1
          axm1f[i,j] = 6 * gtau[i,j]
       end
    end

    pttrs!(spline.axd,spline.axod,axm1f)

    for j=1:n2
       # assemblage du second membre du systeme 2x2
       gamma1 = (6.0/(delta1)^2)*( -axm1f[n1-1,j] + 2*axm1f[n1,j])
       gamma2 = (6.0/(delta1)^2)*(2*axm1f[   1,j] -   axm1f[ 2,j])

       coefnp2 = (gamma1*spline.a4x - gamma2*spline.a2x)/det
       coef1   = (gamma2*spline.a1x - gamma1*spline.a3x)/det
       spline.bcoef[j,n1p2] = coefnp2
       spline.bcoef[j,   1] = coef1

       for  i=2:n1p1
           spline.bcoef[j,i] = axm1f[i-1,j]
                             - spline.axm1gamma1[i-1]*coefnp2
                             - spline.axm1gamma2[i-1]*coef1
       end
    end
end

function nat_y!(spline :: SplineNN)

    n1     = spline.geom.n1
    n2     = spline.geom.n2
    delta1 = spline.geom.delta1
    delta2 = spline.geom.delta2

    n1p2 = n1+2
    n2p2 = n2+2
    det  = spline.a1y*spline.a4y - spline.a2y*spline.a3y

    aym1f = zeros(Float64, (n2,n1p2))
    for  i=1:n1p2
       for j=1:n2
          aym1f[j,i] = 6*spline.bcoef[j,i]
       end
    end

    pttrs!(spline.ayd,spline.ayod,aym1f)

    for i=1:n1p2

       gamma1 = (6.0/(delta2)^2)*( -aym1f[n2-1,i] + 2*aym1f[n2,i])
       gamma2 = (6.0/(delta2)^2)*(2*aym1f[   1,i] -   aym1f[ 2,i])

       coefnp2 = (gamma1*spline.a4y - gamma2*spline.a2y)/det
       coef1   = (gamma2*spline.a1y - gamma1*spline.a3y)/det

       spline.coef[i,n2p2] = coefnp2
       spline.coef[i,   1] = coef1
       for  j=2:n2 + 1
          spline.coef[i,j] = aym1f[j-1,i]
                           - spline.aym1gamma1[j-1]*coefnp2
                           - spline.aym1gamma2[j-1]*coef1
       end

    end
end

function evaltab!(spline :: SplineNN,
                  xd     :: Array{Float64,2},
                  yd     :: Array{Float64,2},
                  f      :: Array{Float64,2})

    n1         = spline.geom.n1
    n2         = spline.geom.n2
    delta1     = spline.geom.delta1
    delta2     = spline.geom.delta2
    delta1x    = delta1*delta1
    delta1xx   = delta1x*delta1
    delta1xx6  = 1/(6*delta1xx)
    delta2y    = delta2*delta2
    delta2yy   = delta2y*delta2
    delta2yy6  = 1/(6*delta2yy)
    idelta1    = 1/delta1
    idelta2    = 1/delta2

    for j = 2:n2-1
       for i = 2:n1-1

          i1 = trunc(Int64,(xd[i,j]-spline.geom.x1min)*idelta1)
          j1 = trunc(Int64,(yd[i,j]-spline.geom.x2min)*idelta2)

          xdp1   = spline.geom.x1grid[i1+2]-xd[i,j]
          bvalx1 = xdp1*xdp1*xdp1
          bvalx2 = delta1xx+3*delta1x*xdp1
                 + 3*delta1*xdp1*xdp1-3*xdp1*xdp1*xdp1
          xd1    = xd[i,j]-spline.geom.x1grid[i1+1]
          bvalx3 = delta1xx+3*delta1x*xd1
                 + 3*delta1*xd1*xd1-3*xd1*xd1*xd1
          bvalx4 = xd1*xd1*xd1

          ydp1   = spline.geom.x2grid[j1+2]-yd[i,j]
          bvaly1 = ydp1*ydp1*ydp1
          bvaly2 = delta2yy+3*ydp1*(delta2y+ydp1*(delta2-ydp1))
          yd1    = yd[i,j]-geom.x2grid[j1+1]
          bvaly3 = delta2yy+3*yd1*(delta2y+yd1*(delta2-yd1))
          bvaly4 = yd1*yd1*yd1

          sval  = 0.
          sval1 = spline.coef[i1+1,j1+1]*bvaly1
          sval1 = sval1+spline.coef[i1+1,j1+2]*bvaly2
          sval1 = sval1+spline.coef[i1+1,j1+3]*bvaly3
          sval1 = sval1+spline.coef[i1+1,j1+4]*bvaly4

          sval  = sval+sval1*bvalx1

          sval2 = spline.coef[i1+2,j1+1]*bvaly1
          sval2 = sval2+spline.coef[i1+2,j1+2]*bvaly2
          sval2 = sval2+spline.coef[i1+2,j1+3]*bvaly3
          sval2 = sval2+spline.coef[i1+2,j1+4]*bvaly4

          sval  = sval+sval2*bvalx2

          sval3 = spline.coef[i1+3,j1+1]*bvaly1
          sval3 = sval3+spline.coef[i1+3,j1+2]*bvaly2
          sval3 = sval3+spline.coef[i1+3,j1+3]*bvaly3
          sval3 = sval3+spline.coef[i1+3,j1+4]*bvaly4

          sval  = sval+sval3*bvalx3

          sval4 = spline.coef[i1+4,j1+1]*bvaly1
          sval4 = sval4+spline.coef[i1+4,j1+2]*bvaly2
          sval4 = sval4+spline.coef[i1+4,j1+3]*bvaly3
          sval4 = sval4+spline.coef[i1+4,j1+4]*bvaly4

          sval  = sval+sval4*bvalx4

          f[i,j] = delta1xx6*delta2yy6*sval

       end
    end

    fill!(view(f, 1, :), 0)
    fill!(view(f,n1, :), 0)

    fill!(view(f,:,  1), 0)
    fill!(view(f,:, n2), 0)

end

function evaldep!(spline::SplineNN,
		  alphax::Float64,
		  alphay::Float64,
		  f::Array{Float64,2})

    n1        = spline.geom.n1
    n2        = spline.geom.n2

    delta1    = spline.geom.delta1
    delta1x   = delta1*delta1
    delta1xx  = delta1x*delta1
    delta1xx6 = 1/(6*delta1xx)

    delta2    = spline.geom.delta2
    delta2y   = delta2*delta2
    delta2yy  = delta2y*delta2
    delta2yy6 = 1/(6*delta2yy)

    if (alphax > 0)
       intaxsdx=trunc(Int64,-alphax/delta1+eps(Float64))-1
       ideb = trunc(Int64,alphax/delta1)+2
       ifin = n1 - 1
    else
       intaxsdx=trunc(Int64,-alphax/delta1)
       ideb = 2
       ifin = -trunc(Int64,-alphax/delta1)+n1-1
    end

    if (alphay>0)
       intaysdy=trunc(Int64,-alphay/delta2+eps(Float64))-1
       jdeb = trunc(Int64,alphay/delta2)+2
       jfin = n2 - 1
    else
       intaysdy=trunc(Int64,-alphay/delta2)
       jdeb = 2
       jfin = -trunc(Int64,-alphay/delta2)+n2-1
    end

    xd1    = -alphax-intaxsdx*delta1
    xdp1   =  delta1-xd1
    yd1    = -alphay-intaysdy*delta2
    ydp1   =  delta2-yd1
    bvalx1 = xdp1*xdp1*xdp1
    bvalx2 = delta1xx+3*delta1x*xdp1+3*delta1*xdp1*xdp1-3*xdp1*xdp1*xdp1
    bvalx3 = delta1xx+3*delta1x*xd1+3*delta1*xd1*xd1-3*xd1*xd1*xd1
    bvalx4 = xd1*xd1*xd1
    bvaly1 = ydp1*ydp1*ydp1
    bvaly2 = delta2yy+3*ydp1*(delta2y+ydp1*(delta2-ydp1))
    bvaly3 = delta2yy+3*yd1*(delta2y+yd1*(delta2-yd1))
    bvaly4 = yd1*yd1*yd1

    for j=jdeb:jfin
       j1=j-1+intaysdy
       for i=ideb:ifin
          i1=i-1+intaxsdx

          f[i,j] = delta1xx6*delta2yy6* (
                 bvalx1 * (
                 spline.coef[i1+1,j1+1] * bvaly1
               + spline.coef[i1+1,j1+2] * bvaly2
               + spline.coef[i1+1,j1+3] * bvaly3
               + spline.coef[i1+1,j1+4] * bvaly4)
               + bvalx2 * (
                 spline.coef[i1+2,j1+1] * bvaly1
               + spline.coef[i1+2,j1+2] * bvaly2
               + spline.coef[i1+2,j1+3] * bvaly3
               + spline.coef[i1+2,j1+4] * bvaly4)
               + bvalx3 * (
                 spline.coef[i1+3,j1+1] * bvaly1
               + spline.coef[i1+3,j1+2] * bvaly2
               + spline.coef[i1+3,j1+3] * bvaly3
               + spline.coef[i1+3,j1+4] * bvaly4)
               + bvalx4 * (
                 spline.coef[i1+4,j1+1] * bvaly1
               + spline.coef[i1+4,j1+2] * bvaly2
               + spline.coef[i1+4,j1+3] * bvaly3
               + spline.coef[i1+4,j1+4] * bvaly4))
       end
    end

    fill!(view(f,1:ideb-1,:),0)
    fill!(view(f,ifin+1:n1,:),0)

    fill!(view(f,:,1:jdeb-1),0)
    fill!(view(f,:,jfin+1:n2),0)

end


text:
# check whether a reaction on a chemical hyperedge is uni- or multimolecular
isunimolecular(che::ChemicalHyperEdge) = nsrcs(che) <=1 && ntgts(che) <= 1
ismultimolecular(che::ChemicalHyperEdge) = !(isunimolecular(che))

# return two reactions (forward and backward) given some species and two rates
function reversible_reaction(s1, s2, rates::AbstractVector)
    length(rates) == 2 || error("too many rates given: expected 2, given $(length(rates))")
    fwd = ChemicalHyperEdge(s1, s2, rates[1])
    bwd = switch(fwd); update_weight!(bwd, rates[2])
    return [fwd, bwd]
end
reversible_reaction(che::ChemicalHyperEdge, rates::AbstractVector) = reversible_reaction(src(che), tgt(che), rates)

# macro to easily use the above
macro reversible(che, rates)
    return esc(:(reversible_reaction($che, $rates)))
end

ext: jl


text:
using JuliaSyntax
using Test

using JuliaSyntax: SourceFile

using JuliaSyntax: GreenNode, SyntaxNode,
    flags, EMPTY_FLAGS, TRIVIA_FLAG, INFIX_FLAG,
    children, child, setchild!, SyntaxHead

module TokenizeTests
    using Test
    @testset "Tokenize" begin
        include("../Tokenize/test/runtests.jl")
    end
end

include("test_utils.jl")
include("parse_stream.jl")
include("parser.jl")
include("parser_api.jl")
include("syntax_tree.jl")

@testset "Parsing values from strings" begin
    include("value_parsing.jl")
end

include("parse_packages.jl")

# Prototypes
#include("syntax_interpolation.jl")
#include("simple_parser.jl")


text: derivative(f, point, h) = imag(f(point + im*h))/h


text:
@enum ExpressionOrder LOWEST EQUALS LESSGREATER SUM PRODUCT PREFIX CALL INDEX

const ORDERS = Dict{TokenType,ExpressionOrder}(
    EQ => EQUALS,
    NOT_EQ => EQUALS,
    LT => LESSGREATER,
    GT => LESSGREATER,
    PLUS => SUM,
    MINUS => SUM,
    SLASH => PRODUCT,
    ASTERISK => PRODUCT,
    LPAREN => CALL,
    LBRACKET => INDEX,
)

mutable struct Parser
    l::Lexer
    errors::Vector{ErrorObj}
    cur_token::Token
    peek_token::Token
    prefix_parse_functions::Dict{TokenType,Function}
    infix_parse_functions::Dict{TokenType,Function}
end

function Parser(l::Lexer)
    cur_token = next_token!(l)
    peek_token = next_token!(l)
    p = Parser(l, [], cur_token, peek_token, Dict(), Dict())

    register_prefix!(p, IDENT, parse_identifier!)
    register_prefix!(p, INT, parse_integer_literal!)
    register_prefix!(p, BANG, parse_prefix_expression!)
    register_prefix!(p, MINUS, parse_prefix_expression!)
    register_prefix!(p, TRUE, parse_boolean!)
    register_prefix!(p, FALSE, parse_boolean!)
    register_prefix!(p, NULL, parse_null!)
    register_prefix!(p, LPAREN, parse_grouped_expression!)
    register_prefix!(p, IF, parse_if_expression!)
    register_prefix!(p, FUNCTION, parse_function_literal!)
    register_prefix!(p, STRING, parse_string_literal!)
    register_prefix!(p, LBRACKET, parse_array_literal!)
    register_prefix!(p, LBRACE, parse_hash_literal!)
    register_prefix!(p, MACRO, parse_macro_literal!)

    register_infix!(p, PLUS, parse_infix_expression!)
    register_infix!(p, MINUS, parse_infix_expression!)
    register_infix!(p, SLASH, parse_infix_expression!)
    register_infix!(p, ASTERISK, parse_infix_expression!)
    register_infix!(p, EQ, parse_infix_expression!)
    register_infix!(p, NOT_EQ, parse_infix_expression!)
    register_infix!(p, LT, parse_infix_expression!)
    register_infix!(p, GT, parse_infix_expression!)
    register_infix!(p, LPAREN, parse_call_expression!)
    register_infix!(p, LBRACKET, parse_index_expression!)

    return p
end

register_prefix!(p::Parser, t::TokenType, fn::Function) = p.prefix_parse_functions[t] = fn
register_infix!(p::Parser, t::TokenType, fn::Function) = p.infix_parse_functions[t] = fn

function next_token!(p::Parser)
    p.cur_token = p.peek_token
    p.peek_token = next_token!(p.l)
end

function parse_statement!(p::Parser)
    if p.cur_token.type == LET || (p.cur_token.type == IDENT && p.peek_token.type == ASSIGN)
        return parse_let_statement!(p)
    elseif p.cur_token.type == RETURN
        return parse_return_statement!(p)
    elseif p.cur_token.type == WHILE
        return parse_while_statement!(p)
    elseif p.cur_token.type == BREAK
        return parse_break_statement!(p)
    elseif p.cur_token.type == CONTINUE
        return parse_continue_statement!(p)
    else
        return parse_expression_statement!(p)
    end
end

function parse_let_statement!(p::Parser)
    token = p.cur_token
    reassign = token.type == IDENT

    if !reassign && !expect_peek!(p, IDENT)
        return nothing
    end

    name = Identifier(p.cur_token, p.cur_token.literal)

    if !expect_peek!(p, ASSIGN)
        return nothing
    end

    next_token!(p)

    value = parse_expression!(p, LOWEST)

    if isa(value, FunctionLiteral)
        value =
            FunctionLiteral(value.token, value.parameters, value.body; name = name.value)
    end

    if p.peek_token.type == SEMICOLON
        next_token!(p)
    end

    return LetStatement(token, name, value, reassign)
end

function parse_return_statement!(p::Parser)
    token = p.cur_token

    next_token!(p)

    value = parse_expression!(p, LOWEST)

    if p.peek_token.type == SEMICOLON
        next_token!(p)
    end

    return ReturnStatement(token, value)
end

function parse_while_statement!(p::Parser)
    token = p.cur_token

    if !expect_peek!(p, LPAREN)
        return nothing
    end

    next_token!(p)
    condition = parse_expression!(p, LOWEST)

    if !expect_peek!(p, RPAREN)
        return nothing
    end

    if !expect_peek!(p, LBRACE)
        return nothing
    end

    body = parse_block_statement!(p; must_close = true)

    return WhileStatement(token, condition, body)
end

function parse_break_statement!(p::Parser)
    token = p.cur_token

    if p.peek_token.type == SEMICOLON
        next_token!(p)
    end

    return BreakStatement(token)
end

function parse_continue_statement!(p::Parser)
    token = p.cur_token

    if p.peek_token.type == SEMICOLON
        next_token!(p)
    end

    return ContinueStatement(token)
end

function parse_expression_statement!(p::Parser)
    token = p.cur_token
    expression = parse_expression!(p, LOWEST)
    if p.peek_token.type == SEMICOLON
        next_token!(p)
    end

    return ExpressionStatement(token, expression)
end

function parse_block_statement!(p::Parser; must_close = false)
    token = p.cur_token
    statements = Statement[]

    next_token!(p)

    while p.cur_token.type != RBRACE && p.cur_token.type != EOF
        statement = parse_statement!(p)

        if !isnothing(statement)
            push!(statements, statement)
        end

        next_token!(p)
    end

    if must_close && p.cur_token.type != RBRACE
        push!(p.errors, ErrorObj("parse error: braces must be closed"))
        return nothing
    end

    return BlockStatement(token, statements)
end

function parse_expression!(p::Parser, order::ExpressionOrder)
    if p.cur_token.type ∉ keys(p.prefix_parse_functions)
        push!(
            p.errors,
            ErrorObj("parse error: no prefix parse function for $(p.cur_token.type) found"),
        )
        return nothing
    else
        prefix_fn = p.prefix_parse_functions[p.cur_token.type]
        left = prefix_fn(p)

        while p.peek_token.type != SEMICOLON && order < peek_order(p)
            if p.peek_token.type ∉ keys(p.infix_parse_functions)
                return left
            end

            infix_fn = p.infix_parse_functions[p.peek_token.type]
            next_token!(p)
            left = infix_fn(p, left)
        end

        return left
    end
end

parse_boolean!(p::Parser) = BooleanLiteral(p.cur_token, p.cur_token.type == TRUE)

parse_null!(p::Parser) = NullLiteral(p.cur_token)

parse_identifier!(p::Parser) = Identifier(p.cur_token, p.cur_token.literal)

function parse_integer_literal!(p::Parser)
    token = p.cur_token
    try
        value = Base.parse(Int64, p.cur_token.literal)
        return IntegerLiteral(token, value)
    catch
        push!(
            p.errors,
            ErrorObj("parse error: could not parse $(p.cur_token.literal) as integer"),
        )
        return nothing
    end
end

function parse_string_literal!(p::Parser)
    return StringLiteral(p.cur_token, p.cur_token.literal)
end

function parse_array_literal!(p::Parser)
    token = p.cur_token
    elements = parse_expression_list!(p, RBRACKET)
    return ArrayLiteral(token, elements)
end

function parse_hash_literal!(p::Parser)
    token = p.cur_token
    pairs = Dict{Expression,Expression}()

    while p.peek_token.type != RBRACE
        next_token!(p)
        key = parse_expression!(p, LOWEST)
        if !expect_peek!(p, COLON)
            return nothing
        end

        next_token!(p)
        val = parse_expression!(p, LOWEST)
        pairs[key] = val

        if p.peek_token.type != RBRACE && !expect_peek!(p, COMMA)
            return nothing
        end
    end

    next_token!(p)

    return HashLiteral(token, pairs)
end

function parse_index_expression!(p::Parser, left::Expression)
    token = p.cur_token
    next_token!(p)
    index = parse_expression!(p, LOWEST)

    if !expect_peek!(p, RBRACKET)
        return nothing
    end

    return IndexExpression(token, left, index)
end

function parse_prefix_expression!(p::Parser)
    token = p.cur_token
    operator = p.cur_token.literal
    next_token!(p)
    right = parse_expression!(p, PREFIX)
    return PrefixExpression(token, operator, right)
end

function parse_infix_expression!(p::Parser, left::Expression)
    token = p.cur_token
    operator = p.cur_token.literal

    order = cur_order(p)
    next_token!(p)
    right = parse_expression!(p, order)

    return InfixExpression(token, left, operator, right)
end

function parse_grouped_expression!(p::Parser)
    next_token!(p)
    expr = parse_expression!(p, LOWEST)
    if !expect_peek!(p, RPAREN)
        return nothing
    end
    return expr
end

function parse_expression_list!(p::Parser, end_token::TokenType)
    expressions = Expression[]

    if p.peek_token.type == end_token
        next_token!(p)
        return expressions
    end

    next_token!(p)
    push!(expressions, parse_expression!(p, LOWEST))

    while p.peek_token.type == COMMA
        next_token!(p)
        next_token!(p)
        push!(expressions, parse_expression!(p, LOWEST))
    end

    if !expect_peek!(p, end_token)
        return nothing
    end

    return expressions
end

function parse_if_expression!(p::Parser)
    token = p.cur_token

    if !expect_peek!(p, LPAREN)
        return nothing
    end

    next_token!(p)
    condition = parse_expression!(p, LOWEST)

    if !expect_peek!(p, RPAREN)
        return nothing
    end

    if !expect_peek!(p, LBRACE)
        return nothing
    end

    consequence = parse_block_statement!(p; must_close = true)

    if p.peek_token.type == ELSE
        next_token!(p)

        if !expect_peek!(p, LBRACE)
            return nothing
        end

        alternative = parse_block_statement!(p; must_close = true)
    else
        alternative = nothing
    end

    return IfExpression(token, condition, consequence, alternative)
end

function parse_function_or_macro_parameters!(p::Parser)
    identifiers = Identifier[]

    if p.peek_token.type == RPAREN
        next_token!(p)
        return identifiers
    end

    next_token!(p)
    ident = Identifier(p.cur_token, p.cur_token.literal)
    push!(identifiers, ident)

    while p.peek_token.type == COMMA
        next_token!(p)
        next_token!(p)
        ident = Identifier(p.cur_token, p.cur_token.literal)
        push!(identifiers, ident)
    end

    if !expect_peek!(p, RPAREN)
        return nothing
    end

    return identifiers
end

function parse_function_or_macro_literal!(p::Parser)
    token = p.cur_token

    if !expect_peek!(p, LPAREN)
        return nothing
    end

    parameters = parse_function_or_macro_parameters!(p)

    if !expect_peek!(p, LBRACE)
        return nothing
    end

    body = parse_block_statement!(p; must_close = true)
    if isnothing(body)
        return nothing
    end

    return (token, parameters, body)
end

parse_function_literal!(p::Parser) = begin
    result = parse_function_or_macro_literal!(p)
    return isnothing(result) ? nothing : FunctionLiteral(result...)
end

parse_macro_literal!(p::Parser) = begin
    result = parse_function_or_macro_literal!(p)
    return isnothing(result) ? nothing : MacroLiteral(result...)
end

function parse_call_expression!(p::Parser, fn::Expression)
    token = p.cur_token
    arguments = parse_expression_list!(p, RPAREN)

    return CallExpression(token, fn, arguments)
end

function expect_peek!(p::Parser, t::TokenType)
    if p.peek_token.type == t
        next_token!(p)
        return true
    else
        peek_error!(p, t)
        return false
    end
end

function peek_error!(p::Parser, t::TokenType)
    push!(
        p.errors,
        ErrorObj(
            "parse error: expected next token to be $t, got $(p.peek_token.type) instead",
        ),
    )
end

function cur_order(p::Parser)
    if p.cur_token.type ∉ keys(ORDERS)
        return LOWEST
    else
        return ORDERS[p.cur_token.type]
    end
end

function peek_order(p::Parser)
    if p.peek_token.type ∉ keys(ORDERS)
        return LOWEST
    else
        return ORDERS[p.peek_token.type]
    end
end

function parse!(p::Parser)
    program = Program([])

    while p.cur_token.type != EOF
        try
            statement = parse_statement!(p)
            if !isnothing(statement)
                push!(program.statements, statement)
            end
        catch
        finally
            next_token!(p)
        end
    end

    return program
end

function parse(code::String; input::IO = stdin, output::IO = stdout)
    if isempty(code)
        return nothing
    end

    lexer = Lexer(code)
    parser = Parser(lexer)
    program = parse!(parser)

    if !isempty(parser.errors)
        println(
            output,
            ErrorObj(
                "parser has $(length(parser.errors)) error$(length(parser.errors) == 1 ? "" : "s")",
            ),
        )
        println(output, join(map(string, parser.errors), "\n"))
        return nothing
    end

    return program
end


text:
for i in 0:9
  println(i)
end

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
function cutting_planes(model::JuMP.Model, VecBin::Vector{Int}, MaxIter::Int64 = 5)
    tol = 1e-5

    Astd, b, cstd, xlb, xub, solver, m, n, flag_sense = extract_data(model)

    iter = 1
    convergence = -2 #not finished
    X = 0
    z = 0

    while convergence == -2
        t_1 , t_2, t_3 = solve_LP(Astd, b, cstd, xlb, xub, solver)

        if flag_sense == 1
            if isnan(t_2)
                return 0, -z
            end
        else
            if isnan(t_2)
                return 0, z
            end
        end

        X, z, status = t_1 , t_2, t_3

        # Finding fractional variable address
        j = find(min.(X[1:n]-floor.(X[1:n]),ceil.(X[1:n])-X[1:n]).>tol)
        # Removing index from continous variables
        tmp = min(size(VecBin)[1], size(j)[1])
        tmp_ind = []
        for i in 1:tmp
            if size(VecBin)[1] <= size(j)[1]
                if VecBin[i] in j
                    push!(tmp_ind,VecBin[i])
                end
            else
                if j[i] in VecBin
                    push!(tmp_ind,j[i])
                end
            end
        end
        j = tmp_ind

        #Test convergence
        if length(j) == 0

            convergence = 1
        elseif iter >= MaxIter
            j = j[1]
            convergence  = -1
        else
            j = j[1] #pick a fractional and binary variable address (the first one)

            NBas = find(X .<= tol)
            Bas = deleteat!(collect(1:n), NBas)

            # u = Linha do dicionario correspondente a variavel x_j
            u = find(Bas.== j)

            DB = X[Bas]
            DN = ((Astd[:,Bas]'*Astd[:,Bas])\Astd[:,Bas]')*Astd[:,NBas]

            a0 = X[Bas]
            a = zeros(size(Bas)[1],n)

            a[:,NBas] = DN
            for i = 1:size(Bas)[1]
              a[i,Bas[i]] = 1
            end
            a[find(abs.(a) .< tol)] = 0

            N1 = intersect(deleteat!(copy(VecBin),find(VecBin .== j)),NBas)
            N2 = intersect(deleteat!(collect(1:n), VecBin),NBas)

            f0 = a0[u] - floor.(a0[u])
            f0 = f0[1]

            f = a[u,union(N1,N2)] - floor.(a[u,union(N1,N2)])

            A_cut = zeros(1,n)

            A_cut[j] = 1
            if size(intersect(find( f .<= f0),N1))[1] != 0
                A_cut[intersect(find( f .<= f0),N1)] = floor.(a[u,intersect(find( f .<= f0),N1)])
            end
            if size(intersect(find( f .>  f0),N1))[1] != 0
                A_cut[intersect(find( f .>  f0),N1)] = floor.(a[u,intersect(find( f .>  f0),N1)]) + (f[intersect(find( f .>  f0),N1)] - f0)/(1 - f0)
            end
            if size(intersect(find(a[u,:] .< 0), N2))[1] != 0
                A_cut[intersect(find(a[u,:] .< 0), N2)] = a[intersect(find(a[u,:] .< 0), N2)]*(1/(1-f0))
            end

            b_cut = floor.(a0[u])

            Astd = [[Astd zeros(m,1)]; [A_cut 1]]
            b = [b ; b_cut]

            cstd = [cstd ; 0]
            xlb = [xlb; 0]
            xub = [xub; +Inf]
            m+=1
            n+=1
            iter+=1

        end
    end
    model_F = compose_model(Astd, b, cstd, xlb, xub, flag_sense, solver, X, z)

    if isnan(getobjectivevalue(model_F))
        return model_F, convergence, getobjectivevalue(model_F), 0

    else
        return model_F, convergence, getobjectivevalue(model_F), 1
    end
end

function extract_data(model::JuMP.Model)

    A = full(JuMP.prepConstrMatrix(model))
    c = JuMP.prepAffObjective(model)
    m, n = size(A)

    xlb = copy(model.colLower)
    xub = copy(model.colUpper)

    rowlb, rowub = JuMP.prepConstrBounds(model)

    folgas = find(rowlb .!= rowub)
    s_f = size(folgas)[1]

    Astd = [A zeros(m,s_f)]

    for i = 1:size(folgas)[1]
        if rowlb[folgas[i]] == -Inf
            Astd[folgas[i],n+i] += 1
        else
            Astd[folgas[i],n+i] -= 1
        end
    end

    cstd = [c;zeros(s_f)]

    xlb = copy(model.colLower)
    xub = copy(model.colUpper)

    nxlb = zeros(n+s_f,1)
    nxub = zeros(n+s_f,1)

    nxlb[1:size(xlb)[1]] = xlb
    nxlb[size(xlb)[1]+1 : s_f+n] = 0

    nxub[1:size(xlb)[1]] = xub
    nxub[size(xub)[1]+1 : s_f+n] = +Inf

    b = zeros(m,1)
    for i in 1:m
        if rowlb[i] == -Inf
            b[i] = rowub[i]
        else
            b[i] = rowlb[i]
        end
    end

    flag_sense = 0

    if model.objSense == :Min
    cstd = -cstd
    flag_sense = 1
    end

    m, n = size(Astd)
    solver = model.solver

    return Astd, b, cstd, nxlb, nxub, solver, m, n, flag_sense
end

function solve_LP(Astd, b, cstd, xlb, xub, solver)
    m, n = size(Astd)
    model_LP = Model(solver=GurobiSolver(OutputFlag = 0))

    @variable(model_LP, xlb[i] <= x[i=1:n] <= xub[i])
    @constraints(model_LP, begin
    constrain[i=1:m], sum(Astd[i,j]*x[j] for j=1:n) == b[i]
    end)

    @objective(model_LP, :Max, sum(cstd[j]*x[j] for j=1:n))
    status = solve(model_LP);
    return model_LP.colVal, model_LP.objVal, status
end

function compose_model(Astd, b, cstd, xlb, xub, flag_sense, solver, X, z)

    m, n = size(Astd)

    model_F = Model(solver=solver)

    @variable(model_F, xlb[i] <= x[i=1:n] <= xub[i])
    @constraints(model_F, begin
    constrain[i=1:m], sum(Astd[i,j]*x[j] for j=1:n) == b[i]
    end)

    if flag_sense == 1
        cstd = -cstd
        @objective(model_F, :Min, sum(cstd[j]*x[j] for j=1:n))
        model_F.objVal = -z
    else
        @objective(model_F, :Max, sum(cstd[j]*x[j] for j=1:n))
        model_F.objVal = z
    end

    model_F.colVal = X

    return model_F
end


text:
function JDKRandomGenerator()
    return JDKRandomGenerator(())
end

function JDKRandomGenerator(arg0::Random)
    return JDKRandomGenerator((Random,), arg0)
end

function JDKRandomGenerator(arg0::jlong)
    return JDKRandomGenerator((jlong,), arg0)
end

function next_boolean(obj::JDKRandomGenerator)
    return jcall(obj, "nextBoolean", jboolean, ())
end

function next_bytes(obj::JDKRandomGenerator, arg0::Vector{jbyte})
    return jcall(obj, "nextBytes", void, (Vector{jbyte},), arg0)
end

function next_double(obj::JDKRandomGenerator)
    return jcall(obj, "nextDouble", jdouble, ())
end

function next_float(obj::JDKRandomGenerator)
    return jcall(obj, "nextFloat", jfloat, ())
end

function next_gaussian(obj::JDKRandomGenerator)
    return jcall(obj, "nextGaussian", jdouble, ())
end

function next_int(obj::JDKRandomGenerator)
    return jcall(obj, "nextInt", jint, ())
end

function next_int(obj::JDKRandomGenerator, arg0::jint)
    return jcall(obj, "nextInt", jint, (jint,), arg0)
end

function next_long(obj::JDKRandomGenerator)
    return jcall(obj, "nextLong", jlong, ())
end

function set_seed(obj::JDKRandomGenerator, arg0::Vector{jint})
    return jcall(obj, "setSeed", void, (Vector{jint},), arg0)
end

function set_seed(obj::JDKRandomGenerator, arg0::jint)
    return jcall(obj, "setSeed", void, (jint,), arg0)
end

function set_seed(obj::JDKRandomGenerator, arg0::jlong)
    return jcall(obj, "setSeed", void, (jlong,), arg0)
end


text:
using Test
using SparseArrays
using LinearAlgebra
using MtSpMV

# Test parallel SpMV
function spmv(m, n)

    # Create random test matrix and random test vector
    A = sprand(m, n, 5/m)
    transA = transpose(sparse(transpose(A)))
    x = rand(n)

    # Compute SpMV in serial for CSC matrix
    y_ser = zeros(m)
    for i = 1:n
        for j = A.colptr[i]:(A.colptr[i+1] - 1)
            y_ser[A.rowval[j]] += A.nzval[j] * x[i]
        end
    end

    # Compute SpMV in parallel with mul! and test
    y_par = zeros(m)
    mul!(y_par, transA, x)
    @test y_ser == y_par

    # Compute SpMV in parallel with * and test
    @test y_ser == transA * x
end

# Tests for quadratic and nonquadratic matrices
@testset "spmv" begin

    @testset "quadratic matrix" begin
        spmv(100, 100)
        spmv(1000, 1000)
    end

    @testset "nonquadratic matrix, nrows < ncols" begin
        spmv(80, 100)
        spmv(800, 1000)
    end

    @testset "nonquadratic matrix, nrows > ncols" begin
        spmv(100, 80)
        spmv(1000, 800)
    end
end


text:
"""
Plots a desnity of a distribution for a selected metric (e.g. effective sample size)
* `df`: dataframe of results
* `metric`: name of metric, such as :ess for effective sample size
* `group`: a tuple of grouping factors, e.g. (:sampler,:Nd)
* `save`: save=true saves each plot
* `figfmt`: figure format
* `dir`: directory of saved plot. Default pwd.
"""
function plotdensity(df::DataFrame, metric::Symbol, group=(:sampler,); save=false,
    figfmt="png", dir="", options...)
    plots = Plots.Plot[]
    layout = SetLayout(df, group)
    grouping = map(x->df[!,x], group)
    for c in names(df)
        !isin(metric,c) ? (continue) : nothing
        p=density(df[!,c], group=grouping, grid=false, xlabel=string(c),
            ylabel="Density", layout=layout, fill=(0,.5), width=1.5; options...)
        push!(plots,p)
        save ? savefig(p,string(dir, "density_", c, ".", figfmt)) : nothing
    end
    return plots
end

"""
Plots a summary of a metric computed according to a function, func. func defaults to the mean.
* `df`: dataframe of results
* `xvar`: variable assigned to x-axis
* `metric`: name of metric, such as :ess for effective sample size
* `group`: a tuple of grouping factors, e.g. (:sampler,:Nd)
* `save`: save=true saves each plot
* `figfmt`: figure format
* func: a function used to summarize results. Default mean
* `dir`: directory of saved plot. Default pwd.
"""
function plotsummary(df::DataFrame, xvar::Symbol, metric::Symbol, group=(:sampler,); save=false,
    figfmt="png", func=mean, dir="", options...)
    plots = Plots.Plot[]
    layout = SetLayout(df, group)
    for c in names(df)
        !isin(metric,c) ? (continue) : nothing
        summary,yvar = summarize(df, c, [xvar,group...], func)
        dfse,yse = summarize(df, c, [xvar,group...], se)
        grouping = map(x->summary[!,x], group)
        p=plot(summary[!,xvar], summary[!,yvar], group=grouping, grid=false, xlabel=string(xvar),
            ylabel=string(yvar), layout=layout, width=1.5, yerror=dfse[!,yse]; options...)
        push!(plots, p)
        display(p)
        save ? savefig(p, string(dir, "summary_", c, ".", figfmt)) : nothing
    end
    return plots
end

"""
Generates a scatter plot to visualize the relationship between two
benchmarking metrics.
* `df`: dataframe of results
* `xvar`: variable assigned to x-axis
* `metric`: name of metric, such as :ess for effective sample size
* `group`: a tuple of grouping factors, e.g. (:sampler,:Nd)
* `save`: save=true saves each plot
* `figfmt`: figure format
* func: a function used to summarize results. Default mean
* `dir`: directory of saved plot. Default pwd.
"""
function plotscatter(df::DataFrame, xvar::Symbol, metric::Symbol, group=(:sampler,); save=false,
    figfmt="png", func=mean,dir="", options...)
    plots = Plots.Plot[]
    layout = SetLayout(df, group)
    grouping = map(x->df[!,x], group)
    for c in names(df)
        !isin(metric,c) ? (continue) : nothing
        p=scatter(df[!,xvar], df[!,c], group=grouping, grid=false, xlabel=string(xvar),
            ylabel=string(c), layout=layout, width=1.5; options...)
        push!(plots, p)
        save ? savefig(p, string(dir, "scatter_", xvar, "_", c, ".", figfmt)) : nothing
    end
    return plots
end

"""
Generates a plot to compare the estimated parameters to true parameter values used to
generate simulated data.
* `df`: dataframe of results
* `parms`: namedtuple of parameter names and true values
* `group`: a tuple of grouping factors, e.g. (:sampler,:Nd)
* `save`: save=true saves each plot
* `figfmt`: figure format
* func: a function used to summarize results. Default mean
* `dir`: directory of saved plot. Default pwd.
"""
function plotrecovery(df::DataFrame, parms, group=(:sampler,); save=false,
    figfmt="png", dir="", options...)
    plots = Plots.Plot[]
    layout = SetLayout(df, group)
    for s in groupby(df, :sampler)
        sampler = s[!,:sampler][1]
        grouping = map(x-> s[!,x], group)
        for (parm,v) in pairs(parms)
            μ = s[!,Symbol(string(parm, "_mean"))]
            lb = μ .- s[!,Symbol(string(parm, "_hdp_lb"))]
            ub = s[!,Symbol(string(parm, "_hdp_ub"))] .- μ
            p=scatter(μ, group=grouping, grid=false, ylabel=string(parm), layout=layout,
                leg=false, yerror=(lb,ub), width=1, title=string(sampler); options...)
            cnt = 0
            [hline!(p, [v], subplot=cnt+=1) for i in 1:layout[1] for j in 1:layout[2]]
            push!(plots, p)
            save ? savefig(p, string(dir, sampler, "_recovery_", parm, ".", figfmt)) : nothing
        end
    end
    return plots
end

"""
Test whether column is a metric, e.g. mu_ess for ess.
"""
function isin(metric, col)
    occursin(string(metric), string(col))
end

"""
Creates a layout = (n,1) where n is the number of factors
in the last grouping variable.
"""
function SetLayout(df, group)
    isempty(group) ? (return (1,1)) : nothing
    length(group) == 1 ? (return (1,1)) : nothing
    col = group[end]
    n = length(unique(df[!,col]))
    return (n,1)
end

se(x) = std(x)/sqrt(length(x))

function summarize(df, metric, grouping, func)
    g = groupby(df, grouping)
    newDF= combine(g, metric=>func)
    yvar = names(newDF)[end]
    return newDF,yvar
end


text:
using Test
using Quiqbox

@testset "Differentiation.jl" begin

pb1 = ParamBox(1, :a)
@test inSymOf(pb1) == :a == outSymOf(pb1)
@test dataOf(pb1)[] == pb1[] == inValOf(pb1)
@test pb1[] == pb1() == 1 == outValOf(pb1)
@test pb1.map == Quiqbox.itself == mapOf(pb1)
@test pb1.canDiff[] == true == isDiffParam(pb1)
toggleDiff!(pb1)
@test pb1.canDiff[] == false
disableDiff!(pb1)
@test pb1.canDiff[] == false
enableDiff!(pb1)
@test pb1.canDiff[] == true
pb2 = inVarCopy(pb1)
@test dataOf(pb1) === dataOf(pb2) === pb1.data
pb2_2 = outValCopy(pb1)
@test dataOf(pb1) == dataOf(pb2_2)
@test dataOf(pb1) !== dataOf(pb2_2)

pb3 = ParamBox(1, :b, abs)
@test inSymOf(pb3) == :x_b
@test outSymOf(pb3) == :b
@test mapOf(pb3) == abs

pb4 = ParamBox(1.2, :c, x->x^2, :x)
@test inSymOf(pb4) == :x
@test outSymOf(pb4) == :c
@test startswith(nameof(pb4.map) |> string, "f_c")

pb5 = outValCopy(pb4)
@test pb5() == pb5[] == pb4()
@test pb5.map == Quiqbox.itself
end


text:
# This is a simplified version of the primal-dual hybrid gradient
# and extragradient implementation, of the package FirstOrderLp
# (https://github.com/google-research/FirstOrderLp.jl)
# aimed to be used for experiments with restarts.

import Printf
import LinearAlgebra
import Random
import DataFrames
const norm = LinearAlgebra.norm
const dot = LinearAlgebra.dot
const randn = Random.randn
include("restart_scheme.jl")

"""
Information for recording iterates on two dimensional subspace.
"""
struct TwoDimensionalSubspace
  """
  Subspace basis vector 1
  """
  basis_vector1::Vector{Float64}
  """
  Subspace basis vector 2, should be orthognal to vector 1
  """
  basis_vector2::Vector{Float64}
end

@enum PrimalDualMethod PDHG EXTRAGRADIENT

function RecoverPrimalDualMethodFromString(method_name::String)
  if method_name == "PDHG"
    return PDHG
  elseif method_name == "extragradient"
    return EXTRAGRADIENT
  else
    error("unsupported method choosen")
  end
end

"""
A PrimalDualOptimizerParameters struct specifies the parameters for solving the saddle
point formulation of an problem using primal-dual hybrid gradient.

It solves a problem of the form (see quadratic_programmming.jl in FirstOrderLp)
minimize objective_vector' * x

s.t. constraint_matrix[1:num_equalities, :] * x =
     right_hand_side[1:num_equalities]

     constraint_matrix[(num_equalities + 1):end, :] * x >=
     right_hand_side[(num_equalities + 1):end, :]

     variable_lower_bound <= x <= variable_upper_bound

We use notation from Chambolle and Pock, "On the ergodic convergence rates of a
first-order primal-dual algorithm"
(http://www.optimization-online.org/DB_FILE/2014/09/4532.pdf).
That paper doesn't explicitly use the terminology "primal-dual hybrid gradient"
but their Theorem 1 is analyzing PDHG. In this file "Theorem 1" without further
reference refers to that paper.

Our problem is equivalent to the saddle point problem:
    min_x max_y L(x, y)
where
    L(x, y) = y' K x + g(x) - h*(y)
    K = -constraint_matrix
    g(x) = objective_vector' x if variable_lower_bound <= x <=
                                                            variable_upper_bound
                               otherwise infinity
    h*(y) = -right_hand_side' y if y[(num_equalities + 1):end] >= 0
                                otherwise infinity

Note that the places where g(x) and h*(y) are infinite effectively limits the
domain of the min and max. Therefore there's no infinity in the code.

We parametrize the primal and dual step sizes (tau and sigma in Chambolle and
Pock) as:
    primal_step_size = step_size / primal_weight
    dual_step_size = step_size * primal_weight.
The algoritm converges if
    primal_stepsize * dual_stepsize * norm(contraint_matrix)^2 < 1.
"""
mutable struct PrimalDualOptimizerParameters
  """
  method = PDHG or extragradient
  """
  method::PrimalDualMethod

  """
  Constant step size used in the algorithm. If nothing is specified, the solver
  computes a provably correct step size.
  """
  step_size::Union{Float64,Nothing}
  """
  Weight relating primal and dual step sizes.
  """
  primal_weight::Float64
  """
  Records iteration stats to csv every record_every iterations.
  During these iteration restarts may also be performed.
  """
  record_every::Int64
  """
  If verbose is true then prints iteration stats every print_every time that
  iteration information is recorded.
  """
  print_every::Int64
  """
  If true a line of debugging info is printed every printing_every
  iterations.
  """
  verbosity::Bool
  """
  Number of loop iterations to run. Must be postive.
  """
  iteration_limit::Int64
  """
  Which restart scheme should we use
  """
  restart_scheme::RestartScheme
  """
  A subspace where projections of iterates are recorded.
  If the value is nothing then the values will be missing from the iteration
  stats.
  """
  two_dimensional_subspace::Union{TwoDimensionalSubspace,Nothing}

  kkt_tolerance::Float64
end

function create_stats_data_frame()
  return DataFrames.DataFrame(
    # Record the current iteration.
    iteration = Int64[],
    # Primal objectives of the iterates; the ith entry corresponds to the primal
    # objective of the ith iterate.
    primal_objectives = Float64[],
    # Dual objectives of the iterates; the ith entry corresponds to the dual
    # objective of the ith iterate.
    dual_objectives = Float64[],
    # Primal norms of the iterates; the ith entry corresponds to the primal
    # norm of the ith iterate.
    primal_solution_norms = Float64[],
    # Dual norms of the iterates; the ith entry corresponds to the dual
    # norm of the ith iterate.
    dual_solution_norms = Float64[],
    # Primal delta norms of the iterates; the ith entry corresponds to the
    # primal delta norm of the ith iterate.
    primal_delta_norms = Float64[],
    # Dual delta norms of the iterates; the ith entry corresponds to the dual
    # delta norm of the ith iterate.
    dual_delta_norms = Float64[],
    # First coordinate of a subspace that the current iterates are projected
    # onto.
    current_subspace_coordinate1 = Float64[],
    # Second coordinate of a subspace that the current iterates are projected
    # onto.
    current_subspace_coordinate2 = Float64[],
    # First coordinate of a subspace that the average iterates are projected
    # onto.
    average_subspace_coordinate1 = Float64[],
    # Second coordinate of a subspace that the average iterates are projected
    # onto.
    average_subspace_coordinate2 = Float64[],
    # Did a restart occur at this iteration
    restart_occurred = Bool[],
    # Number of variables that moved in or out of the active set at this
    # iteration
    number_of_active_set_changes = Int64[],
    kkt_error_current_iterate = Float64[],
    kkt_error_average_iterate = Float64[],
    average_normalized_gap = Float64[],
    current_normalized_gap = Float64[],
  )
end

"""
Output of the solver.
"""
struct PrimalDualOutput
  primal_solution::Vector{Float64}
  dual_solution::Vector{Float64}
  iteration_stats::DataFrames.DataFrame
end

function compute_convergence_information(
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  primal_solution::Vector{Float64},
  dual_solution::Vector{Float64},
)
  dual_stats =
    FirstOrderLp.compute_dual_stats(problem, primal_solution, dual_solution)
  primal_objective = FirstOrderLp.primal_obj(problem, primal_solution)
  dual_objective = dual_stats.dual_objective
  primal_residual =
    FirstOrderLp.compute_primal_residual(problem, primal_solution)
  kkt_error = sqrt(
    sum(primal_residual .^ 2) +
    sum(dual_stats.dual_residual .^ 2) +
    (primal_objective - dual_objective)^2,
  )
  return (
    primal_objective = primal_objective,
    dual_objective = dual_objective,
    kkt_error = kkt_error,
  )
end

"""
Computes statistics for the current iteration. The arguments primal_delta and
dual_delta correspond to the difference between the last two primal and dual
iterates, respectively.
"""
function compute_stats(
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  iteration::Int64,
  two_dimensional_subspace::Union{TwoDimensionalSubspace,Nothing},
  primal_solution::Vector{Float64},
  dual_solution::Vector{Float64},
  solution_weighted_avg::FirstOrderLp.SolutionWeightedAverage,
  primal_delta::Vector{Float64},
  dual_delta::Vector{Float64},
  restart_occurred::Bool,
  number_of_active_set_changes::Int64,
  normalized_gaps::NormalizedDualityGaps,
)
  if restart_occurred
    primal_avg, dual_avg = primal_solution, dual_solution
  else
    primal_avg, dual_avg = FirstOrderLp.compute_average(solution_weighted_avg)
  end
  current_iterate_convergence_info =
    compute_convergence_information(problem, primal_solution, dual_solution)
  average_iterate_convergence_info =
    compute_convergence_information(problem, primal_avg, dual_avg)

  kkt_error_average_iterate = Inf

  if nothing != two_dimensional_subspace
    z = [primal_solution; dual_solution]
    current_subspace_coordinate1 =
      dot(z, two_dimensional_subspace.basis_vector1)
    current_subspace_coordinate2 =
      dot(z, two_dimensional_subspace.basis_vector2)
    z_avg = [primal_avg; dual_avg]
    average_subspace_coordinate1 =
      dot(z_avg, two_dimensional_subspace.basis_vector1)
    average_subspace_coordinate2 =
      dot(z_avg, two_dimensional_subspace.basis_vector2)
  else
    current_subspace_coordinate1 = NaN
    current_subspace_coordinate2 = NaN
    average_subspace_coordinate1 = NaN
    average_subspace_coordinate2 = NaN
  end

  return DataFrames.DataFrame(
    iteration = iteration,
    primal_objectives = current_iterate_convergence_info.primal_objective,
    dual_objectives = current_iterate_convergence_info.dual_objective,
    primal_solution_norms = norm(primal_solution),
    dual_solution_norms = norm(dual_solution),
    primal_delta_norms = norm(primal_delta),
    dual_delta_norms = norm(dual_delta),
    current_subspace_coordinate1 = current_subspace_coordinate1,
    current_subspace_coordinate2 = current_subspace_coordinate2,
    average_subspace_coordinate1 = average_subspace_coordinate1,
    average_subspace_coordinate2 = average_subspace_coordinate2,
    restart_occurred = restart_occurred,
    number_of_active_set_changes = number_of_active_set_changes,
    kkt_error_current_iterate = current_iterate_convergence_info.kkt_error,
    kkt_error_average_iterate = average_iterate_convergence_info.kkt_error,
    average_normalized_gap = normalized_gaps.average_gap,
    current_normalized_gap = normalized_gaps.current_gap,
  )
end

"""
Logging while the algorithm is running.
"""
function log_iteration(
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  iteration_stats::DataFrames.DataFrame,
)

  Printf.@printf(
    "   %5d objectives=(%9g, %9g) norms=(%9g, %9g) res_norm=(%9g, %9g)\n",
    iteration_stats[:iteration][end],
    iteration_stats[:primal_objectives][end],
    iteration_stats[:dual_objectives][end],
    iteration_stats[:primal_solution_norms][end],
    iteration_stats[:dual_solution_norms][end],
    iteration_stats[:primal_delta_norms][end],
    iteration_stats[:dual_delta_norms][end],
  )
end

function take_pdhg_step(
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  current_primal_solution::Vector{Float64},
  current_dual_solution::Vector{Float64},
  primal_weight::Float64,
  step_size::Float64,
)
  # The next lines compute the primal portion of the PDHG algorithm:
  # argmin_x [g(x) + current_dual_solution' K x
  #          + (0.5 * primal_weight / step_size)
  #             * norm(x - current_primal_solution)^2]
  # See Sections 2-3 of Chambolle and Pock and the comment above
  # PrimalDualOptimizerParameters.
  # This minimization is easy to do in closed form since it can be separated
  # into independent problems for each of the primal variables. The
  # projection onto the primal feasibility set comes from the closed form
  # for the above minimization and the cases where g(x) is infinite - there
  # isn't officially any projection step in the algorithm.
  primal_gradient = FirstOrderLp.compute_primal_gradient(
    problem,
    current_primal_solution,
    current_dual_solution,
  )
  next_primal =
    current_primal_solution .- primal_gradient * (step_size / primal_weight)
  FirstOrderLp.project_primal!(next_primal, problem)

  # The next two lines compute the dual portion:
  # argmin_y [H*(y) - y' K (2.0*next_primal - current_primal_solution)
  #           + (0.5 / (primal_weight * step_size))
  #.             * norm(y-current_dual_solution)^2]
  dual_gradient = FirstOrderLp.compute_dual_gradient(
    problem,
    2.0 * next_primal - current_primal_solution,
  )
  next_dual =
    current_dual_solution .+ dual_gradient * (step_size * primal_weight)
  FirstOrderLp.project_dual!(next_dual, problem)

  return next_primal, next_dual
end

function proximal_step(
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  current_primal_solution::Vector{Float64},
  current_dual_solution::Vector{Float64},
  current_primal_gradient::Vector{Float64},
  current_dual_gradient::Vector{Float64},
  primal_weight::Float64,
  step_size::Float64,
)
  prox_primal =
    current_primal_solution -
    (step_size / primal_weight) * current_primal_gradient
  FirstOrderLp.project_primal!(prox_primal, problem)
  prox_dual =
    current_dual_solution + (step_size * primal_weight) * current_dual_gradient
  FirstOrderLp.project_dual!(prox_dual, problem)
  return prox_primal, prox_dual
end

function take_extragradient_step(
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  current_primal_solution::Vector{Float64},
  current_dual_solution::Vector{Float64},
  primal_weight::Float64,
  step_size::Float64,
  solution_weighted_avg::FirstOrderLp.SolutionWeightedAverage,
)
  current_primal_gradient = FirstOrderLp.compute_primal_gradient(
    problem,
    current_primal_solution,
    current_dual_solution,
  )
  current_dual_gradient =
    FirstOrderLp.compute_dual_gradient(problem, current_primal_solution)
  test_primal_solution, test_dual_solution = proximal_step(
    problem,
    current_primal_solution,
    current_dual_solution,
    current_primal_gradient,
    current_dual_gradient,
    primal_weight,
    step_size,
  )
  test_primal_gradient = FirstOrderLp.compute_primal_gradient(
    problem,
    test_primal_solution,
    test_dual_solution,
  )
  test_dual_gradient =
    FirstOrderLp.compute_dual_gradient(problem, test_primal_solution)
  FirstOrderLp.add_to_solution_weighted_average(
    solution_weighted_avg,
    test_primal_solution,
    test_dual_solution,
    1.0,
  )
  return proximal_step(
    problem,
    current_primal_solution,
    current_dual_solution,
    test_primal_gradient,
    test_dual_gradient,
    primal_weight,
    step_size,
  )
end

function compute_number_of_active_set_changes(
  problem,
  current_primal_solution::Vector{Float64},
  current_dual_solution::Vector{Float64},
  primal_delta::Vector{Float64},
  dual_delta::Vector{Float64},
)
  number_of_active_set_changes = 0
  for idx in 1:length(current_primal_solution)
    if primal_delta[idx] != 0.0 && (
      current_primal_solution[idx] == problem.variable_lower_bound[idx] ||
      current_primal_solution[idx] == problem.variable_upper_bound[idx]
    )
      number_of_active_set_changes += 1
    end
  end
  for idx in (problem.num_equalities+1):length(current_dual_solution)
    if dual_delta[idx] != 0.0 && current_dual_solution[idx] == 0.0
      number_of_active_set_changes += 1
    end
  end
  return number_of_active_set_changes
end

"""
`optimize(params::PrimalDualOptimizerParameters,
          problem::QuadraticProgrammingProblem)`

Solves a linear program using primal-dual hybrid gradient or
extragradient. If the step_size
specified in params is negative, picks a step size that ensures
step_size^2 * norm(constraint_matrix)^2 < 1,
a condition that guarantees provable convergence.

# Arguments
- `params::PrimalDualOptimizerParameters`: parameters.
- `original_problem::QuadraticProgrammingProblem`: the QP to solve.

# Returns
A PrimalDualOutput struct containing the solution found.
"""
function optimize(
  params::PrimalDualOptimizerParameters,
  problem::FirstOrderLp.QuadraticProgrammingProblem,
  initial_primal_solution::Vector{Float64},
  initial_dual_solution::Vector{Float64},
)
  primal_size = length(problem.variable_lower_bound)
  dual_size = length(problem.right_hand_side)

  primal_weight = params.primal_weight

  if isnothing(params.step_size)
    desired_relative_error = 0.2
    maximum_singular_value, number_of_power_iterations =
      FirstOrderLp.estimate_maximum_singular_value(
        problem.constraint_matrix,
        probability_of_failure = 0.001,
        desired_relative_error = desired_relative_error,
      )
    step_size = (1 - desired_relative_error) / maximum_singular_value
  else
    step_size = params.step_size
  end
  if params.verbosity
    Printf.@printf(
      "Step size: %9g, Primal weight: %9g, method: %s\n",
      step_size,
      primal_weight,
      params.method
    )
  end

  restart_scheme = params.restart_scheme
  restart_info = initialize_restart_info(primal_size, dual_size)
  solution_weighted_avg =
    FirstOrderLp.initialize_solution_weighted_average(primal_size, dual_size)
  # Difference between current_primal_solution and last primal iterate, i.e.
  # x_k - x_{k-1}.
  primal_delta = zeros(primal_size)
  # Difference between current_dual_solution and last dual iterate, i.e. y_k -
  # y_{k-1}
  dual_delta = zeros(dual_size)
  iteration_limit = params.iteration_limit
  stats = create_stats_data_frame()

  current_primal_solution = initial_primal_solution
  current_dual_solution = initial_dual_solution

  restart_occurred = false
  number_of_active_set_changes = 0
  iteration = 0
  while true
    # store stats and log
    terminate = false
    if iteration >= iteration_limit
      if params.verbosity
        println("Iteration limit reached")
      end
      terminate = true
    end
    do_fixed_frequency_restart =
      isa(restart_scheme, FixedFrequencyRestarts) &&
      solution_weighted_avg.sum_primal_solutions_count >=
      restart_scheme.restart_length
    store_stats =
      mod(iteration, params.record_every) == 0 ||
      terminate ||
      do_fixed_frequency_restart
    print_stats =
      params.verbosity && (
        mod(iteration, params.record_every * params.print_every) == 0 ||
        terminate
      )
    if store_stats
      if iteration > 0
        restart_occurred, normalized_gap_info = run_restart_scheme(
          restart_scheme,
          problem,
          current_primal_solution,
          current_dual_solution,
          primal_weight,
          solution_weighted_avg,
          restart_info,
          params.verbosity,
        )
      else
        restart_occurred = false
        normalized_gap_info = NormalizedDualityGaps(NaN, NaN, NaN, NaN)
      end
      this_iteration_stats = compute_stats(
        problem,
        iteration,
        params.two_dimensional_subspace,
        current_primal_solution,
        current_dual_solution,
        solution_weighted_avg,
        primal_delta,
        dual_delta,
        restart_occurred,
        number_of_active_set_changes,
        normalized_gap_info,
      )
      kkt_error = min(
        this_iteration_stats.kkt_error_current_iterate[1],
        this_iteration_stats.kkt_error_average_iterate[1],
      )
      if kkt_error <= params.kkt_tolerance
        if params.verbosity
          println("Found optimal solution")
        end
        terminate = true
      end
    end
    if store_stats
      append!(stats, this_iteration_stats)
    end
    if print_stats
      log_iteration(problem, this_iteration_stats)
    end

    if terminate
      return PrimalDualOutput(
        current_primal_solution,
        current_dual_solution,
        stats,
      )
    end
    iteration += 1

    if params.method == PDHG
      next_primal, next_dual = take_pdhg_step(
        problem,
        current_primal_solution,
        current_dual_solution,
        primal_weight,
        step_size,
      )
      FirstOrderLp.add_to_solution_weighted_average(
        solution_weighted_avg,
        next_primal,
        next_dual,
        1.0,
      )
    elseif params.method == EXTRAGRADIENT
      next_primal, next_dual = take_extragradient_step(
        problem,
        current_primal_solution,
        current_dual_solution,
        primal_weight,
        step_size,
        solution_weighted_avg,
      )
    else
      error("unknown method")
    end
    # Update deltas and iterates
    primal_delta = next_primal .- current_primal_solution
    dual_delta = next_dual .- current_dual_solution

    number_of_active_set_changes = compute_number_of_active_set_changes(
      problem,
      current_primal_solution,
      current_dual_solution,
      primal_delta,
      dual_delta,
    )

    # update iterates
    current_primal_solution = next_primal
    current_dual_solution = next_dual
  end
end

function optimize(
  params::PrimalDualOptimizerParameters,
  problem::FirstOrderLp.QuadraticProgrammingProblem,
)
  current_primal_solution = zeros(length(problem.variable_lower_bound))
  current_dual_solution = zeros(length(problem.right_hand_side))

  return optimize(
    params,
    problem,
    current_primal_solution,
    current_dual_solution,
  )
end

max_stars_repo_name: pedersor/google-research

max_issues_repo_name: pedersor/google-research

max_forks_repo_name: pedersor/google-research


text:
fetchvcat(J::judiVector) = J
fetchvcat(x) = vcat(x...)

function JUDI.time_modeling(model::Model, srcGeometry, srcData, recGeometry, recData, dm, srcnum::UnitRange{Int64}, op::Char, mode::Int64, options)

    # Broadcast common parameters
    _model = @bcast model
    _dm = isnothing(dm) ? dm : @bcast dm
    opts = make_opt([model, dm, op, srcGeometry])
    iter = make_parts(srcnum)
    # Run on azure
    results = @batchexec pmap(sx -> time_modeling_azure(_model, subsample(srcGeometry, sx), subsample(srcData, sx),
                                                        subsample(recGeometry, sx), subsample(recData, sx), _dm,
                                                        op, mode, subsample(options, sx)), iter) opts
    # Gather results
    if op=='F' || (op=='J' && mode==1)
        temp = fetch(results)
        argout1 = fetchvcat(temp)
    elseif op=='J' && mode==-1
        argout1 = fetchreduce(results; op=+, remote=true)
    else
        error("operation no defined")
    end
    return argout1
end

# FWI
function JUDI.fwi_objective(model::Model, source::judiVector, dObs::judiVector; options=JUDI.Options())
    # fwi_objective function for multiple sources. The function distributes the sources and the input data amongst the available workers.
    # Broadcast common parameters
    _model = @bcast model
    opts = make_opt([model, dObs, source])
    iter = make_parts(1:dObs.nsrc)
    results = @batchexec pmap(j -> fwi_objective_azure(_model, source[j], dObs[j], subsample(options, j)), iter) opts

    # Collect and reduce gradients
    obj, gradient = fetchreduce(results; op=+, remote=true)
    # first value corresponds to function value, the rest to the gradient
    return obj, gradient
end

# lsrtm

function JUDI.lsrtm_objective(model::Model, source::judiVector, dObs::judiVector, dm; options=JUDI.Options(), nlind=false)
    # lsrtm_objective function for multiple sources. The function distributes the sources and the input data amongst the available workers.

    # Broadcast common parameters
    _model = @bcast model
    _dm = isnothing(dm) ? dm : @bcast dm
    opts = make_opt([model, dObs, source, dm])
    iter = make_parts(1:dObs.nsrc)
    results = @batchexec pmap(j -> lsrtm_objective_azure(_model, source[j], dObs[j], _dm, subsample(options, j); nlind=nlind), iter) opts
    # Collect and reduce gradients
    obj, gradient = fetchreduce(results; op=+, remote=true)
    # first value corresponds to function value, the rest to the gradient
    return obj, gradient
end

# Make options for unique batch ids
make_opt(a::Array) = AzureClusterlessHPC.Options(;task_name="task_$(objectid(a))")

# Batchexec function
slicec(i::Integer, j::Integer) = i==j ? i : (i:j)
chunk(arr, n) = [arr[slicec(i, min(i + n - 1, length(arr)))] for i in 1:n:length(arr)]
make_parts(iter) = _njpi > 1 ? chunk(iter, _njpi) : iter


text:
"""
    initialize_hbv_model(config::Config)

Initial part of the SBM model concept. Reads the input settings and data as defined in the
Config object. Will return a Model that is ready to run.
"""
function initialize_hbv_model(config::Config)
    # unpack the paths to the NetCDF files
    tomldir = dirname(config)
    static_path = joinpath(tomldir, config.input.path_static)
    dynamic_path = joinpath(tomldir, config.input.path_forcing)

    reader = prepare_reader(dynamic_path, static_path, config)
    clock = Clock(config, reader)
    Δt = clock.Δt

    reinit = get(config.model, "reinit", true)::Bool
    do_reservoirs = get(config.model, "reservoirs", false)::Bool
    do_lakes = get(config.model, "lakes", false)::Bool
    do_pits = get(config.model, "pits", false)::Bool
    set_kquickflow = get(config.model, "set_kquickflow", false)::Bool

    kw_river_tstep = get(config.model, "kw_river_tstep", 0)
    kw_land_tstep = get(config.model, "kw_land_tstep", 0)

    nc = NCDataset(static_path)
    dims = dimnames(nc[param(config, "input.subcatchment")])

    subcatch_2d = ncread(nc, param(config, "input.subcatchment"); allow_missing = true)
    # indices based on catchment
    inds, rev_inds = active_indices(subcatch_2d, missing)
    n = length(inds)

    cfmax =
        ncread(
            nc,
            param(config, "input.vertical.cfmax", nothing);
            sel = inds,
            defaults = 3.75653,
            type = Float,
        ) .* (Δt / basetimestep)
    tt = ncread(
        nc,
        param(config, "input.vertical.tt", nothing);
        sel = inds,
        defaults = -1.41934,
        type = Float,
    )
    tti = ncread(
        nc,
        param(config, "input.vertical.tti", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float,
    )
    ttm = ncread(
        nc,
        param(config, "input.vertical.ttm", nothing);
        sel = inds,
        defaults = -1.41934,
        type = Float,
    )
    whc = ncread(
        nc,
        param(config, "input.vertical.whc", nothing);
        sel = inds,
        defaults = 0.1,
        type = Float,
    )
    # glacier parameters
    g_tt = ncread(
        nc,
        param(config, "input.vertical.g_tt", nothing);
        sel = inds,
        defaults = 0.0,
        type = Float,
        fill = 0.0,
    )
    g_cfmax =
        ncread(
            nc,
            param(config, "input.vertical.g_cfmax", nothing);
            sel = inds,
            defaults = 3.0,
            type = Float,
            fill = 0.0,
        ) .* (Δt / basetimestep)
    g_sifrac = ncread(
        nc,
        param(config, "input.vertical.g_sifrac", nothing);
        sel = inds,
        defaults = 0.001,
        type = Float,
        fill = 0.0,
    )
    glacierfrac = ncread(
        nc,
        param(config, "input.vertical.glacierfrac", nothing);
        sel = inds,
        defaults = 0.0,
        type = Float,
        fill = 0.0,
    )
    glacierstore = ncread(
        nc,
        param(config, "input.vertical.glacierstore", nothing);
        sel = inds,
        defaults = 5500.0,
        type = Float,
        fill = 0.0,
    )
    fc = ncread(
        nc,
        param(config, "input.vertical.fc", nothing);
        sel = inds,
        defaults = 260.0,
        type = Float,
    )
    betaseepage = ncread(
        nc,
        param(config, "input.vertical.betaseepage", nothing);
        sel = inds,
        defaults = 1.8,
        type = Float,
    )
    lp = ncread(
        nc,
        param(config, "input.vertical.lp", nothing);
        sel = inds,
        defaults = 0.53,
        type = Float,
    )
    k4 =
        ncread(
            nc,
            param(config, "input.vertical.k4", nothing);
            sel = inds,
            defaults = 0.02307,
            type = Float,
        ) .* (Δt / basetimestep)
    kquickflow =
        ncread(
            nc,
            param(config, "input.vertical.kquickflow", nothing);
            sel = inds,
            defaults = 0.09880,
            type = Float,
        ) .* (Δt / basetimestep)
    suz = ncread(
        nc,
        param(config, "input.vertical.suz", nothing);
        sel = inds,
        defaults = 100.0,
        type = Float,
    )
    k0 =
        ncread(
            nc,
            param(config, "input.vertical.k0", nothing);
            sel = inds,
            defaults = 0.30,
            type = Float,
        ) .* (Δt / basetimestep)
    khq =
        ncread(
            nc,
            param(config, "input.vertical.khq", nothing);
            sel = inds,
            defaults = 0.09880,
            type = Float,
        ) .* (Δt / basetimestep)
    hq =
        ncread(
            nc,
            param(config, "input.vertical.hq", nothing);
            sel = inds,
            defaults = 3.27,
            type = Float,
        ) .* (Δt / basetimestep)
    alphanl = ncread(
        nc,
        param(config, "input.vertical.alphanl", nothing);
        sel = inds,
        defaults = 1.1,
        type = Float,
    )
    perc =
        ncread(
            nc,
            param(config, "input.vertical.perc", nothing);
            sel = inds,
            defaults = 0.4,
            type = Float,
        ) .* (Δt / basetimestep)
    cfr = ncread(
        nc,
        param(config, "input.vertical.cfr", nothing);
        sel = inds,
        defaults = 0.05,
        type = Float,
    )
    pcorr = ncread(
        nc,
        param(config, "input.vertical.pcorr", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float,
    )
    rfcf = ncread(
        nc,
        param(config, "input.vertical.rfcf", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float,
    )
    sfcf = ncread(
        nc,
        param(config, "input.vertical.sfcf", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float,
    )
    cflux =
        ncread(
            nc,
            param(config, "input.vertical.cflux", nothing);
            sel = inds,
            defaults = 2.0,
            type = Float,
        ) .* (Δt / basetimestep)
    icf = ncread(
        nc,
        param(config, "input.vertical.icf", nothing);
        sel = inds,
        defaults = 2.0,
        type = Float,
    )
    cevpf = ncread(
        nc,
        param(config, "input.vertical.cevpf", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float,
    )
    epf = ncread(
        nc,
        param(config, "input.vertical.epf", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float,
    )
    ecorr = ncread(
        nc,
        param(config, "input.vertical.ecorr", nothing);
        sel = inds,
        defaults = 1.0,
        type = Float,
    )

    # read x, y coordinates and calculate cell length [m]
    y_nc = read_y_axis(nc)
    x_nc = read_x_axis(nc)
    y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
    cellength = abs(mean(diff(x_nc)))

    sizeinmetres = get(config.model, "sizeinmetres", false)::Bool
    xl, yl = cell_lengths(y, cellength, sizeinmetres)

    threshold = fc .* lp

    altitude =
        ncread(nc, param(config, "input.vertical.altitude"); sel = inds, type = Float)

    hbv = HBV{Float}(
        Δt = Float(tosecond(Δt)),
        n = n,
        fc = fc,
        betaseepage = betaseepage,
        lp = lp,
        threshold = threshold,
        k4 = k4,
        kquickflow = set_kquickflow ? kquickflow :
                     pow.(khq, 1.0 .+ alphanl) .* pow.(hq, -alphanl),
        suz = suz,
        k0 = k0,
        khq = khq,
        hq = hq,
        alphanl = alphanl,
        perc = perc,
        cfr = cfr,
        pcorr = pcorr,
        rfcf = rfcf,
        sfcf = sfcf,
        cflux = cflux,
        icf = icf,
        cevpf = cevpf,
        epf = epf,
        ecorr = ecorr,
        tti = tti,
        tt = tt,
        ttm = ttm,
        cfmax = cfmax,
        whc = whc,
        # glacier parameters
        g_tt = g_tt,
        g_sifrac = g_sifrac,
        g_cfmax = g_cfmax,
        glacierstore = glacierstore,
        glacierfrac = glacierfrac,
        # default (cold) states:
        interceptionstorage = zeros(Float, n),
        snow = zeros(Float, n),
        snowwater = zeros(Float, n),
        soilmoisture = copy(fc),
        upperzonestorage = 0.2 .* fc,
        lowerzonestorage = 1.0 ./ (3.0 .* k4),
        # variables:
        precipitation = fill(mv, n),
        temperature = fill(mv, n),
        potential_evaporation = fill(mv, n),
        potsoilevap = fill(mv, n),
        soilevap = fill(mv, n),
        intevap = fill(mv, n),
        actevap = fill(mv, n),
        rainfallplusmelt = fill(mv, n),
        directrunoff = fill(mv, n),
        hbv_seepage = fill(mv, n),
        in_upperzone = fill(mv, n),
        quickflow = fill(mv, n),
        real_quickflow = fill(mv, n),
        percolation = fill(mv, n),
        capflux = fill(mv, n),
        baseflow = fill(mv, n),
        runoff = fill(mv, n),
    )

    modelsize_2d = size(subcatch_2d)
    river_2d = ncread(nc, param(config, "input.river_location"); type = Bool, fill = false)
    river = river_2d[inds]
    riverwidth_2d =
        ncread(nc, param(config, "input.lateral.river.width"); type = Float, fill = 0)
    riverwidth = riverwidth_2d[inds]
    riverlength_2d =
        ncread(nc, param(config, "input.lateral.river.length"); type = Float, fill = 0)
    riverlength = riverlength_2d[inds]

    inds_riv, rev_inds_riv = active_indices(river_2d, 0)
    nriv = length(inds_riv)

    # reservoirs
    pits = zeros(Bool, modelsize_2d)
    if do_reservoirs
        reservoirs, resindex, reservoir, pits =
            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits)
    else
        reservoir = ()
    end

    # lakes
    if do_lakes
        lakes, lakeindex, lake, pits =
            initialize_natural_lake(config, static_path, nc, inds_riv, nriv, pits)
    else
        lake = ()
    end

    ldd_2d = ncread(nc, param(config, "input.ldd"); allow_missing = true)
    ldd = ldd_2d[inds]
    if do_pits
        pits_2d = ncread(nc, param(config, "input.pits"); type = Bool, fill = false)
        ldd = set_pit_ldd(pits_2d, ldd, inds)
    end

    βₗ = ncread(nc, param(config, "input.lateral.land.slope"); sel = inds, type = Float)
    clamp!(βₗ, 0.00001, Inf)
    dl = fill(mv, n)
    dw = fill(mv, n)

    for i = 1:n
        dl[i] = detdrainlength(ldd[i], xl[i], yl[i])
        dw[i] = detdrainwidth(ldd[i], xl[i], yl[i])
    end

    n_land = ncread(
        nc,
        param(config, "input.lateral.land.n", nothing);
        sel = inds,
        defaults = 0.072,
        type = Float,
    )

    alpha_pow = Float((2.0 / 3.0) * 0.6)
    β = Float(0.6)
    olf = SurfaceFlow(
        β = β,
        sl = βₗ,
        n = n_land,
        dl = dl,
        q = zeros(Float, n),
        qin = zeros(Float, n),
        q_av = zeros(Float, n),
        qlat = zeros(Float, n),
        inwater = zeros(Float, n),
        volume = zeros(Float, n),
        h = zeros(Float, n),
        h_av = zeros(Float, n),
        Δt = Float(tosecond(Δt)),
        its = kw_land_tstep > 0 ? Int(cld(tosecond(Δt), kw_land_tstep)) : kw_land_tstep,
        width = dw,
        wb_pit = pits[inds],
        alpha_pow = alpha_pow,
        alpha_term = fill(mv, n),
        α = fill(mv, n),
        eps = Float(1e-3),
        cel = zeros(Float, n),
        to_river = zeros(Float, n),
        rivercells = fill(false, n),
        reservoir_index = fill(0, n),
        lake_index = fill(0, n),
        reservoir = nothing,
        lake = nothing,
    )

    graph = flowgraph(ldd, inds, pcr_dir)

    riverslope =
        ncread(nc, param(config, "input.lateral.river.slope"); sel = inds_riv, type = Float)
    clamp!(riverslope, 0.00001, Inf)
    riverlength = riverlength_2d[inds_riv]
    riverwidth = riverwidth_2d[inds_riv]
    n_river = ncread(
        nc,
        param(config, "input.lateral.river.n", nothing);
        sel = inds_riv,
        defaults = 0.036,
        type = Float,
    )
    ldd_riv = ldd_2d[inds_riv]
    if do_pits
        ldd_riv = set_pit_ldd(pits_2d, ldd_riv, inds_riv)
    end
    graph_riv = flowgraph(ldd_riv, inds_riv, pcr_dir)

    # the indices of the river cells in the land(+river) cell vector
    index_river = filter(i -> !isequal(river[i], 0), 1:n)
    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, βₗ, n)

    rf = SurfaceFlow(
        β = β,
        sl = riverslope,
        n = n_river,
        dl = riverlength,
        q = zeros(Float, nriv),
        qin = zeros(Float, nriv),
        q_av = zeros(Float, nriv),
        qlat = zeros(Float, nriv),
        inwater = zeros(Float, nriv),
        volume = zeros(Float, nriv),
        h = zeros(Float, nriv),
        h_av = zeros(Float, nriv),
        Δt = Float(tosecond(Δt)),
        its = kw_river_tstep > 0 ? ceil(Int(tosecond(Δt) / kw_river_tstep)) :
              kw_river_tstep,
        width = riverwidth,
        wb_pit = pits[inds_riv],
        alpha_pow = alpha_pow,
        alpha_term = fill(mv, nriv),
        α = fill(mv, nriv),
        eps = Float(1e-3),
        cel = zeros(Float, nriv),
        to_river = zeros(Float, nriv),
        reservoir_index = do_reservoirs ? resindex : fill(0, nriv),
        lake_index = do_lakes ? lakeindex : fill(0, nriv),
        reservoir = do_reservoirs ? reservoirs : nothing,
        lake = do_lakes ? lakes : nothing,
        rivercells = river,
    )

    # setup subdomains for the land and river kinematic wave domain, if nthreads = 1
    # subdomain is equal to the complete domain
    toposort = topological_sort_by_dfs(graph)
    toposort_riv = topological_sort_by_dfs(graph_riv)
    index_pit_land = findall(x -> x == 5, ldd)
    index_pit_river = findall(x -> x == 5, ldd_riv)
    subbas_order, indices_subbas, topo_subbas =
        kinwave_set_subdomains(config, graph, toposort, index_pit_land)
    subriv_order, indices_subriv, topo_subriv =
        kinwave_set_subdomains(config, graph_riv, toposort_riv, index_pit_river)

    modelmap = (vertical = hbv, lateral = (land = olf, river = rf))
    indices_reverse = (
        land = rev_inds,
        river = rev_inds_riv,
        reservoir = isempty(reservoir) ? nothing : reservoir.reverse_indices,
        lake = isempty(lake) ? nothing : lake.reverse_indices,
    )
    writer = prepare_writer(config, reader, modelmap, indices_reverse, x_nc, y_nc, nc)
    close(nc)

    # for each domain save:
    # - the directed acyclic graph (graph),
    # - the traversion order (order),
    # - upstream_nodes,
    # - subdomains for the kinematic wave domains for parallel execution (execution order of
    #   subbasins (subdomain_order), traversion order per subbasin (topo_subdomain) and
    #   Vector indices per subbasin matching the traversion order of the complete domain
    #   (indices_subdomain))
    # - the indices that map it back to the two dimensional grid (indices)

    # for the land domain the x and y length [m] of the grid cells are stored
    # for reservoirs and lakes indices information is available from the initialization
    # functions
    land = (
        graph = graph,
        upstream_nodes = filter_upsteam_nodes(graph, olf.wb_pit),
        subdomain_order = subbas_order,
        topo_subdomain = topo_subbas,
        indices_subdomain = indices_subbas,
        order = toposort,
        indices = inds,
        reverse_indices = rev_inds,
        xl = xl,
        yl = yl,
    )
    river = (
        graph = graph_riv,
        upstream_nodes = filter_upsteam_nodes(graph_riv, rf.wb_pit),
        subdomain_order = subriv_order,
        topo_subdomain = topo_subriv,
        indices_subdomain = indices_subriv,
        order = toposort_riv,
        indices = inds_riv,
        reverse_indices = rev_inds_riv,
    )

    model = Model(
        config,
        (; land, river, reservoir, lake, index_river, frac_toriver),
        (land = olf, river = rf),
        hbv,
        clock,
        reader,
        writer,
    )

    # read and set states in model object if reinit=true
    if reinit == false
        instate_path = joinpath(tomldir, config.state.path_input)
        state_ncnames = ncnames(config.state)
        set_states(instate_path, model, state_ncnames; type = Float)
        # update kinematic wave volume for river and land domain
        @unpack lateral = model
        lateral.land.volume .= lateral.land.h .* lateral.land.width .* lateral.land.dl
        lateral.river.volume .= lateral.river.h .* lateral.river.width .* lateral.river.dl
    end

    # make sure the forcing is already loaded
    # it's fine to run twice, and may help catching errors earlier
    update_forcing!(model)
    if haskey(config.input, "cyclic")
        update_cyclic!(model)
    end
    return model
end

function update(model::Model{N,L,V,R,W}) where {N,L,V<:HBV,R,W}
    @unpack lateral, vertical, network, clock, config = model

    inds_riv = network.index_river
    kinwave_it = get(config.model, "kin_wave_iteration", false)::Bool

    update_forcing!(model)
    if haskey(config.input, "cyclic")
        update_cyclic!(model)
    end

    # vertical hbv concept is updated until snow state, after that (optional)
    # snow transport is possible
    update_until_snow(vertical, config)

    # lateral snow transport
    if get(config.model, "masswasting", false)::Bool
        lateral_snow_transport!(
            vertical.snow,
            vertical.snowwater,
            lateral.land.sl,
            network.land,
        )
    end

    # update vertical hbv concept
    update_after_snow(vertical, config)

    # determine lateral inflow for overland flow based on vertical runoff [mm] from vertical
    # hbv concept
    lateral.land.inwater .=
        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.Δt
    lateral.land.qlat .= lateral.land.inwater ./ lateral.land.dl

    # run kinematic wave for overland flow
    update(
        lateral.land,
        network.land,
        frac_toriver = network.frac_toriver,
        do_iter = kinwave_it,
    )

    # determine lateral inflow (from overland flow) for river flow
    lateral.river.inwater .= copy(lateral.land.to_river[inds_riv])
    lateral.river.qlat .= lateral.river.inwater ./ lateral.river.dl

    # run kinematic wave for river flow
    # check if reservoirs or lakes are defined, the inflow from overland flow is required
    if !isnothing(lateral.river.reservoir) || !isnothing(lateral.river.lake)
        update(
            lateral.river,
            network.river,
            inflow_wb = lateral.land.q_av[inds_riv],
            do_iter = kinwave_it,
            doy = dayofyear(clock.time),
        )
    else
        update(
            lateral.river,
            network.river,
            do_iter = kinwave_it,
            doy = dayofyear(clock.time),
        )
    end

    write_output(model)

    # update the clock
    advance!(clock)

    return model
end


text:
using PortHamiltonian
using LinearAlgebra
using Test

tests = ["ph_discretization","ph_discretization2","ph_discretization_closed",
         "ph_constraint_elim", "ph_discretization_weak",
        "ph_constrained", "ph_coupled_transf"]
for t in tests
    println(t)
    include("$(t).jl")
end

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
#=
-------------------------------------------------------------------------------------------------------------
This file contains the Student template for solving Project 5: the 2D Poisson Equation using the CG method with GGP storage
introduced in F.X. Giraldo's Introduction to Element-based Galerkin Methods using
Tensor-Product Bases: Analysis, Algorithms, and Applications.

The approached follows Algorithm 12.18 in the book.

Written by F.X. Giraldo on July 10, 2021.
Department of Applied Mathematics
Naval Postgraduate School
Monterey, CA 93943

The interpolation points used are the following:
ipoints=1: Lobatto
ipoints=2: Legendre
ipoints=3: Chebyshev
ipoints=4: Equi-spaced

The integration points used are the following:
qpoints=1: Lobatto
qpoints=2: Legendre
-------------------------------------------------------------------------------------------------------------
=#

using Plots
using LinearAlgebra

include("QuadraturePoints.jl")
include("element_matrices.jl")
include("global_matrices.jl")
include("create_grid.jl")
include("compute_metrics.jl")
include("global_matrices.jl")
include("exact_solution.jl")
include("compute_norm.jl")

#Some Constants
DFloat = Float64
machine_zero=eps(DFloat)

#Main Function
function main()

    #-----------------------------Only Change these Input parameters---------------------------------#
    Nex=8
    Ney=8
    N=4
    Q=N+1
    ipoints=1
    qpoints=1
    c=DFloat(1) #Constant in Exact solution
    plot_grid=true
    plot_solution=true
    warp_grid=false
    #-----------------------------Only Change these Input parameters---------------------------------#

    space_method="CG" #only capable of doing CG
    Np=N+1
    Nq=Q+1
    Npts=Np^2
    Nqs=Nq^2
    Ne=Nex*Ney
    Nx=Nex*N+1
    Ny=Ney*N+1
    Npoin=Nx*Ny
    Nboun=2*Nx + 2*(Ny-2)
    @show(N,Q,Ne,Npoin,Nboun,c)

    #Select Interpolation Points
    if ipoints == 1
        (ξ,ω) = QuadraturePoints.lobatto_gauss(Np,DFloat) #Lobatto
    elseif ipoints == 2
        (ξ,ω) = QuadraturePoints.legendre_gauss(Np,DFloat) #Legendre
    elseif ipoints == 3
        (ξ,ω) = QuadraturePoints.chebyshev_gauss(Np,DFloat) #Chebyshev
    elseif ipoints == 4
        (ξ,ω) = QuadraturePoints.equispaced_points(Np,DFloat) #Equi-spaced
    end

    #Select Quadrature Points
    if qpoints == 1
        (ξq,ωq) = QuadraturePoints.lobatto_gauss(Nq,DFloat) #Lobatto
    elseif qpoints == 2
        (ξq,ωq) = QuadraturePoints.legendre_gauss(Nq,DFloat) #Legendre
    end

    #Construct Lagrange Polynomials
    (ψ,dψ) = QuadraturePoints.lagrange_basis(Np,Nq,ξ,ξq,DFloat)

    #Construct Grid and (i,e)->I map
    (coord,intma,iboun) = create_grid(Np,Npoin,Ne,Nboun,Nex,Ney,ξ,plot_grid,warp_grid,DFloat)

    #Construct Metric Terms
    (ξ_x,ξ_y,η_x,η_y,jac) = compute_metrics(coord,intma,ψ,dψ,Ne,Np,Nq,DFloat)

    #-----------------------------Students Add their Functions inside ti_LSRK--------------------------#
    #-----------------------------Students Add their Functions inside ti_LSRK--------------------------#
    #Construct Element Matrices
    (Me,Le) = element_matrices(ψ,dψ,ξ_x,ξ_y,η_x,η_y,jac,ωq,Ne,Np,Nq,DFloat)

    #Construct Global Matrices
    (M,L) = global_matrices(Me,Le,intma,Ne,Np,Npoin,DFloat)
    #-----------------------------Students Add their Functions inside ti_LSRK--------------------------#
    #-----------------------------Students Add their Functions inside ti_LSRK--------------------------#

    #Compute Exact Solution
    (qe,fe) = exact_solution(coord,Npoin,c,DFloat)

    #Apply Dirichlet BCs: Impose exact solution
    R=M*fe
    for i=1:Nboun
        I=iboun[i]
        L[I,:] = zeros(DFloat,Npoin)
        L[I,I] = DFloat(-1)
        R[I]=qe[I]
    end#i

    #Solve Matrix problem
    q=-L\R

    #Compute L2 Norm
    norm2 = compute_norm(q,qe,Npoin,DFloat)
    println("norm.L2 = ",norm2)

    #Plot Solution
    if (plot_solution)
        x=reshape(coord[1,:],(Nx,Ny))
        y=reshape(coord[2,:],(Nx,Ny))
        qq=reshape(q,(Nx,Ny))
        xx=x[:,1]
        yy=y[1,:]

        p1 = contourf(xx,yy,qq)
        title!(p1,"Numerical Solution")
        xlabel!(p1,"x")
        ylabel!(p1,"y")

        qq=reshape(qe,(Nx,Ny))
        p2 = contourf(xx,yy,qq)
        title!(p2,"Exact Solution")
        xlabel!(p2,"x")
        ylabel!(p2,"y")

        figure_1 = plot(p1,p2,layout = (1,2))
        display(figure_1)
    end

    println("Done") #output

end

#----------------------------------#
# Run the main function
#----------------------------------#
main()


text:
module Radon

using ..NESSie
using ..NESSie: _etol, ddot
using Distances: euclidean

export regularyukawacoll, regularyukawacoll!

# =========================================================================================
"""
    regularyukawapot{T}(
        x     ::DenseArray{T,1},
        ξ     ::Vector{T},
        yukawa::T,
              ::Vector{T}=T[]
    )

Computes the regular part of the Yukawa potential, that is, Yukawa minus Laplace:
```math
\\mathcal{G}^Y-\\mathcal{G}^L = \\frac{1}{4π}\\frac{e^{-\\frac{|x - ξ|}{Λ}} - 1}{|x-ξ|}
```

!!! note
    The result is premultiplied by 4π.

# Arguments
 * `yukawa` [Exponent](@ref int-constants) of the Yukawa operator's fundamental solution

# Return type
`T`
"""
function regularyukawapot(
        x     ::DenseArray{T,1},
        ξ     ::Vector{T},
        yukawa::T,
              ::Vector{T}=T[]
    ) where T
    #=== TIME- AND MEMORY-CRITICAL CODE! ===#
    rnorm = euclidean(x, ξ)

    # limit for |x-ξ| → 0
    rnorm <= _etol(T) && return -yukawa

    scalednorm = yukawa * rnorm

    # guard against cancellation
    if scalednorm < .1
        # use alternating series to approximate
        # e^(-c) - 1 = Σ((-c)^i / i!) for i=1 to ∞
        term = -scalednorm
        tolerance = _etol(T) * abs(term)
        tsum = zero(T)
        for i in 1:15
            abs(term) <= tolerance && break

            tsum += term
            term *= -scalednorm / (i+1)
        end
        return tsum / rnorm
    end

    # no danger of cancellation
    (exp(-scalednorm) - 1) / rnorm
end

# =========================================================================================
"""
    ∂ₙregularyukawapot{T}(
        x     ::Vector{T},
        ξ     ::Vector{T},
        yukawa::T,
        normal::Vector{T}
    )

Computes the normal derivative of the regular part of the Yukawa potential, that is, Yukawa
minus Laplace:
```math
\\frac{∂}{∂n} \\frac{1}{4π} \\frac{e^{-\\frac{|x - ξ|}{Λ}} - 1}{|x-ξ|}
= \\frac{1}{4π} \\frac{1 - (1 + Λ^{-1} |x - ξ|)e^{-\\frac{|x - ξ|}{Λ}}}{|x-ξ|²}
\\frac{(x - ξ) ⋅ n}{|x - ξ|}
```

!!! note
    The result is premultiplied by 4π.

# Arguments
 * `yukawa` [Exponent](@ref int-constants) of the Yukawa operator's fundamental solution

# Return type
`T`
"""
function ∂ₙregularyukawapot(
        x     ::Vector{T},
        ξ     ::Vector{T},
        yukawa::T,
        normal::Vector{T}
    ) where T
    #=== TIME- AND MEMORY-CRITICAL CODE! ===#
    rnorm = euclidean(x, ξ)

    # limit for |x-ξ| → 0
    rnorm <= _etol(T) && return yukawa^2 / 2 / T(√3)

    cosovernorm2 = ddot(x, ξ, normal) / rnorm^3
    scalednorm = yukawa * rnorm

    # guard against cancellation
    if scalednorm < .1
        # use alternating series to approximate
        # 1 - (c+1)e^(-c) = Σ((-c)^i * (i-1) / i!) for i=2 to ∞
        term = scalednorm * scalednorm / 2
        tolerance = _etol(T) * abs(term)
        tsum = zero(T)
        for i in 2:16
            abs(term) <= tolerance && break

            tsum += term * (i-1)
            term *= -scalednorm / (i+1)
        end
        return tsum * cosovernorm2
    end

    # no danger of cancellation
    (1 - (1 + scalednorm) * exp(-scalednorm)) * cosovernorm2
end

# =========================================================================================
"""
    setcubpts!{T}(
        dest::Vector{Vector{T}},
        qpts::QuadPts2D{T},
        elem::Triangle{T}
    )

Prepare cubature points for one surface element.

# Return type
`Void`
"""
function setcubpts!(dest::Vector{Vector{T}}, qpts::QuadPts2D{T}, elem::Triangle{T}) where T
    u = elem.v2 - elem.v1
    v = elem.v3 - elem.v1

    # devectorized version of
    # cubpts = [u * qpts.x[i] + v * qpts.y[i] + elem.v1 for i in 1:7]
    for i in 1:qpts.num, j in 1:3
        dest[i][j] = qpts.x[i] * u[j] + qpts.y[i] * v[j] + elem.v1[j]
    end
    nothing
end

# =========================================================================================
"""
    radoncoll!{T}(
            dest    ::DenseArray{T,1},
            elements::Vector{Triangle{T}},
            Ξ       ::Vector{Vector{T}},
            solution::Function,
            yukawa  ::T,
            fvals   ::DenseArray{T,1}
    )

    radoncoll!{T}(
            dest    ::DenseArray{T,2},
            elements::Vector{Triangle{T}},
            Ξ       ::Vector{Vector{T}},
            solution::Function,
            yukawa  ::T
    )

Seven-point Radon cubature [[Rad48]](@ref Bibliography) for a given function and a list of
triangles and observation points `Ξ`. If `dest` is a vector, the function values f for each
surface triangle have to be specified, since each element of the vector represents the dot
product of the corresponding coefficient matrix row and the `fvals` vector.

If you intend computing single/double layer potentials with this function, you might want
to use the shorthand signature `regularyukawacoll!` instead.

!!! note
    The result is premultiplied by 4π.

# Arguments
 * `solution` Fundamental solution; supported functions: `regularyukawapot`,
   `∂ₙregularyukawapot`
 * `yukawa` [Exponent](@ref int-constants) of the Yukawa operator's fundamental solution

# Return type
`Void`
"""
function radoncoll!(
        dest    ::DenseArray{T,1},
        elements::Vector{Triangle{T}},
        Ξ       ::Vector{Vector{T}},
        solution::Function,
        yukawa  ::T,
        fvals   ::DenseArray{T,1}
    ) where T
    #=== MEMORY-CRITICAL CODE! ===#
    @assert length(fvals) == length(elements)
    @assert length(dest) == length(Ξ)

    # pre-allocate memory for cubature points
    qpts = quadraturepoints(Triangle{T})
    cubpts = [zeros(T, 3) for _ in 1:qpts.num]

    @inbounds for eidx in 1:length(elements)
        elem = elements[eidx]
        area = 2 * elem.area
        setcubpts!(cubpts, qpts, elem)

        for oidx in 1:length(Ξ)
            ξ = Ξ[oidx]
            value = zero(T)
            for i in 1:qpts.num
                value += solution(cubpts[i], ξ, yukawa, elem.normal) * qpts.weight[i]
            end
            dest[oidx] += value * area * fvals[eidx]
        end
    end
    nothing
end

function radoncoll!(
        dest    ::DenseArray{T,2},
        elements::Vector{Triangle{T}},
        Ξ       ::Vector{Vector{T}},
        solution::Function,
        yukawa  ::T
    ) where T
    #=== MEMORY-CRITICAL CODE! ===#
    @assert size(dest) == (length(Ξ), length(elements))

    # pre-allocate memory for cubature points
    qpts = quadraturepoints(Triangle{T})
    cubpts = [zeros(T, 3) for _ in 1:qpts.num]

    @inbounds for eidx in 1:length(elements)
        elem = elements[eidx]
        area = 2 * elem.area
        setcubpts!(cubpts, qpts, elem)

        for oidx in 1:length(Ξ)
            ξ = Ξ[oidx]
            value = zero(T)
            for i in 1:qpts.num
                value += solution(cubpts[i], ξ, yukawa, elem.normal) * qpts.weight[i]
            end
            dest[oidx, eidx] = value * area
        end
    end
    nothing
end

# =========================================================================================
"""
    radoncoll{T}(
            ξ       ::Vector{T},
            elem    ::Triangle{T},
            yukawa  ::T,
            solution::Function
    )

Seven-point Radon cubature [[Rad48]](@ref Bibliography) for a given function and a pair of
surface element and observations point.

If you intend computing single/double layer potentials with this function, you might want
to use the shorthand signature `regularyukawacoll` instead.

!!! note
    The result is premultiplied by 4π.

# Arguments
 * `solution` Fundamental solution; supported functions: `regularyukawapot`,
   `∂ₙregularyukawapot`
 * `yukawa` [Exponent](@ref int-constants) of the Yukawa operator's fundamental solution

# Return type
`Void`
"""
function radoncoll(
        ξ       ::Vector{T},
        elem    ::Triangle{T},
        yukawa  ::T,
        solution::Function
    ) where T

    qpts = quadraturepoints(Triangle{T})
    cubpts = [zeros(T, 3) for _ in 1:qpts.num]
    area = 2 * elem.area
    setcubpts!(cubpts, qpts, elem)

    value = zero(T)
    for i in 1:qpts.num
        value += solution(cubpts[i], ξ, yukawa, elem.normal)::T * qpts.weight[i]
    end
    value * area
end

# ========================================================================================
"""
    regularyukawacoll!{T, P <: PotentialType}(
                ::Type{P},
        dest    ::DenseArray{T,1},
        elements::Vector{Triangle{T}},
        Ξ       ::Vector{Vector{T}},
        yukawa  ::T,
        fvals   ::Union{DenseArray{T,1},SubArray{T,1}}
    )

    regularyukawacoll!{T, P <: PotentialType}(
                ::Type{P},
        dest    ::DenseArray{T,2},
        elements::Vector{Triangle{T}},
        Ξ       ::Vector{Vector{T}},
        yukawa  ::T
    )

Computes the regular part of the single or double layer Yukawa potential (that is, Yukawa
minus Laplace) using a seven-point Radon cubature [[Rad48]](@ref Bibliography) for a given
list of triangles and observation points `Ξ`.

The first version of this function uses a vector as destination `dest`, where each element
represents the dot product of the corresponding coefficient matrix row and the `fvals`
vector.

!!! note
    The result is premultiplied by 4π.

# Arguments
 * `yukawa` [Exponent](@ref int-constants) of the Yukawa operator's fundamental solution

# Return type
`Void`
"""
regularyukawacoll!(
            ::Type{SingleLayer},
    dest    ::DenseArray{T,1},
    elements::Vector{Triangle{T}},
    Ξ       ::Vector{Vector{T}},
    yukawa  ::T,
    fvals   ::Union{DenseArray{T,1},SubArray{T,1}}
) where T = radoncoll!(dest, elements, Ξ, regularyukawapot, yukawa, fvals)

regularyukawacoll!(
            ::Type{DoubleLayer},
    dest    ::DenseArray{T,1},
    elements::Vector{Triangle{T}},
    Ξ       ::Vector{Vector{T}},
    yukawa  ::T,
    fvals   ::Union{DenseArray{T,1},SubArray{T,1}}
) where T = radoncoll!(dest, elements, Ξ, ∂ₙregularyukawapot, yukawa, fvals)

regularyukawacoll!(
            ::Type{SingleLayer},
    dest    ::DenseArray{T,2},
    elements::Vector{Triangle{T}},
    Ξ       ::Vector{Vector{T}},
    yukawa  ::T
) where T = radoncoll!(dest, elements, Ξ, regularyukawapot, yukawa)

regularyukawacoll!(
            ::Type{DoubleLayer},
    dest    ::DenseArray{T,2},
    elements::Vector{Triangle{T}},
    Ξ       ::Vector{Vector{T}},
    yukawa  ::T
) where T = radoncoll!(dest, elements, Ξ, ∂ₙregularyukawapot, yukawa)

# ========================================================================================
"""
    regularyukawacoll{T, P <: PotentialType}(
              ::Type{P},
        ξ     ::Vector{T},
        elem  ::Triangle{T},
        yukawa::T
    )

Computes the regular part of the single or double layer Yukawa potential (that is, Yukawa
minus Laplace) using a seven-point Radon cubature [[Rad48]](@ref Bibliography) for a given
triangle and observation point `ξ`.

!!! note
    The result is premultiplied by 4π.

# Arguments
 * `yukawa` [Exponent](@ref int-constants) of the Yukawa operator's fundamental solution

# Return type
`Void`
"""
regularyukawacoll(
          ::Type{SingleLayer},
    ξ     ::Vector{T},
    elem  ::Triangle{T},
    yukawa::T
) where T = radoncoll(ξ, elem, yukawa, regularyukawapot)

regularyukawacoll(
          ::Type{DoubleLayer},
    ξ     ::Vector{T},
    elem  ::Triangle{T},
    yukawa::T
) where T = radoncoll(ξ, elem, yukawa, ∂ₙregularyukawapot)

end # module


text:
function swap(side :: HorizonSide)
    side = HorizonSide((Int(side)+2)%4)
    return side
end

function corner!(r)
    while isborder(r,West)==0
            move!(r,West)
    end
    while isborder(r,Sud)==0
        move!(r,Sud)
    end
end

function GoHome!(r, x, y)
    corner!(r)
    for i = 1:y
        move!(r,Nord)
    end
    for i = 1:x
        move!(r,Ost)
    end
end

function mindhome(r)
    x = 0
    y = 0
    while isborder(r, Sud) == 0
        move!(r,Sud)
        y = y + 1
    end
    while isborder(r, West) == 0
        move!(r,West)
        x = x + 1
    end
    return x, y
end

function put!(r, x, y, n)
    if ((x ÷ n) + (y ÷ n)) % 2 == 0
        putmarker!(r)
    end
end

function chess!(r,side, n, x, y)
    while isborder(r, side) == 0
        move!(r, side)
        if side == Ost
            x = x + 1;
        else
            x = x - 1;
        end
        put!(r, x, y, n)
    end
    return x, y
end

function task_12(r, n)
    z = 0
    m = 0
    z, m = mindhome(r)
    side = Ost;
    x = 0;
    y = 0;
    putmarker!(r)
    while isborder(r,Nord) == 0
        x, y = chess!(r,side, n, x, y)
        move!(r,Nord);
        y = y + 1;
        side = swap(side);
        put!(r, x, y, n)
    end
    chess!(r,side, n, x, y)
    go_home!(r, z, m)
end

n = 5

task_12(r,n)