text:
import Pluto: Pluto, Cell
import Pluto.MoreAnalysis

using Test

@testset "MoreAnalysis" begin

    file = joinpath(@__DIR__, "parallelpaths4.jl")

    newpath = tempname()
    Pluto.readwrite(file, newpath)

    notebook = Pluto.load_notebook(newpath)
    # Run pluto's analysis. This is like opening the notebook, without actually running it
    # s = Pluto.ServerSession()
    # Pluto.update_save_run!(s, notebook, notebook.cells; run)
    notebook.topology = Pluto.updated_topology(notebook.topology, notebook, notebook.cells)

    ind(c::Cell) = findfirst(isequal(c), notebook.cells)
    cell_indices(cs) = sort(ind.(cs))

    @testset "Recursive stuff" begin

        @test cell_indices(MoreAnalysis.downstream_recursive(notebook, notebook.topology, notebook.cells[1:1])) == [1,3,4,5,6,7]
        @test cell_indices(MoreAnalysis.upstream_recursive(notebook, notebook.topology, notebook.cells[1:1])) == [1]

        @test cell_indices(MoreAnalysis.downstream_recursive(notebook, notebook.topology, notebook.cells[6:6])) == [6,7]
        @test cell_indices(MoreAnalysis.upstream_recursive(notebook, notebook.topology, notebook.cells[6:6])) == [1,2,3,4,6]

        @test cell_indices(MoreAnalysis.downstream_recursive(notebook, notebook.topology, notebook.cells[7:7])) == [7]
        @test cell_indices(MoreAnalysis.upstream_recursive(notebook, notebook.topology, notebook.cells[7:7])) == 1:7

        @test cell_indices(MoreAnalysis.downstream_recursive(notebook, notebook.topology, notebook.cells[4:5])) == 4:7
        @test cell_indices(MoreAnalysis.upstream_recursive(notebook, notebook.topology, notebook.cells[4:5])) == 1:5

        @test cell_indices(MoreAnalysis.downstream_recursive(notebook, notebook.topology, notebook.cells[1:2])) == 1:7
        @test cell_indices(MoreAnalysis.upstream_recursive(notebook, notebook.topology, notebook.cells[1:7])) == 1:7

        @test MoreAnalysis.upstream_recursive(notebook, notebook.topology, Cell[]) == Set{Cell}()
        @test MoreAnalysis.downstream_recursive(notebook, notebook.topology, Cell[]) == Set{Cell}()

    end

    @testset "Bond connections" begin

        # bound_variables = (map(notebook.cells) do cell
        #     MoreAnalysis.find_bound_variables(cell.parsedcode)
        # end)

        # @show bound_variables

        connections = MoreAnalysis.bound_variable_connections_graph(notebook)
        # @show connections

        @test !isempty(connections)
        wanted_connections = Dict(
            :x         => [:y, :x],
            :y         => [:y, :x],
            :show_dogs => [:show_dogs],
            :b         => [:b],
            :c         => [:c],
            :five1     => [:five1],
            :five2     => [:five2],
            :six1      => [:six2, :six1],
            :six2      => [:six3, :six2, :six1],
            :six3      => [:six3, :six2],
            :cool1     => [:cool1, :cool2],
            :cool2     => [:cool1, :cool2],
            :world     => [:world],
            :boring    => [:boring],
        )

        transform(d) = Dict(k => sort(v) for (k, v) in d)

        @test transform(connections) == transform(wanted_connections)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 00f8bf47-47d2-4a99-86ba-d7cd0b7f268b
md""" ## First steps

Use this [template](https://github.com/sapal6/nbdev.jl-template.git) to create a new repository. This will create a repository with the directory structure required for nbdev."""

# ╔═╡ fbb6eadc-80a1-4cd9-ad27-bebcdd6e5981
md"""
!!! note

*Right now nbdev is highly opinionated about the folder structure like presence of "src" and "docs" directory.*"""

# ╔═╡ 170cd996-3652-4157-bcd5-cf782ce74896
md"
!!! tip

*Refer to this [link](https://docs.github.com/en/repositories/creating-and-managing-repositories/creating-a-repository-from-a-template) to understand how to use a template repository.*"

# ╔═╡ b0c1fe8f-e796-4835-95f8-17990c7421a7
md"## First notebook"

# ╔═╡ e0694cd4-7d07-494d-a010-f705d1cbeb2c
md"""As a first steps clone your repository to your local machine.

Once you have cloned the repository, start Pluto the usual way as you use to do and create a notebook inside the "nbs" directory.

Each Pluto notebook maintains it's own environment. You can directly import a package and Pluto will install it from the registry. However, in our case Nbdev.jl is not yet registered. So, you can install Nbdev in the following way"""

# ╔═╡ 15a5c4fd-301e-4f7c-9ad3-28ae2d214224
md"""In the first cell of your notebook type in the following code 👇"""

# ╔═╡ 5c86332a-db85-4326-9f17-215a8a0188c8
md"""
```
using Nbdev
```
"""

# ╔═╡ e89e5c16-845e-4d0a-a346-80586c100be0
md"This will install Nbdev in our notebook environment."

# ╔═╡ 2ac2c281-ad4a-4364-8f62-5c8202d4e9c1
md"!!! note

*If you don't want to install Nbdev as part of every notebook that you create in your project, then you can follow the other options to install packages for your notebook projects by following the instructions in the [Pluto wiki](https://github.com/fonsp/Pluto.jl/wiki/%F0%9F%8E%81-Package-management).*"

# ╔═╡ 9258500a-0a33-41cd-b32a-e09e820f84ee
md"## Start building your project"

# ╔═╡ dedef829-e43c-41c9-8238-b9d631dccb35
md"There are some anchors which tell Nbdev about which portion of your notebook should be exported as code, which portion should be part of documentation and which portion should not be touched."

# ╔═╡ 99d43364-5ac6-43bb-a70c-20a79a463fdf
md"### #export"

# ╔═╡ fb11c61e-2b26-4b2a-9a34-835163ac7f87
md"when you put the comment #export at the begining of a cell, it tells nbdev to export that cell as a code into the source file."

# ╔═╡ c9f87989-522c-4adc-9a50-af95e3fa9cfa
md"
```
#export
function somefunc(x)
	x+1
end
```"

# ╔═╡ 0d4fb029-bfc8-4208-a6f4-bd5ea981ff82
md"Nbdev will now put the above cell as part of the source file (with the same name as the current notebook) inside the src directory."

# ╔═╡ ac65836e-6f4f-4212-b7f4-a648de9c761b
md"### #hide"

# ╔═╡ 789ebdb8-0017-4662-9120-a780584f8368
md"Suppose that you have included soem cells as par of experment while developing your package but you want that the cell should neither be exported as a code nor that should be exported as part of documentation.

In such a case use the `#hide` anchor at the begining of the cell."

# ╔═╡ 809102dd-b91a-458f-b5c7-2a49c7536d3b
md"
```
#hide
myexperiment = 1+1
```"

# ╔═╡ f8391a08-edf6-499c-8e8d-8f731dad8fd4
md"such a cell won't be exported."

# ╔═╡ 4d12cf7f-54d7-4a06-9749-a97246e7564e
md"Any other cell which neither has the `#hide` nor the `#export` anchor, is exported as an example block. For such a block the cell is executed and the output is captured.

The output and the code of an example cell appears in the final document."

# ╔═╡ c40eeb9c-dbda-4dc0-9fc3-e91be447ac25
md"### #noop"

# ╔═╡ 929d68db-fac5-4507-8316-a0dcf71df403
md"Nbdev captures the output for cells which are not marked with `#export` by default. but there might be times when you have an example for which you just want to show the code and not the output. Such cases can be handled by the `#noop` anchor. Just start such cells with `#noop` and nbdev will not capture output for such cells."

# ╔═╡ fac3e16d-76ff-476e-a2c4-ff34f7c3f1cc
md"""
```
begin
#noop
"
some experiment
"
function myexperiment
    1+1
end
myexperiment
end
```"""

# ╔═╡ 9fcf1320-0ddf-471c-a584-092586122347
md"In the above case even though Pluto displays the docstrign as an output, that output won't be displayed in your document.

As another example let's consider something like the code below"

# ╔═╡ a3ea3f66-a66c-41d1-8c3c-c1cd3a879871
md"
```
#noop
add = 1+1
add
````
"

# ╔═╡ 77b75b82-c684-4377-882a-ff0a4283c7f5
md"The output for such code cell won't be captured in the documentation. Only the code will be displayed once you build the document."

# ╔═╡ a7714f66-3eff-481a-af7f-ae26c2b8e8e2
md"## Docstrings"

# ╔═╡ 926b30a2-1a86-4952-a991-d75b5ce6e545
md"Define a docstring for functions, structs etc. and use the special function `showDoc` to document that function or struct"

# ╔═╡ 5bd889ce-9522-462b-a843-508061926306
"""
> This is a function docstring
"""
function somefunc()
	1+1
end

# ╔═╡ 8c9f1b35-8f88-44ba-8a37-57b92b42980b
md"`showdoc(somefunc)` would show the function like this 👇"

# ╔═╡ 52e3492a-c27d-476f-9010-7d009fd0eeaf
md"`showdoc(somefunc)`"

# ╔═╡ 306239c9-7fd6-442e-8b97-db674398b307
md"## Exporting code"

# ╔═╡ a8955cbd-2530-405f-b876-18bb35e0b008
md"""Once you are done with your notebook and you can run `notebook2script()`. This reads the current notebook in the "nbs" folder, scrubs out all the cells which are do not start with `#export` and creates a module file in the src folder.

The module file in the src folder will have the same name as the current file with the first letter in uppercase.

The exported code file automatically has all the code inside the module with the same name as the file.

For example, the exported source code file would have the following structure 👇

```
module Example

function somefunc(x)
    x+1
end

end
```"""

# ╔═╡ 84f7040c-0ca9-4d61-abfe-789fc446ebb2
md"""
Once you are happy with the code you can export the source code by calling

`notebook2script(nbsdir, srcdir)`

Here "nbs_dir" is the directory where your notebooks are stored and the "src_dir" is the directory where you want your source code files to be exported. This is usually the `src` directory."""

# ╔═╡ 7b1fdf81-870a-4977-8a9d-c867ee2861ff
md"## Tests"

# ╔═╡ 6ace0489-21ae-488f-8ea5-d3087651394c
md"""Nbdev comes with support for [PlutoTests(in alpha)](https://github.com/JuliaPluto/PlutoTest.jl) out of the box.

What this means is that you can write reactive tests within the same notebook where you are writing your code. when you run that particular notebook, the affected tests are executed. So, you can write tests as you write your code.

In addition to reactivity you can also take advantage of other goodies in PlutoTests like "time travel" to see the different stages of your code (follow the link to PlutoTests for further details)."""

# ╔═╡ 50ef31b0-b5f2-4a8e-b123-1681645a031f
md"!!! note
Currently PlutoTests is in **alpha** so use it with caution. However, you can even use the native test suite in Julia to write your tests but won't get features like time travel while using the native test suite.

Automatic execution of tests using github actions or Travis CI is planned for future versions of PlutoTest"

# ╔═╡ bc79b045-8969-4014-b7ac-01282d969c8c
md"## Building documentation"

# ╔═╡ eb26aca4-4b12-435d-941f-edb18bda0c3c
md"""
Do not forget to edit the "index" notebook available in the "nbs" folder. This notebook would be used as your homepage and the contents of this notebook are also used in the README file."""

# ╔═╡ 547f7e70-71bf-44a6-ba62-0e35edf688fc
md"!!! note

*You will get the index notebook in the nbs directory when you create a new project using the nbdev template.*
"

# ╔═╡ b1e0024c-aa32-46a3-baea-534b1187fb7a
md"""To build documentation, goto the project root and launch the Julia prompt and then import Nbdev with `import Nbdev` and then call the following function `Nbdev.build(<your notebook directory)`. Your code at the Julia prompt should look like this

```
import Nbdev
Nbdev.build("nbs")
```

This will run all your notebooks in the notebook directory (/nbs usually) and pick up the cells which are not marked as `#export` or #hide and then create markdown files in the `/docs` directory."""

# ╔═╡ fa7161ef-219f-47b5-b4a0-b18e0d3b5b54
md"## Creating the Table of contents"

# ╔═╡ 54817dac-18b6-40cc-9ef6-b51370164839
md"""The default configuration in the "mkdocs.yml" file (available to you when you create your repository from the template) contains all the necessary default configuration required to create you project page.

In order to test how your website looks like, go to the project root and type `mkdocs serve` . This would serve the project documentation to view locally (localhost).

The table of contents in the default case would be created by mkdocs on the basis of the documents available under the "docs" directory.

If you want to give your own spin to the table of contents then add a config similar to the following in the mkdocs.yml file-->

```
nav:
    - Home: index.md
    - Docs:
      - Nav1: nav1.md
      - Nav2: nav2.md
```
"""

# ╔═╡ 2b4992b3-f1a4-4587-b0d6-cb30efb1dfa1
md"""
!!! note
*Make sure that the "Home" navigation should always be the index.md file*"""

# ╔═╡ 394597e3-42f8-4bd3-8f10-206b5ea5f225
md"For more details about setting up navigation refer [this](https://www.mkdocs.org/user-guide/writing-your-docs/) and [this](https://squidfunk.github.io/mkdocs-material/setup/extensions/python-markdown/?h=table+con+of#table-of-contents) "

# ╔═╡ 094e214e-6681-4de8-b9a4-e73c4b23ff1c
md"## Other required edits"

# ╔═╡ 23ac38d5-c44a-4b8b-86c4-0bc3833add14
md"""In the mkdocs.yml file, edit the "site_name" config to the name that you want to give to your site. The "repo_url" and "repo_name" should be your git repository url and the name of your git repository respectively."""

# ╔═╡ 33eaee48-fc83-4c71-949d-1c9b9a2b32aa
md"## Optional edits"

# ╔═╡ 952bc125-e3b6-474c-90c4-786fdf300dcc
md"There are loads f other optional edits which you can make to the config file in order to customize your project website. Refer the documentation [here](https://squidfunk.github.io/mkdocs-material/getting-started/) to know about the other options."

# ╔═╡ 7429a0be-4f40-4d36-89cd-df56377c4502
md"## Publishing your project website"

# ╔═╡ 8fa0648a-599b-4ed5-870f-6b3d74e25fcd
md"""The first step is to checkin all your project contents to your repoitory. As soon as you push your changes to the remote branch, the github action (placed under "/.github/workflows" creates a new "gh-pages" branch under your project repository, builds the project webpages and pushes the created pages to the ""gh-pages" branch.

If you go to the "Actions" tab in your github repository then you can view the Action under the "All workflows" section. Here, a green tick mark before the action name suggests that all went fine.

Now, go to the Settings tab in your repository, click on Pages menu item from the left menu and in the "Source" drop-down select "gh-pages" as the source and root as the source folder, Click on save.

A detailed step by step description about configuring github pages can be read [here](https://docs.github.com/en/pages/getting-started-with-github-pages/configuring-a-publishing-source-for-your-github-pages-site).

You project website would be published at the following url --> "<username>.github.io/<project reporsitory name>"
"""

# ╔═╡ Cell order:
# ╠═00f8bf47-47d2-4a99-86ba-d7cd0b7f268b
# ╠═fbb6eadc-80a1-4cd9-ad27-bebcdd6e5981
# ╠═170cd996-3652-4157-bcd5-cf782ce74896
# ╠═b0c1fe8f-e796-4835-95f8-17990c7421a7
# ╠═e0694cd4-7d07-494d-a010-f705d1cbeb2c
# ╠═15a5c4fd-301e-4f7c-9ad3-28ae2d214224
# ╠═5c86332a-db85-4326-9f17-215a8a0188c8
# ╠═e89e5c16-845e-4d0a-a346-80586c100be0
# ╠═2ac2c281-ad4a-4364-8f62-5c8202d4e9c1
# ╠═9258500a-0a33-41cd-b32a-e09e820f84ee
# ╠═dedef829-e43c-41c9-8238-b9d631dccb35
# ╠═99d43364-5ac6-43bb-a70c-20a79a463fdf
# ╠═fb11c61e-2b26-4b2a-9a34-835163ac7f87
# ╠═c9f87989-522c-4adc-9a50-af95e3fa9cfa
# ╠═0d4fb029-bfc8-4208-a6f4-bd5ea981ff82
# ╠═ac65836e-6f4f-4212-b7f4-a648de9c761b
# ╠═789ebdb8-0017-4662-9120-a780584f8368
# ╠═809102dd-b91a-458f-b5c7-2a49c7536d3b
# ╠═f8391a08-edf6-499c-8e8d-8f731dad8fd4
# ╠═4d12cf7f-54d7-4a06-9749-a97246e7564e
# ╠═c40eeb9c-dbda-4dc0-9fc3-e91be447ac25
# ╠═929d68db-fac5-4507-8316-a0dcf71df403
# ╠═fac3e16d-76ff-476e-a2c4-ff34f7c3f1cc
# ╠═9fcf1320-0ddf-471c-a584-092586122347
# ╠═a3ea3f66-a66c-41d1-8c3c-c1cd3a879871
# ╠═77b75b82-c684-4377-882a-ff0a4283c7f5
# ╠═a7714f66-3eff-481a-af7f-ae26c2b8e8e2
# ╠═926b30a2-1a86-4952-a991-d75b5ce6e545
# ╠═5bd889ce-9522-462b-a843-508061926306
# ╠═8c9f1b35-8f88-44ba-8a37-57b92b42980b
# ╠═52e3492a-c27d-476f-9010-7d009fd0eeaf
# ╠═306239c9-7fd6-442e-8b97-db674398b307
# ╠═a8955cbd-2530-405f-b876-18bb35e0b008
# ╠═84f7040c-0ca9-4d61-abfe-789fc446ebb2
# ╠═7b1fdf81-870a-4977-8a9d-c867ee2861ff
# ╠═6ace0489-21ae-488f-8ea5-d3087651394c
# ╠═50ef31b0-b5f2-4a8e-b123-1681645a031f
# ╠═bc79b045-8969-4014-b7ac-01282d969c8c
# ╠═eb26aca4-4b12-435d-941f-edb18bda0c3c
# ╠═547f7e70-71bf-44a6-ba62-0e35edf688fc
# ╠═b1e0024c-aa32-46a3-baea-534b1187fb7a
# ╠═fa7161ef-219f-47b5-b4a0-b18e0d3b5b54
# ╠═54817dac-18b6-40cc-9ef6-b51370164839
# ╠═2b4992b3-f1a4-4587-b0d6-cb30efb1dfa1
# ╠═394597e3-42f8-4bd3-8f10-206b5ea5f225
# ╠═094e214e-6681-4de8-b9a4-e73c4b23ff1c
# ╠═23ac38d5-c44a-4b8b-86c4-0bc3833add14
# ╠═33eaee48-fc83-4c71-949d-1c9b9a2b32aa
# ╠═952bc125-e3b6-474c-90c4-786fdf300dcc
# ╠═7429a0be-4f40-4d36-89cd-df56377c4502
# ╠═8fa0648a-599b-4ed5-870f-6b3d74e25fcd

ext: jl


text:
using Statistics
using CUDA.CUDNN:
    cudnnReduceTensor,
    cudnnReduceTensor!,
    cudnnGetReductionIndicesSize,
    cudnnGetReductionWorkspaceSize,
    cudnnReduceTensorDescriptor,
        cudnnReduceTensorDescriptor_t,
        cudnnCreateReduceTensorDescriptor,
        cudnnSetReduceTensorDescriptor,
        cudnnGetReduceTensorDescriptor,
        cudnnDestroyReduceTensorDescriptor,
    cudnnReduceTensorOp_t,
        CUDNN_REDUCE_TENSOR_ADD,          # 0,
        CUDNN_REDUCE_TENSOR_MUL,          # 1,
        CUDNN_REDUCE_TENSOR_MIN,          # 2,
        CUDNN_REDUCE_TENSOR_MAX,          # 3,
        CUDNN_REDUCE_TENSOR_AMAX,         # 4,
        CUDNN_REDUCE_TENSOR_AVG,          # 5,
        CUDNN_REDUCE_TENSOR_NORM1,        # 6,
        CUDNN_REDUCE_TENSOR_NORM2,        # 7,
        CUDNN_REDUCE_TENSOR_MUL_NO_ZEROS, # 8,
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN, # 0
        CUDNN_PROPAGATE_NAN,     # 1
    cudnnReduceTensorIndices,
    cudnnReduceTensorIndices_t,
        CUDNN_REDUCE_TENSOR_NO_INDICES,        # 0,
        CUDNN_REDUCE_TENSOR_FLATTENED_INDICES, # 1,
    cudnnIndicesType,
    cudnnIndicesType_t,
        CUDNN_32BIT_INDICES, # 0,
        CUDNN_64BIT_INDICES, # 1,
        CUDNN_16BIT_INDICES, # 2,
        CUDNN_8BIT_INDICES,  # 3,
    cudnnDataType

@testset "cudnn/reduce" begin
    @test cudnnReduceTensorDescriptor(C_NULL) isa cudnnReduceTensorDescriptor
    @test Base.unsafe_convert(Ptr, cudnnReduceTensorDescriptor(C_NULL)) isa Ptr
    @test cudnnReduceTensorDescriptor(CUDNN_REDUCE_TENSOR_ADD,cudnnDataType(Float32),CUDNN_NOT_PROPAGATE_NAN,CUDNN_REDUCE_TENSOR_NO_INDICES,CUDNN_32BIT_INDICES) isa cudnnReduceTensorDescriptor

    (ax,ay) = randn(Float32,10,10), randn(Float32,10,1)
    (cx,cy) = CuArray.((ax,ay))

    function reducetensortest(;
        op::cudnnReduceTensorOp_t = CUDNN_REDUCE_TENSOR_ADD,
        compType::DataType = (eltype(ax) <: Float64 ? Float64 : Float32),
        nanOpt::cudnnNanPropagation_t = CUDNN_NOT_PROPAGATE_NAN,
        indices::Union{Vector{<:Unsigned},Nothing} = nothing,
        d::cudnnReduceTensorDescriptor = cudnnReduceTensorDescriptor(op, cudnnDataType(compType), nanOpt, cudnnReduceTensorIndices(op, indices), cudnnIndicesType(indices)),
        alpha::Real = 1,
        beta::Real = 0,
    )
        f0 = (op === CUDNN_REDUCE_TENSOR_ADD          ? sum(ax, dims=2) :
              op === CUDNN_REDUCE_TENSOR_MUL          ? prod(ax, dims=2) :
              op === CUDNN_REDUCE_TENSOR_MIN          ? minimum(ax, dims=2) :
              op === CUDNN_REDUCE_TENSOR_MAX          ? maximum(ax, dims=2) :
              op === CUDNN_REDUCE_TENSOR_AMAX         ? maximum(abs, ax, dims=2) :
              op === CUDNN_REDUCE_TENSOR_AVG          ? mean(ax, dims=2) :
              op === CUDNN_REDUCE_TENSOR_NORM1        ? sum(abs, ax, dims=2) :
              op === CUDNN_REDUCE_TENSOR_NORM2        ? sqrt.(sum(abs2, ax, dims=2)) :
              op === CUDNN_REDUCE_TENSOR_MUL_NO_ZEROS ? (ax1=copy(ax);ax1[ax.==0].=1;prod(ax1,dims=2)) :
              error("Unknown reducetensor"))
        f1 = alpha * f0
        f2 = f1 + beta * ay
        dims = size(ay)
        @test f1 ≈ cudnnReduceTensor(cx; dims, op, compType, nanOpt, indices, alpha) |> Array
        @test f1 ≈ cudnnReduceTensor(cx, d; dims, indices, alpha) |> Array
        @test f2 ≈ cudnnReduceTensor!(copy(cy), cx; op, compType, nanOpt, indices, alpha, beta) |> Array
        @test f2 ≈ cudnnReduceTensor!(copy(cy), cx, d; indices, alpha, beta) |> Array
    end

    reducetensortest()
    reducetensortest(op = CUDNN_REDUCE_TENSOR_MUL)
    reducetensortest(op = CUDNN_REDUCE_TENSOR_MIN)
    reducetensortest(op = CUDNN_REDUCE_TENSOR_MAX)
    reducetensortest(op = CUDNN_REDUCE_TENSOR_AMAX)
    reducetensortest(op = CUDNN_REDUCE_TENSOR_AVG)
    reducetensortest(op = CUDNN_REDUCE_TENSOR_NORM1)
    reducetensortest(op = CUDNN_REDUCE_TENSOR_NORM2)
    reducetensortest(op = CUDNN_REDUCE_TENSOR_MUL_NO_ZEROS)
    reducetensortest(nanOpt = CUDNN_PROPAGATE_NAN)
    reducetensortest(alpha = 2)
    reducetensortest(beta = 2)
end


text:
using JuMP, MathProgBase, CPLEX, Random, LinearAlgebra, JLD, Mosek, LaTeXStrings, DataFrames, Test, Suppressor, DelimitedFiles, CSV, StatsBase

include("types.jl")
include("miop_formulation.jl")
include("inner_problem.jl")
include("hillclimb.jl")
include("socp_relaxation.jl")

include("kelley_primal.jl")
include("raw.jl")


text:
using RomanNumerals
using Test
using Random

const MT = MersenneTwister

using RomanNumerals: InvalidRomanNumeral

@testset "Construction" begin
    @test RomanNumeral("I") == RomanNumeral(1)
    @test RomanNumeral("V") == RomanNumeral(5)
    @test RomanNumeral("X") == RomanNumeral(10)
    @test RomanNumeral("L") == RomanNumeral(50)
    @test RomanNumeral("C") == RomanNumeral(100)
    @test RomanNumeral("D") == RomanNumeral(500)
    @test RomanNumeral("M") == RomanNumeral(1000)

    @test RomanNumeral("II") == RomanNumeral(2)
    @test RomanNumeral("VII") == RomanNumeral(7)
    @test RomanNumeral("XXX") == RomanNumeral(30)
    @test RomanNumeral("LXIV") == RomanNumeral(64)
    @test RomanNumeral("CCL") == RomanNumeral(250)
    @test RomanNumeral("CD") == RomanNumeral(400)
    @test RomanNumeral("M") == RomanNumeral(1000)

    @test rn"I" == RomanNumeral(1)
    @test rn"V" == RomanNumeral(5)
    @test rn"X" == RomanNumeral(10)
    @test rn"L" == RomanNumeral(50)
    @test rn"C" == RomanNumeral(100)
    @test rn"D" == RomanNumeral(500)
    @test rn"M" == RomanNumeral(1000)

    @test rn"MMMCDLIII" == RomanNumeral(3453)
    @test rn"MCDXLIX" == RomanNumeral(1449)
    @test rn"MMMCMXCIV" == RomanNumeral(3994)
    @test rn"CCCCXXXXIIII" == rn"CDXLIV" == RomanNumeral(444)

    @test_throws InvalidRomanNumeral RomanNumeral("Y")
    @test_throws InvalidRomanNumeral RomanNumeral("XY")
end

@testset "Consts" begin
    @test I == RomanNumeral(1)
    @test V == RomanNumeral(5)
    @test X == RomanNumeral(10)
    @test L == RomanNumeral(50)
    @test C == RomanNumeral(100)
    @test D == RomanNumeral(500)
    @test M == RomanNumeral(1000)
end

@testset "Conversion" begin
    @testset "Integer" begin
        @test Int(rn"I") == 1
        @test Int8(rn"V") == 5
        @test Int16(rn"X") == 10
        @test Int32(rn"L") == 50
        @test Int64(rn"C") == 100
        @test Int16(rn"D") == 500
        @test Int32(rn"M") == 1000
    end
    @testset "String" begin
        @test string(rn"I") == "I"
        @test string(rn"V") == "V"
        @test string(rn"X") == "X"
        @test string(rn"L") == "L"
        @test string(rn"C") == "C"
        @test string(rn"D") == "D"
        @test string(rn"M") == "M"
    end
end

@testset "Arithmetic" begin
    @testset "Multiplication/Division" begin
        @test rn"II" == 2I == RomanNumeral(2)
        @test rn"III" == 3I == RomanNumeral(3)
        @test rn"IV" == 4I == RomanNumeral(4)
        @test rn"IX" == 9I == RomanNumeral(9)
        @test rn"XX" == 2X == RomanNumeral(20)
        @test rn"CCC" == 3C == RomanNumeral(300)
        @test rn"MMMM" == 4M == RomanNumeral(4000)

        @test I ÷ I == RomanNumeral(1)
        @test X ÷ V == RomanNumeral(2)
        @test L ÷ X == V
        @test C ÷ L == 2I
        @test C ÷ V == 2I * X
        @test M ÷ L == X + X
    end
    @testset "Addition/Subtraction" begin
        @test I + I == RomanNumeral(2)
        @test I + I + I == RomanNumeral(3)
        @test V - I == RomanNumeral(4)
        @test X - I == RomanNumeral(9)
        @test X + X == RomanNumeral(20)
        @test D - 2C == RomanNumeral(300)
        @test M + M + M + M == RomanNumeral(4000)
        @test V - 2 == RomanNumeral(3)
        @test 5 - rn"II" == RomanNumeral(3)
    end
end

@testset "Random" begin
    @test rand(MT(1), RomanNumeral) == rn"MCXCIX"
    @test rand(MT(1), RomanNumeral, 3) == RomanNumeral.(["MCXCIX", "MCLV", "CCCLX"])
    @test rand(MT(2), RomanNumeral, 3) ==
        RomanNumeral.(["DCCCXCV", "DCCXXXVI", "MDCCXXXVIII"])
    @test rand(RomanNumeral) isa RomanNumeral{Int}
    for T = [Int8, Int16, Int32, Int64, Int128]
        @test rand(RomanNumeral{T}) isa RomanNumeral{T}
    end
end


text:
module TestBasic
using Telegram
using Telegram.API
using Test

@testset "basic api call" begin
    tg = TelegramClient("xxx"; chat_id = "yyy")
    useglobally!(tg)
    @test Telegram.DEFAULT_OPTS.client === tg
end

end # module


text:
using Test
using Bijectors
using ForwardDiff: derivative, jacobian
using LinearAlgebra: logabsdet, I
using Random

Random.seed!(123)

# logabsdet doesn't handle scalars.
_logabsdet(x::AbstractArray) = logabsdet(x)[1]
_logabsdet(x::Real) = log(abs(x))

# Generate a (vector / matrix of) random number(s).
_rand_real(::Real) = randn()
_rand_real(x) = randn(size(x))
_rand_real(x, e) = (y = randn(size(x)); y[end] = e; y)

# Standard tests for all distributions involving a single-sample.
function single_sample_tests(dist, jacobian)
    ϵ = eps(Float64)
    # Do the regular single-sample tests.
    single_sample_tests(dist)

    # Check that the implementation of the logpdf agrees with the AD version.
    x = rand(dist)
    if dist isa SimplexDistribution
        logpdf_ad = logpdf(dist, x .+ ϵ) - _logabsdet(jacobian(x->link(dist, x, Val{false}), x))
    else
        logpdf_ad = logpdf(dist, x) - _logabsdet(jacobian(x->link(dist, x), x))
    end
    @test logpdf_ad ≈ logpdf_with_trans(dist, x, true)
end

# Standard tests for all distributions involving a single-sample. Doesn't check that the
# logpdf implementation is consistent with the link function for technical reasons.
function single_sample_tests(dist)
    ϵ = eps(Float64)

    # Check that invlink is inverse of link.
    x = rand(dist)
    @test invlink(dist, link(dist, copy(x))) ≈ x atol=1e-9

    # Check that link is inverse of invlink. Hopefully this just holds given the above...
    y = link(dist, x)
    if dist isa Dirichlet
        # `logit` and `logistic` are not perfect inverses. This leads to a diversion.
        # Example:
        # julia> logit(logistic(0.9999999999999998))
        #    1.0
        # julia> logistic(logit(0.9999999999999998))
        # 0.9999999999999998
        @test link(dist, invlink(dist, copy(y))) ≈ y atol=0.5
    else
        @test link(dist, invlink(dist, copy(y))) ≈ y atol=1e-9
    end
    if dist isa SimplexDistribution
        # This should probably be exact.
        @test logpdf(dist, x .+ ϵ) == logpdf_with_trans(dist, x, false)
        # Check that invlink maps back to the apppropriate constrained domain.
        @test all(isfinite, logpdf.(Ref(dist), [invlink(dist, _rand_real(x, 0)) .+ ϵ for _ in 1:100]))
    else
        # This should probably be exact.
        @test logpdf(dist, x) == logpdf_with_trans(dist, x, false)
        @test all(isfinite, logpdf.(Ref(dist), [invlink(dist, _rand_real(x)) for _ in 1:100]))
    end
    # This is a quirk of the current implementation, of which it would be nice to be rid.
    @test typeof(x) == typeof(y)
end

# Standard tests for all distributions involving multiple samples. xs should be whatever
# the appropriate repeated version of x is for the distribution in question. ie. for
# univariate distributions, just a vector of identical values. For vector-valued
# distributions, a matrix whose columns are identical.
function multi_sample_tests(dist, x, xs, N)
    ys = link(dist, copy(xs))
    @test invlink(dist, link(dist, copy(xs))) ≈ xs atol=1e-9
    @test link(dist, invlink(dist, copy(ys))) ≈ ys atol=1e-9
    @test logpdf_with_trans(dist, xs, true) == fill(logpdf_with_trans(dist, x, true), N)
    @test logpdf_with_trans(dist, xs, false) == fill(logpdf_with_trans(dist, x, false), N)

    # This is a quirk of the current implementation, of which it would be nice to be rid.
    @test typeof(xs) == typeof(ys)
end

# Scalar tests
@testset "scalar" begin
let
    # Tests with scalar-valued distributions.
    uni_dists = [
        Arcsine(2, 4),
        Beta(2,2),
        BetaPrime(),
        Biweight(),
        Cauchy(),
        Chi(3),
        Chisq(2),
        Cosine(),
        Epanechnikov(),
        Erlang(),
        Exponential(),
        FDist(1, 1),
        Frechet(),
        Gamma(),
        InverseGamma(),
        InverseGaussian(),
        # Kolmogorov(),
        Laplace(),
        Levy(),
        Logistic(),
        LogNormal(1.0, 2.5),
        Normal(0.1, 2.5),
        Pareto(),
        Rayleigh(1.0),
        TDist(2),
        TruncatedNormal(0, 1, -Inf, 2),
    ]
    for dist in uni_dists

        single_sample_tests(dist, derivative)

        # specialised multi-sample tests.
        N = 10
        x = rand(dist)
        xs = fill(x, N)
        multi_sample_tests(dist, x, xs, N)
    end
end
end

# Tests with vector-valued distributions.
@testset "vector" begin
let ϵ = eps(Float64)
    vector_dists = [
        Dirichlet(2, 3),
        Dirichlet([1000 * one(Float64), eps(Float64)]),
        Dirichlet([eps(Float64), 1000 * one(Float64)]),
        MvNormal(randn(10), exp.(randn(10))),
        MvLogNormal(MvNormal(randn(10), exp.(randn(10)))),
        Dirichlet([1000 * one(Float64), eps(Float64)]),
        Dirichlet([eps(Float64), 1000 * one(Float64)]),
    ]
    for dist in vector_dists

        if dist isa Dirichlet
            single_sample_tests(dist)

            # This should fail at the minute. Not sure what the correct way to test this is.
            x = rand(dist)
            logpdf_turing = logpdf_with_trans(dist, x, true)
            J = jacobian(x->link(dist, x, Val{false}), x)
            @test logpdf(dist, x .+ ϵ) - _logabsdet(J) ≈ logpdf_turing

            # Issue #12
            stepsize = 1e10
            dim = length(dist)
            x = [logpdf_with_trans(dist, invlink(dist, link(dist, rand(dist)) .+ randn(dim) .* stepsize), true) for _ in 1:1_000]
            @test !any(isinf, x) && !any(isnan, x)
        else
            single_sample_tests(dist, jacobian)
        end

        # Multi-sample tests. Columns are observations due to Distributions.jl conventions.
        N = 10
        x = rand(dist)
        xs = repeat(x, 1, N)
        multi_sample_tests(dist, x, xs, N)
    end
end
end

# Tests with matrix-valued distributions.
@testset "matrix" begin
let
    matrix_dists = [
        Wishart(7, [1 0.5; 0.5 1]),
        InverseWishart(2, [1 0.5; 0.5 1]),
    ]
    for dist in matrix_dists

        single_sample_tests(dist)

        x = rand(dist); x = x + x' + 2I
        lowerinds = [LinearIndices(size(x))[I] for I in CartesianIndices(size(x)) if I[1] >= I[2]]
        upperinds = [LinearIndices(size(x))[I] for I in CartesianIndices(size(x)) if I[2] >= I[1]]
        logpdf_turing = logpdf_with_trans(dist, x, true)
        J = jacobian(x->link(dist, x), x)
        J = J[lowerinds, upperinds]
        @test logpdf(dist, x) - _logabsdet(J) ≈ logpdf_turing

        # Multi-sample tests comprising vectors of matrices.
        N = 10
        x = rand(dist)
        xs = [x for _ in 1:N]
        multi_sample_tests(dist, x, xs, N)
    end
end
end

################################## Miscelaneous old tests ##################################

# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), exp.([-1000., -1000., -1000.]), true)
# NaN
# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), [-1000., -1000., -1000.], true, true)
# -1999.30685281944
#
# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), exp.([-1., -2., -3.]), true)
# -3.006450206744678
# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), [-1., -2., -3.], true, true)
# -3.006450206744678
d  = Dirichlet([1., 1., 1.])
r  = [-1000., -1000., 0.0]
r2 = [-1., -2., 0.0]

# test link
#link(d, r)

# test invlink
@test invlink(d, r) ≈ [0., 0., 1.] atol=1e-9

# test logpdf_with_trans
#@test logpdf_with_trans(d, invlink(d, r), true) -1999.30685281944 1e-9 ≈ # atol=NaN
@test logpdf_with_trans(d, invlink(d, r2), true) ≈ -3.760398892580863 atol=1e-9


text:
abstract type Dynamics end
abstract type Population end
abstract type CaseModelFitResult end


text:
# Training a normalizing flow for 2D Rosenbrock example
# Authors: Ali Siahkoohi, alisk@gatech.edu
#          Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
#          Philipp Witte, pwitte3@gatech.edu
# Date: September 2020
"""
Example:
julia train_supervised_hint.jl
"""

using DrWatson
@quickactivate :FastApproximateInference

using InvertibleNetworks
using Random
using Statistics
using ArgParse
using BSON
using Printf
using Flux: gpu, cpu
using Flux.Data: DataLoader
using Flux.Optimise: Optimiser, ExpDecay, update!, ADAM

# Random seed
Random.seed!(19)

s = ArgParseSettings()
@add_arg_table! s begin
    "--max_epoch"
        help = "Maximum # of epochs"
        arg_type = Int
        default = 25
    "--lr"
        help = "starting learning rate"
        arg_type = Float32
        default = 1f-3
    "--lr_step"
        help = "lr scheduler applied at evey lr_step epoch"
        arg_type = Int
        default = 1
    "--batchsize"
        help = "batch size"
        arg_type = Int
        default = 64
    "--n_hidden"
        help = "# hidden channels"
        arg_type = Int
        default = 64
    "--depth"
        help = "depth of the network"
        arg_type = Int
        default = 8
    "--sigma"
        help = "noise std"
        arg_type = Float32
        default = 4f-1
    "--sim_name"
        help = "simulation name"
        arg_type = String
        default = "supervised-hint-rosenbrock"
end
parsed_args = parse_args(s)

max_epoch = parsed_args["max_epoch"]
lr = parsed_args["lr"]
lr_step = parsed_args["lr_step"]
batchsize = parsed_args["batchsize"]
n_hidden = parsed_args["n_hidden"]
depth = parsed_args["depth"]
sigma = parsed_args["sigma"]
sim_name = parsed_args["sim_name"]

# Define network
nx = 1
ny = 1
n_in = 2

CH = NetworkConditionalHINT(n_in, n_hidden, depth) |> gpu

# Training data
ntrain = 5000
X_train = rand(RosenbrockDistribution(0f0, 5f-1, ones(Float32, 2, 1)), ntrain)
Y_train = X_train + sigma*randn(Float32, nx, ny, n_in, ntrain)
X_train = X_train |> gpu
Y_train = Y_train |> gpu

# Training
# Batch extractor
num_batches = cld(ntrain, batchsize)
train_loader = DataLoader((X_train, Y_train), batchsize=batchsize, shuffle=true)

# Optimizer
opt = Optimiser(ExpDecay(lr, .9f0, num_batches*lr_step, 0f0), ADAM(lr))

# Training log keeper
fval = zeros(Float32, num_batches*max_epoch)

for j=1:max_epoch
    for (itr, (X, Y)) in enumerate(train_loader)

        fval[(j-1)*num_batches + itr] = loss_supervised(CH, X, Y)[1]

        @printf(
            "Epoch: [%d/%d] | Itr: [%d/%d] | Average negative log-likelihood: %4.2f\r",
            j, max_epoch, itr, num_batches, fval[(j-1)*num_batches + itr]
        )

        # Update params
        for p in get_params(CH)
            update!(opt, p.data, p.grad)
        end
        clear_grad!(CH)
    end
end

# Saving parameters and logs
CH = CH |> cpu
Params = get_params(CH)
save_dict = @strdict max_epoch lr lr_step batchsize n_hidden depth sigma sim_name Params fval
@tagsave(datadir(sim_name, savename(save_dict, "bson")), save_dict; safe=true)


text:
using Pkg
Pkg.activate(@__DIR__)

using Downloads, Tar, ArtifactUtils, Artifacts, Dates, Pkg

#-----------------------------------------------------------------------------# Start from scratch
for file in readdir(@__DIR__)
    if endswith(file, ".tar") || endswith(file, ".gz")
        rm(joinpath(@__DIR__, file), force=true)
    end
end

#-----------------------------------------------------------------------------# Templates
# Need to install plotly.py in order to get themes because they are generated in python
# run(`conda create -n conda_jl python conda`)
# ENV["CONDA_JL_HOME"] = "/opt/homebrew/Caskroom/miniconda/base/envs/conda_jl"
# Pkg.build("Conda")
# using Conda
# Conda.add("plotly")

dir = "/opt/homebrew/Caskroom/miniconda/base/envs/conda_jl/pkgs"
templates = joinpath(dir, filter(x -> startswith(x, "plotly"), readdir(dir))[1], "site-packages/plotly/package_data/templates")

run(`gzip $(Tar.create(templates, joinpath(@__DIR__, "templates.tar")))`)

#-----------------------------------------------------------------------------# Plotly.js
url = "https://cdn.plot.ly/plotly-2.11.0.min.js"
file = basename(url)
dir = mkpath(joinpath(@__DIR__, "plotlyjs"))
Downloads.download(url, joinpath(dir, file))

run(`gzip $(Tar.create(dir, joinpath(@__DIR__, "plotly.tar")))`)

#-----------------------------------------------------------------------------# upload
try
    artifacts_today = "artifacts_$(today())"

    run(`gh release create $artifacts_today $(joinpath(@__DIR__, "templates.tar.gz")) $(joinpath(@__DIR__, "plotly.tar.gz")) --title $artifacts_today --notes ""`)

    @info "Sleeping so artifacts are ready on GitHub..."
    sleep(10)
    add_artifact!(
        "Artifacts.toml",
        "plotly.min.js",
        "https://github.com/joshday/PlotlyLight.jl/releases/download/$artifacts_today/plotly.tar.gz",
        force=true,
    )
    add_artifact!(
        "Artifacts.toml",
        "plotly_templates",
        "https://github.com/joshday/PlotlyLight.jl/releases/download/$artifacts_today/templates.tar.gz",
        force=true,
    )
catch ex
    @error "Error (probably the release already exists): $ex"
end


text:
import Bridge: _scale, endpoint, samplepath, ProcessOrCoefficients, _b#, SamplePath
using StaticArrays

function forcedSolve!(em, Y, u, W::SamplePath,
                      P::ProcessOrCoefficients)
    forcedSolve!(em, Y, u, W, P, domain(P))
end

"""
    forcedSolve!(em, Y, u::T, W::SamplePath,
                 P::ProcessOrCoefficients, ::DiffusionDomain=domain(P))

Solve stochastic differential equation ``dX_t = b(t,X_t)dt + σ(t,X_t)dW_t``
using the Euler-Maruyama scheme in place. `forcedSolve!` as opposed to `solve!`
enforces adherence to the diffusion's domain (which numerical schemes are prone
to violate). By default, no restrictions are made, so calls `solve!`.
"""
function forcedSolve!(em::EulerMaruyamaBounded, Y, u::T, W::SamplePath,
                      P::ProcessOrCoefficients, ::UnboundedDomain) where {T}
    solve!(em, Y, u, W, P)
end

"""
    forcedSolve!(::EulerMaruyamaBounded, Y, u::T,
                 W::Union{SamplePath{SVector{D,S}},SamplePath{S}},
                 P::ProcessOrCoefficients, d::LowerBoundedDomain=domain(P))

Solve stochastic differential equation ``dX_t = b(t,X_t)dt + σ(t,X_t)dW_t``
using the Euler-Maruyama scheme in place. `forcedSolve!` as opposed to `solve!`
enforces adherence to the diffusion's domain (which numerical schemes are prone
to violate). This function enforces lower bounds by modifying `W` in place.
"""
function forcedSolve!(::EulerMaruyamaBounded, Y, u::T, W::SamplePath{S},
                      P::ProcessOrCoefficients, d::DiffusionDomain) where {S,T}
    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    ww = W.yy
    tt = Y.tt
    yy = Y.yy
    tt[:] = W.tt
    y::T = u

    offset = zero(S)
    for i in 1:N-1
        ww[..,i+1] += offset
        yy[.., i] = y
        dWt = ww[.., i+1]-ww[.., i]
        increm = _b((i,tt[i]), y, P)*(tt[i+1]-tt[i]) + _scale(dWt, σ(tt[i], y, P))
        y_new = y + increm
        offset_addon = zero(S)
        while !bound_satisfied(d, y_new)
            rootdt = √(tt[i+1]-tt[i])
            dWt = rootdt*randn(S)
            increm = _b((i,tt[i]), y, P)*(tt[i+1]-tt[i]) + _scale(dWt, σ(tt[i], y, P))
            y_new = y + increm
            offset_addon = ww[.., i] + dWt - ww[.., i+1]
        end
        offset += offset_addon
        ww[.., i+1] = ww[.., i] + dWt
        y = y_new
    end
    yy[.., N] = endpoint(y, P)
    Y
end

function forcedSolve(em, u::T, W::SamplePath,
                     P::ProcessOrCoefficients) where T
    WW = deepcopy(W)
    XX = samplepath(W.tt, zero(T))
    forcedSolve!(em, XX, u, WW, P)
    WW, XX
end

function forcedSolve(em, u::T, W::SamplePath,
                     P::GuidPropBridge) where T
    WW = deepcopy(W)
    XX = samplepath(W.tt, zero(T))
    forcedSolve!(em, XX, u, WW, P, domain(P.Target))
    WW, XX
end

import Bridge.solve!

function solve!(::Union{EulerMaruyamaBounded,Bridge.EulerMaruyama}, Y, u::T, W::SamplePath, P::ProcessOrCoefficients) where {T}
    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    ww = W.yy
    tt = Y.tt
    yy = Y.yy
    tt[:] = W.tt

    y::T = u
    dom = domain(P)
    for i in 1:N-1
        yy[.., i] = y
        y = y + _b((i,tt[i]), y, P)*(tt[i+1]-tt[i]) + _scale((ww[.., i+1]-ww[..,i]), σ(tt[i], y, P))
        !bound_satisfied(dom, y) && return false, nothing
    end
    yy[.., N] = endpoint(y, P)
    true, Y
end

#NOTE ↓↓↓↓↓↓↓↓↓↓↓ not the right way to go for now ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓

function solveAndll(::EulerMaruyama, u::T, W::SamplePath, P::GuidPropBridge, θ
                    ) where T
    WW = deepcopy(W)
    XX = samplepath(W.tt, zero(T))
    ll = solveAndll!(EulerMaruyamaBounded(), XX, u, WW, P, domain(P.Target), θ)
    WW, XX, ll
end

function solveAndll!(::EulerMaruyama, u::T, W::SamplePath, P::GuidPropBridge, θ
                     ) where T
    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    ww = W.yy
    tt = Y.tt
    yy = Y.yy
    tt[:] = W.tt
    y::T = u

    ll = 0.0
    for i in 1:N-1
        yy[.., i] = y
        dWt = ww[.., i+1]-ww[.., i]
        s = tt[i]
        dt = tt[i+1]-tt[i]
        b_prop = _b((i,s), y, P, θ)
        y = y + b_prop*dt + _scale(dWt, σ(s, y, P, θ))

        b_trgt = _b((i,s), y, target(P), θ)
        b_aux = _b((i,s), y, auxiliary(P), θ)
        rₜₓ = r((i,s), y, P, θ)
        ll += dot(b_trgt-b_aux, rₜₓ) * dt

        if !constdiff(P)
            Hₜₓ = H((i,s), y, P, θ)
            aₜₓ = a((i,s), y, target(P), θ)
            ãₜ = ã((i,s), y, P, θ)
            ll -= 0.5*sum( (aₜₓ - ãₜ).*Hₜₓ ) * dt
            ll += 0.5*( rₜₓ'*(aₜₓ - ãₜ)*rₜₓ ) * dt
        end
    end
    yy[.., N] = endpoint(y, P)
    ll
end


text:
"""
 Type: GenericStopping

 Methods: start!, stop!, update\\_and\\_start!, update\\_and\\_stop!, fill\\_in!, reinit!, status

 A generic Stopping to solve instances with respect to some
 optimality conditions. Optimality is decided by computing a score, which is then
 tested to zero.

 Tracked data include:
- pb         : A problem
- state      : The information relative to the problem, see `GenericState`.
- (opt) meta : Metadata relative to a stopping criterion, see `StoppingMeta`.
- (opt) main_stp : Stopping of the main loop in case we consider a Stopping
                       of a subproblem.
                       If not a subproblem, then `nothing`.
- (opt) listofstates : ListofStates designed to store the history of States.
- (opt) stopping_user_struct : Contains any structure designed by the user.

 Constructor: `GenericStopping(:: Any, :: AbstractState; meta :: AbstractStoppingMeta = StoppingMeta(), main_stp :: Union{AbstractStopping, Nothing} = nothing, stopping_user_struct :: Any = nothing, kwargs...)`

 Note: Metadata can be provided by the user or created with the Stopping
       constructor via kwargs. If a specific StoppingMeta is given and
       kwargs are provided, the kwargs have priority.

 Examples:
 GenericStopping(pb, GenericState(ones(2)), rtol = 1e-1)

 Besides optimality conditions, we consider classical emergency exit:
       - domain error        (for instance: NaN in x)
       - unbounded problem   (not implemented)
       - unbounded x         (x is too large)
       - tired problem       (time limit attained)
       - resources exhausted (not implemented)
       - stalled problem     (not implemented)
       - iteration limit     (maximum number of iteration (i.e. nb of stop) attained)
       - main_pb limit       (tired or resources of main problem exhausted)

 There is an additional default constructor which creates a Stopping with a default State.

 `GenericStopping(:: Any, :: Union{Number, AbstractVector}; kwargs...)`

 Note: Keywords arguments are forwarded to the classical constructor.

 Examples:
 GenericStopping(pb, x0, rtol = 1e-1)
"""
mutable struct GenericStopping{Pb, M, SRC, T, MStp, LoS} <: AbstractStopping{Pb, M, SRC, T, MStp, LoS}

  # Problem
  pb                   :: Pb

  # Problem stopping criterion
  meta                 :: M
  stop_remote          :: SRC

  # Current information on the problem
  current_state        :: T

  # Stopping of the main problem, or nothing
  main_stp             :: MStp

  # History of states
  listofstates         :: LoS

  # User-specific structure
  stopping_user_struct :: AbstractDict

end

function GenericStopping(pb            :: Pb,
                         meta          :: M,
                         stop_remote   :: SRC,
                         current_state :: T;
                         main_stp      :: AbstractStopping = VoidStopping(),
                         list          :: AbstractListofStates = VoidListofStates(),
                         user_struct   :: AbstractDict = Dict(),
                         kwargs...
                         ) where {Pb  <: Any,
                                  M   <: AbstractStoppingMeta,
                                  SRC <: AbstractStopRemoteControl,
                                  T   <: AbstractState}

  return GenericStopping(pb, meta, stop_remote, current_state,
                        main_stp, list, user_struct)
end

function GenericStopping(pb            :: Pb,
                         meta          :: M,
                         current_state :: T;
                         main_stp      :: AbstractStopping = VoidStopping(),
                         list          :: AbstractListofStates = VoidListofStates(),
                         user_struct   :: AbstractDict = Dict(),
                         kwargs...
                         ) where {Pb <: Any,
                                  M  <: AbstractStoppingMeta,
                                  T  <: AbstractState}

  stop_remote = StopRemoteControl() #main_stp == VoidStopping() ? StopRemoteControl() : cheap_stop_remote_control()

  return GenericStopping(pb, meta, stop_remote, current_state,
                        main_stp, list, user_struct)
end

function GenericStopping(pb            :: Pb,
                         current_state :: T;
                         main_stp      :: AbstractStopping = VoidStopping(),
                         list          :: AbstractListofStates = VoidListofStates(),
                         user_struct   :: AbstractDict = Dict(),
                         kwargs...
                         ) where {Pb <: Any, T <: AbstractState}

  meta = StoppingMeta(; kwargs...)
  stop_remote = StopRemoteControl() #main_stp == VoidStopping() ? StopRemoteControl() : cheap_stop_remote_control()

  return GenericStopping(pb, meta, stop_remote, current_state,
                        main_stp, list, user_struct)
end

function GenericStopping(pb            :: Pb,
                         stop_remote   :: SRC,
                         current_state :: T;
                         main_stp      :: AbstractStopping = VoidStopping(),
                         list          :: AbstractListofStates = VoidListofStates(),
                         user_struct   :: AbstractDict = Dict(),
                         kwargs...
                         ) where {Pb  <: Any,
                                  SRC <: AbstractStopRemoteControl,
                                  T   <: AbstractState}

  meta = StoppingMeta(; kwargs...)

  return GenericStopping(pb, meta, stop_remote, current_state,
                        main_stp, list, user_struct)
end

"""
`GenericStopping(pb :: Any, x :: T; n_listofstates :: Int = 0, kwargs...)`

Setting the keyword argument `n_listofstates > 0` initialize a ListofStates of length `n_listofstates`.
"""
function GenericStopping(pb :: Any, x :: T; n_listofstates :: Int = 0, kwargs...) where T
  state = GenericState(x)
  if n_listofstates > 0 && :list ∉ keys(kwargs)
    list = ListofStates(n_listofstates, Val{typeof(state)}())
    return GenericStopping(pb, state, list = list; kwargs...)
  end
  return GenericStopping(pb, state; kwargs...)
end

"""
update!: generic update function for the Stopping

`update!(:: AbstractStopping; kwargs...)`

Shortcut for update!(stp.current_state; kwargs...)
"""
function update!(stp :: AbstractStopping; kwargs...)
  return update!(stp.current_state; kwargs...)
end

"""
fill_in!: fill in the unspecified values of the AbstractState.

`fill_in!(:: AbstractStopping, x :: Union{Number, AbstractVector})`

Note: NotImplemented for Abstract/Generic-Stopping.
"""
function fill_in!(stp :: AbstractStopping, x :: T) where T
  return throw(error("NotImplemented function"))
end

"""
update\\_and\\_start!: update values in the State and initialize the Stopping.
Returns the optimality status of the problem as a boolean.

`update_and_start!(:: AbstractStopping; kwargs...)`

 Note:
  - Kwargs are forwarded to the `update!` call.
  - `no_opt_check` skip optimality check in `start!` (`false` by default).
"""
function update_and_start!(stp          :: AbstractStopping;
                           no_opt_check :: Bool = false,
                           kwargs...)

  if stp.stop_remote.cheap_check
    _smart_update!(stp.current_state; kwargs...)
  else
    update!(stp; kwargs...)
  end
  OK = start!(stp, no_opt_check = no_opt_check)

  return OK
end

"""
 Update the Stopping and return `true` if we must stop.

 `start!(:: AbstractStopping; no_opt_check :: Bool = false, kwargs...)`

 Purpose is to know if there is a need to even perform an optimization algorithm
 or if we are at an optimal solution from the beginning.
 Set `no_opt_check` to `true` avoid checking optimality and domain errors.

 The function `start!` successively calls: `_domain_check(stp, x)`,
 `_optimality_check!(stp, x)`, `_null_test(stp, x)` and
 `_user_check!(stp, x, true)`.

 Note: - `start!` initializes `stp.meta.start_time` (if not done before),
 `stp.current_state.current_time` and `stp.meta.optimality0`
 (if `no_opt_check` is false).
       - Keywords argument are passed to the `_optimality_check!` call.
       - Compatible with the `StopRemoteControl`.
"""
function start!(stp :: AbstractStopping;
                no_opt_check :: Bool = false,
                kwargs...)

  state = stp.current_state
  src   = stp.stop_remote

  #Initialize the time counter
  if src.tired_check && isnan(stp.meta.start_time)
    stp.meta.start_time = time()
  end
  #and synchornize with the State
  if src.tired_check && isnan(state.current_time)
    _update_time!(state, stp.meta.start_time)
  end

  if !no_opt_check
    stp.meta.domainerror = if src.domain_check
                              #don't check current_score
                              _domain_check(stp.current_state,
                                            current_score = true)
                            else
                              stp.meta.domainerror
                            end
    if !stp.meta.domainerror && src.optimality_check
      optimality0          = _optimality_check!(stp; kwargs...)
      norm_optimality0     = norm(optimality0, Inf)
      if src.domain_check && isnan(norm_optimality0)
        stp.meta.domainerror = true
      elseif norm_optimality0 == Inf
        stp.meta.optimality0 = one(typeof(norm_optimality0))
      else
        stp.meta.optimality0 = norm_optimality0
      end

      if _null_test(stp, optimality0) stp.meta.optimal = true end
    end
  end

  src.user_start_check && _user_check!(stp, state.x, true)

  OK = OK_check(stp.meta)

  #do nothing if typeof(stp.listofstates) == VoidListofStates
  add_to_list!(stp.listofstates, stp.current_state)

  return OK
end

"""
 reinit!: reinitialize the meta-data in the Stopping.

 `reinit!(:: AbstractStopping; rstate :: Bool = false, kwargs...)`

 Note:
- If `rstate` is set as `true` it reinitializes the current State
(with the kwargs).
- If `rlist` is set as true the list of states is also reinitialized, either
set as a `VoidListofStates` if `rstate` is `true` or a list containing only the current
state otherwise.
"""
function reinit!(stp    :: AbstractStopping;
                 rstate :: Bool = false,
                 rlist  :: Bool = true,
                 kwargs...)

  stp.meta.start_time  = NaN
  stp.meta.optimality0 = 1.0

  #reinitialize the boolean status
  reinit!(stp.meta)

  #reinitialize the counter of stop
  stp.meta.nb_of_stop = 0

  #reinitialize the list of states
  if rlist && (typeof(stp.listofstates) != VoidListofStates)
    #TODO: Warning we cannot change the type of ListofStates
    stp.listofstates = rstate ? VoidListofStates() : ListofStates(stp.current_state)
  end

  #reinitialize the state
  if rstate
    reinit!(stp.current_state; kwargs...)
  end

  return stp
end

"""
update\\_and\\_stop!: update the values in the State and
return the optimality status of the problem as a boolean.

`update_and_stop!(stp :: AbstractStopping; kwargs...)`

Note: Kwargs are forwarded to the `update!` call.
"""
function update_and_stop!(stp :: AbstractStopping; kwargs...)

  if stp.stop_remote.cheap_check
    _smart_update!(stp.current_state; kwargs...)
    OK = cheap_stop!(stp)
  else
    update!(stp; kwargs...)
    OK = stop!(stp)
  end

 return OK
end

"""
stop!: update the Stopping and return a boolean true if we must stop.

`stop!(:: AbstractStopping; kwargs...)`

It serves the same purpose as `start!` in an algorithm; telling us if we
stop the algorithm (because we have reached optimality or we loop infinitely,
etc).

The function `stop!` successively calls: `_domain_check`, `_optimality_check`,
`_null_test`, `_unbounded_check!`, `_tired_check!`, `_resources_check!`,
`_stalled_check!`, `_iteration_check!`, `_main_pb_check!`, `add_to_list!`

Note:
- kwargs are sent to the `_optimality_check!` call.
- If `listofstates != VoidListofStates`, call `add_to_list!`.
"""
function stop!(stp          :: AbstractStopping;
               no_opt_check :: Bool = false,
               kwargs...)

  x   = stp.current_state.x
  src = stp.stop_remote

  src.unbounded_and_domain_x_check && _unbounded_and_domain_x_check!(stp, x)
  stp.meta.domainerror = if src.domain_check
                            #don't check x and current_score
                            _domain_check(stp.current_state, x = true,
                                                          current_score = true)
                          else
                            stp.meta.domainerror
                          end
  if !no_opt_check && !stp.meta.domainerror
    # Optimality check
    if src.optimality_check
      score = _optimality_check!(stp; kwargs...)
      if any(isnan, score)
        stp.meta.domainerror = true
      end
      if _null_test(stp, score) stp.meta.optimal = true end
    end

    src.infeasibility_check     && _infeasibility_check!(stp, x)
    src.unbounded_problem_check && _unbounded_problem_check!(stp, x)
    src.tired_check             && _tired_check!(stp, x)
    src.resources_check         && _resources_check!(stp, x)
    src.stalled_check           && _stalled_check!(stp, x)
    src.iteration_check         && _iteration_check!(stp, x)
    src.main_pb_check           && _main_pb_check!(stp, x)
    src.user_check              && _user_check!(stp, x)
  end

  OK = OK_check(stp.meta)

  _add_stop!(stp)

  #do nothing if typeof(stp.listofstates) == VoidListofStates
  add_to_list!(stp.listofstates, stp.current_state)

  return OK
end

"""
stop!: update the Stopping and return a boolean true if we must stop.

`cheap_stop!(:: AbstractStopping; kwargs...)`

It serves the same purpose as `stop!`, but avoids any potentially expensive checks.
We no longer browse `x` and `res` in the State, and no check on the `main_stp`.
Check only the updated entries in the meta.

The function `cheap_stop!` successively calls:
`_null_test`, `_unbounded_check!`, `_tired_check!`, `_resources_check!`,
`_stalled_check!`, `_iteration_check!`, `add_to_list!`

Note:
- kwargs are sent to the `_optimality_check!` call.
- If `listofstates != VoidListofStates`, call `add_to_list!`.
"""
function cheap_stop!(stp :: AbstractStopping; kwargs...)

  x   = stp.current_state.x
  src = stp.stop_remote

  # Optimality check
  if src.optimality_check
    score = _optimality_check!(stp; kwargs...) #update state.current_score
    if _null_test(stp, score) stp.meta.optimal = true end
  end
  OK = stp.meta.optimal

  OK = OK || (src.infeasibility_check     && _infeasibility_check!(stp, x))
  OK = OK || (src.unbounded_problem_check && _unbounded_problem_check!(stp, x))
  OK = OK || (src.tired_check             && _tired_check!(stp, x))
  OK = OK || (src.resources_check         && _resources_check!(stp, x))
  OK = OK || (src.iteration_check         && _iteration_check!(stp, x))
  OK = OK || (src.user_check              && _user_check!(stp, x))

  _add_stop!(stp)

  #do nothing if typeof(stp.listofstates) == VoidListofStates
  add_to_list!(stp.listofstates, stp.current_state)

  return OK
end

"""
\\_add\\_stop!: increment a counter of stop.

`_add_stop!(:: AbstractStopping)`

Fonction called everytime `stop!` is called. In theory should be called once
every iteration of an algorithm.

Note: update `meta.nb_of_stop`.
"""
function _add_stop!(stp :: AbstractStopping)

  stp.meta.nb_of_stop += 1

  return stp
end

"""
\\_iteration\\_check!: check if the optimization algorithm has reached the
iteration limit.

`_iteration_check!(:: AbstractStopping,  :: T)`

Note: Compare `meta.iteration_limit` with `meta.nb_of_stop`.
"""
function _iteration_check!(stp :: AbstractStopping,
                           x   :: T) where T

  max_iter = stp.meta.nb_of_stop >= stp.meta.max_iter
  if max_iter stp.meta.iteration_limit = true end

  return stp.meta.iteration_limit
end

"""
\\_stalled\\_check!: check if the optimization algorithm is stalling.

`_stalled_check!(:: AbstractStopping, :: T)`

Note: Do nothing by default for AbstractStopping.
"""
function _stalled_check!(stp :: AbstractStopping,
                         x   :: T) where T
  return stp.meta.stalled
end

"""
\\_tired\\_check!: check if the optimization algorithm has been running
 for too long.

`_tired_check!(:: AbstractStopping, :: T)`

Note:
  - Return `false` if `meta.start_time` is `NaN` (by default).
  - Update `meta.tired`.
"""
function _tired_check!(stp :: AbstractStopping,
                       x   :: T) where T

  stime = stp.meta.start_time #can be NaN
  ctime = time()

  #Keep the current_state updated
  _update_time!(stp.current_state, ctime)

  elapsed_time = ctime - stime
  max_time     = elapsed_time > stp.meta.max_time #NaN > 1. is false

  if max_time stp.meta.tired = true end

  return stp.meta.tired
end

function _tired_check!(stp :: VoidStopping,
                       x   :: T) where T
  return false
end

"""
\\_resources\\_check!: check if the optimization algorithm has exhausted
 the resources.

`_resources_check!(:: AbstractStopping, :: T)`

Note: Do nothing by default `meta.resources` for AbstractStopping.
"""
function _resources_check!(stp :: AbstractStopping,
                           x   :: T) where T
  return stp.meta.resources
end

function _resources_check!(stp :: VoidStopping,
                           x   :: T) where T
  return false
end

"""
\\_main\\_pb\\_check!: check the resources and the time of the upper problem
 if `main_stp != VoidStopping`.

`_main_pb_check!(:: AbstractStopping, :: T)`

Note: - Modify the meta of the `main_stp`.
      - return `false` for `VoidStopping`.
"""
function _main_pb_check!(stp :: AbstractStopping,
                         x   :: T) where T

  max_time     = _tired_check!(stp.main_stp, x)
  resources    = _resources_check!(stp.main_stp, x)
  main_main_pb = _main_pb_check!(stp.main_stp, x)

  check = max_time || resources || main_main_pb

  if check stp.meta.main_pb = true end

  return stp.meta.main_pb
end

function _main_pb_check!(stp :: VoidStopping,
                         x   :: T) where T
 return false
end

"""
\\_unbounded\\_and\\_domain\\_x\\_check!: check if x gets too big,
and if it has NaN or missing values.

`_unbounded_and_domain_x_check!(:: AbstractStopping, :: T)`

Note:
- compare `||x||_∞` with `meta.unbounded_x` and update `meta.unbounded`.
- it also checks `NaN` and `missing` and update `meta.domainerror`.
- short-circuiting if one of the two is `true`.
"""
function _unbounded_and_domain_x_check!(stp :: AbstractStopping,
                                        x   :: T) where T

  bigX(z :: eltype(T)) = (abs(z) >= stp.meta.unbounded_x)
  (stp.meta.unbounded, stp.meta.domainerror) = _large_and_domain_check(bigX, x)
  return stp.meta.unbounded || stp.meta.domainerror
end

function _large_and_domain_check(f, itr)
    for x in itr
        v   = f(x)
        w   = ismissing(x) || isnan(x)
        if w
          return (false, true)
        elseif v
          return (true, false)
        end
    end
    return (false, false)
end

"""
\\_unbounded\\_problem!: check if problem relative informations are unbounded

`_unbounded_problem_check!(:: AbstractStopping, :: T)`

Note: Do nothing by default.
"""
function _unbounded_problem_check!(stp :: AbstractStopping,
                                   x   :: T) where T

 return stp.meta.unbounded_pb
end

"""
\\_infeasibility\\_check!: check if problem is infeasible

`_infeasibility_check!(:: AbstractStopping, :: T)`

Note: `meta.infeasible` is `false` by default.
"""
function _infeasibility_check!(stp :: AbstractStopping,
                               x   :: T) where T

 return stp.meta.infeasible
end

"""
\\_optimality\\_check: compute the optimality score.

`_optimality_check!(:: AbstractStopping; kwargs...)`

"""
function _optimality_check!(stp :: AbstractStopping{Pb, M, SRC, T,
                                                    MStp, LoS};
                           kwargs...) where {Pb, M, SRC, T, MStp, LoS}

  setfield!(stp.current_state, :current_score,
            stp.meta.optimality_check(stp.pb, stp.current_state; kwargs...))

  return stp.current_state.current_score
end

"""
\\_null\\_test: check if the score is close enough to zero
(up to some precisions found in the meta).

`_null_test(:: AbstractStopping, :: T)`

Note:
- the second argument is compared with
`meta.tol_check(meta.atol, meta.rtol, meta.optimality0)`,
and `meta.tol_check_neg(meta.atol, meta.rtol, meta.optimality0)`.
- Compatible sizes is not verified.
"""
function _null_test(stp  :: AbstractStopping, optimality :: T) where T

  check_pos, check_neg = tol_check(stp.meta)
  optimal = _inequality_check(optimality, check_pos, check_neg)

  return optimal
end

#remove the Missing option here
_inequality_check(opt :: Number, check_pos :: Number, check_neg :: Number) = (opt <= check_pos) && (opt >= check_neg)
_inequality_check(opt, check_pos :: Number, check_neg :: Number) :: Bool = !any(z->(ismissing(z) || (z > check_pos) || (z < check_neg)), opt)
function _inequality_check(opt :: T, check_pos :: T, check_neg :: T) where T

  size_check = try
                 n = size(opt)
                 ncp, ncn = size(check_pos), size(check_neg)
                 n != ncp || n != ncn
               catch
                 false
               end

  if size_check
    throw(ErrorException("Error: incompatible size in _null_test wrong size of optimality, tol_check and tol_check_neg"))
  end

  for (o, cp, cn) in zip(opt, check_pos, check_neg)
    v = o > cp || o < cn
    if v
      return false
    end
  end

  return true
end

"""
\\_user\\_check: nothing by default.

`_user_check!( :: AbstractStopping, x :: T, start :: Bool)`

Call the `user_check_func!(:: AbstractStopping, :: Bool)` from the meta.
The boolean `start` is `true` when called from the `start!` function.
"""
function _user_check!(stp :: AbstractStopping, x :: T, start :: Bool) where T
  #callback function furnished by the user
  stp.meta.user_check_func!(stp, start)

  return stp.meta.stopbyuser
end

function _user_check!(stp :: AbstractStopping, x :: T) where T
  return _user_check!(stp, x, false)
end

const status_meta_list = Dict([(:Optimal, :optimal),
             (:SubProblemFailure, :fail_sub_pb),
             (:SubOptimal, :suboptimal),
             (:Unbounded, :unbounded),
             (:UnboundedPb, :unbounded_pb),
             (:Stalled, :stalled),
             (:IterationLimit, :iteration_limit),
             (:TimeLimit, :tired),
             (:EvaluationLimit, :resources),
             (:ResourcesOfMainProblemExhausted, :main_pb),
             (:Infeasible, :infeasible),
             (:StopByUser, :stopbyuser),
             (:Exception, :exception),
             (:DomainError, :domainerror)])

"""
status: returns the status of the algorithm:

`status(:: AbstractStopping; list = false)`

The different statuses are:
- Optimal: reached an optimal solution.
- SubProblemFailure
- SubOptimal: reached an acceptable solution.
- Unbounded: current iterate too large in norm.
- UnboundedPb: unbouned problem.
- Stalled: stalled algorithm.
- IterationLimit: too many iterations of the algorithm.
- TimeLimit: time limit.
- EvaluationLimit: too many ressources used,
                          i.e. too many functions evaluations.
- ResourcesOfMainProblemExhausted: in the case of a substopping, EvaluationLimit or TimeLimit
  for the main stopping.
- Infeasible: default return value, if nothing is done the problem is
               considered feasible.
- StopByUser: stopped by the user.
- DomainError: there is a NaN somewhere.
- Exception: unhandled exception
- Unknwon: if stopped for reasons unknown by Stopping.

Note:
- Set keyword argument `list` to true, to get an `Array` with all the statuses.
- The different statuses correspond to boolean values in the meta.
"""
function status(stp :: AbstractStopping; list = false)

  if list
    list_status = findall(x -> getfield(stp.meta, x), status_meta_list)
    if list_status == zeros(0) list_status = [:Unknown] end
  else
    list_status = findfirst(x -> getfield(stp.meta, x), status_meta_list)
    if isnothing(list_status) list_status = :Unknown end
  end

  return list_status
end

"""
elapsed_time: returns the elapsed time.

`elapsed_time(:: AbstractStopping)`

`current_time` and `start_time` are NaN if not initialized.
"""
function elapsed_time(stp :: AbstractStopping)
  return stp.current_state.current_time - stp.meta.start_time
end


text:
# Example with 3 independent variables (constant, dummy, and continuous variable)
using StatsPlots

N = 1000;
X = [ones(N) rand(0:1, N) randn(N)];
Time, Indicator, λ, β, α = simulate([0.5, 1.0, -1.0], .5, X);
b, a = estimate(Time, Indicator, X, MvNormal([0, 0, 0], 100I), Exponential(1), 2_500);

# Plotting
gr(size = (800, 600))
plot(layout = (3, 1), label = false)

plot!(b, color = [:blue :red :green], subplot = 1, label = false)
hline!(β', color = [:blue :red :green], label = false, subplot = 1)

plot!(a, color = :grey, subplot = 2, label = false)
hline!([α], color = :black, label = false, linewidth = 2, subplot = 2)

Fit = lambda(mean(b[501:2500, :]; dims = 1)[:], mean(a[501:2500]), X)
#scatter!(Time, Fit, subplot = 3, label = false)
#plot!(0:10, 0:10, subplot = 3, label = false)
plot!([Time Fit], subplot = 3, label = false)


text:
function estimateBM(xs::Vector{Float64}, b::Int64)
  lxs::Int64 = length(xs)
  a = floor(Int64, lxs/b)
  @assert a > 1
  n::Int64 = a*b
  start::Int64 = length(xs) - n
  overallMean::Float64 = 0.0
  for i = 1:n
    overallMean += xs[start+i]
  end
  overallMean /= n
  acc::Float64 = 0.0
  for i = 1:a
    batchAcc::Float64 = 0.0
    batchStart::Int64 = start + (i-1)*b
    for j = 1:b
      batchAcc += xs[batchStart + j]
    end
    tmp::Float64 = batchAcc/b
    tmp -= overallMean
    acc += tmp * tmp
  end
  return b/(a-1)*acc
end

## Basic batch means estimation of the asymptotic variance
function estimateBM(xs::Vector{Float64})
  return estimateBM(xs, floor(Int64, sqrt(length(xs))))
end


text:
# Passing Julia callback functions to Python routines.
#
# Note that this will typically involve two functions: the
# desired Julia function/closure, and a top-level C-callable
# wrapper function used with PyCFunction_NewEx -- the latter
# is called from Python and calls the former as needed.

################################################################

# To pass an arbitrary Julia object to Python, we wrap it
# in a jlwrap Python class, where jlwrap.__call__
# executes the pyjlwrap_call function in Julia, which
# in turn fetches the actual Julia value from the python wrapper
# and calls it.

# convert Python args to Julia; overridden below for a FuncWrapper type
# that allows the user to specify the argument types.
julia_args(f, args) = convert(PyAny, args)
julia_kwarg(f, kw, arg) = convert(PyAny, arg)

################################################################
# allow the user to convert a Julia function into a Python
# function with specified argument types, both to give more control
# and to avoid the overhead of type introspection.

struct FuncWrapper{T,F}
    f::F
    kwargs::Dict{Symbol,Any}
end
(f::FuncWrapper)(args...; kws...) = f.f(args...; kws...)
julia_args(f::FuncWrapper{T}, args) where {T} = convert(T, args)
julia_kwarg(f::FuncWrapper, kw, arg) = convert(get(f.kwargs, kw, PyAny), arg)

FuncWrapper(f, kwargs, argtypes::Type) = FuncWrapper{argtypes,typeof(f)}(f, kwargs)

"""
    pyfunction(f, argtypes...; kwtypes...)

Create a Python object that wraps around the Julia function (or callable
object) `f`.   Unlike `PyObject(f)`, this allows you to specify the argument
types that the Julia function expects — giving you more control and potentially
better performance.

`kwtypes...` should be a set of `somekeyword=SomeType` arguments giving
the desired Julia type for each keyword `somekeyword`.  Unspecified keywords
are converted to `PyAny` (i.e. auto-detected) by default.

The return value `ret = f(...)` is still converted back to a Python object by
`PyObject(ret)`.   If you want a different return-type conversion than the default
of `PyObject(ret)`, you can instead call `pyfunctionret`.
"""
function pyfunction(f, argtypes...; kwtypes...)
    kwdict = Dict{Symbol,Any}(k => v for (k,v) in kwtypes)
    return pyjlwrap_new(FuncWrapper(f, kwdict, Tuple{argtypes...}))
end

"""
    pyfunctionret(f, returntype, argtypes...; kwtypes...)

Like `pyfunction`, but also lets you specify the `returntype` for
conversion back to Python.   In particular, if `ret = f(...)` is
the return value of `f`, then it is converted to Python via
`PyObject(returntype(ret))`.

If `returntype` is `Any`, then `ret` is not converted to a "native"
Python type at all, and is instead converted to a "wrapped" Julia
object in Python.  If `returntype` is `nothing`, then the return value
is discarded and `nothing` is returned to Python.
"""
function pyfunctionret(f, returntype, argtypes...; kwtypes...)
    if returntype === Any
        return pyfunction((args...; kws...) -> pyjlwrap_new(f(args...; kws...)),
                          argtypes...; kwtypes...)
    elseif returntype === nothing
        return pyfunction((args...; kws...) -> begin f(args...; kws...); nothing; end,
                          argtypes...; kwtypes...)
    else
        return pyfunction((args...; kws...) -> returntype(f(args...; kws...)),
                          argtypes...; kwtypes...)
    end
end


text:
# A simple 2D transmission full waveform inversion experiment
# to showcase how this framework operates. Much of the
#

using JOLI   # Linear operators
using PyPlot # For plotting
using WAVEFORM

# Turn off annoying JOLI type checks
JOLI.jo_type_mismatch_error_set(false)
JOLI.jo_type_mismatch_warn_set(false)

# number of model points in each direction (z,x)
n = 101 .* [1;1];

# grid spacing in metres
d = 10.0 .*[1;1];

# location of the coordinate axis origin
o = 0.0 .* [1;1];

# Frequencies
freqs = [5.0;10.0;15.0];

# Source wavelet time shift
t0 = 0.0;

# Peak frequency of Ricker Wavelet
f0 = 10.0;

# Units of velocity parameter vector
unit = "m/s";

# Background velocity (m/s)
vel_background = 2000;

# Model domain range
L = n.*d;

# Source grid definition
xsrc = 0.0:100.0:L[2];
ysrc = [0.0];
zsrc = [10.0];

# Receiver grid definition
xrec = 0:10.0:L[2];
yrec = [0.0];
zrec = [950.0];

# Model type containing the domain geometry
model = Model{Int64,Float64}(n,d,o,t0,f0,unit,freqs,xsrc,ysrc,zsrc,xrec,yrec,zrec);
nsrc = length(model.xsrc); nfreq = length(model.freq);

comp_n = n;
comp_d = d;
comp_o = o;
npml = 20*ones(Int64,2,2);

# Stencil used to discretize the PDE
pde_scheme = helm2d_chen9p;

# If true,
cut_pml = true;

implicit_matrix = false;

# Objective function
misfit = least_squares;

# Binary mask to indicate which sources/frequencies to use (default: all of them)
srcfreqmask = trues(nsrc,nfreq);

# Use a direct solver for 2D problems
lsopts = LinSolveOpts(solver=:lufact);

# Options for discretizing + solving the Helmholtz PDEs
opts = PDEopts{Int64,Float64}(pde_scheme,comp_n,comp_d,comp_o,cut_pml,implicit_matrix,npml,misfit,srcfreqmask,lsopts);

v0 = vel_background*ones(n...);
v = copy(v0);
v[div(n[1],3):2*div(n[1],3),div(n[2],3):2*div(n[2],3)] = 1.25.*vel_background;
v = vec(v);
v0 = vec(v0);

# True velocity
imshow(v);

Q = diagm(ones(nsrc));
D = forw_model(v,Q,model,opts);

# Plot a frequency slice
imshow(real(D[:,1:nsrc]),aspect="auto");

obj! = construct_pde_misfit(v,Q,D,model,opts,batch_mode=false)

using OptimPackNextGen
proj! = (xproj,x)->project_bounds!(x,minimum(v),maximum(v),xproj);
vest = spg(obj!,proj!,v0,3,maxfc=50,verb=true)
imshow(reshape(vest,model.n...),vmin=minimum(v),vmax=maximum(v))


text:
module HilbertSchmidtIndependenceCriterion

  using StatsFuns, Distances

  # load Base modules
  using Statistics, LinearAlgebra
  import StatsBase: sample
	# includes
  include("common.jl")
  include("gammaHSIC.jl")

	# function exports
  export gammaHSIC, estimateKernelSize, rbfDotProduct

end # module


text:
image_distance(img1::AbstractMatrix, img2::AbstractMatrix) = mse(img1, img2)
image_distance(
    img1::AbstractMatrix{<:Gray},
    img2::AbstractMatrix{<:Gray},
    weights::AbstractMatrix,
    weights_sum = sum(weights)
) = sum(((float(img1) .- img2) .^ 2) .* weights) / weights_sum
image_distance(
    img1::AbstractMatrix{<:RGB},
    img2::AbstractMatrix,
    weights::AbstractMatrix,
    weights_sum = sum(weights)
) = sum((norm.(float(img1) .- img2, (2f0,)) .^ 2) .* weights) / weights_sum
image_distance(
    img1::AbstractMatrix{<:Gray},
    img2::AbstractMatrix{<:RGB},
    weights::AbstractMatrix,
    weights_sum = sum(weights)
) = image_distance(img1, Gray.(img2), weights, weights_sum)

function template_match_all(
    template, img;
    indices = axes(img),
    σ = 0.5f0,
    mask = nothing,
    border = Fill(1)
)
    if σ > 0
        template, img = blur(template, σ), blur(img, σ)
    end
    h, w = (size(template) .+ 1) .÷ 2 .* 2 .- 1
    template = @view template[1:h, 1:w]
    f = if isnothing(mask)
        (window) -> image_distance(window, template)
    else
        mask = @view mask[1:h, 1:w]
        (window) -> image_distance(window, template, mask)
    end
    window_size = size(template)
    # `mapwindow` is not stable for one line of pixels:
    # https://github.com/JuliaImages/ImageFiltering.jl/issues/169
    dists = mapwindow(
        f, img, window_size,
        indices = indices,
        border = border
    )
    perms = sortperm(reshape(dists, :))
    perms, dists
end

"""
    template_match(template, image; indices = axes(image) σ = 0.5f0, mask = nothing, border = Fill(1))

Match `template` in `image`. `indices` are the points of the centers of windows in the image.

If `σ ≤ 0`, no blur is applied.
"""
function template_match(
    template, img;
    indices = axes(img),
    σ = 0.5f0,
    mask = nothing,
    border = Fill(1)
)
    perms, dists = template_match_all(
        template, img;
        indices = indices,
        σ = σ,
        mask = mask,
        border = border
    )
    perms[1], dists[perms[1]]
end

function get_table_size(img, block_size)
    th, tw = size(img)
    bh, bw = block_size
    nrows = th ÷ bh
    ncols = tw ÷ bw
    (nrows, ncols)
end

"""
    table_search(
        template, table;
        block_size = size(template),
        rect = ((1, 1), block_size),
        indices = 1:block_count
    )

Find `template` in `table`. `rect` = (topleft, size) is the region of the block to search.
"""
function table_search(
    template, table;
    block_size = size(template),
    table_size = get_table_size(table, block_size),
    rect = Rect((1, 1), block_size),
    indices = nothing,
    mask = nothing
)
    if isnothing(indices)
        indices = 1:(table_size[1] * table_size[2])
    end

    rect_size = size(rect)
    if size(template) != rect_size
        template = imresize(template, rect_size)
    end

    _, closest_i = find_closest(indices) do i
        block_img = block(table, i, table_size, block_size)
        subsection = subimage(block_img, rect)
        dist = if isnothing(mask)
            image_distance(subsection, template)
        else
            block_mask = block(mask, i, table_size, block_size)
            submask = subimage(block_mask, rect)
            image_distance(subsection, template, submask)
        end
        dist, i
    end
    closest_i
end

struct SpriteSheet{T <: Union{Gray{Float32}, RGB{Float32}}, TM <: Union{Matrix{Float32}, Nothing}}
    image::Matrix{T}
    mask::TM
    block_size::Tuple{Int, Int}
    # (nrows, ncols)
    table_size::Tuple{Int, Int}
    data::OrderedDict{Int, String}
    filename::String
end

function SpriteSheet(filename; gray = false)
    filename, ext = splitext(filename)
    if lowercase(ext) ∉ ("", ".png")
        throw("Should be a png file: $filename")
    end
    img = load_image("$filename.png", gray = gray)
    mask = try
        reinterpret(Float32, float(load_image("$filename.mask.png", gray = true))) |> copy
    catch
        nothing
    end
    raw_data, block_size, indices = open("$filename.yaml") do fi
        data = YAML.load(fi)
        data["data"], data["block-size"], get(data, "indices", nothing)
    end
    block_size = (block_size[1], block_size[2])
    table_size = get_table_size(img, block_size)
    data = OrderedDict{Int, String}()
    if isnothing(indices)
        indices = 1:length(raw_data)
    end
    @assert length(indices) ≡ length(raw_data)
    for (i, s) in zip(indices, raw_data)
        data[i] = string(s)
    end
    SpriteSheet(float(img), mask, block_size, table_size, data, filename)
end

@forward SpriteSheet.data Base.length, Base.getindex
image(s::SpriteSheet) = s.image
Base.show(io::IO, mime::MIME"image/png", s::SpriteSheet; kwargs...) = show(io, mime, s.image; kwargs...)

block(img::AbstractMatrix, i, table_size, block_size) = let (_, ncols) = table_size, (bh, bw) = block_size
    subimage(
        img,
        (
            (i - 1) ÷ ncols * bh + 1,
            ((i - 1) % ncols) * bw + 1
        ),
        (bh, bw)
    )
end
block(s::SpriteSheet, i) = block(s.image, i, s.table_size, s.block_size)

function table_search(
    template, sheet::SpriteSheet;
    rect = Rect((1, 1), sheet.block_size),
    indices = keys(sheet.data)
)
    i = table_search(
        template, sheet.image,
        block_size = sheet.block_size,
        table_size = sheet.table_size,
        rect = rect,
        indices = indices,
        mask = sheet.mask
    )
    sheet[i]
end


text:
function test_api_endpoint(request::HTTP.Request)::HTTP.Response

    # initialize -
    response_json_dictionary = Dict{String,Any}()
    response_code = 200 # I'm an optimist ... default: all is good in the world

    try

        # First, lets get the body of the request, and then build the request_body_dictionary -
        # note: if this call fails, it will be captured in the catch - this is an insanely
        # huge malfunction, so a 500 is going back to the user ...
        request_body_dictionary = JSON.parse(String(request.body))

        # debug: did we we get anyting from the server?
        with_logger(simpleLogger) do
            @info("rcvd: $(request_body_dictionary)")
            flush(io)
        end

        # test logic: takes a string, reverses it and creates a nee field in the respone
        original_string = request_body_dictionary["to"];
        new_string = reverse(original_string)

        # package -
        response_json_dictionary["string_sent_to_server"] = original_string
        response_json_dictionary["string_sent_from_server"] = new_string

        # wrap and send -
        # encode -
        buffer = JSON.json(response_json_dictionary)
        response = HTTP.Response(response_code, buffer)

        # return -
        return response
    catch error

        # ooops ... if we get here, then something super bad happend.
        # send 500 back to user, and log
        response_code = 500

        # can we get some error information please?
        msg = sprint(showerror, error, catch_backtrace())

        # log -
        with_logger(simpleLogger) do
            @error(msg)
            flush(io)
        end

        # check: give specific error messages per error type
        response_json_dictionary["error"] = "Hmmm. that's not right"
        response_json_dictionary["status"] = 0
        if (isa(error, KeyError) == true)
            response_json_dictionary["error"] = "Key error: check your request for the correct keys"
        end

        # encode -
        buffer = JSON.json(response_json_dictionary)
        response = HTTP.Response(response_code, buffer)

        # return -
        return response
    end
end


text:
using GraphQLClient
using GraphQLParser
using HTTP
using JSON3
using StructTypes
using Test

include("utils.jl")

@testset "GraphQLClient" begin
    @testset "Client" begin include("client.jl") end
    @testset "Introspection" begin include("introspection.jl") end
    @testset "Schema Utils" begin include("schema_utils.jl") end
    @testset "Output Fields" begin include("output_fields.jl") end
    @testset "Args" begin include("args.jl") end
    @testset "Variables" begin include("variables.jl") end
    @testset "GQLResponse" begin include("gqlresponse.jl") end
    @testset "Subscriptions" begin include("subscriptions.jl") end
    @testset "HTTP Execution" begin include("http_execution.jl") end
    @testset "Type Construction" begin include("type_construction.jl") end
    @testset "GQL String" begin include("gql_string.jl") end
end


text:
# Algorithm 7.1
# Jacobi Iteration

using LinearAlgebra

function jacobi(mat, approx, known, maxIter, tol)
	n, m = size(mat)
	sol = Vector(undef, n)
	a = mat
	for k in 1:maxIter
		for i in 1:n
			sum = 0
			for j in 1:n
				if j != i
					sum += (a[i, j]*approx[j]) + known[i]
				end
			end
			sol[i] = (-sum)/a[i, i]
		end
		if LinearAlgebra.norm(sol - approx) < tol
			println("Success")
			break
		end
		for i in 1:n
			approx[i] = sol[i]
		end
	end
	return sol
end

A = [[1.0 1.33]; [2.0 -2.55];]
x = [0.0, 0.0]
b = [2.0, 5.0]
println(jacobi(A, x, b, 50, 0.01))


text:
function append_tag!(tags, kind, info)
    isnothing(info) && return
    tag = ElementNode("tag")
    link!(tag, AttributeNode("name", kind))
    link!(tag, AttributeNode("value", "$info"))
    link!(tags, tag)
end

clean_artist_name(artist_name) = replace(artist_name, r" \(\d\)$" => "")

function generate_xml(release::Release)
    docs = []
    position = 1
    for track in release.tracklist
        isnothing(track.position) && continue
        doc = XMLDocument()
        tags = ElementNode("tags")
        append_tag!(tags, "ALBUM", release.title)
        append_tag!(tags, "YEAR", release.year)
        append_tag!(tags, "TITLE", track.title)
        append_tag!(tags, "TRACKNUMBER", position)
        isnothing(release.genres) || append_tag!(tags, "GENRE", release.genres[1])
        if isnothing(track.artists)
            append_tag!(tags, "ARTIST", clean_artist_name(release.artists[1].name))
        else
            append_tag!(tags, "ARTIST", clean_artist_name(track.artists[1].name))
        end
        setroot!(doc, tags)
        name = "$(position) - $(track.title)"
        push!(docs, (name, doc))
        position += 1
    end

    return docs
end

function save_xml(release::Release, dir::AbstractString = "")
    docs = generate_xml(release::Release)
    for (name, doc) in docs
        filename = replace(name, r"[/\:*?\"<>|]" => "") * ".xml"
        write(joinpath(dir, filename), doc)
    end
end