text:
# # Cubical Persistent Homology

# In this example, we will demonstrate computing sublevel set persistent homology of time
# series and image data. We will need the following packages.

using Ripserer
using Images
using Plots
gr() # hide
using Random # hide
Random.seed!(1337) # hide
nothing # hide

# We will use two data sets in this example. The first will be a curve:

n = 1000
x = range(0, 1; length=n)
curve = sin.(2π * 5x) .* x

curve_plot = plot(curve; legend=false, title="Curve")

# The other will be [the Event Horizon Telescope picture of a black
# hole](https://en.wikipedia.org/wiki/File:Black_hole_-_Messier_87_crop_max_res.jpg). We
# will use a small, 240×240 pixel version of the image. Ripserer should have no problems
# with processing larger images, but this will work well enough for this tutorial.

blackhole_image = load(
    joinpath(@__DIR__, "../assets/data/240px-Black_hole_-_Messier_87_crop_max_res.jpg")
)
blackhole_plot = plot(blackhole_image; title="Black Hole")

# To use the image with Ripserer, we have to convert it to grayscale.

blackhole = Gray.(blackhole_image)
nothing # hide

# ## One-dimensional Case

# Sublevel set persistent homology provides a stable description of the critical points of a
# function. The zeroth persistent homology group ``H_0`` corresponds to its local minima. To
# compute this with Ripserer, we use cubical persistent homology. Note that there is no
# information in ``H_1``, since the function is one-dimensional, so we only grab the first
# part of the result.

result, _ = ripserer(Cubical(curve))
plot(curve_plot, plot(result))

# To see which minimum each interval corresponds to, we compute representatives by setting
# `reps=true`.

result, _ = ripserer(Cubical(curve); reps=true)

# The infinite interval's representative will include the whole function. To plot a
# representative, simply pass the interval along with the data to `plot`. An alternative
# way to plot the same thing is to get the representative of the interval first and then
# plot that in the same way. Let's use both ways in the following example.

infinite_interval = only(filter(!isfinite, result))
plt = plot(
    representative(infinite_interval),
    curve;
    legend=false,
    title="Representatives",
    seriestype=:path,
)

for interval in filter(isfinite, result)
    plot!(plt, interval, curve; seriestype=:path)
end

# To get the locations of the minima, extract the critical simplices from intervals. Since
# simplices act like collections of vertex indices, and zero-dimensional representatives
# have a single point each, we can use `only` to extract them. Cubical simplices store
# vertices as `CartesianIndex`, so we need to index into the x-values with them to plot
# them.

min_indices = [only(birth_simplex(int)) for int in result]
min_x = eachindex(curve)[min_indices]

scatter!(plt, min_x, curve[min_x]; color=1:6, markershape=:star)
plot(plt, plot(result; markercolor=1:6, markeralpha=1))

# Note that each interval's birth is equal to the value of the corresponding local minimum
# and its death is equal to the higher of the two adjacent maxima. An intuitive way of
# thinking about the result is imagining you pour water in the curve. Water is collected in
# a valley and once it reaches a local maximum, it starts pouring in the adjacent valley.

# ### Dealing With Noise

# Let's add some noise to the curve from earlier and try to repeat the example.

noisy_curve = curve .+ rand(length(curve))
plot(noisy_curve)

# Compute the persistence diagram and plot it. The `persistence` argument to `plot` plots
# the diagram in birth, persistence coordinates, which will make the level of noise easier
# to see.

noisy_result, _ = ripserer(Cubical(noisy_curve); reps=true)
plot(noisy_result; persistence=true)

# We can see a lot of noise in the diagram, but notice how all of it is below 1, which is
# exactly the amount we added. We can filter this out by supplying the `cutoff` argument to
# `ripserer`. This will ignore all intervals below this cutoff.

filtered_result, _ = ripserer(Cubical(noisy_curve); reps=true, cutoff=1)
plot(filtered_result; persistence=true)

# We could achieve the same result by simply filtering the resulting diagram.

@assert filter(int -> persistence(int) > 1, noisy_result) == filtered_result # hide
filter(int -> persistence(int) > 1, noisy_result) == filtered_result

# Keep in mind, however, that for very large diagrams, collecting representatives on all of
# the noise could make `ripserer` much run slower. In that case, it is recommended to use
# `cutoff`, rather than filtering.

# Now we can repeat what we did earlier.

infinite_interval = only(filter(!isfinite, filtered_result))
plt = plot(
    representative(infinite_interval),
    noisy_curve;
    legend=false,
    title="Representatives",
    seriestype=:path,
)

for interval in sort(filter(isfinite, filtered_result); by=birth)
    plot!(plt, interval, noisy_curve; seriestype=:path)
end

min_indices = [first(birth_simplex(int)) for int in filtered_result]
min_x = eachindex(curve)[min_indices]

scatter!(
    plt, min_x, noisy_curve[min_x]; color=eachindex(filtered_result), markershape=:star
)
plot(plt, plot(filtered_result; markercolor=eachindex(filtered_result), markeralpha=1))

# The result we got is very similar even though the data set was very noisy.

# ## Two-dimensional Case

# Now let's do a similar thing for a 2d example. Nothing is stopping us from going
# into higher dimensions, but we will skip those. The principles are the same.
# Instead of looking for local minima, let's look for local maxima. To do that, we have to
# invert the image.

result = ripserer(Cubical(-blackhole))
plot(blackhole_plot, plot(result))

# We notice there are quite a lot of intervals along the diagonal. These correspond to the
# local geometry of the image, so we are not interested in them right now. To filter them
# out, we set a cutoff.

result = ripserer(Cubical(-blackhole); cutoff=0.1)
plot(blackhole_plot, plot(result))

# Like earlier, we can show the local extrema in the image. We will show a different way to
# plot them. Let's use the `threshold` argument with `plot`, which only keeps parts of the
# representative with a diameter equal to or lower than `threshold`. If we needed a strict
# `<`, we could use `threshold_strict`.

result = ripserer(Cubical(-blackhole); cutoff=0.1, reps=true)
plt = plot(blackhole_image; title="Black Hole")
for interval in result[1]
    plot!(plt, interval, blackhole; threshold=birth(interval))
end
plot(plt, plot(result[1]; markercolor=2:3, markeralpha=1))

# Some maxima have multiple values because more than one pixel in the image has the same
# value. Using `birth_simplex` instead of plotting the filtered representative would only
# show one of each.

# Unlike in the previous example, we now also have access to ``H_1``, which corresponds to
# the cycles in the image. Let's try to plot the representative of ``H_1``.

plt = plot(blackhole_image; title="Black Hole")
plot!(plt, only(result[2]), blackhole; label="H₁ cocycle", color=1)

# Notice that the result is not a cycle, but rather a collection of pixels that
# would destroy the cycle if removed. The reason is that Ripserer computes persistent
# _co_homology by default. The persistence diagrams of persistent homology and persistent
# cohomology are the same, but persistent cohomology is much more efficient to compute. The
# representatives it finds, however, tend to not be as informative. Keep this in mind when
# trying persistent homology out for larger datasets; it can take a very very long
# time. This case is quite small and computing persistent homology should pose no problem.
# We compute persistent homology with the argument `alg=:homology`.

# !!! warning "Infinite Intervals in Persistent Homology"
#     Ripserer currently can't compute infinite intervals in dimensions
#     higher than zero with persistent homology.

result = ripserer(Cubical(-blackhole); cutoff=0.1, reps=true, alg=:homology)
plot!(plt, only(result[2]), blackhole; label="H₁ cycle", color=3)

# We have successfully found the hole in a black hole.


text:
Base.@kwdef struct DensityAnalysis{D, K, B}
    datalimits::D=automatic
    npoints::Int=200
    kernel::K=automatic
    bandwidth::B=automatic
end

# Work around lack of length 1 tuple method
# TODO: also add weights support here
_kde(data::NTuple{1, Any}; kwargs...) = kde(data...; kwargs...)
_kde(data::Tuple; kwargs...) = kde(data; kwargs...)

function compute_datalimits(positional, datalimits)
    return if datalimits === automatic
        map(v -> mapreduce(extrema, extend_extrema, v), Tuple(positional))
    else
        datalimits
    end
end

applydatalimits(f::Function, d) = map(f, d)
applydatalimits(limits::Union{AbstractArray, Tuple}, _) = limits

function _density(data...; datalimits, npoints, kwargs...)
    k = _kde(data; kwargs...)
    es = applydatalimits(datalimits, data)
    rgs = map(e -> range(e...; length=npoints), es)
    res = pdf(k, rgs...)
    return (rgs..., res)
end

function (d::DensityAnalysis)(input::ProcessedLayer)
    datalimits = compute_datalimits(input.positional, d.datalimits)
    options = valid_options(; datalimits, d.npoints, d.kernel, d.bandwidth)
    output = map(input) do p, n
        return _density(p...; pairs(n)..., pairs(options)...), (;)
    end
    N = length(input.positional)
    labels = set(input.labels, N+1 => "pdf")
    plottypes = [LinesFill, Heatmap, Volume]
    default_plottype = plottypes[N]
    plottype = Makie.plottype(input.plottype, default_plottype)
    return ProcessedLayer(output; plottype, labels)
end

"""
    density(; datalimits=automatic, kernel=automatic, bandwidth=automatic, npoints=200)

Fit a kernel density estimation of `data`.
Here, `datalimits` specifies the range for which the density should be calculated,
and `kernel` and `bandwidth` are forwarded to `KernelDensity.kde`.
`npoints` is the number of points used by Makie to draw the line
"""
density(; options...) = transformation(DensityAnalysis(; options...))

ext: jl


text:
"""
This is Tutorial_netCDF.jl. It contains all the necessary commands to import a netCDF file from a seismic tomography,
to convert it to the GeoData format and to export it to a Paraview format. The following steps are performed:
1. Read data from this file.
2. Put the data in a GeoData format (this is the format that is used internally in the GMG).
3. Export that data to a format readable by Paraview.

The data file used in this example can be downloaded here:
[https://ds.iris.edu/files/products/emc/emc-files/El-Sharkawy-etal-G3.2020-MeRE2020-Mediterranean-0.0.nc](https://ds.iris.edu/files/products/emc/emc-files/El-Sharkawy-etal-G3.2020-MeRE2020-Mediterranean-0.0.nc)
"""

using NetCDF, GeophysicalModelGenerator

# 1. define where the file is located on your computer
filename = "./El-Sharkawy-etal-G3.2020-MeRE2020-Mediterranean-0.0.nc" # this only works if you are in the directory where the data is located

# 2. load desired data

# Now check with ncinfo(filename), what the variables are called exactly and what the contents of your netCDF file are

lat = ncread(filename,"latitude");
lon = ncread(filename,"longitude");
depth = ncread(filename,"depth");
vs    = ncread(filename,"Vs");
depth = -1 .* depth # CAREFUL: MAKE SURE DEPTH IS NEGATIVE, AS THIS IS THE ASSUMTPION IN GeoData

# For netCDF data, 3D coordinates of a regular grid are only given as 1D vectors. As we need to compute Cartesian coordinates for
# Paraview, we need the full matrix of grid coordinates
Lon3D,Lat3D,Depth3D = LonLatDepthGrid(lon, lat, depth);

# Set up the Data structure
Data_set       =   GeoData(Lon3D,Lat3D,Depth3D,(VS=vs,))

# Export the data structure to Paraview format
Write_Paraview(Data_set, "test_netcdf_3D")


text:
function test_bcs()

  fname = "input_vals_channel_bc.jl"
  mesh, sbp, eqn, opts = run_solver(fname)
  params = eqn.params
  gamma = params.gamma
  gamma_1 = params.gamma_1

  #----------------------------------------------------------------------------
  # Subsonic Inflow
  # construct a state such that the BC function will compute the boundary state
  # equal to the input state

  Ma = 0.5
  # these must match the values in the BC function
  pt = 102010.0/params.p_free
  Tt = 288.6/params.T_free

  println("pt = ", pt)
  println("Tt = ", Tt)

#  Ub = 1.0  #TODO: make sure this is subsonic
#  a = Ub/Ma

  operand = 1 + 0.5*gamma_1*Ma*Ma
  p = pt*(operand^(-gamma/gamma_1))
  T = Tt*(operand^(-1))

  rho = p/(params.R_ND*T)
  a = sqrt(gamma*p/rho)
  Ub = Ma*a
  E = p/params.gamma_1 + 0.5*rho*Ub*Ub

  q = zeros(mesh.numDofPerNode)
  q[1] = rho
  q[2] = 0
  q[3] = rho*Ub
  q[4] = E

  F = zeros(mesh.numDofPerNode)
  F2 = zeros(mesh.numDofPerNode)
  aux_vars = Float64[EulerEquationMod.calcPressure(params, q)]
  coords = mesh.coords_interface[:, 1, 1]
#  nrm_xy = mesh.nrm_face[:, 1, 1]
  nrm_xy = -Float64[0.0, 2.0]

  EulerEquationMod.calcEulerFlux(params, q, aux_vars, nrm_xy, F)
  obj = EulerEquationMod.BCDict["subsonicInflowBC"](mesh, eqn)
  obj(params, q, aux_vars, coords, nrm_xy, F2)

  @test isapprox( norm((F - F2)/norm(q)), 0.0) atol=1e-12

  #----------------------------------------------------------------------------
  # Subsonic Outflow

  p = 101300.0/params.p_free # must match the value in the BC function
  q[1] = 2
  q[2] = 0.2
  q[3] = 0.2
  q[4] = p/gamma_1 + 0.5*( q[2]*q[2] + q[3]*q[3])/q[1]

  nrm_xy = Float64[0.0, 2.0]

  EulerEquationMod.calcEulerFlux(params, q, aux_vars, nrm_xy, F)
  obj = EulerEquationMod.BCDict["subsonicOutflowBC"](mesh, eqn)
  obj(params, q, aux_vars, coords, nrm_xy, F2)

  @test isapprox( norm((F - F2)/norm(q)), 0.0) atol=1e-12

  return nothing
end

add_func1!(EulerTests, test_bcs, [TAG_SHORTTEST])


text:
# ----------------------------------------------------------------------------
# This file was created by manually inspecting 2D Brillouin zones.
# ----------------------------------------------------------------------------

const pathsd_2d = Dict(
    :hp1 => [[:Γ, :M, :K, :Γ]],
    # ------------------------------------------------------------------------------------ #
    :tp1 => [[:Γ, :M, :X, :Γ]],
    # ------------------------------------------------------------------------------------ #
    :op1 => [[:Γ, :Y, :S, :X, :Γ]], # op1 and op2 are really the same paths; they differ
    :op2 => [[:Γ, :X, :S, :Y, :Γ]], # only in the definition of :X and :Y points (swapped)
    # ------------------------------------------------------------------------------------ #
    :oc1 => [[:C₀, :Y, :Γ, :Σ₀], [:Γ, :S]], # the k-point labels picked here are plane group
    :oc2 => [[:Δ₀, :Γ, :Y, :F₀], [:Γ, :S]], # 9's; they may not agree with plane group 5
    # TODO: think about if we can do better for plane group 5
    # ------------------------------------------------------------------------------------ #
    :mp1 => [[:Γ, :B], [:Γ, :Y], [:Γ, :A]], # the two settings have different definitions of
    :mp2 => [[:Γ, :B], [:Γ, :Y], [:Γ, :A]]  # :A [(1/2, +1/2) in mp1 and (1/2, -1/2) in mp2]
)

const pointsd_2d = Dict(
    :hp1 => [:Γ => :([0, 0]), :M => :([1/2, 0]), :K => :([1/3, 1/3])],
    # ------------------------------------------------------------------------------------ #
    :tp1 => [:Γ => :([0, 0]), :X => :([0, 1/2]), :M => :([1/2, 1/2])],
    # ------------------------------------------------------------------------------------ #
    :op1 => [:Γ => :([0, 0]), :X => :([0, 1/2]), :S => :([1/2, 1/2]), :Y => :([1/2, 0])],
    :op2 => [:Γ => :([0, 0]), :X => :([1/2, 0]), :S => :([1/2, 1/2]), :Y => :([0, 1/2])],
    # ------------------------------------------------------------------------------------ #
    :oc1 => [:Γ => :([0, 0]), :Y => :([-1/2, 1/2]), :C₀ => :([-δ, 1-δ]), :S => :([0, 1/2]), :Σ₀ => :([δ,  δ])],
    :oc2 => [:Γ => :([0, 0]), :Y => :([1/2, 1/2]),  :F₀ => :([δ,  1-δ]), :S => :([0, 1/2]), :Δ₀ => :([-δ, δ])],
    # ------------------------------------------------------------------------------------ #
    :mp1 => [:Γ => :([0, 0]), :B => :([1/2, 0]), :Y => :([0, 1/2]), :A => :([1/2, +1/2])],
    :mp2 => [:Γ => :([0, 0]), :B => :([1/2, 0]), :Y => :([0, 1/2]), :A => :([1/2, -1/2])]
)

const paramsd_2d = Dict(
   :oc1 => [:δ => :(0.25 * ( 1 + (a*a) / (b*b) ) )],
   :oc2 => [:δ => :(0.25 * ( 1 + (b*b) / (a*a) ) )]
)


text:
using LossFunctions, Plots
pyplot()

_mkpath(fname) = joinpath(dirname(@__FILE__()), fname)
srand(123)

# ---------------------------------------------------------------
# README Plots

xld  = "\$\\hat{y} - y\$"
xlm  = "\$y \\cdot \\hat{y}\$"
yld  = "\$L(\\hat{y} - y)\$"
yldd = "\$L'(\\hat{y} - y)\$"
ylm  = "\$L(y \\cdot \\hat{y})\$"
ylmd = "\$L'(y \\cdot \\hat{y})\$"
default(bg_outside=RGBA(1,1,1,0))

plot([L2DistLoss(), L1DistLoss(), LPDistLoss(.5), LogitDistLoss(),
      HuberLoss(), L1EpsilonInsLoss(1), L2EpsilonInsLoss(1),
      QuantileLoss(0.7)],
     -4:0.01:4,
     xlab=xld, ylab=yld,
     xlim=(-4,4), ylim=(0,5), palette=:darktest, grid=false,
     linestyle=[:dash :solid :solid :dashdot :dash :solid :solid :solid])
Plots.svg(_mkpath("distance"))

plot([ZeroOneLoss(), L1HingeLoss(),L2HingeLoss(), LogitMarginLoss(),
      ModifiedHuberLoss(), PerceptronLoss(), SmoothedL1HingeLoss(1),
      L2MarginLoss(), ExpLoss(), SigmoidLoss(), DWDMarginLoss(1)],
     -2:0.001:2,
     xlab=xlm, ylab=ylm,
     ylim=(0,7), palette=:darktest, grid=false,
     linestyle=[:dash :solid :solid :dashdot :dash :solid :solid :dash :solid :dot :dot])
Plots.svg(_mkpath("margin"))

# ---------------------------------------------------------------
# Readthedocs Distance Plots

default(bg_outside=RGBA(1,1,1,0), aspect_ratio = 1, palette=:darktest, size = (350, 350))

# LP Dist Loss
plot([LPDistLoss(0.5), LPDistLoss(1.5), LPDistLoss(3)],
     xlab=xld, ylab=yld,
     -2:0.01:2, ylim=(0,4))
Plots.svg(_mkpath("LPDistLoss1"))
plot([LPDistLoss(0.5)', LPDistLoss(1.5)', LPDistLoss(3)'],
     xlab=xld, ylab=yldd,
     -2:0.005:2, ylim=(-2,2))
Plots.svg(_mkpath("LPDistLoss2"))

# L1 Dist Loss
plot(L1DistLoss(), -2:0.05:2, ylim=(0,4), xlab=xld, ylab=yld)
Plots.svg(_mkpath("L1DistLoss1"))
plot(L1DistLoss()', -2:0.01:2, ylim=(-2,2), xlab=xld, ylab=yldd)
Plots.svg(_mkpath("L1DistLoss2"))

# L2 Dist Loss
plot(L2DistLoss(), -2:0.05:2, ylim=(0,4), xlab=xld, ylab=yld)
Plots.svg(_mkpath("L2DistLoss1"))
plot(L2DistLoss()', -2:0.05:2, ylim=(-2,2), xlab=xld, ylab=yldd)
Plots.svg(_mkpath("L2DistLoss2"))

# Logit Dist Loss
plot(LogitDistLoss(), -2:0.05:2, ylim=(0,4), xlab=xld, ylab=yld)
Plots.svg(_mkpath("LogitDistLoss1"))
plot(LogitDistLoss()', -2:0.05:2, ylim=(-2,2), xlab=xld, ylab=yldd)
Plots.svg(_mkpath("LogitDistLoss2"))

# HuberLoss
plot([HuberLoss(0.5), HuberLoss(1), HuberLoss(1.5)],
     xlab=xld, ylab=yld,
     -2:0.01:2, ylim=(0,4))
Plots.svg(_mkpath("HuberLoss1"))
plot([HuberLoss(0.5)', HuberLoss(1)', HuberLoss(1.5)'],
     xlab=xld, ylab=yldd,
     -2:0.005:2, ylim=(-2,2))
Plots.svg(_mkpath("HuberLoss2"))

# L1EpsilinInsLoss
plot([L1EpsilonInsLoss(0.5), L1EpsilonInsLoss(1), L1EpsilonInsLoss(1.5)],
     xlab=xld, ylab=yld,
     -2:0.01:2, ylim=(-0.05,4))
Plots.svg(_mkpath("L1EpsilonInsLoss1"))
plot([L1EpsilonInsLoss(0.5)', L1EpsilonInsLoss(1)', L1EpsilonInsLoss(1.5)'],
     xlab=xld, ylab=yldd,
     -2:0.005:2, ylim=(-2.05,2))
Plots.svg(_mkpath("L1EpsilonInsLoss2"))

# L2EpsilonInsLoss
plot([L2EpsilonInsLoss(0.5), L2EpsilonInsLoss(1), L2EpsilonInsLoss(1.5)],
     xlab=xld, ylab=yld,
     -2:0.01:2, ylim=(-0.05,4))
Plots.svg(_mkpath("L2EpsilonInsLoss1"))
plot([L2EpsilonInsLoss(0.5)', L2EpsilonInsLoss(1)', L2EpsilonInsLoss(1.5)'],
     xlab=xld, ylab=yldd,
     -2:0.005:2, ylim=(-2.05,2))
Plots.svg(_mkpath("L2EpsilonInsLoss2"))

# PeriodicLoss
plot([PeriodicLoss(.5), PeriodicLoss(1), PeriodicLoss(1.5)],
     xlab=xld, ylab=yld,
     -1:0.01:1, ylim=(0,2))
Plots.svg(_mkpath("PeriodicLoss1"))
plot([PeriodicLoss(.5)', PeriodicLoss(1)', PeriodicLoss(1.5)'],
     xlab=xld, ylab=yldd,
     -8:0.01:8, ylim = (-8,8))
Plots.svg(_mkpath("PeriodicLoss2"))

# QuantileLoss
plot([QuantileLoss(0.25), QuantileLoss(0.5), QuantileLoss(.75)],
     xlab=xld, ylab=yld,
     -2:0.01:2, ylim=(0,4))
Plots.svg(_mkpath("QuantileLoss1"))
plot([QuantileLoss(0.25)', QuantileLoss(0.5)', QuantileLoss(.75)'],
     xlab=xld, ylab=yldd,
     -2:0.005:2, ylim=(-2,2))
Plots.svg(_mkpath("QuantileLoss2"))

# ---------------------------------------------------------------
# Readthedocs Margin Plots

# Zero One Loss
plot(ZeroOneLoss(), -2:0.01:2, ylim=(-.05,4), xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("ZeroOneLoss1"))
plot(ZeroOneLoss()', -2:0.05:2, ylim=(-2.05,2), xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("ZeroOneLoss2"))

# Perceptron Loss
plot(PerceptronLoss(), -2:0.01:2, ylim=(-.05,4), xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("PerceptronLoss1"))
plot(PerceptronLoss()', -2:0.005:2, ylim=(-2.05,2), xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("PerceptronLoss2"))

# Hinge Loss
plot(L1HingeLoss(), -2:0.01:2, ylim=(-.05,4), xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("L1HingeLoss1"))
plot(L1HingeLoss()', -2:0.005:2, ylim=(-2.05,2), xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("L1HingeLoss2"))

# L2 Hinge Loss
plot(L2HingeLoss(), -2:0.01:2, ylim=(-.05,4), xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("L2HingeLoss1"))
plot(L2HingeLoss()', -2:0.005:2, ylim=(-2.05,2), xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("L2HingeLoss2"))

# Hinge Loss
plot([SmoothedL1HingeLoss(0.5),SmoothedL1HingeLoss(1.0),SmoothedL1HingeLoss(1.5),SmoothedL1HingeLoss(2.0)],
      -2:0.01:2, ylim=(-.05,4),
      xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("SmoothedL1HingeLoss1"))
plot([SmoothedL1HingeLoss(0.5)',SmoothedL1HingeLoss(1.0)',SmoothedL1HingeLoss(1.5)',SmoothedL1HingeLoss(2.0)'],
     -2:0.005:2, ylim=(-2.05,2),
     xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("SmoothedL1HingeLoss2"))

# Modified Huber Loss
plot(ModifiedHuberLoss(), -4:0.01:4, xticks=[-4:4...], yticks=[0:8...], ylim=(-.05,8), xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("ModifiedHuberLoss1"))
plot(ModifiedHuberLoss()', -4:0.005:4, xticks=[-4:4...], yticks=[-4:4...], ylim=(-4.05,4), xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("ModifiedHuberLoss2"))

# DWD Margin Loss
plot([DWDMarginLoss(0.5),DWDMarginLoss(1.0),DWDMarginLoss(1.5)],
      -2:0.01:2, ylim=(0,4),
      xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("DWDMarginLoss1"))
plot([DWDMarginLoss(0.5)',DWDMarginLoss(1.0)',DWDMarginLoss(1.5)'],
     -2:0.005:2, ylim=(-2,2),
     xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("DWDMarginLoss2"))

# L2 Margin Loss
plot(L2MarginLoss(), -2:0.01:2, ylim=(0,4), xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("L2MarginLoss1"))
plot(L2MarginLoss()', -2:0.005:2, ylim=(-2,2), xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("L2MarginLoss2"))

# Logit Margin Loss
plot(LogitMarginLoss(), -2:0.01:2, ylim=(0,4), xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("LogitMarginLoss1"))
plot(LogitMarginLoss()', -2:0.005:2, ylim=(-2,2), xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("LogitMarginLoss2"))

# Exp Loss
plot(ExpLoss(), -2:0.01:2, ylim=(0,4), xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("ExpLoss1"))
plot(ExpLoss()', -2:0.005:2, ylim=(-2,2), xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("ExpLoss2"))

# Sigmoid Loss
plot(SigmoidLoss(), -2:0.01:2, ylim=(0,4), xlab=xlm, ylab=ylm)
Plots.svg(_mkpath("SigmoidLoss1"))
plot(SigmoidLoss()', -2:0.005:2, ylim=(-2,2), xlab=xlm, ylab=ylmd)
Plots.svg(_mkpath("SigmoidLoss2"))

max_stars_repo_name: JuliaTagBot/FileStorage

max_issues_repo_name: JuliaTagBot/FileStorage

max_forks_repo_name: JuliaTagBot/FileStorage


text:
module Headings

import Genie, Stipple

export heading

function heading( title::String = "";
                  class::String = "",
                  img::Union{AbstractString,Nothing} = nothing,
                  h = Genie.Renderer.Html.h1(title, class="st-header__title text-h3"),
                  content::Union{AbstractString,Nothing,Vector{AbstractString}} = nothing,
                  kwargs...)

  Genie.Renderer.Html.header(; class="$(class) st-header q-pa-sm", kwargs...) do; [
    img
    h
    content
  ]end
end

end


text:
const SKVec = [
"AdaBoostClassifier",
"AdaBoostRegressor",
"BaggingClassifier",
"ColumnTransformer",
"DecisionTreeClassifier",
"DecisionTreeRegressor",
"DummyClassifier",
"DummyRegressor",
"ExtraTreesClassifier",
"ExtraTreesRegressor",
"FeatureAgglomeration",
"FunctionTransformer",
"GaussianNB",
"GradientBoostingClassifier",
"GradientBoostingRegressor",
"KNeighborsClassifier",
"KNeighborsRegressor",
"LinearRegression",
"LinearSVC",
"LogisticRegression",
"MLPClassifier",
"MinMaxScaler",
"MissingIndicator",
"MultinomialNB",
"NMF",
"Normalizer",
"Nystroem",
"OneHotEncoder",
"OrdinalEncoder",
"PCA",
"PassiveAggressiveClassifier",
"Pipeline",
"PolynomialFeatures",
"QuadraticDiscriminantAnalysis",
"QuantileTransformer",
"RFE",
"RandomForestClassifier",
"RandomForestRegressor",
"Ridge",
"RidgeClassifier",
"RobustScaler",
"SGDClassifier",
"SGDRegressor",
"SVC",
"SVR",
"SelectKBest",
"SimpleImputer",
"StandardScaler",
"TfidfVectorizer",
"VarianceThreshold",
"VotingClassifier"]


text:
ts_input_data = load_timeseries_data(:CEP_GER18)

 ev1 = SimpleExtremeValueDescr("wind-dena42","max","absolute") #idx 39
 ev2 = SimpleExtremeValueDescr("solar-dena42","min","integral") # idx 359
 ev3 = SimpleExtremeValueDescr("el_demand-dena21","max","integral") # idx 19
 ev4 = SimpleExtremeValueDescr("el_demand-dena21","min","absolute") # idx 185
 ev = [ev1, ev2, ev3]

# simple extreme day selection
@testset "simple extr day selection" begin
    @testset "single day" begin
        #max absolute"# max absolute
        ts_input_data_mod,extr_vals,extr_idcs = simple_extr_val_sel(ts_input_data,ev1;rep_mod_method="feasibility")
        test_ClustData(ts_input_data_mod, ts_input_data)
        for (k,v) in extr_vals.data
            @test all(extr_vals.data[k] .≈ ts_input_data.data[k][:,39])
        end
        @test extr_idcs == [39]

        # min integral
        ts_input_data_mod,extr_vals,extr_idcs = simple_extr_val_sel(ts_input_data,ev2;rep_mod_method="feasibility")
        test_ClustData(ts_input_data_mod, ts_input_data)
        for (k,v) in extr_vals.data
            @test all(extr_vals.data[k] .≈ ts_input_data.data[k][:,359])
        end
        @test extr_idcs == [359]

        # max integral
        ts_input_data_mod,extr_vals,extr_idcs = simple_extr_val_sel(ts_input_data,ev3;rep_mod_method="feasibility")
        test_ClustData(ts_input_data_mod, ts_input_data)
        for (k,v) in extr_vals.data
            @test all(extr_vals.data[k] .≈ ts_input_data.data[k][:,19])
        end
        @test extr_idcs == [19]

        # min absolute
        ts_input_data_mod,extr_vals,extr_idcs = simple_extr_val_sel(ts_input_data,ev4;rep_mod_method="feasibility")
        test_ClustData(ts_input_data_mod, ts_input_data)
        for (k,v) in extr_vals.data
            @test all(extr_vals.data[k] .≈ ts_input_data.data[k][:,185])
        end
        @test extr_idcs == [185]
    end

    @testset "multiple days" begin
        ts_input_data_mod,extr_vals,extr_idcs = simple_extr_val_sel(ts_input_data,ev;rep_mod_method="feasibility")
        test_ClustData(ts_input_data_mod, ts_input_data)
        for (k,v) in extr_vals.data
            @test all(extr_vals.data[k] .≈ ts_input_data.data[k][:,[39,359,19]])
        end
        @test extr_idcs == [39,359,19]

    end
end

@testset "representation modification" begin
     ev1 = SimpleExtremeValueDescr("wind-dena42","max","absolute") #idx 39
     ev2 = SimpleExtremeValueDescr("solar-dena42","min","integral") # idx 359
     ev3 = SimpleExtremeValueDescr("el_demand-dena21","max","integral") # idx 19
     ev = [ev1, ev2, ev3]
     mod_methods = ["feasibility","append"]
     @testset "$mod_method" for mod_method in mod_methods begin
         @testset "single day" begin
            ts_input_data_mod,extr_vals,extr_idcs = simple_extr_val_sel(ts_input_data,ev1;rep_mod_method=mod_method)
            for (k,v) in extr_vals.data
                @test all(extr_vals.data[k] .≈ ts_input_data.data[k][:,39])
            end
            if mod_method=="feasibility"
                @test all(extr_vals.weights .==0.)
            elseif mod_method=="append"
                @test all(extr_vals.weights .==1.)
            end
            @test extr_vals.T==24
            @test extr_vals.K==1
            @test extr_idcs == [39]

            ts_clust_res = run_clust(ts_input_data_mod;method="kmeans",representation="centroid",n_init=10,n_clust=5) # default k-means
            ts_clust_extr = representation_modification(extr_vals,ts_clust_res.clust_data)

            @test ts_clust_extr.T == ts_input_data.T
            @test ts_clust_extr.K == ts_clust_res.clust_data.K + 1
            for (k,v) in ts_clust_extr.data
                @test all(ts_clust_extr.data[k][:,1:ts_clust_res.clust_data.K] .≈ ts_clust_res.clust_data.data[k])
                @test all(ts_clust_extr.data[k][:,ts_clust_res.clust_data.K+1] .≈ extr_vals.data[k])
            end
            @test all(ts_clust_extr.weights[1:ts_clust_res.clust_data.K] .≈ ts_clust_res.clust_data.weights)
            @test all(ts_clust_extr.weights[ts_clust_res.clust_data.K+1] .≈ extr_vals.weights)
         end
         @testset "multiple days" begin
            ts_input_data_mod,extr_vals,extr_idcs = simple_extr_val_sel(ts_input_data,ev;rep_mod_method=mod_method)
            for (k,v) in extr_vals.data
                @test all(extr_vals.data[k] .≈ ts_input_data.data[k][:,[39,359,19]])
            end
            if mod_method=="feasibility"
                @test all(extr_vals.weights .==0.)
            elseif mod_method=="append"
                @test all(extr_vals.weights .==1.)
            end
            @test extr_vals.T==24
            @test extr_vals.K==3
            @test extr_idcs == [39,359,19]

            ts_clust_res = run_clust(ts_input_data_mod;method="kmeans",representation="centroid",n_init=10,n_clust=5) # default k-means
            ts_clust_extr = representation_modification(extr_vals,ts_clust_res.clust_data)

            @test ts_clust_extr.T == ts_input_data.T
            @test ts_clust_extr.K == ts_clust_res.clust_data.K + 3
            for (k,v) in ts_clust_extr.data
                @test all(ts_clust_extr.data[k][:,1:ts_clust_res.clust_data.K] .≈ ts_clust_res.clust_data.data[k])
                @test all(ts_clust_extr.data[k][:,ts_clust_res.clust_data.K+1:ts_clust_res.clust_data.K+3] .≈ extr_vals.data[k])
            end
            @test all(ts_clust_extr.weights[1:ts_clust_res.clust_data.K] .≈ ts_clust_res.clust_data.weights)
            @test all(ts_clust_extr.weights[ts_clust_res.clust_data.K+1:ts_clust_res.clust_data.K+3] .≈ extr_vals.weights)
         end

        end
    end

end


text:
using Test
using Logging

using Mozi

const PATH=pwd()

include("./parser.jl")

macro showbanner(word,total=99)
    n=length(word)
    m=(total-n)÷2
    for i in 1:m
        print("-")
    end
    print(word)
    for i in 1:total-m-n
        print("-")
    end
    println()
end

include("./sparse/sparse_test.jl")

include("./beams/element_test.jl")
include("./quads/element_test.jl")
# include("./trias/element_test.jl")

include("./static_cantilever.jl")
include("./dynamic_cantilever.jl")
include("./buckling_cantilever.jl")

# include("./static_quad.jl")
# include("./dynamic_quad.jl")
# include("./buckling_quad.jl")

# include("./static_tria.jl")
# include("./dynamic_tria.jl")
# include("./buckling_tria.jl")

# include("./static_p_delta.jl")
# include("./dynamic_p_delta.jl")
# include("./buckling_p_delta.jl")
#
# include("./static_large_deform.jl")
# include("./dynamic_large_deform.jl")
# include("./buckling_large_deform.jl")

# include("./dynamic_cantilever.jl")
# cdiff=result_nodal_time_history(assembly,"DIFF",2,0,3)
# newmark=result_nodal_time_history(assembly,"NEWMARK",2,0,3)
# wilson=result_nodal_time_history(assembly,"WILSON",2,0,3)
# modaldecomp=result_nodal_time_history(assembly,"MODALDECOMP",2,0,3)
# using PyPlot
# cla()
# plot(t,cdiff)
# plot(t,newmark)
# # plot(t,wilson)
# plot(t,modaldecomp)
# display(gcf())


text:
using NMEA

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

nmeas = NMEAData()
f = open("testdata.txt", "r")
for line = readlines(f)
    mtype = parse_msg!(nmeas, line)
    if (mtype == "GGA")
        @test nmeas.last_GGA.valid
    elseif (mtype == "RMC")
        @test nmeas.last_RMC.valid
    elseif (mtype == "GSA")
        @test nmeas.last_GSA.valid
    elseif (mtype == "GSV")
        @test nmeas.last_GSV.valid
    elseif (mtype == "GBS")
        @test nmeas.last_GBS.valid
    elseif (mtype == "VTG")
        @test nmeas.last_VTG.valid
    elseif (mtype == "GLL")
        @test nmeas.last_GLL.valid
    elseif (mtype == "ZDA")
        @test nmeas.last_ZDA.valid
    elseif (mtype == "DTM")
        @test nmeas.last_DTM.valid
    elseif (mtype == "PASHR")
        @test nmeas.last_PASHR.valid
    else
        continue
    end
end

example = NMEA.parse(raw"$GPGGA,134740.000,5540.3248,N,01231.2992,E,1,09,0.9,20.2,M,41.5,M,,0000*61")
@test(example.latitude == 55.67208)

example = NMEA.parse("\$GNGGA,033158.447,4532.5164,N,12258.1463,W,1,9,1.25,161.5,M,-19.6,M,,*7E")
@test(example.longitude == -122.969105)


text:
Base.log(x::SReal{X}) where X = SFloat64{log(X)}()

"""
# StaticValues Examples
```jldoctest
julia> prevpow(2, 7)
4
julia> prevpow(2, 9)
8
julia> prevpow(5, 20)
5
julia> prevpow(4, 16)
16
```
SUInt(0x0000000000000004)
"""
function Base.prevpow(a::SReal, x::SReal)
    x < 1 && throw(DomainError(x, "`x` must be ≥ 1."))
    # See comment in nextpos() for a == special case.
    a == 2 && isa(x, Integer) && return _prevpow2(x)
    a <= 1 && throw(DomainError(a, "`a` must be greater than 1."))
    n = floor(Integer, log(a, x))
    p = a^(n+SOne)
    p <= x ? p : a^n
end

_prevpow2(x::Unsigned) = one(x) << unsigned((SInt(sizeof(eltype(x))) << SInt(3))-leading_zeros(x)-SOne)
_prevpow2(x::Integer) = convert(typeof(x), x < 0 ? -_prevpow2(unsigned(-x)) :
                                                        _prevpow2(unsigned(x)))
#=
    # only iterate over <:Integer
    for (ST2,BT2) in zip(static_integer, base_integer)
        @eval begin
            Base.round(::Type{$BT2}, ::$ST{V}) where V = $ST2{round($BT2, V::$BT)}()
        end
    end
end

Base.isfinite(s::SNumber) = isfinite(values(s))

#=
FIXME:
Weird inference error on gcdx

a = SInt(9)
b = SInt(-9)

julia> @inferred(gcdx(a,b))
(9, 0, -1)

a> a = SInt(4)
4

julia> b = SInt(-8)
-8

julia> @inferred(gcdx(a,b))
(4, 1, 0)

a = SInt(6)
b = SInt(-9)
julia> @inferred(gcdx(a,b))
ERROR: return type Tuple{SInt64{3},SInt64{-1},SInt64{-1}} does not match inferred return type Tuple{SInt64{_A} where _A,SInt64{_A} where _A,SInt64{_A} where _A}
=#

for (ST,BT) in zip(static_integer, base_integer)
    @eval begin
        function Base.invmod(n::$ST, m::$ST)
            g, x, y = gcdx(n, m)
            g != 1 && throw(DomainError((n, m), "Greatest common divisor is $g."))
            m == 0 && throw(DomainError(m, "`m` must not be 0."))
            # Note that m might be negative here.
            # For unsigned T, x might be close to typemax; add m to force a wrap-around.
            r = mod(x + m, m)
            # The postcondition is: mod(r * n, m) == mod(T(1), m) && div(r, m) == 0
            r
        end
    end
end

>>(x::SInteger, y::SInteger) = SInteger(>>(values(x), values(y)))

nbits16 = SInt(cld(precision(Float16), 2))
nbits32 = SInt(cld(precision(Float32), 2))
nbits64 = SInt(cld(precision(Float64), 2))

# lack of specificity in base requires that these be more verbosely written out
Base.nbitslen(::Type{Float16}, l::SReal{V1}, f::SReal{V2}) where {V1,V2} =
    min(nbits16, Base.nbitslen(l, f))
Base.nbitslen(::Type{Float32}, l::SReal{V1}, f::SReal{V2}) where {V1,V2} =
    min(nbits32, Base.nbitslen(l, f))
Base.nbitslen(::Type{Float64}, l::SReal{V1}, f::SReal{V2}) where {V1,V2} =
    min(nbits64, Base.nbitslen(l, f))

Base.nbitslen(::Type{<:SFloat16}, l::SReal{V1}, f::SReal{V2}) where {V1,V2} =
    min(nbits16, Base.nbitslen(l, f))
Base.nbitslen(::Type{<:SFloat32}, l::SReal{V1}, f::SReal{V2}) where {V1,V2} =
    min(nbits32, Base.nbitslen(l, f))
Base.nbitslen(::Type{<:SFloat64}, l::SReal{V1}, f::SReal{V2}) where {V1,V2} =
    min(nbits64, Base.nbitslen(l, f))

Base.AbstractFloat(x::SInteger) = SFloat64(x)

=#
# TODO
# - ispow2
# - isqrt
# - factorial
# - binomial
# - all string things
# - >> is type unstable if not same exact type
#
# we get the following for free
# - nextpow
# - prevpow
# - ndigits


text:
module TestRecursiveBrim
   using Base.Test
   using Brim
   #=testdir = dirname(@__FILE__)=#
   @test eye(Int64, 50) |> partition_random |> recursive_brim! |> Q >= 0.90
   @test ones(Int64, (50, 50)) |> partition_random |> recursive_brim! |> Q <= 0.1
end


text:
using Distributions, Roots

function BlackOverK(moneyness, stdDev, callOrPut)
    d1 = log(moneyness) / stdDev + stdDev / 2.0
    d2 = d1 - stdDev
    norm = Normal()
    return callOrPut * (moneyness*cdf(norm,callOrPut*d1)-cdf(norm,callOrPut*d2))
end

function Black(strike, forward, sigma, T, callOrPut)
    nu = sigma*sqrt(T)
    if nu<1.0e-12   # assume zero
        return max(callOrPut*(forward-strike),0.0)  # intrinsic value
    end
    return strike * BlackOverK(forward/strike,nu,callOrPut)
end

function BlackVega(strike, forward, sigma, T)
    stdDev = sigma*sqrt(T)
    d1 = log(forward/strike) / stdDev + stdDev / 2.0
    norm = Normal()
    return forward * cdf(norm,d1) * sqrt(T)
end

function BlackImpliedVol(price, strike, forward, T, callOrPut)
    objective(sigma) = Black(strike, forward, sigma, T, callOrPut) - price
    return find_zero(objective,(0.01,1.0),Roots.Brent(),xatol=1.0e-8)
end

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
"""
assembler is the functions in charge of the construction of the matrix and RHS of the main FEM system
"""
function assembler( nodalCoords, elemNodalConnec, MEBIValsMat, MEBIVec, materialsData, elementsData, BCsData, currentSolution )

    nElems = size( elemNodalConnec, 1)
    nNodes = size( nodalCoords,     1)

    # The global stiffness equations are computed for the unknown displacements (free dofs)
    KG = zeros( 6*nNodes, 6*nNodes )

    Kmatrices = [] ;

    for elem in (1:nElems);

        # extract the MEBI entry
        elemMEBIEntry = MEBIVec[elem]

        elemMat = MEBIValsMat[ elemMEBIEntry, 1 ]
        elemEle = MEBIValsMat[ elemMEBIEntry, 2 ]

        print( "elemMebi entry", elemMEBIEntry )

        if elemMat > 0  # if it has a material

            nodesOfThisElem = elemNodalConnec[ elem ] ;
            dofsElem = nodes2dofs( nodesOfThisElem , 6 ) ;

            if cmp( elementsData[elemEle].type, "truss" ) == 0

                dofsElem = dofsElem[ [1,5,7,11]]

                E = materialsData[ MEBIValsMat[ elemMEBIEntry, 1 ] ].youngModulus
                A = elementsData[  MEBIValsMat[ elemMEBIEntry, 2 ] ].crossSectionArea

                nodesOfElemCoords = nodalCoords[ nodesOfThisElem, : ] ;
                L = norm( nodalCoords[ nodesOfThisElem[1], : ] - nodalCoords[ nodesOfThisElem[2], : ] );

                KGelem = linearTrussStiffnessMatrix( E, A, nodesOfElemCoords );

            elseif cmp( elementsData[elemEle].type, "tetrahedron" ) == 0

                dofsElem = dofsElem[ Vector(1:2:6*4) ]

                E = materialsData[ MEBIValsMat[ elemMEBIEntry, 1 ] ].youngModulus

                KGelem = linearTetraSolid( E );

            end

            KG[ dofsElem, dofsElem ] = KG[ dofsElem, dofsElem ] .+ KGelem

            push!( Kmatrices, KGelem )
        end
    end

    FG, neumDofs = assembleFext( nodalCoords, elemNodalConnec, MEBIValsMat, MEBIVec, BCsData )

    return KG, FG, Kmatrices, neumDofs
end

function assembleFext( nodalCoords, elemNodalConnec, MEBIValsMat, MEBIVec, BCsData )

    nNodes     = size( nodalCoords,     1)

    neumDofs = computeNeumDofs( nNodes )

    print("neumdofs", neumDofs)
    FG = zeros( 6*nNodes ) ;
    elementsBC = computeElemsByMEBI( MEBIValsMat[:,3], MEBIVec )
    for indBC in (1:length(BCsData))
        print("indBC", indBC)
        for elem in elementsBC[indBC]
            print("elem:", elem)
            nodesElem = elemNodalConnec[ elem ]
            dofs      = nodes2dofs( nodesElem, 6)

            if length( BCsData[indBC].NeumannNodalDOFs) >0
                dofsNeu = dofs[ BCsData[indBC].NeumannNodalDOFs ]
                FG[ dofsNeu ] = FG[ dofsNeu ] + BCsData[indBC].NeumannNodalVals
            end

            if length( BCsData[indBC].DirichletNodalDOFs) >0
                @assert sum( BCsData[indBC].DirichletNodalVals )==0 "only homogeneous dirichlet BC are allowed by now!"
                print("dofs:", dofs)
                neumDofs[ dofs[ BCsData[indBC].DirichletNodalDOFs ] ] .= 0
            end
        end
    end

    print("neumdofs", neumDofs)

    neumDofs = sort( unique(neumDofs) )
    neumDofs[1] == 0 && popfirst!(neumDofs)

    print("neumdofs", neumDofs)
    print("listo\n")

    return FG, neumDofs
end


text:
using Dates
using ArchGDAL
const AG=ArchGDAL
using Glob

using ESDL

"""
   getdate(x,reg = r"[0-9]{8}T[0-9]{6}", df = dateformat"yyyymmddTHHMMSS")
Return a DateTime object from a string where the time stamp is found by `reg`
and it is parsed according to the `df` dateformat.
"""
function getdate(x,reg = r"[0-9]{8}T[0-9]{6}", df = dateformat"yyyymmddTHHMMSS")
   m = match(reg,x).match
   date =DateTime(m,df)
end

"""
grouptimes(times, timediff=200000)
Group a sorted vector of time stamps into subgroups
where the difference between neighbouring elements are less than `timediff` milliseconds.
This returns the indices of the subgroups as a vector of vectors.
"""
function grouptimes(times, timediff=200000)
   @assert sort(times) == times
   group = [1]
   groups = [group]

   for i in 2:length(times)
      t = times[i]
      period = t - times[group[end]]
      if period.value < timediff
         push!(group, i)
      else
         push!(groups, [i])
         group = groups[end]
      end
   end
   return groups
end

#indir = "/home/crem_fe/Daten/bonds/jurua/tifs"
#pol="VV"
"""
   gdalcube(indir, pol)

Load the datasets in `indir` with a polarisation `pol` as a ESDLArray.
We assume, that `indir` is a folder with geotiffs in the same CRS which are mosaicked into timesteps and then stacked as a threedimensional array.

"""
function gdalcube(indir, pol)

filenames = glob("*$(pol)*.tif", indir)
dates = getdate.(filenames)
# Sort the dates and files by DateTime
p = sortperm(dates)
sdates = dates[p]
sfiles = filenames[p]
@show sfiles
# Put the dates which are 200 seconds apart into groups
groupinds = grouptimes(sdates, 200000)

datasets = AG.read.(sfiles)
datasetgroups = [datasets[group] for group in groupinds]
#We have to save the vrts because the usage of nested vrts is not working as a rasterdataset

temp = tempdir()
outpaths = [joinpath(temp, splitext(basename(sfiles[group][1]))[1] * ".vrt") for group in groupinds]
@show outpaths
@show length(outpaths), length(datasetgroups)
vrt_grouped = AG.unsafe_gdalbuildvrt.(datasetgroups)
AG.write.(vrt_grouped, outpaths)
vrt_grouped = AG.read.(outpaths)
vrt_vv = AG.unsafe_gdalbuildvrt(vrt_grouped, ["-separate"])
rvrt_vv = AG.RasterDataset(vrt_vv)
cube=ESDLArray(rvrt_vv)
#bandnames = AG.GDAL.gdalgetfilelist(vrt_vv.ptr)

# Set the timestaps from the bandnames as time axis
dates_grouped = [sdates[group[begin]] for group in groupinds]

taxis = RangeAxis(:Time, dates_grouped)
ESDL.Cubes.renameaxis!(cube, "Band"=>taxis)
return cube
end
#Tryout chunking the data

function getbandnames(dataset)
   meta = AG.metadata(dataset)
   meta_split = map(x->split(x, "="),meta)
   #str = collect(meta_split[1])
   #print(typeof(meta_split[1]))
   #by_func = x -> parse(Int,replace(match(r"_\d",x[1]).match,"_",""))
   filter!(x->occursin(r"^Band",x[1]),meta_split)
   by_func =  x-> parse(Int, replace(match(r"_(?<number>\d+)", x[1])[:number], "_" => ""))
   sort!(meta_split, by = by_func)
   #map(x->x[1],meta_split)
   map(x->x[2],meta_split)
end

"""
extracts the dates from the bandnames in format "yyyy-mm-dd"

Parameters:
-----------
    input: an array of the 'full' band names
    reg: a reqular expression to find a pattern matching the dates in a 'full' band name
    df: format of the date to be converted

Returns:
-------
    array with the dates
"""
function getdates(input::Array{T} where T<:AbstractString, reg = r"[0-9]{8}", df = dateformat"yyyymmdd")
            dates = map(x->Dates.Date(match(reg, x).match,df), input)
            #strdates = map(x-> Dates.format(x, "yyyy-mm-dd"), dates)
            #@show dates
end

#=
function getbandnames(dataset)
   metadata = AG.metadata(dataset)
   metaarray = split.(bandnames, "=")
   bandnames = String[]
   for meta in metaarray
      if "Band" in meta[1]
         bandnumber = split(meta[1], "_")[2]
         push!()

=#


text:
"""
    ScalarizeBridge{T, F, S}

Transforms a constraint `AbstractVectorFunction`-in-`vector_set_type(S)` where
`S <: LPCone{T}` to `F`-in-`S`.
"""
mutable struct ScalarizeBridge{T,F,S} <: AbstractBridge
    scalar_constraints::Vector{MOI.ConstraintIndex{F,S}}
    constants::Vector{T}
end

function bridge_constraint(
    ::Type{ScalarizeBridge{T,F,S}},
    model::MOI.ModelLike,
    f::MOI.AbstractVectorFunction,
    set::MOI.Utilities.VectorLinearSet,
) where {T,F,S}
    dimension = MOI.output_dimension(f)
    constants = MOI.constant(f, T)
    new_f = MOI.Utilities.scalarize(f, true)
    constraints = Vector{MOI.ConstraintIndex{F,S}}(undef, dimension)
    for i in 1:dimension
        constraints[i] = MOI.add_constraint(model, new_f[i], S(-constants[i]))
    end
    return ScalarizeBridge{T,F,S}(constraints, constants)
end

function MOI.supports_constraint(
    ::Type{ScalarizeBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:MOI.Utilities.VectorLinearSet},
) where {T}
    # If `F` is `MOI.VectorAffineFunction{Complex{Float64}}`, `S` is `MOI.Zeros`
    # and `T` is `Float64`, it would create a set `MOI.EqualTo{Float64}` which
    # is incorrect hence we say we only support it if the coefficient type of
    # `F` is `T`.
    return MOI.Utilities.is_coefficient_type(F, T)
end

function MOI.Bridges.added_constrained_variable_types(::Type{<:ScalarizeBridge})
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(
    ::Type{ScalarizeBridge{T,F,S}},
) where {T,F,S}
    return Tuple{Type,Type}[(F, S)]
end

function concrete_bridge_type(
    ::Type{<:ScalarizeBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    S::Type{<:MOI.Utilities.VectorLinearSet},
) where {T}
    return ScalarizeBridge{
        T,
        MOI.Utilities.scalar_type(F),
        MOI.Utilities.scalar_set_type(S, T),
    }
end

function MOI.get(
    bridge::ScalarizeBridge{T,F,S},
    ::MOI.NumberOfConstraints{F,S},
)::Int64 where {T,F,S}
    return length(bridge.scalar_constraints)
end

function MOI.get(
    bridge::ScalarizeBridge{T,F,S},
    ::MOI.ListOfConstraintIndices{F,S},
) where {T,F,S}
    return copy(bridge.scalar_constraints)
end

function MOI.delete(model::MOI.ModelLike, bridge::ScalarizeBridge)
    for ci in bridge.scalar_constraints
        MOI.delete(model, ci)
    end
    return
end

function MOI.delete(
    model::MOI.ModelLike,
    bridge::ScalarizeBridge,
    i::MOI.Bridges.IndexInVector,
)
    MOI.delete(model, bridge.scalar_constraints[i.value])
    deleteat!(bridge.scalar_constraints, i.value)
    deleteat!(bridge.constants, i.value)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::ScalarizeBridge{T},
) where {T}
    func = MOI.Utilities.vectorize(
        MOI.get.(model, attr, bridge.scalar_constraints),
    )
    if !(func isa MOI.VectorOfVariables)
        # `func` is in terms of bridged variables here while in
        # `bridge_constraint` it was in terms of the solver variables so
        # `MOI.constant(set)` might be different than `bridge.constants[i]`.
        for i in eachindex(bridge.scalar_constraints)
            set = MOI.get(
                model,
                MOI.ConstraintSet(),
                bridge.scalar_constraints[i],
            )
            func = MOI.Utilities.operate_output_index!(
                -,
                T,
                i,
                func,
                MOI.constant(set),
            )
        end
    end
    return func
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::ScalarizeBridge{T,F,S},
) where {T,F,S}
    return MOI.Utilities.vector_set_type(S)(length(bridge.scalar_constraints))
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{ScalarizeBridge{T,F,S}},
) where {T,F,S}
    return MOI.supports(model, attr, MOI.ConstraintIndex{F,S})
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimal,
    bridge::ScalarizeBridge,
)
    return MOI.get.(model, attr, bridge.scalar_constraints) .+ bridge.constants
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::ScalarizeBridge,
)
    values = MOI.get.(model, attr, bridge.scalar_constraints)
    if any(value -> value === nothing, values)
        return
    end
    return values .+ bridge.constants
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::ScalarizeBridge,
    ::Nothing,
)
    for ci in bridge.scalar_constraints
        MOI.set(model, attr, ci, nothing)
    end
    return
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::ScalarizeBridge,
    value,
)
    MOI.set.(model, attr, bridge.scalar_constraints, value .- bridge.constants)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDual,
    bridge::ScalarizeBridge,
)
    return MOI.get.(model, attr, bridge.scalar_constraints)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDualStart,
    bridge::ScalarizeBridge,
)
    values = MOI.get.(model, attr, bridge.scalar_constraints)
    if any(value -> value === nothing, values)
        return
    end
    return values
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDualStart,
    bridge::ScalarizeBridge,
    value,
)
    MOI.set.(model, attr, bridge.scalar_constraints, value)
    return
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::ScalarizeBridge{T,F,S},
    change::MOI.VectorConstantChange{T},
) where {T,F,S}
    bridge.constants = change.new_constant
    MOI.set.(
        model,
        MOI.ConstraintSet(),
        bridge.scalar_constraints,
        S.(-change.new_constant),
    )
    return
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::ScalarizeBridge,
    change::MOI.MultirowChange{T},
) where {T}
    for (index, value) in change.new_coefficients
        MOI.modify(
            model,
            bridge.scalar_constraints[index],
            MOI.ScalarCoefficientChange{T}(change.variable, value),
        )
    end
    return
end

function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintFunction,
    bridge::ScalarizeBridge{T,F,S},
    func,
) where {T,F,S}
    old_constants = bridge.constants
    bridge.constants = MOI.constant(func, T)
    new_func = MOI.Utilities.scalarize(func, true)
    MOI.set.(
        model,
        MOI.ConstraintFunction(),
        bridge.scalar_constraints,
        new_func,
    )
    for i in eachindex(bridge.constants)
        if bridge.constants[i] != old_constants[i]
            MOI.set(
                model,
                MOI.ConstraintSet(),
                bridge.scalar_constraints[i],
                S(-bridge.constants[i]),
            )
        end
    end
    return
end

# TODO implement transform


text:
mutable struct MCVBCallback <: AbstractCallback
  every::Int64
  dkl::Float64
  aval::Float64
  # biasing observables
  Afunc::Function
  #Bfunc::Function
  # force Malliavin weights
  mvy::Array{Float64}
  mvydot::Array{Float64}
  gradf::Array{Float64}
  # value baseline Malliavin weights
  mvz::Array{Float64}
  mvzdot::Array{Float64}
  gradv::Array{Float64}
  # value baseline function
  vbl::AbstractValueBaseline
  # extra memory for computing some stuff
  em::Array{Float64}
  # IO stuff for writing things to files
  #avgsio::Union{IO,nothing}
  function MCVBCallback(Afunc::Function,model::AbstractPotential,vbl::AbstractValueBaseline)
    mvy::Array{Float64} = zeros(Float64,gradshape(model))
    mvydot::Array{Float64} = zeros(Float64,jacshape(model))
    gradf::Array{Float64} = zeros(Float64,gradshape(model))
    mvz::Array{Float64} = zeros(Float64,gradshape(vbl))
    mvzdot::Array{Float64} = zeros(Float64,gradshape(vbl))
    gradv::Array{Float64} = zeros(Float64,gradshape(vbl))
    em = zeros(Float64,getdimensionality(model))
    new(1,0.0,0.0,Afunc,mvy,mvydot,gradf,mvz,mvzdot,gradv,vbl,em)
  end
end

"""
Function that zeros all the quantities for restarting trajectories
"""
function initialize!(cb::MCVBCallback)
  cb.mvy .= 0.0
  cb.mvydot .= 0.0
  cb.mvz .= 0.0
  cb.mvzdot .= 0.0
end

function initializeavgs!(cb::MCVBCallback)
  cb.dkl = 0.0
  cb.aval = 0.0
  cb.gradf .= 0.0
  cb.gradv .= 0.0
end

function initializeall!(cb::MCVBCallback)
  cb.dkl = 0.0
  cb.aval = 0.0
  cb.mvy .= 0.0
  cb.mvydot .= 0.0
  cb.gradf .= 0.0
  cb.mvz .= 0.0
  cb.mvzdot .= 0.0
  cb.gradv .= 0.0
end

"""
Function that averages all the quantities by the number of trajectories
"""
function average!(ntraj::Int64,cb::MCVBCallback)
  cb.dkl /= ntraj
  cb.aval /= ntraj
  cb.gradf /= ntraj
  cb.gradv /= ntraj
end

"""
Implementation of the Monte Carlo Value Baseline algorithm
"""
function callback(cb::MCVBCallback,system::AbstractThermostattedSystem,mm::MixedModel,args...)
  # extract some arguments
  ind::Int64 = args[1]
  dt::Float64 = args[2]

  # compute gradient of variable force and update ydot
  jacobian!(cb.mvydot,system,mm.potentials[length(mm.potentials)])
  cb.mvydot .*= (sqrt(dt)/system.thermostat.scale)
  cb.mvy .+= (cb.mvydot*system.thermostat.rands)

  # calculate value baseline function
  vval::Float64 = callvbl(system,cb.vbl)

  # compute gradient of value baseline
  gradient!(cb.mvzdot,system,cb.vbl)
  cb.mvz .+= (dt .* cb.mvzdot)

  ### compute the instantaneous return
  # action difference - positive term
  cb.em .= ((system.thermostat.scale/sqrt(dt)) .* system.thermostat.rands).^2
  # action difference - negative term
  cb.em .-= (mm.potentials[length(mm.potentials)].f .+ (system.thermostat.scale/sqrt(dt)) .* system.thermostat.rands).^2
  # compute return
  rval::Float64 = sum(cb.em) / (2*system.thermostat.scale^2)

  # add return to KL divergence
  cb.dkl += rval*dt

  # observable bias
  aval::Float64 = cb.Afunc(system,mm,ind,dt)
  rval += aval# + Bfunc(system,mm,ind,dt)
  cb.aval += aval

  # update the gradients of the force parameters and the value baseline parameters
  cb.gradf .+= (rval .* cb.mvy) .- ((vval/dt) .* (cb.mvydot*system.thermostat.rands))
  cb.gradv .+= (rval .* cb.mvz) .- (vval .* cb.mvzdot)
end

#"""
#Implementation of the Monte Carlo Value Baseline algorithm for an active particle simulation
#"""
#function callback(cb::MCVBCallback,system::AbstractActiveSystem,mm::MixedModel,args...)
#  # extract some arguments
#  ind::Int64 = args[1]
#  dt::Float64 = args[2]
#
#  # compute gradient of variable force and update ydot
#  jacobian!(cb.mvydot,system,mm.potentials[length(mm.potentials)])
#  for i in 1:system.dim
#    @. @views cb.mvydot[:,i] .*= (sqrt(dt) / system.thermostat.scale[i])
#  end
#  cb.mvy .+= (cb.mvydot*system.thermostat.rands)
#
#  # calculate value baseline function
#  vval::Float64 = callvbl(system,cb.vbl)
#
#  # compute gradient of value baseline
#  gradient!(cb.mvzdot,system,cb.vbl)
#  cb.mvz .+= (dt .* cb.mvzdot)
#
#  ### compute the instantaneous return
#  # action difference - positive term
#  #cb.em .= ((system.thermostat.scale/sqrt(dt)) .* system.thermostat.rands).^2
#  # action difference - negative term
#  cb.em .= -(mm.potentials[length(mm.potentials)].f .+ (system.thermostat.scale/sqrt(dt)) .* system.thermostat.rands).^2
#  #cb.em .-= (mm.potentials[length(mm.potentials)].f .+ (system.thermostat.scale/sqrt(dt)) .* system.thermostat.rands).^2
#  # compute return
#  rval::Float64 = sum(cb.em ./ (2 .* system.thermostat.scale.^2))
#
#  # add return to KL divergence
#  cb.dkl += rval*dt
#
#  # observable bias
#  aval::Float64 = cb.Afunc(system,mm,ind,dt)
#  rval += aval# + Bfunc(system,mm,ind,dt)
#  cb.aval += aval
#
#  # update the gradients of the force parameters and the value baseline parameters
#  cb.gradf .+= (rval .* cb.mvy) .- ((vval/dt) .* (cb.mvydot*system.thermostat.rands))
#  cb.gradv .+= (rval .* cb.mvz) .- (vval .* cb.mvzdot)
#
#end
#


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Solving Scalar Equations
#
# - https://docs.sciml.ai/v5.0.0/tutorials/ode_example.html#Example-1-:-Solving-Scalar-Equations-1

using DifferentialEquations
using Plots

# ## simple form

# +
f(u,p,t) = 1.01*u
u0=1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")
# -

# ## a little generic form

# +
# function form and `du` argument is added
function f(du, u, p, t)
    du[1] = 1.01 * u[1]
end

u0 = [1 / 2] # 1/2 -> [1/2] scalar to arrays
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")
# -

# # Lotka-Voltera equation
#
# - Code and other topics is taken from https://julialang.org/blog/2019/01/fluxdiffeq/#what_is_the_neural_ordinary_differential_equation_ode
#
# $$
# \begin{align}
# x' &= \alpha x + \beta xy \\
# y' &= -\delta y + \gamma xy
# \end{align}
# $$

# +
using DifferentialEquations

# define Differenial Equation
function lotka_volterra(du, u, p, t)
    x, y = u
    α, β, δ, γ = p
    du[1] = α*x - β*x*y
    du[2] = -δ*y + γ*x*y
end

# set initial
u0 = [1.0, 1.0]

tspan = (0.0, 10.0)
p = [1.5, 1.0, 4.0, 1.0]

# set ODE solver
prob = ODEProblem(lotka_volterra, u0, tspan, p)

sol = solve(prob)
using Plots
plot(sol)


text:
using PyPlot

function circle(p; radius=1, scale=1, dth=0.1)
  scale*hcat(map(x->p+[radius*cos(x);radius*sin(x)], 0:dth:(2*pi)+dth)...)
end

function plot_circle(p; scale=1, radius=1, color="k", linestyle="-", linewidth=1.0)
  c = circle(p[1:2], scale=scale, radius=radius)
  plot3D(c[1,:][:], c[2,:][:], p[3]*ones(size(c,2)), color=color, linestyle=linestyle, linewidth=linewidth)
end

so2(theta) = [cos(theta) -sin(theta); sin(theta) cos(theta)]

function fill3d(ps; color = "k", alpha=1.0)
  ax = gca()
  tuples = map(x->tuple(ps[:,x]...), 1:size(ps,2))
  poly = art3D[:Poly3DCollection](Array[tuples])
  poly[:set_color](color)
  poly[:set_alpha](alpha)

  ax[:add_collection3d](poly)
end

function plot_quadrotor(p; scale=1, color="k", linewidth=2.0, theta=0.0, alpha=1.0)
  r = 0.6
  rot = so2(theta)
  circle_ps = rot*scale*[[1;0] [0;1] [-1;0] [0;-1]]

  p3 = 0.0
  if length(p) == 3
    p3 = p[3]
  end

  p2 = p[1:2]

  l1 = [p2+scale*rot[:,1] p2-scale*rot[:,1]]
  l2 = [p2+scale*rot[:,2] p2-scale*rot[:,2]]
  plot3D(l1[1,:][:], l1[2,:][:], p3*ones(2), color=color, linewidth=linewidth,
          alpha=alpha)
  plot3D(l2[1,:][:], l2[2,:][:], p3*ones(2), color=color, linewidth=linewidth,
          alpha = alpha)

  for ii = 1:size(circle_ps,2)
    c = circle(p2+circle_ps[:,ii], radius=scale*r)
    plot3D(c[1,:][:], c[2,:][:], p3*ones(size(c,2)), color=color,
            linestyle="--", linewidth=linewidth, alpha=alpha)
    #fill(c[1,:][:], c[2,:][:], p3*ones(size(c,2)), color=color, alpha=0.5)
    fill3d([c;p3*ones(1,size(c,2))], color=color, alpha=alpha*0.5)
  end
end

# iterations in 2, trials in 1
function standard_error(data, color = "k")
  stds = std(data, 1)[:] / sqrt(size(data,1))
  means = mean(data,1)[:]

  indices = collect(1:size(data,2))
  fill([indices;reverse(indices)], [means+stds;reverse(means-stds)],
    color = color, alpha=0.2, linewidth=0.0)
end

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
#!/usr/bin/env julia

# this script runs each runprofile.* script, with different values for each argument,
# collecting GPU and CUDA API traces in CSV files for further processing.

using Glob
using DataFrames

# run a benchmark once, returning the measurement data
function run_nvprof(callback, cmd, args)
    output_file = joinpath(tempdir(), "nvprof.csv")
    output_pattern = "nvprof.csv.*"

    # delete old profile output
    rm.(glob(output_pattern, tempdir()))

    # run and measure
    out = Pipe()
    nvprof = ```
        nvprof
        $args
        --profile-from-start off
        --concurrent-kernels off
        --profile-child-processes
        --unified-memory-profiling off
        --normalized-time-unit us
        --csv
        --log-file $output_file.%p
        $cmd
    ```
    cmd_success = success(pipeline(ignorestatus(nvprof), stdout=out, stderr=out))
    close(out.in)
    output_files = glob(output_pattern, tempdir())

    # read all data
    if !cmd_success
        println(readstring(out))
        error("benchmark did not succeed")
    else
        # when precompiling, an additional process will be spawned,
        # but the resulting CSV should not contain any data.
        output_data = []
        for output_file in output_files
            data = callback(output_file)
            if data != nothing
                push!(output_data, data)
            end
        end

        if length(output_data) == 0
            error("no output files")
        elseif length(output_data) > 1
            error("too many output files")
        else
            return output_data[1]
        end
    end
end

function save(path, data)
    open(path, "w") do io
        serialize(io, data)
    end
end

function collect_trace(cmd)
    @info "Collecting trace" cmd
    return run_nvprof(cmd, ```
        --print-api-trace
        --print-gpu-trace
    ```) do input_path
        table = mktemp() do temp_path,io
            for (line, contents) in enumerate(eachline(input_path; chomp=false))
                contains(contents, "No kernels were profiled.") && return nothing
                line < 4 && continue                        # skip nvprof banner
                line == 5 && continue                       # skip units header
                contains(contents, "cuEvent") && continue   # skip verbose TF event handling
                write(io, contents)
            end
            flush(io)
            readtable(temp_path)
        end

        return table
    end
end

function collect_metrics(cmd)
    @info "Collecting metrics" cmd
    return run_nvprof(cmd, ```
        --metrics achieved_occupancy,branch_efficiency,warp_execution_efficiency,warp_nonpred_execution_efficiency
    ```) do input_path
        table = mktemp() do temp_path,io
            for (line, contents) in enumerate(eachline(input_path; chomp=false))
                contains(contents, "No events/metrics were profiled.") && return nothing
                startswith(contents, "==") && continue                  # skip nvprof banner
                write(io, contents)
            end
            flush(io)
            readtable(temp_path)
        end

        return table
    end
end

function collect(cmd_trace, cmd_metrics, tag, dir)
    let fn = joinpath(dir, "$(tag)_trace.jls")
        if !isfile(fn)
            data = collect_trace(cmd_trace)
            save(fn, data)
        end
    end

    let fn = joinpath(dir, "$(tag)_metrics.jls")
        if !isfile(fn)
            data = collect_metrics(cmd_metrics)
            save(fn, data)
        end
    end
end

const ITERATIONS = length(ARGS) >= 1 ? parse(Int, ARGS[1]) : 128

dir = if length(ARGS) >= 2
    abspath(ARGS[2])
else
    @__DIR__
end
mkpath(dir)

cd(@__DIR__) do
    for problem_size in [512,1024,2048]
        let
            collect(`python3 runprofile.py $problem_size $ITERATIONS`,
                    `python3 runprofile.py $problem_size 1`,
                    "python_$(problem_size)",
                    dir)
        end

        for tfstyle in [true, false], uniform in [true, false]
            # uniform control doesn't impact tfstyle/non-fused execution
            (tfstyle && uniform) && continue

            collect(`julia --depwarn=no runprofile.jl $tfstyle $problem_size $ITERATIONS 0 $uniform`,
                    `julia --depwarn=no runprofile.jl $tfstyle $problem_size 1 0 $uniform`,
                    "julia_$(tfstyle ? "tfstyle" : "fused")_$(uniform ? "uniform" : "random")_$(problem_size)",
                    dir)
        end

        for arity in 1:10
            collect(`julia --depwarn=no runprofile.jl false $problem_size $ITERATIONS $arity`,
                    `julia --depwarn=no runprofile.jl false $problem_size 1 $arity`,
                    "julia_arity$(arity)_$(problem_size)",
                    dir)
        end
    end
end


text:
import Base: setindex!, getindex

include("params.jl")

"""
    LindaEnv

Optimizer environment. Holds parameters' values.
"""
mutable struct LindaEnv

    num_cols_rmp_max::LindaIntParam  # Maximum number of columns in RMP

    num_columns_max::LindaIntParam  # Maximum number of columns before pricing stops

    num_cgiter_max::LindaIntParam  # Maximum number of Column-Generation iterations

    time_limit::LindaFloatParam  # Time limit, in seconds

    tol_reduced_cost::LindaFloatParam  # Numerical tolerance for reduced costs

    verbose::LindaIntParam  # Verbosity level

    function LindaEnv()
        env = new()

        env.num_cols_rmp_max = LindaRealParam(:num_cols_rmp_max,
            typemax(Int64), zero(Int64), typemax(Int64)
        )

        env.num_columns_max = LindaRealParam(:num_columns_max,
            typemax(Int64), zero(Int64), typemax(Int64)
        )

        env.num_cgiter_max = LindaRealParam(:num_cgiter_max,
            typemax(Int64), zero(Int64), typemax(Int64)
        )

        env.time_limit = LindaRealParam(:time_limit, Inf, 0.0, Inf)

        env.tol_reduced_cost = LindaRealParam(:tol_reduced_cost, 10.0^-6, 0.0, Inf)

        env.verbose = LindaRealParam(:verbose, 0, 0, 1 )

        return env
    end

end

# TODO: use `getproperty` and `setproperty!` instead (> Julia 0.7)
getindex(env::LindaEnv, name::Symbol) = get_param_value(Core.getfield(env, name))

getindex(env::LindaEnv, ::Type{Val{:num_cols_rmp_max}}) = env.num_cols_rmp_max.val
getindex(env::LindaEnv, ::Type{Val{:num_columns_max}}) = env.num_columns_max.val
getindex(env::LindaEnv, ::Type{Val{:num_cgiter_max}}) = env.num_cgiter_max.val
getindex(env::LindaEnv, ::Type{Val{:time_limit}}) = env.time_limit.val
getindex(env::LindaEnv, ::Type{Val{:tol_reduced_cost}}) = env.tol_reduced_cost.val
getindex(env::LindaEnv, ::Type{Val{:verbose}}) = env.verbose.val

setindex!(env::LindaEnv, v, name::Symbol) = set_param_value!(Core.getfield(env, name), v)

"""
    setparam!(env; kwargs...)

Set parameters to specified values.
"""
function setparam!(env; kwargs...)
    for (s, v) in kwargs
        set_param_value!(Core.getfield(env, s), v)
    end
end

"""
    reset!(env)

Reset all parameters to their default values.
"""
function reset!(env)
    for s in fieldnames(LindaEnv)
        set_param_default!(Core.getfield(env, s))
    end
    return nothing
end


text:
macro c(ret_type, func, arg_types, lib)
  local args_in = Any[ symbol(string('a',x)) for x in 1:length(arg_types.args) ]
  quote
    $(esc(func))($(args_in...)) = ccall( ($(string(func)), $(esc(lib))), $ret_type, $arg_types, $(args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  quote
    typealias $fake_t $real_t
  end
end

@ctypedef XML_Char Uint8
@ctypedef XML_LChar Uint8
@ctypedef XML_Index Int32
@ctypedef XML_Size Uint32
@ctypedef XML_Parser Ptr{Void}
@ctypedef XML_Bool Uint8
# enum XML_Status
const XML_STATUS_ERROR = 0
const XML_STATUS_OK = 1
const XML_STATUS_SUSPENDED = 2
# end
# enum XML_Error
const XML_ERROR_NONE = 0
const XML_ERROR_NO_MEMORY = 1
const XML_ERROR_SYNTAX = 2
const XML_ERROR_NO_ELEMENTS = 3
const XML_ERROR_INVALID_TOKEN = 4
const XML_ERROR_UNCLOSED_TOKEN = 5
const XML_ERROR_PARTIAL_CHAR = 6
const XML_ERROR_TAG_MISMATCH = 7
const XML_ERROR_DUPLICATE_ATTRIBUTE = 8
const XML_ERROR_JUNK_AFTER_DOC_ELEMENT = 9
const XML_ERROR_PARAM_ENTITY_REF = 10
const XML_ERROR_UNDEFINED_ENTITY = 11
const XML_ERROR_RECURSIVE_ENTITY_REF = 12
const XML_ERROR_ASYNC_ENTITY = 13
const XML_ERROR_BAD_CHAR_REF = 14
const XML_ERROR_BINARY_ENTITY_REF = 15
const XML_ERROR_ATTRIBUTE_EXTERNAL_ENTITY_REF = 16
const XML_ERROR_MISPLACED_XML_PI = 17
const XML_ERROR_UNKNOWN_ENCODING = 18
const XML_ERROR_INCORRECT_ENCODING = 19
const XML_ERROR_UNCLOSED_CDATA_SECTION = 20
const XML_ERROR_EXTERNAL_ENTITY_HANDLING = 21
const XML_ERROR_NOT_STANDALONE = 22
const XML_ERROR_UNEXPECTED_STATE = 23
const XML_ERROR_ENTITY_DECLARED_IN_PE = 24
const XML_ERROR_FEATURE_REQUIRES_XML_DTD = 25
const XML_ERROR_CANT_CHANGE_FEATURE_ONCE_PARSING = 26
const XML_ERROR_UNBOUND_PREFIX = 27
const XML_ERROR_UNDECLARING_PREFIX = 28
const XML_ERROR_INCOMPLETE_PE = 29
const XML_ERROR_XML_DECL = 30
const XML_ERROR_TEXT_DECL = 31
const XML_ERROR_PUBLICID = 32
const XML_ERROR_SUSPENDED = 33
const XML_ERROR_NOT_SUSPENDED = 34
const XML_ERROR_ABORTED = 35
const XML_ERROR_FINISHED = 36
const XML_ERROR_SUSPEND_PE = 37
const XML_ERROR_RESERVED_PREFIX_XML = 38
const XML_ERROR_RESERVED_PREFIX_XMLNS = 39
const XML_ERROR_RESERVED_NAMESPACE_URI = 40
# end
# enum XML_Content_Type
const XML_CTYPE_EMPTY = 1
const XML_CTYPE_ANY = 2
const XML_CTYPE_MIXED = 3
const XML_CTYPE_NAME = 4
const XML_CTYPE_CHOICE = 5
const XML_CTYPE_SEQ = 6
# end
# enum XML_Content_Quant
const XML_CQUANT_NONE = 0
const XML_CQUANT_OPT = 1
const XML_CQUANT_REP = 2
const XML_CQUANT_PLUS = 3
# end
@ctypedef XML_Content Void
@ctypedef XML_ElementDeclHandler Ptr{Void}
@ctypedef XML_AttlistDeclHandler Ptr{Void}
@ctypedef XML_XmlDeclHandler Ptr{Void}
type XML_Memory_Handling_Suite
  malloc_fcn::Ptr{Void}
  realloc_fcn::Ptr{Void}
  free_fcn::Ptr{Void}
end
@ctypedef XML_StartElementHandler Ptr{Void}
@ctypedef XML_EndElementHandler Ptr{Void}
@ctypedef XML_CharacterDataHandler Ptr{Void}
@ctypedef XML_ProcessingInstructionHandler Ptr{Void}
@ctypedef XML_CommentHandler Ptr{Void}
@ctypedef XML_StartCdataSectionHandler Ptr{Void}
@ctypedef XML_EndCdataSectionHandler Ptr{Void}
@ctypedef XML_DefaultHandler Ptr{Void}
@ctypedef XML_StartDoctypeDeclHandler Ptr{Void}
@ctypedef XML_EndDoctypeDeclHandler Ptr{Void}
@ctypedef XML_EntityDeclHandler Ptr{Void}
@ctypedef XML_UnparsedEntityDeclHandler Ptr{Void}
@ctypedef XML_NotationDeclHandler Ptr{Void}
@ctypedef XML_StartNamespaceDeclHandler Ptr{Void}
@ctypedef XML_EndNamespaceDeclHandler Ptr{Void}
@ctypedef XML_NotStandaloneHandler Ptr{Void}
@ctypedef XML_ExternalEntityRefHandler Ptr{Void}
@ctypedef XML_SkippedEntityHandler Ptr{Void}
type XML_Encoding
  map::Void
  data::Ptr{Void}
  convert::Ptr{Void}
  release::Ptr{Void}
end
@ctypedef XML_UnknownEncodingHandler Ptr{Void}
# enum XML_Parsing
const XML_INITIALIZED = 0
const XML_PARSING = 1
const XML_FINISHED = 2
const XML_SUSPENDED = 3
# end
type XML_ParsingStatus
  parsing::Int32
  finalBuffer::XML_Bool
end
# enum XML_ParamEntityParsing
const XML_PARAM_ENTITY_PARSING_NEVER = 0
const XML_PARAM_ENTITY_PARSING_UNLESS_STANDALONE = 1
const XML_PARAM_ENTITY_PARSING_ALWAYS = 2
# end
type XML_Expat_Version
  major::Int32
  minor::Int32
  micro::Int32
end
# enum XML_FeatureEnum
const XML_FEATURE_END = 0
const XML_FEATURE_UNICODE = 1
const XML_FEATURE_UNICODE_WCHAR_T = 2
const XML_FEATURE_DTD = 3
const XML_FEATURE_CONTEXT_BYTES = 4
const XML_FEATURE_MIN_SIZE = 5
const XML_FEATURE_SIZEOF_XML_CHAR = 6
const XML_FEATURE_SIZEOF_XML_LCHAR = 7
const XML_FEATURE_NS = 8
const XML_FEATURE_LARGE_SIZE = 9
const XML_FEATURE_ATTR_INFO = 10
# end
type XML_Feature
  feature::Int32
  name::Ptr{XML_LChar}
  value::Int32
end


text:
"""

    plot_density_check(y, yrep, plot_legend, args...; kwargs...)

Draw a density plot of a random sample of draws from the response variable posterior distribution against the density estimation of the actual data to visualise model fit.

Usage:
```julia-repl
plot_density_check(y, yrep, plot_legend)
```

Details:

The `yrep` matrix should have individual draws from the posterior distribution in rows and unique values (matching the length of the `y` vector in length) in columns.

Arguments:

- `y` : The vector of response variable values.
- `yrep` : The Draws x Values matrix of posterior predictions.
- `plot_legend` : Boolean of whether to add a legend to the plot or not.
"""
function plot_density_check(y::Array, yrep, plot_legend::Bool, args...; kwargs...)

    # Check object sizes

    size(y, 2) == 1 || error("`y` should have a dimension of N x 1.")
    length1 = size(y, 1)

    if isa(yrep, Array)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    elseif isa(yrep, DataFrame)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    else
        error("`yrep` should be an object of class `Array` or `DataFrame` containing Draws x Values matrix information.")
    end

    length1 == length2 || error("Number of columns in `yrep` should match the length of vector `y`.")

    #-------- Draw plot --------

    gr() # gr backend for graphics
    drawscolour = cgrad(:blues)[1]
    actualcolour = cgrad(:blues)[2]
    Random.seed!(123) # Fix seed for reproducibility
    myPlot = plot()

    # Wrangle from wide to long

    if isa(yrep, Array)
        tmp = DataFrame(yrep)
    else
        tmp = yrep
    end

    # Plot posterior draws

    for n in 1:nrows

        tmp2 = DataFrame(tmp[n, :])
        tmp2 = stack(tmp2, 1:length2)

        # Add iteration to the plot

        if n == nrows
            plot!(tmp2[!, :value], linealpha = 0.4, xlabel = "", ylabel = "", label = "yrep",
            seriestype = :density, color = drawscolour, size = (400, 400))
        else
            plot!(tmp2[!, :value], linealpha = 0.4, xlabel = "", ylabel = "", label = "",
            seriestype = :density, color = drawscolour, size = (400, 400))
        end
    end

    # Plot actual data

    plot!(y, linealpha = 1, xlabel = "Value", ylabel = "Density", label = "y",
        seriestype = :density, color = actualcolour, legend = plot_legend,
        title = "Posterior Predictive Check", linewidth = 2)

    return myPlot
end

"""

    plot_hist_check(y, yrep, plot_legend, args...; kwargs...)

Draw a plot with a binned histogram of posterior-predicted response variable values against a measure of centrality of the actual data to visualise model fit.

Usage:
```julia-repl
plot_hist_check(y, yrep, plot_legend)
```

Arguments:

- `y` : The vector of response variable values.
- `yrep` : The Draws x Values matrix of posterior predictions.
- `plot_legend` : Boolean of whether to add a legend to the plot or not.
"""
function plot_hist_check(y::Array, yrep, plot_legend::Bool, args...; kwargs...)

    # Check object sizes

    size(y, 2) == 1 || error("`y` should have a dimension of N x 1.")
    length1 = size(y, 1)

    if isa(yrep, Array)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    elseif isa(yrep, DataFrame)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    else
        error("`yrep` should be an object of class `Array` or `DataFrame` containing Draws x Values matrix information.")
    end

    length1 == length2 || error("Number of columns in `yrep` should match the length of vector `y`.")

    # Compute median for real data vector

    m = median(y)

    # Wrangle posterior draws from wide to long and compute median for each value

    tmp = DataFrame(yrep)
    tmp = stack(tmp, 1:length2)
    tmp = combine(groupby(tmp, :variable), :value => median)

    # Set up graphics helpers

    drawscolour = cgrad(:blues)[1]
    actualcolour = cgrad(:blues)[2]

    # Draw plot

    myPlot = plot(tmp[!, :value_median], seriestype = :histogram, fillalpha = 0.5,
                  xlabel = "Value", ylabel = "", label = "yrep",
                  color = drawscolour, title = "Posterior Predictive Check", size = (400, 400),
                  legend = plot_legend)

    plot!([m], seriestype = "vline", color = actualcolour, label = "y", linewidth = 2.5)

    return myPlot
end

"""

    plot_ecdf_check(y, yrep, plot_legend, args...; kwargs...)

Draw an empirical cumulative density function plot of a random sample of draws from the response variable posterior distribution against the density estimation of the actual data to visualise model fit.

Usage:
```julia-repl
plot_ecdf_check(y, yrep, plot_legend)
```

Details:

The `yrep` matrix should have individual draws from the posterior distribution in rows and unique values (matching the length of the `y` vector in length) in columns.

Arguments:

- `y` : The vector of response variable values.
- `yrep` : The Draws x Values matrix of posterior predictions.
- `plot_legend` : Boolean of whether to add a legend to the plot or not.
"""
function plot_ecdf_check(y::Array, yrep, plot_legend::Bool, args...; kwargs...)

    # Check object sizes

    size(y, 2) == 1 || error("`y` should have a dimension of N x 1.")
    length1 = size(y, 1)

    if isa(yrep, Array)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    elseif isa(yrep, DataFrame)
        length2 = size(yrep, 2)
        nrows = size(yrep, 1)
    else
        error("`yrep` should be an object of class `Array` or `DataFrame` containing Draws x Values matrix information.")
    end

    length1 == length2 || error("Number of columns in `yrep` should match the length of vector `y`.")

    #-------- Draw plot --------

    gr() # gr backend for graphics
    drawscolour = cgrad(:blues)[1]
    actualcolour = cgrad(:blues)[2]
    Random.seed!(123) # Fix seed for reproducibility
    myPlot = plot()

    # Wrangle from wide to long

    if isa(yrep, Array)
        tmp = DataFrame(yrep)
    else
        tmp = yrep
    end

    # Compute ECDF and plot posterior draws

    for n in 1:nrows

        tmp2 = DataFrame(tmp[n, :])
        tmp2 = stack(tmp2, 1:length2)

        # ECDF calculation

        myECDF = StatsBase.ecdf(tmp2[!, :value])

        # Add iteration to the plot

        if n == nrows
            plot!(myECDF, linealpha = 0.4, xlabel = "", ylabel = "", label = "yrep",
            color = drawscolour, size = (400, 400))
        else
            plot!(myECDF, linealpha = 0.4, xlabel = "", ylabel = "", label = "",
            color = drawscolour, size = (400, 400))
        end
    end

    # Compute ECDF of actual data and add to plot

    yECDF = StatsBase.ecdf(y)

    if plot_legend == true
        plot!(yECDF, linealpha = 1, xlabel = "Value", ylabel = "Density", label = "y",
        color = actualcolour, legend = :bottomright, title = "Posterior Predictive Check",
        linewidth = 1.5)
    else
        plot!(yECDF, linealpha = 1, xlabel = "Value", ylabel = "Density", label = "y",
        color = actualcolour, legend = false, title = "Posterior Predictive Check", linewidth = 1.5)
    end

    return myPlot
end