text:
using PlotlyJS, RDatasets, DataFrames

function violin_box_overlay()
    y = abs.(100 .* randn(300))
    data = [
        violin(x0="sample 1", name="violin", y=y, points="all"),
        box(x0="sample 1", name="box", y=y, boxpoints=false)
    ]
    plot(data, Layout(legend=attr(x=0.95, xanchor="right")))
end

function violin_grouped()
    days = repeat(["day 1", "day 2"], inner=5)
    y_kale = [0.2, 0.2, 0.6, 1, 0.5, 0.4, 0.2, 0.7, 0.9, 0.1, 0.5, 0.3]
    y_radish = [0.6, 0.7, 0.3, 0.6, 0, 0.5, 0.7, 0.9, 0.5, 0.8, 0.7, 0.2]
    y_carrot = [0.1, 0.3, 0.1, 0.9, 0.6, 0.6, 0.9, 1, 0.3, 0.6, 0.8, 0.5]
    colors = ["#3D9970", "#FF4136", "#FF851B"]
    names = ["kale", "radishes", "carrots"]
    ys = (y_kale, y_radish, y_carrot)

    data = [
        violin(
            y=y, name=name, x=days, span=[0, nothing], jitter=0, points="all",
            marker=attr(symbol="line-ew", color=color, line=attr(color=color, width=2))
        ) for (y, name, color) in zip(ys, names, colors)
    ]
    layout = Layout(
        yaxis=attr(zeroline=false, title="normalized moisture"),
        violinmode="group"
    )
    plot(data, layout)
end

function violin_nonlinear()
    p1 = plot(
        violin(
            x=["1798-01-01", "1798-04-04", "1798-05-05",
               "1798-05-05", "1798-07-05", "1798-07-22", "1799-01-01"],
            orientation="h", box_visible=true, xcalendar="discworld",
            name="discworld dates"
       )
    )
    p2 = plot(
        violin(
            x=["A", "B", "C", "C", "C", "D", "G"],
            orientation="h", box_visible=true, xcalendar="discworld",
            name="categories"
       ), Layout(xaxis_categoryarray='A':'G')
    )
    p = [p1; p2]
    relayout!(p, showlegend=false, margin_l=100)
    p
end

function violin_old_faithful()
    y = [79, 54, 74, 62, 85, 55, 88, 85, 51, 85, 54, 84, 78, 47, 83, 52, 62,
    84, 52, 79, 51, 47, 78, 69, 74, 83, 55, 76, 78, 79, 73, 77, 66, 80, 74, 52,
    48, 80, 59, 90, 80, 58, 84, 58, 73, 83, 64, 53, 82, 59, 75, 90, 54, 80, 54,
    83, 71, 64, 77, 81, 59, 84, 48, 82, 60, 92, 78, 78, 65, 73, 82, 56, 79, 71,
    62, 76, 60, 78, 76, 83, 75, 82, 70, 65, 73, 88, 76, 80, 48, 86, 60, 90, 50,
    78, 63, 72, 84, 75, 51, 82, 62, 88, 49, 83, 81, 47, 84, 52, 86, 81, 75, 59,
    89, 79, 59, 81, 50, 85, 59, 87, 53, 69, 77, 56, 88, 81, 45, 82, 55, 90, 45,
    83, 56, 89, 46, 82, 51, 86, 53, 79, 81, 60, 82, 77, 76, 59, 80, 49, 96, 53,
    77, 77, 65, 81, 71, 70, 81, 93, 53, 89, 45, 86, 58, 78, 66, 76, 63, 88, 52,
    93, 49, 57, 77, 68, 81, 81, 73, 50, 85, 74, 55, 77, 83, 83, 51, 78, 84, 46,
    83, 55, 81, 57, 76, 84, 77, 81, 87, 77, 51, 78, 60, 82, 91, 53, 78, 46, 77,
    84, 49, 83, 71, 80, 49, 75, 64, 76, 53, 94, 55, 76, 50, 82, 54, 75, 78, 79,
    78, 78, 70, 79, 70, 54, 86, 50, 90, 54, 54, 77, 79, 64, 75, 47, 86, 63, 85,
    82, 57, 82, 67, 74, 54, 83, 73, 73, 88, 80, 71, 83, 56, 79, 78, 84, 58, 83,
    43, 60, 75, 81, 46, 90, 46, 74]
    plot(violin(y=y, points="all", name="Old Faithful", meanline_visible=true))
end

function violin_side_by_side()
    # requires RDatasets and DataFrames
    tips = RDatasets.dataset("reshape2", "tips")
    parts = zip(
        ("Female", "Male"),
        ("positive", "negative"),
        ("#bebada", "#8dd3c7"),
        (1.0, -0.6)
    )
    traces = GenericTrace[]
    for (sex, side, color, pointpos) in parts
        sub_tips = tips[tips[!, :Sex] .== sex, :]
        sub_traces = violin(sub_tips,
            group=:Day,
            x=:TotalBill, y0=(df) -> df[1, :Day],
            side=side, orientation="h",
            marker=attr(line=attr(width=2, color=color), symbol="line-ns"),
            line_color=color,
            hoveron="points+kde", text=(df) -> "Sample length $(size(df, 1))",
            scalemode="count", scalegroup=sex, legendgroup=sex, name=sex,
            points="all", jitter=0, pointpos=pointpos,
            span=[0],
            box_visible=true, meanline_visible=true,
            showlegend=false,
        )
        sub_traces[1][:showlegend] = true
        append!(traces, sub_traces)
    end
    # TODO: make the layout
    layout = Layout(
        hovermode="closest", violinmode="overlay",
        title="Total bill distribution<br><i>scaled by number of bills per gender",
        legend_tracegroupgap=0, violingap=0, violingroupgap=0,
        yaxis=attr(showgrid=true, categoryarray=["Thur", "Fri", "Sat", "Sun"]),
    )
    plot(traces, layout)
end

function violin_style()
    y1 = vcat(abs.(20 .* rand(100)), rand(UInt16, 300) .* 500 ./ typemax(UInt16))
    y2 = [25.261999999999997, 66.5419, 98.2114, 0.09070629 ]
    box = attr(fillcolor="black", line_color="black", width=0.01)
    span = [0, nothing]
    trace1 = violin(
        bandwidth=5, points=false, y=y1, name="Radial Velocity",
        span=span, line_color="#67353E", box=box
    )
    trace2 = violin(
        bandwidth=5, points=false, y=y2, name="Pulsar Timing",
        span=span, line_color="#34ABA2", box=box
    )

    layout = Layout(
        paper_bgcolor="#d3d3d3", plot_bgcolor="#d3d3d3",
        showlegend=false, violingap=0, xaxis_side="top",
        yaxis=attr(
            showline=false, showticklabels=false, range=(-5, 550),
            zeroline=false, visible=false, showgrid=false,
        ),
        annotations=[attr(
            text="Orbital Period Estimations", font_size=20,
            xref="paper", yref="paper", showarrow=false,
        )]
    )
    plot([trace1, trace2], layout)
end


text:
""" Data Migration functors
"""
module DataMigration
export Functor, Delta, Sigma, migrate!

using ...Present, ...Theories
using ...Theories: SchemaDesc, ob, hom, dom, codom, attr
using ..FinSets, ..CSets, ..Limits, ...Graphs, ..FreeDiagrams

import ...CategoricalAlgebra.FreeDiagrams: FreeDiagram
import ...Present: Presentation
import ...Theories: dom, codom, Ob, Hom

""" Abstract type for a functor.
"""
abstract type AbstractFunctor end

"""   Functor

A functor ``F: \\mathcal{C} \\to \\mathcal{D}`` consists of a map of objects and a
map of homomorphisms from a domain category ``\\mathcal{C}``  to a codomain
category ``\\mathcal{D}``.
"""
struct Functor{O,H,C,D} <: AbstractFunctor
  FOb::O
  FHom::H
  dom::C
  codom::D
end

# interface for abstract functor
dom(F::Functor) = F.dom
codom(F::Functor) = F.codom
Ob(F::Functor) = F.FOb
Hom(F::Functor) = F.FHom

# To construct a functor where FOb and FHom are keyed by the symbols of
# the domain and codomain presentations. Turn them into a dictionaries
# whose keys and values are homs.
function Functor(FOb::Dict{Symbol, Symbol}, FHom::Dict{Symbol, H},
                 dom::Presentation, codom::Presentation) where H
  Functor(
    Dict(dom[c_name] => codom[Fc_name] for (c_name, Fc_name) in FOb),
    Dict(dom[f_name] => get_hom(codom, Ff_name) for (f_name, Ff_name) in FHom),
    dom, codom
  )
end

get_hom(presY::Presentation, Ff::FreeSchema.Hom) = Ff
get_hom(presY::Presentation, Ff::Symbol) = presY[Ff]
get_hom(presY::Presentation, Ff::Array) =
  reduce(compose, map(h -> get_hom(presY, h), Ff))

""" Abstract type for migration functor which migrates data from one ACSet to another
"""
abstract type MigrationFunctor <: AbstractFunctor end

(F::MigrationFunctor)(X::ACSet) = F(codom(F)(),X)

dom(MF::MigrationFunctor) = MF.dom
codom(MF::MigrationFunctor) = MF.codom

#### Delta Migration
###################

""" Pullback functorial data migration from one ACSet to another.

Note that this operation is contravariant: the data is transferred from `Y` to
`X` but the functor, represented by two dictionaries, maps the schema for `X`
to the schema for `Y`.

When the functor is the identity, this function is equivalent to
[`copy_parts!`](@ref).
"""
struct DeltaMigration{D,C} <: MigrationFunctor
  F::Functor # on the schemas
  dom::D
  codom::C
end

"""   Delta(F::Functor)

Given a functor ``F: \\mathcal{C} \\to \\mathcal{D}`` produces a `MigrationFunctor`
which maps a ``\\mathcal{D}``-set ``X`` to the ``\\mathcal{C}``-set
``\\Delta_F(X) := X \\circ F``.

See (Spivak, 2014, *Category Theory for the Sciences*) for details.
"""
Delta(F::Functor, dom, codom) = DeltaMigration(F, dom, codom)

function (ΔF::DeltaMigration)(X::StructACSet{S}, Y::ACSet) where {S}
  FOb = Ob(ΔF.F)
  FHom = Hom(ΔF.F)

  partsX = Dict(map(ob(S)) do c_name
    c = dom(ΔF.F)[c_name]
    c => add_parts!(X, c_name, nparts(Y, nameof(FOb[c])))
  end)

  for (f, Ff) in collect(FHom)
    doms = partsX[dom(f)]
    subpt = f isa FreeSchema.Hom ? partsX[codom(f)][subpart(Y, Ff)] : subpart(Y, Ff)
    set_subpart!(X, doms, nameof(f), subpt)
  end

  return X
end

"""   migrate!(X::ACSet, Y::ACSet, FOb, FHom)

Migrates the data from `Y` to `X` via the pullback
data migration induced by the functor defined on objects by `FOb` and
on morphisms by `FHom`.
"""
migrate!(X::ACSet, Y::ACSet, F::Functor) = Delta(F, typeof(Y), typeof(X))(X,Y)

migrate!(X::ACSet, Y::ACSet, FOb, FHom) =
  migrate!(X,Y, Functor(FOb, FHom, Presentation(X), Presentation(Y)))

function (::Type{T})(Y::ACSet, FOb::AbstractDict, FHom::AbstractDict) where T <: ACSet
  X = T()
  migrate!(X, Y, FOb, FHom)
end

function (::Type{T})(Y::ACSet, F::Functor) where T <: ACSet
  X = T()
  migrate!(X, Y, F)
end

### SigmaMigration
###################
""" Left pushforward functorial data migration from one ACSet to another.
"""
struct SigmaMigration{D, C, CC} <: MigrationFunctor
  F::Functor # on the schemas
  dom::D
  codom::C
  comma_cats::CC

  function SigmaMigration(F::Functor, dom::D, codom::C) where {D, C}
    cc = get_comma_cats(F)
    new{D, C, typeof(cc)}(F, dom, codom, cc)
  end
end

"""   Sigma(F::Functor)

Given a functor ``F: \\mathcal{C} \\to \\mathcal{D}`` produces a `MigrationFunctor`
which maps a ``\\mathcal{C}``-set ``X`` to the ``\\mathcal{D}``-set ``\\Sigma_F(X)``.
``\\Sigma_F`` is left-adjoint to ``\\Delta_F``.

For ``d \\in \\mathsf{ob}\\mathcal{D}``,
``\\Sigma_F(X)(d) := \\mathsf{colim}_{F \\downarrow d} X \\circ \\pi`` where
``\\pi: F\\downarrow d \\to \\mathcal{C}`` is the projection.

See (Spivak, 2014, *Category Theory for the Sciences*) for details.
"""
Sigma(F::Functor, dom, codom) = SigmaMigration(F, dom, codom)

function (ΣF::SigmaMigration)(Y::ACSet, X::ACSet)
  comma_cats = ΣF.comma_cats
  diagramD = FreeDiagram(codom(ΣF.F))

  # define Y on objects by taking colimits
  colimX = map(parts(diagramD, :V)) do i
    F∇d = ob(comma_cats, i)
    Xobs = map(ob(F∇d)) do (c,_)
      FinSet{Int, Int}(nparts(X, nameof(c)))
    end
    Xhoms = map(parts(F∇d, :E)) do g
      FinFunction(X[nameof(hom(F∇d, g))], Xobs[src(F∇d, g)], Xobs[tgt(F∇d, g)])
    end
    colimit(FreeDiagram(Xobs, collect(zip(Xhoms, src(F∇d), tgt(F∇d)))))
  end

  for d in parts(diagramD, :V)
    add_parts!(Y, nameof(ob(diagramD, d)), length(apex(colimX[d])))
  end

  # Define Y on morphisms by the universal property
  for g in parts(diagramD, :E)
    if nparts(Y, nameof(dom(hom(diagramD, g)))) == 0
      continue
    end
    set_subpart!(Y, nameof(hom(diagramD, g)),
      universal(colimX[src(diagramD, g)],
        Multicospan(legs(colimX[tgt(diagramD, g)])[hom(comma_cats, g)[:V].func])
      ).func
    )
  end
  return Y
end

"""   inv(D::AbstractDict)

Given a dictionary `D` returns a function that maps a proposed value
`y` to the array of all keys with value `y`.
"""
inv(D::AbstractDict) = y -> [k for (k,v) in D if v == y]

""" add_hom!(d::FreeDiagram{Ob, Hom}, src_ob::Ob, tgt_ob::Ob, hom::Hom)

Adds a hom to `d` from the vertex with object `src_ob` to the vertex with object `tgt_ob`.
"""
function add_hom!(d::FreeDiagram, src_ob, tgt_ob, hom)
  src_idx = first(incident(d, src_ob, :ob))
  tgt_idx = first(incident(d, tgt_ob, :ob))
  return add_edge!(d, src_idx, tgt_idx, hom = hom)
end

"""   comma_cats(diagramD::FreeDiagram{FreeSchema.Ob, FreeSchema.Hom}, FOb, FHom)

Given a free category ``\\mathcal{D}`` with no cycles and
a functor represented by the pair `(FOb, FHom)`, returns a diagram
``\\mathcal{D} \\to \\mathsf{Cat}`` which on objects takes ``d \\in \\mathcal{D}`` to the
comma category ``F \\downarrow d``.
"""
function get_comma_cats(F::Functor)
  diagramD = FreeDiagram(codom(F))
  FObInv = inv(Ob(F)); FHomInv = inv(Hom(F))
  comma_cats = FreeDiagram{FreeDiagram, ACSetTransformation}()
  add_vertices!(comma_cats,
    nparts(diagramD, :V),
    ob = map(parts(diagramD, :V)) do _
      FreeDiagram{Pair{FreeSchema.Ob, FreeSchema.Hom}, FreeSchema.Hom}()
    end
  )

  for d in topological_sort(diagramD)
    F∇d = ob(comma_cats, d)
    id_d = id(ob(diagramD, d))

    # If FOb[c] = d, add an objects (c, id(d)) to F∇d
    preimage_d = FObInv(ob(diagramD, d))
    id_obs = add_vertices!(F∇d, length(preimage_d), ob = map(c->c=>id_d, preimage_d))

    # if FHom[h] = id(d), add a hom h: (dom(h), id(d)) -> (codom(h), id(d)) to F∇d
    id_homs = map(FHomInv(id_d)) do h
      add_hom!(F∇d, dom(h) => id_d, codom(h) => id_d, h)
    end

    # for g: d' -> d in D (note that F∇d' must have already been constructed)
    #     populate F∇d with obs and homs coming from F∇d′
    for g in incident(diagramD, d, :tgt)
      d′ = src(diagramD, g)
      F∇g = comma_cat_hom!(F∇d, ob(comma_cats, d′), id_d, hom(diagramD, g), FHomInv)
      add_edge!(comma_cats, d′, d, hom=F∇g)
    end
  end

  return comma_cats
end

function comma_cat_hom!(F∇d, F∇d′, id_d, g, FHomInv)
  # for (c,f) in F∇d' add ob (c, compose(f,g)) to F∇d
  vs = add_vertices!(F∇d, nparts(F∇d′, :V), ob = map(ob(F∇d′)) do (c,f)
    c => compose(f, g)
  end)

  # for h: (c, f) -> (c', f') in diagram_d', add hom h to F∇d
  es = add_edges!(F∇d, vs[src(F∇d′)], vs[tgt(F∇d′)], hom = hom(F∇d′))

  # for every (c,f) in the recently added objects. If FHom[h] == f, then add the hom
  #       h: (c,f) -> (codom(h), idv)
  # relies on D being free
  for (c, f) in ob(F∇d, vs)
    for h in FHomInv(f)
      dom(h) == c && add_hom!(F∇d, c => f, codom(h) => id_d, h)
    end
  end

  # return the inclusion from F∇d′ to F∇d
  return ACSetTransformation((V = collect(vs), E = collect(es)), F∇d′, F∇d)
end

### Translate between ACSets, Presentations, and FreeDiagrams
###############################################################
"""Get the Schema from an ACSet
"""
function Presentation(::StructACSet{S}) where S
  return Presentation(S)
end

"""   FreeDiagram(pres::Presentation{FreeSchema, Symbol})

Returns a `FreeDiagram` whose objects are the generating objects of `pres` and
whose homs are the generating homs of `pres`.
"""
function FreeDiagram(pres::Presentation{Schema, Symbol}) where Schema
  obs = Array{FreeSchema.Ob}(generators(pres, :Ob))
  homs = Array{FreeSchema.Hom}(generators(pres, :Hom))
  doms = map(h -> generator_index(pres, nameof(dom(h))), homs)
  codoms = map(h -> generator_index(pres, nameof(codom(h))), homs)
  return FreeDiagram(obs, collect(zip(homs, doms, codoms)))
end

end #module

ext: jl


text:
@testset "dataSelection.jl" begin
    data_surv = CSO.dataset("cso_raw")
    data_trudeau = CSO.dataset("Montreal-Trudeau_15min_precipitations")
    data = innerjoin(data_surv, data_trudeau, on=:Date)

    dropmissing!(data, :Duration)
    addColOverflow!(data, "Duration", "Surverse", 3)

    @testset "pctOnly(df, ColPct)" begin
        newDF = pctOnly(data, "d24h")

        #Only rows where precipitations are positive are kept
        @test sum(newDF.d24h) == sum(data.d24h)

        #Ancient DF has not been changed
        @test nrow(data) != nrow(newDF)
    end

    @testset "pctOnly!(df, ColPct)" begin
        n = sum(data.d24h)
        transformedDF = pctOnly!(data, "d24h")

        #Only rows where precipitations are positive are kept
        @test sum(transformedDF.d24h) == n

        #Ancient DF has been changed
        @test nrow(data) == nrow(transformedDF)
    end

    @testset "overflowOnly(df, ColOverflow)" begin
        data = data
        n = sum(data.Surverse)
        newDF = overflowOnly(data, "Surverse")

        #Only rows where precipitations are positive are kept
        @test sum(newDF.Surverse) == n

        #Ancient DF has not been changed
        @test nrow(data) != nrow(newDF)
    end

    @testset "overflowOnly!(df, ColOverflow)" begin
        data = data
        n = sum(data.Surverse)
        transformedDF = overflowOnly!(data, "Surverse")

        #Only rows where precipitations are positive are kept
        @test sum(transformedDF.Surverse) == n

        #Ancient DF has been changed
        @test nrow(data) == nrow(transformedDF)
    end

end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
"""
    NCBITaxon

An `NCBITaxon` is minimally defined by its name and its id. Note that across the
entire package the `name` field is the valid scientific name according to NCBI.
These objects can be passed to `authority`, `synonyms`, and `vernacular` to get
other names as strings only.
"""
struct NCBITaxon
    name::String
    id::Int
end

Base.show(io::IO, t::NCBITaxon) = print(io, "$(t.name) (ncbi:$(t.id))")

@enum NCBINameClass begin
    class_scientific_name=1
    class_equivalent_name=2
    class_common_name=3
    class_authority=4
    class_synonym=5
    class_in_part=6
    class_includes=7
    class_acronym=8
    class_genbank_common_name=9
    class_genbank_synonym=10
    class_genbank_acronym=11
    class_blast_name=12
end


text:
module Visualise

using PyCall, JSON
using JuLIP: set_pbc!, positions, neighbourlist
using JuLIP.ASE: AbstractAtoms, ASEAtoms, write, rnn, chemical_symbols

export draw

@pyimport IPython.display as ipydisp
@pyimport imolecule

"""

### KW-args (copied from imolecule)

* `size`: Starting dimensions of visualization, in pixels, e.g., `(500,500)`
* `drawing_type`: Specifies the molecular representation. Can be 'ball and
    stick', "wireframe", or 'space filling'.
* `camera_type`: Can be 'perspective' or 'orthographic'.
* `shader`: Specifies shading algorithm to use. Can be 'toon', 'basic',
    'phong', or 'lambert'.
* `display_html`: If True (default), embed the html in a IPython display.
    If False, return the html as a string.
* `element_properties`: A dictionary providing color and radius information
    for custom elements or overriding the defaults in imolecule.js
* `show_save`: If True, displays a save icon for rendering molecule as an
    image.

## More kwargs:

* bonds: :babel, :auto
* cutoff: :auto or a number
"""
function draw(at::AbstractAtoms; bonds=:babel, box=:auto,
            camera_type="perspective", size=(500,500), display_html=false,
            cutoff=:auto,
            kwargs...)

   # TODO: implement box (ignore for now)
   # TODO: check whether we are in a notebook or REPL
   if bonds == :babel
      # in this case we assume that OpenBabel and in particular `pybel` is
      # installed; we write the atoms object as xyz and let `pybel`
      # do the rest
      fn = "$(tempname()).xyz"
      write(fn, at)
      out = ipydisp.HTML(
           imolecule.draw( fn, format="xyz", camera_type=camera_type,
                                size=size, display_html=false, kwargs...) )
      rm(fn)
      return out
   elseif bonds == :auto
      # TODO: bonds = autobondlenghts(at)
   elseif !isa(bonds, Dict)
      error("""draw(at ...) : at the momement, `bonds` must be
               `:babel`, `:auto` or a `Dict` containing the bonds-length
               information""")
   end

   # we have a `Dict` that contains all the bond-length information
   # (which we guessed really; this needs more work)
   # we can therefore write a JSON file with atom positions + connectivity
   # and let `imolecule` load that; this never needs to import `pybel` so it
   # will work even if OpenBabel is not installed.

   # TODO: bond-length guesses for multiple species
   # TODO: bounding box especially with PBCs
   # TODO: check which other information imolecule accepts

   # temporarily assume there is a single species!
   # and make the bond length 120% of the NN-length
   if cutoff == :auto
      bondlength = 1.2 * rnn( chemical_symbols(at)[1] )
   else
      bondlength = cutoff
   end

   # we need to turn of PBC otherwise we get weird bonds
   set_pbc!(at, (false, false, false))
   atoms = [Dict("element" => e, "location" => Vector(x))
            for (e, x) in zip(chemical_symbols(at), positions(at)) ]
   nlist = neighbourlist(at, bondlength)
   # we are assuming that nlist is the Matscipy nlist
   ij = unique(sort([nlist.i nlist.j], 2), 1) - 1
   bondlist = [Dict("atoms" => [ij[n,1];ij[n,2]], "order" => 1)
               for n = 1:Base.size(ij, 1)]
   molecule_data = Dict("atoms" => atoms, "bonds" => bondlist)

   # write to a temporary file
   fn = "$(tempname()).json"
   fio = open(fn, "w")
   print(fio, JSON.json(molecule_data))
   close(fio)

   # plot
   out = ipydisp.HTML(
        imolecule.draw(fn, format="json", camera_type=camera_type,
                             size=size, display_html=false, kwargs...) )
   rm(fn)
   return out
end

end # module Visualise


text:
input = joinpath(@__DIR__, "input")
lines = readlines(input)

numbers = map(n -> parse(Int, n), split(lines[1], ','))
boards = lines[2:end]

function genboards(boards)
    filter!(l -> l != "", boards)
    rows = map(l -> map(r -> parse(Int, r), split(l)), boards)
    acc = []
    while length(rows) != 0
        append!(acc, [rows[1:5]])
        rows = rows[6:end]
    end
    return acc
end

function score(numbers, board)
    for i in 1:length(numbers)
        drawn = numbers[1:i]
        b = copy(board)
        # rows (b) and columns (map(collect, zip(b...)))
        append!(b, map(collect, zip(b...)))
        for s in b
            if issubset(s, drawn)
                rem = setdiff(Iterators.flatten(board), drawn)
                # return "round" board won in and score
                return (length(drawn), sum(rem) * last(drawn))
            end
        end
    end
end

boards = genboards(boards)
scores = map(b -> score(numbers, b), boards)

(_, fastest) = minimum(scores)
(_, slowest) = maximum(scores)

p1 = fastest
p2 = slowest

println("-----------------------------------------------------------------------")
println("giant squid -- part one :: $p1")
println("giant squid -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 44736)
@assert(p2 == 1827)


text:
module TasksController
using Genie.Renderer
using Genie.Router
using SearchLight
using Tasks

function index()
  html!(:tasks, :tasks, tasks = SearchLight.all(Tasks.Task))
end

function new()
  html!(:tasks, :new)
end

function create()
  Tasks.Task(content = @params(:task_content), done = false) |> save && redirect_to(:get_tasks)
end

function edit()
  # html!(:tasks, :edit, task = SearchLight.find_one(Tasks.Task, @params:(:task_id)))
end

function update()
  Tasks.Task(id = @params(:task_id), done = @params:(:task_done)) |> save && redirect_to(:get_tasks)
end

end

max_stars_repo_name: yshutaro/GenieTodoList

max_issues_repo_name: yshutaro/GenieTodoList

max_forks_repo_name: yshutaro/GenieTodoList


text:
module Demo

include("function_1.jl")

export my_f, my_ff

end

max_stars_repo_name: aivodji/Yg

max_issues_repo_name: aivodji/Yg

max_forks_repo_name: aivodji/Yg


text:
# TODO: Should this also return false (or error) if `a` or `b` don't have unique elements?
issetequal(a,b) = a ⊆ b && b ⊆ a


text:
### A Pluto.jl notebook ###
# v0.17.2

using Markdown
using InteractiveUtils

# ╔═╡ c24044bc-2fc3-4e34-b99c-7dc4ec6a7515
import Pkg

# ╔═╡ 42da96ce-8122-4aa0-99d6-635f68db6e66
Pkg.activate(".")

# ╔═╡ 88e46392-1fbc-4079-8207-dd5b3ac01d7d
using DeltaRCWA

# ╔═╡ 16316725-c608-49a5-b90c-8caff973683a
using Plots

# ╔═╡ 83d9578c-3960-486c-ad3d-b4449f7bae6f
md"
# Fresnel scattering at a dielectric interface

Note that our GSTCs simplify to continuity of the tangential fields when
the surface conductivities vanish.
Therefore, our scattering matrix should recover the Fresnel coefficients
when we allow the medium on each side of the interface to differ.

The analytical solution for this problem is known and I will just cite it [here](https://homerreid.github.io/scuff-em-documentation/tests/FresnelScattering/FresnelScattering/)
``\newcommand{\bm}{\boldsymbol}``
"

# ╔═╡ 6cbc5f34-fe14-45e8-9fee-9787d76ad0fe
struct EmptySheet <: Sheet end

# ╔═╡ fb8df44c-a354-42de-9e99-df64dd4c5727
DeltaRCWA.ElectricResponseStyle(::Type{EmptySheet}) = Admittance()

# ╔═╡ f3d36c35-c09b-483f-bb41-046b6ecd0773
const Water = UniformMedium{1.7, 1.0}()

# ╔═╡ e899dbc0-03bb-4a13-9c09-bfad33e5c1d3
const stack = SheetStack(EmptySheet(), (Vacuum, Water))

# ╔═╡ 43d29086-510d-4870-94b7-2ff5ff34324b
n = 100 # number of grid points

# ╔═╡ 5377ae36-32fc-41c5-bd03-c737422b17f4
prob = DeltaRCWAProblem(stack, ((n, 1.0), ), 10.0, [i==2 ? 1 : 0 for i in 1:n], zeros(n))

# ╔═╡ 65c2e966-8cb0-4d76-96db-18f9045695a3
sol = solve(prob)

# ╔═╡ bf1213cd-2ec1-4d1e-8ebc-5be5af26621c
plot(sol; combine=true)

# ╔═╡ ba0c4060-f9d4-43b6-bc31-b48dcd81dd91
md"""
``k_z = \sqrt{\omega^2 \epsilon \mu - k_x^2 - k_y^2}``
"""

# ╔═╡ 687627f9-6938-4d7a-9395-9f1783b8a449
md"""
## Scattering matrix in 2d

We restrict ourselves to the decoupled TE/TM mode case (2d with no off-diagonal impedance matrix elements) for simplicity (we can make the 3d case work too).
Starting from the GSTCs for the TM modes (when ``k_y=0`` and the TE mode zero, ``E_y=H_x=0``)
```math
\begin{align}
-(H_y^{(2,+)} + H_y^{(2,-)} - H_y^{(1,+)} - H_y^{(1,-)})
	&= \sigma^e_{xx} (E_x^{(2,+)} + E_x^{(2,-)} + E_x^{(1,+)} + E_x^{(1,-)})/2
\\
(E_x^{(2,+)} + E_x^{(2,-)} - E_x^{(1,+)} - E_x^{(1,-)})
	&= -\sigma^m_{yy} (H_y^{(2,+)} + H_y^{(2,-)} + H_y^{(1,+)} + H_y^{(1,-)})/2
\end{align}
```
and the corresponding planewave relations (indexed by medium ``i``)
```math
\begin{align}
% -\omega^2 \epsilon^i \mu^i H_x^i &= \omega \epsilon^i k_z^i E_y^i
% \\
(k_x^2 -\omega^2 \epsilon^i \mu^i) H_y^i =  -(k_z^i)^2 H_y^i &= -\omega \epsilon^i k_z^i E_x^i
\\
-\omega^2 \epsilon^i \mu^i E_x^i &= -\omega \mu^i k_z^i H_y^i
% \\
% (k_x^2 -\omega^2 \epsilon^i \mu^i) E_y^i &= \omega \mu^i k_z^i H_x^i
\end{align}
```
which both are the same and let us eliminate the electric field using
```math
\begin{align}
E_x^i &= \frac{k_z^i}{\omega \epsilon^i} H_y^i
\end{align}
```
Taking into account the direction of propagation means
```math
\begin{align}
E_x^{(i,\pm)} &= \pm\frac{k_z^i}{\omega \epsilon^i} H_y^{(i,\pm)}
\end{align}
```
Note that this is the point where the 2d case is advantageous because the decoupled modes make the inversion in this case just the inverse of a diagonal matrix.
Hence the GSTCs are
```math
\begin{align}
&- \left( H_y^{(2,+)} + H_y^{(2,-)} - H_y^{(1,+)} - H_y^{(1,-)} \right)
\\ &\qquad = \sigma^e_{xx} \left( \frac{k_z^2}{2\omega\epsilon^2} \left( H_y^{(2,+)} - H_y^{(2,-)} \right) + \frac{k_z^1}{2\omega\epsilon^1} \left( H_y^{(1,+)} - H_y^{(1,-)} \right) \right)
\\
&\frac{2k_z^2}{\omega\epsilon^2} \left( H_y^{(2,+)} - H_y^{(2,-)} \right) - \frac{2k_z^1}{\omega\epsilon^1} \left( H_y^{(1,+)} - H_y^{(1,-)} \right)
\\ &\qquad = -\sigma^m_{yy} \left( H_y^{(2,+)} + H_y^{(2,-)} + H_y^{(1,+)} + H_y^{(1,-)} \right)
\end{align}
```
which gives the following scattering matrix
```math
\begin{align}
&\begin{pmatrix}
I + \sigma^e_{xx} \frac{k_z^1}{2\omega\epsilon^1}
&
-I - \sigma^e_{xx} \frac{k_z^2}{2\omega\epsilon^2}
\\
\sigma^m_{yy} + \frac{2k_z^1}{\omega\epsilon^1}
&
\sigma^m_{yy} + \frac{2k_z^2}{\omega\epsilon^2}
\end{pmatrix}
\begin{pmatrix}
H_y^{(1,-)}
\\
H_y^{(2,+)}
\end{pmatrix}
\\ =
&\begin{pmatrix}
-I + \sigma^e_{xx} \frac{k_z^1}{2\omega\epsilon^1}
&
I - \sigma^e_{xx} \frac{k_z^2}{2\omega\epsilon^2}
\\
-\sigma^m_{yy} + \frac{2k_z^1}{\omega\epsilon^1}
&
-\sigma^m_{yy} + \frac{2k_z^2}{\omega\epsilon^2}
\end{pmatrix}
\begin{pmatrix}
H_y^{(1,+)}
\\
H_y^{(2,-)}
\end{pmatrix}
\end{align}
```
"""

# ╔═╡ 860dc600-c86d-43a2-a88a-eb876eecbd6e
md"""
## Recovering Fresnel scattering

We can recover Fresnel scattering by setting the metasurface conductivities to zero, inverting the left hand side and reading off the scattered amplitudes.
To illustrate this, let's calculate this in this special case, which is just requiring the continuity of the tangential fields
```math
\begin{pmatrix}
I & -I
\\
\epsilon^2 k_z^1 & \epsilon^1 k_z^2
\end{pmatrix}
\begin{pmatrix}
H_y^{(1,-)}
\\
H_y^{(2,+)}
\end{pmatrix}
=
\begin{pmatrix}
-I & I
\\
\epsilon^2 k_z^1
& \epsilon^1 k_z^2
\end{pmatrix}
\begin{pmatrix}
H_y^{(1,+)}
\\
H_y^{(2,-)}
\end{pmatrix}
```
Let's invert the left hand side
```math
\begin{align}
\begin{pmatrix}
H_y^{(1,-)}
\\
H_y^{(2,+)}
\end{pmatrix}
&=
\begin{pmatrix}
(1 + \frac{\epsilon^2 k_z^1}{\epsilon^1 k_z^2})^{-1}
&
(\epsilon^1 k_z^2 + \epsilon^1 k_z^2)^{-1}
\\
-(1 + \frac{\epsilon^1 k_z^2}{\epsilon^2 k_z^1})^{-1}
&
(\epsilon^1 k_z^2 + \epsilon^2 k_z^1)^{-1}
\end{pmatrix}
\begin{pmatrix}
-I & I
\\
\epsilon^2 k_z^1 & \epsilon^1 k_z^2
\end{pmatrix}
\begin{pmatrix}
H_y^{(1,+)}
\\
H_y^{(2,-)}
\end{pmatrix}
\\ &=
(\epsilon^1 k_z^2 + \epsilon^2 k_z^1)^{-1}
\begin{pmatrix}
\epsilon^2 k_z^1 - \epsilon^1 k_z^2 & 2 \epsilon^1 k_z^2
\\
2 \epsilon^2 k_z^1 & \epsilon^1 k_z^2 - \epsilon^2 k_z^1
\end{pmatrix}
\begin{pmatrix}
H_y^{(1,+)}
\\
H_y^{(2,-)}
\end{pmatrix}
\end{align}
```
Note how the scattering depends only on the angle of incidence from the normal axis (the ratio of ``\omega`` and ``k_z``) and the impedances of the materials.
An analogous calculation can be done for the TE mode
"""

# ╔═╡ 271422af-db75-4373-be1c-104745f941bf
md"""
## Scattering matrix in 3d

The planewave form of Ampère's law we will use within dielectric medium ``i`` with permittivity ``\epsilon^i`` and permeability ``\mu^i`` is
```math
\bm{R} \bm{\Omega}^i \bm{R}^{-1} \bm{E}_\parallel^i = \omega \bm{R} \mu^i k_z^i \bm{H}_\parallel^i
```
where
```math
\bm{\Omega}^i = \bm{K} - \omega^2 \epsilon^i \mu^i \bm{I}
, \qquad
\bm{K} =
\begin{pmatrix}
k_x k_x & k_x k_y
\\
k_y k_x & k_y k_y
\end{pmatrix}
, \qquad
\bm{R} =
\begin{pmatrix}
0 & -1
\\
1 & 0
\end{pmatrix}
.
```
The GSTCs (with an ``\bm{R}`` and a 2 moved about) are
```math
\begin{align}
2 \bm{\rho}^e \bm{R} (\bm{H}_\parallel^2 - \bm{H}_\parallel^1)
&= \bm{E}_\parallel^2 + \bm{E}_\parallel^1
\\
\bm{E}_\parallel^2 - \bm{E}_\parallel^1
&= \bm{R} \bm{\sigma}^m (\bm{H}_\parallel^2 + \bm{H}_\parallel^1)/2
\end{align}
```
It is basic algebra to show that ``\bm{\Omega}^i`` commutes with ``\bm{\Omega}^j`` so multiplying the GSTCs by ``\bm{\Omega}^1 \bm{\Omega}^2 \bm{R}^{-1} = \bm{\Omega}^2 \bm{\Omega}^1 \bm{R}^{-1}`` gives
```math
\begin{align}
2 \bm{\Omega}^1 \bm{\Omega}^2 \bm{R}^{-1}
\bm{\rho}^e \bm{R} (\bm{H}_\parallel^2 - \bm{H}_\parallel^1)
&= \bm{\Omega}^1 \bm{R}^{-1} \bm{R} \bm{\Omega}^2 \bm{R}^{-1} \bm{E}_\parallel^2
+ \bm{\Omega}^2 \bm{R}^{-1} \bm{R} \bm{\Omega}^1 \bm{R}^{-1} \bm{E}_\parallel^1
\\
\bm{\Omega}^1 \bm{R}^{-1} \bm{R} \bm{\Omega}^2 \bm{R}^{-1} \bm{E}_\parallel^2
- \bm{\Omega}^2 \bm{R}^{-1} \bm{R} \bm{\Omega}^1 \bm{R}^{-1} \bm{E}_\parallel^1
&= \bm{\Omega}^1 \bm{\Omega}^2 \bm{\sigma}^m (\bm{H}_\parallel^2 + \bm{H}_\parallel^1)/2
\end{align}
```
Next we will do 3 operations in the same step:
1. Break each mode into forward and backward propagating components: ``\bm{H}_\parallel^i \to \bm{H}_\parallel^{(i,+)} + \bm{H}_\parallel^{(i,-)},\ \bm{E}_\parallel^i \to \bm{E}_\parallel^{(i,+)} + \bm{E}_\parallel^{(i,-)}``
2. Eliminate the electric field with Ampère's law: ``\bm{R} \bm{\Omega}^i \bm{R}^{-1} \bm{E}_\parallel^{(i,\pm)} = \pm \omega \bm{R} \mu^i k_z^i \bm{H}_\parallel^{(i,\pm)}``
3. Cancel any pairs of ``\bm{R}`` and ``\bm{R}^{-1}``
```math
\begin{align}
&2 \bm{\Omega}^1 \bm{\Omega}^2 \bm{R}^{-1} \bm{\rho}^e \bm{R} \left( \bm{H}_\parallel^{(2,+)} + \bm{H}_\parallel^{(2,-)} - \bm{H}_\parallel^{(1,+)} - \bm{H}_\parallel^{(1,-)} \right)
\\ &\qquad = \omega \left(
\bm{\Omega}^1 \mu^2 k_z^2 \left( \bm{H}_\parallel^{(2,+)} - \bm{H}_\parallel^{(2,-)} \right)
+ \bm{\Omega}^2 \mu^1 k_z^1 \left( \bm{H}_\parallel^{(1,+)} - \bm{H}_\parallel^{(1,-)} \right) \right)
\\
&\omega \left(
\bm{\Omega}^1 \mu^2 k_z^2 \left( \bm{H}_\parallel^{(2,+)} - \bm{H}_\parallel^{(2,-)} \right)
- \bm{\Omega}^2 \mu^1 k_z^1 \left( \bm{H}_\parallel^{(1,+)} - \bm{H}_\parallel^{(1,-)} \right) \right)
\\ &\qquad = \bm{\Omega}^1 \bm{\Omega}^2 \bm{\sigma}^m \left( \bm{H}_\parallel^{(2,+)} + \bm{H}_\parallel^{(2,-)} + \bm{H}_\parallel^{(1,+)} + \bm{H}_\parallel^{(1,-)} \right)/2
\end{align}
```
Next if we move the scattered components, (1,-) and (2,+), to the left and the incident components to the right, (1,+) and (2,-), we can express the scattering matrix as
```math
\begin{align}
&\begin{pmatrix}
-\bm{\Omega}^1 \bm{\Omega}^2 \bm{R}^{-1} \bm{\rho}^e \bm{R} + \frac{1}{2} \omega \bm{\Omega}^2 \mu^1 k_z^1
&
\bm{\Omega}^1 \bm{\Omega}^2 \bm{R}^{-1} \bm{\rho}^e \bm{R} - \frac{1}{2} \omega \bm{\Omega}^1 \mu^2 k_z^2
\\
-\bm{\Omega}^1 \bm{\Omega}^2 \bm{\sigma}^m + 2 \omega \bm{\Omega}^2 \mu^1 k_z^1
&
-\bm{\Omega}^1 \bm{\Omega}^2 \bm{\sigma}^m + 2 \omega \bm{\Omega}^1 \mu^2 k_z^2
\end{pmatrix}
\begin{pmatrix}
\bm{H}_\parallel^{(1,-)}
\\
\bm{H}_\parallel^{(2,+)}
\end{pmatrix}
\\ =
&\begin{pmatrix}
\bm{\Omega}^1 \bm{\Omega}^2 \bm{R}^{-1} \bm{\rho}^e \bm{R} + \frac{1}{2} \omega \bm{\Omega}^2 \mu^1 k_z^1
&
-\bm{\Omega}^1 \bm{\Omega}^2 \bm{R}^{-1} \bm{\rho}^e \bm{R} - \frac{1}{2} \omega \bm{\Omega}^1 \mu^2 k_z^2
\\
\bm{\Omega}^1 \bm{\Omega}^2 \bm{\sigma}^m + 2 \omega \bm{\Omega}^2 \mu^1 k_z^1
&
\bm{\Omega}^1 \bm{\Omega}^2 \bm{\sigma}^m + 2 \omega \bm{\Omega}^1 \mu^2 k_z^2
\end{pmatrix}
\begin{pmatrix}
\bm{H}_\parallel^{(1,+)}
\\
\bm{H}_\parallel^{(2,-)}
\end{pmatrix}
\end{align}
```
Note that if the two media are the same (i.e. ``\epsilon^1=\epsilon^2, \mu^1=\mu^2``) then ``\bm{\Omega}^1 = \bm{\Omega}^2`` and the scattering matrix has an overall factor of ``\Omega`` that can be removed, which recovers the result for the scattering matrix in a homogenous medium.
In fact we can manipulate this expression further using the specific form of the inverse
```math
(\bm{\Omega}^i)^{-1} = \bm{R}^{-1} \bm{\Omega}^i \bm{R} \frac{1}{\omega^2 \epsilon^i \mu^i (k_z^i)^2}
```
due to the fact ``\bm{\Omega}^i`` is a square block matrix whose blocks are diagonal matrices.
"""

# ╔═╡ 325b8359-ebac-434b-bba1-31aba05c2c39
md"
## Higher order corrections for a metasurface at a dielectric interface

The GSTCs we use are published by Epstein (DOI: 10.1109/TAP.2014.2354419)

These are based on ones published by Kuester (DOI: 10.1109/TAP.2003.817560), which
make the following assumptions.
'In this paper, we will derive what we believe to be the most flexible form of a
generalized sheet transition condition (GSTC) applicable to metafilms located in a
homogeneous medium. The case of a metafilm on a dielectric interface will be
considered in a separate paper.'

Kuester published a follow-up article (DOI: 10.1109/APS.2010.5562250)
with a first order correction to the previous one which should generalize to the
case where the metafilm has different media on either side.
The form of the GSTC is the same as before, except the impedance matrices are modified functions of the actual material parameters
"

# ╔═╡ Cell order:
# ╟─83d9578c-3960-486c-ad3d-b4449f7bae6f
# ╠═c24044bc-2fc3-4e34-b99c-7dc4ec6a7515
# ╠═42da96ce-8122-4aa0-99d6-635f68db6e66
# ╠═88e46392-1fbc-4079-8207-dd5b3ac01d7d
# ╠═6cbc5f34-fe14-45e8-9fee-9787d76ad0fe
# ╠═fb8df44c-a354-42de-9e99-df64dd4c5727
# ╠═f3d36c35-c09b-483f-bb41-046b6ecd0773
# ╠═e899dbc0-03bb-4a13-9c09-bfad33e5c1d3
# ╠═43d29086-510d-4870-94b7-2ff5ff34324b
# ╠═5377ae36-32fc-41c5-bd03-c737422b17f4
# ╠═65c2e966-8cb0-4d76-96db-18f9045695a3
# ╠═16316725-c608-49a5-b90c-8caff973683a
# ╠═bf1213cd-2ec1-4d1e-8ebc-5be5af26621c
# ╟─ba0c4060-f9d4-43b6-bc31-b48dcd81dd91
# ╟─687627f9-6938-4d7a-9395-9f1783b8a449
# ╟─860dc600-c86d-43a2-a88a-eb876eecbd6e
# ╟─271422af-db75-4373-be1c-104745f941bf
# ╟─325b8359-ebac-434b-bba1-31aba05c2c39


text:
"""
---
Apply Chebyshev polynom of order L to input x
# --- Syntax
	y = cheb(L,x)
# --- Input parameters
- L	  : Chebyshev order [Int]
- x	  : Input [Any]
# --- Output parameters
- y	  : Output
# ---
# v 1.0 - Robin Gerzaguet.
"""
function cheb(L,x)
	## Apply Chebyshev polynom of order m to input x
	 x[abs.(x).<=1] = cos.(L .* acos.(x[abs.(x).<=1]));
	 x[abs.(x).>1]  = cosh.(L .* acosh.(abs.(x[abs.(x).>1])));
	 return x;
end

"""
---
Returns the Dolp Chebyshev filter of order L with desired attenation att

See Peter Lynch, "The Dolph-Chebyshev Window: A Simple Optimal Filter", Monthly Weather Review, Vol. 125, pp. 655-660, April 1997. (http://www.maths.tcd.ie/~plynch/Publications/Dolph.pdf)
Dolph, "A current distribution for broadside arrays which optimizes the relationship between beam width and side-lobe level", Proc. IEEE, 34, pp. 335-348.
		cheb(m-1, beta * cos(pi * k/m))
W(k) =	-------------------------------
			  cheb(m-1, beta)
# --- Syntax
 filterDC = dolphChebyshev(n,at)
# --- Input parameters
- n		  : Size of desired filter
- at	  : Attenation in dBB
# --- Output parameters
- filterDC : Filter impulse response [Array{Float64}]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function dolphChebyshev(n,at);
	# --- Filter in frequency domain
	# --- DC constant
    gamma = 10^(-at/20);
    beta  = cosh(1/(n-1) * acosh(1/gamma));
	k     = collect(0:n-1);
    x     = beta*cos.(pi*k/n);
	p	  = cheb(n-1, x);
	if mod(n,2) == 1
		w = real(fft(p));
		M = Int((n+1)/2);
		w = w[1:M]/w[1];
		w = [w[M:-1:2];w];
	else
		off = exp.(1im*pi/n.*collect(0:n-1));
		# half-sample delay (even order)
		p .= p * transpose(off);
		w = real(fft(p));
		M = Int(n/2+1);
		w = w./w[2];
		w .= [w[M:-1:2];w[2:M]];
	end
	return w
end


text:
module BenchSum

using BenchmarkTools
using CUDA
using Folds
using FoldsCUDA

const CACHE = Ref{Any}(nothing)

fact(::Val{0}) = 1
fact(::Val{N}) where {N} = N * fact(Val(N - 1))

""" Some function with controllable complexity. """
function exp_approx(x, ::Val{N}) where {N}
    coefficients = ntuple(p -> 1.0 / fact(Val(p - 1)), Val(N))
    return evalpoly(x, coefficients)
end

function setup(; test = false, n = test ? 10 : 2^30, exp_approx_n = 10)
    CACHE[] = CUDA.randn(n)

    suite = BenchmarkGroup()

    s = suite["identity"] = BenchmarkGroup()
    s["base"] = @benchmarkable sum(CACHE[])
    s["folds"] = @benchmarkable Folds.sum(CACHE[])
    s["coalesced"] = @benchmarkable Folds.sum(CACHE[], CoalescedCUDAEx())

    ea = let ea, n = Val(exp_approx_n)
        ea(x) = exp_approx(x, n)
    end
    s = suite["exp_approx"] = BenchmarkGroup()
    s["base"] = @benchmarkable sum($ea, CACHE[])
    s["folds"] = @benchmarkable Folds.sum($ea, CACHE[])
    s["coalesced"] = @benchmarkable Folds.sum($ea, CACHE[], CoalescedCUDAEx())

    s = suite["sin"] = BenchmarkGroup()
    s["base"] = @benchmarkable sum(sin, CACHE[])
    s["folds"] = @benchmarkable Folds.sum(sin, CACHE[])
    s["coalesced"] = @benchmarkable Folds.sum(sin, CACHE[], CoalescedCUDAEx())

    return suite
end

function clear()
    xs = CACHE[]
    CACHE[] = nothing
    if xs isa CuArray
        CUDA.unsafe_free!(xs)
    end
end

end  # module BenchSum


text:
### Process examples

# Always rerun examples
const EXAMPLES_OUT = joinpath(@__DIR__, "src", "examples")

# State explicitly which examples are to be included, rather than including everything by
# default. This can make the docs development cycle faster -- it means that if you're not
# working on the examples then you can just make `examples` an empty vector after building
# all of the examples in a first run of the docs.
# Similarly, when developing only a subset of the examples, only a subset need to be
# re-built every time. To achieve this, comment out those examples that you don't wish to
# repeatedly re-build. Alternatively, if you're just working on a single example, either
# work on the source directly, or work with the literate.jl file to avoid needing to start
# a fresh Julia session each time you want to run the example.
examples = [
    "getting_started"
]

example_locations = map(example -> joinpath(@__DIR__, "..", "examples", example), examples)

if ispath(EXAMPLES_OUT)
    map(examples) do example
        path = joinpath(EXAMPLES_OUT, example)
        println(path)
        isdir(path) && rm(path; recursive=true)
    end
else
    mkpath(EXAMPLES_OUT)
end

let script = "using Pkg; Pkg.activate(ARGS[1]); Pkg.instantiate()"
    for example in examples
        if !success(`$(Base.julia_cmd()) -e $script $example`)
            error(
                "project environment of example ",
                basename(example),
                " could not be instantiated",
            )
        end
    end
end

# Run examples asynchronously
literate_path = joinpath(@__DIR__, "literate.jl")
processes = map(examples) do example
    return run(
        pipeline(
            `$(Base.julia_cmd()) $literate_path $(basename(example)) $EXAMPLES_OUT`;
            stdin=devnull,
            stdout=devnull,
            stderr=stderr,
        );
        wait=true,
    )::Base.Process
end

# Check that all examples were run successfully
isempty(processes) || success(processes) || error("some examples were not run successfully")

### Build documentation
using Pkg
Pkg.develop(path=joinpath(@__DIR__, ".."))
using Documenter, Stheno

DocMeta.setdocmeta!(
    Stheno,
    :DocTestSetup,
    :(using Stheno.AbstractGPs, Stheno, Random, LinearAlgebra);
    recursive=true,
)

makedocs(
    modules=[Stheno],
    format=Documenter.HTML(),
    sitename="Stheno.jl",
    pages=[
        "Home" => "index.md",
        "Getting Started" => joinpath("examples", "getting_started.md"),
        "Input Types" => "input_types.md",
        "Kernel Design" => "kernel_design.md",
        "Internals" => "internals.md",
        "API" => "api.md",
    ],
    doctestfilters=[
        r"{([a-zA-Z0-9]+,\s?)+[a-zA-Z0-9]+}",
        r"(Array{[a-zA-Z0-9]+,\s?1}|Vector{[a-zA-Z0-9]+})",
        r"(Array{[a-zA-Z0-9]+,\s?2}|Matrix{[a-zA-Z0-9]+})",
    ],
)

deploydocs(; repo="github.com/JuliaGaussianProcesses/Stheno.jl.git", push_preview=true)


text:
using Brainstorm

Brainstorm._Test.test()


text:
using Base.Test
using FastGaussQuadrature

include("test_gausslegendre.jl")
include("test_gaussjacobi.jl")
include("test_gaussradau.jl")
include("test_gausslobatto.jl")


text:
"""
    zipandgutzwiller(mps1, mps2, mode)

given two MPS, with the same number of sites and two physical
dimensions perform gutzwiller projection of the two by making the
matrices using iterative SVD. The final result is an MPS with 2
physical dimension where the final matrices are the tensor product of
the given ones according to the `mode` which can be :B14 for bosonic
or :F23 for fermionic.

"""

function zipandgutzwiller!(mps1::MPState{Y},
                           mps2::MPState{Y};
                           mode::Symbol=:B14,
                           maxdim::Int64=200) where {Y<:Array}
    if mode != :B14
        error("only Bosinic ↑↑ and ↓↓ is possible for no-symmetry MPSs!")
    end
    @assert mps1.d == mps2.d == 2
    lx = mps1.lx
    @assert mps2.lx == lx

    center_at!(mps1, 1)
    center_at!(mps2, 1)

    dims = ones(Int64, lx+1)
    matrices = Array{T, 3}[]

    # Gutzwiller projector for mode :B14
    gutzp = zeros(T, 2,2,2)
    gutzp[1,1,1] = 1
    gutzp[2,2,2] = 1

    E = ones(T, 1,1,1)
    for site=1:lx-1
        mat1 = mps1.matrices[site]
        mat2 = mps2.matrices[site]

        diml1, dimr1 = size(mat1, 1), size(mat1, 3)
        diml2, dimr2 = size(mat2, 1), size(mat2, 3)
        dimle = size(E, 1)

        @tensor C[le, dg, r2, r1] := ((E[le, l2, l1] * mat1[l1,d1,r1]) *
                                      gutzp[dg,d1,d2]) * mat2[l2,d2, r2]

        fact = svd(reshape(C, dimle*2, dimr2*dimr1), full=false)
        S, n, ratio = MPStateTools.truncate(fact.S, maxdim=maxdim)
        dims[site+1] = n
        U = fact.U[:, 1:n]
        push!(matrices, reshape(U, dimle, 2, n))
        E = reshape(Diagonal(S) * fact.Vt[1:n, :], n, dimr2, dimr1)
    end
    mat1 = mps1.matrices[lx]
    mat2 = mps2.matrices[lx]
    @tensor C[le, dg, r2, r1] := ((E[le, l2, l1] * mat1[l1,d1,r1]) *
                                  gutzp[dg,d1,d2]) * mat2[l2,d2,r2]
    push!(matrices, reshape(C, size(E, 1), 2, 1))

    return MPState{T}(lx, 2, dims, matrices, lx)
end

function zipandgutzwiller!(mps1::MPState{Y},
                           mps2::MPState{Y};
                           mode::Symbol=:F23,
                           maxdim::Int64=200) where {Y<:SymTensor}
    if mode==:B14
        return _zipandgutzwiller_B14!(mps1, mps2, maxdim=maxdim)
    elseif mode==:F23
        return _zipandgutzwiller_F23!(mps1, mps2, maxdim=maxdim)
    else
        error("mode not defined : ", mode)
    end
end

#bosonic version
function _zipandgutzwiller_B14!(mps1::MPState{Y},
                                mps2::MPState{Y};
                                maxdim::Int64=200) where {Y}
    #@assert mps1.d == mps2.d == 2
    lx = length(mps1)
    @assert length(mps2) == lx

    center_at!(mps1, 1)
    center_at!(mps2, 1)

    mps = MPState{Y}()

    ## NOTE: in order to make the gutzwiller projector respect the U1
    ## symmetry we need to do the follwoing. Assume that ↑↑
    ## corresponds to ↑ or 2 and ↓↓ corresponds or ↓ or 0.
    G = fill(one(T), 0,
             (STLeg(+1, [0,2], [1,1]),
              STLeg(-1, [0,1], [1,1]),
              STLeg(-1, [0,1], [1,1])))

    E = fill(one(T), 0, (STLeg(+1, [0], [1]),
                         STLeg(-1, [0], [1]),
                         STLeg(-1, [0], [1])))
    for l=1:lx-1
        A = mps1.matrices[l]
        B = mps2.matrices[l]

        ##NOTE:
        # contracting (E2, A1)
        # First tensor has EA = E1 E3 A2 A3
        # contracting (A2, G2)
        # Scond tensor has EAG = E1 E3 G1 A3 G3
        # contracting (E3, B1) and (G3, B2)
        # Final tensor has E1 G1 A3 B3
        C = contract(contract(contract(E, (1,-1, 2), A, (-1, 3, 4)),
                              (1, 2, -1, 4), G, (3, -1, 5)),
                     (1,-1, 2, 3,-2), B, (-1,-2, 4))

        u,s,v = svdtrunc(fuselegs(fuselegs(C, 1, 2), 2, 2), maxdim=maxdim)
        normalize!(s)
        push!(matrices,
              mapcharges(x->div(x,2),
                         unfuseleg(u, 1, (E.legs[1], G.legs[1]))))
        E = unfuseleg(s*v, 2, (A.legs[3], B.legs[3]))
    end
    A = mps1.As[lx]
    B = mps2.As[lx]
    C = contract(contract(contract(E, (1,-1, 2), A, (-1, 3, 4)),
                          (1, 2, -1, 4), G, (3, -1, 5)),
                 (1,-1, 2, 3,-2), B, (-1,-2, 4))

    push!(matrices, mapcharges(x->div(x,2), fuselegs(C, 3, 2)))

    return MPState{Y}(lx, 2, dims, matrices, lx)
end

function _zipandgutzwiller_F23!(mps1::MPState{Y},
                                mps2::MPState{Y};
                                truncation::Bool=true,
                                maxdim::Int64=200) where {Y}
    #@assert mps1.d == mps2.d == 2
    lx = length(mps1)
    @assert length(mps2) == lx

    center_at!(mps1, 1)
    center_at!(mps2, 1)

    mps = MPState{Y}()

    ## NOTE: in order to make the gutzwiller projector respect the U1
    ## symmetry we need to do the follwoing. Assume the first mps
    ## corresponds to ↑ or 1 and second mps to ↓ or -1.
    T = eltype(Y)
    S = vtype(Y)

    Vd = VectorSpace{S}(0=>1, 1=>1)
    G = fill(one(T), (dual(Vd), mapcharges(x->2*x-1, Vd), Vd))

    Vdummy = VectorSpace{S}(0=>1)
    E = fill(one(T), (Vdummy, dual(Vdummy), Vdummy))

    u,s,v = _svd_(eye(T, Vdummy))
    for l=1:lx-1
        A = mps1.As[l]
        B = dual(mps2.As[l], conjugate=false)

        fswap = fermionswapgate(A.space[2], B.space[1])
        ##NOTE:
        # contracting (E3, A1)
        # First tensor has EA = E1 E2 A2 A3
        # contracting with swap
        # second tensor has EAX = E1 XA2 E3 A3
        # contracting (A2, G2)
        # Third tensor has EAG = E1 G1 G3 E3 A3
        # contracting (E3, B1) and (G3, B2)
        # Final tensor has E1 G1 A3 B3
        C = contract(contract(contract(contract(E, (1, 2, -1), A, (-1, 3, 4)),
                                       (1, -1, -2, 4), fswap, (-1, 2, 3, -2)),
                              (1, -1, 4, 5), G, (-1, 2, 3)),
                     (1,2, -2, -1, 4), B, (-1,-2, 3))

        if truncation
            u,s,v = svdtrunc(SymMatrix(C, [1,2], [3,4]), maxdim=maxdim)
        else
            u,s,v = _svd_(SymMatrix(C, [1,2], [3,4]))
        end
        normalize!(s)

        fnl = x->div(x+l-1, 2)
        fnd = x->div(x+1, 2)
        fnr = x->div(x+l, 2)
        push!(mps,
              mapcharges((fnl,fnd,fnr),
                         splitleg(u, 1, (E.space[1], G.space[2]))))
        E = splitleg(s*v, 2, (B.space[3], A.space[3]))
    end
    A = mps1.As[lx]
    B =  dual(mps2.As[lx], conjugate=false)
    fswap = fermionswapgate(A.space[2], B.space[1])
    C = contract(contract(contract(contract(E, (1, 2, -1), A, (-1, 3, 4)),
                                   (1, -1, -2, 4), fswap, (-1, 2, 3, -2)),
                          (1, -1, 4, 5), G, (-1, 2, 3)),
                 (1,2, -2, -1, 4), B, (-1,-2, 3))

    if truncation
        u,s,v = svdtrunc(SymMatrix(C, [1,2], [3,4]), maxdim=maxdim)
    else
        u,s,v = _svd_(SymMatrix(C, [1,2], [3,4]))
    end
    normalize!(s)

    C = splitleg(u*s*v, 1, (E.space[1], G.space[2]))
    fnl = x->div(x+lx-1, 2)
    fnd = x->div(x+1, 2)
    fnr = x->div(x+lx, 2)
    push!(mps, mapcharges((fnl,fnd,fnr), C))

    # could the below alone be the issue?!
    # is this correct (need explanation and stuff)
    mps.center = lx
    mps
end

function gutzwillerexact(mps1::MPState{Y},
                         mps2::MPState{Y};
                         mode::Symbol=:F23) where {Y}
    _zipandgutzwiller_F23(mps1, mps2, truncation=false)
end

function _tensorproductzip!(mps1::MPState{Y},
                            mps2::MPState{Y};
                            truncation::Bool=true,
                            maxdim::Int64=200,
                            verbose::Bool=true) where {Y}

    #@assert mps1.d == mps2.d == 2
    lx = length(mps1)
    @assert length(mps2) == lx

    center_at!(mps1, 1)
    center_at!(mps2, 1)

    mps = MPState{Y}()

    ## NOTE: in order to make the gutzwiller projector respect the U1
    ## symmetry we need to do the follwoing. Assume the first mps
    ## corresponds to ↑ or 1 and second mps to ↓ or -1.
    T = eltype(Y)
    S = vtype(Y)

    Vdummy = VectorSpace{S}(0=>1)
    E = fill(one(T), (Vdummy, dual(Vdummy), Vdummy))

    u,s,v = _svd_(eye(T, Vdummy))
    for l=1:lx-1
        A = mps1.As[l]
        B = dual(mps2.As[l], conjugate=false)

        fswap = fermionswapgate(A.space[2], B.space[1])
        ##NOTE:
        # contracting (E3, A1)
        # First tensor has EA = E1 E2 A2 A3
        # contracting with swap
        # second tensor has EAX = E1 XA2 E2 A3
        # contracting (EX2, B1)
        # Final tensor has E1 XA2 B2 B3 A3
        C = contract(contract(contract(E, (1, 2, -1), A, (-1, 3, 4)),
                              (1, -1, -2, 4), fswap, (-1, 2, 3, -2)),
                     (1, 2, -1, 5), B, (-1, 3, 4))

        if truncation
            u,s,v = svdtrunc(SymMatrix(C, [1,2,3], [4,5]), maxdim=maxdim)
        else
            u,s,v = _svd_(SymMatrix(C, [1,2,3], [4,5]))
        end
        normalize!(s)

        push!(mps, fuselegs(
            splitleg(u, 1, (E.space[1], A.space[2], B.space[2])), 2, 2))

        E = splitleg(s*v, 2, (B.space[3], A.space[3]))
    end
    A = mps1.As[lx]
    B = dual(mps2.As[lx], conjugate=false)
    fswap = fermionswapgate(A.space[2], B.space[1])
    C = contract(contract(contract(E, (1, 2, -1), A, (-1, 3, 4)),
                          (1, -1, -2, 4), fswap, (-1, 2, 3, -2)),
                 (1, 2, -1, 5), B, (-1, 3, 4))

    push!(mps, fuselegs(fuselegs(C, 2, 2), 3, 2))
    mps.center = lx
    mps
end

function _applygutzwillerzip!(mps::MPState{Y}) where{Y}
    T = eltype(mps)
    G = Tensor(ones, zero(U1),
               (U1Space([-1=>1, 1=>1]),
                dual(U1Space([-1=>1, 0=>2, 1=>1]))
                ))

    mpsgutz = MPState{Y}()

    center_at!(mps, 1)
    C = contract(G, (2, -1), mps.As[1], (1, -1, 3))
    for l=1:lx-1
        u,s,v = svdtrunc(fuselegs(C, 1, 2))

        fnl = x->div(x+lx-1, 2)
        fnd = x->div(x+1, 2)
        fnr = x->div(x+lx, 2)
        push!(mpsgutz, mapcharges((fnl,fnd,fnr),
                                  splitleg(u, 1, (C.space[1], G.space[1]))))
        E = s*v
        C = contract(E, (1, -1),
                     contract(G, (2, -1), mps.As[l+1], (1, -1, 3)),
                     (-1,2,3))
    end
    push!(mpsgutz, C)

    mpsgutz
end

function _zipandgutzwiller_F23_analysis!(mps1::MPState{Y},
                                         mps2::MPState{Y};
                                         maxdim::Int64=200) where {Y}

    lx = length(mps1)
    @assert length(mps2) == lx

    center_at!(mps1, 1)
    center_at!(mps2, 1)

    mps = MPState{Y}()
    ## NOTE: in order to make the gutzwiller projector respect the U1
    ## symmetry we need to do the follwoing. Assume the first mps
    ## corresponds to ↑ or 1 and second mps to ↓ or -1.
    T = eltype(Y)
    S = vtype(Y)

    Vd = VectorSpace{S}(0=>1, 1=>1)
    G = fill(one(T), (dual(Vd), mapcharges(x->2*x-1, Vd), Vd))

    Vdummy = VectorSpace{S}(0=>1)

    ### Making the F tensors
    Vend1 = rightspace(mps1)
    Vend2 = rightspace(mps2)
    Fs = Vector{SymTensor{S,T,3}}(undef, lx)
    F = fill(one(T), (Vend1, dual(Vend2), fuse(true, dual(Vend1), Vend2)))
    for l=lx:-1:2
        A = mps1.As[l]
        B = dual(mps2.As[l], conjugate=false)
        fswap = fermionswapgate(A.space[2], B.space[1])
        C = contract(contract(contract(contract(A, (1, 2, -1), F, (-1, 3, 4)),
                                       (1, 2, -1, 5), B, (3, 4, -1)),
                              (1, -2, -1, 4, 5), fswap, (2, 3, -1, -2)),
                     (1, 2, -1, -2, 4), G, (-1, 3, -2))
        u,s,v = _svd_(SymMatrix(C, [1,2], [3,4]))
        normalize!(s)
        F = splitleg(u*s, 1, (C.space[1], C.space[2]))
        Fs[l] = F
    end

    error_ortho = Vector{Float64}()
    error_ets = Vector{Float64}()
    E = fill(one(T), (Vdummy, Vdummy, dual(Vdummy)))
    for l=1:lx-1
        A = mps1.As[l]
        B = dual(mps2.As[l], conjugate=false)

        fswap = fermionswapgate(A.space[2], B.space[1])
        ##NOTE:
        # contracting (E3, A1)
        # First tensor has EA = E1 E2 A2 A3
        # contracting with swap
        # second tensor has EAX = E1 XA2 E3 A3
        # contracting (A2, G2)
        # Third tensor has EAG = E1 G1 G3 E3 A3
        # contracting (E3, B1) and (G3, B2)
        # Final tensor has E1 G1 A3 B3
        C = contract(contract(contract(contract(E, (1, 2, -1), A, (-1, 3, 4)),
                                       (1, -1, -2, 4), fswap, (-1, 2, 3, -2)),
                              (1, -1, 4, 5), G, (-1, 2, 3)),
                     (1,2, -2, -1, 4), B, (-1,-2, 3))

        # actual zipgutz process
        u,s,v = svdtrunc(SymMatrix(C, [1,2], [3,4]), maxdim=maxdim)
        normalize!(s)
        println(norm(SymMatrix(C, [1,2], [3,4])))
        # The ortho center matrix
        CFmat = SymMatrix(contract(C, (1,2,-1,-2), Fs[l+1], (-2,-1, 3)), [1, 2], [3])
        println(norm(SymMatrix(Fs[l+1], [2,1], [3])))
        norm_CFmat = norm(CFmat)
        println(norm_CFmat)
        println()
        # othonormal truncation error
        _u,_s,_v = svdtrunc(CFmat, maxdim=maxdim)
        normalize!(_s)
        push!(error_ortho, dot(_u*_s*_v, CFmat) / norm_CFmat)

        # zipgutz truncation error
        CtrFmat = (u*s*v) * SymMatrix(Fs[l+1], [2,1], [3])
        norm_CtrFmat = norm(CtrFmat)
        push!(error_ets, dot(CtrFmat, CFmat) / (norm_CFmat*norm_CtrFmat))

        fnl = x->div(x+l-1, 2)
        fnd = x->div(x+1, 2)
        fnr = x->div(x+l, 2)
        push!(mps,
              mapcharges((fnl,fnd,fnr),
                         splitleg(u, 1, (E.space[1], G.space[2]))))
        E = splitleg(s*v, 2, (B.space[3], A.space[3]))
    end
    A = mps1.As[lx]
    B =  dual(mps2.As[lx], conjugate=false)
    fswap = fermionswapgate(A.space[2], B.space[1])
    C = contract(contract(contract(contract(E, (1, 2, -1), A, (-1, 3, 4)),
                                   (1, -1, -2, 4), fswap, (-1, 2, 3, -2)),
                          (1, -1, 4, 5), G, (-1, 2, 3)),
                 (1,2,-2,-1,4), B, (-1,-2, 3))

    u,s,v = svdtrunc(SymMatrix(C, [1,2], [3,4]), maxdim=maxdim)
    normalize!(s)

    push!(error_ortho, 1.0)
    push!(error_ets, 1.0)

    C = splitleg(u*s*v, 1, (E.space[1], G.space[2]))
    fnl = x->div(x+lx-1, 2)
    fnd = x->div(x+1, 2)
    fnr = x->div(x+lx, 2)
    push!(mps, mapcharges((fnl,fnd,fnr), C))

    mps.center = lx
    mps, error_ortho, error_ets
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
@testset "Validate" begin
    import ONNXNaiveNASflux: modelproto, graphproto
    import ONNXNaiveNASflux: validate, uniqueoutput, optypedefined, outputused, inputused, hasname, ONNX.NodeProto

    np(ins,outs,op="A") = ONNX.NodeProto(input=ins, output=outs, op_type=op)
    vip(name,shape=(2,:A)) = ONNX.ValueInfoProto(name, shape)

    @testset "Duplicate output name" begin
        gp = graphproto()

        push!(gp.node, np(["in"], ["n1"]))
        push!(gp.node, np(["n1"], ["n2"]))
        push!(gp.node, np(["n2"], ["n1", "n3"]))

        mp = modelproto()
        mp.graph = gp
        @test_throws ErrorException uniqueoutput(mp)
        @test_throws ErrorException validate(mp)
        @test_logs (:warn, r"Duplicate output name: n1 found in NodeProto.*in.*n1.* and NodeProto.*n2.*n1.*n3") uniqueoutput(mp, s -> @warn replace(s, "\n" => ""))
    end

    @testset "No OP specified" begin
        gp = graphproto()

        push!(gp.node, ONNX.NodeProto(output=["n1"], input=["in"]))
        mp = modelproto()
        mp.graph = gp
        @test_throws ErrorException optypedefined(mp)
        @test_throws ErrorException validate(mp)
        @test_logs (:warn, r"No op_type defined.*in.*n1") optypedefined(mp, s -> @warn replace(s, "\n" => ""))
    end

    @testset "All outputs not used" begin
        gp = graphproto()

        push!(gp.node, np(["in"], ["n1", "nf1"]))
        push!(gp.node, np(["n1"], ["n2", "nf2"]))
        push!(gp.node, np(["n2"], ["n3"]))
        gp.output = vip.(["n3", "n1"])
        gp.input = vip.(["in", "nf3"])

        mp = modelproto()
        mp.graph = gp
        @test_throws ErrorException outputused(mp)
        @test_throws ErrorException validate(mp)
        @test_logs (:warn, "Found unused outputs: nf1, nf2, nf3") outputused(mp, s -> @warn s)
    end

    @testset "All inputs not used" begin
        gp = graphproto()

        push!(gp.node, np(["in"], ["n1"]))
        push!(gp.node, np(["n1", "nf1"], ["n2"]))
        push!(gp.node, np(["n2"], ["n3"]))
        gp.output = vip.(["n3", "n1", "nf2"])
        gp.input = vip.(["in"])

        mp = modelproto()
        mp.graph = gp
        @test_throws ErrorException inputused(mp)
        @test_throws ErrorException validate(mp)
        @test_logs (:warn, "Found unused inputs: nf1, nf2") inputused(mp, s -> @warn s)
    end

    @testset "Graph has a name" begin
        gp = ONNX.GraphProto(;node=ONNX.NodeProto[],input=ONNX.ValueInfoProto[], output=ONNX.ValueInfoProto[],initializer=ONNX.TensorProto[],)
        mp = modelproto()
        mp.graph = gp
        @test_throws ErrorException hasname(mp)
        @test_throws ErrorException validate(mp)
        @test_logs (:warn, "Graph name not defined!") hasname(mp, s -> @warn s)

        # Empty string is not valid either
        gp.name = ""
        @test_throws ErrorException hasname(mp)
        @test_throws ErrorException validate(mp)
        @test_logs (:warn, "Graph name is empty string!") hasname(mp, s -> @warn s)
    end
end


text:
@testset "Estimator" begin

D = Mod1Dynamic(x->2*x)
B = Ulam(32)
P = assemble(B, D)

# C = InvariantMeasures.boundnorm(B, P, 10)

# @test C == [1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

end


text:
import Base: copy, copyto!, ==
import Random: rand, seed!
import RandomNumbers: AbstractRNG, gen_seed, seed_type

"""
```julia
SplitMix64 <: AbstractRNG{UInt64}
SplitMix64([seed])
```

Used for initializing a random seed.
"""
mutable struct SplitMix64 <: AbstractRNG{UInt64}
    x::UInt64
end
function SplitMix64(seed::Integer=gen_seed(UInt64))
    r = SplitMix64(zero(UInt64))
    seed!(r, seed)
    r
end

@inline seed_type(::Type{SplitMix64}) = UInt64

function copyto!(dest::SplitMix64, src::SplitMix64)
    dest.x = src.x
    dest
end

copy(src::SplitMix64) = copyto!(SplitMix64(zero(UInt64)), src)

==(r1::SplitMix64, r2::SplitMix64) = r1.x == r2.x

function seed!(r::SplitMix64, seed::Integer=gen_seed(UInt64))
    r.x = splitmix64(seed % UInt64)
    r
end

@inline function rand(r::SplitMix64, ::Type{UInt64})
    r.x = splitmix64(r.x)
end

@inline function splitmix64(x::UInt64)
    x += 0x9e3779b97f4a7c15
    x = (x ⊻ (x >> 30)) * 0xbf58476d1ce4e5b9
    x = (x ⊻ (x >> 27)) * 0x94d049bb133111eb
    x ⊻ (x >> 31)
end

init_seed(seed, ::Type{UInt64}) = splitmix64(seed % UInt64)
function init_seed(seed, ::Type{UInt64}, N::Int)
    x = seed % UInt64
    NTuple{N, UInt64}(
        x = splitmix64(x)
        for _ in 1:N
    )
end
function init_seed(seed, ::Type{UInt32}, N::Int)
    x = seed % UInt64
    NTuple{N, UInt32}(begin
        x = splitmix64(x)
        x % UInt32
    end for _ in 1:N)
end


text:
function Base.cat(a::NamedDimsArray{L}; dims) where L
    newL = expand_dimnames(L, dims)
    numerical_dims = dim(newL, dims)
    data = Base.cat(parent(a); dims=numerical_dims) # Base.cat is type unstable
    T = eltype(a)  # therefore the element type has to be inferred manually
    N = length(newL)  # as must the size of the array
    return NamedDimsArray{newL, T, N, Array{T,N}}(data)
end

# While the following two functions are covered by the general case below where splatted c
# is empty, they are included because they reduce allocations to that of regular arrays.
function Base.cat(a::NamedDimsArray{L}, b::AbstractArray; dims) where L
    newL = expand_dimnames(L, dims)
    numerical_dims = dim(newL, dims)
    data = Base.cat(parent(a), b; dims=numerical_dims)
    T = promote_type(eltype(a), eltype(b))
    N = length(newL)
    return NamedDimsArray{newL, T, N, Array{T,N}}(data)
end

function Base.cat(a::AbstractArray, b::NamedDimsArray{L}; dims) where L
    newL = expand_dimnames(L, dims)
    numerical_dims = dim(newL, dims)
    data = Base.cat(a, parent(b); dims=numerical_dims)
    T = promote_type(eltype(a), eltype(b))
    N = length(newL)
    return NamedDimsArray{newL, T, N, Array{T,N}}(data)
end

# to dispatch on the first _or the second_ argument being the NDA.
for (T, S) in [
    (:NamedDimsArray, :AbstractArray),
    (:AbstractArray, :NamedDimsArray),
    (:NamedDimsArray, :NamedDimsArray)
    ]

    @eval function Base.cat(a::$T, b::$S, c::AbstractArray...; dims)
        combL = unify_names_shortest(dimnames(a), dimnames(b), ntuple(i->dimnames(c[i]), length(c))...)
        newL = expand_dimnames(combL, dims)
        numerical_dims = dim(newL, dims)
        data = Base.cat(parent(a), parent(b), ntuple(i->parent(c[i]), length(c))...; dims=numerical_dims)
        T = promote_type(eltype(a), eltype(b), ntuple(i->eltype(c[i]), length(c))...)
        N = length(newL)
        return NamedDimsArray{newL, T, N, Array{T,N}}(data)
    end
end

# vcat and hcat
function Base.hcat(a::NamedDimsArray{L}) where L
    newL = expand_dimnames(L, 2)
    data = Base.hcat(parent(a))
    T = eltype(a)
    N = length(newL)
    return NamedDimsArray{newL, T, N, Array{T,N}}(data)
end

Base.vcat(a::NamedDimsArray{L}) where L = a

# Base.hcat and Base.vcat specialise on this Union
const AbsVecOrMat = Union{AbstractVector, AbstractMatrix}
for (T, S) in [
    (:NamedDimsArray, :AbsVecOrMat),
    (:NamedDimsArray, :AbstractArray),
    (:AbsVecOrMat, :NamedDimsArray),
    (:AbstractArray, :NamedDimsArray),
    (:NamedDimsArray, :NamedDimsArray)
    ]

    for (fun, d) in zip((:vcat, :hcat), (1, 2))
        # this special case is needed to resolve ambiguities
        @eval function Base.$fun(a::$T, b::$S)
            combL = unify_names_shortest(dimnames(a), dimnames(b))
            newL = expand_dimnames(combL, $d)
            data = Base.$fun(parent(a), parent(b))
            T = promote_type(eltype(a), eltype(b))
            N = length(newL)
            return NamedDimsArray{newL, T, N, Array{T,N}}(data)
        end

        @eval function Base.$fun(a::$T, b::$S, c::NamedDimsArray...)
            combL = unify_names_shortest(dimnames(a), dimnames(b), ntuple(i->dimnames(c[i]), length(c))...)
            newL = expand_dimnames(combL, $d)
            data = Base.$fun(parent(a), parent(b), ntuple(i->parent(c[i]), length(c))...)
            T = promote_type(eltype(a), eltype(b), ntuple(i->eltype(c[i]), length(c))...)
            N = length(newL)
            return NamedDimsArray{newL, T, N, Array{T,N}}(data)
        end
    end
end


text:
using ImplicitBias
using Test

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

max_stars_repo_name: brechetp/ImplicitBias

max_issues_repo_name: brechetp/ImplicitBias

max_forks_repo_name: brechetp/ImplicitBias


text:
using Test
using OrdinaryDiffEq
using SparseArrays
using LinearAlgebra

## in-place
#https://github.com/JuliaDiffEq/SparseDiffTools.jl/blob/master/test/test_integration.jl
function f_ip(dx,x,p,t)
  for i in 2:length(x)-1
    dx[i] = x[i-1] - 2x[i] + x[i+1]
  end
  dx[1] = -2x[1] + x[2]
  dx[end] = x[end-1] - 2x[end]
  nothing
end

## out-of-place
function f_oop(x,p,t)
  dx = similar(x)
  for i in 2:length(x)-1
    dx[i] = x[i-1] - 2x[i] + x[i+1]
  end
  dx[1] = -2x[1] + x[2]
  dx[end] = x[end-1] - 2x[end]
  return dx
end

function second_derivative_stencil(N)
  A = zeros(N,N)
  for i in 1:N, j in 1:N
    (j-i==-1 || j-i==1) && (A[i,j]=1)
    j-i==0 && (A[i,j]=-2)
  end
  A
end

function generate_sparsity_pattern(N::Integer)
  dl = repeat([1.0],N-1)
  du = repeat([1.0],N-1)
  d = repeat([-2.0],N)
  return Tridiagonal(dl,d,du)
end

jac_sp = sparse(generate_sparsity_pattern(10))
#jac = second_derivative_stencil(10)
colors = repeat(1:3,10)[1:10]
u0=[1.,2.,3,4,5,5,4,3,2,1]
tspan=(0.,10.)

for f in [f_oop, f_ip]
  odefun_std = ODEFunction(f)
  prob_std = ODEProblem(odefun_std,u0,tspan)

  for ad in [true, false]
    for Solver in [Rodas5, Trapezoid, KenCarp4]
      for tol in [nothing, 1e-10]
        sol_std=solve(prob_std,Solver(autodiff=ad),reltol=tol,abstol=tol)
        @test sol_std.retcode==:Success
        for (i,prob) in enumerate(map(f->ODEProblem(f,u0,tspan),
                        [ODEFunction(f,colorvec=colors,jac_prototype=jac_sp),
                         ODEFunction(f,jac_prototype=jac_sp),
                         ODEFunction(f,colorvec=colors,sparsity=jac_sp)
                         ]))
          sol=solve(prob,Solver(autodiff=ad),reltol=tol,abstol=tol)
          @test sol.retcode==:Success
          if tol !=nothing
            @test sol_std.u[end]≈sol.u[end] atol=tol
          else
            @test sol_std.u[end]≈sol.u[end]
          end
          @test length(sol_std.t)==length(sol.t)
        end
      end
    end
  end
end


text:
# Unwrap vecs into tuples
data(x) = x
abstract type AbstractVec{T,N} end

datatype(x) = typeof(x)
datatype(::AbstractVec{T}) where T = T

data(x::AbstractVec) = error("`data` not implemented for $(typeof(x))")

Base.length(xs::AbstractVec{T,N}) where {T,N} = N
Base.getindex(xs::AbstractVec, i::Integer) = data(xs)[i]

Base.getindex(xs::AbstractVec, is::AbstractVec) = vect(map(i -> xs[i], is)...)

function Base.show(io::IO, v::AbstractVec)
  print(io, summary(v), "{")
  join(io, sprint.(show, v), ", ")
  print(io, "}")
end

Base.iterate(v::AbstractVec, i = 1) =
  i > length(v) ? nothing : (v[i], i+1)

const Mask{N} = AbstractVec{Bool,N}

struct Vec{T,N} <: AbstractVec{T,N}
  data::NTuple{N,T}
  Vec{T,N}(data::NTuple{N,T}) where {T,N} = new(data)
end

data(vec::Vec) = vec.data

Vec(xs::Tuple) = Vec{Union{map(typeof,xs)...},length(xs)}(xs)

vect(xs...) = Vec(xs)
vect(xs::Nothing...) = Vec(xs)

Base.summary(::Vec{T}) where T = "Vec{$T}"

struct BitVec{N,T<:Unsigned} <: AbstractVec{Bool,N}
  data::T
end

BitVec{N}(data::T) where {T<:Unsigned,N} = BitVec{N,T}(data)

Base.summary(::BitVec) = "BitVec"

Base.getindex(v::BitVec, i) = Bool(v.data >> (i-1) & 0x01)

@generated function bitpack(::Type{T}, xs::NTuple{N,Any}) where {T<:Unsigned,N}
  y = zero(T)
  for i = 1:N
    y = :($y | (xs[$i] << $(i-1)))
  end
  return y
end

# Define some mask operations directly

import Base: &, |, ~

(a::Mask{N} & b::Mask{N}) where N = spmd(mask(N), &, a, b)
(a::Mask{N} | b::Mask{N}) where N = spmd(mask(N), |, a, b)
(~a::Mask{N}) where N = spmd(mask(N), ~, a)

Base.sum(x::AbstractVec) = sum(data(x))

# TODO use smaller words where possible
# vect(xs::Bool...) = BitVec{length(xs)}(bitpack(UInt64, xs))

import Base.convert, Base.promote_rule

convert(::Type{Vec{T,N}}, xs::Vec{S,N}) where {T,S,N} = Vec(convert.(T, xs.data))

promote_rule(::Type{Vec{S,N}}, ::Type{Vec{T,N}}) where {S,T,N} = Vec{promote_type(T,S),N}

#HACK
spmd(mask, ::typeof(===), x, y) = vect((x .=== y)...)


text:
using Base.Dates
using Base.Test

using MarketData

using TimeSeries

@testset "meta" begin

@testset "construction with and without meta field" begin
    nometa = TimeArray(cl.timestamp, cl.values, cl.colnames)

    @testset "default meta field to nothing" begin
        @test nometa.meta == nothing
    end

    @testset "allow objects in meta field" begin
        @test mdata.meta == "Apple"
    end
end

@testset "get index operations preserve meta" begin
    @testset "index by integer row" begin
        @test mdata[1].meta == "Apple"
    end

    @testset "index by integer range" begin
        @test mdata[1:2].meta == "Apple"
    end

    @testset "index by column name" begin
        @test mdata["Close"].meta == "Apple"
    end

    @testset "index by date range" begin
        @test mdata[[Date(2000,1,3), Date(2000,1,14)]].meta == "Apple"
    end
end

@testset "split operations preserve meta" begin
    @testset "when" begin
        @test when(mdata, dayofweek, 1).meta == "Apple"
    end

    @testset "from" begin
        @test from(mdata, Date(2000,1,1)).meta == "Apple"
    end

    @testset "to" begin
        @test to(mdata, Date(2000,1,1)).meta == "Apple"
    end
end

@testset "apply operations preserve meta" begin
    @testset "lag" begin
        @test lag(mdata).meta == "Apple"
    end

    @testset "lead" begin
        @test lead(mdata).meta == "Apple"
    end

    @testset "percentchange" begin
        @test percentchange(mdata).meta == "Apple"
    end

    @testset "moving" begin
        @test moving(mdata,mean,10).meta == "Apple"
    end

    @testset "upto" begin
        @test upto(mdata,sum).meta == "Apple"
    end
end

@testset "combine operations preserve meta" begin
    @testset "merge when both have identical meta" begin
        @test merge(cl, op).meta         == "AAPL"
        @test merge(cl, op, :left).meta  == "AAPL"
        @test merge(cl, op, :right).meta == "AAPL"
        @test merge(cl, op, :outer).meta == "AAPL"
    end

    @testset "merged meta field value concatenates when both objects' meta field values are strings" begin
        @test merge(mdata, cl).meta         == "Apple_AAPL"
        @test merge(mdata, cl, :left).meta  == "Apple_AAPL"
        @test merge(mdata, cl, :right).meta == "Apple_AAPL"
        @test merge(mdata, cl, :outer).meta == "Apple_AAPL"
    end

    @testset "merge when supplied with meta" begin
        @test merge(mdata, mdata, meta=47).meta         == 47
        @test merge(mdata, mdata, :left, meta=47).meta  == 47
        @test merge(mdata, mdata, :right, meta=47).meta == 47
        @test merge(mdata, mdata, :outer, meta=47).meta == 47
        @test merge(mdata, cl, meta=47).meta            == 47
        @test merge(mdata, cl, :left, meta=47).meta     == 47
        @test merge(mdata, cl, :right, meta=47).meta    == 47
        @test merge(mdata, cl, :outer, meta=47).meta    == 47
    end

    @testset "merged meta field value for disparate types in meta field defaults to Void" begin
        @test merge(mdata, merge(cl, op, meta=47)).meta         == Void
        @test merge(mdata, merge(cl, op, meta=47), :left).meta  == Void
        @test merge(mdata, merge(cl, op, meta=47), :right).meta == Void
        @test merge(mdata, merge(cl, op, meta=47), :outer).meta == Void
    end

    @testset "collapse" begin
        @test collapse(mdata, week, first).meta == "Apple"
    end
end

@testset "basecall operations preserve meta" begin
    @testset "basecall" begin
        @test basecall(mdata, cumsum).meta == "Apple"
    end
end

@testset "mathematical and comparison operations preserve meta" begin
    @testset ".+" begin
        @test (mdata .+ mdata).meta == "Apple"
        @test (mdata .+ cl).meta == Void
    end

    @testset ".<" begin
        @test (mdata .< mdata).meta == "Apple"
        @test (mdata .< cl).meta == Void
    end
end

@testset "readwrite accepts meta argument" begin
    @testset "Apple is present" begin
        @test mdata.meta == "Apple"
    end
end

end  # @testset "meta"


text:
let
    lrec = ListRecord(1.0, Float64, Bool, Int)
    append!(lrec.frames, [RecordFrame(1,2), RecordFrame(3,4)])
    append!(lrec.states, [RecordState(1.0,1),RecordState(2.0,2),RecordState(3.0,1),RecordState(4.0,3)])
    lrec.defs[1] = true
    lrec.defs[2] = true
    lrec.defs[3] = false

    file = tempname()
    open(file, "w") do io
        write(io, MIME"text/plain"(), lrec)
    end

    lrec2 = open(file, "r") do io
        read(io, MIME"text/plain"(), ListRecord{Float64, Bool, Int})
    end
    @test nframes(lrec2) == 2
    @test lrec2.frames == [RecordFrame(1,2), RecordFrame(3,4)]
    @test lrec2.states == [RecordState(1.0,1),RecordState(2.0,2),RecordState(3.0,1),RecordState(4.0,3)]
    @test lrec2.defs[1] == true
    @test lrec2.defs[2] == true
    @test lrec2.defs[3] == false
    @test length(lrec2.defs) == 3

    rm(file)
end

let
    frames = [
            Frame([Entity(1.0,true,1), Entity(2.0, true,2)]),
            Frame([Entity(3.0,true,1), Entity(4.0,false,3)]),
            ]

    file = tempname()
    open(file, "w") do io
        write(io, MIME"text/plain"(), frames)
    end

    frames2 = open(file, "r") do io
        read(io, MIME"text/plain"(), Vector{Frame{Entity{Float64, Bool, Int}}})
    end
    @test frames2[1][1] == Entity(1.0,true,1)
    @test frames2[1][2] == Entity(2.0,true,2)
    @test frames2[2][1] == Entity(3.0,true,1)
    @test frames2[2][2] == Entity(4.0,false,3)
    @test length(frames2) == 2

    rm(file)
end

# let
#     qrec = QueueRecord(Entity{Float64, Bool, Int}, 2, 0.1)
#     @test nframes(qrec) == 0

#     update!(qrec, Frame([Entity(1.0,true,1), Entity(2.0,true,2)]))
#     update!(qrec, Frame([Entity(3.0,true,1), Entity(4.0,false,3)]))

#     file = tempname()
#     open(file, "w") do io
#         write(io, MIME"text/plain"(), qrec)
#     end

#     qrec2 = open(file, "r") do io
#         read(io, MIME"text/plain"(), QueueRecord{Entity{Float64, Bool, Int}})
#     end
#     @test qrec2[-1][1] == Entity(1.0,true,1)
#     @test qrec2[-1][2] == Entity(2.0,true,2)
#     @test qrec2[ 0][1] == Entity(3.0,true,1)
#     @test qrec2[ 0][2] == Entity(4.0,false,3)
#     @test nframes(qrec2) == 2

#     rm(file)
# end

max_stars_repo_name: tawheeler/Records

max_issues_repo_name: tawheeler/Records

max_forks_repo_name: tawheeler/Records


text:
"""
    emd(μ, ν, C; kwargs...)

Compute the optimal transport plan for the Monge-Kantorovich problem with source and target
marginals `μ` and `ν` and cost matrix `C` of size `(length(μ), length(ν))`.

The optimal transport plan `γ` is of the same size as `C` and solves
```math
\\inf_{\\gamma \\in \\Pi(\\mu, \\nu)} \\langle \\gamma, C \\rangle.
```

This function is a wrapper of the function
[`emd`](https://pythonot.github.io/all.html#ot.emd) in the Python Optimal Transport
package. Keyword arguments are listed in the documentation of the Python function.

# Examples

```jldoctest
julia> μ = [0.5, 0.2, 0.3];

julia> ν = [0.0, 1.0];

julia> C = [0.0 1.0; 2.0 0.0; 0.5 1.5];

julia> emd(μ, ν, C)
3×2 Matrix{Float64}:
 0.0  0.5
 0.0  0.2
 0.0  0.3
```

See also: [`emd2`](@ref)
"""
function emd(μ, ν, C; kwargs...)
    return pot.lp.emd(μ, ν, PyCall.PyReverseDims(permutedims(C)); kwargs...)
end

"""
    emd2(μ, ν, C; kwargs...)

Compute the optimal transport cost for the Monge-Kantorovich problem with source and target
marginals `μ` and `ν` and cost matrix `C` of size `(length(μ), length(ν))`.

The optimal transport cost is the scalar value
```math
\\inf_{\\gamma \\in \\Pi(\\mu, \\nu)} \\langle \\gamma, C \\rangle.
```

This function is a wrapper of the function
[`emd2`](https://pythonot.github.io/all.html#ot.emd2) in the Python Optimal Transport
package. Keyword arguments are listed in the documentation of the Python function.

# Examples

```jldoctest
julia> μ = [0.5, 0.2, 0.3];

julia> ν = [0.0, 1.0];

julia> C = [0.0 1.0; 2.0 0.0; 0.5 1.5];

julia> emd2(μ, ν, C)
0.95
```

See also: [`emd`](@ref)
"""
function emd2(μ, ν, C; kwargs...)
    return pot.lp.emd2(μ, ν, PyCall.PyReverseDims(permutedims(C)); kwargs...)
end

"""
    emd_1d(xsource, xtarget; kwargs...)

Compute the optimal transport plan for the Monge-Kantorovich problem with univariate
discrete measures with support `xsource` and `xtarget` as source and target marginals.

This function is a wrapper of the function
[`emd_1d`](https://pythonot.github.io/all.html#ot.emd_1d) in the Python Optimal Transport
package. Keyword arguments are listed in the documentation of the Python function.

# Examples

```jldoctest
julia> xsource = [0.2, 0.5];

julia> xtarget = [0.8, 0.3];

julia> emd_1d(xsource, xtarget)
2×2 Matrix{Float64}:
 0.0  0.5
 0.5  0.0

julia> histogram_source = [0.8, 0.2];

julia> histogram_target = [0.7, 0.3];

julia> emd_1d(xsource, xtarget; a=histogram_source, b=histogram_target)
2×2 Matrix{Float64}:
 0.5  0.3
 0.2  0.0
```

See also: [`emd`](@ref), [`emd2_1d`](@ref)
"""
function emd_1d(xsource, xtarget; kwargs...)
    return pot.lp.emd_1d(xsource, xtarget; kwargs...)
end

"""
    emd2_1d(xsource, xtarget; kwargs...)

Compute the optimal transport cost for the Monge-Kantorovich problem with univariate
discrete measures with support `xsource` and `xtarget` as source and target marginals.

This function is a wrapper of the function
[`emd2_1d`](https://pythonot.github.io/all.html#ot.emd2_1d) in the Python Optimal Transport
package. Keyword arguments are listed in the documentation of the Python function.

# Examples

```jldoctest
julia> xsource = [0.2, 0.5];

julia> xtarget = [0.8, 0.3];

julia> round(emd2_1d(xsource, xtarget); sigdigits=6)
0.05

julia> histogram_source = [0.8, 0.2];

julia> histogram_target = [0.7, 0.3];

julia> round(emd2_1d(xsource, xtarget; a=histogram_source, b=histogram_target); sigdigits=6)
0.201
```

See also: [`emd2`](@ref), [`emd2_1d`](@ref)
"""
function emd2_1d(xsource, xtarget; kwargs...)
    return pot.lp.emd2_1d(xsource, xtarget; kwargs...)
end

"""
    sinkhorn(μ, ν, C, ε; kwargs...)

Compute the optimal transport plan for the entropic regularization optimal transport problem
with source and target marginals `μ` and `ν`, cost matrix `C` of size
`(length(μ), length(ν))`, and entropic regularization parameter `ε`.

The optimal transport plan `γ` is of the same size as `C` and solves
```math
\\inf_{\\gamma \\in \\Pi(\\mu, \\nu)} \\langle \\gamma, C \\rangle
+ \\varepsilon \\Omega(\\gamma),
```
where ``\\Omega(\\gamma) = \\sum_{i,j} \\gamma_{i,j} \\log \\gamma_{i,j}`` is the entropic
regularization term.

This function is a wrapper of the function
[`sinkhorn`](https://pythonot.github.io/all.html#ot.sinkhorn) in the Python Optimal
Transport package. Keyword arguments are listed in the documentation of the Python function.

# Examples

```jldoctest sinkhorn
julia> μ = [0.5, 0.2, 0.3];

julia> ν = [0.0, 1.0];

julia> C = [0.0 1.0; 2.0 0.0; 0.5 1.5];

julia> sinkhorn(μ, ν, C, 0.01)
3×2 Matrix{Float64}:
 0.0  0.5
 0.0  0.2
 0.0  0.3
```

It is possible to provide multiple target marginals as columns of a matrix. In this case the
optimal transport costs are returned:

```jldoctest sinkhorn
julia> ν = [0.0 0.5; 1.0 0.5];

julia> round.(sinkhorn(μ, ν, C, 0.01); sigdigits=6)
2-element Vector{Float64}:
 0.95
 0.45
```

See also: [`sinkhorn2`](@ref)
"""
function sinkhorn(μ, ν, C, ε; kwargs...)
    return pot.sinkhorn(μ, ν, PyCall.PyReverseDims(permutedims(C)), ε; kwargs...)
end

"""
    sinkhorn2(μ, ν, C, ε; kwargs...)

Compute the optimal transport cost for the entropic regularization optimal transport problem
with source and target marginals `μ` and `ν`, cost matrix `C` of size
`(length(μ), length(ν))`, and entropic regularization parameter `ε`.

The optimal transport cost is the scalar value
```math
\\inf_{\\gamma \\in \\Pi(\\mu, \\nu)} \\langle \\gamma, C \\rangle
+ \\varepsilon \\Omega(\\gamma),
```
where ``\\Omega(\\gamma) = \\sum_{i,j} \\gamma_{i,j} \\log \\gamma_{i,j}`` is the entropic
regularization term.

This function is a wrapper of the function
[`sinkhorn2`](https://pythonot.github.io/all.html#ot.sinkhorn2) in the Python Optimal
Transport package. Keyword arguments are listed in the documentation of the Python function.

# Examples

```jldoctest sinkhorn2
julia> μ = [0.5, 0.2, 0.3];

julia> ν = [0.0, 1.0];

julia> C = [0.0 1.0; 2.0 0.0; 0.5 1.5];

julia> round.(sinkhorn2(μ, ν, C, 0.01); sigdigits=6)
1-element Vector{Float64}:
 0.95
```

It is possible to provide multiple target marginals as columns of a matrix.

```jldoctest sinkhorn2
julia> ν = [0.0 0.5; 1.0 0.5];

julia> round.(sinkhorn2(μ, ν, C, 0.01); sigdigits=6)
2-element Vector{Float64}:
 0.95
 0.45
```

See also: [`sinkhorn`](@ref)
"""
function sinkhorn2(μ, ν, C, ε; kwargs...)
    return pot.sinkhorn2(μ, ν, PyCall.PyReverseDims(permutedims(C)), ε; kwargs...)
end

"""
    empirical_sinkhorn_divergence(xsource, xtarget, ε; kwargs...)

Compute the Sinkhorn divergence from empirical data, where `xsource` and `xtarget` are
arrays representing samples in the source domain and target domain, respectively, and  `ε`
is the regularization term.

This function is a wrapper of the function
[`ot.bregman.empirical_sinkhorn_divergence`](https://pythonot.github.io/gen_modules/ot.bregman.html#ot.bregman.empirical_sinkhorn_divergence)
in the Python Optimal Transport package. Keyword arguments are listed in the documentation of the Python function.

# Examples

```jldoctest
julia> xsource = [1];

julia> xtarget = [2, 3];

julia> ε = 0.01;

julia> empirical_sinkhorn_divergence(xsource, xtarget, ε) ≈
       sinkhorn2([1], [0.5, 0.5], [1 4], ε) -
       (
           sinkhorn2([1], [1], zeros(1, 1), ε) +
           sinkhorn2([0.5, 0.5], [0.5, 0.5], [0 1; 1 0], ε)
       ) / 2
true
```

See also: [`sinkhorn2`](@ref)
"""
function empirical_sinkhorn_divergence(xsource, xtarget, ε; kwargs...)
    return pot.bregman.empirical_sinkhorn_divergence(
        reshape(xsource, Val(2)), reshape(xtarget, Val(2)), ε; kwargs...
    )
end

"""
    sinkhorn_unbalanced(μ, ν, C, ε, λ; kwargs...)

Compute the optimal transport plan for the unbalanced entropic regularization optimal
transport problem with source and target marginals `μ` and `ν`, cost matrix `C` of size
`(length(μ), length(ν))`, entropic regularization parameter `ε`, and marginal relaxation
term `λ`.

The optimal transport plan `γ` is of the same size as `C` and solves
```math
\\inf_{\\gamma} \\langle \\gamma, C \\rangle
+ \\varepsilon \\Omega(\\gamma)
+ \\lambda \\mathrm{KL}(\\gamma 1, \\mu)
+ \\lambda \\mathrm{KL}(\\gamma^{\\mathsf{T}} 1, \\nu),
```
where ``\\Omega(\\gamma) = \\sum_{i,j} \\gamma_{i,j} \\log \\gamma_{i,j}`` is the entropic
regularization term and ``\\mathrm{KL}`` is the Kullback-Leibler divergence.

This function is a wrapper of the function
[`sinkhorn_unbalanced`](https://pythonot.github.io/all.html#ot.sinkhorn_unbalanced) in the
Python Optimal Transport package. Keyword arguments are listed in the documentation of the
Python function.

# Examples

```jldoctest sinkhorn_unbalanced
julia> μ = [0.5, 0.2, 0.3];

julia> ν = [0.0, 1.0];

julia> C = [0.0 1.0; 2.0 0.0; 0.5 1.5];

julia> sinkhorn_unbalanced(μ, ν, C, 0.01, 1_000)
3×2 Matrix{Float64}:
 0.0  0.499964
 0.0  0.200188
 0.0  0.29983
```

It is possible to provide multiple target marginals as columns of a matrix. In this case the
optimal transport costs are returned:

```jldoctest sinkhorn_unbalanced
julia> ν = [0.0 0.5; 1.0 0.5];

julia> round.(sinkhorn_unbalanced(μ, ν, C, 0.01, 1_000); sigdigits=6)
2-element Vector{Float64}:
 0.949709
 0.449411
```

See also: [`sinkhorn_unbalanced2`](@ref)
"""
function sinkhorn_unbalanced(μ, ν, C, ε, λ; kwargs...)
    return pot.sinkhorn_unbalanced(
        μ, ν, PyCall.PyReverseDims(permutedims(C)), ε, λ; kwargs...
    )
end

"""
    sinkhorn_unbalanced2(μ, ν, C, ε, λ; kwargs...)

Compute the optimal transport cost for the unbalanced entropic regularization optimal
transport problem with source and target marginals `μ` and `ν`, cost matrix `C` of size
`(length(μ), length(ν))`, entropic regularization parameter `ε`, and marginal relaxation
term `λ`.

The optimal transport cost is the scalar value
```math
\\inf_{\\gamma} \\langle \\gamma, C \\rangle
+ \\varepsilon \\Omega(\\gamma)
+ \\lambda \\mathrm{KL}(\\gamma 1, \\mu)
+ \\lambda \\mathrm{KL}(\\gamma^{\\mathsf{T}} 1, \\nu),
```
where ``\\Omega(\\gamma) = \\sum_{i,j} \\gamma_{i,j} \\log \\gamma_{i,j}`` is the entropic
regularization term and ``\\mathrm{KL}`` is the Kullback-Leibler divergence.

This function is a wrapper of the function
[`sinkhorn_unbalanced2`](https://pythonot.github.io/all.html#ot.sinkhorn_unbalanced2) in the
Python Optimal Transport package. Keyword arguments are listed in the documentation of the
Python function.

# Examples

```jldoctest sinkhorn_unbalanced2
julia> μ = [0.5, 0.2, 0.3];

julia> ν = [0.0, 1.0];

julia> C = [0.0 1.0; 2.0 0.0; 0.5 1.5];

julia> round.(sinkhorn_unbalanced2(μ, ν, C, 0.01, 1_000); sigdigits=6)
1-element Vector{Float64}:
 0.949709
```

It is possible to provide multiple target marginals as columns of a matrix:

```jldoctest sinkhorn_unbalanced2
julia> ν = [0.0 0.5; 1.0 0.5];

julia> round.(sinkhorn_unbalanced2(μ, ν, C, 0.01, 1_000); sigdigits=6)
2-element Vector{Float64}:
 0.949709
 0.449411
```

See also: [`sinkhorn_unbalanced`](@ref)
"""
function sinkhorn_unbalanced2(μ, ν, C, ε, λ; kwargs...)
    return pot.sinkhorn_unbalanced2(
        μ, ν, PyCall.PyReverseDims(permutedims(C)), ε, λ; kwargs...
    )
end

"""
    barycenter(A, C, ε; kwargs...)

Compute the entropically regularized Wasserstein barycenter with histograms `A`, cost matrix
`C`, and entropic regularization parameter `ε`.

The Wasserstein barycenter is a histogram and solves
```math
\\inf_{a} \\sum_{i} W_{\\varepsilon,C}(a, a_i),
```
where the histograms ``a_i`` are columns of matrix `A` and ``W_{\\varepsilon,C}(a, a_i)``
is the optimal transport cost for the entropically regularized optimal transport problem
with marginals ``a`` and ``a_i``, cost matrix ``C``, and entropic regularization parameter
``\\varepsilon``. Optionally, weights of the histograms ``a_i`` can be provided with the
keyword argument `weights`.

This function is a wrapper of the function
[`barycenter`](https://pythonot.github.io/all.html#ot.barycenter) in the
Python Optimal Transport package. Keyword arguments are listed in the documentation of the
Python function.

# Examples

```jldoctest
julia> A = rand(10, 3);

julia> A ./= sum(A; dims=1);

julia> C = rand(10, 10);

julia> isapprox(sum(barycenter(A, C, 0.01; method="sinkhorn_stabilized")), 1; atol=1e-4)
true
```
"""
barycenter(A, C, ε; kwargs...) = pot.barycenter(A, C, ε; kwargs...)

"""
    barycenter_unbalanced(A, C, ε, λ; kwargs...)

Compute the entropically regularized unbalanced Wasserstein barycenter with histograms `A`, cost matrix
`C`, entropic regularization parameter `ε` and marginal relaxation parameter `λ`.

The Wasserstein barycenter is a histogram and solves
```math
\\inf_{a} \\sum_{i} W_{\\varepsilon,C,\\lambda}(a, a_i),
```
where the histograms ``a_i`` are columns of matrix `A` and ``W_{\\varepsilon,C,\\lambda}(a, a_i)}``
is the optimal transport cost for the entropically regularized optimal transport problem
with marginals ``a`` and ``a_i``, cost matrix ``C``, entropic regularization parameter
``\\varepsilon`` and marginal relaxation parameter ``\\lambda``. Optionally, weights of the histograms ``a_i`` can be provided with the
keyword argument `weights`.

This function is a wrapper of the function
[`barycenter_unbalanced`](https://pythonot.github.io/gen_modules/ot.unbalanced.html#ot.unbalanced.barycenter_unbalanced) in the
Python Optimal Transport package. Keyword arguments are listed in the documentation of the
Python function.

# Examples

```jldoctest
julia> A = rand(10, 3);

julia> A ./= sum(A; dims=1);

julia> C = rand(10, 10);

julia> isapprox(sum(barycenter_unbalanced(A, C, 0.01, 1; method="sinkhorn_stabilized")), 1; atol=1e-4)
false

julia> isapprox(sum(barycenter_unbalanced(
           A, C, 0.01, 10_000; method="sinkhorn_stabilized", numItermax=5_000
       )), 1; atol=1e-4)
true
```

See also: [`barycenter`](@ref)
"""
function barycenter_unbalanced(A, C, ε, λ; kwargs...)
    return pot.barycenter_unbalanced(A, C, ε, λ; kwargs...)
end


text:
mutable struct SourceTerm
    λ::Float64
    μ::Float64
end

function (st::SourceTerm)(x::Float64,p::Float64)
    return st.λ*tanh(st.μ*(p-x))
end