text:
module GroundDataV
"""
Module for the download and processing of atmospheric data gathered by measuring stations located in Veneto, Italy.
"""

using CombinedParsers
using CombinedParsers.Regexp
using CSV
using DataFrames
using Dates
using HTTP
using JSON
using JSONTables
using Revise

export getData,
       getRegionAttributes, getRegionIds, getRegionStationInfo

@syntax informations = Sequence(
    :info => Sequence(
        Sequence( "<", "PERIODO", ">" ),
        "<![CDATA[",
        re"[^\[\]]+",
        "]]>",
        "</PERIODO>"
    ),
    :begin => Sequence(
        Sequence( "<", "INIZIO", ">" ),
        Numeric(Int64),
        "</INIZIO>"
    ),
    :end => Sequence(
        Sequence( "<", "FINE", ">" ),
        Numeric(Int64),
        "</FINE>"
    ),
    :proj => Sequence(
        Sequence( "<", "PROJECTION", ">" ),
        re"[^<:>]+",
        ":",
        Numeric(Int64),
        "</PROJECTION>"
    )
)

@syntax stations = Repeat(
    Sequence(
        re"<[^<>]+>",
        Either(
            Numeric(Int64),
            Numeric(Float64),
            Sequence( "<![CDATA[", re"[^\[\]]+", "]]>" ),
            re"[^<>]+"
        ),
        re"</[^<>]+>"
    )
)

@syntax sensor = Sequence(
    Repeat(
        re"<(?!D)[^>]+>",
        Either(
            Numeric(Int64),
            Sequence( "<![CDATA[", re"[^\[\]]+", "]]>" ),
            re"[^>]+"
        ),
        re"</[^>]+>"
    ),
    Repeat(
        Sequence( "<DATI ISTANTE=\"", Numeric(Int64), "\">" ),
        Sequence( "<VM>", Numeric(Float64), "</VM>" ),
        "</DATI>"
    )
)

"""
    getRegionAttributes( [ type::Symbol=:METEO ] )

Obtain the names of the columns of the region's dataframe required by `GroundData.createMap`'s `attributes` parameter to create `GroundData.standardize`'s `map` parameter.
"""
function getRegionAttributes( type::Symbol=:METEO )
    return type == :METEO ? [ :paramnm, :unitnm, :value, nothing, :instant, :x, :y, :quota, nothing, nothing, :rmh ] :
        type == :AIRQUALITY ? [ :param, nothing, :value, nothing, :date, :lon, :lat, nothing, nothing, :tipozona ] :
            throw(DomainError(type, "`type` must be either `:METEO` OR `:AIRQUALITY`."))
end

"""
    getRegionAttributes( [ type::Symbol=:METEO ] )

Obtain the names of the columns of the dataframe required for `GroundData.standardize`'s `bridge` parameter.
"""
function getRegionIds( type::Symbol=:METEO )
    return type == :METEO ? :idstaz :
        type == :AIRQUALITY ? :codseqst : throw(DomainError(type, "`type` must be either `:METEO` OR `:AIRQUALITY`."))
end

"""
    getRegionStationInfo( [ type::Symbol=:METEO  ] )

Obtain the names of the columns of the region's stations dataframe required by `GroundData.createMap`'s `attributes` parameter to be used in `GroundData.generateUuidsTable`.
"""
function getRegionStationsInfo( type::Symbol=:METEO )
    return type == :METEO ? [ :idstaz, :nome, :x, :y ] :
        type == :AIRQUALITY ? [ :codseqst, :nome, :lon, :lat ] :
            throw(DomainError(type, "`type` must be either `:METEO` OR `:AIRQUALITY`."))
end

"""
    getStationsInfo()

Obtain the dataframe containing informations on all available ARPAV stations.
"""
function getMeteoStationsData()
 #Get the information on the available stations
    # Get the page containing informations on all the available stations
    page = String( HTTP.get( "https://www.arpa.veneto.it/bollettini/meteo/h24/img07/stazioni.xml" ).body )
    # Keep only the useful portion of the page
    useful = split( page, "</PERIODO>" )[2][1:end-15]
    # Split the downloaded string in substrings one for each station
    arr = split( useful, r"</?STAZIONE>", keepempty=false )
    #info = informations(arr[1])
    stats = stations.(arr[2:end])
    #Create the dictionary of the stations, checking the type of `attribute[2]` (contains the value of the attribute)
    df = DataFrame([
        Dict(
            Symbol( lowercase( String( attribute[1][2] ) ) )
            =>
            isa(attribute[2], Number) ? attribute[2] : # If it's a number leave it as is
                isa(attribute[2], Array) ? lowercase( String( attribute[2] ) ) : #If it is an array (of strings) convert it to String
                    titlecase( String( attribute[2][2] ) ) # If it is a Tuple ( Es. `("<![CDATA[", "Arabba", "]]>")` ) take the actual value ( `"Arabba"` )
            for attribute in station
        ) for station in stats
    ])
    insertcols!(df, :rmh => "0m")
    return df
end

# res = getMeteoStationData()

"""
    getStationsData( stats::AbstractVector{String} )

Obtain the dataframe containing the data of all the stations described by the elements of `stats`.
"""
function getMeteoData( ids::AbstractVector{Int64} )
    pages_vect = [ String( HTTP.get("https://www.arpa.veneto.it/bollettini/meteo/h24/img07/$(lpad(id, 4, "0")).xml").body ) for id in ids ]
    # For each of the station remove the first part of the string and the closing tags
    stat_strings_vect = [ split(page, "</ATTIVAZIONE>")[2][1:end-26] for page in pages_vect ]
    # From each station generate the corresponding array of sensors
    sensor_vect = split.( stat_strings_vect, r"</?SENSORE>", keepempty=false )
    # vect[i][j] i-th station, j-th sensor
    # vect[i][j][1] j-th sensor's info (Vector)
    # vect[i][j][2] j-th sensor's data (Vector)
    # vect[i][j][2][l] j-th sensor'data, l-th measurement (Vector)
    # vect[i][j][2][l][4] l-th measurement's values (Vector 1 => mean, 2 => min, 3 => max)
    vect = [ sensor.(sensor_group) for sensor_group in sensor_vect ]
    df = DataFrame([
        push!(
            Dict(
                Symbol( lowercase( String( attribute[3] ) ) )
                =>
                isa( attribute[5], Number ) ? attribute[5] :
                    isa( attribute[5], Array ) ? String( attribute[5] ) :
                        titlecase( String( attribute[5][2] ) )
              for attribute in sensor[1]
            ),
            :idstaz => id,
            :instant => entry[2],
            :value => entry[5]
        )
        for (id, station) in zip(ids, vect)
            for sensor in station
                for entry in sensor[2]
    ])
    transform!(
        df,
        [:unitnm] => ByRow( x -> x = replace( replace( x, "\xb0" => "°" ), "2" => "²" ) ) => :unitnm,
        [:instant] => ByRow( x -> DateTime( string(x), "yyyymmddHHMM" ) ) => :instant
    )
    return df
end

#   df = getMeteoStationsData()
#   data = getMeteoSensorsInfo( df[1:3, :idstaz] )

"""
    getAqStationsData()

Obtain a dataframe containing informations on the measuring stations in the area.
"""
function getAqStationsData()
    page = String(HTTP.get("http://213.217.132.81/aria-json/exported/aria/stats.json").body)
    df = DataFrame(jsontable(page[14:end]))
    return df
end

#   res = getAqStationsData()

"""
    getAqData()

Obtain the data gathered by the stations in the area.
"""
function getAqData()
    page = String(HTTP.get("http://213.217.132.81/aria-json/exported/aria/data.json").body)
    # js["stazioni"][i]                                            i-th station
    # js["stazioni"][i]["codseqst"]                                i-th station's id
    # js["stazioni"][i]["misurazioni"][j]                          i-th station's j-th measurement
    # js["stazioni"][i]["misurazioni"][j]["pm10"/"ozono"][l]       j-th measurement's l-th entry
    js = JSON.parse(page)["stazioni"]
    arr = Dict[]
    for station in js
        if !isempty(station["misurazioni"])
            for measurement in station["misurazioni"]
                for key in collect(keys(measurement))
                    for entry in measurement[key]
                        push!(
                            arr,
                            Dict(
                                :codseqst => station["codseqst"],
                                :param => key,
                                :date => entry["data"],
                                :value => entry["mis"]
                            )
                        )
                    end
                end
            end
        else
            push!( arr, Dict( :codseqst => station["codseqst"], :param => missing, :date => missing, :value => missing ) )
        end
    end
    return DataFrame(arr)
end

#   df = getAqSensorsData()

"""
    getData(; <keyword arguments> )

Obtain data of category `type` and source of category `kind`.

# Arguments
- `type::Symbol=:METEO`: defines the type of data to be downloaded, it must either be `:METEO` or `:AIRQUALITY`.
- `kind::Symbol=:STATIONS`: defines if the data to be downloaded has to regard the stations or their actual measurements, it must either be `:STATIONS` or `:SENSORS`.
"""
function getData(; type::Symbol=:METEO, kind::Symbol=:STATIONS )
    if type == :METEO
        stations = getMeteoStationsData()
        if kind == :STATIONS
            return stations
        elseif kind == :SENSORS
            return getMeteoData( stations[:, :idstaz] )
        else
            throw(DomainError(kind, "`kind` must either be `:STATIONS` or `:SENSORS`."))
        end
    elseif type == :AIRQUALITY
        if kind == :STATIONS
            return getAqStationsData()
        elseif kind == :SENSORS
            return getAqData()
        else
            throw(DomainError(kind, "`kind` must either be `:STATIONS` or `:SENSORS`."))
        end
    else
        throw(DomainError(type, "`type` must either be `:METEO` or `:AIRQUALITY`."))
    end
end

#   ressta = getData()
#   ressen = getData( source=:SENSORS )
#   ressta = getData( type=:AIRQUALITY, source=:STATIONS )
#   ressen = getData( type=:AIRQUALITY, source=:SENSORS )

end # moduled

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
### Right-hand side and solution vectors

const SaddleVector = ArrayPartition

#SaddleVector(u::TU,f::TF) where {TU,TF} = ArrayPartition(u,f)

"""
    SaddleVector(u,f)

Construct a vector of a state part `u` and constraint part `f` of a
saddle-point vector, to be associated with a [`SaddleSystem`](@ref).
"""
function SaddleVector end

"""
    state(x::SaddleVector)

Provide the state part of the given saddle vector `x`
"""
state(u::SaddleVector) = u.x[1]

"""
    constraint(x::SaddleVector)

Provide the constraint part of the given saddle vector `x`
"""
constraint(u::SaddleVector) = u.x[2]

ext: jl


text:
# indicator of a point

export IndPoint

"""
**Indicator of a singleton**

    IndPoint(p=0.0)

Returns the indicator of the set
```math
C = \\{p \\}.
```
Parameter `p` can be a scalar, in which case the unique element of `S` has uniform coefficients.
"""
struct IndPoint{R <: Real, C <: Union{R, Complex{R}}, T <: Union{C, AbstractArray{C}}} <: ProximableFunction
    p::T
    function IndPoint{R, C, T}(p::T) where {R, C, T}
        new(p)
    end
end

is_separable(f::IndPoint) = true
is_convex(f::IndPoint) = true
is_singleton(f::IndPoint) = true
is_cone(f::IndPoint) = norm(f.p) == 0
is_affine(f::IndPoint) = true

IndPoint(p::T=0.0) where {R <: Real, C <: Union{R, Complex{R}}, T <: Union{C, AbstractArray{C}}} = IndPoint{R, C, T}(p)

function (f::IndPoint{R, C})(x::AbstractArray{C}) where {R, C}
    if all(x .≈ f.p)
        return R(0)
    end
    return R(Inf)
end

function prox!(y::AbstractArray{C}, f::IndPoint{R, C}, x::AbstractArray{C}, gamma=R(1)) where {R, C}
    y .= f.p
    return R(0)
end

fun_name(f::IndPoint) = "indicator of a point"
fun_dom(f::IndPoint) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndPoint) = "x ↦ 0 if x = p, +∞ otherwise"
fun_params(f::IndPoint) =
    string( "p = ", typeof(f.p) <: AbstractArray ? string(typeof(f.p), " of size ", size(f.p)) : f.p, ", ")

function prox_naive(f::IndPoint{R, C}, x::AbstractArray{C}, gamma=R(1)) where {R, C}
    y = similar(x)
    y .= f.p
    return y, R(0)
end


text:
# EXCLUDE FROM TESTING
# NOTE: This is temporarily not working until we merge the DMSTAG routines with the new Clang branch
#
#
# This is an example of a 1D viscoelastic porosity wave as described in
# Vasyliev et al. Geophysical Research Letters (25), 17. p. 3239-3242
# https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/98GL52358
#
# It simulates how a pulse of magma migrates upwards in the Earth, which
# can be described by a set of coupled nonlinear PDE's
#
# This example only requires the specification of a residual routine; automatic
# differentiation is used to generate the jacobian.

using PETSc, MPI

using SparseArrays, SparseDiffTools, ForwardDiff

PETSc.initialize()

CreatePlots = false;
if CreatePlots==true
    using Plots
end

function FormRes!(cfx_g, cx_g, user_ctx)

    # Note that in PETSc, cx and cfx are pointers to global vectors.

    # Copy global to local vectors
    PETSc.DMGlobalToLocal(user_ctx.dm, cx_g,  PETSc.INSERT_VALUES,  user_ctx.x_l)
    PETSc.DMGlobalToLocal(user_ctx.dm, cfx_g, PETSc.INSERT_VALUES,  user_ctx.f_l)

    # Retrieve arrays from the local vectors
    ArrayLocal_x     =   PETSc.DMStagVecGetArrayRead(user_ctx.dm,   user_ctx.x_l);      # array with all local x-data
    ArrayLocal_f     =   PETSc.DMStagVecGetArray(user_ctx.dm,       user_ctx.f_l);      # array with all local residual

    # Compute local residual
    ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx)

    # Finalize local arrays
    Base.finalize(ArrayLocal_x)
    Base.finalize(ArrayLocal_f)

    # Copy local into global residual vector
    PETSc.DMLocalToGlobal(user_ctx.dm,user_ctx.f_l, PETSc.INSERT_VALUES, cfx_g)

end

function  ForwardDiff_res(x, user_ctx)
    f   = zero(x)               # vector of zeros, of same type as x (local vector)

    ArrayLocal_x     =   PETSc.DMStagVecGetArray(user_ctx.dm, x);        # array with all local x-data
    ArrayLocal_f     =   PETSc.DMStagVecGetArray(user_ctx.dm, f);        # array with all local residual

    ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx);

    # As the residual vector f is linked with ArrayLocal_f, we don't need to
    # pass ArrayLocal_f back to f

    return f;
end

function  f(out, x, user_ctx)

    ArrayLocal_x     =   PETSc.DMStagVecGetArray(user_ctx.dm, x);        # array with all local x-data
    ArrayLocal_f     =   PETSc.DMStagVecGetArray(user_ctx.dm, out);        # array with all local residual

    ComputeLocalResidual(user_ctx.dm, ArrayLocal_x, ArrayLocal_f, user_ctx);

    return nothing
end

function FormJacobian!(cx_g, J, P, user_ctx)
    # This requires several steps:
    #
    #   1) Extract local vector from global solution (x) vector
    #   2) Compute local jacobian from the residual routine (note that
    #       this routine requires julia vectors as input)

    # Extract the local vector
    PETSc.DMGlobalToLocal(user_ctx.dm, cx_g,  PETSc.INSERT_VALUES,  user_ctx.x_l)
    x               =   PETSc.unsafe_localarray(Float64, user_ctx.x_l.ptr;  write=false, read=true)

    if isnothing(user_ctx.jac)
        # Compute sparsity pattern of jacobian. This is relatvely slow, but only has to be done once.
        # Theoretically, more efficient tools for this exists (jacobian_sparsity in the SparsityDetection.jl package),
        # but they don't seem to work with the PETSc approach we use. Therefore we employ
        f_Residual  =   (x -> ForwardDiff_res(x, user_ctx));        # pass additional arguments into the routine
        J_julia     =   ForwardDiff.jacobian(f_Residual,x);

        # employ sparse structure to compute jacobian - to be moved inside routines
        jac         =   sparse(J_julia);
        colors      =   matrix_colors(jac)          # the number of nonzeros per row

    else
        jac     =   user_ctx.jac;
        colors  =   user_ctx.colors;

    end
    out         =   similar(x);

   f_Res           =   ((out,x)->f(out, x, user_ctx));        # pass additional arguments into the routine
   forwarddiff_color_jacobian!(jac, f_Res, x, colorvec = colors)

    ind             =   PETSc.LocalInGlobalIndices(user_ctx.dm);    # extr
    J              .=   jac[ind,ind];

    user_ctx.jac    =   jac;
    user_ctx.colors =   colors;

   return jac, ind
end

# Define a struct that holds data we need in the local SNES routines below
mutable struct Data_PorWav1D
    dm
    x_l
    xold_l
    xold_g
    f_l
    dt
    dz
    De
    jac
    colors
end
user_ctx = Data_PorWav1D(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing);  # holds data we need in the local

function ComputeLocalResidual(dm, ArrayLocal_x, ArrayLocal_f, user_ctx)
    # Compute the local residual. The vectors include ghost points
    n           =   3.0
    m           =   2.0
    dt          =   user_ctx.dt;
    dz          =   user_ctx.dz;
    De          =   user_ctx.De;

    Phi         =   PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x   ,   PETSc.DMSTAG_ELEMENT,   0);
    Phi_old     =   PETSc.DMStagGetGhostArrayLocationSlot(dm,user_ctx.xold_l,   PETSc.DMSTAG_ELEMENT,   0);
    Pe          =   PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_x   ,   PETSc.DMSTAG_ELEMENT,   1);
    Pe_old      =   PETSc.DMStagGetGhostArrayLocationSlot(dm,user_ctx.xold_l,   PETSc.DMSTAG_ELEMENT,   1);

    res_Phi     =   PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f,      PETSc.DMSTAG_ELEMENT,   0);
    res_Pe      =   PETSc.DMStagGetGhostArrayLocationSlot(dm,ArrayLocal_f,      PETSc.DMSTAG_ELEMENT,   1);

    # compute the FD stencil
    ind         =     PETSc.DMStagGetCentralNodes(dm);          # indices of (center/element) points, not including ghost values.

    # Porosity residual @ center points
    iz                  =   ind.Center.lower[1]+1:ind.Center.upper[1]-1;            # Phi is on center points
    res_Phi[iz[1]-1]    =   Phi[iz[1]-1]   - 1.0;       # Bottom BC
    res_Phi[iz[end]+1]  =   Phi[iz[end]+1] - 1.0;       # Top BC
    res_Phi[iz]         =   (Phi[iz] - Phi_old[iz])/dt + De.*(Pe[iz]-Pe_old[iz])/dt + (Phi[iz].^m)   .* Pe[iz]

    # Pressure update @ nodal points
    iz                  =   ind.Vertex.lower[1]+1:ind.Vertex.upper[1]-1;
    #iz                  =   sx[1]+1:sn[1]-1;            # Pe is on center points as well (dof=2)
    res_Pe[iz[1]-1]     =   Pe[iz[1]-1]   - 0.;         # Bottom BC
    res_Pe[iz[end]+1]   =   Pe[iz[end]+1] - 0.;         # Top BC
    res_Pe[iz]          =   De.*(Pe[iz]-Pe_old[iz])/dt -   ( ((0.5*(Phi[iz .+ 1] + Phi[iz .+ 0])).^n) .* ( (Pe[iz .+ 1] - Pe[iz     ])/dz .+ 1.0)  -
                                                             ((0.5*(Phi[iz .- 1] + Phi[iz .+ 0])).^n) .* ( (Pe[iz     ] - Pe[iz .- 1])/dz .+ 1.0))/dz  +
                                                             (Phi[iz].^m)   .* Pe[iz];

    # Cleanup
    Base.finalize(Phi);    Base.finalize(Phi_old);
    Base.finalize(Pe);     Base.finalize(Pe_old);

end

function SetInitialPerturbations(user_ctx, x_g)
    # Computes the initial perturbations as in the paper

    # Retrieve coordinates from DMStag
    DMcoord     =   PETSc.DMGetCoordinateDM(user_ctx.dm)
    vec_coord   =   PETSc.DMGetCoordinatesLocal(user_ctx.dm);
    Coord       =   PETSc.DMStagVecGetArray(DMcoord, vec_coord);
    Z_cen       =   PETSc.DMStagGetGhostArrayLocationSlot(user_ctx.dm,Coord, PETSc.DMSTAG_ELEMENT,    0); # center (has 1 extra)
    user_ctx.dz =   Z_cen[2]-Z_cen[1];

    Phi         =   PETSc.DMStagGetGhostArrayLocationSlot(user_ctx.dm,user_ctx.x_l,     PETSc.DMSTAG_ELEMENT, 0);
    Pe          =   PETSc.DMStagGetGhostArrayLocationSlot(user_ctx.dm,user_ctx.x_l,     PETSc.DMSTAG_ELEMENT, 1);
    Phi_old     =   PETSc.DMStagGetGhostArrayLocationSlot(user_ctx.dm,user_ctx.xold_l,  PETSc.DMSTAG_ELEMENT, 0);
    Pe_old      =   PETSc.DMStagGetGhostArrayLocationSlot(user_ctx.dm,user_ctx.xold_l,  PETSc.DMSTAG_ELEMENT, 1);

    Phi0 =1.0; dPhi1=8.0; dPhi2=1.0; z1=0.0; z2=40.0; lambda=1.0
    dPe1        =   dPhi1/user_ctx.De;
    dPe2        =   dPhi2/user_ctx.De;

    Phi        .=   Phi0 .+ dPhi1.*exp.( -((Z_cen .- z1).^2.0)/lambda^2) +  dPhi2.*exp.( -((Z_cen .- z2).^2.0)/lambda^2);
    Pe         .=          -dPe1 .*exp.( -((Z_cen .- z1).^2.0)/lambda^2) -   dPe2.*exp.( -((Z_cen .- z2).^2.0)/lambda^2);

    Phi_old    .=   Phi0 .+ dPhi1.*exp.( -((Z_cen .- z1).^2.0)/lambda^2) +  dPhi2.*exp.( -((Z_cen .- z2).^2.0)/lambda^2);
    Pe_old     .=          -dPe1 .*exp.( -((Z_cen .- z1).^2.0)/lambda^2) -   dPe2.*exp.( -((Z_cen .- z2).^2.0)/lambda^2);

    # Copy local into global residual vector
    PETSc.DMLocalToGlobal(user_ctx.dm,user_ctx.x_l, PETSc.INSERT_VALUES, x_g)

    # send back coordinates
    return Z_cen
end

# Main Solver
nx              =   1001;
L               =   150;
#user_ctx.De     =   1e-2;        # Deborah number
#user_ctx.dt     =   5e-5;        # Note that the timestep has to be tuned a bit depending on De in order to obtain convergence
user_ctx.De     =   1e2;
user_ctx.dt     =   1e-2;        # Note that the timestep has to be tuned a bit depending on De in order to obtain convergence

user_ctx.dm     =   PETSc.DMStagCreate1d(MPI.COMM_SELF,PETSc.DM_BOUNDARY_NONE,nx,0,2);  # both Phi and Pe are on center points
PETSc.DMStagSetUniformCoordinatesExplicit(user_ctx.dm, -20, L)            # set coordinates
x_g             =   PETSc.DMCreateGlobalVector(user_ctx.dm)

f_g             =   PETSc.DMCreateGlobalVector(user_ctx.dm)
user_ctx.x_l    =   PETSc.DMCreateLocalVector(user_ctx.dm)
user_ctx.xold_l =   PETSc.DMCreateLocalVector(user_ctx.dm)
user_ctx.xold_g =   PETSc.DMCreateGlobalVector(user_ctx.dm)
user_ctx.f_l    =   PETSc.DMCreateLocalVector(user_ctx.dm)
J               =   PETSc.DMCreateMatrix(user_ctx.dm);                  # Jacobian from DMStag

# initial non-zero structure of jacobian
Z_cen           =   SetInitialPerturbations(user_ctx, x_g)

x0              =   PETSc.VecSeq(rand(size(x_g,1)));
J_julia,ind     =   FormJacobian!(x0.ptr, J, J, user_ctx)

S = PETSc.SNES{Float64}(MPI.COMM_SELF;
        snes_rtol=1e-12,
        snes_monitor=true,
        snes_max_it = 500,
        snes_monitor_true_residual=true,
        snes_converged_reason=true);
S.user_ctx  =       user_ctx;

SetInitialPerturbations(user_ctx, x_g)

PETSc.setfunction!(S, FormRes!, f_g)
PETSc.setjacobian!(S, FormJacobian!, J, J)

# Preparation of visualisation
if CreatePlots
    ENV["GKSwstype"]="nul";
    if isdir("viz_out")==true
        rm("viz_out", recursive=true)
    end
    mkdir("viz_out")
    loadpath = "./viz_out/"; anim = Animation(loadpath,String[])
end

time = 0.0;
it   = 1;
while time<2.5
    global time, Z, Z_cen, it

    # Solve one (nonlinear) timestep
    PETSc.solve!(x_g, S);

    # Update old local values
    user_ctx.xold_g  =  x_g;
    PETSc.DMGlobalToLocal(user_ctx.dm, x_g,  PETSc.INSERT_VALUES,  user_ctx.x_l)
    PETSc.DMGlobalToLocal(user_ctx.dm, user_ctx.xold_g,  PETSc.INSERT_VALUES,  user_ctx.xold_l)

    # Update time
    time += user_ctx.dt;
    it   += 1;

    if (mod(it,200)==0) & (CreatePlots==true) # Visualisation
        # Extract values and plot
        Phi         =   PETSc.DMStagGetGhostArrayLocationSlot(user_ctx.dm,user_ctx.x_l,     PETSc.DMSTAG_ELEMENT, 0);
        Pe          =   PETSc.DMStagGetGhostArrayLocationSlot(user_ctx.dm,user_ctx.x_l,     PETSc.DMSTAG_ELEMENT, 1);

        p1 = plot(Phi[1:end-1], Z_cen[1:end-1],  ylabel="Z", xlabel="ϕ",  xlims=( 0.0, 2.0), label=:none, title="De=$(user_ctx.De)");
        #p2 = plot(Pe[1:end-1],  Z_cen[1:end-1]    ,                       xlabel="Pe", xlims=(-1.25, 1.25), label=:none, title="$(round(time;sigdigits=3))");
        p2 = plot(Pe[1:end-1],  Z_cen[1:end-1]    ,                       xlabel="Pe", xlims=(-0.025, 0.01), label=:none, title="$(round(time;sigdigits=3))");

       # p1 = plot(Phi[1:end-1], Z_cen[1:end-1],  ylabel="Z", xlabel="ϕ",   label=:none, title="De=$(user_ctx.De)");
       # p2 = plot(Pe[1:end-1],  Z_cen[1:end-1]    ,                       xlabel="Pe",  label=:none, title="$(round(time;sigdigits=3))");

        plot(p1, p2, layout=(1,2)); frame(anim)

        Base.finalize(Pe);  Base.finalize(Phi)
    end

    println("Timestep $it, time=$time")

end

if CreatePlots==true
    gif(anim, "Example_1D.gif", fps = 15)   # create a gif animation
end


text:
using HTTP, JSON
const headers = (("content-type", "application/json"), ("Accept", "application/json"))

struct IBMQUser #store other info here
    id::String
end

struct Qobj
    data::Dict{String, Any}
end

struct IBMQReg{B} <: AbstractRegister{B}
    id::String
    device::String
    nactive::B
    nqubits::B
end

struct Job
    reg::IBMQReg
    data::Qobj
    jobid::String
end

Yao.nactive(reg::IBMQReg) = reg.nactive
Yao.nqubits(reg::IBMQReg) = reg.nqubits

function YaoBlocks.apply!(reg::IBMQReg, qc::Array{ChainBlock{N}}) where N
    qobj = yaotoqobj(qc, reg.device)
    job = run(reg, qobj)
end

"""
Login method
"""
function authenticate(token::String = "") # todo : save the token
    if length(token) == 0
        println("IBM QE token > ")
        token = readline()
    end
    conf = (readtimeout = 1,     # todo handle timeouts better
        pipeline_limit = 4,
        retry = false,
        redirect = false)
    url = "https://auth.quantum-computing.ibm.com/api/users/loginWithToken"
    req = Dict("apiToken" => token)
    print("Logging You in...")
    response = HTTP.post(url, headers, JSON.json(req); conf...)

    response_json = String(response.body)
    response_parsed = JSON.parse(response_json)
    if response.status == 200
        println("✔")
    else
        println("❌")
    end
    IBMQUser(response_parsed["id"])
end

"""
Get backends info
"""
function createreg(user::IBMQUser)
    id = user.id
    url_back = "https://api.quantum-computing.ibm.com/api/Network/ibm-q/Groups/open/Projects/main/devices/v/1?access_token=$(id)"
    print("Fetching Backends...")
    response_back = HTTP.get(url_back)
    if response_back.status == 200
        println("✔")
    else
        println("❌")
        return
    end
    response_back_json = String(response_back.body)
    response_back_parsed = JSON.parse(response_back_json)
    id_back = [i["backend_name"] for i in response_back_parsed]
    confirm = false
    while !confirm
        println("The following backends are available > ")
        for i in enumerate(id_back)
            println(i)
        end
        println("Enter the serial number of the backend you wish to use")
        backend_no = readline()
        device_info = response_back_parsed[parse(Int64, backend_no)]
        n_qubits=device_info["n_qubits"]
        println("n_qubits = $(n_qubits)")
        println("basis_gates=$(device_info["basis_gates"])")
        println("Confirm? (Y/N)")
        conf = readline()
        if conf !== "N" || conf !== "n"
            return IBMQReg(id, id_back[parse(Int, backend_no)], n_qubits, n_qubits)
        else
            confirm = false
        end
    end
end

function run(reg::IBMQReg, qobj::Qobj)
    url = "https://api.quantum-computing.ibm.com/api/Network/ibm-q/Groups/open/Projects/main/Jobs?access_token=$(reg.id)"

    req = Dict("backend" => Dict("name" => reg.device), "allowObjectStorage" => true, "shareLevel"=> "none")
    print("Connecting to $(reg.device)...")
    request = HTTP.post(url , headers, JSON.json(req))
    if request.status == 200
        println("✔")
    else
        println("❌")
    end
    response_json = String(request.body)
    response_parsed = JSON.parse(response_json)

    objectinfo = response_parsed["objectStorageInfo"]
    upload_url = objectinfo["uploadUrl"]
    jobid = response_parsed["id"]

    print("Preparing Data...")
    json = JSON.json(qobj.data)
    println("✔")

    print("Uploading circuit to $(reg.device)...")
    ckt_upload = HTTP.put(upload_url, [], json)
    if ckt_upload.status == 200
        println("✔")
        print("Notifying backend...")
        # Notify the backend that the job has been uploaded
        url = "https://api.quantum-computing.ibm.com/api/Network/ibm-q/Groups/open/Projects/main/Jobs/$(jobid)/jobDataUploaded?access_token=$(reg.id)"
        json_step4 ="""{
                "data": "none",
                "json": "none"
                }"""
        request = HTTP.post(url, headers, json_step4)
        println("✔")
        return Job(reg, qobj, jobid)
    else
        println("❌")
    end
end

function status(job::Job)
    jobid = job.jobid
    id = job.reg.id
    url = "https://api.quantum-computing.ibm.com/api/Network/ibm-q/Groups/open/Projects/main/Jobs/$(jobid)?access_token=$(id)"
    result = HTTP.get(url)
    response_json = String(result.body)
    response_parsed = JSON.parse(response_json)
    return response_parsed["status"]
end

# get result
function getresult(job::Job)
    jobid = job.jobid
    id = job.reg.id
    url = "https://api.quantum-computing.ibm.com/api/Network/ibm-q/Groups/open/Projects/main/Jobs/$(jobid)/resultDownloadUrl?access_token=$(id)"
    print("Connecting to $(job.reg.device)...")
    result = HTTP.get(url)
    if result.status == 200
        println("✔")
    else
        println("❌")
    end
    response_json = String(result.body)
    response_parsed = JSON.parse(response_json)
    downloadUrl = response_parsed["url"]
    final_res = HTTP.get(downloadUrl)

    if final_res.status == 200
        print("Fetching result from $(job.reg.device)...")
        # STEP8: Confirm the data was downloaded
        url = "https://api.quantum-computing.ibm.com/api/Network/ibm-q/Groups/open/Projects/main/Jobs/$(jobid)/resultDownloaded?access_token=$(id)"
        json_step8 = """{
            "data": "none",
            "json": "none"
            }"""
        result = HTTP.post(url, headers, json_step8)
        if result.status == 200
            println("✔")
        else
            println("❌")
        end
        response_json = String(final_res.body)
        response_parsed = JSON.parse(response_json)
        res = response_parsed["results"]
    end
end


text:
using OrdinaryDiffEq
using Trixi
using Printf
using PrettyTables

######### EINSTELLUNGEN ############
N_vec = [3, 4]
Nq_vec = [2, 4, 8, 16, 32]
CFL = 0.8
# timespan for periodic solution
tspan = (0.0, 2.0)
latex = false

surface_flux = FluxPlusDissipation(flux_chandrashekar, DissipationLocalLaxFriedrichs(max_abs_speed_naive)) #flux_lax_friedrichs
volume_flux  = flux_chandrashekar
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
# volume_integral = Trixi.VolumeIntegralPseudoStrongForm()

# mapping as described in the worksheet
function mapping(xi_, eta_)

    xi = xi_
    eta = eta_

    x = xi + 0.15 * cos(0.5 * pi * xi) * cos((3/2) * pi * eta)
    y = eta + 0.15 * cos(2 * pi * xi) * cos(0.5 * pi * eta)

    return SVector(x, y)
  end

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

# copied from EUlerEquation2D to do modifications
function WR2_initial_condition_convergence_test(x, t, equations::CompressibleEulerEquations2D)
    c = 2
    A = 0.1
    L = 2
    f = 1 / L
    ω = 2 * pi * f
    ini = c + A * sin(ω * (x[1] + x[2] - t))

    rho = ini
    rho_v1 = ini
    rho_v2 = ini
    rho_e = ini^2

    return SVector(rho, rho_v1, rho_v2, rho_e)
end
initial_condition = WR2_initial_condition_convergence_test

# coordinates neccessary for periodic solution with periodic boundaries
coordinates_min = (0.0, 0.0)
coordinates_max = (2.0,  2.0)

####

    # Table preparation
kwargs = Dict{Symbol,Any}(
        :formatters => (
            ft_printf("%.5e", [2]),
            ft_printf("%.2f", [3])
        )
    )

if latex
    kwargs[:backend] = :latex
    kwargs[:highlighters] = LatexHighlighter(
            (data, i, j) -> j > 1,
            (data, i, j, str) -> "\$\\num{$str}\$"
        )
end

error = Vector{Float64}(undef, length(Nq_vec))
for N in N_vec

    basis = LobattoLegendreBasis(N)
    solver = DGSEM(basis, surface_flux, volume_integral)
    println(typeof(solver.volume_integral))
    # solver = DGSEM(polydeg=N, surface_flux=surface_flux, volume_integral=volume_integral)
    for j in eachindex(Nq_vec)
        Nq = Nq_vec[j]

        cells_per_dimension = (Nq, Nq)

        mesh = CurvedMesh(cells_per_dimension, mapping)#(-1.0, -1.0), (1.0, 1.0))
        # mesh = CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max)

        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms=source_terms_convergence_test)

            ###############################################################################
            # ODE solvers, callbacks etc.

        ode = semidiscretize(semi, tspan)

        summary_callback = SummaryCallback()
        analysis_interval = 1000

        analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
                                                #  extra_analysis_errors=(:conservation_error,),
                                                #  extra_analysis_integrals=(entropy, energy_total,
                                                #                            energy_kinetic, energy_internal)

                                                )

            # alive_callback = AliveCallback(analysis_interval=analysis_interval)

            # save_solution = SaveSolutionCallback(interval=100,
            #                                      save_initial_solution=true,
            #                                      save_final_solution=true,
            #                                      solution_variables=cons2prim)

        stepsize_callback = StepsizeCallback(cfl= CFL)

        callbacks = CallbackSet(
                                    # summary_callback,
                                    analysis_callback,
                                    # alive_callback,
                                    # save_solution,
                                    stepsize_callback)

            ###############################################################################
            # run the simulation

        sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
                        dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
                        save_everystep=false, callback=callbacks);
        summary_callback() # print the timer summary

        error[j] = maximum(abs.(sol.u[1] - sol.u[2])) # maximum(abs.(exact - sol))
    end

        # Compute EOC
    eoc = vcat("", [log(error[j + 1] / error[j]) /
            log(Nq_vec[j] / Nq_vec[j + 1]) for j in 1:length(Nq_vec) - 1])
    println("N = $N")

    pretty_table(hcat(Nq_vec, error, eoc), ["Nq", "Error", "EOC"]; kwargs...)

    mean = sum(eoc[2:end]) / (length(Nq_vec) - 1)
    @printf("Mean EOC: %.2f\n", mean)
end
####

# built in convergence analysis doesnt work yet
# default_example() = joinpath(examples_dir(), "2d", "elixir_advection_basic.jl")
# convergence_test(default_example(), 4)

#   "c:\\Users\\jmbr2\\git\\TrixiFork\\Trixi.jl\\src\\WR2\\milestone1.jl"
# convergence_test("c:\\Users\\jmbr2\\git\\TrixiFork\\Trixi.jl\\src\\WR2\\milestone1.jl", 4)

max_stars_repo_path: src/WR2/MS3 high order curved.jl

max_issues_repo_path: src/WR2/MS3 high order curved.jl

max_forks_repo_path: src/WR2/MS3 high order curved.jl


text:
input = open("../inputs/day06.txt") do f
    read(f, String)
end

groups = [
    [
        Set([c for c in ln])
        for ln in split(group, "\n")
    ]
    for group in split(input, "\n\n")
]

part1 = sum([length(union(group...)) for group in groups])
println("Part 1: ", part1)

part2 = sum([length(intersect(group...)) for group in groups])
println("Part 2: ", part2)


text:
"""
# module ShootingMethod

# Examples

```jldoctest
julia>
```
"""
module ShootingMethod

using Roots: find_zero, AbstractUnivariateZeroMethod

using NumericalMethodsInQuantumMechanics.EigenvalueProblems:
    InitialCondition, BoundaryCondition
using NumericalMethodsInQuantumMechanics.EigenvalueProblems.NumerovMethod: integrate

export shoot

# last means ``y(t1)``, `ic` contains a guess, `bc[2]` is ``y1``.
# `last(integrate(ic_guess, r, g, s)) - bc[2]` means ``y(t1; yd0) - y1``.
function shoot(g, s, bc::BoundaryCondition, init, args...)
    f = guess -> last(integrate(g, s, InitialCondition(bc.y0, guess), args...)) - bc.y1
    # `x` of the root finding is `y′0`.
    return find_zero(f, init)
end
shoot(g::AbstractVector, bc::BoundaryCondition, init, args...) =
    shoot(g, zeros(length(g)), bc, init, args...)
shoot(g::Function, bc::BoundaryCondition, init, args...) = shoot(g, zero, bc, init, args...)

end


text:
Double64(x::Float64) = Double64((x, 0.0))

for T in (:Float32, :Float16, :Int32, :Int16, :Int8)
    @eval Double64(x::$T) = ((Float64(x), 0.0))
end

for T in (:Int64, :Int128)
  @eval function Double64(x::$T)
    if abs(x) <= maxintfloat(Float64)
        hi = Float64(x)
        lo = 0.0
    else
        bf = Float128(x)
        hi = Float64(x)
        lo = Float64(bf - hi)
    end
    return Double64((hi, lo))
  end
end

function Double64(x::BigFloat)
    hi = Float64(x)
    lo = Float64(x - hi)
    return Double64((hi, lo))
end

function Double64(x::Float128)
    hi = Float64(x)
    lo = Float64(x - hi)
    return Double64((hi, lo))
end

function Double64(x::BigInt)
    if abs(x) <= typemax(Int128)
        Double64(Int128(x))
    else
         Double64(BigFloat(x))
    end
end

Double64(x::Real, y::Real) = Double64(Float64(x), Float64(y))

Double64(x::ComplexD64) = real(x)
Double64(x::ComplexF64) = Double64(real(x))

ComplexD64(x::Float64) = ComplexD64(Double64(x))
ComplexD64(x::Float64, y::Float64) = ComplexD64( Double64(x), Double64(y) )
ComplexD64(x::ComplexF64) = ComplexD64(x.re, x.im)

# inverse constructors
Base.Float64(x::Double64) = Hi(x)
Base.Float64(x::ComplexD64) = real(Hi(x))
Base.ComplexF64(x::Double64) = ComplexF64(Hi(x), 0.0)
Base.ComplexF64(x::ComplexD64) = ComplexF64(Hi(x)...)
Base.BigFloat(x::Double64) = BigFloat(Hi(x)) + BigFloat(Lo(x))
Base.BigFloat(x::ComplexD64) = BigFloat(real(Hi(x))) + BigFloat(real(Lo(x)))
Base.BigInt(x::Double64) = BigInt(BigFloat(x))
Base.BigInt(x::ComplexD64) = BigInt(BigFloat(x))
Base.Int128(x::Double64) = Int128(BigInt(x))
Base.Int128(x::ComplexD64) = Int128(BigInt(x))

function Base.Int64(x::Double64)
    Int64(Hi(x)) + Int64(Lo(x))
end

for T in (:Float32, :Float16, :Int32, :Int16, :Int8)
    @eval Base.$T(x::Double64) = $T(Hi(x))
end

Float128(x::Double64) = Float128(Hi(x)) + Float128(Lo(x))
Double64(x::Float128) = Double64((Float64(x), Float64(x-Float64(x))))
ComplexD64(x::ComplexF128) = (Double64(real(x)), Double64(imag(x)))
ComplexF128(x::ComplexD64) = ComplexF128(real(x), imag(x))


text:
import Base: *,        Ac_mul_B,  A_mul_Bc,  Ac_mul_Bc,  At_mul_B,  A_mul_Bt,  At_mul_Bt
import Base: A_mul_B!, Ac_mul_B!, A_mul_Bc!, Ac_mul_Bc!, At_mul_B!, A_mul_Bt!, At_mul_Bt!

typealias BlasEltypes Union{Float64, Float32, Complex{Float64}, Complex{Float32}}

# TODO size-inferrable products with AbstractArray (such as StaticMatrix * AbstractVector)
# TODO Potentially a loop version for rather large arrays? Or try and figure out inference problems?

# TODO make faster versions of A*_mul_B*
@generated function A_mul_Bc(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        A * ctranspose(B)
    end
end
@generated function Ac_mul_B(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        ctranspose(A) * B
    end
end
@generated function Ac_mul_Bc(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        ctranspose(B*A) # is this always safe?
    end
end

@generated function A_mul_Bt(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        A * transpose(B)
    end
end
@generated function At_mul_B(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        transpose(A) * B
    end
end
@generated function At_mul_Bt(A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        transpose(B*A) # is this always safe?
    end
end

# mutating

@generated function A_mul_Bc!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        A_mul_B!(C, A, ctranspose(B))
    end
end
@generated function Ac_mul_B!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        A_mul_B!(C, ctranspose(A), B)
    end
end
@generated function Ac_mul_Bc!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        A_mul_B!(C, ctranspose(A), ctranspose(B))
    end
end

@generated function A_mul_Bt!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        A_mul_B!(C, A, transpose(B))
    end
end
@generated function At_mul_B!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        A_mul_B!(C, transpose(A), B)
    end
end
@generated function At_mul_Bt!(C::Union{StaticMatrix, StaticVector}, A::Union{StaticMatrix, StaticVector}, B::Union{StaticMatrix, StaticVector})
    return quote
        $(Expr(:meta, :inline))
        A_mul_B!(C, transpose(A), transpose(B))
    end
end

@generated function *(A::StaticMatrix, b::StaticVector)
    TA = eltype(A)
    Tb = eltype(b)
    sA = size(A)
    sb = size(b)

    s = (sA[1],)
    T = typeof(zero(TA)*zero(Tb))

    if sb[1] != sA[2]
        error("Dimension mismatch")
    end

    if s == sb
        if T == Tb
            newtype = b
        else
            newtype = similar_type(b, T)
        end
    else
        if T == Tb
            newtype = similar_type(b, s)
        else
            newtype = similar_type(b, T, s)
        end
    end

    if sA[2] != 0
        exprs = [reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k, j))]*b[$j]) for j = 1:sA[2]]) for k = 1:sA[1]]
    else
        exprs = [zero(T) for k = 1:sA[1]]
    end

    return quote
        $(Expr(:meta,:inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function *(a::StaticVector, B::StaticMatrix)
    Ta = eltype(a)
    TB = eltype(B)
    sa = size(a)
    sB = size(B)

    s = (sa[1],sB[2])
    T = typeof(zero(Ta)*zero(TB))

    if sB[1] != 1
        error("Dimension mismatch")
    end

    if s == sB
        if T == TB
            newtype = B
        else
            newtype = similar_type(B, T)
        end
    else
        if T == TB
            newtype = similar_type(B, s)
        else
            newtype = similar_type(B, T, s)
        end
    end

    exprs = [:(a[$j] * B[$k]) for j = 1:s[1], k = 1:s[2]]

    return quote
        $(Expr(:meta,:inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

#@inline *{S1,S2,S3}(A::MMatrix{S1,S2}, B::MMatrix{S2,S3}) = MMatrix{S1,S3}(SMatrix{S1,S2}(A)*SMatrix{S2,S3}(B))

@generated function *(A::StaticMatrix, B::StaticMatrix)
    TA = eltype(A)
    TB = eltype(B)

    T = typeof(zero(TA)*zero(TB))

    can_mutate = !isbits(A) || !isbits(B) # !isbits implies can get a persistent pointer (to pass to BLAS). Probably will change to !isimmutable in a future version of Julia.
    can_blas = T == TA && T == TB && T <: Union{Float64, Float32, Complex{Float64}, Complex{Float32}}

    if can_mutate
        sA = size(A)
        sB = size(B)
        s = (sA[1], sB[2])

        # Heuristic choice between BLAS and explicit unrolling (or chunk-based unrolling)
        if can_blas && size(A,1)*size(A,2)*size(B,2) >= 14*14*14
            return quote
                $(Expr(:meta, :inline))
                C = similar(A, $T, $s)
                A_mul_B_blas!(C, A, B)
                return C
            end
        elseif size(A,1)*size(A,2)*size(B,2) < 8*8*8
            return quote
                $(Expr(:meta, :inline))
                return A_mul_B_unrolled(A, B)
            end
        elseif size(A,1) <= 14 && size(A,2) <= 14 && size(B,2) <= 14
            return quote
                $(Expr(:meta, :inline))
                return $(similar_type(A, T, s))(A_mul_B_unrolled_chunks(A, B))
            end
        else
            return quote
                $(Expr(:meta, :inline))
                return A_mul_B_loop(A, B)
            end
        end
    else # both are isbits type...
        # Heuristic choice for amount of codegen
        if size(A,1)*size(A,2)*size(B,2) <= 8*8*8
            return quote
                $(Expr(:meta, :inline))
                return A_mul_B_unrolled(A, B)
            end
        elseif size(A,1) <= 14 && size(A,2) <= 14 && size(B,2) <= 14
            return quote
                $(Expr(:meta, :inline))
                return A_mul_B_unrolled_chunks(A, B)
            end
        else
            return quote
                $(Expr(:meta, :inline))
                return A_mul_B_loop(A, B)
            end
        end
    end
end

@generated function A_mul_B_unrolled(A::StaticMatrix, B::StaticMatrix)
    sA = size(A)
    sB = size(B)
    TA = eltype(A)
    TB = eltype(B)

    s = (sA[1], sB[2])
    T = typeof(zero(TA)*zero(TB))

    if sB[1] != sA[2]
        error("Dimension mismatch")
    end

    # TODO think about which to be similar to
    if s == sB
        if T == TB
            newtype = B
        else
            newtype = similar_type(B, T)
        end
    else
        if T == TB
            newtype = similar_type(B, s)
        else
            newtype = similar_type(B, T, s)
        end
    end

    if sA[2] != 0
        exprs = [reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k1, j))]*B[$(sub2ind(sB, j, k2))]) for j = 1:sA[2]]) for k1 = 1:sA[1], k2 = 1:sB[2]]
    else
        exprs = [zero(T) for k1 = 1:sA[1], k2 = 1:sB[2]]
    end

    return quote
        $(Expr(:meta,:inline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

@generated function A_mul_B_loop(A::StaticMatrix, B::StaticMatrix)
    sA = size(A)
    sB = size(B)
    TA = eltype(A)
    TB = eltype(B)

    s = (sA[1], sB[2])
    T = typeof(zero(TA)*zero(TB))

    if sB[1] != sA[2]
        error("Dimension mismatch")
    end

    # TODO think about which to be similar to
    if s == sB
        if T == TB
            newtype = B
        else
            newtype = similar_type(B, T)
        end
    else
        if T == TB
            newtype = similar_type(B, s)
        else
            newtype = similar_type(B, T, s)
        end
    end

    tmps = [Symbol("tmp_$(k1)_($k2)") for k1 = 1:sA[1], k2 = 1:sB[2]]
    exprs_init = [:($(tmps[k1,k2])  = A[$k1] * B[1 + $((k2-1) * sB[1])]) for k1 = 1:sA[1], k2 = 1:sB[2]]
    exprs_loop = [:($(tmps[k1,k2]) += A[$(k1-sA[1]) + $(sA[1])*j] * B[j + $((k2-1) * sB[1])]) for k1 = 1:sA[1], k2 = 1:sB[2]]

    return quote
        $(Expr(:meta,:inline))

        @inbounds $(Expr(:block, exprs_init...))
        for j = 2:$(sA[2])
            @inbounds $(Expr(:block, exprs_loop...))
        end
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, tmps...)))
    end
end

# Concatenate a series of matrix-vector multiplications
# Each function is N^2 not N^3 - aids in compile time.
@generated function A_mul_B_unrolled_chunks(A::StaticMatrix, B::StaticMatrix)
    sA = size(A)
    sB = size(B)
    TA = eltype(A)
    TB = eltype(B)

    s = (sA[1], sB[2])
    T = typeof(zero(TA)*zero(TB))

    if sB[1] != sA[2]
        error("Dimension mismatch")
    end

    # TODO think about which to be similar to
    if s == sB
        if T == TB
            newtype = B
        else
            newtype = similar_type(B, T)
        end
    else
        if T == TB
            newtype = similar_type(B, s)
        else
            newtype = similar_type(B, T, s)
        end
    end

    #vect_exprs = [:($(Symbol("tmp_$k2")) = partly_unrolled_multiply(A, B[:, $k2])) for k2 = 1:sB[2]]

    # Do a custom B[:, k2] to return a SVector (an isbits type) rather than (possibly) a mutable type. Avoids allocation == faster
    tmp_type_in = SVector{sB[1], T}
    tmp_type_out = SVector{sA[1], T}
    vect_exprs = [:($(Symbol("tmp_$k2"))::$tmp_type_out = partly_unrolled_multiply(A, $(Expr(:call, tmp_type_in, [Expr(:ref, :B, sub2ind(s, i, k2)) for i = 1:sB[1]]...)))::$tmp_type_out) for k2 = 1:sB[2]]

    exprs = [:($(Symbol("tmp_$k2"))[$k1]) for k1 = 1:sA[1], k2 = 1:sB[2]]

    return Expr(:block,
        Expr(:meta,:inline),
        vect_exprs...,
        :(@inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...))))
    )
end

@generated function partly_unrolled_multiply(A::StaticMatrix, b::StaticVector)
    TA = eltype(A)
    Tb = eltype(b)
    sA = size(A)
    sb = size(b)

    s = (sA[1],)
    T = typeof(zero(TA)*zero(Tb))

    if sb[1] != sA[2]
        error("Dimension mismatch")
    end

    if s == sb
        if T == Tb
            newtype = b
        else
            newtype = similar_type(b, T)
        end
    else
        if T == Tb
            newtype = similar_type(b, s)
        else
            newtype = similar_type(b, T, s)
        end
    end

    if sA[2] != 0
        exprs = [reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k, j))]*b[$j]) for j = 1:sA[2]]) for k = 1:sA[1]]
    else
        exprs = [zero(T) for k = 1:sA[1]]
    end

    return quote
        $(Expr(:meta,:noinline))
        @inbounds return $(Expr(:call, newtype, Expr(:tuple, exprs...)))
    end
end

#function A_mul_B_blas(a, b, c, A, B)
#q
#end

# The idea here is to get pointers to stack variables and call BLAS.
# This saves an aweful lot of time compared to copying SArray's to Ref{SArray{...}}
# and should be fastest for (very) large SArrays

# Here is an LLVM function that gets the pointer to its input, %x
# After this we would make the ccall above.
#
# define i8* @f(i32 %x) #0 {
#     %1 = alloca i32, align 4
#     store i32 %x, i32* %1, align 4
#     ret i32* %1
# }

@generated function A_mul_B!(C::StaticMatrix, A::StaticMatrix, B::StaticMatrix)
    if isbits(C)
        error("Cannot mutate $C")
    end

    TA = eltype(A)
    TB = eltype(B)
    T = typeof(zero(TA)*zero(TB))

    can_blas = T == TA && T == TB && T <: Union{Float64, Float32, Complex{Float64}, Complex{Float32}}

    if can_blas
        if size(A,1) * size(A,2) * size(A,3) < 4*4*4
            return quote
                $(Expr(:meta, :inline))
                A_mul_B_unrolled!(C, A, B)
                return C
            end
        elseif size(A,1) * size(A,2) * size(A,3) < 14*14*14 # Something seems broken for this one with large matrices (becomes allocating)
            return quote
                $(Expr(:meta, :inline))
                A_mul_B_unrolled_chunks!(C, A, B)
                return C
            end
        else
            return quote
                $(Expr(:meta, :inline))
                A_mul_B_blas!(C, A, B)
                return C
            end
        end
    else
        if size(A,1) * size(A,2) * size(A,3) < 4*4*4
            return quote
                $(Expr(:meta, :inline))
                A_mul_B_unrolled!(C, A, B)
                return C
            end
        else
            return quote
                $(Expr(:meta, :inline))
                A_mul_B_unrolled_chunks!(C, A, B)
                return C
            end
        end
    end
end

@generated function A_mul_B_blas!(C::StaticMatrix, A::StaticMatrix, B::StaticMatrix)
    sA = size(A)
    sB = size(B)
    TA = eltype(A)
    TB = eltype(B)

    s = size(C)
    T = eltype(C)

    if sB[1] != sA[2] || sA[1] != s[1] || sB[2] != s[2]
        error("Dimension mismatch")
    end

    if sA[1] > 0 && sA[2] > 0 && sB[2] > 0 && T == TA && T == TB && T <: Union{Float64, Float32, Complex{Float64}, Complex{Float32}}
        # This code adapted from `gemm!()` in base/linalg/blas.jl

        if T == Float64
            gemm = :dgemm_
        elseif T == Float32
            gemm = :sgemm_
        elseif T == Complex{Float64}
            gemm = :zgemm_
        else # T == Complex{Float32}
            gemm = :cgemm_
        end

        return quote
            alpha = $(one(T))
            beta = $(zero(T))
            transA = 'N'
            transB = 'N'
            m = $(sA[1])
            ka = $(sA[2])
            kb = $(sB[1])
            n = $(sB[2])
            strideA = $(sA[1])
            strideB = $(sB[1])
            strideC = $(s[1])

            ccall((Base.BLAS.@blasfunc($gemm), Base.BLAS.libblas), Void,
                (Ptr{UInt8}, Ptr{UInt8}, Ptr{Base.BLAS.BlasInt}, Ptr{Base.BLAS.BlasInt},
                 Ptr{Base.BLAS.BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{Base.BLAS.BlasInt},
                 Ptr{$T}, Ptr{Base.BLAS.BlasInt}, Ptr{$T}, Ptr{$T},
                 Ptr{Base.BLAS.BlasInt}),
                 &transA, &transB, &m, &n,
                 &ka, &alpha, A, &strideA,
                 B, &strideB, &beta, C,
                 &strideC)
            return C
        end
    else
        error("Cannot call BLAS gemm with $C = $A * $B")
    end
end

@generated function A_mul_B_unrolled!(C::StaticMatrix, A::StaticMatrix, B::StaticMatrix)
    sA = size(A)
    sB = size(B)
    TA = eltype(A)
    TB = eltype(B)
    T = eltype(C)

    s = (sA[1], sB[2])

    if sB[1] != sA[2]
        error("Dimension mismatch")
    end

    if s != size(C)
        error("Dimension mismatch")
    end

    if sA[2] != 0
        exprs = [:(C[$(sub2ind(s, k1, k2))] = $(reduce((ex1,ex2) -> :(+($ex1,$ex2)), [:(A[$(sub2ind(sA, k1, j))]*B[$(sub2ind(sB, j, k2))]) for j = 1:sA[2]]))) for k1 = 1:sA[1], k2 = 1:sB[2]]
    else
        exprs = [:(C[$(sub2ind(s, k1, k2))] = $(zero(T))) for k1 = 1:sA[1], k2 = 1:sB[2]]
    end

    return quote
        $(Expr(:meta,:inline))
        @inbounds $(Expr(:block, exprs...))
    end
end

@generated function A_mul_B_unrolled_chunks!(C::StaticMatrix, A::StaticMatrix, B::StaticMatrix)
    sA = size(A)
    sB = size(B)
    TA = eltype(A)
    TB = eltype(B)

    s = size(C)
    T = eltype(C)

    if sB[1] != sA[2] || sA[1] != s[1] || sB[2] != s[2]
        error("Dimension mismatch")
    end

    #vect_exprs = [:($(Symbol("tmp_$k2")) = partly_unrolled_multiply(A, B[:, $k2])) for k2 = 1:sB[2]]

    # Do a custom B[:, k2] to return a SVector (an isbits type) rather than a mutable type. Avoids allocation == faster
    tmp_type = SVector{sB[1], T}
    vect_exprs = [:($(Symbol("tmp_$k2")) = partly_unrolled_multiply(A, $(Expr(:call, tmp_type, [Expr(:ref, :B, sub2ind(s, i, k2)) for i = 1:sB[1]]...)))) for k2 = 1:sB[2]]

    exprs = [:(C[$(sub2ind(s, k1, k2))] = $(Symbol("tmp_$k2"))[$k1]) for k1 = 1:sA[1], k2 = 1:sB[2]]

    return Expr(:block,
        Expr(:meta,:inline),
        vect_exprs...,
        exprs...
    )
end


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

"""
    Polytope{K,Dim,T}

We say that a geometry is a K-polytope when it is a collection of "flat" sides
that constitue a `K`-dimensional subspace. They are called polygon and polyhedron
respectively for 2D (`K=2`) and 3D (`K=3`) subspaces, embedded in a `Dim`-dimensional
space. The parameter `K` is also known as the rank or parametric dimension of the
polytope: https://en.wikipedia.org/wiki/Abstract_polytope.

The term polytope expresses a particular combinatorial structure. A polyhedron,
for example, can be decomposed into faces. Each face can then be decomposed into
edges, and edges into vertices. Some conventions act as a mapping between vertices
and higher dimensional features (edges, faces, cells...), removing the need to
store all features.

Additionally, the following property must hold in order for a geometry to be considered
a polytope: the boundary of a (K+1)-polytope is a collection of K-polytopes, which may
have (K-1)-polytopes in common. See https://en.wikipedia.org/wiki/Polytope.

### Notes

- Type aliases are `Polygon`, `Polyhedron`.
"""
abstract type Polytope{K,Dim,T} <: Geometry{Dim,T} end

(::Type{PL})(vertices::Vararg{P}) where {PL<:Polytope,P<:Point} = PL(SVector(vertices))
(::Type{PL})(vertices::AbstractVector{TP}) where {PL<:Polytope,TP<:Tuple} = PL(Point.(vertices))
(::Type{PL})(vertices::Vararg{TP}) where {PL<:Polytope,TP<:Tuple} = PL(collect(vertices))

"""
    paramdim(polytope)

Return the parametric dimension or rank of the polytope.
"""
paramdim(::Type{<:Polytope{K}}) where {K} = K

"""
    vertices(polytope)

Return the vertices of a `polytope`.
"""
vertices(p::Polytope) = p.vertices

"""
    nvertices(polytope)

Return the number of vertices in the `polytope`.
"""
nvertices(p::Polytope) = length(vertices(p))

"""
    p1 == p2

Tells whether or not polytopes `p1` and `p2` are equal.
"""
==(p1::Polytope, p2::Polytope) = vertices(p1) == vertices(p2)

"""
    centroid(polytope)

Return the centroid of the `polytope`.
"""
centroid(p::Polytope) = Point(sum(coordinates, vertices(p)) / length(vertices(p)))

"""
    unique(polytope)

Return a new `polytope` without duplicate vertices.
"""
Base.unique(p::Polytope) = unique!(deepcopy(p))

"""
    unique!(polytope)

Remove duplicate vertices in `polytope`.
"""
function Base.unique!(::Polytope) end

function Base.show(io::IO, p::Polytope)
  kind = prettyname(p)
  vert = join(vertices(p), ", ")
  print(io, "$kind($vert)")
end

function Base.show(io::IO, ::MIME"text/plain", p::Polytope{K,Dim,T}) where {K,Dim,T}
  kind = prettyname(p)
  lines = ["  └─$v" for v in vertices(p)]
  println(io, "$kind{$Dim,$T}")
  print(io, join(lines, "\n"))
end

prettyname(p::Polytope) = prettyname(typeof(p))
function prettyname(PL::Type{<:Polytope})
  n = string(PL)
  i = findfirst('{', n)
  isnothing(i) ? n : n[1:i-1]
end

# --------
# POLYGON
# --------

"""
    Polygon{Dim,T}

A polygon is a 2-polytope, i.e. a polytope with parametric dimension 2.
"""
const Polygon = Polytope{2}

"""
    area(polygon)

Return the area of the `polygon`.
"""
area(p::Polygon) = measure(p)

"""
    chains(polygon)

Return the outer and inner chains of the polygon.
"""
function chains(::Polygon) end

"""
    hasholes(polygon)

Tells whether or not the `polygon` contains holes.
"""
function hasholes(::Polygon) end

"""
    issimple(polygon)

Tells whether or not the `polygon` is simple.
See [https://en.wikipedia.org/wiki/Simple_polygon]
(https://en.wikipedia.org/wiki/Simple_polygon).
"""
issimple(p::Polygon) = issimple(typeof(p))

"""
    windingnumber(point, polygon)

Winding number of `point` with respect to the `polygon`.
"""
function windingnumber(::Point, ::Polygon) end

"""
    orientation(polygon)

Returns the orientation of the rings of the `polygon`
as either counter-clockwise (CCW) or clockwise (CW).
"""
orientation(p::Polygon) = orientation(p, WindingOrientation())

function orientation(p::Polygon, algo)
  o = [orientation(c, algo) for c in chains(p)]
  hasholes(p) ? o : first(o)
end

"""
    boundary(polygon)

Returns the boundary of the `polygon`.
"""
boundary(p::Polygon) = hasholes(p) ? Multi(chains(p)) : first(chains(p))

"""
    isconvex(polygon)

Tells whether or not the `polygon` is convex.
"""
isconvex(p::Polygon{Dim,T}) where {Dim,T} =
  issimple(p) && all(≤(T(π)), innerangles(boundary(p)))

"""
    bridge(polygon; width=0)

Transform `polygon` with holes into a single outer chain
via bridges of given `width` as described in Held 1998.
Return the outer chain and a vector with pairs of indices
for duplicate vertices. These indices can be used to undo
the bridges.

## References

* Held. 1998. [FIST: Fast Industrial-Strength Triangulation of Polygons]
  (https://link.springer.com/article/10.1007/s00453-001-0028-4)
"""
function bridge(p::Polygon{Dim,T}; width=zero(T)) where {Dim,T}
  # polygons without holes are trivial
  if !hasholes(p)
    outerchain = first(chains(p))
    duplicates = Tuple{Int,Int}[]
    return outerchain, duplicates
  end

  # retrieve chains as vectors of coordinates
  pchains = [coordinates.(vertices(open(c))) for c in chains(p)]

  # sort vertices lexicographically
  coords  = [coord for pchain in pchains for coord in pchain]
  indices = sortperm(sortperm(coords))

  # each chain has its own set of indices
  pinds = Vector{Int}[]; offset = 0
  for nvertex in length.(pchains)
    push!(pinds, indices[offset+1:offset+nvertex])
    offset += nvertex
  end

  # sort chains based on leftmost vertex
  leftmost = argmin.(pinds)
  minimums = getindex.(pinds, leftmost)
  reorder  = sortperm(minimums)
  leftmost = leftmost[reorder]
  minimums = minimums[reorder]
  pchains  = pchains[reorder]
  pinds    = pinds[reorder]

  # initialize outer boundary
  outer = first(pchains)
  oinds = first(pinds)

  # merge holes into outer boundary
  for i in 2:length(pchains)
    inner = pchains[i]
    iinds = pinds[i]
    l = leftmost[i]
    m = minimums[i]

    # find closest vertex in boundary
    dmin, jmin = Inf, 0
    for j in findall(oinds .≤ m)
      d = sum(abs, outer[j] - inner[l])
      if d < dmin
        dmin, jmin = d, j
      end
    end

    # create a bridge of given width δ
    # from line segment A--B. The point
    # A is split into A′ and A′′ and the
    # point B is split into B′ and B′′
    A = outer[jmin]
    B = inner[l]
    δ = width
    v = B - A
    u = Vec(-v[2], v[1])
    n = u / norm(u)
    A′  = A + δ/2 * n
    A′′ = A - δ/2 * n
    B′  = B + δ/2 * n
    B′′ = B - δ/2 * n

    # insert hole at closest vertex
    outer = [
      outer[begin:jmin-1]; [A′, B′];
      circshift(inner, -l+1)[2:end];
      [B′′, A′′]; outer[jmin+1:end]
    ]
    oinds = [
      oinds[begin:jmin];
      circshift(iinds, -l+1);
      [iinds[l]];
      oinds[jmin:end]
    ]
  end

  # find duplicate vertices
  duplicates = Tuple{Int,Int}[]
  occurred = Dict{Int,Int}()
  for (i, ind) in enumerate(oinds)
    if haskey(occurred, ind)
      push!(duplicates, (occurred[ind], i))
    else
      occurred[ind] = i
    end
  end

  # close outer boundary
  push!(outer, first(outer))

  outerchain = Chain(Point.(outer))

  outerchain, duplicates
end

# -----------
# POLYHEDRON
# -----------

"""
    Polyhedron{Dim,T}

A polyhedron is a 3-polytope, i.e. a polytope with parametric dimension 3.
"""
const Polyhedron = Polytope{3}

"""
   volume(polyhedron)

Return the volume of the `polyhedron`.
"""
volume(p::Polyhedron) = measure(p)

# ----------------
# IMPLEMENTATIONS
# ----------------

include("polytopes/segment.jl")
include("polytopes/ngon.jl")
include("polytopes/chain.jl")
include("polytopes/polyarea.jl")
include("polytopes/pyramid.jl")
include("polytopes/tetrahedron.jl")
include("polytopes/hexahedron.jl")


text:
validtestvalues = [
    # Structure: MTrk, length in bytes (4 bytes long), data
    # Perhaps longer than it needs to be, but covers meta events, MIDI events, and running statuses.
    (
        [0x4d, 0x54, 0x72, 0x6b, 0x00, 0x00, 0x00, 0x52, 0x60, 0x90, 0x3c, 0x7f, 0x30, 0x43, 0x7f, 0x30, 0x80, 0x3c, 0x7f, 0x00, 0x90, 0x40, 0x7f, 0x30, 0x80, 0x43, 0x7f, 0x30, 0xc0, 0x3f, 0x00, 0x80, 0x40, 0x7f, 0x00, 0x90, 0x3c, 0x7f, 0x30, 0x43, 0x7f, 0x30, 0x80, 0x3c, 0x7f, 0x00, 0x90, 0x40, 0x7f, 0x30, 0x80, 0x43, 0x7f, 0x30, 0xc0, 0x22, 0x00, 0x80, 0x40, 0x7f, 0x00, 0x90, 0x3c, 0x7f, 0x30, 0x43, 0x7f, 0x30, 0x80, 0x3c, 0x7f, 0x00, 0x90, 0x40, 0x7f, 0x30, 0x80, 0x43, 0x7f, 0x30, 0xc0, 0x52, 0x00, 0x80, 0x40, 0x7f, 0x00, 0xff, 0x2f, 0x00],
        MIDI.MIDITrack(MIDI.TrackEvent[MIDI.MIDIEvent(96,0x90,UInt8[60,127]),MIDI.MIDIEvent(48,0x90,UInt8[67,127]),MIDI.MIDIEvent(48,0x80,UInt8[60,127]),MIDI.MIDIEvent(0,0x90,UInt8[64,127]),MIDI.MIDIEvent(48,0x80,UInt8[67,127]),MIDI.MIDIEvent(48,0xc0,UInt8[63]),MIDI.MIDIEvent(0,0x80,UInt8[64,127]),MIDI.MIDIEvent(0,0x90,UInt8[60,127]),MIDI.MIDIEvent(48,0x90,UInt8[67,127]),MIDI.MIDIEvent(48,0x80,UInt8[60,127]),MIDI.MIDIEvent(0,0x90,UInt8[64,127]),MIDI.MIDIEvent(48,0x80,UInt8[67,127]),MIDI.MIDIEvent(48,0xc0,UInt8[34]),MIDI.MIDIEvent(0,0x80,UInt8[64,127]),MIDI.MIDIEvent(0,0x90,UInt8[60,127]),MIDI.MIDIEvent(48,0x90,UInt8[67,127]),MIDI.MIDIEvent(48,0x80,UInt8[60,127]),MIDI.MIDIEvent(0,0x90,UInt8[64,127]),MIDI.MIDIEvent(48,0x80,UInt8[67,127]),MIDI.MIDIEvent(48,0xc0,UInt8[82]),MIDI.MIDIEvent(0,0x80,UInt8[64,127])])
    ),
]

invalidtestvalues = [
    ([0x00], EOFError),
    ([0x4d, 0x54, 0x72, 0x6b, 0x00, 0x00, 0x00, 0xFF, 0x00], EOFError),
    ([0x4d, 0x54, 0x72, 0x6b], EOFError),
    ([0x4d, 0x54, 0x72, 0x6c], ErrorException),
]

for (input, output) in validtestvalues
    result = MIDI.readtrack(IOBuffer(input))
    @test length(result.events) == length(output.events)
    for (e1, e2) in zip(result.events, output.events)
        e1 == e2
    end
end

for (output, input) in validtestvalues
    buf = IOBuffer()
    MIDI.writetrack(buf, input)
    @test takebuf_array(buf) == output
end

for (input, errtype) in invalidtestvalues
    @test_throws errtype MIDI.readtrack(IOBuffer(input))
end

# Test writing notes and program change events to a track
C = MIDI.Note(60, 96, 0, 0)
G = MIDI.Note(67, 96, 48, 0)
E = MIDI.Note(64, 96, 96, 0)
inc = 96

track = MIDI.MIDITrack()
notes = MIDI.Note[]
for v in UInt8[1,2,3]
    push!(notes, C)
    push!(notes, E)
    push!(notes, G)
    MIDI.programchange(track, E.position + inc + inc, UInt8(0), v)
    C.position += inc
    E.position += inc
    G.position += inc
    C = MIDI.Note(60, 96, C.position+inc, 0)
    E = MIDI.Note(64, 96, E.position+inc, 0)
    G = MIDI.Note(67, 96, G.position+inc, 0)
end

MIDI.addnotes(track, notes)

buf = IOBuffer()
MIDI.writetrack(buf, track)
@test takebuf_array(buf) == [0x4d, 0x54, 0x72, 0x6b, 0x00, 0x00, 0x00, 0x52, 0x60, 0x90, 0x3c, 0x7f, 0x30, 0x43, 0x7f, 0x30, 0x80, 0x3c, 0x7f, 0x00, 0x90, 0x40, 0x7f, 0x30, 0x80, 0x43, 0x7f, 0x30, 0xc0, 0x00, 0x00, 0x80, 0x40, 0x7f, 0x00, 0x90, 0x3c, 0x7f, 0x30, 0x43, 0x7f, 0x30, 0x80, 0x3c, 0x7f, 0x00, 0x90, 0x40, 0x7f, 0x30, 0x80, 0x43, 0x7f, 0x30, 0xc0, 0x01, 0x00, 0x80, 0x40, 0x7f, 0x00, 0x90, 0x3c, 0x7f, 0x30, 0x43, 0x7f, 0x30, 0x80, 0x3c, 0x7f, 0x00, 0x90, 0x40, 0x7f, 0x30, 0x80, 0x43, 0x7f, 0x30, 0xc0, 0x02, 0x00, 0x80, 0x40, 0x7f, 0x00, 0xff, 0x2f, 0x00]

sort!(notes, lt=((x, y)->x.position<y.position))

# Test getting notes from a track
for (n1, n2) in zip(notes, MIDI.getnotes(track))
    @test n1 == n2
end


text:
function kendallsTau{T<:AbstractFloat}(Y::Matrix{T})

  n, p = size(Y)
  covM = zeros(T, (p, p))
  denom = 0

  for a=1:p, b=a:p
    if a == b
      covM[a,b] = 1.
      continue
    end
    t = zero(T)
    for i=1:n, j=i+1:n
      @inbounds t += sign((Y[i, a] - Y[j, a])*(Y[i, b] - Y[j, b]))
    end
    t = 2. * t / ((n-1)*n)
    t = sin(π/2*T)
    covM[a,b] = t
    covM[b,a] = t
  end
  covM
end

function nonparanormalCorrelation{T<:AbstractFloat}(Y::Matrix{T})
  n, p = size(Y)
  stdNormal = Normal()
  δn = 1. / (4 * n^(1./4.) * sqrt(π*log(n)))
  IX = zeros(Y)
  for j=1:p
    sortperm!(view(IX, :, j), view(Y, :, j))
  end
  scale!(IX, 1./n)
  for j=1:p, i=1:n
    if IX[i,j] < δn
      IX[i,j] < δn
    elseif IX[i,j] >= 1. - δn
      IX[i,j] = 1. - δn
    end
    IX[i,j] = quantile(stdNormal, IX[i,j])
  end
  return cor(IX)
end


text:
#####################################################
# Set-based coverage.
#
# This is an example of another data structure to
# support generation of test cases. It's a set of sets
# of the remaining covers.
#####################################################

mutable struct SetCoverage{T <: Integer}
    cover::Dict{Tuple{Vararg{T}},Set{Tuple{Vararg{T}}}}
    arity::Array{T, 1}
end

function SetCoverage(arity::Array{T}, n_way::Integer) where {T <: Integer}
    param_cnt = length(arity)
    param_combo = [Tuple(x) for x in combinations(1:param_cnt, n_way)]
    # The dictionary with keys = parameter indices and dictionary
    # value is set of tuples of parameter values.
    cover = Dict(x => Set{NTuple{n_way,T}}() for x in param_combo)
    SetCoverage{T}(cover, arity)
end

eltype(sc::SetCoverage{T}) where T = T
remaining(sc::SetCoverage) = sum([length(aset) for aset in values(sc.cover)])

function build_all_combinations!(sc::SetCoverage, n_way)
    allc = all_combinations(sc.arity, n_way)
    param_cnt = size(allc, 1)
    for cover_idx in 1:size(allc, 2)
        acover = allc[:, cover_idx]
        nonzero = acover .!= 0
        indices = tuple((1:param_cnt)[nonzero]...)
        values = tuple(acover[nonzero]...)
        if indices ∉ keys(sc.cover)
            sc.cover[indices] = Set{NTuple{n_way, Int}}()
        end
        push!(sc.cover[indices], values)
    end
end

function add_coverage!(sc::SetCoverage, entry::Array{T, 1}) where T
    add_cnt = 0  # Record the number removed for visibility.
    for indices in keys(sc.cover)
        value = tuple([entry[x] for x in indices]...)
        if value ∈ sc.cover[indices]
            pop!(sc.cover[indices], value)
            add_cnt += 1
        end
    end
    add_cnt
end

function add_coverage!(sc::SetCoverage, entries::Array{T, 2}) where T
    add_cnt = 0  # Record the number removed for visibility.
    for entry_idx in 1:size(entries, 2)
        add_cnt += add_coverage!(sc, entries[:, entry_idx])
    end
    add_cnt
end

function test_coverage(test_cases, arity, n_way)
    sc = SetCoverage(arity, n_way)
    build_all_combinations!(sc, n_way)
    start = remaining(sc)
    add_coverage!(sc, test_cases)
    finish = remaining(sc)
    (start = start, finish = finish)
end

"""
    tuples_in_trials(trials, n_way)

Find the number of distinct tuples in a set of trials. This uses a different method
to count and store those tuples, as a check on other methods. It returns a dictionary
where each key is the index of non-zero parameters and the value is the set of all
ways those values appear.
"""
function tuples_in_trials(trials, n_way)
    param_cnt = length(trials[1])
    param_combo = [Tuple(x) for x in combinations(1:param_cnt, n_way)]
    # The dictionary with keys = parameter indices and dictionary
    # value is set of tuples of parameter values.
    seen = Dict(x => Set{NTuple{n_way,Int}}() for x in param_combo)
    for trial_idx in 1:length(trials)
        entry = trials[trial_idx]
        for params in param_combo
            values = entry[[x for x in params]]
            # values can be missing, in which case this isn't a tuple.
            if !any(values .== 0)
                push!(seen[params], Tuple(values))
            end
        end
    end
    seen
end

function coverage_by_tuple(trials, n_way)
    accumulated = tuples_in_trials(trials, n_way)
    sum([length(aset) for aset in values(accumulated)])
end


text:
#=
# Based on the file JacobiP.m from the book:
# Nodal Discontinuous Galerkin Method - Hesthaven, Warburton
#  https://link.springer.com/book/10.1007/978-0-387-72067-8
=#

# Evaluate Nth order Jacobi Polynomial of type (alpha,beta) at points x.
# Note   : They are normalized to be orthonormal.
function jacobi_polynomial(x, alpha::Int, beta::Int, N::Int)
    PL = zeros(N+1,length(x));

    # Initial values P_0(x) and P_1(x)
    gamma0 = 2^(alpha+beta+1)/(alpha+beta+1)*factorial(alpha)*
                factorial(beta)/factorial(alpha+beta);
    PL[1,:] = 1.0/sqrt(gamma0) .* ones(length(x));
    if N==0
        return PL[1,:];
    end
    gamma1 = (alpha+1)*(beta+1)/(alpha+beta+3) * gamma0;
    PL[2,:] = ((alpha+beta+2)/(2*sqrt(gamma1))) .* x .* ones(length(x)) .+ (alpha-beta)/(2*sqrt(gamma1));
    if N==1
        return PL[2,:];
    end

    # Higher orders are generated by this recurrence.
    # P(i+1) = -(a(i-1)/a(i))P(i-1) + ((x-b(i))/a(i))P(i)
    aold = 2/(2+alpha+beta)*sqrt((alpha+1)*(beta+1)/(alpha+beta+3));
    for i=1:N-1
      h1 = 2*i+alpha+beta;
      anew = 2/(h1+2)*sqrt( (i+1)*(i+1+alpha+beta)*(i+1+alpha)*
          (i+1+beta)/((h1+1)*(h1+3)));
      bnew = - (alpha^2-beta^2)/(h1*(h1+2));
      PL[i+2,:] = (-aold/anew) .* PL[i,:] + ((x.-bnew)./anew) .* PL[i+1,:];
      aold =anew;
    end;

    return PL[N+1,:];
end

# sqrt(N*(N+alpha+beta+1))*JacobiP(r(:),alpha+1,beta+1, N-1);
function grad_jacobi_polynomial(x, alpha::Int, beta::Int, N::Int)
    if N < 1
        return zeros(size(x));
    end

    return sqrt(N*(N+alpha+beta+1)) .* jacobi_polynomial(x, alpha+1, beta+1, N-1);
end


text:
immutable Geometric <: DiscreteUnivariateDistribution
    prob::Float64
    function Geometric(p::Real)
        zero(p) < p < one(p) || error("prob must be in (0, 1)")
    	new(float64(p))
    end
end

Geometric() = Geometric(0.5) # Flips of a fair coin

## Support
insupport(::Geometric, x::Real) = isinteger(x) && x >= zero(x)
insupport(::Type{Geometric}, x::Real) = isinteger(x) && x >= zero(x)

isupperbounded(d::Union(Geometric, Type{Geometric})) = false
islowerbounded(d::Union(Geometric, Type{Geometric})) = true
isbounded(d::Union(Geometric, Type{Geometric})) = false

minimum(d::Union(Geometric, Type{Geometric})) = 0
maximum(d::Geometric) = Inf

## Properties
mean(d::Geometric) = (1.0 - d.prob) / d.prob

median(d::Geometric) = -fld(logtwo,log1p(-d.prob)) - 1.0

mode(d::Geometric) = 0

var(d::Geometric) = (1.0 - d.prob) / d.prob^2
skewness(d::Geometric) = (2.0 - d.prob) / sqrt(1.0 - d.prob)
kurtosis(d::Geometric) = 6.0 + d.prob^2 / (1.0 - d.prob)

entropy(d::Geometric) = (-xlogx(1.0 - d.prob) - xlogx(d.prob)) / d.prob

## Functions
function pdf(d::Geometric, x::Real)
    !insupport(d,x) && return 0.0
    p = d.prob
    p < 0.5 ? p*exp(log1p(-p)*x) : p*(1.0-p)^x
end

logpdf(d::Geometric, x::Real) = insupport(d,x) ? log(d.prob) + log1p(-d.prob)*x : -Inf

function cdf(d::Geometric, x::Real)
    x < zero(x) && return 0.0
    p = d.prob
    n = floor(x) + 1.0
    p < 0.5 ? -expm1(log1p(-p)*n) : 1.0-(1.0-p)^n
end
function ccdf(d::Geometric, x::Real)
    x < zero(x) && return 1.0
    p = d.prob
    n = floor(x) + 1.0
    p < 0.5 ? exp(log1p(-p)*n) : (1.0-p)^n
end
logcdf(d::Geometric, q::Real) = q < zero(q) ? -Inf : log1mexp(log1p(-d.prob) * (floor(q) + 1.0))
logccdf(d::Geometric, q::Real) =  q < zero(q) ? 0.0 : log1p(-d.prob) * (floor(q) + 1.0)

quantile(d::Geometric, p::Real) = invlogccdf(d,log1p(-p))
cquantile(d::Geometric, p::Real) = invlogccdf(d,log(p))
invlogcdf(d::Geometric, lp::Real) = invlogccdf(d,log1mexp(lp))

function invlogccdf(d::Geometric, lp::Real)
    if (lp > zero(lp)) || isnan(lp)
        return NaN
    elseif isinf(lp)
        return Inf
    elseif lp == zero(lp)
        return 0.0
    end
    max(ceil(lp/log1p(-d.prob))-1.0,0.0)
end

function mgf(d::Geometric, t::Real)
    p = d.prob
    p / (expm1(-t) + p)
end

function cf(d::Geometric, t::Real)
    p = d.prob
    # replace with expm1 when complex version available
    p / (exp(-t*im) - 1.0 + p)
end

## Sampling
function rand(d::Geometric)
    e = Base.Random.randmtzig_exprnd()
    floor(-e/log1p(-d.prob))
end

## Fitting
immutable GeometricStats <: SufficientStats
    sx::Float64
    tw::Float64

    GeometricStats(sx::Real, tw::Real) = new(float64(sx), float64(tw))
end

suffstats{T<:Integer}(::Type{Geometric}, x::Array{T}) = GeometricStats(sum(x), length(x))

function suffstats{T<:Integer}(::Type{Geometric}, x::Array{T}, w::Array{Float64})
    n = length(x)
    if length(w) != n
        throw(ArgumentError("Inconsistent argument dimensions."))
    end
    sx = 0.
    tw = 0.
    for i = 1:n
        wi = w[i]
        sx += wi * x[i]
        tw += wi
    end
    GeometricStats(sx, tw)
end

fit_mle(::Type{Geometric}, ss::GeometricStats) = Geometric(1.0 / (ss.sx / ss.tw + 1.0))