text:
using .Networks

using Dates
using DataStructures

function compute_prod(schedule, ts, scenario)
    prod = 0.
    for (gen_id, _) in schedule.generator_schedules
        gen_prod = safeget_prod_value(schedule, gen_id, ts, scenario)
        prod += gen_prod
    end
    return prod
end

"""
Computes EOD (i.e. supply-demand) without specifically considering the capping or the LoL :
sum(prod_schedule(gen))) - sum(load_uncertainties(bus)))
"""
function compute_eod(uncertainties_at_ech,
                    schedule,
                    network,
                    ts, scenario)
    prod = compute_prod(schedule, ts, scenario)
    load = compute_load(uncertainties_at_ech, network, ts, scenario)
    return (prod - load)
end
function compute_eod(uncertainties,
                    schedule,
                    network,
                    ech, ts, scenario)
    return compute_eod(get_uncertainties(uncertainties, ech),
                        schedule,
                        network,
                        ts, scenario)
end

function compute_flow(branch_id::String,
                uncertainties_at_ech::UncertaintiesAtEch,
                schedule::Schedule,
                network::Networks.Network,
                ts, scenario)
    #NOTE : this is not the function called in the models
    flow = 0.
    for bus in Networks.get_buses(network)
        bus_id = Networks.get_id(bus)
        ptdf = Networks.safeget_ptdf_elt(network, branch_id, bus_id)

        load = get_uncertainties(uncertainties_at_ech, bus_id, ts, scenario)
        cut_load = get_loss_of_load(schedule, bus_id, ts, scenario)
        load = ismissing(cut_load) ? load : ( load - cut_load )

        prod = 0.
        for generator in Networks.get_generators(bus)
            gen_id = Networks.get_id(generator)
            prod_l = safeget_prod_value(schedule, gen_id, ts, scenario)
            capped_prod = get_capping(schedule, gen_id, ts, scenario)
            prod_l = ismissing(capped_prod) ? prod_l : ( prod_l - capped_prod )

            prod += prod_l
        end

        flow -= ptdf * load
        flow += ptdf * prod
    end
    return flow
end
function compute_flow(branch_id::String,
                uncertainties::Uncertainties,
                schedule::Schedule,
                network::Networks.Network,
                ech, ts, scenario)
    return compute_flow(branch_id,
                        get_uncertainties(uncertainties, ech),
                        schedule,
                        network,
                        ts, scenario)
end

function compute_flows(uncertainties_at_ech::UncertaintiesAtEch,
                        schedule::Schedule,
                        network::Networks.Network,
                        TS, scenarios)
    #branch, ts, s
    flows = SortedDict{Tuple{String, DateTime, String}, Float64}()
    for branch in Networks.get_branches(network)
        branch_id = Networks.get_id(branch)
        for ts in TS
            for scenario in scenarios
                val = compute_flow(branch_id,
                                uncertainties_at_ech,
                                schedule,
                                network,
                                ts, scenario)
                flows[branch_id, ts, scenario] = val
            end
        end
    end
    return flows
end

function compute_flows(context::PSCOPFContext,
                        schedule::Schedule)
    ech = schedule.decision_time
    return compute_flows(get_uncertainties(context, ech), schedule, get_network(context),
                        get_target_timepoints(context), get_scenarios(context))
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# these tests spawn subprocesses, so reset the current context to conserve memory
CUDA.async_alloc[] || CUDA.device_reset!()

@testset "stack traces at different debug levels" begin

script = """
    function kernel(arr, val)
        arr[1] = val
        return
    end

    cpu = zeros(Int)
    gpu = CuArray(cpu)
    @cuda kernel(gpu, 1.2)
    synchronize()

    # FIXME: on some platforms (Windows...), for some users, the exception flag change
    # doesn't immediately propagate to the host, and gets caught during finalization.
    # this looks like a driver bug, since we threadfence_system() after setting the flag.
    # https://stackoverflow.com/questions/16417346/cuda-pinned-memory-flushing-from-the-device
    sleep(1)
    synchronize()
"""

# NOTE: kernel exceptions aren't always caught on the CPU as a KernelException.
#       on older devices, we emit a `trap` which causes a CUDA error...
#

let (code, out, err) = julia_script(script, `-g0`)
    @test code == 1
    @test  occursin("ERROR: KernelException: exception thrown during kernel execution on device", err) ||
           occursin("ERROR: CUDA error: an illegal instruction was encountered", err) ||
           occursin("ERROR: CUDA error: unspecified launch failure", err)
    @test !occursin(r"ERROR: a \w+ was thrown during kernel execution", out)
    # NOTE: stdout sometimes contain a failure to free the CuArray with ILLEGAL_ACCESS
end

let (code, out, err) = julia_script(script, `-g1`)
    @test code == 1
    @test occursin("ERROR: KernelException: exception thrown during kernel execution on device", err) ||
          occursin("ERROR: CUDA error: an illegal instruction was encountered", err) ||
          occursin("ERROR: CUDA error: unspecified launch failure", err)
    @test occursin(r"ERROR: a \w+ was thrown during kernel execution", out)
    @test occursin("Run Julia on debug level 2 for device stack traces", out)
end

let (code, out, err) = julia_script(script, `-g2`)
    @test code == 1
    @test occursin("ERROR: KernelException: exception thrown during kernel execution on device", err) ||
          occursin("ERROR: CUDA error: an illegal instruction was encountered", err) ||
          occursin("ERROR: CUDA error: unspecified launch failure", err)
    @test occursin(r"ERROR: a \w+ was thrown during kernel execution", out)
    if VERSION < v"1.3.0-DEV.270"
        @test occursin("[1] Type at float.jl", out)
    else
        @test occursin("[1] Int64 at float.jl", out)
    end
    @test occursin("[4] kernel at none:5", out)
end

end

@testset "#329" begin

script = """
    @noinline foo(a, i) = a[1] = i
    bar(a) = (foo(a, 42); nothing)

    ptr = reinterpret(Core.LLVMPtr{Int,AS.Global}, C_NULL)
    arr = CuDeviceArray{Int,1,AS.Global}((0,), ptr)

    CUDA.@sync @cuda bar(arr)
"""

let (code, out, err) = julia_script(script, `-g2`)
    @test code == 1
    @test occursin("ERROR: KernelException: exception thrown during kernel execution on device", err) ||
          occursin("ERROR: CUDA error: an illegal instruction was encountered", err) ||
          occursin("ERROR: CUDA error: unspecified launch failure", err)
    @test occursin(r"ERROR: a \w+ was thrown during kernel execution", out)
    @test occursin("foo at none:4", out)
    @test occursin("bar at none:5", out)
end

end

ext: jl


text:
module NodeModule
import Base.show

    struct Node
        id::Int16
        x::Float32
        y::Float32
    end

    function show(io::IO, n::Node)
        print(string("Node:", n.id))
    end
end

max_stars_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_issues_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_forks_repo_name: zahidkizmaz/firefly-algorithm-tsp


text:
"""
````
function simplexgrid(coord::Array{Tc,2},
                     cellnodes::Array{Ti,2},
                     cellregions::Array{Ti,1},
                     bfacenodes::Array{Ti,2},
                     bfaceregions::Array{Ti,1}
                     ) where {Tc,Ti}
````

    Create simplex grid from five arrays.
"""
function simplexgrid(coord::Array{Tc,2},
                     cellnodes::Array{Ti,2},
                     cellregions::Array{Ti,1},
                     bfacenodes::Array{Ti,2},
                     bfaceregions::Array{Ti,1}
                     ) where {Tc,Ti}
    @assert size(coord,2)>0
    dim=size(coord,1)
    if dim==1
        eltype=Edge1D
        btype=Vertex0D
        csys=Cartesian1D
    elseif dim==2
        eltype=Triangle2D
        btype=Edge1D
        csys=Cartesian2D
    elseif dim==3
        eltype=Tetrahedron3D
        btype=Triangle2D
        csys=Cartesian3D
    end

    grid=ExtendableGrid{Tc,Ti}()
    grid[Coordinates]=coord
    grid[CellNodes]=cellnodes
    grid[CellRegions]=cellregions
    grid[CellGeometries]=VectorOfConstants(eltype,length(cellregions))
    grid[BFaceNodes]=bfacenodes
    grid[BFaceRegions]=bfaceregions
    grid[BFaceGeometries]=VectorOfConstants(btype,length(bfaceregions))
    grid[CoordinateSystem]=csys
    return grid
end

"""
$(TYPEDSIGNATURES)

(Try to) create a subdivision of interval (a,b) stored in the
returned array X such that
  - `X[1]==a, X[end]==b`
  - `(X[2]-X[1])<=ha+tol*(b-a)`
  - `(X[end]-X[end-1])<=hb+tol*(b-a)`
  - There is a number q such that  `X[i+1]-X[i] == q*(X[i]-X[i-1])`
  - X is the array with the minimal possible number of points with the above property

Caveat: the algorithm behind this is  well tested but unproven.

Returns an Array containing the points of the subdivision.
"""
function geomspace(a, b, ha, hb ; tol=1.0e-10)

    function _geomspace0(l,h0, hl, tol=1.0e-10)
        @assert (l>0.0)
        @assert (h0>0.0)
        @assert (hl>=h0)
        @assert((hl+h0)<l)

        #  We need to  adjust two things:

        # The sum of the geometric progression must
        # match the length of the interval, so lmismatch
        # should be zero:
        function lmismatch(q,k)
            return l - h0*(1-q^k)/(1-q)
        end

        # The claim from experimenral evidence (Wolfram)
        # is that, if written as a polynomial,
        # it has two real zeros: one and the value searched for which
        # is slightly larger than one. All other zeros are on one circle.

        # The size of the last interval should be close to
        # to hl, so hmismatch should be close to one and not larger than one
        function  hmismatch(q,k)
            return  h0*q^(k-1)/hl
        end

        # define initial number of intervals from
        # average of minmal and maximal h
        n=Int(ceil((2.0*l/(h0+hl))))

        if n==1
            n=2
        end

        # define initial q such that hmismatch is one.
        q=(hl/h0)^(1.0/(n-1.0))

        # Iteration until both mismatches are satisfactory
        # Outer loop runs until hmismatch is less than 1
        hmiss=10.0 # some initial value >1 just to run the loop at least once
        if abs(q-1.0)<tol
            hmiss=1.0
        end

        while  hmiss>1.0
            # increase number of intervals until
            # lmismatch becomes less than zero
            while  lmismatch(q,n)>0.0
                n+=1
            end

            # find initial interval for q containing
            # value with zero lmismatch
            ns=0
            nsmax=100

            while lmismatch(q,n)<0 &&  ns<nsmax
                q*=0.99
                ns+=1
            end

            xl=q
            xr=q/0.99
            @assert ns<nsmax

            # bisection to define q with zero lmismatch
            ns=0
            xm=0.5*(xl+xr)
            while !(xr≈xl) && ns<nsmax
                ns+=1
                mmm=lmismatch(xm,n)
                if mmm≈0.0
                    break
                elseif   lmismatch(xl,n)*mmm<0
                    xr=xm
                else
                    xl=xm
                end
                xm=0.5*(xl+xr)
            end
            q=xm
            @assert ns<nsmax
            hmiss=hmismatch(q,n)
            if hmiss>1.0+tol
                n=n+1
            end
        end
        #  printf("%d %g %g %g\n",n,q,lmismatch(q,n),hmismatch(q,n))

        X = Array{Float64,1}(undef,n+1)
        X[1]=0
        h=h0
        for i=1:n
            X[i+1]=X[i]+h
            h*=q
        end
        X[n+1]=l
        return X
    end

    # Map things to [0,b-a]
    @assert (ha>0.0)
    @assert (hb>0.0)
    @assert (a<b)

    tol=tol*(b-a)
    if ha<=hb
        X=_geomspace0(b-a,ha,hb,tol)
        X.+=a
    else
        X=-reverse(_geomspace0(b-a,hb,ha,tol))
        X.+=b
    end

    @assert (X[2]-X[1])<=ha+tol
    @assert (X[end]-X[end-1])<=hb+tol

    return X
end

collect_or_assign(X)=X
collect_or_assign(X::AbstractRange)=collect(X)
is_monotone(X)=all(X[2:end]-X[1:end-1].>0)

"""
$(TYPEDSIGNATURES)

Glue together two vectors a and b resulting in a vector c. They last element
of a shall be equal (up to tol) to the first element of b.
The result fulfills `length(c)=length(a)+length(b)-1`
"""
function glue(_a::AbstractVector, _b::AbstractVector; tol=1.0e-10)
    a=collect_or_assign(_a)
    b=collect_or_assign(_b)

    is_monotone(a) || error("a not monotone")
    is_monotone(b) || error("b not monotone")
    Tv=promote_type(eltype(a),eltype(b))
    na=length(a)
    nb=length(b)

    d=b[1]-a[na-1]
    @assert(d>0)
    d=b[1]-a[na]
    @assert(d>-tol)
    @assert(d<tol)

    c=Vector{Tv}(undef,na+nb-1)
    ic=0
    for ia=1:na
        ic+=1
        c[ic]=a[ia]
    end
    for ib=2:nb
        ic+=1
        c[ic]=b[ib]
    end
    return c
end

##########################################################
abstract type XCoordinates <: AbstractGridFloatArray1D end
abstract type YCoordinates <: AbstractGridFloatArray1D end
abstract type ZCoordinates <: AbstractGridFloatArray1D end

"""
$(TYPEDSIGNATURES)

Constructor for 1D grid.

Construct 1D grid from an array of node cordinates.
It creates two boundary regions with index 1 at the left end and
index 2 at the right end.

Primal grid holding unknowns: marked by `o`, dual
grid marking control volumes: marked by `|`.

```@raw html
 o-----o-----o-----o-----o-----o-----o-----o-----o
 |--|-----|-----|-----|-----|-----|-----|-----|--|
```
"""
function simplexgrid(_X::AbstractVector)
    X=collect_or_assign(_X)
    #    is_monotone(X) || error("X not monotone")
    coord=reshape(X,1,length(X))
    cellnodes=zeros(Int32,2,length(X)-1)
    cellregions=zeros(Int32,length(X)-1)
    for i=1:length(X)-1
        cellnodes[1,i]=i
        cellnodes[2,i]=i+1
        cellregions[i]=1
    end
    bfacenodes=Array{Int32}(undef,1,2)
    bfaceregions=zeros(Int32,2)
    bfacenodes[1,1]=1
    bfacenodes[1,2]=length(X)
    bfaceregions[1]=1
    bfaceregions[2]=2
    grid=simplexgrid(coord,
                     cellnodes,
                     cellregions,
                     bfacenodes,
                     bfaceregions)
    grid[XCoordinates]=X
    grid
end

##########################################################
"""
$(TYPEDSIGNATURES)

Constructor for 2D grid
from coordinate arrays.
Boundary region numbers count counterclockwise:

| location  |  number |
| --------- | ------- |
| south     |       1 |
| east      |       2 |
| north     |       3 |
| west      |       4 |

"""
function  simplexgrid(_X::AbstractVector, _Y::AbstractVector)
    X=collect_or_assign(_X)
    Y=collect_or_assign(_Y)
    is_monotone(X) || error("X not monotone")
    is_monotone(Y) || error("Y not monotone")

    Tc=promote_type(eltype(X),eltype(Y))

    function leq(x, x1, x2)
        if (x>x1)
            return false
        end
        if (x>x2)
            return false
        end
        return true
    end

    function geq(x, x1, x2)
        if (x<x1)
            return false
        end
        if (x<x2)
            return false
        end
        return true
    end

    nx=length(X)
    ny=length(Y)

    hmin=X[2]-X[1]
    for i=1:nx-1
        h=X[i+1]-X[i]
        if h <hmin
            hmin=h
        end
    end
    for i=1:ny-1
        h=Y[i+1]-Y[i]
        if h <hmin
            hmin=h
        end
    end

    @assert(hmin>0.0)
    eps=1.0e-5*hmin

    x1=X[1]+eps
    xn=X[nx]-eps
    y1=Y[1]+eps
    yn=Y[ny]-eps

    function  check_insert_bface(ibface,coord,n1,n2)
        if (geq(x1,coord[1,n1],coord[1,n2]))
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
	    bfaceregions[ibface]=4
        elseif (leq(xn,coord[1,n1],coord[1,n2]))
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
	    bfaceregions[ibface]=2
        elseif (geq(y1,coord[2,n1],coord[2,n2]))
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
	    bfaceregions[ibface]=1
        elseif (leq(yn,coord[2,n1],coord[2,n2]))
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
	    bfaceregions[ibface]=3
        end
        ibface
    end

    num_nodes=nx*ny
    num_cells=2*(nx-1)*(ny-1)
    num_bfacenodes=2*(nx-1)+2*(ny-1)

    coord=zeros(Tc,2,num_nodes)
    cellnodes=zeros(Int32,3,num_cells)
    cellregions=zeros(Int32,num_cells)
    bfacenodes=Array{Int32}(undef,2,num_bfacenodes)
    bfaceregions=zeros(Int32,num_bfacenodes)

    ipoint=0
    for iy=1:ny
        for ix=1:nx
            ipoint=ipoint+1
            coord[1,ipoint]=X[ix]
            coord[2,ipoint]=Y[iy]
        end
    end
    @assert(ipoint==num_nodes)

    icell=0
    for iy=1:ny-1
        for ix=1:nx-1
	    ip=ix+(iy-1)*nx
	    p00 = ip
	    p10 = ip+1
	    p01 = ip  +nx
	    p11 = ip+1+nx

            icell=icell+1
            cellnodes[1,icell]=p00
            cellnodes[2,icell]=p10
            cellnodes[3,icell]=p11
            cellregions[icell]=1

            icell=icell+1
            cellnodes[1,icell]=p11
            cellnodes[2,icell]=p01
            cellnodes[3,icell]=p00
            cellregions[icell]=1
        end
    end
    @assert(icell==num_cells)

    #lazy way to  create boundary grid

    ibface=0
    for icell=1:num_cells
        n1=cellnodes[1,icell]
	n2=cellnodes[2,icell]
	n3=cellnodes[3,icell]
        ibface=check_insert_bface(ibface,coord,n1,n2)
	ibface=check_insert_bface(ibface,coord,n1,n3)
	ibface=check_insert_bface(ibface,coord,n2,n3)
    end
    @assert(ibface==num_bfacenodes)

    grid=simplexgrid(coord,
                     cellnodes,
                     cellregions,
                     bfacenodes,
                     bfaceregions)

    grid[XCoordinates]=X
    grid[YCoordinates]=Y
    grid
end

##########################################################
"""
$(TYPEDSIGNATURES)

Constructor for 3D grid
from coordinate arrays.
Boundary region numbers:

| location  |  number |
| --------- | ------- |
| south     |       1 |
| east      |       2 |
| north     |       3 |
| west      |       4 |
| bottom    |       5 |
| top       |       6 |

"""
function  simplexgrid(_X::AbstractVector,_Y::AbstractVector,_Z::AbstractVector)
    X=collect_or_assign(_X)
    Y=collect_or_assign(_Y)
    Z=collect_or_assign(_Z)

    is_monotone(X) || error("X not monotone")
    is_monotone(Y) || error("Y not monotone")
    is_monotone(Z) || error("Z not monotone")

    Tc=promote_type(eltype(X),eltype(Y),eltype(Z))

    leq(x, x1, x2, x3)=x≤x1 && x≤x2 && x≤x3
    geq(x, x1, x2, x3)=x≥x1 && x≥x2 && x≥x3

    nx=length(X)
    ny=length(Y)
    nz=length(Z)

    hmin=X[2]-X[1]
    for i=1:nx-1
        h=X[i+1]-X[i]
        if h <hmin
            hmin=h
        end
    end
    for i=1:ny-1
        h=Y[i+1]-Y[i]
        if h <hmin
            hmin=h
        end
    end

    for i=1:nz-1
        h=Z[i+1]-Z[i]
        if h <hmin
            hmin=h
        end
    end

    @assert(hmin>0.0)
    eps=1.0e-5*hmin

    x1=X[1]+eps
    xn=X[end]-eps

    y1=Y[1]+eps
    yn=Y[end]-eps

    z1=Z[1]+eps
    zn=Z[end]-eps

    function  check_insert_bface(ibface,coord,bfacenodes,n1,n2,n3)
        if geq(x1,coord[1,n1],coord[1,n2],coord[1,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=4
        elseif leq(xn,coord[1,n1],coord[1,n2],coord[1,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=2
        elseif geq(y1,coord[2,n1],coord[2,n2],coord[2,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=1
        elseif leq(yn,coord[2,n1],coord[2,n2],coord[2,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=3
        elseif geq(z1,coord[3,n1],coord[3,n2],coord[3,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=5
        elseif leq(zn,coord[3,n1],coord[3,n2],coord[3,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=6
        end
        ibface
    end

    num_nodes=nx*ny*nz
    num_cells=6*(nx-1)*(ny-1)*(nz-1)
    num_bfacenodes=4*(nx-1)*(ny-1)+4*(nx-1)*(nz-1)+4*(ny-1)*(nz-1)

    Ti=Int32
    coord=zeros(Tc,3,num_nodes)
    cellnodes=zeros(Ti,4,num_cells)
    cellregions=zeros(Ti,num_cells)
    bfacenodes=zeros(Ti,3,num_bfacenodes)
    bfaceregions=zeros(Ti,num_bfacenodes)

    ipoint=0
    for iz=1:nz
        for iy=1:ny
            for ix=1:nx
                ipoint=ipoint+1
                coord[1,ipoint]=X[ix]
                coord[2,ipoint]=Y[iy]
                coord[3,ipoint]=Z[iz]
            end
        end
    end

    @assert(ipoint==num_nodes)

    icell=0
    nxy=nx*ny
    for iz=1:nz-1
        for iy=1:ny-1
            for ix=1:nx-1

	        ip=ix+(iy-1)*nx+(iz-1)*nxy;

	        p000 = ip ;
	        p100 = ip+1;
	        p010 = ip  +nx;
	        p110 = ip+1+nx;
	        p001 = ip	+nxy;
	        p101 = ip+1	+nxy;
	        p011 = ip  +nx+nxy;
	        p111 = ip+1+nx+nxy;

                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p100,p110,p111)
                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p100,p101,p111)
                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p010,p011,p111)
                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p010,p110,p111)
                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p001,p101,p111)
                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p001,p011,p111)
            end
        end
    end
    @assert(icell==num_cells)
    #lazy way to  create boundary grid
    # lazy but easy... bad for partitioning !!!

    ibface=0
    for icell=1:num_cells
        n1=cellnodes[1,icell]
	n2=cellnodes[2,icell]
	n3=cellnodes[3,icell]
	n4=cellnodes[4,icell]
        ibface=check_insert_bface(ibface,coord,bfacenodes,n1,n2,n3)
	ibface=check_insert_bface(ibface,coord,bfacenodes,n1,n2,n4)
	ibface=check_insert_bface(ibface,coord,bfacenodes,n1,n3,n4)
	ibface=check_insert_bface(ibface,coord,bfacenodes,n2,n3,n4)
    end
    @assert(ibface==num_bfacenodes)

    grid=simplexgrid(coord,
                     cellnodes,
                     cellregions,
                     bfacenodes,
                     bfaceregions)

    grid[XCoordinates]=X
    grid[YCoordinates]=Y
    grid[ZCoordinates]=Z
    grid
end

######################################################
"""
$(TYPEDSIGNATURES)

Read grid from file. Currently for pdelib sg format only
"""
function simplexgrid(file::String;format="")
    (fbase,fext)=splitext(file)
    if format==""
        format=fext[2:end]
    end
    @assert format=="sg"

    tks=TokenStream(file)
    expecttoken(tks,"SimplexGrid")
    version=parse(Float64,gettoken(tks))
    version20=false;

    if (version==2.0)
        version20=true;
    elseif (version==2.1)
        version20=false;
    else
        error("Read grid: wrong format version: $(version)")
    end

    dim::Int32=0
    coord=Array{Float64,2}(undef,0,0)
    cells=Array{Int32,2}(undef,0,0)
    regions=Array{Int32,1}(undef,0)
    faces=Array{Int32,2}(undef,0,0)
    bregions=Array{Int32,1}(undef,0)
    while(true)
        if (trytoken(tks,"DIMENSION"))
            dim=parse(Int32,gettoken(tks));
        elseif (trytoken(tks,"NODES"))
            nnodes=parse(Int32,gettoken(tks));
            embdim=parse(Int32,gettoken(tks));
            if(dim!=embdim)
                error("Dimension error (DIMENSION $(dim)) in section NODES")
            end
            coord=Array{Float64,2}(undef,dim,nnodes)
            for inode=1:nnodes
                for idim=1:embdim
                    coord[idim,inode]=parse(Float64,gettoken(tks))
                end
            end
        elseif (trytoken(tks,"CELLS"))
            ncells=parse(Int32,gettoken(tks));
            cells=Array{Int32,2}(undef,dim+1,ncells)
            regions=Array{Int32,1}(undef,ncells)
            for icell=1:ncells
                for inode=1:dim+1
                    cells[inode,icell]=parse(Int32,gettoken(tks));
                end
                regions[icell]=parse(Int32,gettoken(tks));
	        if version20
		    for j=1:dim+1
		        gettoken(tks);  # skip file format garbage
                    end
                end
            end
        elseif (trytoken(tks,"FACES"))
            nfaces=parse(Int32,gettoken(tks));
            faces=Array{Int32,2}(undef,dim,nfaces)
            bregions=Array{Int32,1}(undef,nfaces)
            for iface=1:nfaces
                for inode=1:dim
                    faces[inode,iface]=parse(Int32,gettoken(tks));
                end
                bregions[iface]=parse(Int32,gettoken(tks));
	        if (version20)
		    for j=1:dim+2
		        gettoken(tks); #skip file format garbage
                    end
                end
            end
        else
            expecttoken(tks,"END")
            break
        end
    end
    simplexgrid(coord,cells,regions,faces,bregions);
end

"""
$(TYPEDSIGNATURES)

Write grid to file. Currently for pdelib sg format only
"""
function Base.write(fname::String, g::ExtendableGrid; format="")
    (fbase,fext)=splitext(fname)
    if format==""
        format=fext[2:end]
    end
    @assert format=="sg"

    dim_g=dim_grid(g)
    dim_s=dim_space(g)
    nn=num_nodes(g)
    nc=num_cells(g)
    nbf=num_bfaces(g)
    coord=g[Coordinates]
    cellnodes=g[CellNodes]
    bfacenodes=g[BFaceNodes]
    cellregions=g[CellRegions]
    bfaceregions=g[BFaceRegions]

    # TODO: replace @sprintf by someting non-allocating
    open(fname, "w") do file
        write(file,@sprintf("SimplexGrid"))
        write(file,@sprintf(" "))
        write(file,@sprintf("2.1\n"))
        write(file,@sprintf("#created by ExtendableGrids.jl (c) J.Fuhrmann et al\n"))
        write(file,@sprintf("#mailto:{fuhrmann|streckenbach}@wias-berlin.de\n"))
        write(file,@sprintf("#%s\n",Dates.format(Dates.now(),"yyyy-mm-ddTHH-mm-SS")))

        write(file,@sprintf("DIMENSION\n%d\n",dim_g))
        write(file,@sprintf("NODES\n%d %d\n",nn,dim_s))

        for inode=1:nn
            for idim=1:dim_s
	        write(file,@sprintf("%.20e ",coord[idim,inode]))
                write(file,@sprintf("\n"))
            end
        end

        write(file,@sprintf("CELLS\n%d\n",nc))
        for icell=1:nc
            for inode=1:dim_g+1
	        write(file,@sprintf("%d ", cellnodes[inode,icell]))
            end
            write(file,@sprintf("%d\n",cellregions[icell]))
        end

        write(file,@sprintf("FACES\n%d\n",nbf))
        for ibface=1:nbf
            for inode=1:dim_g
	        write(file,@sprintf("%d ", bfacenodes[inode,ibface]))
            end
            write(file,@sprintf("%d\n",bfaceregions[ibface]))
        end
        write(file,@sprintf("END\n"))
        flush(file)
        flush(file)
    end
    nothing
end


text:
export PendulumNonInteractiveEnv

struct PendulumNonInteractiveEnvParams{Fl<:AbstractFloat}
    gravity::Fl
    length::Fl
    mass::Fl
    step_size::Fl
    maximum_time::Fl
end

"""
A non-interactive pendulum environment.

Accepts only `nothing` actions, which result in the system being simulated for one time step.
Sets `env.done` to `true` once `maximum_time` is reached. Resets to a random position and momentum.
Always returns zero rewards.

Useful for debugging and development purposes, particularly in model-based reinforcement learning.
"""
mutable struct PendulumNonInteractiveEnv{
    Fl<:AbstractFloat,
    VFl<:AbstractVector{Fl},
    R<:AbstractRNG,
} <: NonInteractiveEnv
    parameters::PendulumNonInteractiveEnvParams{Fl}
    state::VFl
    done::Bool
    t::Int
    rng::R
end

"""
    PendulumNonInteractiveEnv(;kwargs...)

# Keyword arguments

- `float_type = Float64`
- `gravity = 9.8`
- `length = 2.0`
- `mass = 1.0`
- `step_size = 0.01`
- `maximum_time = 10.0`
- `rng = Random.GLOBAL_RNG`
"""
function PendulumNonInteractiveEnv(;
    float_type = Float64,
    gravity = 9.8,
    length = 2.0,
    mass = 1.0,
    step_size = 0.01,
    maximum_time = 10.0,
    rng = Random.GLOBAL_RNG,
)
    parameters = PendulumNonInteractiveEnvParams{float_type}(
        gravity,
        length,
        mass,
        step_size,
        maximum_time,
    )
    env = PendulumNonInteractiveEnv(parameters, zeros(float_type, 2), false, 0, rng)
    reset!(env)
    env
end

Random.seed!(env::PendulumNonInteractiveEnv, seed) = Random.seed!(env.rng, seed)

RLBase.reward(env::PendulumNonInteractiveEnv) = 0
RLBase.is_terminated(env::PendulumNonInteractiveEnv) = env.done
RLBase.state(env::PendulumNonInteractiveEnv) = env.state
RLBase.state_space(env::PendulumNonInteractiveEnv{T}) where {T} =
    Space([typemin(T)..typemax(T), typemin(T)..typemax(T)])

function RLBase.reset!(env::PendulumNonInteractiveEnv{Fl}) where {Fl}
    env.state .= (Fl(2 * pi) * rand(env.rng, Fl), randn(env.rng, Fl))
    env.t = 0
    env.done = false

    nothing
end

function (env::PendulumNonInteractiveEnv{Fl})(a::Nothing) where {Fl}
    (g, l, m) = (env.parameters.gravity, env.parameters.length, env.parameters.mass)
    (dt, T) = (env.parameters.step_size, env.parameters.maximum_time)
    (theta, p_theta) = env.state

    # H = p_theta^2 / (2*m*l^2) + m*g*l*(1 - cos(theta))
    p_theta -= (dt / 2) * m * g * l * sin(theta)
    theta += dt * p_theta / (m * l^2)
    p_theta -= (dt / 2) * m * g * l * sin(theta)

    theta = mod(theta, 2 * Fl(pi))

    env.t += 1
    env.state .= (theta, p_theta)
    env.done = env.t * dt > T ? true : false

    nothing
end


text:
# Precompose with a linear mapping + translation (= affine)

export Precompose

"""
**Precomposition with linear mapping/translation**

    Precompose(f, L, μ, b)

Returns the function
```math
g(x) = f(Lx + b)
```
where ``f`` is a convex function and ``L`` is a linear mapping: this must satisfy ``LL^* = μI`` for ``μ ⩾ 0``. Furthermore, either ``f`` is separable or parameter `μ` is a scalar, for the `prox` of ``g`` to be computable.

Parameter `L` defines ``L`` through the `A_mul_B!` and `Ac_mul_B!` methods. Therefore `L` can be an `AbstractMatrix` for example, but not necessarily.

In this case, `prox` and `prox!` are computed according to Prop. 24.14 in Bauschke, Combettes "Convex Analisys and Monotone Operator Theory in Hilbert Spaces", 2nd edition, 2016. The same result is Prop. 23.32 in the 1st edition of the same book.
"""
struct Precompose{T <: ProximableFunction, R <: Real, C <: Union{R, Complex{R}}, U <: Union{C, AbstractArray{C}}, V <: Union{C, AbstractArray{C}}, M} <: ProximableFunction
  f::T
  L::M
  mu::U
  b::V
  function Precompose{T, R, C, U, V, M}(f::T, L::M, mu::U, b::V) where {T <: ProximableFunction, R <: Real, C <: Union{R, Complex{R}}, U <: Union{R, AbstractArray{R}}, V <: Union{C, AbstractArray{C}}, M}
    if !is_convex(f)
      error("f must be convex")
    end
    if any(mu .<= 0.0)
      error("elements of μ must be positive")
    end
    new(f, L, mu, b)
  end
end

is_prox_accurate(f::Precompose) = is_prox_accurate(f.f)
is_convex(f::Precompose) = is_convex(f.f)
is_set(f::Precompose) = is_set(f.f)
is_singleton(f::Precompose) = is_singleton(f.f)
is_cone(f::Precompose) = is_cone(f.f)
is_affine(f::Precompose) = is_affine(f.f)
is_smooth(f::Precompose) = is_smooth(f.f)
is_quadratic(f::Precompose) = is_quadratic(f.f)
is_generalized_quadratic(f::Precompose) = is_generalized_quadratic(f.f)
is_strongly_convex(f::Precompose) = is_strongly_convex(f.f)

Precompose(f::T, L::M, mu::U, b::V) where {T <: ProximableFunction, R <: Real, C <: Union{R, Complex{R}}, U <: Union{R, AbstractArray{R}}, V <: Union{C, AbstractArray{C}}, M} = Precompose{T, R, C, U, V, M}(f, L, mu, b)

Precompose(f::T, L::M, mu::U) where {T <: ProximableFunction, R <: Real, U <: Union{R, AbstractArray{R}}, M} = Precompose{T, R, R, U, R, M}(f, L, mu, zero(R))

function (g::Precompose)(x::T) where {T <: Union{Tuple, AbstractArray}}
  return g.f(g.L*x .+ g.b)
end

function gradient!(y::AbstractArray{T}, g::Precompose, x::AbstractArray{T}) where T <: RealOrComplex
  res = g.L*x .+ g.b
  gradres = similar(res)
  v = gradient!(gradres, g.f, res)
  mul!(y, adjoint(g.L), gradres)
  return v
end

function prox!(y::AbstractArray{C}, g::Precompose, x::AbstractArray{C}, gamma::R=one(R)) where {R <: Real, C <: Union{R, Complex{R}}}
  # See Prop. 24.14 in Bauschke, Combettes "Convex Analisys and Monotone Operator Theory in Hilbert Spaces", 2nd ed., 2016.
  # The same result is Prop. 23.32 in the 1st ed. of the same book.
  #
  # This case has an additional translation: if f(x) = h(x + b) then
  #   prox_f(x) = prox_h(x + b) - b
  # Then one can apply the above mentioned result to g(x) = f(Lx).
  #
  res = g.L*x .+ g.b
  proxres = similar(res)
  v = prox!(proxres, g.f, res, g.mu.*gamma)
  proxres .-= res
  proxres ./= g.mu
  mul!(y, adjoint(g.L), proxres)
  y .+= x
  return v
end

function prox_naive(g::Precompose, x::AbstractArray{C}, gamma::R=1.0) where {R <: Real, C <: Union{R, Complex{R}}}
  res = g.L*x .+ g.b
  proxres, v = prox_naive(g.f, res, g.mu .* gamma)
  y = x + g.L'*((proxres .- res)./g.mu)
  return y, v
end

fun_name(f::Precompose) = "Precomposition"
fun_dom(f::Precompose) = fun_dom(f.f)
fun_expr(f::Precompose) = "x ↦ f(L(x) + b)"


text:
module Trajectories

using ..UpdateMethods

import Base:push!, length

export Trajectory, heun_fill, heun_fill!
mutable struct Trajectory
    t::Vector{<:Real}
    x::Vector{<:Real}
    v::Union{Vector{<:Real}, Nothing}
end
Trajectory(t::Vector{<:Real}, x::Vector{<:Real}) = Trajectory(t, x, nothing)
Trajectory(T::Type) = Trajectory(Vector{T}(), Vector{T}(), nothing)
Trajectory() = Trajectory(Float64)

function push!(tj::Trajectory, t::Real, x::Real)
    push!(tj.t, t)
    push!(tj.x, x)
    tj
end

function genrate_v!(tj::Trajectory)
    tj.v = (tj.x[3:end] - tj.x[1:end-2]) ./ (tj.t[3:end] - tj.t[1:end-2])
end

lastx(tj::Trajectory)::Real = tj.x[end]

function length(tj::Trajectory)
    length(tj.x)
end

function heun_fill!(tj::Trajectory, h::Heun, f::Function, trange::Tuple{T, T})::Trajectory where T<:Real
    t0, tend = trange
    for t in t0:h.step:tend
        push!(tj, heun_update(f, h, lastx(tj), t)...)
    end
    tj
end

function heun_fill(x0::Real, h::Heun, f::Function, trange::Tuple{T, T})::Trajectory where T<:Real
    heun_fill!(Trajectory([trange[1]],[x0]), h, f, trange)
end

end


text:
using LPVSpectral, Plots, DSP, ProximalOperators

@testset "lasso" begin
    function generate_signal(f,w,N, modphase=false)
        x = sort(10rand(N)) # Sample points
        v = range(0, stop=1, length=N) # Scheduling variable

        # generate output signal
        dependence_matrix = Float64[f[(i-1)%length(f)+1](v) for v in v, i in eachindex(w)] # N x nw
        frequency_matrix  = [cos(w*x -0.5modphase*(dependence_matrix[i,j])) for (i,x) in enumerate(x), (j,w) in enumerate(w)] # N x nw
        y = sum(dependence_matrix.*frequency_matrix,dims=2)[:] # Sum over all frequencies
        y += 0.1randn(size(y))
        y,v,x,frequency_matrix, dependence_matrix
    end

    N      = 500 # Number of training data points
    f      = [v->2v^2, v->2/(5v+1), v->3exp(-10*(v-0.5)^2),] # Functional dependences on the scheduling variable
    w      = 2π*[2,10,20] # Frequency vector
    w_test = 2π*(2:2:25) # Test Frequency vector, set w_test = w For a nice Function visualization

    Y,V,X,frequency_matrix, dependence_matrix = generate_signal(f,w,N, true)

    λ      = 0.02 # Regularization parameter
    λs     = 5    # Regularization parameter group-lasso
    normal = true # Use normalized basis functions
    Nv     = 50   # Number of basis functions

    zeronan(x) = ifelse(x==0, NaN, x)
    callback(x,z) = ()#plot(max.(abs2.([complex.(x[1:end÷2], x[end÷2+1:end]) complex.(z[1:end÷2], z[end÷2+1:end])]), 1e-20), reuse=true, show=true)

    ses = ls_sparse_spectral_lpv(Y,X,V,w_test,Nv; λ = λs, normalize = normal, tol=1e-8, printerval=10, iters=2000, cb=callback) # Perform LPV spectral estimation
    se  = ls_spectral_lpv(Y,X,V,w_test,Nv; λ = 0.02, normalize = normal)
    xs  = LPVSpectral.ls_sparse_spectral(Y,X,1:0.1:25; λ=20, tol=1e-9, printerval=2, iters=9000, μ=0.00001,cb=callback)
    # xsi = LPVSpectral.ls_sparse_spectral(Y,X,1:0.1:25; proxg=IndBallL0(6), λ=0.5, tol=1e-9, printerval=100, iters=9000, μ=0.0001,cb=callback)
    xsw = ls_windowpsd(Y,X,1:0.5:22; nw=2, estimator=ls_sparse_spectral, λ=0.2, tol=1e-10, printerval=10000, iters=60000, μ=0.0001)

    # plot(X,[Y V], linewidth=[1 2], lab=["\$y_t\$" "\$v_t\$"], xlabel=L"$x$ (sampling points)", title=L"Test signal $y_t$ and scheduling signal $v_t$", legend=true, xlims=(0,10), grid=false, c=[:cyan :blue])
    plot(se; normalization=:none, dims=2, l=:solid, c = :orange, fillalpha=0.5, nMC = 100, fillcolor=:orange, linewidth=2, bounds=true, lab=reshape(["Est. \$\\omega = $(round(w/π))\\pi \$" for w in w_test],1,:), phase = false)
    plot!(ses; normalization=:none, dims=2, l=:solid, c = :green, linewidth=2, lab=reshape(["Est. \$\\omega = $(round(w/π))\\pi \$" for w in w_test],1,:), phase = false)
    plot!(V,dependence_matrix, title="Functional dependencies \$A(\\omega,v)\$", xlabel="\$v\$", ylabel="\$A(\\omega,v)\$", c = [:blue], l=:dot, linewidth=2,lab=reshape(["True \$\\omega = $(round(w/π))\\pi\$" for w in w],1,:), grid=false)

    ## Plot regular spectrum
    # spectrum_lpv   = psd(se) # Calculate power spectral density
    # spectrum_lpvs  = psd(ses) # Calculate sparse power spectral density
    # fs             = N/(X[end]-X[1]) # This is the (approximate) sampling freqency of the generated signal
    # spectrum_per   = DSP.periodogram(Y, fs=fs)
    # spectrum_welch = DSP.welch_pgram(Y, fs=fs)
    # plot(2π*collect(spectrum_per.freq), spectrum_per.power, lab="Periodogram", l=:path, m=:none, yscale=:log10, c=:cyan, legend=:bottomright)
    # plot!(2π*collect(spectrum_welch.freq), spectrum_welch.power, lab="Welch", l=:path, m=:none, yscale=:log10, linewidth=2, c=:blue)
    # plot!(2π*(1:0.1:25), zeronan.(abs2.(xs)), lab="sparse", l=:path, m=:none, yscale=:log10, linewidth=2, c=:magenta)
    # plot!(w_test,spectrum_lpv/fs, xlabel=L"\omega [rad/s]", ylabel="Spectral density", ylims=(-Inf,Inf), grid=false, lab="LPV", l=:scatter, m=:o, yscale=:log10, c=:orange)
    # plot!(w_test,zeronan.(spectrum_lpvs)./fs, lab="Sparse LPV", l=:scatter, m=:x, c=:green)

    # plot!(2π*(1:0.1:25), max.(abs2.(xsi), 1e-15), lab="sparse ind ball", l=:path, m=:none, yscale=:log10, linewidth=2, c=:yellow)
    # plot!(2π*(1:0.5:22), zeronan.(abs2.(xsw)), lab="sparse windowed", l=:path, m=:none, yscale=:log10, linewidth=2, c=:orange)
    # savetikz("/local/home/fredrikb/phdthesis/spectralest/figs/spectrum_gen2.tex")
end


text:
"""
set_ground_transit_graph!(env::CoordinatedMAPFEnv, drone_sp_dists::Vector{Float64})

Generate the road-and-transit graph as the underlying Stage 2 MAPF graph
"""
function set_ground_transit_graph!(env::CoordinatedMAPFEnv, drone_sp_dists::Vector{Float64})

    gtg = env.ground_transit_graph
    weights_dict = env.ground_transit_weights  # Tracks the edge weights as added
    expanded_paths_dict = env.ground_transit_expanded_paths
    drone_sg_idxs = env.gtg_drone_start_goal_idxs    # Notes the start-goal idxs of drones relative to gtg
    ground_path_start_ends = env.gtg_ground_path_start_ends

    v_idx = 1

    # First add drone start and goal vertices
    for (ai, atask) in enumerate(env.aerial_task_list)
        push!(drone_sg_idxs, (v_idx, v_idx+1))

        # Compute road graph distance between drone starts and goals
        weights_dict[(v_idx, v_idx+1)] = drone_sp_dists[ai]

        Graphs.add_vertex!(gtg, GroundTransitState(idx=v_idx, road_vtx_id=atask.origin))
        v_idx += 1
        Graphs.add_vertex!(gtg, GroundTransitState(idx=v_idx, road_vtx_id=atask.dest))
        v_idx += 1
    end

    # Loop over ground paths
    for (gi, gpath) in enumerate(env.ground_paths)

        # Now loop over agent's path
        ground_start = v_idx
        for (path_idx, state) in enumerate(gpath.path_states)

            new_gts = GroundTransitState(idx=v_idx, ground_agent_id=gi, vertex_along_path=path_idx, road_vtx_id=state.road_vtx_id, timeval=state.timeval)
            Graphs.add_vertex!(gtg, new_gts)

            # Compute to and from paths relative to drone starts and goals over ROAD GRAPH
            # And update the gtg weights correctly
            for (ai, (ais, aig)) in enumerate(drone_sg_idxs)

                atask = env.aerial_task_list[ai]

                if atask.origin == state.road_vtx_id
                    weights_dict[(ais, v_idx)] = 0.0
                else
                    bidir_sp_to = compute_bidir_astar_euclidean(env.road_graph, atask.origin, state.road_vtx_id, env.road_graph_wts, env.location_list)
                    weights_dict[(ais, v_idx)] = bidir_sp_to.dist
                    expanded_paths_dict[(ais, v_idx)] = bidir_sp_to.path
                end

                if state.road_vtx_id == atask.dest
                    weights_dict[(v_idx, aig)] = 0.0
                else
                    bidir_sp_from = compute_bidir_astar_euclidean(env.road_graph, state.road_vtx_id, atask.dest, env.road_graph_wts, env.location_list)
                    weights_dict[(v_idx, aig)] = bidir_sp_from.dist
                    expanded_paths_dict[(v_idx, aig)] = bidir_sp_from.path
                end
            end

            v_idx += 1
        end # enumerate(gpath.path_states)
        ground_end = v_idx-1
        push!(ground_path_start_ends, (ground_start, ground_end))
    end

end # function generate_ground_transit_graph

## Stage 2 MAPF functions
function MultiAgentPathFinding.add_constraint!(cbase::GroundTransitConstraint, cadd::GroundTransitConstraint)
    union!(cbase.avoid_gid_vertex_set, cadd.avoid_gid_vertex_set)
end

function MultiAgentPathFinding.overlap_between_constraints(cbase::GroundTransitConstraint, cadd::GroundTransitConstraint)
    return ~(isempty(intersect(cbase.avoid_gid_vertex_set, cadd.avoid_gid_vertex_set)))
end

MultiAgentPathFinding.get_mapf_state_from_idx(env::CoordinatedMAPFEnv, idx::Int64) = env.aerial_mapf_graph.vertices[idx]

function get_mapf_transit_action(env::CoordinatedMAPFEnv, u::AerialMAPFState, v::AerialMAPFState)
    ugtg = env.ground_transit_graph.vertices[u.ground_transit_idx]
    vgtg = env.ground_transit_graph.vertices[v.ground_transit_idx]
    if ugtg.ground_agent_id == vgtg.ground_agent_id
        return GroundTransitAction(Stay)
    else
        return GroundTransitAction(Fly)
    end
end

function MultiAgentPathFinding.get_mapf_action(env::CoordinatedMAPFEnv, u::Int64, v::Int64)
    vtx_u = env.aerial_mapf_graph.vertices[u]
    vtx_v = env.aerial_mapf_graph.vertices[v]

    return get_mapf_transit_action(env, vtx_u, vtx_v)
end

function MultiAgentPathFinding.set_low_level_context!(::CoordinatedMAPFEnv, ::Int64, ::GroundTransitConstraint)
end

MultiAgentPathFinding.get_empty_constraint(::Type{GroundTransitConstraint}) = GroundTransitConstraint()

function MultiAgentPathFinding.get_first_conflict(env::CoordinatedMAPFEnv, solution::Vector{PlanResult{AerialMAPFState,GroundTransitAction,Float64}})

    gtg_idx_to_aerial_ids = Dict{Int64,Vector{Int64}}()

    for (i, sol_i) in enumerate(solution)
        for ((si, _), (sip,_)) in zip(sol_i.states[1:end-1], sol_i.states[2:end])
            si_gtg = env.ground_transit_graph.vertices[si.ground_transit_idx]
            sip_gtg = env.ground_transit_graph.vertices[sip.ground_transit_idx]
            if si_gtg.ground_agent_id !=0 && si_gtg.ground_agent_id == sip_gtg.ground_agent_id
                if ~haskey(gtg_idx_to_aerial_ids, si_gtg.idx)
                    gtg_idx_to_aerial_ids[si_gtg.idx] = Int64[]
                end
                push!(gtg_idx_to_aerial_ids[si_gtg.idx], i)
            end
        end
    end

    # Whenever we hit a key with more than capacity drones, continue until fewer than capacity OR non-consec IDXs
    gid_to_aerial_ids = Dict{Int64,Vector{Int64}}()
    for (gtgidx, aerial_ids) in gtg_idx_to_aerial_ids

        if length(aerial_ids) > env.car_capacity

            gid = env.ground_transit_graph.vertices[gtgidx].ground_agent_id

            # Now just continue until different running agent set
            if ~haskey(gid_to_aerial_ids, gid)
                gid_to_aerial_ids[gid] = aerial_ids
            end
        end
    end

    if ~isempty(gid_to_aerial_ids)

        env.num_global_conflicts += 1
        println("Aerial CBS conflict number $(env.num_global_conflicts)")

        if env.num_global_conflicts > env.threshold_global_conflicts
            throw(DomainError("Too many conflicts in Stage 2!"))
        end

        overlap_gids = collect(keys(gid_to_aerial_ids))
        aerial_agent_ids = collect(values(gid_to_aerial_ids))

        # @infiltrate
        return GroundTransitConflict(overlap_gids, aerial_agent_ids)
    end

    return nothing
end # function get_first_conflict

function MultiAgentPathFinding.create_constraints_from_conflict(env::CoordinatedMAPFEnv, conflict::GroundTransitConflict)

    # Generate all NC2
    aerial_gids_to_avoid = Dict{Int64,Set{Int64}}()

    # Get all car IDs from overlap_gtg_vertices
    for (gid, aids) in zip(conflict.overlap_gids, conflict.aerial_agent_ids)
        for aid in aids
            if ~haskey(aerial_gids_to_avoid, aid)
                aerial_gids_to_avoid[aid] = Set{Int64}(gid)
            else
                push!(aerial_gids_to_avoid[aid], gid)
            end
        end
    end

    constraint = Dict{Int64,GroundTransitConstraint}(aid => GroundTransitConstraint(gids) for (aid,gids) in aerial_gids_to_avoid)

    return [constraint]
end

function MultiAgentPathFinding.low_level_search!(solver::ECBSSolver, aerial_agent_id::Int64, s::AerialMAPFState, constraint::GroundTransitConstraint, solution::Vector{PlanResult{AerialMAPFState,GroundTransitAction,Float64}})

    env = solver.env

    env.aerial_mapf_graph = SimpleVListGraph{AerialMAPFState}()
    Graphs.add_vertex!(env.aerial_mapf_graph, s)

    env.gtg_idx_gval =  Dict{Int64,Float64}()

    edge_wt_fn(u, v) = amapfg_edge_wt_fn(env, u, v)
    admissible_heuristic(u) = direct_flight_time_heuristic(env, u)

    # Set start and goal
    start_idx = 1   # Always one since graph reinitialized
    env.next_gtg_goal_idx = env.gtg_drone_start_goal_idxs[aerial_agent_id][2]
    env.next_amapfg_goal_idx = 0

    vis = GroundTransitGoalVisitor(env, constraint.avoid_gid_vertex_set)

    a_star_states = a_star_implicit_shortest_path!(env.aerial_mapf_graph,
                                                    edge_wt_fn,
                                                    start_idx,
                                                    vis,
                                                    admissible_heuristic
                                                    )

    # Get solution from a_star_eps_states
    # NOTE: It should always have a solution because drone can just fly
    @assert haskey(a_star_states.parent_indices, env.next_amapfg_goal_idx)

    # Get path cost
    states = Tuple{AerialMAPFState,Float64}[]
    actions = Tuple{GroundTransitAction, Float64}[]

    # Get path cost
    sp_cost = a_star_states.dists[env.next_amapfg_goal_idx]

    # Get path
    sp_idxs = shortest_path_indices(a_star_states.parent_indices, env.aerial_mapf_graph, start_idx, env.next_amapfg_goal_idx)

    push!(states, (s, 0.0))

    # Traverse sp_idxs and set up states and actions
    for i = 2:length(sp_idxs)

        u, v =  (sp_idxs[i-1], sp_idxs[i])

        # NOTE: Create new state because we are going to delete graph later
        new_state = deepcopy(env.aerial_mapf_graph.vertices[v])
        push!(states, (new_state, a_star_states.dists[v]))

        # All derived info but can be useful
        push!(actions, (get_mapf_action(env, u, v), a_star_states.dists[v] - a_star_states.dists[u]))
    end

    return PlanResult{AerialMAPFState,GroundTransitAction,Float64}(states, actions, sp_cost, sp_cost)
end # function low_level_search

function MultiAgentPathFinding.low_level_search!(solver::CBSSolver, aerial_agent_id::Int64, s::AerialMAPFState, constraint::GroundTransitConstraint)

    env = solver.env

    env.aerial_mapf_graph = SimpleVListGraph{AerialMAPFState}()
    Graphs.add_vertex!(env.aerial_mapf_graph, s)

    env.gtg_idx_gval =  Dict{Int64,Float64}()

    edge_wt_fn(u, v) = amapfg_edge_wt_fn(env, u, v)
    admissible_heuristic(u) = direct_flight_time_heuristic(env, u)

    # Set start and goal
    start_idx = 1   # Always one since graph reinitialized
    env.next_gtg_goal_idx = env.gtg_drone_start_goal_idxs[aerial_agent_id][2]
    env.next_amapfg_goal_idx = 0

    vis = GroundTransitGoalVisitor(env, constraint.avoid_gid_vertex_set)

    a_star_states = a_star_implicit_shortest_path!(env.aerial_mapf_graph,
                                                    edge_wt_fn,
                                                    start_idx,
                                                    vis,
                                                    admissible_heuristic
                                                    )

    # Get solution from a_star_eps_states
    # NOTE: It should always have a solution because drone can just fly
    @assert haskey(a_star_states.parent_indices, env.next_amapfg_goal_idx)

    # Get path cost
    states = Tuple{AerialMAPFState,Float64}[]
    actions = Tuple{GroundTransitAction, Float64}[]

    # Get path cost
    sp_cost = a_star_states.dists[env.next_amapfg_goal_idx]

    # Get path
    sp_idxs = shortest_path_indices(a_star_states.parent_indices, env.aerial_mapf_graph, start_idx, env.next_amapfg_goal_idx)

    push!(states, (s, 0.0))

    # Traverse sp_idxs and set up states and actions
    for i = 2:length(sp_idxs)

        u, v =  (sp_idxs[i-1], sp_idxs[i])

        # NOTE: Create new state because we are going to delete graph later
        new_state = deepcopy(env.aerial_mapf_graph.vertices[v])
        push!(states, (new_state, a_star_states.dists[v]))

        # All derived info but can be useful
        push!(actions, (get_mapf_action(env, u, v), a_star_states.dists[v] - a_star_states.dists[u]))
    end

    return PlanResult{AerialMAPFState,GroundTransitAction,Float64}(states, actions, sp_cost, sp_cost)
end # function low_level_search

"""
Supports a scaled average of distance and time, but we use only distance in practice, i.e., alpha_weight_distance = 1.0
"""
function amapfg_edge_wt_fn(env::CoordinatedMAPFEnv, u::AerialMAPFState, v::AerialMAPFState)

    ddiff = v.distval - u.distval
    tdiff = v.timeval - u.timeval

    return env.alpha_weight_distance*ddiff + (1.0 - env.alpha_weight_distance)*tdiff
end

function direct_flight_time_heuristic(env::CoordinatedMAPFEnv, u::AerialMAPFState)

    if u.ground_transit_idx == env.next_gtg_goal_idx
        return 0.0
    else
        ddiff = env.ground_transit_weights[(u.ground_transit_idx, env.next_gtg_goal_idx)]
        return (1.0 - env.alpha_weight_distance)*ddiff/AvgSpeed
    end
end

struct GroundTransitGoalVisitor <: Graphs.AbstractDijkstraVisitor
    env::CoordinatedMAPFEnv
    avoid_gid_vertex_set::Set{Int64}
end

function Graphs.include_vertex!(vis::GroundTransitGoalVisitor, u::AerialMAPFState, v::AerialMAPFState, d::Float64, nbrs::Vector{Int64})

    env = vis.env
    avoid_gid_vertex_set = vis.avoid_gid_vertex_set
    gtg_weights = env.ground_transit_weights

    # When you expand goal, return
    if v.ground_transit_idx == env.next_gtg_goal_idx
        env.next_amapfg_goal_idx = v.idx
        return false
    end

    # Now lookup where drone is on the GTG
    ugtg = env.ground_transit_graph.vertices[u.ground_transit_idx]
    vgtg = env.ground_transit_graph.vertices[v.ground_transit_idx]
    idx_ctr = Graphs.num_vertices(env.aerial_mapf_graph)

    # Case 1: Drone start
    if vgtg.ground_agent_id == 0

        ddiff = gtg_weights[(v.ground_transit_idx, env.next_gtg_goal_idx)]
        tdiff = ddiff / AvgSpeed

        idx_ctr += 1
        new_goal_vert = AerialMAPFState(idx=idx_ctr, ground_transit_idx=env.next_gtg_goal_idx, timeval=v.timeval+tdiff, distval=v.distval+ddiff)

        # Add to graph and add neighbor
        Graphs.add_vertex!(env.aerial_mapf_graph, new_goal_vert)
        push!(nbrs, idx_ctr)

        env.gtg_idx_gval[env.next_gtg_goal_idx] = env.alpha_weight_distance*new_goal_vert.distval + (1.0 - env.alpha_weight_distance)*new_goal_vert.timeval

        for (gid, (gstart, gend)) in enumerate(env.gtg_ground_path_start_ends)

            if ~(gid in avoid_gid_vertex_set)

                gts_tup = get_best_transit_connection(env, gstart, gend, v, false)

                if ~isnothing(gts_tup)

                    idx_ctr += 1
                    new_vert = AerialMAPFState(idx=idx_ctr, ground_transit_idx=gts_tup.idx, timeval=gts_tup.timeval, distval=gts_tup.distval)

                    Graphs.add_vertex!(env.aerial_mapf_graph, new_vert)
                    push!(nbrs, idx_ctr)

                    env.gtg_idx_gval[gts_tup.idx] = env.alpha_weight_distance*new_vert.distval + (1.0 - env.alpha_weight_distance)*new_vert.timeval

                end # isnothing(gts_tup)
            end # ~(gid in avoid_gid_vertex_set)
        end

    else
        # Either first transfer to new ground point OR along the way
        if vgtg.ground_agent_id != ugtg.ground_agent_id

            new_considered_cars = union(v.considered_cars, Set{Int64}([vgtg.ground_agent_id]))

            next_gtg = env.ground_transit_graph.vertices[vgtg.idx + 1]

            idx_ctr += 1

            # IMP NOTE: Timeval may be greater than vgtg.timeval (probably will be)
            # So the new timeval should be v.timeval + car_time_diff (next_gtg.timeval - vgtg.timeval)
            tval = v.timeval + (next_gtg.timeval - vgtg.timeval)
            new_vert = AerialMAPFState(idx=idx_ctr, ground_transit_idx=next_gtg.idx, timeval=tval, distval=v.distval, considered_cars=new_considered_cars)

            Graphs.add_vertex!(env.aerial_mapf_graph, new_vert)
            push!(nbrs, idx_ctr)

        else
            # Ongoing car route vertex. First try to add goal as usual
            ddiff = gtg_weights[(v.ground_transit_idx, env.next_gtg_goal_idx)]
            tdiff = ddiff / AvgSpeed
            goal_gval = env.alpha_weight_distance*(v.distval+ddiff) + (1.0 - env.alpha_weight_distance)*(v.timeval + tdiff)

            if goal_gval < env.gtg_idx_gval[env.next_gtg_goal_idx]

                idx_ctr += 1
                new_goal_vert = AerialMAPFState(idx=idx_ctr, ground_transit_idx=env.next_gtg_goal_idx, timeval=v.timeval+tdiff, distval=v.distval+ddiff, considered_cars=v.considered_cars)

                # Add to graph and add neighbor
                Graphs.add_vertex!(env.aerial_mapf_graph, new_goal_vert)
                push!(nbrs, idx_ctr)

                env.gtg_idx_gval[env.next_gtg_goal_idx] = goal_gval
            end

            # Then add next along car route IF not last
            if vgtg.idx < env.gtg_ground_path_start_ends[vgtg.ground_agent_id][2]
                next_gtg = env.ground_transit_graph.vertices[vgtg.idx + 1]

                idx_ctr += 1
                tval = v.timeval + (next_gtg.timeval - vgtg.timeval)
                new_vert = AerialMAPFState(idx=idx_ctr, ground_transit_idx=next_gtg.idx, timeval=tval, distval=v.distval, considered_cars=v.considered_cars)

                Graphs.add_vertex!(env.aerial_mapf_graph, new_vert)
                push!(nbrs, idx_ctr)

            end

            for (gid, (gstart, gend)) in enumerate(env.gtg_ground_path_start_ends)

                if gid != vgtg.ground_agent_id && ~(gid in v.considered_cars) && ~(gid in avoid_gid_vertex_set)

                    gts_tup = get_best_transit_connection(env, gstart, gend, v, true)

                    # Add AerialMAPFState versions of them
                    if ~isnothing(gts_tup)

                        true_timeval = gts_tup.timeval
                        true_distval = gts_tup.distval

                        temp_gval = env.alpha_weight_distance*true_distval + (1.0 - env.alpha_weight_distance)*true_timeval

                        if ~haskey(env.gtg_idx_gval, gts_tup.idx) || temp_gval < env.gtg_idx_gval[gts_tup.idx]

                            idx_ctr += 1
                            new_vert = AerialMAPFState(idx=idx_ctr, ground_transit_idx=gts_tup.idx, timeval=true_timeval, distval=true_distval, considered_cars=v.considered_cars)

                            Graphs.add_vertex!(env.aerial_mapf_graph, new_vert)
                            push!(nbrs, idx_ctr)

                            env.gtg_idx_gval[gts_tup.idx] = temp_gval
                        end
                    end # isnothing(gts_tup)
                end
            end # (gstart, gend) in gtg_ground_path_start_ends
        end # vgtg.ground_agent_id != ugtg.ground_agent_id
    end # vgtg.ground == or != 0

    return true
end

"""
    update_aerial_ground_paths_cbs_solution!(env::CoordinatedMAPFEnv, solution::Vector{PlanResult}, drone_flight_paths::Vector{AgentPathInfo})

Reconcile the time and distance value of all the aerial-ground paths at the end of the Stage 2 MAPF.
"""
function update_aerial_ground_paths_cbs_solution!(env::CoordinatedMAPFEnv, solution::Vector{PlanResult{AerialMAPFState,GroundTransitAction,Float64}}, drone_flight_paths::Vector{AgentPathInfo})

    final_aerial_paths = Vector{AgentPathInfo}(undef, length(solution))

    # Iterate over solution
    for (aerial_agent_id, aerial_soln) in enumerate(solution)

        # Direct flight path
        if length(aerial_soln.states) == 2
            # Just copy over drone flight paths
            final_aerial_paths[aerial_agent_id] = drone_flight_paths[aerial_agent_id]
        else
            # Iterate over drone car usage
            # If time val of boarding is greater than drone's timestamp, edit drone path
            # If time val of drone is greater than car vertex then edit car path
            path_states = AgentPathState[]
            running_time = 0.0
            running_dist = 0.0

            for i = 1:length(aerial_soln.states) - 1

                si, sip = (aerial_soln.states[i][1], aerial_soln.states[i+1][1])

                si_gtg = env.ground_transit_graph.vertices[si.ground_transit_idx]
                sip_gtg = env.ground_transit_graph.vertices[sip.ground_transit_idx]

                push!(path_states, AgentPathState(road_vtx_id=si_gtg.road_vtx_id, timeval=running_time, coord_agent_id=Set{Int64}(si_gtg.ground_agent_id)))

                # Update ground path coord_id for drone
                if si_gtg.ground_agent_id != 0
                    push!(env.ground_paths[si_gtg.ground_agent_id].path_states[si_gtg.vertex_along_path].coord_agent_id, aerial_agent_id)
                end

                # Two cases: transit or transfer
                # Case 1 - flight expanded path between vertices (first or last hop)
                if si_gtg.ground_agent_id != sip_gtg.ground_agent_id

                    if haskey(env.ground_transit_expanded_paths, (si_gtg.idx, sip_gtg.idx))

                        running_dist += env.ground_transit_weights[(si_gtg.idx, sip_gtg.idx)]

                        # Add expanded paths up to last but one
                        exp_path = env.ground_transit_expanded_paths[(si_gtg.idx, sip_gtg.idx)]
                        for idx = 2:length(exp_path)-1
                            running_time += env.road_graph_wts[exp_path[idx-1], exp_path[idx]] / AvgSpeed
                            push!(path_states, AgentPathState(road_vtx_id=exp_path[idx], timeval=running_time))
                        end
                        running_time += env.road_graph_wts[exp_path[end-1], exp_path[end]] / AvgSpeed
                    else
                        # Actually run bidir_eucl for transfer and add path states
                        bidir_sp = compute_bidir_astar_euclidean(env.road_graph, si_gtg.road_vtx_id, sip_gtg.road_vtx_id, env.road_graph_wts, env.location_list)

                        if isempty(bidir_sp.path)

                            # Just add flight path; don't need to add any states
                            flight_dist = distance_lat_lon_euclidean(env.location_list[si_gtg.road_vtx_id], env.location_list[sip_gtg.road_vtx_id])
                            running_dist += flight_dist
                            running_time += flight_dist / AvgSpeed

                        elseif length(bidir_sp.path) > 1

                            running_dist += bidir_sp.dist
                            exp_path = bidir_sp.path
                            for idx = 2:length(exp_path)-1
                                running_time += env.road_graph_wts[exp_path[idx-1], exp_path[idx]] / AvgSpeed
                                push!(path_states, AgentPathState(road_vtx_id=exp_path[idx], timeval=running_time))
                            end
                            running_time += env.road_graph_wts[exp_path[end-1], exp_path[end]] / AvgSpeed

                        end
                    end # haskey(env.ground_transit_expanded_paths, (si_gtg.idx, sip_gtg.idx))

                    # Do some extra work if last is
                    if sip_gtg.ground_agent_id != 0

                        # If running time is more than car, update car
                        # If running time is less than car, update running time
                        if sip_gtg.timeval < running_time

                            ground_route_running_time = running_time

                            # Continue until subsumed by existing time...may go all the way to end
                            for gtg_idx = sip_gtg.idx : env.gtg_ground_path_start_ends[sip_gtg.ground_agent_id][2]-1

                                gtg_vtx = env.ground_transit_graph.vertices[gtg_idx]
                                if gtg_vtx.timeval >= ground_route_running_time
                                    break
                                end

                                # Compute difference to next one BEFORE updating ground route running time
                                time_to_incr = env.ground_transit_graph.vertices[gtg_idx+1].timeval - env.ground_transit_graph.vertices[gtg_idx].timeval

                                # Edit timeval for BOTH GTS and ground paths at vertex along path
                                env.ground_transit_graph.vertices[gtg_idx].timeval = ground_route_running_time
                                env.ground_paths[sip_gtg.ground_agent_id].path_states[gtg_vtx.vertex_along_path].timeval = ground_route_running_time

                                # Update ground route running time by difference
                                ground_route_running_time += time_to_incr

                            end

                            env.ground_transit_graph.vertices[env.gtg_ground_path_start_ends[sip_gtg.ground_agent_id][2]].timeval = ground_route_running_time
                            env.ground_paths[sip_gtg.ground_agent_id].path_states[end].timeval = ground_route_running_time

                        else
                            running_time = sip_gtg.timeval
                        end # if sip_gtg.timeval < running_time
                    end # sip_gtg.ground_agent_id != 0
                else
                    # Transfer, just update running time
                    running_time += sip_gtg.timeval - si_gtg.timeval
                end
            end # i = 1: length(states) - 1

            # Add last state
            push!(path_states, AgentPathState(road_vtx_id = env.aerial_task_list[aerial_agent_id].dest, timeval=running_time))

            final_aerial_paths[aerial_agent_id] = AgentPathInfo(path_states=path_states, total_dist=running_dist, total_time=running_time)

        end # length(states) == or != 2
    end # enumerate(solution)

    env.aerial_paths = final_aerial_paths

    # Update total time of ground paths
    for (i, gp) in enumerate(env.ground_paths)
        env.ground_paths[i].total_time = gp.path_states[end].timeval
    end

end


text:
using .LazySets

function (opl::LinearOettliPrager)(A, b)
    n = length(b)
    Ac = mid.(A)
    bc = mid.(b)
    Ar = IntervalArithmetic.radius.(A)
    br = IntervalArithmetic.radius.(b)

    polytopes = Vector{HPolytope}(undef, 2^n)
    orthants = DiagDirections(n)

    @inbounds for (i, d) in enumerate(orthants)
        D = Diagonal(d)
        Ard = -Ar*D

        A1 = [Ac + Ard; -Ac + Ard; -D]
        b1 = [br + bc; br - bc; zeros(n)]

        polytopes[i] =  HPolytope(A1, b1)
    end
    return identity.(filter!(!isempty, polytopes))
end

_default_precondition(_, ::LinearOettliPrager) = NoPrecondition()


text:
function infer(M::Model, X)
  infer(M.alg, M.arc, X)
end

function update!(M::Model, y, ŷ, x, learn_rate::Float64)
    update!(M, M.alg, M.arc, y, ŷ, x, learn_rate)
end

"""
Linear Algorithm
"""
struct Linear <: Algorithm
    link
end

function Linear()::Linear
  Linear(identity)
end

function infer(alg::Linear, arc::SingleLayer, X)
  return arc.wts * X'
end

# function infer(alg::Linear, arc::MultiLayer, X)
# end

function update!(M::Model, alg::Linear, arc::SingleLayer,
                 y, ŷ, x, learn_rate::Float64)
  M.arc.wts .-= learn_rate * assessgrad(M, y, ŷ, x)
end

# function update!(M::Model, alg::Linear, arc::MultiLayer,
#                  y, ŷ, x, learn_rate::Float64)
# end


text:
module Toxcore

include("gen/tox.h.jl")
include("gen/toxav.h.jl")
include("gen/toxdns.h.jl")

module CInterface

	importall Toxcore
	import Toxcore.Array_128_Uint8

	const libtoxcore = @windows? Pkg.dir("Toxcore", "deps", "bin", "libtoxcore") : @linux? "libtoxcore.so" : error("platform not support")

	include("gen/tox.h.functions.jl")
	include("gen/toxav.h.functions.jl")
	include("gen/toxdns.h.functions.jl")

end # module CInterface

include("util.jl")
include("julia_api.jl")

end # module Toxcore


text:
# The prime 41, can be written as the sum of six consecutive primes:
#
# 41 = 2 + 3 + 5 + 7 + 11 + 13
#
# This is the longest sum of consecutive primes that adds to a prime below
# one-hundred.
#
# The longest sum of consecutive primes below one-thousand that adds to a prime,
# contains 21 terms, and is equal to 953.
#
# Which prime, below one-million, can be written as the sum of the most
# consecutive primes?

using ProjectEulerSolutions

# Get the first n primes that sum to less than n, then start looking at the
# longest to the shortest summed sequences to see which sum is also prime.
function p050solution_iterative(n::Integer=1_000)::Integer

    # Bound the max length to be when the sum of the first primes surpasses n
    primes = Array{Integer}([2, 3])
    running_sum = 2
    upper_ind = 1
    test_val = 3
    while running_sum < n
        upper_ind += 1
        running_sum += primes[upper_ind]
        if length(primes) == upper_ind
            test_val += 1
            while !isprime(test_val)
                test_val += 1
            end
            append!(primes, test_val)
        end
    end

    start_val = sum(primes[1:upper_ind+1])
    for seqlength in upper_ind:-1:1 # sequence lengths to consider
        start_val -= primes[seqlength+1]
        # Even lengths can only work if they include the first prime
        if seqlength % 2 == 0
            if isprime(start_val)
                return start_val
            end
            continue
        end
        prev_val = start_val
        i = 2
        while prev_val < n
            s = prev_val - primes[i-1] + primes[seqlength+i-1]
            if isprime(s)
                return s
            end
            prev_val = s
            i += 1
        end
    end
    return -1
end

# Use seive like before, but build up a cumulative sum so that you only need to
# subtract two numbers to get the sum of primes from index 1 to index j.
function p050solution_cumsum(n::Integer=1_000)::Integer

    primes = Array{Integer}([2, 3])

    cumsum = zeros(Integer, 1)
    cumsum[1] = primes[1]

    test_val = 3
    i = 2
    while cumsum[i-1] <= n
        if length(primes) == i
            test_val += 1
            while !isprime(test_val)
                test_val += 1
            end
            append!(primes, test_val)
        end
        append!(cumsum, cumsum[i-1] + primes[i])
        i += 1

    end

    upper_ind = min(i, length(primes)) - 1

    for seqlength in upper_ind:-1:1 # sequence lengths to consider
        # Even lengths can only work if they include the first prime
        s = cumsum[seqlength]
        if seqlength % 2 == 0
            if isprime(s)
                return s
            end
            continue
        end
        i = 1
        while s < n
            s = cumsum[seqlength + i] - cumsum[i]
            if isprime(s)
                return s
            end
            i += 1
        end
    end

    return -1
end

p050 = Problems.Problem(Dict("Iterative" => p050solution_iterative,
                             "Cumulative sum" => p050solution_cumsum))

Problems.benchmark(p050, 1_000_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LocalAnisotropy(rotation, magnitude)

Create a `LocalAnisotropy` object. The `rotation` is a vector of
`ReferenceFrameRotations.Quaternion` and the magnitude is an array of size
(number of dimensions, number of elements)
"""
struct LocalAnisotropy
  rotation::AbstractVector{Quaternion}
  magnitude::AbstractArray
end

Ints = Union{Int,AbstractVector{Int}}
rotation(L::LocalAnisotropy) = L.rotation
rotation(L::LocalAnisotropy, i::Int) = L.rotation[i]
rotation(L::LocalAnisotropy, i::AbstractVector{Int}) = view(L.rotation, i)
rotmat(L::LocalAnisotropy, i::Int) = quat_to_dcm(L.rotation[i])
magnitude(L::LocalAnisotropy) = L.magnitude
magnitude(L::LocalAnisotropy, i::Ints) = view(L.magnitude, :, i)
ratio1(L::LocalAnisotropy) = L.magnitude[2,:] ./ L.magnitude[1,:]
ratio1(L::LocalAnisotropy, i::Ints) = L.magnitude[2,i] / L.magnitude[1,i]
ratio2(L::LocalAnisotropy) = L.magnitude[3,:] ./ L.magnitude[1,:]
ratio2(L::LocalAnisotropy, i::Ints) = L.magnitude[3,i] / L.magnitude[1,i]
slice(L::LocalAnisotropy, i::Ints) = LocalAnisotropy(rotation(L,i),magnitude(L,i))
nvals(L::LocalAnisotropy) = length(L.rotation)
ndims(L::LocalAnisotropy) = size(L.magnitude,1)
iaxis(ax) = (X=1, Y=2, Z=3)[ax]

Base.show(io::IO, lp::LocalAnisotropy) = print(io, "LocalAnisotropy")

function Base.show(io::IO, ::MIME"text/plain", lp::LocalAnisotropy)
	print(io,"LocalAnisotropy $(ndims(lp))-D")
end

SpatialData = Union{GeoData,Domain,DataView,DomainView}

struct LocalGeoData
	data::Union{SpatialData,Nothing}
	domain::SpatialData
	localaniso::LocalAnisotropy
	refvario::Union{Variogram,Nothing}
	datapars::Union{AbstractVector{Int},Nothing}
    graph::Union{SimpleWeightedGraph,Nothing}
end

"""
    LocalGeoData(domain, localaniso)

Association of a spatial domain and its `LocalAnisotropy`.
"""
function LocalGeoData(obj::SpatialData, lpars::LocalAnisotropy)
	LocalGeoData(nothing, obj, lpars, nothing, nothing, nothing)
end

"""
    LocalGeoData(domain, localaniso, refvario)

Association of a spatial domain, its `LocalAnisotropy` and its reference variogram.
"""
function LocalGeoData(obj::SpatialData, lpars::LocalAnisotropy, refvario::Variogram)
	LocalGeoData(nothing, obj, lpars, refvario, nothing, nothing)
end

"""
    LocalGeoData(samples, domain, localaniso)

Association of samples, a spatial domain and their `LocalAnisotropy`. The
samples `LocalAnisotropy` are passed from the spatial data via nearest neighbors.
"""
function LocalGeoData(hd::SpatialData, obj::SpatialData, lpars::LocalAnisotropy)
	hdids = grid2hd_ids(hd,obj)
	LocalGeoData(hd, obj, lpars, nothing, hdids, nothing)
end

"""
    LocalGeoData(samples, domain, localaniso, refvario)

Association of samples, a spatial domain, their `LocalAnisotropy` and the
reference variogram. The samples `LocalAnisotropy` are passed from the spatial data
via nearest neighbors.
"""
function LocalGeoData(hd::SpatialData, obj::SpatialData, lpars::LocalAnisotropy, refvario::Variogram)
	hdids = grid2hd_ids(hd,obj)
	LocalGeoData(hd, obj, lpars, refvario, hdids, nothing)
end

obj(D::LocalGeoData) = D.domain
sobj(D::LocalGeoData) = D.data
ndims(D::LocalGeoData) = embeddim(D.domain)
sdata(D::LocalGeoData) = sobj(D) != nothing
nvals(D::LocalGeoData) = nelements(D.domain)
snvals(D::LocalGeoData) = nelements(D.data)
spars(D::LocalGeoData) = D.datapars
spars(D::LocalGeoData, i::Int) = spars(D)[i-nvals(D)]
spars(D::LocalGeoData, i::AbstractVector{Int}) = view(spars(D), i .- nvals(D))
nall(D::LocalGeoData) = sdata(D) ? nvals(D)+snvals(D) : nvals(D)
centro(D::LocalGeoData, i::Int) = i<=nvals(D) ? centroid(obj(D),i) : centroid(sobj(D),i-nvals(D))
coords(D::LocalGeoData, i::Int) = coordinates(centro(D,i))
coords(D::SpatialData) = reduce(hcat, [coordinates(centroid(D,x)) for x in 1:nelements(D)])
coords(D::SpatialData, i::AbstractVector{Int}) = reduce(hcat, [coordinates(centroid(D,x)) for x in i])

rotation(D::LocalGeoData) = D.localaniso.rotation
srotation(D::LocalGeoData) = view(rotation(D),spars(D))
magnitude(D::LocalGeoData) = D.localaniso.magnitude
smagnitude(D::LocalGeoData) = view(magnitude(D),:,spars(D))

rotation(D::LocalGeoData, i::Int) = i<=nvals(D) ? rotation(D)[i] : rotation(D)[spars(D,i)]
rotation(D::LocalGeoData, i::AbstractVector{Int}) = i<=nvals(D) ? view(rotation(D), i) : view(rotation(D), spars(D,i))
magnitude(D::LocalGeoData, i::Ints) = i<=nvals(D) ? view(magnitude(D), :, i) : view(magnitude(D), :, spars(D,i))

Base.show(io::IO, ld::LocalGeoData) = print(io, "LocalGeoData")

function Base.show(io::IO, ::MIME"text/plain", ld::LocalGeoData)
	print(io,"LocalGeoData $(ndims(ld))-D")
end

abstract type LocalParMethods end
struct Geometric <: LocalParMethods end
struct Gradients <: LocalParMethods end

abstract type LocalMetric end
struct AnisoDistance  <: LocalMetric end
struct KernelVariogram <: LocalMetric end
struct GraphDistance  <: LocalMetric end


text:
module WAdag
using DataFrames, XLSX
import Unicode: normalize

export cleanaddress
export bfunc

include("EmployeeVendorAddress.jl")
include("basicfunc.jl")

end # module


text:
module Trapla

using QuadGK

export  smoothjump, pulsebell
include("gevrey.jl")

end # module


text:
using Pkg; haskey(Pkg.installed(),"CodecZlib") || Pkg.add("CodecZlib")
using CodecZlib

"Where to download mnist from"
mnisturl = "http://yann.lecun.com/exdb/mnist"

"Where to download mnist to"
mnistdir = joinpath(@__DIR__, "mnist")

"""

This utility loads the [MNIST](http://yann.lecun.com/exdb/mnist)
hand-written digits dataset.  There are 60000 training and 10000 test
examples. Each input x consists of 784 pixels representing a 28x28
grayscale image.  The pixel values are converted to Float32 and
normalized to [0,1].  Each output y is a UInt8 indicating the correct
class.  10 is used to represent the digit 0.

```
# Usage:
include("mnist.jl")
xtrn, ytrn, xtst, ytst = mnist()
# xtrn: 28×28×1×60000 Array{Float32,4}
# ytrn: 60000-element Array{UInt8,1}
# xtst: 28×28×1×10000 Array{Float32,4}
# ytst: 10000-element Array{UInt8,1}
```

"""
function mnist()
    global _mnist_xtrn,_mnist_ytrn,_mnist_xtst,_mnist_ytst
    if !@isdefined(_mnist_xtrn)
        @info("Loading MNIST...")
        _mnist_xtrn = _mnist_xdata("train-images-idx3-ubyte.gz")
        _mnist_xtst = _mnist_xdata("t10k-images-idx3-ubyte.gz")
        _mnist_ytrn = _mnist_ydata("train-labels-idx1-ubyte.gz")
        _mnist_ytst = _mnist_ydata("t10k-labels-idx1-ubyte.gz")
    end
    return _mnist_xtrn,_mnist_ytrn,_mnist_xtst,_mnist_ytst
end

"Utility to view a MNIST image, requires Images."
mnistview(x,i)=colorview(Gray,permutedims(x[:,:,1,i],(2,1)))

"Return minibatched mnist data, requires Knet."
function mnistdata(;batchsize=100,xtype=(gpu()>=0 ? KnetArray{Float32} : Array{Float32}), o...)
    if !@isdefined(_mnist_xtrn); mnist(); end
    minibatch(_mnist_xtrn, _mnist_ytrn, batchsize; xtype=xtype, o...),
    minibatch(_mnist_xtst, _mnist_ytst, batchsize; xtype=xtype, o...)
end

function _mnist_xdata(file)
    a = _mnist_gzload(file)[17:end]
    reshape(a ./ 255f0, (28,28,1,div(length(a),784)))
end

function _mnist_ydata(file)
    a = _mnist_gzload(file)[9:end]
    a[a.==0] .= 10
    # full(sparse(a,1:length(a),1f0,10,length(a)))
    return a
end

function _mnist_gzload(file)
    if !isdir(mnistdir)
        mkpath(mnistdir)
    end
    path = joinpath(mnistdir,file)
    if !isfile(path)
        url = "$mnisturl/$file"
        download(url, path)
    end
    f = GzipDecompressorStream(open(path))
    a = read(f)
    close(f)
    return(a)
end

nothing


text:
"""
    ForwardSensFspOutputSparse{NS,IntT<:Integer,RealT<:AbstractFloat}

Struct to store Forward Sensitivity Finite State Projection outputs based on the sparse representation of the FSP solution.

# Fields
- `t`: Array of solution output times.
- `p`: Array of `FspVectorSparse` instances. `p[i]` is the CME solution at time `t[i]`.
- 'S': Array of arrays of `FspVectorSparse` instances. `S[j][i]` is the `j`-th partial derivative of `p[i]`.
- `sinks`: Probability mass accumulated at the sink states over time. `sinks[i]` is the vector of sink state probabilities at time `t[i]`.
- `dsinks`: Partial derivatives of sink probabilties. `dsinks[j][i]` is the `j`-th partial derivative of `sinks[i]`.

# Usage
If `sol` is of type `FspOutputSparse`, `sol[i]` will return a slice, of type `FspOutputSliceSparse` of the solution set at the `i`-th index.

# See also
[`ForwardSensFspOutputSliceSparse`](@ref).

"""
Base.@kwdef mutable struct ForwardSensFspOutputSparse{NS,IntT<:Integer,RealT<:AbstractFloat}
    t::Vector{RealT}
    p::Vector{FspVectorSparse{NS,IntT,RealT}}
    sinks::Vector{AbstractVector{RealT}}
    S::Vector{Vector{FspVectorSparse{NS,IntT,RealT}}}
    dsinks::Vector{Vector{AbstractVector{RealT}}}
end

"""
    ForwardSensFspOutputSliceSparse

Struct to store the Forward Sensitivity FSP solution at a single time. Fields: `t`, `p`, `S`, `sinks`, `dsinks`.
"""
mutable struct ForwardSensFspOutputSliceSparse{NS,IntT<:Integer,RealT<:AbstractFloat}
    t::RealT
    p::FspVectorSparse{NS,IntT,RealT}
    sinks::AbstractVector{RealT}
    S::Vector{FspVectorSparse{NS,IntT,RealT}}
    dsinks::Vector{AbstractVector{RealT}}
end

function Base.getindex(sfspoutput::ForwardSensFspOutputSparse, ind::Integer)
    if ind > length(sfspoutput.t)
        throw(ArgumentError("Requested index exceeds array limit."))
    end
    return ForwardSensFspOutputSliceSparse(
        sfspoutput.t[ind],
        sfspoutput.p[ind],
        sfspoutput.sinks[ind],
        sfspoutput.S[ind],
        sfspoutput.dsinks[ind]
    )
    return
end

function Base.getindex(sensfspoutput::ForwardSensFspOutputSparse, inds::Vector{T}) where {T <: Integer}
    return [sensfspoutput[i] for i in inds]
end

function Base.lastindex(sensfspoutput::ForwardSensFspOutputSparse)
    return length(sensfspoutput.t)
end

function Base.length(sfspoutput::ForwardSensFspOutputSparse)
    return length(sfspoutput.t)
end
function Base.size(sfspoutput::ForwardSensFspOutputSparse)
    return size(sfspoutput.t)
end

using Printf
function Base.show(io::IO, output::ForwardSensFspOutputSparse)
    println(io, "Forward Sensitivity FSP solution set with $(length(output)) snapshots")
end

function Base.show(io::IO, outputslice::ForwardSensFspOutputSliceSparse)
    println(io, "Forward Sensitivity FSP solution output slice:")
    @printf(io, "t = %.2e \n", outputslice.t)
    println(io, "sinks = $(outputslice.sinks)")
    @printf(io, "p = \n")
    show(io, outputslice.p)
    @printf(io, "S = \n")
    for svec in outputslice.S
        show(io, svec)
    end
end