Update the documentation setup

This sets everything up as built doc examples, splits tutorials from examples, adds more informative titles, and sets up buildkite

.buildkite/pipeline.yml

@@ -0,0 +1,18 @@
+steps:
+  - label: "Julia 1"
+    plugins:
+      - JuliaCI/julia#v1:
+          version: "1"
+      - JuliaCI/julia-test#v1:
+           coverage: false # 1000x slowdown
+    agents:
+      queue: "juliagpu"
+      cuda: "*"
+    env:
+      GROUP: 'GPU'
+    timeout_in_minutes: 60
+    # Don't run Buildkite if the commit message includes the text [skip tests]
+    if: build.message !~ /\[skip tests\]/
+
+env:
+  JULIA_PKG_SERVER: "" # it often struggles with our large artifacts

docs/pages.jl

@@ -1,13 +1,17 @@
 pages = [
     "Home" => "index.md",
     "Getting started" => "getting_started.md",
-    "Solver Algorithms" => ["MLP.md",
-        "DeepSplitting.md",
-        "DeepBSDE.md"],
     "Tutorials" => [
         "tutorials/deepsplitting.md",
         "tutorials/deepbsde.md",
         "tutorials/mlp.md",
     ],
+    "Extended Examples" => [
+        "examples/blackscholes",
+    ],
+    "Solver Algorithms" => 
+        ["MLP.md",
+        "DeepSplitting.md",
+        "DeepBSDE.md"],
     "Feynman Kac formula" => "Feynman_Kac.md",
 ]

docs/src/examples/blackscholes.md

@@ -0,0 +1,56 @@
+# Solving the 100-dimensional Black-Scholes-Barenblatt Equation
+
+Black Scholes equation is a model for stock option price.
+In 1973, Black and Scholes transformed their formula on option pricing and corporate liabilities into a PDE model, which is widely used in financing engineering for computing the option price over time. [1.]
+In this example, we will solve a Black-Scholes-Barenblatt equation of 100 dimensions.
+The Black-Scholes-Barenblatt equation is a nonlinear extension to the Black-Scholes
+equation, which models uncertain volatility and interest rates derived from the
+Black-Scholes equation. This model results in a nonlinear PDE whose dimension
+is the number of assets in the portfolio.
+
+To solve it using the `PIDEProblem`, we write:
+
+```julia
+d = 100 # number of dimensions
+X0 = repeat([1.0f0, 0.5f0], div(d,2)) # initial value of stochastic state
+tspan = (0.0f0,1.0f0)
+r = 0.05f0
+sigma = 0.4f0
+f(X,u,σᵀ∇u,p,t) = r * (u - sum(X.*σᵀ∇u))
+g(X) = sum(X.^2)
+μ_f(X,p,t) = zero(X) #Vector d x 1
+σ_f(X,p,t) = Diagonal(sigma*X) #Matrix d x d
+prob = PIDEProblem(g, f, μ_f, σ_f, X0, tspan)
+```
+
+As described in the API docs, we now need to define our `NNPDENS` algorithm
+by giving it the Flux.jl chains we want it to use for the neural networks.
+`u0` needs to be a `d`-dimensional -> 1-dimensional chain, while `σᵀ∇u`
+needs to be `d+1`-dimensional to `d` dimensions. Thus we define the following:
+
+```julia
+hls  = 10 + d #hide layer size
+opt = Flux.Optimise.Adam(0.001)
+u0 = Flux.Chain(Dense(d,hls,relu),
+                Dense(hls,hls,relu),
+                Dense(hls,1))
+σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu),
+                  Dense(hls,hls,relu),
+                  Dense(hls,hls,relu),
+                  Dense(hls,d))
+pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)
+```
+
+And now we solve the PDE. Here, we say we want to solve the underlying neural
+SDE using the Euler-Maruyama SDE solver with our chosen `dt=0.2`, do at most
+150 iterations of the optimizer, 100 SDE solves per loss evaluation (for averaging),
+and stop if the loss ever goes below `1f-6`.
+
+```julia
+ans = solve(prob, pdealg, verbose=true, maxiters=150, trajectories=100,
+                            alg=EM(), dt=0.2, pabstol = 1f-6)
+```
+
+## References
+
+1. Shinde, A. S., and K. C. Takale. "Study of Black-Scholes model and its applications." Procedia Engineering 38 (2012): 270-279.

docs/src/getting_started.md

@@ -117,7 +117,7 @@ sol = solve(prob,
 
 [`DeepSplitting`](@ref deepsplitting) can run on the GPU for (much) improved performance. To do so, just set `use_cuda = true`.
 
-```julia
+```@example DeepSplitting_gpu
 sol = solve(prob, 
             alg, 
             0.1, 

docs/src/tutorials/deepbsde.md

@@ -1,11 +1,9 @@
-# `DeepBSDE`
-<!-- TODO: revise the code, which should not work right here -->
-### Solving a 100-dimensional Hamilton-Jacobi-Bellman Equation
+# Solving a 100-dimensional Hamilton-Jacobi-Bellman Equation with `DeepBSDE`
 
 First, here's a fully working code for the solution of a 100-dimensional
 Hamilton-Jacobi-Bellman equation that takes a few minutes on a laptop:
 
-```julia
+```@example deepbsde
 using NeuralPDE
 using Flux, OptimizationOptimisers
 using DifferentialEquations
@@ -65,7 +63,7 @@ with terminating condition $g(x) = \log(1/2 + 1/2 \|x\|^2))$.
 
 To get the solution above using the `PIDEProblem`, we write:
 
-```julia
+```@example deepbsde2
 d = 100 # number of dimensions
 X0 = fill(0.0f0,d) # initial value of stochastic control process
 tspan = (0.0f0, 1.0f0)
@@ -85,7 +83,7 @@ by giving it the Flux.jl chains we want it to use for the neural networks.
 `u0` needs to be a `d` dimensional -> 1 dimensional chain, while `σᵀ∇u`
 needs to be `d+1` dimensional to `d` dimensions. Thus we define the following:
 
-```julia
+```@example deepbsde2
 hls = 10 + d #hidden layer size
 opt = Flux.Optimise.Adam(0.01)  #optimizer
 #sub-neural network approximating solutions at the desired point
@@ -102,7 +100,7 @@ pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)
 
 #### Solving with Neural Nets
 
-```julia
+```@example deepbsde2
 @time ans = solve(prob, pdealg, verbose=true, maxiters=100, trajectories=100,
                             alg=EM(), dt=0.2, pabstol = 1f-2)
 
@@ -111,61 +109,4 @@ pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)
 Here we want to solve the underlying neural
 SDE using the Euler-Maruyama SDE solver with our chosen `dt=0.2`, do at most
 100 iterations of the optimizer, 100 SDE solves per loss evaluation (for averaging),
-and stop if the loss ever goes below `1f-2`.
-
-## Solving the 100-dimensional Black-Scholes-Barenblatt Equation
-
-Black Scholes equation is a model for stock option price.
-In 1973, Black and Scholes transformed their formula on option pricing and corporate liabilities into a PDE model, which is widely used in financing engineering for computing the option price over time. [1.]
-In this example, we will solve a Black-Scholes-Barenblatt equation of 100 dimensions.
-The Black-Scholes-Barenblatt equation is a nonlinear extension to the Black-Scholes
-equation, which models uncertain volatility and interest rates derived from the
-Black-Scholes equation. This model results in a nonlinear PDE whose dimension
-is the number of assets in the portfolio.
-
-To solve it using the `PIDEProblem`, we write:
-
-```julia
-d = 100 # number of dimensions
-X0 = repeat([1.0f0, 0.5f0], div(d,2)) # initial value of stochastic state
-tspan = (0.0f0,1.0f0)
-r = 0.05f0
-sigma = 0.4f0
-f(X,u,σᵀ∇u,p,t) = r * (u - sum(X.*σᵀ∇u))
-g(X) = sum(X.^2)
-μ_f(X,p,t) = zero(X) #Vector d x 1
-σ_f(X,p,t) = Diagonal(sigma*X) #Matrix d x d
-prob = PIDEProblem(g, f, μ_f, σ_f, X0, tspan)
-```
-
-As described in the API docs, we now need to define our `NNPDENS` algorithm
-by giving it the Flux.jl chains we want it to use for the neural networks.
-`u0` needs to be a `d`-dimensional -> 1-dimensional chain, while `σᵀ∇u`
-needs to be `d+1`-dimensional to `d` dimensions. Thus we define the following:
-
-```julia
-hls  = 10 + d #hide layer size
-opt = Flux.Optimise.Adam(0.001)
-u0 = Flux.Chain(Dense(d,hls,relu),
-                Dense(hls,hls,relu),
-                Dense(hls,1))
-σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu),
-                  Dense(hls,hls,relu),
-                  Dense(hls,hls,relu),
-                  Dense(hls,d))
-pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)
-```
-
-And now we solve the PDE. Here, we say we want to solve the underlying neural
-SDE using the Euler-Maruyama SDE solver with our chosen `dt=0.2`, do at most
-150 iterations of the optimizer, 100 SDE solves per loss evaluation (for averaging),
-and stop if the loss ever goes below `1f-6`.
-
-```julia
-ans = solve(prob, pdealg, verbose=true, maxiters=150, trajectories=100,
-                            alg=EM(), dt=0.2, pabstol = 1f-6)
-```
-
-## References
-
-1. Shinde, A. S., and K. C. Takale. "Study of Black-Scholes model and its applications." Procedia Engineering 38 (2012): 270-279.
+and stop if the loss ever goes below `1f-2`.

docs/src/tutorials/deepsplitting.md

@@ -1,5 +1,5 @@
 
-# `DeepSplitting`
+# Solving the 10-dimensional Fisher-KPP equation with `DeepSplitting`
 Consider the Fisher-KPP equation with non-local competition
 
 ```math
@@ -10,7 +10,7 @@ where $\Omega = [-1/2, 1/2]^d$, and let's assume Neumann Boundary condition on $
 
 Let's solve Eq. (1) with the [`DeepSplitting`](@ref deepsplitting) solver. 
 
-```julia
+```@example deepsplitting
 using HighDimPDE
 
 ## Definition of the problem
@@ -50,10 +50,12 @@ sol = solve(prob,
             maxiters = 1000,
             batch_size = 1000)
 ```
+
 #### Solving on the GPU
+
 `DeepSplitting` can run on the GPU for (much) improved performance. To do so, just set `use_cuda = true`.
 
-```julia
+```@example deepsplitting2
 sol = solve(prob, 
             alg, 
             0.1, 

docs/src/tutorials/mlp.md

@@ -1,12 +1,12 @@
 
-# `MLP`
-## Local PDE
+# Solving the 10-dimensional Fisher-KPP equation with `MLP`
+
 Let's solve the [Fisher KPP](https://en.wikipedia.org/wiki/Fisher%27s_equation) PDE in dimension 10 with [`MLP`](@ref mlp).
 ```math
 \partial_t u = u (1 - u) + \frac{1}{2}\sigma^2\Delta_xu \tag{1}
 ```
 
-```julia
+```@example mlp
 using HighDimPDE
 
 ## Definition of the problem
@@ -32,7 +32,7 @@ Let's include in the previous equation non local competition, i.e.
 \partial_t u = u (1 - \int_\Omega u(t,y)dy) + \frac{1}{2}\sigma^2\Delta_xu \tag{2}
 ```
 where $\Omega = [-1/2, 1/2]^d$, and let's assume Neumann Boundary condition on $\Omega$.
-```julia
+```@example mlp2
 using HighDimPDE
 
 ## Definition of the problem