text:
#=
    The file contains bases searching related code
=#

#------------------------------------------------------------------------------

include("parsik.jl")
include("gizmos.jl")
include("wiedemannchik.jl")
include("../src/structs/dok_sparsik.jl")
include("../src/structs/sparsik.jl")
include("../src/structs/densik.jl")
include("../src/structs/subspacik.jl")

#------------------------------------------------------------------------------

import Primes: nextprime
import Nemo: fmpz

#------------------------------------------------------------------------------

function find_basis_1_γβ(vectors)
    domain = base_ring(first(vectors))
    sz = size(first(vectors))

    # what if..
    hash_vector = random_densik(sz, domain)
    alg = linear_span!(deepcopy(vectors), hash_vector)

    fat_vectors = apply_matrices_inplace!(alg, deepcopy(vectors), ω=0.08)

    for vect in fat_vectors
        eat_sparsik!(alg, vect)
    end

    apply_matrices_inplace!(alg, deepcopy(vectors), ω=1.0)

    return alg
end

#------------------------------------------------------------------------------

function find_basis_1(vectors)
    alg = linear_span!(deepcopy(vectors))
    apply_matrices_inplace!(alg, deepcopy(vectors))
    return alg
end

#------------------------------------------------------------------------------

function find_basis_1_beta(vectors)
    global CNT
    global NONZERO

    # eat all input vectors
    alg = linear_span!(deepcopy(vectors))

    # apply them with the threshold of ω
    fat_vectors = apply_matrices_inplace!(alg, deepcopy(vectors), ω=0.05)

    push!(count_β_before, CNT)
    push!(count_β_skipped, length(fat_vectors))
    CNT = 0

    push!(nonzero_β_before, NONZERO)
    NONZERO = 0

    # eat discarded veced vectors
    for vect in fat_vectors
        eat_sparsik!(alg, vect)
    end

    # apply again, but not no vectors would be discarded

    # CONTROVERSIAL at the best
    apply_matrices_inplace!(alg, deepcopy(vectors), ω=1.0)

    push!(count_β_after, CNT)
    CNT = 0

    return alg
end

#------------------------------------------------------------------------------

#
function find_basis_2(vectors)

    first_vector = first(vectors)
    first_pivot = first_nonzero(first_vector)
    n = dim(first_vector)

    # the size of the current square
    k = 1
    field = base_ring(first_vector)

    # indices of column pivots
    current_pivots = Int[ first_pivot ]

    # long rows, not restricted to current_pivots
    current_vectors = [ deepcopy(first_vector) ]

    # select hashing vector
    # we should adjust the hash_vector density
    # Gleb: I think we should just make it dense, otherwise we can loose some vectors
    hash_vector = random_sparsik(size(first_vector), field, density=1)
    hash_col = [ inner(hash_vector, first_vector) ]

    # this is the rows of the `square` k×k of our big matrix
    # actually, current_rows = current_vectors|current_pivots
    # if `|` is the restriction to some coordinates
    first_row = Sparsik(1, field, Int[1], Dict(1 => first_vector[first_pivot]))
    current_square = from_rows(
        1, 1, field, Int[1],
        Dict(1 => first_row)
    )
    current_rows = Dict(
        1 => deepcopy(first_row)
    )

    first_row_extended = restrict(first_vector, first_pivot, first(hash_col))
    current_rows_extended = Dict(
        1 => first_row_extended
    )

    # ort vector of size 2×1 ( orthogonal for extended rows with hashes )
    ort_vector = find_orthogonal!(current_rows_extended)

    new_vectors = deepcopy(vectors)
    i = 0

    while ! isempty(new_vectors)
        # multiplyiing vectors as long as new ones appear
        # if no new_vectors were produced, we stop

        buffer = deepcopy(new_vectors)
        empty!(new_vectors)

        for new_vect in buffer
            for vect in vectors

                product = apply_vector(new_vect, vect)

                i += 1
                i % 10 == 0 && print(".")

                if !iszero(product)
                    product_hash = inner(hash_vector, product)

                    # product is a square matrix,
                    # product_restricted is a Sparsik vector ( plain , one dimension )
                    # and only columns from `current_pivots` are present in product_restricted
                    product_restricted = restrict(product, current_pivots)
                    product_restricted_extended = restrict(product, current_pivots, product_hash)

                    # (ii) checking for linear dependency
                    if ! iszero(inner(product_restricted_extended, ort_vector))

                        # now we want `y`, such that rows of A|k with coeffs `y`
                        # give us the product|k  ,  so we transpose A|k
                        # ( because otherwise y will be the coeffs of the columns, not of the rows )
                        # (iii) finding dependencies

                        transpose!(current_square)
                        # Gleb: we can do randomized as it is is guaranteed to give a correct result
                        # Alex: seems untrue..., but randomized one is indeed much faster
                        y = square_nonsingular_randomized_wiedemann(current_square, product_restricted)
                        # should be O(1)
                        transpose!(current_square)

                        # new pivot is the column in `product` that is not zero
                        # if we add all `current_rows` with coeffs `y` to `product`
                        # (iv) selecting new pivot
                        # this one can be very slow
                        new_pivot = 0
                        accum = zero(field)
                        col = 0
                        while iszero(accum)
                            # Gleb: I think we should skip the columns that are already pivots, this can be a big speedup
                            # Another option for the future would be to do sort-of binary search by multiplying by some
                            # new hash vectors
                            # Alex: Agreed.
                            #       I tried searching for new pivot by descending the binary tree
                            #       induced by multiple hash-vectors
                            #       It works fine in test and that Algorithm is
                            #       guaranteed to be better asymptotically,
                            #       but I do not want to paste it here just yet
                            col += 1
                            accum = sum(y[j] * current_vectors[j][col]
                                          for j in 1:k) - product[col]
                        end
                        new_pivot = col

                        # Here we are just updating everything: adding new vector,
                        # new pivot, new orthogonal vector and etc.
                        # ...
                        # (v) updating rows
                        k += 1

                        push!(current_vectors, product)
                        current_rows_extended[k] = product_restricted_extended
                        current_rows[k] = product_restricted

                        push!(current_pivots, new_pivot)
                        push!(hash_col, product_hash)

                        # append -> find_orthogonal -> pop
                        # it must be simplified
                        # Gleb: how about making the value of the hash to be
                        # the first coordinate of the restricted vector?
                        # Then there will be non need in this push-pop every time
                        # Done
                        for i in 1 : k
                            append!(current_rows_extended[i],
                                    current_vectors[i][new_pivot])

                            append!(current_rows[i],
                                    current_vectors[i][new_pivot])
                        end
                        ort_vector = find_orthogonal!(current_rows_extended)

                        current_square = from_rows(k, k, field, Array(1 : k), current_rows)

                        push!(new_vectors, product)

                    end
                end
            end
        end

    end

    return linear_span!(current_vectors)
end

#------------------------------------------------------------------------------

# input:
#       vectors - an array-like of AbstractSparsiks over Q
# output:
#       V - a Subspacik consisting of basis vectors
# O(∞)
function find_basis(vectors; used_algorithm=find_basis_1)
    # used_algorithm is assumed not to throw normally
    # and to handle errors by thyself

    V = Subspacik(QQ)
    primes = BigInt[ 2^31 - 1 ]
    i = 0

    @info "generating Algebra with $used_algorithm"

    while true
        prime = last(primes)
        field = GF(fmpz(prime))

        @info "new modulo = $prime"

        xs = [ modular_reduction(x, field)
               for x in vectors ]

        start = time_ns()
        @time V = used_algorithm(xs)
        push!(time_β, timeit(start))

        V = rational_reconstruction(V)

        if check_inclusion!(V, vectors)
            if check_invariance!(V, vectors)
                break
            end
            @warn "invariance check failed.."
        else
            @warn "inclusion check failed.."
        end

        # how to choose the next better?
        push!(primes, nextprime(prime ^ 2 + 1))

        i += 1
        i % 10 == 0 && error("something is wrong")
    end

    @info "constructed a basis for algebra of dim $(dim(V))"

    return V
end

#------------------------------------------------------------------------------

timeit(start) = (time_ns() - start) * 1e-9

times_γβ = []
times_β = []

CNT = 0
count_β_before = []
count_β_skipped = []
count_β_after = []
time_β = []

NONZERO = 0
nonzero_β_before = []

# 100 - 150
function owo()
    for (i, (mfn, mdim, msz, mdata)) in enumerate(load_COO_if(from_dim=45, to_dim=50))

        @info "\n$i-th model : $mfn of dim : $mdim"

        As = map(matr -> from_COO(matr..., QQ), mdata)

        start = time_ns()
        @time V = find_basis(deepcopy(As), used_algorithm=find_basis_1_beta)
        push!(times_γβ, timeit(start))

        #=
        start = time_ns()
        @time V = find_basis(deepcopy(As), used_algorithm=find_basis_1_beta)
        push!(times_β, timeit(start))
        =#

    end
end

# owo()

A1 = from_dense([1 0 0 0; 1 0 0 0; 0 0 0 0; 0 0 0 0], QQ)
A2 = from_dense([1 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0], QQ)
A3 = from_dense([0 0 0 1; 0 0 1 0; 0 1 0 1; 0 1 1 0], QQ)

# A = find_basis([A1, A2, A3],  used_algorithm=find_basis_1_γβ)

# datas = load_COO_if(from_dim=110, to_dim=140)
# As = [map(matr -> from_COO(matr..., QQ), data[4]) for data in datas ]

function test_γβ(idx)
    V = find_basis(deepcopy(As[idx]), used_algorithm=find_basis_1_γβ)
    1
end
function test_β(idx)
    V = find_basis(deepcopy(As[idx]), used_algorithm=find_basis_1_beta)
    1
end

#=
plot(1:length(norm), norm, label="norm")
plot!(1:length(norm), alpha, label="alpha")
plot!(1:length(norm), beta, label="beta")
=#

# Strange SuitSparse
#   ### Big Rationals??
#       bfwa62/bfwa62.mtx, bfwb62/bfwb62.mtx, west0067/west0067.mtx
#   ### What??
#       Journals/Journals.mtx
#
# !! BIOMD0000000259.json
# !! MODEL1604100000.json

# big numbers:
#   Hamrle1/Hamrle1.mtx


text:
"""
    correction_hull(A::IntervalMatrix{T}, t, p) where {T}

Compute the correction term for the convex hull of a point and its linear map
with an interval matrix in order to contain all trajectories of a linear system.

### Input

- `A` -- interval matrix
- `t` -- non-negative time value
- `p` -- order of the approximation

### Output

An interval matrix representing the correction term.

### Algorithm

See Theorem 3 in [1].

[1] M. Althoff, O. Stursberg, M. Buss. Reachability Analysis of Linear Systems
with Uncertain Parameters and Inputs. CDC 2007.
"""
function correction_hull(A::IntervalMatrix{T}, t, p) where {T}
    # initialize interval matrix with zero intervals
    m, n = size(A)
    F = IntervalMatrix(zeros(Interval{T}, m, n))

    A2i = A
    fac_i = 1
    t2i = t
    @inbounds for i in 2:p
        t2i *= t
        left = (one(T) / i^(i/i-1) - one(T) / i^(1/i-1)) * t2i
        itv = Interval(left, zero(T))
        A2i = A2i * A
        fac_i *= i
        F += itv * A2i * (1/fac_i)
    end
    F += _expm_remainder(A, t, p)
    return F
end

ext: jl


text:
int_rules_1_2_2_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.2.2.6*P(x)*(d*x)^m*(a+b*x^2+c*x^4)^p =#
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => Module([q = Expon(~Pq, ~x), k], Antiderivative((~d * ~x) ^ ~m * Sum(Coeff(~Pq, ~x, 2k) * (~x) ^ (2k), [k, 0, q / 2 + 1]) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x) + (1 / ~d) * Antiderivative((~d * ~x) ^ (~m + 1) * Sum(Coeff(~Pq, ~x, 2k + 1) * (~x) ^ (2k), [k, 0, (q - 1) / 2 + 1]) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && Not(PolyQ(~Pq, (~x) ^ 2)))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (1 / 2) * Subst(Antiderivative((~x) ^ ((~m - 1) / 2) * SubstFor((~x) ^ 2, ~Pq, ~x) * (~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~p, ~x), ~x, (~x) ^ 2) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, (~x) ^ 2) && IntegerQ((~m - 1) / 2))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~d * ~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (PolyQ(~Pq, (~x) ^ 2) && IGtQ(~p, -2))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => (1 / (~d) ^ 2) * Antiderivative((~d * ~x) ^ (~m + 2) * ExpandToSum(~Pq / (~x) ^ 2, ~x) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (PolyQ(~Pq, (~x) ^ 2) && EqQ(Coeff(~Pq, ~x, 0), 0))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => With([e = Coeff(~Pq, ~x, 0), f = Coeff(~Pq, ~x, 2), g = Coeff(~Pq, ~x, 4)], (e * (~d * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1)) / (~a * ~d * (~m + 1)) <-- EqQ(~a * f * (~m + 1) - ~b * e * (~m + 2 * ~p + 3), 0) && (EqQ(~a * g * (~m + 1) - ~c * e * (~m + 4 * ~p + 5), 0) && NeQ(~m, -1))) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (PolyQ(~Pq, (~x) ^ 2) && EqQ(Expon(~Pq, ~x), 4))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => ((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ FracPart(~p) / ((4 * ~c) ^ IntPart(~p) * (~b + 2 * ~c * (~x) ^ 2) ^ (2 * FracPart(~p)))) * Antiderivative((~d * ~x) ^ ~m * ~Pq * (~b + 2 * ~c * (~x) ^ 2) ^ (2 * ~p), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && (PolyQ(~Pq, (~x) ^ 2) && (GtQ(Expon(~Pq, (~x) ^ 2), 1) && EqQ((~b) ^ 2 - 4 * ~a * ~c, 0)))
    @apply_utils Antiderivative((~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => With([d = Coeff(PolynomialRemainder((~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x), ~x, 0), e = Coeff(PolynomialRemainder((~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x), ~x, 2)], (~x * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1) * ((~a * ~b * e - d * ((~b) ^ 2 - 2 * ~a * ~c)) - ~c * (~b * d - 2 * ~a * e) * (~x) ^ 2)) / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1) * ExpandToSum((((2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * PolynomialQuotient((~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x) + (~b) ^ 2 * d * (2 * ~p + 3)) - 2 * ~a * ~c * d * (4 * ~p + 5)) - ~a * ~b * e) + ~c * (4 * ~p + 7) * (~b * d - 2 * ~a * e) * (~x) ^ 2, ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (PolyQ(~Pq, (~x) ^ 2) && (GtQ(Expon(~Pq, (~x) ^ 2), 1) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~p, -1) && IGtQ(~m / 2, 0)))))
    @apply_utils Antiderivative((~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~p, ~x) => With([d = Coeff(PolynomialRemainder((~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x), ~x, 0), e = Coeff(PolynomialRemainder((~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x), ~x, 2)], (~x * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1) * ((~a * ~b * e - d * ((~b) ^ 2 - 2 * ~a * ~c)) - ~c * (~b * d - 2 * ~a * e) * (~x) ^ 2)) / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c)) + (1 / (2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c))) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ (~p + 1) * ExpandToSum(2 * ~a * (~p + 1) * ((~b) ^ 2 - 4 * ~a * ~c) * (~x) ^ -(~m) * PolynomialQuotient((~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4, ~x) + (((~b) ^ 2 * d * (2 * ~p + 3) - 2 * ~a * ~c * d * (4 * ~p + 5)) - ~a * ~b * e) * (~x) ^ -(~m) + ~c * (4 * ~p + 7) * (~b * d - 2 * ~a * e) * (~x) ^ (2 - ~m), ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (PolyQ(~Pq, (~x) ^ 2) && (GtQ(Expon(~Pq, (~x) ^ 2), 1) && (NeQ((~b) ^ 2 - 4 * ~a * ~c, 0) && (LtQ(~p, -1) && ILtQ(~m / 2, 0)))))
    #= Antiderivative((~x)^(~m')*(~Pq)*((~a)+(~b')*(~x)^2+(~c')*(~x)^4)^(~p),~x) :=  With([d=Coeff(PolynomialRemainder(x^m*Pq,a+b*x^2+c*x^4,x),x,1), e=Coeff(PolynomialRemainder(x^m*Pq,a+b*x^2+c*x^4,x),x,3)], x^2*(a+b*x^2+c*x^4)^(p+1)*(a*b*e-d*(b^2-2*a*c)-c*(b*d-2*a*e)*x^2)/( 2*a*(p+1)*(b^2-4*a*c)) + 1/(a*(p+1)*(b^2-4*a*c))*Antiderivative(x^m*(a+b*x^2+c*x^4)^(p+1)* ExpandToSum(a*(p+1)*(b^2-4*a*c)*x^(-m)*PolynomialQuotient(x^m*Pq, a+b*x^2+c*x^4,x)+ (b^2*d*(p+2)-2*a*c*d*(2*p+3)-a*b*e)*x^(1-m)+2*c*(p+2)*(b*d-2*a* e)*x^(3-m),x),x)) <-- FreeQ([a,b,c],x) && PolyQ(Pq,x^2) &&  GtQ(Expon(Pq,x^2),1) && NeQ(b^2-4*a*c,0) && LtQ(p,-1) &&  IntegerQ((m-1)/2) =#
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2 + ~(c') * (~x) ^ 4) ^ ~(p'), ~x) => Unintegrable((~d * ~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ 2 + ~c * (~x) ^ 4) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~p], ~x) && PolyQ(~Pq, ~x)
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.2 Quartic/1.2.2.6 P(x) (d x)^m (a+b x^2+c x^4)^p.jl


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
## original example, following here

move_to(cr, 50.0, 75.0);
line_to(cr, 200.0, 75.0);

move_to(cr, 50.0, 125.0);
line_to(cr, 200.0, 125.0);

move_to(cr, 50.0, 175.0);
line_to(cr, 200.0, 175.0);

set_line_width(cr, 30.0);
set_line_cap(cr, Cairo.CAIRO_LINE_CAP_ROUND);
stroke(cr);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_multi_segment_caps.png");


text:
module Immutable
using Toivo
export @immutable, @get!  # shouldn't need to export @get!

macro immutable(ex)
    code_immutable(ex)
end
function code_immutable(ex)
    @expect is_expr(ex, :type, 2)
    sig, body = ex.args
    typename = (is_expr(sig, :(<:), 2) ? sig.args[1] : sig)::Symbol
    fields, types, sigs = Symbol[], {}, {}
    for arg in body.args
        if isa(arg, Symbol)
            push(fields, arg)
            push(types,  quot(Any))
            push(sigs, arg)
        elseif is_expr(arg, :(::), 2)
            push(fields, arg.args[1])
            push(types,  arg.args[2])
            push(sigs, arg)
        end
    end

    instances = Dict()
    esc(quote
        type $sig
            $(body.args...)
            $typename($(sigs...)) = @get!($(quot(instances)), ($(fields...),),
                                          new($(fields...)))
        end
    end)
end

end


text:
# Functions for sparse array view sums - from https://discourse.julialang.org/t/slow-arithmetic-on-views-of-sparse-matrices/3644

rowsum(x) = sum(x,dims=2)
rowsum(x::SubArray{T,2,P}) where {T,P<:SparseMatrixCSC} = (x.parent * sparse(x.indices[2],ones(Int,length(x.indices[2])), ones(Int,length(x.indices[2])), size(x.parent,2),1))[x.indices[1]]
function rowsum(x::SparseMatrixCSC{Bool})
    active = zeros(Int, x.m)
    @inbounds for i in 1:nnz(x)
        active[x.rowval[i]] += 1
    end
    active
end

colsum(x) = sum(x,dims=1)
colsum(x::SparseMatrixCSC{Bool}) = diff(view(x.colptr, 1:x.n+1))
colsum(x::SubArray{T,2,P}) where {T,P<:SparseMatrixCSC} = (sparse(ones(Int,length(x.indices[1])), x.indices[1], ones(Int,length(x.indices[1])),1,size(x.parent,1))*x.parent)[x.indices[2]]

# Functions for finding nonzero rows and columns from Dan Getz, http://stackoverflow.com/questions/43968445/identify-which-rows-or-columns-have-values-in-sparse-matrix

function nzrows(a::SparseMatrixCSC)
    active = falses(a.m)
    @inbounds for i in 1:nnz(a)
        active[a.rowval[i]] = true
    end
    findall(active)
end

function nzcols(a::SparseMatrixCSC)
    res=Vector{Int}()
    foldl((x,y)->(if (x<a.colptr[y]) push!(res,y-1) end; a.colptr[y]),2:a.n+1, init=a.colptr[1])
    res
end

inrange(v,r) = searchsortedlast(v,last(r))>=searchsortedfirst(v,first(r))

function sortedintersecting(v1, v2)
    i,j = firstindex(v1), firstindex(v2)
    @inbounds while i <= lastindex(v1) && j <= lastindex(v2)
        if v1[i] == v2[j] return true
        elseif v1[i] > v2[j] j += 1
        else i += 1
        end
    end
    return false
end

function nzcols(b::SubArray{T,2,P,Tuple{Vector{Int64},Vector{Int64}}} where {T,P<:SparseMatrixCSC}
  )
    brows = sort(unique(b.indices[1]))
    @inbounds return [k
      for (k,i) in enumerate(b.indices[2])
      if b.parent.colptr[i]<b.parent.colptr[i+1] &&
        sortedintersecting(b.parent.rowval[nzrange(b.parent,i)], brows)]
end

function nzcols(b::SubArray{T,2,P,M} where {T,P<:SparseMatrixCSC, M<:Tuple{AbstractUnitRange{Int64},AbstractUnitRange{Int64}}}
  )
    @inbounds return collect(i+1-start(b.indices[2])
      for i in b.indices[2]
      if b.parent.colptr[i]<b.parent.colptr[i+1] &&
        inrange(b.parent.rowval[nzrange(b.parent,i)], b.indices[1]))
end

function nzcols(b::SubArray{T,2,P,M} where {T,P<:SparseMatrixCSC, M<:Tuple{AbstractUnitRange{Int64},Vector{Int64}}}
  )
  @inbounds return [k
    for (k,i) in enumerate(b.indices[2])
    if b.parent.colptr[i]<b.parent.colptr[i+1] &&
        inrange(b.parent.rowval[nzrange(b.parent,i)], b.indices[1])]
end

function nzcols(b::SubArray{T,2,P,M} where {T,P<:SparseMatrixCSC,M<:Tuple{Vector{Int64},AbstractUnitRange{Int64}}}
  )
    brows = sort(unique(b.indices[1]))
    @inbounds return collect(i+1-start(b.indices[2])
      for i in b.indices[2]
      if b.parent.colptr[i]<b.parent.colptr[i+1] &&
        sortedintersecting(b.parent.rowval[nzrange(b.parent,i)], brows))
end

function findin2(inds,v,w)
    i, j = firstindex(v), firstindex(w)
    res = Vector{Int}()
    @inbounds while i<=lastindex(v) && j<=lastindex(w)
        if v[i]==w[j]
            push!(res,inds[i])
            i += 1
        elseif (v[i]<w[j]) i += 1
        else j += 1
        end
    end
    return res
end

function nzrows(b::SubArray{T,2,P,Tuple{Vector{Int64}, U}} where {T,P<:SparseMatrixCSC, U<:Union{M, Vector{Int}}} where {M <: AbstractUnitRange{Int64}}
  )
    active = falses(length(b.indices[1]))
    inds = sortperm(b.indices[1])
    brows = (b.indices[1])[inds]
    @inbounds for c in b.indices[2]
      active[findin2(inds,brows,b.parent.rowval[nzrange(b.parent,c)])] .= true
    end
    return findall(active)
end

function nzrows(b::SubArray{T,2,P,Tuple{M, U}} where {T,P<:SparseMatrixCSC, U<:Union{M, Vector{Int}}} where M <: AbstractUnitRange{Int64}
  )
    active = falses(length(b.indices[1]))
    @inbounds for c in b.indices[2]
        for r in nzrange(b.parent,c)
            if b.parent.rowval[r] in b.indices[1]
                active[b.parent.rowval[r]+1-first(b.indices[1])] = true
            end
        end
    end
    return findall(active)
end

_nnz(x::AbstractMatrix) = nnz(x)
function _nnz(b::SubArray{T,2,P}) where {T,P<:SparseMatrixCSC}
    ret = 0
    @inbounds for c in b.indices[2]
        for r in nzrange(b.parent,c)
            if b.parent.rowval[r] in b.indices[1]
                ret += 1
            end
        end
    end
    return ret
end


text:
"""
    CoolDict([itr])

    CoolDict{V}()

Like a normal `Dict`, but cool. Values are settable and gettable through dot notation as
well as standard `Dict` indexing notation. The only additional restriction is that keys must
be `Symbolic`.

# Examples
```julia-repl
julia> a = CoolDict();

julia> a.a = 1;

julia> a.b = "howdy";

julia> a
CoolDict{Any} with 2 entries:
  :a => 1
  :b => "howdy"

julia> a.b
"howdy"
```

See also: [`Dict`, `NamedTuple`](@ref).
"""
struct CoolDict{V} <: AbstractDict{Symbol,V}
    _data::Dict{Symbol,V}
    CoolDict(dict::Dict{Symbol,V}) where {V} = new{V}(dict)
    CoolDict(dict::Dict{S,V}) where {S,V} = new{V}(Dict{S,V}(dict))
end
CoolDict(; kwargs...) = CoolDict(Dict(kwargs...))
CoolDict(args...) = CoolDict(Dict(args...))
CoolDict(nt::NamedTuple) = CoolDict(pairs(nt))
CoolDict{V}(args...) where {V} = CoolDict(Dict{Symbol,V}(args...))
CoolDict{V}(cd::CoolDict) where {V} = CoolDict(Dict{Symbol,V}(getfield(cd, :_data)))
CoolDict{V}(nt::NamedTuple) where {V} = CoolDict(pairs(map(V, nt)))
CoolDict{V}(; kwargs...) where {V} = CoolDict{V}(pairs(kwargs))


text:
using ProximalOperators
using ProximalAlgorithms
using BenchmarkTools
using LinearAlgebra
using SparseArrays
using Random
using FileIO

const SUITE = BenchmarkGroup()

for (benchmark_name, file_name) in [
    ("Lasso tiny", joinpath(@__DIR__, "data", "lasso_tiny.jld2")),
    ("Lasso small", joinpath(@__DIR__, "data", "lasso_small.jld2")),
    ("Lasso medium", joinpath(@__DIR__, "data", "lasso_medium.jld2")),
]
    for T in [Float64]
        k = "$benchmark_name ($T)"
        SUITE[k] = BenchmarkGroup([k])

        data = load(file_name)

        A = T.(data["A"])
        b = T.(data["b"])
        R = real(T)
        lam = R(data["lambda"])
        m, n = size(A)

        SUITE[k]["ForwardBackward"] = @benchmarkable solver(x0, f=f, g=g) setup=begin
            solver = ProximalAlgorithms.ForwardBackward(tol=1e-6)
            x0 = zeros($T, size($A, 2))
            f = LeastSquares($A, $b)
            g = NormL1($lam)
        end

        SUITE[k]["FastForwardBackward"] = @benchmarkable solver(x0, f=f, g=g) setup=begin
            solver = ProximalAlgorithms.FastForwardBackward(tol=1e-6)
            x0 = zeros($T, size($A, 2))
            f = LeastSquares($A, $b)
            g = NormL1($lam)
        end

        SUITE[k]["ZeroFPR"] = @benchmarkable solver(x0, f=f, A=$A, g=g) setup=begin
            solver = ProximalAlgorithms.ZeroFPR(tol=1e-6)
            x0 = zeros($T, size($A, 2))
            f = Translate(SqrNormL2(), -$b)
            g = NormL1($lam)
        end

        SUITE[k]["PANOC"] = @benchmarkable solver(x0, f=f, A=$A, g=g) setup=begin
            solver = ProximalAlgorithms.PANOC(tol=1e-6)
            x0 = zeros($T, size($A, 2))
            f = Translate(SqrNormL2(), -$b)
            g = NormL1($lam)
        end

        SUITE[k]["DouglasRachford"] = @benchmarkable solver(x0, f=f, g=g, gamma=$R(1)) setup=begin
            solver = ProximalAlgorithms.DouglasRachford(tol=1e-6)
            x0 = zeros($T, size($A, 2))
            f = LeastSquares($A, $b)
            g = NormL1($lam)
        end

        SUITE[k]["DRLS"] = @benchmarkable solver(x0, f=f, g=g, Lf=Lf) setup=begin
            solver = ProximalAlgorithms.DRLS(tol=1e-6)
            x0 = zeros($T, size($A, 2))
            f = LeastSquares($A, $b)
            Lf = opnorm(($A)' * $A)
            g = NormL1($lam)
        end

        SUITE[k]["AFBA-1"] = @benchmarkable solver(x0, y0, f=f, g=g, beta_f=beta_f) setup=begin
            beta_f = opnorm($A)^2
            solver = ProximalAlgorithms.AFBA(theta=$R(1), mu=$R(1), tol=$R(1e-6))
            x0 = zeros($T, size($A, 2))
            y0 = zeros($T, size($A, 2))
            f = LeastSquares($A, $b)
            g = NormL1($lam)
        end

        SUITE[k]["AFBA-2"] = @benchmarkable solver(x0, y0, h=h, L=$A, g=g) setup=begin
            beta_f = opnorm($A)^2
            solver = ProximalAlgorithms.AFBA(theta=$R(1), mu=$R(1), tol=$R(1e-6))
            x0 = zeros($T, size($A, 2))
            y0 = zeros($T, size($A, 1))
            h = Translate(SqrNormL2(), -$b)
            g = NormL1($lam)
        end

        SUITE[k]["SFISTA"] = @benchmarkable solver(y0, f=f, Lf=Lf, h=h) setup=begin
            solver = ProximalAlgorithms.SFISTA(tol=$R(1e-3))
            y0 = zeros($T, size($A, 2))
            f = LeastSquares($A, $b)
            h = NormL1($lam)
            Lf = opnorm($A)^2
        end
    end
end


text:
module Db

# Write your package code here.

export my_db_function

my_db_function() = "I'm a function inside src/Db.jl!"

end


text:
# MathOptInterface results

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
	return model.inner.status
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return string(model.inner.status)
end

function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
    if state_is_optimal(model.inner.status; allow_almost=model.inner.options.allow_almost_solved)
        return MOI.FEASIBLE_POINT
    else
        return MOI.INFEASIBLE_POINT
    end
end

function MOI.get(model::Optimizer, ::MOI.DualStatus)
    if state_is_optimal(model.inner.status; allow_almost=model.inner.options.allow_almost_solved)
        return MOI.FEASIBLE_POINT
    else
        return MOI.INFEASIBLE_POINT
    end
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveValue)
    if model.inner.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    return model.inner.objval
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveBound)
    if model.inner.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    return model.inner.best_bound
end

function MOI.get(model::Optimizer, ::MOI.RelativeGap)
    if model.inner.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    if isnan(model.inner.objval) || isnan(model.inner.best_bound)
        return NaN
    end
    return abs(model.inner.best_bound-model.inner.objval)/abs(model.inner.objval)
end

function MOI.get(model::Optimizer, ::MOI.SolveTime)
    if model.inner.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    return model.inner.soltime
end

function MOI.get(model::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if model.inner.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    check_inbounds(model, vi)
    return model.inner.solution[vi.value]
end


text:
module PairPotentials

include("../ext_imports.jl")

using LinearAlgebra: norm, dot
using JuLIP.Potentials: ZList, SZList, @pot, @D,
                        PairPotential, SimplePairPotential
using StaticArrays: SMatrix

include("pair_basis.jl")

include("pair_pot.jl")

include("repulsion.jl")

end


text:
module VTUFileHandler

using Base64, CodecZlib
#include(joinpath(".","XMLParser.jl"))
using XMLParser

include(joinpath(".","VTUFileHandler","defs.jl"))

struct VTUHeader{T<:Union{UInt32,UInt64}}
	num_blocks::T
	blocksize::T
	last_blocksize::T
	compressed_blocksizes::Vector{T}
	function VTUHeader(::Type{T}) where {T<:Union{UInt32,UInt64}}
		return new{T}(zero(T),zero(T),zero(T),T[])
	end
	function VTUHeader(::Type{T},input::Vector{UInt8}) where {T<:Union{UInt32,UInt64}}
		dat = reinterpret(T,input)
		num_blocks = dat[1]
		blocksize = dat[2]
		last_blocksize = dat[3]
		compressed_blocksizes=Vector{T}(undef,num_blocks)
		for (j,i) = enumerate(4:length(dat))
			compressed_blocksizes[j] = dat[i]
		end
		return new{T}(num_blocks,blocksize,last_blocksize,compressed_blocksizes)
	end
	function VTUHeader(a::T,b::T,c::T,d::Vector{T}) where {T<:Union{UInt32,UInt64}}
		return new{T}(a,b,c,d)
	end
end

include(joinpath(".","VTUFileHandler","vtuheader_utils.jl"))
include(joinpath(".","VTUFileHandler","io.jl"))

struct VTUDataField{T}
	dat::Vector{T}
	VTUDataField(x::Vector{T}) where {T} = new{T}(x)
	VTUDataField(::Type{T}) where {T} = new{T}()
	VTUDataField(x::Base.ReinterpretArray{T, A, B, Vector{C}, D}) where {T,A,B,C,D} = new{T}(x)
end

mutable struct VTUData
	names::Vector{String}
	header::Vector{VTUHeader}
	data::Vector{VTUDataField}
	interp_data::Vector{VTUDataField{Float64}}
	idat::Vector{Int}
	VTUData(names,header,data,interp_data,idat) = new(names,header,data,interp_data,idat)
	function VTUData(dataarrays::Vector{XMLElement},appendeddata::Vector{XMLElement},headertype::Union{Type{UInt32},Type{UInt64}},offsets::Vector{Int},compressed_dat::Bool)
		names = Vector{String}()
		data = Vector{VTUDataField}()
		header = Vector{VTUHeader}()
		for (i,el) in enumerate(dataarrays)
			_type = replace(getAttribute(el,"type"),"\""=>"")
			type = eval(Meta.parse(_type))
			_format = getAttribute(el,"format")
			_name = getAttribute(el,"Name")
			if findfirst(x->x==_name,vtukeywords.interpolation_keywords) != nothing	|| findfirst(x->x==_name,vtukeywords.uncompress_keywords) != nothing
				if _format == "appended" || _format == "\"appended\""
					vtuheader = readappendeddata!(data,i,appendeddata,offsets,type,headertype,compressed_dat)
				else
					vtuheader = readdataarray!(data,el,type,headertype,compressed_dat)
				end
				push!(names,_name)
				push!(header,vtuheader)
			end
		end
		idat = findall(map(y->findfirst(x->x==y,vtukeywords.interpolation_keywords)!=nothing,names))
		interp_data = data[idat]
		return new(names,header,data,interp_data,idat)
	end
end

include(joinpath(".","VTUFileHandler","vtudata_utils.jl"))
include(joinpath(".","VTUFileHandler","vtudata_math.jl"))

function deletefieldata!(xmlroot)
	dataarrays = getElements(xmlroot,"DataArray")
	fielddata = getElements(xmlroot,"FieldData")[1]
	_offset = parse(Int,replace(getAttribute(dataarrays[length(fielddata.content)+1],"offset"),"\""=>""))
	empty!(fielddata.content)
	els = getElements(xmlroot,"DataArray")
	for el in els
		setAttribute(el,"offset",string(parse(Int,replace(getAttribute(el,"offset"),"\""=>""))-_offset))
	end
	appendeddata = getElements(xmlroot,"AppendedData")
	appendeddata[1].content[1] = String(deleteat!(collect(appendeddata[1].content[1]),2:_offset+1))
end

mutable struct VTUFile
	name::String
	xmlroot::XMLElement
	dataarrays::Vector{XMLElement}
	appendeddata::Vector{XMLElement}
	headertype::Union{Type{UInt32},Type{UInt64}}
	offsets::Vector{Int}
	data::VTUData
	compressed_dat::Bool
	VTUFile(name,xmlroot,dataarrays,appendeddata,headertype,offsets,data) = new(name,xmlroot,dataarrays,appendeddata,headertype,offsets,data,true)
	VTUFile(name,xmlroot,dataarrays,appendeddata,headertype,offsets,data,compr_dat) = new(name,xmlroot,dataarrays,appendeddata,headertype,offsets,data,compr_dat)
	function VTUFile(name::String)
		state = IOState(name);
		xmlroot = readXMLElement(state);
		if !isempty(getElements(xmlroot,"FieldData"))
			deletefieldata!(xmlroot)
		end
		dataarrays = getElements(xmlroot,"DataArray")
		appendeddata = getElements(xmlroot,"AppendedData")
		vtkfile = getElements(xmlroot,"VTKFile")
		@assert length(vtkfile) == 1
		compressed_dat = false
		if hasAttributekey(vtkfile[1],"compressor")
			compressed_dat = true
		end
		#headertype = eval(Symbol(getAttribute(first(vtkfile),"header_type")))
		attr = replace(getAttribute(first(vtkfile),"header_type"),"\""=>"")
		headertype = eval(Meta.parse(attr))
		offsets = Vector{Int}()
		for el in dataarrays
			if hasAttributekey(el,"offset")
				_offset = replace(getAttribute(el,"offset"),"\""=>"")
				offset = parse(Int,_offset)
				push!(offsets,offset)
			else
				push!(offsets,0)
			end
		end
		retval = new(name,xmlroot,dataarrays,appendeddata,headertype,offsets,VTUData(dataarrays,appendeddata,headertype,offsets,compressed_dat),compressed_dat)
		return retval
	end
end

include(joinpath(".","VTUFileHandler","utils.jl"))
include(joinpath(".","VTUFileHandler","vtufile_math.jl"))

end #module VTUHandler

max_stars_repo_name: nagelt/VTUFileHandler

max_issues_repo_name: nagelt/VTUFileHandler

max_forks_repo_name: nagelt/VTUFileHandler


text:
using Plots
pgfplotsx()
# using LaTeXStrings

using JSON
using FileIO
using DelimitedFiles

# --- NeuralNet realated modules --------
# using Flux
# using CUDA
# using BSON:@load

include("./digit_mods.jl")
using .digit_mods

# EPIC COLOR PALLETTE
col = [RGB(i) for i in [0x00ab9f,0x7409b7,0xd2cf00,0xff6f00] ] # mah precious colors

# DATA LOADING

data = readdlm("./digitalization_folder/exp_points.csv",',')

# println(data)
# println(typeof(data["exp_up"]),"\n",data["exp_up"])
# println(digit_mods.interpolate([0.1 0.5;0.3 0.25], exp_up"]),"\n",data["exp_up"])[0.15,0.2]))

exp_down = str_el([data[2:end,1],data[2:end,2]],concat=true)
exp_up = str_el([data[2:end,3],data[2:end,4]],concat=true)
x = str_el([data[2:end,5]],concat=false)

int_down = interpolate(exp_down,x)
int_up = interpolate(exp_up,x)

meshlm = Array{Float32}(undef,(2,49))
mesh1  = Array{Float32}(undef,(2,49))
nasa   = Array{Float32}(undef,(2,49))

popul("/home/freshstart/Dropbox/DTC/testRAE/",meshlm)
popul("/home/freshstart/DiplomaThesisData/RAE_2822_baseline1/",mesh1)
popul("/home/freshstart/DiplomaThesisData/RAE_2822_baseline/",nasa)

m_delta_up =   [[100*abs(int_up[i]-meshlm[1,i])/maximum([abs(i) for i in int_up]) for i in 1:length(x)],
                [100*abs(int_up[i]-mesh1[1,i])/maximum([abs(i) for i in int_up]) for i in 1:length(x)],
                [100*abs(int_up[i]-nasa[1,i])/maximum([abs(i) for i in int_up]) for i in 1:length(x)]]
m_delta_down = [[100*abs(int_down[i]-meshlm[2,i])/maximum([abs(i) for i in int_down]) for i in 1:length(x)],
                [100*abs(int_down[i]-mesh1[2,i])/maximum([abs(i) for i in int_down]) for i in 1:length(x)],
                [100*abs(int_down[i]-nasa[2,i])/maximum([abs(i) for i in int_down]) for i in 1:length(x)]]
delta_up =     [[100*(int_up[i]-meshlm[1,i])/int_up[i] for i in 1:length(x)],
                [100*(int_up[i]-mesh1[1,i])/int_up[i] for i in 1:length(x)],
                [100*(int_up[i]-nasa[1,i])/int_up[i] for i in 1:length(x)]]
delta_down =   [[100*abs(int_down[i]-meshlm[2,i])/int_down[i] for i in 1:length(x)],
                [100*abs(int_down[i]-mesh1[2,i])/int_down[i] for i in 1:length(x)],
                [100*abs(int_down[i]-nasa[2,i])/int_down[i] for i in 1:length(x)]]

up=plot(title="Cp distribution divegence from the experiment for the top face [%]",
        ylabel="%",yguidefontrotation=-90,xlabel="x-axis normalized",minorgrid=true)
# plot!(up,x[2:end],delta_up[1][2:end],c=col[2],label="Mesh #1")
plot!(up,x[2:end],m_delta_up[2][2:end],c=col[1],label = "Mesh #2")
plot!(up,x[2:end],m_delta_up[3][2:end],c=col[4],label = "NASA Mesh")
plot!(up,[0],[0],c=RGB(0xffffff),label = L"\frac{|Cp_{Mesh\ i} - Cp_{exp\ i}|}{max(|Cp_{exp}|)}")
plot!(up,[0],[0],c=RGB(0xffffff),label = " ")

bot=plot(title="Cp distribution divegence from the experiment for the bottom face [%]",
        ylabel="%",yguidefontrotation=-90,xlabel="x-axis normalized",minorgrid=true)
# plot!(bot,x[2:end],delta_down[1][2:end],c=col[2],label="Mesh #1")
plot!(bot,x[2:end],m_delta_down[2][2:end],c=col[1],label = "Mesh #2")
plot!(bot,x[2:end],m_delta_down[3][2:end],c=col[4],label = "NASA Mesh")
plot!(bot,[0],[0],c=RGB(0xffffff),label = L"\frac{|Cp_{Mesh\ i} - Cp_{exp\ i}|}{max(|Cp_{exp}|)}")
plot!(bot,[0],[0],c=RGB(0xffffff),label = " ")

savefig(up,"/home/freshstart/DiplomaThesisData/images/Cp_delta_u.pdf")
savefig(bot,"/home/freshstart/DiplomaThesisData/images/Cp_delta_b.pdf")

_lw=1
s=(800,600)
plt=plot(size=s,legend=:topright,yflip=true,title="Comparison between the different meshes",
        ylabel="Cp",yguidefontrotation=-90,xlabel="x-axis normalized",background=:gray99,
        thickness_scaling=1.5,minorgrid=true)

plot!(plt,exp_up[:,1],exp_up[:,2],color=col[4],lw=_lw,line=:dot ,label="exp top side")
plot!(exp_down[:,1],exp_down[:,2],color=col[4],line=:dash,lw=_lw,label="exp bottom side")
#----------------
plot!(x,meshlm[1,:],color=col[2],lw=_lw,line=:dashdot ,label="Mesh #1 top side")
plot!(x,meshlm[2,:],color=col[2],lw=_lw,label="Mesh #1 bottom side")
#---------------
plot!(x,nasa[1,:],color=col[1],lw=_lw,line=:dashdot ,label="NASA mesh top side")
plot!(x,nasa[2,:],color=col[1],lw=_lw,label="NASA mesh top side")
#---------------
plot!(x,mesh1[1,:],color=col[3],lw=_lw,line=:dashdot ,label="Mesh #2 top side")
plot!(x,mesh1[2,:],color=col[3],lw=_lw,label="Mesh #2 bottom side")
savefig(plot!(),"/home/freshstart/DiplomaThesisData/images/Meshes_versus.pdf")

# plot!(x,mesh1[1,:],c=col[2])
# plot!(x,mesh1[2,:],c=col[2])
# plot!(x,int_down,c=col[3])
# plot!(x,int_up,c=col[3])


text:
module LLVM

using Unicode
using Printf
using Libdl

## source code includes

include("util.jl")
include("base.jl")

const libllvm = Ref{String}()
module API
using CEnum
using ..LLVM
using ..LLVM: libllvm, @runtime_ccall
const off_t = Csize_t
libdir = joinpath(@__DIR__, "..", "lib")
include(joinpath(libdir, "libLLVM_common.jl"))
include(joinpath(libdir, "libLLVM_h.jl"))
include(joinpath(libdir, "libLLVM_extra.jl"))
end

# LLVM API wrappers
include("support.jl")
include("types.jl")
include("passregistry.jl")
include("init.jl")
include("core.jl")
include("linker.jl")
include("irbuilder.jl")
include("analysis.jl")
include("moduleprovider.jl")
include("pass.jl")
include("passmanager.jl")
include("execution.jl")
include("buffer.jl")
include("target.jl")
include("targetmachine.jl")
include("datalayout.jl")
include("ir.jl")
include("bitcode.jl")
include("transform.jl")
include("debuginfo.jl")

has_orc_v1() = v"8" <= LLVM.version() < v"12"
include("orc.jl")

include("interop.jl")

include("deprecated.jl")

## initialization

function __init__()
    # find the libLLVM loaded by Julia
    #
    # we only support working with the copy of LLVM that Julia uses, because we have
    # additional library calls compiled in the Julia binary which cannot be used with
    # another copy of LLVM. loading multiple copies of LLVM typically breaks anyhow.
    libllvm[] = if VERSION >= v"1.6.0-DEV.1429"
        path = Base.libllvm_path()
        if path === nothing
            error("""Cannot find the LLVM library loaded by Julia.
                     Please use a version of Julia that has been built with USE_LLVM_SHLIB=1 (like the official binaries).
                     If you are, please file an issue and attach the output of `Libdl.dllist()`.""")
        end
        String(path)
    else
        libllvm_paths = filter(Libdl.dllist()) do lib
            occursin(r"LLVM\b", basename(lib))
        end
        if isempty(libllvm_paths)
            error("""Cannot find the LLVM library loaded by Julia.
                     Please use a version of Julia that has been built with USE_LLVM_SHLIB=1 (like the official binaries).
                     If you are, please file an issue and attach the output of `Libdl.dllist()`.""")
        end
        if length(libllvm_paths) > 1
            error("""Multiple LLVM libraries loaded by Julia.
                     Please file an issue and attach the output of `Libdl.dllist()`.""")
        end
        first(libllvm_paths)
    end

    @debug "Using LLVM $(version()) at $(Libdl.dlpath(libllvm[]))"
    if version() !== runtime_version()
        @error "Using a different version of LLVM ($(runtime_version())) than the one shipped with Julia ($(version())); this is unsupported"
    end

    _install_handlers()
    _install_handlers(GlobalContext())
end

end


text:
using Test
using DataStructures: OrderedDict
using Printf: @sprintf
using QuantumLattices.Mathematics.QuantumNumbers
using QuantumLattices.Interfaces: ⊕,⊗,dimension,expand,permute,decompose,regularize
import Base.Iterators: Iterators

@testset "SQN" begin
    @test SQN|>fieldnames==(:Sz,)
    @test SQN|>periods==(Inf,)
end

@testset "PQN" begin
    @test PQN|>fieldnames==(:N,)
    @test PQN|>periods==(Inf,)
end

@testset "SPQN" begin
    @test SPQN|>fieldnames==(:N,:Sz)
    @test SPQN|>periods==(Inf,Inf)
end

@testset "Z2QN" begin
    @test Z2QN|>fieldnames==(:N,)
    @test Z2QN|>periods==(2,)
end

@testset "SQNS" begin
    @test SQNS(0.5)==AbelianNumbers('C',[SQN(-0.5),SQN(0.5)],[0,1,2],qnsindptr)
end

@testset "PQNS" begin
    @test PQNS(1.0)==AbelianNumbers('C',[PQN(0.0),PQN(1.0)],[0,1,2],qnsindptr)
end

@testset "SzPQNS" begin
    @test SzPQNS(0.5)==AbelianNumbers('C',[SPQN(0.0,0.0),SPQN(1.0,0.5)],[0,1,2],qnsindptr)
end

@testset "SPQNS" begin
    @test SPQNS(0.5)==AbelianNumbers('C',[SPQN(0.0,0.0),SPQN(1.0,-0.5),SPQN(1.0,0.5),SPQN(2.0,0.0)],[0,1,2,3,4],qnsindptr)
end

@testset "Z2QNS" begin
    @test Z2QNS()==AbelianNumbers('C',[Z2QN(0.0),Z2QN(1.0)],[0,1,2],qnsindptr)
end

@abeliannumber "CN" (:N,) (Inf,)
@abeliannumber "Z4" (:Z,) (4,)
@abeliannumber "CNZ4" (:N,:Z) (Inf,4)

@testset "regularize" begin
    @test regularize(CNZ4,[1.5,5.0])==[1.5,1.0]
    @test regularize(CNZ4,[4.0,-1.0])==[4.0,3.0]
    @test regularize(CNZ4,[1.5 4.0; 5.0 -1.0])==[1.5 4.0; 1.0 3.0]
end

@testset "AbelianNumbers" begin
    qn1,qn2=CNZ4(+1.0,+3.0),CNZ4(-1.0,+1.0)
    qns1=AbelianNumbers('U',[qn1,qn2],[0,2,5],qnsindptr)
    qns2=AbelianNumbers('U',[qn1,qn2],[2,3],qnscounts)
    @test isequal(qns1,qns2)
    qns3=AbelianNumbers(OrderedDict(qn1=>2,qn2=>3))
    qns4=AbelianNumbers(OrderedDict(qn1=>1:2,qn2=>3:5))
    @test qns1==qns2==qns3==qns4

    qns=AbelianNumbers('U',[qn1,qn2],[0,3,5],qnsindptr)
    @test qns|>dimension==5
    @test qns|>string=="QNS(2,5)"
    @test qns|>length==2
    @test qns|>eltype==CNZ4
    @test qns|>typeof|>eltype==CNZ4
    @test @sprintf("%s",qns)=="QNS(CNZ4(1.0,3.0)=>1:3,CNZ4(-1.0,1.0)=>4:5)"
    @test qns[1]==qn1
    @test qns[2]==qn2
    @test qns[1:2]==AbelianNumbers('G',[qn1,qn2],[0,3,5],qnsindptr)
    @test qns[[2,1]]==AbelianNumbers('G',[qn2,qn1],[0,2,5],qnsindptr)
    @test qns|>collect==[qn1,qn2]
    @test qns|>Iterators.reverse|>collect==[qn2,qn1]
    @test qns|>keys|>collect==[qn1,qn2]
    @test values(qns,qnsindptr)|>collect==[1:3,4:5]
    @test values(qns,qnscounts)|>collect==[3,2]
    @test pairs(qns,qnsindptr)|>collect==[qn1=>1:3,qn2=>4:5]
    @test pairs(qns,qnscounts)|>collect==[qn1=>3,qn2=>2]
end

@testset "toordereddict" begin
    qn1,qn2=CNZ4(1.0,2.0),CNZ4(2.0,3.0)
    qns=AbelianNumbers('U',[qn1,qn2],[2,3],qnscounts)
    @test toordereddict(qns,qnsindptr)==OrderedDict(qn1=>1:2,qn2=>3:5)
    @test toordereddict(qns,qnscounts)==OrderedDict(qn1=>2,qn2=>3)
end

@testset "arithmetic" begin
    qn=CNZ4(1.0,3.0)
    @test qn|>dimension==1
    @test qn|>typeof|>dimension==1
    @test +qn==CNZ4(+1.0,+3.0)
    @test -qn==CNZ4(-1.0,+1.0)
    @test qn*4==4*qn==CNZ4(4.0,0.0)
    @test qn^3==CNZ4(3.0,1.0)
    @test kron(CNZ4,CN(1.0),Z4(3.0))==qn

    qn1,qn2=CNZ4(1.0,2.0),CNZ4(2.0,3.0)
    @test qn1+qn2==CNZ4(3.0,1.0)
    @test qn1-qn2==CNZ4(-1.0,3.0)
    @test kron(qn1,qn2,signs=(+1,+1))==CNZ4(+3.0,1.0)
    @test kron(qn1,qn2,signs=(+1,-1))==CNZ4(-1.0,3.0)
    @test kron(qn1,qn2,signs=(-1,+1))==CNZ4(+1.0,1.0)
    @test kron(qn1,qn2,signs=(-1,-1))==CNZ4(-3.0,3.0)
    @test union(qn1,qn2,signs=(+1,+1))==AbelianNumbers('G',[+qn1,+qn2],[0,1,2],qnsindptr)
    @test union(qn1,qn2,signs=(+1,-1))==AbelianNumbers('G',[+qn1,-qn2],[0,1,2],qnsindptr)
    @test union(qn1,qn2,signs=(-1,+1))==AbelianNumbers('G',[-qn1,+qn2],[0,1,2],qnsindptr)
    @test union(qn1,qn2,signs=(-1,-1))==AbelianNumbers('G',[-qn1,-qn2],[0,1,2],qnsindptr)
    @test ⊗(qn1,qn2)==kron(qn1,qn2)
    @test ⊕(qn1,qn2)==union(qn1,qn2)

    qns=AbelianNumbers('U',[qn1,qn2],[0,2,4],qnsindptr)
    @test +qns==AbelianNumbers('U',[+qn1,+qn2],[0,2,4],qnsindptr)
    @test -qns==AbelianNumbers('U',[-qn1,-qn2],[0,2,4],qnsindptr)
    @test qns*3==3*qns==AbelianNumbers('U',[qn1*3,qn2*3],[0,2,4],qnsindptr)
    @test qns^2==AbelianNumbers('G',[qn1+qn1,qn1+qn2,qn1+qn1,qn1+qn2,qn2+qn1,qn2+qn2,qn2+qn1,qn2+qn2],[2,2,2,2,2,2,2,2],qnscounts)

    @test qns+qn==qn+qns==AbelianNumbers('U',[qn1+qn,qn2+qn],[0,2,4],qnsindptr)
    @test qns-qn==AbelianNumbers('U',[qn1-qn,qn2-qn],[0,2,4],qnsindptr)
    @test qn-qns==AbelianNumbers('U',[qn-qn1,qn-qn2],[0,2,4],qnsindptr)

    qns1=AbelianNumbers('U',[+qn1,-qn2],[0,2,4],qnsindptr)
    qns2=AbelianNumbers('U',[-qn1,+qn2],[0,3,4],qnsindptr)
    @test union(qns1,qns2,signs=(+1,+1))==AbelianNumbers('G',[+qn1,-qn2,-qn1,+qn2],[2,2,3,1],qnscounts)
    @test union(qns1,qns2,signs=(+1,-1))==AbelianNumbers('G',[+qn1,-qn2,+qn1,-qn2],[2,2,3,1],qnscounts)
    @test union(qns1,qns2,signs=(-1,+1))==AbelianNumbers('G',[-qn1,+qn2,-qn1,+qn2],[2,2,3,1],qnscounts)
    @test union(qns1,qns2,signs=(-1,-1))==AbelianNumbers('G',[-qn1,+qn2,+qn1,-qn2],[2,2,3,1],qnscounts)
    @test kron(qns1,qns2,signs=(+1,+1))==AbelianNumbers('G',[+qn1-qn1,+qn1+qn2,+qn1-qn1,+qn1+qn2,-qn2-qn1,-qn2+qn2,-qn2-qn1,-qn2+qn2],[3,1,3,1,3,1,3,1],qnscounts)
    @test kron(qns1,qns2,signs=(+1,-1))==AbelianNumbers('G',[+qn1+qn1,+qn1-qn2,+qn1+qn1,+qn1-qn2,-qn2+qn1,-qn2-qn2,-qn2+qn1,-qn2-qn2],[3,1,3,1,3,1,3,1],qnscounts)
    @test kron(qns1,qns2,signs=(-1,+1))==AbelianNumbers('G',[-qn1-qn1,-qn1+qn2,-qn1-qn1,-qn1+qn2,+qn2-qn1,+qn2+qn2,+qn2-qn1,+qn2+qn2],[3,1,3,1,3,1,3,1],qnscounts)
    @test kron(qns1,qns2,signs=(-1,-1))==AbelianNumbers('G',[-qn1+qn1,-qn1-qn2,-qn1+qn1,-qn1-qn2,+qn2+qn1,+qn2-qn2,+qn2+qn1,+qn2-qn2],[3,1,3,1,3,1,3,1],qnscounts)
    @test ⊕(qns1,qns2)==union(qns1,qns2)
    @test ⊗(qns1,qns2)==kron(qns1,qns2)
end

@testset "sort" begin
    qn1,qn2,qn3,qn4=CNZ4(3.0,2.0),CNZ4(4.0,1.0),CNZ4(4.0,2.0),CNZ4(3.0,3.0)
    oldqns=AbelianNumbers('G',[qn1,qn2,qn3,qn4,qn2,qn3,qn1,qn4],[2,1,3,4,1,1,2,3],qnscounts)
    newqns,permutation=sort(oldqns)
    @test newqns==AbelianNumbers('C',[qn1,qn4,qn2,qn3],[4,7,2,4],qnscounts)
    @test permutation==[1,2,13,14,7,8,9,10,15,16,17,3,11,4,5,6,12]
end

@testset "findall" begin
    qn1,qn2,qn3=CNZ4(2.0,3.0),CNZ4(1.0,2.0),CNZ4(0.0,0.0)
    qns=AbelianNumbers('C',[qn2,qn1],[2,3],qnscounts)
    @test findall(qn1,qns,qnscompression)==[2]
    @test findall(qn2,qns,qnscompression)==[1]
    @test findall(qn3,qns,qnscompression)==[]
    @test findall(qn1,qns,qnsexpansion)==[3,4,5]
    @test findall(qn2,qns,qnsexpansion)==[1,2]
    @test findall(qn3,qns,qnsexpansion)==[]

    qns=AbelianNumbers('G',[qn1,qn2,qn1],[2,2,1],qnscounts)
    @test findall(qn1,qns,qnscompression)==[1,3]
    @test findall(qn2,qns,qnscompression)==[2]
    @test findall(qn3,qns,qnscompression)==[]
    @test findall(qn1,qns,qnsexpansion)==[1,2,5]
    @test findall(qn2,qns,qnsexpansion)==[3,4]
    @test findall(qn3,qns,qnsexpansion)==[]
end

@testset "filter" begin
    qn1,qn2=CNZ4(1.0,2.0),CNZ4(3.0,0.0)
    qns=AbelianNumbers('G',[qn1,qn2,qn1,qn2],[1,2,3,4],qnscounts)
    @test filter((qn1,qn2),qns)==AbelianNumbers('G',[qn1,qn2,qn1,qn2],[1,2,3,4],qnscounts)
    @test filter(qn2,qns)==AbelianNumbers('G',[qn2,qn2],[2,4],qnscounts)
end

@testset "ukron" begin
    qn1,qn2=CNZ4(1.0,1.0),CNZ4(2.0,3.0)
    qns1=AbelianNumbers('U',[qn2,qn1],[4,5],qnscounts)
    qns2=AbelianNumbers('U',[qn1,qn2],[2,3],qnscounts)
    qns,records=ukron(qns1,qns2,signs=(+1,-1))
    @test qns==AbelianNumbers('C',[qn1-qn2,CNZ4|>zero,qn2-qn1],[15,22,8],qnscounts)
    @test records==Dict(
        (qn1-qn2) => Dict((qn1,qn2)=>1:15),
        CNZ4|>zero => Dict((qn2,qn2)=>1:12,(qn1,qn1)=>13:22),
        (qn2-qn1) => Dict((qn2,qn1)=>1:8)
    )
end

@testset "expand" begin
    qn=CNZ4(3.0,2.0)
    qns=AbelianNumbers(qn,4)
    @test expand(qns,qnsindices)==[1,1,1,1]
    @test expand(qns,qnscontents)==[qn,qn,qn,qn]
end

@testset "decompose" begin
    qn1,qn2=CNZ4(0.0,0.0),CNZ4(1.0,1.0)
    qns=AbelianNumbers('U',[qn1,qn2],[0,1,2],qnsindptr)
    target=CNZ4(2.0,2.0)
    result=Set(((1,1,2,2),(1,2,1,2),(1,2,2,1),(2,1,1,2),(2,1,2,1),(2,2,1,1)))
    @test ⊆(Set(decompose((qns,qns,qns,qns),target,(1,1,1,1),qnsbruteforce,nmax=10)),result)
    @test ⊆(Set(decompose((qns,qns,qns,qns),target,(1,1,1,1),qnsmontecarlo,nmax=10)),result)
end

@testset "permute" begin
    qn1,qn2,qn3=CNZ4(1.0,2.0),CNZ4(3.0,0.0),CNZ4(4.0,1.0)
    qns=AbelianNumbers('G',[qn1,qn2,qn3],[2,3,4],qnscounts)
    @test permute(qns,[3,2,1],qnscompression)==AbelianNumbers('G',[qn3,qn2,qn1],[4,3,2],qnscounts)
    @test permute(qns,[4,6,9,8],qnsexpansion)==AbelianNumbers('G',[qn2,qn3,qn3,qn3],[1,1,1,1],qnscounts)
end


text:
#--------------------------------------
## LOAD PACKAGES
using RegNeuralDE, StochasticDiffEq, Flux, DiffEqFlux, Tracker
using YAML, Dates, BSON, Random, Statistics, Printf
using CUDA
using RegNeuralDE: accuracy
using Flux: @functor
using Tracker: TrackedReal
import Base.show

# This code will not work on GPU. Please trick your system to not detect GPUs
# using `CUDA_VISIBLE_DEVICES=""`. Will get this fixed soon
CUDA.allowscalar(false)
#--------------------------------------

#--------------------------------------
## CONFIGURATION
config_file = joinpath(pwd(), "experiments", "configs", "mnist_nsde.yml")
config = YAML.load_file(config_file)

Random.seed!(config["seed"])

hparams = config["hyperparameters"]
BATCH_SIZE = hparams["batch_size"]
EPOCHS = hparams["epochs"]
REGULARIZE = hparams["regularize"]
REG_TYPE = hparams["type"]
identifier =
    REGULARIZE ? "$(string(now()))_$(REGULARIZE)_$(REG_TYPE)" : "$(string(now()))_vanilla"
EXPERIMENT_LOGDIR = joinpath(pwd(), "results", "mnist_nsde", identifier)
MODEL_WEIGHTS = joinpath(EXPERIMENT_LOGDIR, "weights.bson")
FILENAME = joinpath(EXPERIMENT_LOGDIR, "results.yml")

# Create a directory to store the results
isdir(EXPERIMENT_LOGDIR) || mkpath(EXPERIMENT_LOGDIR)
cp(config_file, joinpath(EXPERIMENT_LOGDIR, "config.yml"))
#--------------------------------------

#--------------------------------------
## SETUP THE MODELS + DATASET + TRAINING UTILS
# Get the dataset
train_dataloader, test_dataloader = load_mnist(BATCH_SIZE, x -> flatten(Float32.(x)))

agg = mean
if REG_TYPE == "error_est"
    λ₀ = 1.0f1
    λ₁ = 1.0f1
    save_func(u, t, integrator) = integrator.EEst * integrator.dt
    solver = SOSRI()
    global agg = mean
elseif REG_TYPE == "stiff_est"
    λ₀ = 0.1f0
    λ₁ = 0.1f0
    const stability_size =
        Tracker.TrackedReal(1 / Float32(StochasticDiffEq.alg_stability_size(SOSRI2())))
    function save_func(u, t, integrator)
        stiff_est = abs(integrator.eigen_est)
        return stability_size * ((iszero(stiff_est) || isnan(stiff_est)) ? 0 : stiff_est)
    end
    global agg = mean
    solver = AutoSOSRI2(SOSRI2())
else
    global agg = mean
    solver = SOSRI()
end
k = log(λ₀ / λ₁) / EPOCHS
# Exponential Decay
λ_func(t) = λ₀ * exp(-k * t)

nsde = ClassifierNSDE(
    Dense(784, 32) |> track,
    TrackedNeuralDSDE(
        Chain(Dense(32, 64, tanh), Dense(64, 32)) |> track,
        Dense(32, 32) |> track,
        [0.0f0, 1.0f0],
        REGULARIZE,
        solver,
        save_everystep = false,
        reltol = 1.4f-1,
        abstol = 1.4f-1,
        save_start = false,
    ),
    Dense(32, 10) |> track,
)
ps = Flux.trainable(nsde)

opt = Flux.Optimise.Optimiser(InvDecay(1.0e-5), ADAM(0.01))

function loss_function(
    x,
    y,
    model,
    p1,
    p2,
    p3;
    trajectories = 1,
    λ = 1.0f2,
    notrack = false,
)
    pred, _, _, sv = model(x, p1, p2, p3; trajectories = trajectories, func = save_func)
    cross_entropy = Flux.Losses.logitcrossentropy(pred, y)
    reg = REGULARIZE ? λ * agg(sv.saveval) : zero(eltype(pred))
    total_loss = cross_entropy + reg
    if !notrack
        ce_un = cross_entropy |> untrack
        reg_un = reg |> untrack
        total_loss_un = total_loss |> untrack
        logger(
            false,
            Dict(
                "Total Loss" => total_loss_un,
                "Cross Entropy Loss" => ce_un,
                "Regularization" => reg_un,
            ),
        )
    end
    return total_loss
end
#--------------------------------------

#--------------------------------------
## LOGGING UTILITIES
nfe1_counts = Vector{Float64}(undef, EPOCHS + 1)
nfe2_counts = Vector{Float64}(undef, EPOCHS + 1)
train_accuracies = Vector{Float64}(undef, EPOCHS + 1)
test_accuracies = Vector{Float64}(undef, EPOCHS + 1)
train_runtimes = Vector{Float64}(undef, EPOCHS + 1)
inference_runtimes = Vector{Float64}(undef, EPOCHS + 1)
train_runtimes[1] = 0

logger = table_logger(
    [
        "Epoch Number",
        "NFE1 Count",
        "NFE2 Count",
        "Train Accuracy",
        "Test Accuracy",
        "Train Runtime",
        "Inference Runtime",
    ],
    ["Total Loss", "Cross Entropy Loss", "Regularization"],
)
#--------------------------------------

#--------------------------------------
## RECORD DETAILS BEFORE TRAINING STARTS
dummy_data = train_dataloader.data[1][:, 1:BATCH_SIZE]
stime = time()
_, _nfe1, _nfe2, _ = nsde(dummy_data; func = save_func)
inference_runtimes[1] = time() - stime
train_runtimes[1] = 0.0
nfe1_counts[1] = _nfe1
nfe2_counts[1] = _nfe2
train_accuracies[1] = accuracy(nsde, train_dataloader; no_gpu = true, trajectories = 10)
test_accuracies[1] = accuracy(nsde, test_dataloader; no_gpu = true, trajectories = 10)

logger(
    false,
    Dict(),
    0.0,
    nfe1_counts[1],
    nfe2_counts[1],
    train_accuracies[1],
    test_accuracies[1],
    train_runtimes[1],
    inference_runtimes[1],
)
#--------------------------------------

#--------------------------------------
## WARMSTART THE GRADIENT COMPUTATION
y_ = zeros(Float32, 10, 1)
y_[1, :] .= 1.0
_ = Tracker.gradient(
    (p1, p2, p3) ->
        loss_function(rand(Float32, 784, 1), y_, nsde, p1, p2, p3; notrack = true),
    ps...,
)
#--------------------------------------

#--------------------------------------
## TRAINING
for epoch = 1:EPOCHS
    λ = λ_func(epoch - 1)
    timing = 0

    for (i, (x_, y_)) in enumerate(train_dataloader)
        x = x_
        y = y_

        start_time = time()
        gs = Tracker.gradient(
            (p1, p2, p3) -> loss_function(x, y, nsde, p1, p2, p3; λ = λ),
            ps...,
        )
        update_parameters!(ps, gs, opt)
        timing += time() - start_time

        x = y = nothing
        GC.gc(true)
    end

    # Record the time per epoch
    train_runtimes[epoch+1] = timing

    # Record the NFE count
    start_time = time()
    _, nfe1, nfe2, _ = nsde(dummy_data)
    inference_runtimes[epoch+1] = time() - start_time
    nfe1_counts[epoch+1] = nfe1
    nfe2_counts[epoch+1] = nfe2

    # Test and Train Accuracy
    train_accuracies[epoch+1] =
        accuracy(nsde, train_dataloader; no_gpu = true, trajectories = 10)
    test_accuracies[epoch+1] =
        accuracy(nsde, test_dataloader; no_gpu = true, trajectories = 10)

    logger(
        false,
        Dict(),
        epoch,
        nfe1_counts[epoch+1],
        nfe2_counts[epoch+1],
        train_accuracies[epoch+1],
        test_accuracies[epoch+1],
        train_runtimes[epoch+1],
        inference_runtimes[epoch+1],
    )
end
logger(true, Dict())
#--------------------------------------

#--------------------------------------
## STORE THE RESULTS
results = Dict(
    :nfe1_counts => nfe1_counts,
    :nfe2_counts => nfe2_counts,
    :train_accuracies => train_accuracies,
    :test_accuracies => test_accuracies,
    :train_runtimes => train_runtimes,
    :inference_runtimes => inference_runtimes,
)

weights = Flux.params(nsde) .|> cpu .|> untrack
BSON.@save MODEL_WEIGHTS weights

YAML.write_file(FILENAME, results)
#--------------------------------------


text:
#=
= Simulation state structures and methods.
=#

#= Misc. ======================================================================#

mutable struct Counts
    lprops::Int
    hprops::Int
    laccepts::Int
    haccepts::Int
    Counts() = new(0,0,0,0)
end

struct Constants
    k::Int
    R::Float64
    α::Float64
    β::Float64
    ll::Float64 # leapfrog step sizes for
    lh::Float64 # location and height moves.
    Ll::Int     # leapfrog step counts for
    Lh::Int     # location and height moves.
end

#= State ======================================================================#

struct State
    locs::Vector{Float64}
    hts::Vector{Float64}
end

"""
Initialises a new state with evenly spaced step locations and average rates.
"""
function State(Y::Vector{Float64}, k::Int, R::Float64)
    S = State(collect(linspace(0.,R,k+2)), Vector{Float64}(k+1))
    s,y = 0., 1 # step location and data index.
    for i = 1:k+1
        sn = S.locs[i+1] # next step location.
        n = searchsortedfirst(@view(Y[y:end]), sn) - 1 # step data count.
        S.hts[i] = n/(sn-s)
        s,y = sn, y+n
    end
    return S
end

function string(S::State)
    io = IOBuffer()
    for s in S.locs
        write(io, string(s,", "))
    end
    for h in @view(S.hts[1:end-1])
        write(io, string(h,", "))
    end
    write(io, string(S.hts[end]))
    return String(take!(io))
end

#= Extended state =============================================================#

struct Auxiliary
    locs::Vector{Float64}
    hts::Vector{Float64}
    widths::Vector{Float64}
    counts::Vector{Int}
    tlocs::Vector{Float64} # t = logit(s/R).
    lhts::Vector{Float64}  # j = ln(h).
    ∇tlocs::Vector{Float64}
    ∇lhts::Vector{Float64}
    ptlocs::Vector{Float64}
    plhts::Vector{Float64}
    ix::Vector{Int} # index vector used for sorting.
    Auxiliary(k) = new(
        Vector{Float64}(k+2),
        Vector{Float64}(k+1),
        Vector{Float64}(k+1),
        Vector{Int}(k+1),
        Vector{Float64}(k),
        Vector{Float64}(k+1),
        Vector{Float64}(k),
        Vector{Float64}(k+1),
        Vector{Float64}(k),
        Vector{Float64}(k+1),
        collect(1:k)
    )
end

max_stars_repo_name: kevdur/rjcp

max_issues_repo_name: kevdur/rjcp

max_forks_repo_name: kevdur/rjcp


text:
using Plots

f(x, y) = (1/2pi)*exp(-(x^2+y^2)/2)

function main()
    plotlyjs()
    xs = -5.:0.01:5.
    ys = -1.:0.01:1.
    zs = [f(x, y) for y in ys, x in xs]
    plot(xs, ys, zs, linetype=:surface)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
## This code implements the algorithm in
## Rudolfo Gorenflo, Joulia Loutchko and Yuri Loutchko, *Computation of the Mittag-Leffler function and its derivative*,  Fract. Calc. Appl. Anal, **(2002)**
## Changes to the algorithm are:
## * added some special cases.

if VERSION >= v"0.5-"
    import QuadGK: quadgk
end

#using Cubature

macro br(n)
#    esc(:(println(STDERR, "branch ", $n)))
    nothing
end

function P(α,β,ϵ,ϕ,z)
    ω = ϕ * (1+(1-β)/α) + ϵ^(1/α) * sin(ϕ/α)
    res = (1/(2*α*pi)) * ϵ^(1+(1-β)/α)*exp(ϵ^(1/α) * cos(ϕ/α)) * (cos(ω) + im * sin(ω))/(ϵ*exp(im*ϕ)-z)
end

ourquadgk(f,a,b) = quadgk(f,a,b; order=7)[1]
#ourquadgk(f,a,b) = hquadrature(f,a,b)[1]
#ourquadgk(f,a,b) = pquadrature(f,a,b)[1]

Pint(α,β,z) = Pint(α,β,1,z)
Pint(α,β,ϵ,z) = ourquadgk( ϕ -> P(α,β,ϵ,ϕ,z), -α*pi, α*pi)

function K(α,β,χ,z)
    den = (χ^2-2*χ*z*cos(α*pi)+z^2)
    res = (1/(α*pi)) * χ^((1-β)/α) * exp(-χ^(1/α))*( χ * sin(pi*(1-β)) - z * sin(pi*(1-β+α)))/den
end

Kint(α,β,χ0,z) = Kint(α,β,0,χ0,z)

function Kint(α,β,a,χ0,z)
    ourquadgk(χ -> K(α,β,χ,z), a, χ0)
end

mpow(x::Complex,y) = x^y
mpow(x::Real,y) = x >= 0 ? x^y : Complex(x,0)^y

function mittleffsum(α,β,z)
    @br 1
    k0 = floor(Int,α) + 1
    s = zero(z)
    for k=0:(k0-1)
        s += mittleff(α/k0,β,mpow(z,(1/k0))*exp(2*pi*im*k/k0))
    end
    s / k0
end

function mittleffsum2(α,β,z,ρ)
    @br 2
    k0 = max(ceil(Int,(1-β)/α), ceil(Int, log(ρ*(1-abs(z)))/log(abs(z))))
    s = zero(z)
    for k=0:k0
        s += z^k/gamma(β+α*k)
    end
    s
end

function sum2(α,β,z,k0)
    s = zero(z)
    for k=1:k0
        arg = β-α * k
        if !( round(arg) == arg && arg < 0)
            s += mpow(z,-k)/gamma(arg)
        end
    end
    s
end

function choosesum(α,β,z,ρ)
    k0 = floor(Int, -log(ρ)/log(abs(z)))
    if abs(angle(z)) < pi*α/4 + 1//2 * min(pi,pi*α)
        @br 3
        return 1/α * z^((1-β)/α) * exp(z^(1/α)) - sum2(α,β,z,k0)
    else
        @br 4
        return - sum2(α,β,z,k0)
    end
end

function mittleffints(α,β,z,ρ)
    az = abs(z)
    ab = abs(β)
    χ0 = β >= 0 ?
      max(1,2*az,(-log(pi*ρ/6))^α) :
      max((ab+1)^α, 2*az,(-2*log( pi*ρ/(6*(ab+2)*(2*ab)^ab)))^α)
    aaz = abs(angle(z))
    if aaz > α * pi
        if β <= 1
            @br 5
            return Kint(α,β,χ0,z)
        else
            @br 6
            return Kint(α,β,1,χ0,z) + Pint(α,β,z)
        end
    elseif aaz < pi*α
        if β <= 1
            @br 7
            return Kint(α,β,χ0,z) + (1/α)*z^((1-β)/α) * exp(z^(1/α))
        else
            @br 8
            return Kint(α,β,az/2,χ0,z) + Pint(α,β,(az)/2,z) + (1/α)*z^((1-β)/α) * exp(z^(1/α))
        end
    else
        @br 9
        return Kint(α,β,(az+1)/2,χ0,z) + Pint(α,β,(az+1)/2,z)
    end
end

"""
    mittlefferr(α,z,ρ)

Compute mittlefferr(α,1,z,ρ).
"""
mittlefferr(α,z,ρ) = mittlefferr(α,1,z,ρ)

"""
    mittlefferr(α,β,z,ρ)

Compute the Mittag-Leffler function at `z` for parameters `α,β` with
accuracy `ρ`.
"""
function mittlefferr(α,β,z,ρ)
    ρ > 0 || throw(DomainError())
    _mittlefferr(α,β,z,ρ)
end

_mittlefferr(α::Real,β::Real,z::Real,ρ::Real) = real(_mittleff(α,β,z,ρ))
_mittlefferr(α::Real,β::Real,z::Complex,ρ::Real) = _mittleff(α,β,z,ρ)

# The second definition would work for both complex and real
myeps(x) = x |> one |> float |> eps
myeps(x::Complex) =  x |> real |> myeps

"""
    mittleff(α,β,z)

Compute the Mittag-Leffler function at `z` for parameters `α,β`.
"""
mittleff(α,β,z) = _mittlefferr(α,β,z,myeps(z))

"""
    mittleff(α,z)

Compute `mittleff(α,1,z)`.
"""
mittleff(α,z) = _mittlefferr(α,1,z,myeps(z))

function _mittleff(α,β,z,ρ)
#    if β == 1
        # if α == 1/2   # disable this. There is never an error here. But, this triggers a mysterious, untraceable bug in quadgk
        #     res = try
        #         exp(z^2)*erfc(-z)
        #     catch
        #         error("Failed in exp. erfc")
        #     end
        #     return res
        # end
        # Disable these, because we would have to take care of domain errors, etc.
        # α == 0 && return 1/(1-z)
        # α == 1 && return exp(z)
        # α == 2 && return cosh(sqrt(z))
        # α == 3 && return (1//3)*(exp(z^(1//3)) + 2*exp(-z^(1//3)/2) * cos(sqrt(convert(typeof(z),3))/2 * z^(1//3)))
        # α == 4 && return (1//2)*(cosh(z^(1//4)) + cos(z^(1//4)))
#    end
    z == 0 && return 1/gamma(β)
    α == 1 && β == 1 && return(exp(z))
    α < 0  && throw(DomainError())
    az = abs(z)
    1 < α && return mittleffsum(α,β,z)
    az < 1 && return mittleffsum2(α,β,z,ρ)
    az > floor(10+5*α) && return choosesum(α,β,z,ρ)
    mittleffints(α,β,z,ρ)
end

_mittleff(α,β,z) = mittleff(α,β,z,myeps(z))

"""
    mittleffderiv(α,β,z)

Compute the derivative of the Mittag-Leffler function at `z` for parameters `α,β`.
"""
function mittleffderiv(α, β, z)
    #derivative of Mittag Leffler function WRT to main argument Z.
    #take q = 0.5. Paper requires |z| <= q < 1.
    q = 1//2

    #case 1, small z
    if abs(z) <= q

        ω = α + β - 3//2
        D = α^2 - 4*α*β + 6*α + 1

        #k1
        if α>1
            k₁ = ((2-α-β)/(α-1)) + 1
        elseif α>0 && α<=1 && D<=0
            k₁ = ((3-α-β)/α) + 1
        else
            k₁ = maximum([((3-α-β)/α) + 1, ((1-2*ω*α+sqrt(D))/(2*(α^2)))+1])
        end

        k₀ = maximum([k₁, log(myeps(z)*(1-abs(z)))/log(abs(z))])
        k₀ = ceil(Int,k₀) #take ceiling (not specified in paper whether floor or ceiling)

        out = zero(z)
        for k in 0:k₀
            out = out + ((k+1)*z^k)/(gamma(α+β+α*k))
        end

    #case 2, larger z
    else

        out = (mittleff(α,β-1,z) - (β-1)*mittleff(α,β,z))/(α*z)

    end

    return out
end

"""
    mittleffderiv(α,z)

Compute mittleffderiv(α,1,z)
"""
mittleffderiv(α, z) = mittleffderiv(α,1,z)


text:
import JSON2

function validate_serialization(sys::IS.SystemData)
    #path, io = mktemp()
    # For some reason files aren't getting deleted when written to /tmp. Using current dir.
    path = "test_system_serialization.json"
    io = open(path, "w")
    @info "Serializing to $path"

    try
        IS.prepare_for_serialization!(sys, path)
        IS.to_json(io, sys)
    catch
        close(io)
        rm(path)
        rethrow()
    end
    close(io)

    ts_file = nothing
    try
        ts_file = open(path) do file
            JSON2.read(file).time_series_storage_file
        end
        sys2 = IS.SystemData(path)
        return IS.compare_values(sys, sys2)
    finally
        @debug "delete temp file" path
        rm(path)
        rm(ts_file)
    end
end

@testset "Test JSON serialization of system data" begin
    for in_memory in (true, false)
        sys = create_system_data(; with_forecasts=true, time_series_in_memory=in_memory)
        @test validate_serialization(sys)
        text = JSON2.write(sys)
        @test length(text) > 0
    end
end

@testset "Test prepare_for_serialization" begin
    sys = create_system_data(; with_forecasts=true)
    IS.prepare_for_serialization!(sys, joinpath("dir1", "dir2", "sys.json"))
    @test sys.time_series_storage_file == joinpath("dir1", "dir2",
                                                   "sys_" * IS.TIME_SERIES_STORAGE_FILE)
end


text:
__precompile__()

"""
*MultipleTesting* package

* Adjusting p-values for multiple testing
* Estimating the fraction π0 of tests under the null hypothesis
"""
module MultipleTesting

using StatsBase
import StatsBase: fit

using Distributions

export
    storey_pi0,
    bootstrap_pi0,
    lsl_pi0,
    twostep_pi0,
    rightboundary_pi0,
    flat_grenander_pi0,
    adjust,
    PValueAdjustmentMethod,
    Bonferroni,
    BenjaminiHochberg,
    BenjaminiHochbergAdaptive,
    BenjaminiYekutieli,
    BenjaminiLiu,
    Hochberg,
    Holm,
    Hommel,
    Sidak,
    ForwardStop,
    bonferroni,
    benjamini_hochberg,
    benjamini_yekutieli,
    benjamini_liu,
    holm,
    hommel,
    hochberg,
    sidak,
    forwardstop,
    qValues,
    estimate_pi0,
    Pi0Estimator,
    Storey,
    StoreyBootstrap,
    LeastSlope,
    Oracle,
    TwoStep,
    RightBoundary,
    CensoredBUM,
    CensoredBUMFit,
    BUM,
    BUMFit,
    FlatGrenander,
    isin,
    fit,
    BetaUniformMixtureModel

include("types.jl")
include("utils.jl")
include("pval-adjustment.jl")
include("pi0-estimators.jl")
include("qvalue.jl")
include("model.jl")

end


text:
using GeoData, Test, Statistics, Dates, Plots, DiskArrays, RasterDataSources
using GeoData.LookupArrays, GeoData.Dimensions
import ArchGDAL, NCDatasets
using GeoData: FileArray, GDALfile

include(joinpath(dirname(pathof(GeoData)), "../test/test_utils.jl"))
url = "https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif"
gdalpath = maybedownload(url)

@testset "array" begin

    @time gdalarray = GeoArray(gdalpath; mappedcrs=EPSG(4326), name=:test)

    @testset "load from url" begin
        A = GeoArray("/vsicurl/" * url)
        B = GeoArray(url)
        @test read(A) == read(B) == gdalarray
    end

    @testset "open" begin
        @test open(A -> A[Y=1], gdalarray) == gdalarray[:, 1, :]
        tempfile = tempname() * ".tif"
        cp(gdalpath, tempfile)
        gdalwritearray = GeoArray(tempfile)
        @test_throws ArchGDAL.GDAL.GDALError gdalwritearray .*= UInt8(2)
        open(gdalwritearray; write=true) do A
            A .*= UInt8(2)
            nothing
        end
        @test GeoArray(tempfile) == gdalarray .* UInt8(2)
    end

    @testset "read" begin
        @time A = read(gdalarray);
        @test A isa GeoArray
        @test parent(A) isa Array
        A2 = zero(A)
        @time read!(gdalarray, A2);
        A3 = zero(A)
        @time read!(gdalpath, A3);
        @test A == A2 == A3
    end

    @testset "view" begin
        A = view(gdalarray, 1:10, 1:10, 1)
        @test A isa GeoArray
        @test parent(A) isa DiskArrays.SubDiskArray
        @test parent(parent(A)) isa GeoData.FileArray
    end

    @testset "array properties" begin
        @test size(gdalarray) == (514, 515, 1)
        @test gdalarray isa GeoArray{UInt8,3}
    end

    @testset "dimensions" begin
        @test length(dims(gdalarray, X)) == 514
        @test ndims(gdalarray) == 3
        @test dims(gdalarray) isa Tuple{<:X,<:Y,<:Band}
        @test lookup(gdalarray, Band) isa DimensionalData.Categorical;
        # @test span(gdalarray, (Y, X)) ==
            # (Regular(-60.02213698319351), Regular(60.02213698319374))
        @test sampling(gdalarray, (Y, X)) ==
            (Intervals(Start()), Intervals(Start()))
        @test refdims(gdalarray) == ()
        # Bounds calculated in python using rasterio
        @test all(bounds(gdalarray, Y) .≈ (4224973.143255847, 4255884.5438021915))
        @test all(bounds(gdalarray, X) .≈ (-28493.166784412522, 2358.211624949061))
    end

    @testset "other fields" begin
        # This file has an inorrect missing value
        @test missingval(gdalarray) == nothing
        @test metadata(gdalarray) isa Metadata{GDALfile}
        @test basename(metadata(gdalarray).val[:filepath]) == "cea.tif"
        @test name(gdalarray) == :test
        @test label(gdalarray) == "test"
        @test units(gdalarray) == nothing
        @test mappedcrs(dims(gdalarray, Y)) == EPSG(4326)
        @test mappedcrs(dims(gdalarray, X)) == EPSG(4326)
        @test mappedcrs(gdalarray) == EPSG(4326)
        @test mappedcrs(gdalarray[Y(1)]) == EPSG(4326)
        @test_throws ErrorException mappedcrs(gdalarray[Y(1), X(1)])
        wkt = WellKnownText(GeoFormatTypes.CRS(),
          "PROJCS[\"unnamed\",GEOGCS[\"NAD27\",DATUM[\"North_American_Datum_1927\",SPHEROID[\"Clarke 1866\",6378206.4,294.978698213898,AUTHORITY[\"EPSG\",\"7008\"]],AUTHORITY[\"EPSG\",\"6267\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4267\"]],PROJECTION[\"Cylindrical_Equal_Area\"],PARAMETER[\"standard_parallel_1\",33.75],PARAMETER[\"central_meridian\",-117.333333333333],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")
        @test crs(dims(gdalarray, Y)) == wkt
        @test crs(dims(gdalarray, X)) == wkt
        @test crs(gdalarray) == wkt
        @test crs(gdalarray[Y(1)]) == wkt
        @test_throws ErrorException crs(gdalarray[Y(1), X(1)])
    end

    @testset "indexing" begin
        @test gdalarray[Band(1)] isa GeoArray{UInt8,2}
        @test gdalarray[Y(1), Band(1)] isa GeoArray{UInt8,1}
        @test gdalarray[X(1), Band(1)] isa GeoArray{UInt8,1}
        @test gdalarray[X(1), Y(1), Band(1)] isa UInt8
        @test gdalarray[1, 1, 1] isa UInt8
        # Value indexed in python with rasterio
        @test gdalarray[Y(21), X(21), Band(1)] == 82
    end

    @testset "selectors" begin
        # TODO verify the value with R/gdal etc
        @test gdalarray[Y(Contains(33.8)), X(Contains(-117.5)), Band(1)] isa UInt8
        @test gdalarray[Y(Between(33.7, 33.9)), Band(1)] isa GeoArray
    end

    @testset "methods" begin
        @testset "mean" begin
            @test all(mean(gdalarray; dims=Y) .=== mean(parent(gdalarray); dims=2))
        end

        @testset "trim, crop, extend" begin
            a = replace_missing(gdalarray, zero(eltype(gdalarray)))
            a[X(1:100)] .= missingval(a)
            trimmed = trim(a)
            @test size(trimmed) == (414, 514, 1)
            cropped = GeoData.crop(a; to=trimmed)
            @test size(cropped) == (414, 514, 1)
            @test all(collect(cropped .=== trimmed))
            extended = extend(cropped; to=a)
            @test all(collect(extended .=== a))
        end

        @testset "mask and mask! to disk" begin
            msk = replace_missing(gdalarray, missing)
            msk[X(1:100), Y([1, 5, 95])] .= missingval(msk)
            @test !any(gdalarray[X(1:100)] .=== missingval(msk))
            masked = mask(gdalarray; to=msk)
            @test all(masked[X(1:100), Y([1, 5, 95])] .=== missingval(msk))
            tempfile = tempname() * ".tif"
            cp(gdalpath, tempfile)
            @test !all(GeoArray(tempfile)[X(1:100), Y([1, 5, 95])] .=== 0x00)
            open(GeoArray(tempfile); write=true) do A
                mask!(A; to=msk, missingval=0x00)
            end
            @test all(GeoArray(tempfile)[X(1:100), Y([1, 5, 95])] .=== 0x00)
            rm(tempfile)
        end

        @testset "polygon mask" begin
            A = read(GeoArray(gdalpath; name=:test))
            ds = map(dims(A)) do d
                DimensionalData.maybeshiftlocus(Center(), d)
            end
            A = set(A, ds)
            A = set(A, X => Points, Y=> Points)
            polymask = [[-20000, 4.23e6],
                        [-20000, 4.245e6],
                        [0.0, 4.245e6],
                        [0.0, 4.23e6],
                        [-20000, 4.23e6]]
            rastermask = replace_missing(copy(A), missing)
            # mask or `Between` is not exactly accurate yet...
            section = X(Between(-20000.0, 0.0)), Y(Between(4.23005e6, 4.24495e6)), Band(1)
            rastermask .= missing
            rastermask[section...] .= A[section...]
            pmasked = mask(A; to=polymask) |> read
            rmasked = mask(A; to=rastermask) |> read
            @test all(rmasked .=== pmasked)
        end

        @testset "classify! to disk" begin
            tempfile = tempname() * ".tif"
            cp(gdalpath, tempfile)
            open(GeoArray(tempfile); write=true) do A
                classify!(A, [0x01 0xcf 0x00; 0xd0 0xff 0xff])
            end
            @test count(==(0x00), GeoArray(tempfile)) + count(==(0xff), GeoArray(tempfile)) == length(GeoArray(tempfile))
        end

        @testset "aggregate" begin
            ag = aggregate(mean, gdalarray, 4)
            @test ag == aggregate(mean, gdalarray, (X(4), Y(4), Band(1)))
            tempfile = tempname() * ".tif"
            write(tempfile, ag)
            open(GeoArray(tempfile); write=true) do dst
                aggregate!(mean, dst, gdalarray, 4)
            end
            @test GeoArray(tempfile) == ag
        end

        @testset "mosaic" begin
            @time gdalarray = GeoArray(gdalpath; name=:test)
            A1 = gdalarray[X(1:300), Y(1:200)]
            A2 = gdalarray[X(57:500), Y(101:301)]
            tempfile = tempname() * ".tif"
            Afile = mosaic(first, A1, A2; missingval=0x00, atol=1e-8, filename=tempfile)
            Amem = mosaic(first, A1, A2; missingval=0x00, atol=1e-8)
            Atest = gdalarray[X(1:500), Y(1:301)]
            Atest[X(1:56), Y(201:301)] .= 0x00
            Atest[X(301:500), Y(1:100)] .= 0x00
            @test all(Atest .=== Amem .=== Afile)
        end

        @testset "rasterize" begin
            A = rebuild(read(gdalarray); missingval=0x00)
            # We round-trip rasterise the Tables.jl form of A
            R = rasterize(A; to=A)
            @test all(A .=== R .== gdalarray)
            R = rasterize(A; to=A, name=:test)
            @test all(A .=== R .== gdalarray)
            B = rebuild(read(gdalarray) .= 0x00; missingval=0x00)
            rasterize!(B, read(gdalarray))
            @test all(B .=== gdalarray |> collect)
        end

        @testset "chunk_series" begin
            @test GeoData.chunk_series(gdalarray) isa GeoSeries
            @test size(GeoData.chunk_series(gdalarray)) == (1, 1, 1)
        end

    end

    @testset "conversion to GeoArray" begin
        geoA = gdalarray[X(1:50), Y(1:1), Band(1)]
        @test size(geoA) == (50, 1)
        @test eltype(geoA) <: UInt8
        @time geoA isa GeoArray{UInt8,1}
        @test dims(geoA) isa Tuple{<:X,Y}
        @test refdims(geoA) isa Tuple{<:Band}
        @test metadata(geoA) == metadata(gdalarray)
        @test missingval(geoA) == nothing
        @test name(geoA) == :test
    end

    @testset "write" begin
        gdalarray = GeoArray(gdalpath; mappedcrs=EPSG(4326), name=:test);

        @testset "2d" begin
            geoA = view(gdalarray, Band(1))
            filename = tempname() * ".asc"
            @time write(filename, geoA)
            saved1 = GeoArray(filename; mappedcrs=EPSG(4326))[Band(1)];
            eltype(saved1)
            eltype(geoA)
            @test all(saved1 .== geoA)
            # @test typeof(saved1) == typeof(geoA)
            @test val(dims(saved1, X)) ≈ val(dims(geoA, X))
            @test val(dims(saved1, Y)) ≈ val(dims(geoA, Y))
            @test all(metadata.(dims(saved1)) .== metadata.(dims(geoA)))
            @test metadata(dims(saved1)[1]) == metadata(dims(geoA)[1])
            @test missingval(saved1) === missingval(geoA)
            @test refdims(saved1) == refdims(geoA)
        end

        @testset "3d, with subsetting" begin
            geoA2 = gdalarray[Y(Between(33.7, 33.9)), X(Between(-117.6, -117.4))]
            filename2 = tempname() * ".tif"
            write(filename2, geoA2)
            saved2 = read(GeoArray(filename2; name=:test, mappedcrs=EPSG(4326)))
            @test size(saved2) == size(geoA2) == length.(dims(saved2)) == length.(dims(geoA2))
            @test refdims(saved2) == refdims(geoA2)
            #TODO test a file with more metadata
            @test val(metadata(saved2))[:filepath] == filename2
            @test missingval(saved2) === missingval(geoA2)
            @test GeoData.name(saved2) == GeoData.name(geoA2)
            @test step(lookup(dims(saved2, Y))) ≈ step(lookup(dims(geoA2, Y)))
            @test step(lookup(dims(saved2, X))) ≈ step(lookup(dims(geoA2, X)))
            @test typeof(dims(saved2)) == typeof(dims(geoA2))
            @test all(val(dims(saved2, Band)) .≈ val(dims(geoA2, Band)))
            @test all(val(dims(saved2, X)) .≈ val(dims(geoA2, X)))
            @test all(val(dims(saved2, Y)) .≈ val(dims(geoA2, Y)))
            @test all(metadata.(dims(saved2)) .== metadata.(dims(geoA2)))
            @test parent(saved2) == parent(geoA2)
            @test typeof(saved2) == typeof(geoA2)
            filename3 = tempname() * ".tif"
            geoA3 = cat(gdalarray[Band(1)], gdalarray[Band(1)], gdalarray[Band(1)]; dims=Band(1:3))
            write(filename3, geoA3)
            saved3 = read(GeoArray(filename3; mappedcrs=EPSG(4326)))
            @test all(saved3 .== geoA3)
            @test val(dims(saved3, Band)) == 1:3
        end

        @testset "resave current" begin
            filename = tempname() * ".rst"
            write(filename, gdalarray)
            gdalarray2 = GeoArray(filename)
            write(gdalarray2)
            @test read(GeoArray(filename)) == read(gdalarray2)
        end

        @testset "to grd" begin
            write("testgrd.gri", gdalarray)
            grdarray = GeoArray("testgrd.gri")
            @test crs(grdarray) == convert(ProjString, crs(gdalarray))
            @test all(map((a, b) -> all(a .≈ b), bounds(grdarray), bounds(gdalarray)))
            @test index(grdarray, Y) ≈ index(gdalarray, Y)
            @test val(dims(grdarray, X)) ≈ val(dims(gdalarray, X))
            @test grdarray == gdalarray
        end

        @testset "from GeoArray" begin
            filename = tempname() * ".tiff"
            ga = GeoArray(rand(100, 200), (X, Y))
            write(filename, ga)
            saved = GeoArray(filename)
            @test all(reverse(saved[Band(1)]; dims=Y) .=== ga)
            @test saved[1, end, 1] == saved[At(1.0), At(1.0), At(1.0)]
            @test saved[100, 1, 1] == saved[At(100), At(200), At(1)]
            filename2 = tempname() * ".tif"
            ga2 = GeoArray(rand(100, 200), (X(Sampled(101:200)), Y(Sampled(1:200))))
            write(filename2, ga2)
            @test reverse(GeoArray(filename2)[Band(1)]; dims=Y) == ga2
        end

        @testset "to netcdf" begin
            filename2 = tempname() * ".nc"
            write(filename2, gdalarray[Band(1)])
            saved = GeoArray(filename2; crs=crs(gdalarray))
            @test size(saved) == size(gdalarray[Band(1)])
            @test saved ≈ gdalarray[Band(1)]
            clat, clon = DimensionalData.shiftlocus.(Ref(Center()), dims(gdalarray, (Y, X)))
            @test mappedindex(clat) ≈ mappedindex(saved, Y)
            @test mappedindex(clon) ≈ mappedindex(saved, X)
            @test all(mappedbounds(saved, X) .≈ mappedbounds(clon))
            @test all(mappedbounds(saved, Y) .≈ mappedbounds(clat))
            @test projectedindex(clon) ≈ projectedindex(saved, X)
            @test all(projectedbounds(clon) .≈ projectedbounds(saved, X))
            # reason lat crs conversion is less accrurate than lon TODO investigate further
            @test all(map((a, b) -> isapprox(a, b; rtol=1e-6),
                projectedindex(gdalarray, Y),
                projectedindex(DimensionalData.shiftlocus(Start(), dims(saved, Y)))
            ))
            @test all(map((a, b) -> isapprox(a, b; rtol=1e-6), projectedbounds(saved, Y),  projectedbounds(gdalarray, Y)))
        end

        @testset "write missing" begin
            A = replace_missing(gdalarray, missing)
            filename = tempname() * ".tif"
            write(filename, A)
            @test missingval(GeoArray(filename)) === typemin(UInt8)
            rm(filename)
        end

    end

    @testset "show" begin
        sh = sprint(show, MIME("text/plain"), gdalarray)
        # Test but don't lock this down too much
        @test occursin("GeoArray", sh)
        @test occursin("Y", sh)
        @test occursin("X", sh)
        @test occursin("Band", sh)
    end

    @testset "plot" begin # TODO write some tests for this
        gdalarray |> plot
        gdalarray[Y(1)] |> plot
    end

    @testset "nodatavalue type matches the array type" begin
        A = GeoArray(WorldClim{Climate}, :tavg; res="10m", month=1)
        @test typeof(missingval(A)) === eltype(A)
        @test missingval(A) === -3.4f38
    end

end

@testset "stack" begin
    @time gdalstack = GeoStack((a=gdalpath, b=gdalpath))

    @test length(gdalstack) == 2
    @test dims(gdalstack) isa Tuple{<:X,<:Y,<:Band}

    @testset "read" begin
        st = read(gdalstack)
        read!((a=gdalpath, b=gdalpath), st)
        @test st isa GeoStack
        @test st.data isa NamedTuple
        @test first(st.data) isa Array
        @test all(st[:a] .=== gdalstack[:a])
    end

    @testset "child array properties" begin
        @test size(gdalstack[:a]) == (514, 515, 1)
        @test gdalstack[:a] isa GeoArray{UInt8,3}
    end

    @testset "indexing" begin
        @test gdalstack[:a][Y(2:3), X(1), Band(1)] == [0x00, 0x6b]
        @test gdalstack[:a][Y(1), X(1), Band(1)] == 0x00
        @test gdalstack[:b, Band(1)] == gdalstack[:b][Band(1)]
        @test typeof(gdalstack[:b, Band(1)]) == typeof(gdalstack[:b][Band(1)])
        @test view(gdalstack, Y(2:3), X(1), Band(1))[:a] == [0x00, 0x6b]
    end

    @testset "methods" begin
        @testset "mean" begin
            means = map(A -> mean(parent(A); dims=2), gdalstack)
            @test map((a, b) -> all(a .== b), mean(gdalstack; dims=Y), means) |> all
        end
        @testset "trim, crop, extend" begin
            mv = zero(eltype(gdalstack[:a]))
            st = replace_missing(gdalstack, mv)
            st = map(A -> (view(A, X(1:100)) .= mv; A), st)
            trimmed = trim(st)
            @test size(trimmed) == (414, 514, 1)
            cropped = crop(st; to=trimmed)
            @test size(cropped) == (414, 514, 1)
            @test map((c, t) -> all(collect(c .=== t)), cropped, trimmed) |> all
            extended = extend(cropped; to=st)
            @test all(collect(extended .== st))
        end
        @testset "mask and mask!" begin
            st = read(gdalstack)
            msk = replace_missing(gdalstack[:a], missing)
            msk[X(1:100), Y([1, 5, 95])] .= missingval(msk)
            @test !any(st[:b][X(1:100)] .=== missingval(msk))
            masked = mask(st; to=msk)
            masked[:b][X(1:100), Y([1, 5, 95])]
            @test all(masked[:b][X(1:100), Y([1, 5, 95])] .=== missing)
            mask!(st; to=msk, missingval=0x00)
            @test all(st[:a][X(1:100), Y([1, 5, 95])] .=== 0x00)
            @test all(st[:b][X(1:100), Y([1, 5, 95])] .=== 0x00)
        end

        @testset "rasterize" begin
            st = map(A -> rebuild(A; missingval=0x00), gdalstack) |> read
            # We round-trip rasterise the Tables.jl form of A
            r_st = rasterize(read(gdalstack); to=st)
            @test all(map((a, b, c) -> all(a .=== b .=== c), st, r_st, gdalstack))
            r_st = rasterize(read(gdalstack); to=st, name=(:a, :b))
            @test all(map((a, b, c) -> all(a .=== b .=== c), st, r_st, gdalstack))
            st = map(A -> rebuild(A .* 0x00; missingval=0x00), gdalstack) |> read
            rasterize!(st, read(gdalstack))
            @test all(map((a, b) -> all(a .=== b), st, gdalstack))
        end
    end

    if VERSION > v"1.1-"
        @testset "copy" begin
            geoA = zero(GeoArray(gdalstack[:a]))
            copy!(geoA, gdalstack, :a)
            # First wrap with GeoArray() here or == loads from disk for each cell.
            # we need a general way of avoiding this in all disk-based sources
            @test all(geoA .== GeoArray(gdalstack[:a]))
        end
    end

    @testset "write" begin
        geoA = read(gdalstack[:a])
        @testset "write multiple files" begin
            filename = tempname() * ".tif"
            write(filename, gdalstack)
            base, ext = splitext(filename)
            filename_b = string(base, "_b", ext)
            saved = read(GeoArray(filename_b))
            @test all(saved .== geoA)
        end

        @testset "write multiple files with custom suffix" begin
            filename = tempname() * ".tif"
            write(filename, gdalstack; suffix=("_first", "_second"))
            base, ext = splitext(filename)
            filename_b = string(base, "_second", ext)
            saved = read(GeoArray(filename_b))
            @test all(saved .== geoA)
        end

        @testset "write netcdf" begin
            filename = tempname() * ".nc"
            write(filename, gdalstack);
            saved = GeoStack(filename);
            @test all(read(saved[:a]) .== geoA)
            rm(filename)
        end
    end

    @testset "show" begin
        sh = sprint(show, MIME("text/plain"), gdalstack)
        # Test but don't lock this down too much
        @test occursin("GeoStack", sh)
        @test occursin("Y", sh)
        @test occursin("X", sh)
        @test occursin("Band", sh)
        @test occursin(":a", sh)
        @test occursin(":b", sh)
    end

end

@testset "series" begin
    gdalser = GeoSeries([gdalpath, gdalpath], (Ti(),); mappedcrs=EPSG(4326), name=:test)
    @test read(gdalser[Ti(1)]) == read(GeoArray(gdalpath; mappedcrs=EPSG(4326), name=:test))
    @test read(gdalser[Ti(1)]) == read(GeoArray(gdalpath; mappedcrs=EPSG(4326), name=:test))

    gdalser = GeoSeries((a=[gdalpath, gdalpath], b=[gdalpath, gdalpath]), (Ti,))
    # Rebuild the ser by wrapping the disk data in Array.
    # `modify` forces `rebuild` on all containers as in-Memory variants
    modified_ser = modify(Array, gdalser)
    @test typeof(modified_ser) <: GeoSeries{<:GeoStack{<:NamedTuple{(:a,:b),<:Tuple{<:Array{UInt8,3},Vararg}}}}

    @testset "read" begin
        ser1 = read(gdalser)
        @test ser1 isa GeoSeries{<:GeoStack}
        @test ser1.data isa Vector{<:GeoStack}
        @test first(ser1.data[1].data) isa Array
        ser2 = modify(A -> A .* 0, ser1)
        ser3 = modify(A -> A .* 0, ser1)
        read!([(a=gdalpath, b=gdalpath), (a=gdalpath, b=gdalpath)], ser2)
        read!(ser1, ser3)
        @test map(ser1, ser2, ser3) do st1, st2, st3
            map(st1, st2, st3) do A1, A2, A3
                (A2 .=== A2 .=== A3) |> all
            end |> all
        end |> all
    end

    @testset "show" begin
        sh = sprint(show, MIME("text/plain"), gdalser)
        # Test but don't lock this down too much
        @test occursin("GeoSeries", sh)
        @test occursin("GeoStack", sh)
        @test occursin("Ti", sh)
    end
end


text:
n = rand(1:10)
@test matrixdepot("hilb", n) == matrixdepot("hilb", n, n)
@test_approx_eq matrixdepot("hilb" , Float16, n) matrixdepot("hilb", Float32, n)
x = Float64[1:n;]
y = x -1
@test_approx_eq matrixdepot("hilb", n) matrixdepot("cauchy", x, y)
println("'hilb' passed test...")


text:
# solver settings
using GeometricIntegrators.Config

set_config(:nls_atol_break, 1E3)
set_config(:nls_rtol_break, 1E3)
set_config(:nls_stol_break, 1E3)

set_config(:nls_atol, 1E-14)
set_config(:nls_rtol, 1E-14)
set_config(:nls_stol, 1E-14)

include("charged_particle_3d_tests.jl")
include("guiding_center_4d_tests.jl")
include("gyro_kinetics_4d_tests.jl")
include("pauli_particle_3d_tests.jl")


text:
using DemoLoopVectorization
using Test

@testset "DemoLoopVectorization.jl" begin
    # Write your tests here.
end


text:
"""
Submodule of `SirenSeq`.
Defines some shorcut names for the GM drum mapping.
"""
module GmDrums

const aBass 	= 35
const bass 		= 36
const aSnare  	= 38
const eSnare 	= 40

const flTom 	= 41
const fhTom 	= 43
const lTom		= 45
const mTom 		= 47
const mhTom 	= 48
const hTom 		= 50

const cHh 		= 42
const pHh 		= 44
const oHh 		= 46
const cCym1 	= 49
const rCym1 	= 51
const tamborine = 54
const cow 		= 56
const vibraslap = 58

const chCym 	= 52
const ride1 	= 53
const sCym 		= 55
const cCym2 	= 57
const ride2 	= 59

const lWood 	= 77
const hWood 	= 76

const cTri 		= 80
const oTri 		= 81

end


text:
export MetaAgentCBS_Solver

"""
    MetaAgentCBS_Solver

Path planner that employs Meta Agent Conflict-Based Search
"""
@with_kw struct MetaAgentCBS_Solver{L,C} <: AbstractCBSSolver
    low_level_planner::L    = VanillaAStar()
    logger::SolverLogger{C} = SolverLogger{cost_type(low_level_planner)}()
    beta::Int               = 1
end
MetaAgentCBS_Solver(planner) = MetaAgentCBS_Solver(low_level_planner=planner)

export MetaSolution

"""
    MetaSolution{S}

Wrapper for a LowLevelSolution that allows for keeping track of groups.
"""
struct MetaSolution{S}
    solution::S
    group_idxs::Vector{Vector{Int}}
end
MetaSolution(mapf::AbstractMAPF) = MetaSolution(
    get_initial_solution(mapf),
    map(i->[i],1:num_agents(mapf))
    )
for op in [
    :state_type,
    :action_type,
    :cost_type,
    :get_paths,
    :get_path_costs,
    :get_cost,
    :get_cost_model,
    :set_solution_path!,
    :set_path_cost!,
    :set_cost!,
    :is_consistent,
    :is_valid,
    ]
    @eval $op(s::MetaSolution,args...) = $op(s.solution,args...)
end
Base.copy(solution::L) where {L<:MetaSolution} = L(
    copy(solution.solution),
    deepcopy(solution.group_idxs)
    )

function initialize_root_node(solver::MetaAgentCBS_Solver,mapf::AbstractMAPF,solution=MetaSolution(mapf))
    initialize_root_node(mapf,solution)
end

"""
    `combine_agents(conflict_table, groups::Vector{Vector{Int}})`

    Helper for merging two (meta) agents into a meta-agent
"""
function combine_agents!(solver, node)
    groups = node.solution.group_idxs
    N = length(groups)
    conflict_counts = zeros(Int,N,N)
    for (i,idxs1) in enumerate(groups)
        for (j,idxs2) in enumerate(groups[i:end])
            conflict_counts[i,j] += count_conflicts(node.conflict_table,idxs1,idxs2)
        end
    end
    idx = argmax(conflict_counts).I
    i = minimum(idx)
    j = maximum(idx)
    if conflict_counts[i,j] > solver.beta
        # @show i,j,groups[i],groups[j],conflict_counts[i,j]
        @assert i != j
        @log_info(1,verbosity(solver),"Conflict Limit exceeded. Merging agents ", groups[i], " and ", groups[j],"\n")
        # groups = deepcopy(groups)
        groups[i] = [groups[i]..., groups[j]...]
        deleteat!(groups, j)
        # node.groups = groups
        # node.solution.group_idxs = groups
        return groups, i
    end
    return groups, -1
end

function get_group_index(groups, agent_idx)
    group_idx = -1
    for i in 1:length(groups)
        if agent_idx in groups[i]
            group_idx = i
            break
        end
    end
    return group_idx
end
function get_group_index(solution::MetaSolution, agent_idx)
    get_group_index(solution.group_idxs,agent_idx)
end

function cbs_bypass!(solver::MetaAgentCBS_Solver,mapf,node,priority_queue)
    groups, group_idx = combine_agents!(solver, node)
    @assert groups == node.solution.group_idxs
    if group_idx > 0 # New Meta Agent has been constructed
        new_node = initialize_child_search_node(solver, mapf, node)
        consistent_flag = low_level_search!(solver, mapf, new_node, [group_idx])
        if consistent_flag # is_valid(new_node.solution, mapf)
            # for agent_id in groups[group_idx]
            #     detect_conflicts!(new_node.conflict_table,new_node.solution,agent_id)
            # end
            detect_conflicts!(new_node,[group_idx])
            for i in groups[group_idx]
                for j in groups[group_idx]
                    @assert count_conflicts(new_node,i,j) == 0 "$i,$j,$(count_conflicts(node,i,j))"
                end
            end
            enqueue!(priority_queue, new_node => get_cost(new_node))
        end
        return true
    end
    return false
end

# function CRCBS.low_level_search!(solver::MetaAgentCBS_Solver,mapf,node,group_idxs=collect(1:length(node.solution.group_idxs)))
#     # group_idxs = union(map(i->get_group_index(node.solution,i), idxs))
#     low_level_search!(low_level(solver),mapf,node,group_idxs)
# end

function detect_conflicts!(conflict_table,solution::MetaSolution,
        group_idxs=collect(1:length(solution.group_idxs)),
        args...;
        kwargs...
    )
    for group_idx in group_idxs
        for idx in solution.group_idxs[group_idx]
            detect_conflicts!(conflict_table,solution.solution,idx,args...;kwargs...)
        end
    end
end

solve!(solver::MetaAgentCBS_Solver, mapf) = cbs!(solver,mapf)


max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia