From 91488e0d3dddba95a9305a6c9a3d020b511bd4b3 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 18 Jan 2024 22:14:01 -0500 Subject: [PATCH 01/11] 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 | 18 ++++++++++ docs/pages.jl | 7 ++++ docs/src/examples/blackscholes.md | 56 +++++++++++++++++++++++++++++ docs/src/getting_started.md | 2 +- docs/src/tutorials/deepbsde.md | 14 ++++---- docs/src/tutorials/deepsplitting.md | 8 +++-- docs/src/tutorials/mlp.md | 8 ++--- 7 files changed, 97 insertions(+), 16 deletions(-) create mode 100644 .buildkite/pipeline.yml create mode 100644 docs/src/examples/blackscholes.md diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml new file mode 100644 index 0000000..6a675e4 --- /dev/null +++ b/.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 \ No newline at end of file diff --git a/docs/pages.jl b/docs/pages.jl index e733586..92d62ce 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -16,5 +16,12 @@ pages = [ "tutorials/nnkolmogorov.md", "tutorials/nnparamkolmogorov.md", ], + "Extended Examples" => [ + "examples/blackscholes", + ], + "Solver Algorithms" => + ["MLP.md", + "DeepSplitting.md", + "DeepBSDE.md"], "Feynman Kac formula" => "Feynman_Kac.md", ] diff --git a/docs/src/examples/blackscholes.md b/docs/src/examples/blackscholes.md new file mode 100644 index 0000000..59ff3f1 --- /dev/null +++ b/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. \ No newline at end of file diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index 8cbe510..dfaff3c 100644 --- a/docs/src/getting_started.md +++ b/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, diff --git a/docs/src/tutorials/deepbsde.md b/docs/src/tutorials/deepbsde.md index 76bac9f..1aadb8d 100644 --- a/docs/src/tutorials/deepbsde.md +++ b/docs/src/tutorials/deepbsde.md @@ -1,11 +1,9 @@ -# `DeepBSDE` - -### 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) @@ -168,4 +166,4 @@ ans = solve(prob, pdealg, verbose=true, maxiters=150, trajectories=100, ## References -1. Shinde, A. S., and K. C. Takale. "Study of Black-Scholes model and its applications." Procedia Engineering 38 (2012): 270-279. \ No newline at end of file +1. Shinde, A. S., and K. C. Takale. "Study of Black-Scholes model and its applications." Procedia Engineering 38 (2012): 270-279. diff --git a/docs/src/tutorials/deepsplitting.md b/docs/src/tutorials/deepsplitting.md index 60ac77c..11f47ea 100644 --- a/docs/src/tutorials/deepsplitting.md +++ b/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 @@ -48,10 +48,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, diff --git a/docs/src/tutorials/mlp.md b/docs/src/tutorials/mlp.md index 3774808..4756949 100644 --- a/docs/src/tutorials/mlp.md +++ b/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 From 3a9bca188f680b1266105e926f9d556d3cc8efd7 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 22 Jan 2024 05:22:32 -0500 Subject: [PATCH 02/11] Update docs/pages.jl --- docs/pages.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/pages.jl b/docs/pages.jl index 92d62ce..365b884 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -17,7 +17,7 @@ pages = [ "tutorials/nnparamkolmogorov.md", ], "Extended Examples" => [ - "examples/blackscholes", + "examples/blackscholes.md", ], "Solver Algorithms" => ["MLP.md", From 0be3a53c9784d2db3c9955343f191738b948f302 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 22 Jan 2024 06:19:56 -0500 Subject: [PATCH 03/11] Update docs/src/tutorials/deepbsde.md --- docs/src/tutorials/deepbsde.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/tutorials/deepbsde.md b/docs/src/tutorials/deepbsde.md index 1aadb8d..1e11fa5 100644 --- a/docs/src/tutorials/deepbsde.md +++ b/docs/src/tutorials/deepbsde.md @@ -4,7 +4,7 @@ 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: ```@example deepbsde -using NeuralPDE +using HighDimPDE using Flux, OptimizationOptimisers using DifferentialEquations using LinearAlgebra From c12fc24cd73242596dba430612acaf5305b740be Mon Sep 17 00:00:00 2001 From: Anant Thazhemadam Date: Tue, 23 Jan 2024 14:42:59 +0530 Subject: [PATCH 04/11] ci(buildkite): run documentation builds on buildkite --- .buildkite/pipeline.yml | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 6a675e4..ccde26e 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -14,5 +14,23 @@ steps: # Don't run Buildkite if the commit message includes the text [skip tests] if: build.message !~ /\[skip tests\]/ + - label: "Documentation" + plugins: + - JuliaCI/julia#v1: + version: "1" + command: | + julia --project=docs -e ' + println("--- :julia: Instantiating project") + using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate() + println("+++ :julia: Building documentation") + include("docs/make.jl")' + agents: + queue: "juliagpu" + cuda: "*" + if: build.message !~ /\[skip docs\]/ && !build.pull_request.draft + timeout_in_minutes: 1000 + env: + SECRET_DOCUMENTER_KEY: "v/2iYv+rW3iCgDU5Cp4DN1P51xQvXkC4aTO28pDLysKaNTVopjzsQNNtgJaswckg5ahdusdnAmlDRd8oMEobmyPM7PY7wppGdr6MHfUhwakJ7AAWAyAEIb+e+fNFuu+EIHS43bU+vkgzCgKojFEGwa5AFSW69rucnuwROXzIrQQ6GgXF4LayvDMDbUgYPjjnS2zAgfkWOR1L5UJ2ODuy3zK+VESS5YVFRNSmhAnT0EdS2AsenBwa25IgHsYkGO1fR/CL+1epFg7yzSOA+bExDWD0L4WkjKzU8qTlMAI2BAw+clJLhcpZyy4MdUG8hLmic0vvYjBOdjgqtl1f9xtMwQ==;U2FsdGVkX1+MBH0JHcFRDZ6tuHvLUrc8cxidXL6EZJZYrOllAtI7WHCEqx57m0T49q/ALas8vzZb9FJBGcjcrhCNu6jmbKNsBGNb7ZNRbex9qZCiZltv9g0cFU6xA7562CuYDhLxwIF4bXwszpYg/I+AhSCCIfaeaA1couJ0YTUxvKmW/VCR4/anH4SZD9HqblfMHa7N3FMo1VIEN91n+/xgyPNx3unQnre/UPDK0+Z7qjsYoXhqqQokJENdrsDWifjvdpWKh/1viGlMVzN0KFRYR0HIa8SN9wBuJAjVwgiz4SG7oGansuE1sowHon85iNsgIYr4k2WFpWY/f813EHeGvdFAu8/+/EJ42YvSHKJX2yOah1kw43OsBz/nBzQnffQNgG1/xz3WlqH5nXChlmfamzIaN4eERK2jQDIiPMTkbyPBNlDZ0TRHjzpWbpma2gMvnN2L0nyw6d7+yiop1DdYF/GTKuGXXR4I6dA1q5NdYvMQGmaXkDwqz2c6XLLGDPSbzfy8AOsqmAFfzCQTCVva96ynneBlpClEm2A2cW5l8N8/2Y/SfOguBKfFzjncPd+IK7BdHb0FFOhYmZVvcDglGp/Fy2snUe1JiSrtNNKA/CiRC5UQKcffTBP7qboSIOlamFKTR6GaNCYM3teLCodPiaV7mQuidCiFsR7WwORLRMpxkctMj08YZ+da8riUSrfXndXppZ5X1lNxSOQSmj4BaKVNkWb9h2pcoV3P6gD2PzvKiVsLRpqiLfIEuzfQxco2rU16DZj51zh28Gb8tJSbZWQB+pT2kOyVjBfI6AuWEJ3wkHWxyLntzXTUqc4WWnCu1wLPLWZWBA0rq28jh/jJE2NpGmYL5z8/+9T8i4RrSRwZWMbrSekiruSwrk12mWmbshVkArdSxR/Lk3AvvNetu4SccdDnD8CLxeTkdGYV1tU5OklWpxpPCjG8lI1oCCHCrEMTFu5rc/YeLbPcPCg4oKC5rqLpcJI/bmF+9fwqbhgULGfmHABqMZhs1fjvWEqkURCira0WsPHRVPrqRVqXaejpXbwWZbhdYHe9OoViBebj+6gm3l9KUthZBF+/acp5RvCk2KIH+GIhglVYdp9R3JZelRGTLZa1QGha2+RsnNhfE9Wp8ynjTQBknWUbCSCxOM+bjOTLyp1nZD+AjHtu56N30v8rGmBeWjae7t0e0lSxhdQmt5388R1TbeSxhPLGxTVW3X2WmAaIEnxzBWroTVELjchHq6CJIT7vF1UHFVmr4e3e8XH8o+9Mu3M36d5wAT9wT+ggWJTp+Oo48uZF+OBQ0erkNm37Qz307jPrLRO0UNsmuxawEGKATRgsJCWU1FE6YA+4CDXbcHBtj81flO+7MahO39/y1BLpHDt15sfhM94R8LLmjwsTQhk5kiiFHkxG1JT8I87fq3spO6hFXUFO04QSPYu6GeH6uca4COXV/lGbqv7ahTSZwKQR8EVoaiYHcAJpfK0gXXo9G4pe2xaEU6VP6TIB6wv6Vdag2JAOYt7v0pP08y3UiUMgfdPcdRwa/S7F68YT12gcKXgj/I9DAImyX01lfnMYNsE2w7vIPFgvDJ003uDiD+k5Dg6QPOBuEOryYpdaaDVpD8BKbKEMon7mxV3xJTQEXc84WXk+cSgR2a5fa8yyPn2AuZvc/ICyhxYNriMqd55b5ZP70jeBI0FJpZw00CvGroMuqUDCq7Kmf3tdKB8XHU7sJ/Q2OQEYdJkLSO9lFFu4iU856Xu4XWCwdFwH6CQh04/RFx8Sm18OCuczyHH/P4y930qQEqsjmjyJyU8+dQyGTZV08wPCcv23XAGOfBpnIJyaF44wW/rM6a+C+7dXZYFb1/viMpPxZg9C1ILl7gdo4l0QZxLsv/eVXthrDhSiG1AFQtG8QdS+JiQ2icxA9qpxrsdeCCxjhcZgkqV6MofMAnMr+mSZOS6FdVZcHskoi1M7Dq6iyVjNE9jrWnWWINiEwZ4PnulebI20xXsp27EcI9paHVmzfAKXs2+4P35GHPzGX7Foyx6rYkZkP/emhDjqInq2wUArkV21u7VL4AL4Rl01lAOyWdxfTGFVcrGDylXpKSZVTcJgmeCNdbq8uVHykjabgV+XytKWFf50m23bW5T7GZsBk+wuXe6vSrUfL3V5h29VbPorsMZV4p5a5Yn37aROD3pvb/Zex1WXSaxZyIn0z1hMyaElsiYvQ9T3jIwcjHrm/WEbcCPojNbJTGvXKAjcKNFuvMUAlpSNIWDcouKRD6ekOKlNhuJiDclLYKEj7Hj8n/btttruz1Z0G/Znpz7k6dzV7TlLU5bvmefX/yn0l5cMG70cb5111/Zv6twkWT7vXe2rX0DWzVst15DxplGID0VM4tDJO5lXTYNrvoVrK9rVCtAJtzlMLIF6zQxQWMxqhU7Sc7YC2nQJgl99igQFjqB1/BF2ZQZ+EZDFDuKSIrTsscN4xJTtjFSQkZWJT10RyIO0IamoBP2K6SVlRFG91+fC+vgyauqdvqBAWEQbQHsy3vyi3vxS3AeZFDLqeBuRU/L/aYiGde4ok9YXr6ab5d6Aj57bReDBt5r50GY1M/Aaz9gxLvfWBMp8EsZVHovItHL6BE8Ejk6C/nr2uDklktvbVPzE8YSh5BUfeH/PCymPTE6iCcu20tZOVxU1dhGgOkstRHJWgKi3RydyzyHihOsPZXpTsWzTmU/sc9rgbz4ypVRfk3s6hc2ClC2gxQEpVUMacaluvazXn8G24pvlhpEoX6rgjaILjkGiCTIuP4k9e3AQjIW6704Hj9o9j0E5GYpw0ed60ALkoCNoEW9X3f1DXN9x4lVefX2Ix522EOs6dKfbOJsJcI4C56MFlJkKbwETBt4YCumfqDafO98OKkSAHWQn+UvJraVE/64C5b7rJQOb00G9bxovDcyoGXZ28AyFA1abzvhWj7L9AFiOwe6sdouz7aQC2yTNZxNXbwZG5hJuMWv1TVM+TPxzQm/fvy60brAfPhE+vdJVbhzZ+LIwHVGwdfH8XlYJYcHj7sUM01ax0154Y6+V6/T8sJjDrXrHxUYhxJBJDDDsEDy6ZW891r/3mbdfk8faCxtFqD8YCJJaO7C43S/f4ELOtBq7VhFRLeNdLy4yfGcJLSVNDi2Mb26rsdxRIVo7ppekCQNNRSoNGWJUwmuiarbbti2+KdpJYib7KQH862VU1ki64GhlD6ZAO3S2PZ5RnPfBCZcsJ5fyMI8u637wk6kxoIGqsuH0A64gXp6qQo+z3FqRz2X35B9hrXTmshpewsUNMoFIz3WWijUbHDfmMyPAziX5ZDzEAX+C29CmaY6TACtN9KcYCQ+7/MQe8Ye8mHwFqCdNmpCNL9RWjweRykzPDO7M8nsCsnwsi2jOHy1u9C+KQaw+sN9yiJKxpITXWdJTwPwEKhTv0lL2AUW9h+ue77nh6eqrQy+p7FJELEWQePRKAf0ryk+GOX7uSiDuUqHdyPOhWGOlWdQ0IEurcNjVxC8NG4Md5wD2V+iQe49kCPeUDgjkm8YPt6HALSLIdRe7EfpJ9QAe8Z2yCYZuk0ckHcyHcPotSvCyMAY4zmke/GjjbBR3jpcFmGZT0PyD76+yEy5te86dH1wVT7gPozKAwH4Wt3l9xzBYEDx5WSGnC/AEe8bkHAQSUFdeFZIvt7/poonr5/DNBP33z2DEw675ym+Jx9wtYrgQK4HUy9mSf7BYGmt0hWy8s4+t8lOohxcIYytShJXybJA1L4904PA13pd5VJgDFGB2xrHQ/UywsAIGFYsv9fEkkH83pvRDAufmxcrB5DTf0cCul3qv7gI38Gh+FT+By98k9ucBCAzH0APthuM7ERrno0oEnASXeGNAkN/1vVcVSors6CJjFdB5I3+e2yub4ZLa99qXoyeTr4aR9oIuMcvcRagqpVzAwxirHu9mNlMQ1MNt41BbLPn36gfx5jWoUxZCwtNIvKXUfz8kAnjryUL0qCjrr1ZvFp0CC0P6tprpUiLfp2OUbz36PvGnbK3SaDkj05X/BUnLtRObl2o5YOdXZCuzTCGYP3GSbzI3ot3ps8N6RxhMCsN8xyn6yzojLuzm2hWhkgLH626KcZt0fvxWMUinJKIShnjJveWcX1FPQpYRv2k/EwF5lidmBI51DDtU+N9c34FMA3bgYyN2LwNP6HesAZAEtQ0GHYDXPJzjS2t01gVnb3ei7Gdm6GY4Tc0XimM/IIf68qfsESwMYGG2X8siGtM/kFJSxGXwbIAmwLq3wYO4TYvL0ZD4z0ObMpIOiQfmVJngitZgsFCrImrMpDV4nhGePS7nlu2SkLnTKN1CyQvoCrwfwSKGeqNtGsFeUY26zTS+c3h6X0pudKSt2zfIEl/yJBTfotKYtDT+GnYdXxsCo/RixxOTS7HpfNvarCFLnU6p" + env: JULIA_PKG_SERVER: "" # it often struggles with our large artifacts \ No newline at end of file From a5bc5d7188f922d0770559a5504b7b9f00cbffe5 Mon Sep 17 00:00:00 2001 From: Ashutosh Bharambe Date: Sat, 22 Jun 2024 20:37:32 +0530 Subject: [PATCH 05/11] docs: fix docs --- docs/Project.toml | 2 + docs/make.jl | 2 +- docs/src/NNKolmogorov.md | 2 +- docs/src/NNParamKolmogorov.md | 2 +- docs/src/assets/Project.toml | 11 ++++ docs/src/getting_started.md | 6 +- docs/src/index.md | 13 +++-- docs/src/tutorials/deepbsde.md | 76 +++++++++++++++---------- docs/src/tutorials/deepsplitting.md | 57 +++++++++++++++++-- docs/src/tutorials/mlp.md | 4 +- docs/src/tutorials/nnkolmogorov.md | 8 ++- docs/src/tutorials/nnparamkolmogorov.md | 7 ++- docs/src/tutorials/nnstopping.md | 16 +++--- src/DeepSplitting.jl | 2 +- 14 files changed, 147 insertions(+), 61 deletions(-) create mode 100644 docs/src/assets/Project.toml diff --git a/docs/Project.toml b/docs/Project.toml index 45b25ee..f02b906 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,6 +3,8 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" HighDimPDE = "57c578d5-59d4-4db8-a490-a9fc372d19d2" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" [compat] CUDA = "5" diff --git a/docs/make.jl b/docs/make.jl index 896ce68..ea0bd67 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -6,7 +6,7 @@ cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true) include("pages.jl") makedocs(sitename = "HighDimPDE.jl", - authors = "Victor Boussange", + authors = "#", pages = pages, clean = true, doctest = false, linkcheck = true, format = Documenter.HTML(assets = ["assets/favicon.ico"], diff --git a/docs/src/NNKolmogorov.md b/docs/src/NNKolmogorov.md index 73d554b..669f01a 100644 --- a/docs/src/NNKolmogorov.md +++ b/docs/src/NNKolmogorov.md @@ -1,4 +1,4 @@ -# [The `NNKolmogorov` algorithm](@id nn_komogorov) +# [The `NNKolmogorov` algorithm](@id nn_kolmogorov) ### Problems Supported: 1. [`ParabolicPDEProblem`](@ref) diff --git a/docs/src/NNParamKolmogorov.md b/docs/src/NNParamKolmogorov.md index 5635f8c..0b36d93 100644 --- a/docs/src/NNParamKolmogorov.md +++ b/docs/src/NNParamKolmogorov.md @@ -1,4 +1,4 @@ -# [The `NNParamKolmogorov` algorithm](@id nn_komogorov) +# [The `NNParamKolmogorov` algorithm](@id nn_paramkolmogorov) ### Problems Supported: 1. [`ParabolicPDEProblem`](@ref) diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml new file mode 100644 index 0000000..4ef18fb --- /dev/null +++ b/docs/src/assets/Project.toml @@ -0,0 +1,11 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" +HighDimPDE = "57c578d5-59d4-4db8-a490-a9fc372d19d2" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" + +[compat] +Documenter = "1" +Flux = "0.13, 0.14" +HighDimPDE = "2" diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index dfaff3c..f7da87f 100644 --- a/docs/src/getting_started.md +++ b/docs/src/getting_started.md @@ -32,8 +32,8 @@ x0 = fill(0., d) # initial point g(x) = exp(-sum(x.^2)) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients -f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the -prob = PIDEProblem(μ, σ, x0, tspan, g, f) # defining the problem +f(x, v_x, ∇v_x, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the +prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f) # defining the problem ## Definition of the algorithm alg = MLP() # defining the algorithm. We use the Multi Level Picard algorithm @@ -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`. -```@example DeepSplitting_gpu +```@example DeepSplitting_non_local_PDE sol = solve(prob, alg, 0.1, diff --git a/docs/src/index.md b/docs/src/index.md index 32c70e0..7c50a08 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -3,7 +3,7 @@ **HighDimPDE.jl** is a Julia package to **solve Highly Dimensional non-linear, non-local PDEs** of the forms: -1. Partial Integro Differential Equations: +### 1. Partial Integro Differential Equations: ```math \begin{aligned} (\partial_t u)(t,x) &= \int_{\Omega} f\big(t,x,{\bf x}, u(t,x),u(t,{\bf x}), ( \nabla_x u )(t,x ),( \nabla_x u )(t,{\bf x} ) \big) \, d{\bf x} \\ @@ -11,9 +11,9 @@ \end{aligned} ``` -where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject to initial and boundary conditions, and where $d$ is large. +where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject to initial and boundary conditions, and where $d$ is large. These equations are defined using the [`PIDEProblem`](@ref) -2. Parabolic Partial Differential Equations: +### 2. Parabolic Partial Differential Equations: ```math \begin{aligned} (\partial_t u)(t,x) &= f\big(t,x, u(t,x), ( \nabla_x u )(t,x )\big) @@ -21,7 +21,7 @@ where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject \end{aligned} ``` -where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject to initial and boundary conditions, and where $d$ is large. +where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject to initial and boundary conditions, and where $d$ is large. These equations are defined using the [`ParabolicPDEProblem`](@ref) !!! note The difference between the two problems is that in Partial Integro Differential Equations, the integrand is integrated over **x**, while in Parabolic Integro Differential Equations, the function `f` is just evaluated for `x`. @@ -32,8 +32,11 @@ where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject * the [Multi-Level Picard iterations scheme](@ref mlp) -* the Deep BSDE scheme (@ref deepbsde). +* [the Deep BSDE scheme](@ref deepbsde). +* [NNKolmogorov](@ref nn_kolmogorov) and [NNParamKolmogorov](@ref nn_paramkolmogorov) schemes for Kolmogorov PDEs. + +* [NNStopping] scheme for solving optimal stopping time problem. To make the most out of **HighDimPDE.jl**, we advise to first have a look at the diff --git a/docs/src/tutorials/deepbsde.md b/docs/src/tutorials/deepbsde.md index 1e11fa5..f439e18 100644 --- a/docs/src/tutorials/deepbsde.md +++ b/docs/src/tutorials/deepbsde.md @@ -5,33 +5,46 @@ Hamilton-Jacobi-Bellman equation that takes a few minutes on a laptop: ```@example deepbsde using HighDimPDE -using Flux, OptimizationOptimisers -using DifferentialEquations +using Flux +using StochasticDiffEq using LinearAlgebra + d = 100 # number of dimensions -X0 = fill(0.0f0, d) # initial value of stochastic control process +x0 = fill(0.0f0, d) tspan = (0.0f0, 1.0f0) +dt = 0.2f0 λ = 1.0f0 +# +g(X) = log(0.5f0 + 0.5f0 * sum(X .^ 2)) +f(X, u, σᵀ∇u, p, t) = -λ * sum(σᵀ∇u .^ 2) +μ_f(X, p, t) = zero(X) #Vector d x 1 λ +σ_f(X, p, t) = Diagonal(sqrt(2.0f0) * ones(Float32, d)) #Matrix d x d +prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan; g, f) -g(X) = log(0.5f0 + 0.5f0 * sum(X.^2)) -f(X,u,σᵀ∇u,p,t) = -λ * sum(σᵀ∇u.^2) -μ_f(X,p,t) = zero(X) # Vector d x 1 λ -σ_f(X,p,t) = Diagonal(sqrt(2.0f0) * ones(Float32, d)) # Matrix d x d -prob = PIDEProblem(μ_f, σ_f, X0, tspan, g, f) -hls = 10 + d # hidden layer size -opt = Optimisers.Adam(0.01) # optimizer -# sub-neural network approximating solutions at the desired point +hls = 10 + d #hidden layer size +opt = Flux.Optimise.Adam(0.1) #optimizer +#sub-neural network approximating solutions at the desired point u0 = Flux.Chain(Dense(d, hls, relu), - Dense(hls, hls, relu), - Dense(hls, 1)) + Dense(hls, hls, relu), + Dense(hls, hls, relu), + Dense(hls, 1)) # sub-neural network approximating the spatial gradients at time point σᵀ∇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) -@time ans = solve(prob, pdealg, verbose=true, maxiters=100, trajectories=100, - alg=EM(), dt=1.2, pabstol=1f-2) + Dense(hls, hls, relu), + Dense(hls, hls, relu), + Dense(hls, hls, relu), + Dense(hls, d)) +pdealg = DeepBSDE(u0, σᵀ∇u, opt = opt) + +@time sol = solve(prob, + pdealg, + StochasticDiffEq.EM(), + verbose = true, + maxiters = 150, + trajectories = 30, + dt = 1.2f0, + pabstol = 1.0f-4) + ``` Now, let's explain the details! @@ -61,9 +74,14 @@ with terminating condition $g(x) = \log(1/2 + 1/2 \|x\|^2))$. #### Define the Problem -To get the solution above using the `PIDEProblem`, we write: +To get the solution above using the [`ParabolicPDEProblem`](@ref), we write: ```@example deepbsde2 +using HighDimPDE +using Flux +using StochasticDiffEq +using LinearAlgebra + d = 100 # number of dimensions X0 = fill(0.0f0,d) # initial value of stochastic control process tspan = (0.0f0, 1.0f0) @@ -73,12 +91,12 @@ g(X) = log(0.5f0 + 0.5f0*sum(X.^2)) f(X,u,σᵀ∇u,p,t) = -λ*sum(σᵀ∇u.^2) μ_f(X,p,t) = zero(X) #Vector d x 1 λ σ_f(X,p,t) = Diagonal(sqrt(2.0f0)*ones(Float32,d)) #Matrix d x d -prob = PIDEProblem(μ_f, σ_f, X0, tspan, g, f) +prob = ParabolicPDEProblem(μ_f, σ_f, X0, tspan; g, f) ``` #### Define the Solver Algorithm -As described in the API docs, we now need to define our `NNPDENS` algorithm +As described in the API docs, we now need to define our `DeepBSDE` 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: @@ -95,14 +113,13 @@ u0 = Flux.Chain(Dense(d,hls,relu), Dense(hls,hls,relu), Dense(hls,hls,relu), Dense(hls,d)) -pdealg = NNPDENS(u0, σᵀ∇u, opt=opt) +pdealg = DeepBSDE(u0, σᵀ∇u, opt = opt) ``` #### Solving with Neural Nets ```@example deepbsde2 -@time ans = solve(prob, pdealg, verbose=true, maxiters=100, trajectories=100, - alg=EM(), dt=0.2, pabstol = 1f-2) +@time ans = solve(prob, pdealg, EM(), verbose=true, maxiters=100, trajectories=100, dt=0.2f0, pabstol = 1f-2) ``` @@ -121,7 +138,7 @@ 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: +To solve it using the `ParabolicPDEProblem`, we write: ```julia d = 100 # number of dimensions @@ -133,7 +150,7 @@ 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(μ_f, σ_f, X0, tspan, g, f) +prob = ParabolicPDEProblem(μ_f, σ_f, X0, tspan; g, f) ``` As described in the API docs, we now need to define our `NNPDENS` algorithm @@ -151,7 +168,7 @@ u0 = Flux.Chain(Dense(d,hls,relu), Dense(hls,hls,relu), Dense(hls,hls,relu), Dense(hls,d)) -pdealg = NNPDENS(u0, σᵀ∇u, opt=opt) +pdealg = DeepBSDE(u0, σᵀ∇u, opt=opt) ``` And now we solve the PDE. Here, we say we want to solve the underlying neural @@ -160,8 +177,7 @@ SDE using the Euler-Maruyama SDE solver with our chosen `dt=0.2`, do at most 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) +ans = solve(prob, pdealg, EM(), verbose=true, maxiters=150, trajectories=100, dt=0.2f0) ``` ## References diff --git a/docs/src/tutorials/deepsplitting.md b/docs/src/tutorials/deepsplitting.md index 11f47ea..138c7c9 100644 --- a/docs/src/tutorials/deepsplitting.md +++ b/docs/src/tutorials/deepsplitting.md @@ -11,7 +11,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. ```@example deepsplitting -using HighDimPDE +using HighDimPDE, Flux ## Definition of the problem d = 10 # dimension of the problem @@ -22,10 +22,13 @@ g(x) = exp.(- sum(x.^2, dims=1) ) # initial condition σ(x, p, t) = 0.1f0 # diffusion coefficients x0_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = v_x .* (1f0 .- v_y) +``` +Since this is a non-local equation, we will define our problem as a [`PIDEProblem`](@ref) +```@example deepsplitting prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample) - +``` +```@example deepsplitting ## Definition of the neural network to use -using Flux # needed to define the neural network hls = d + 50 #hidden layer size @@ -53,7 +56,7 @@ sol = solve(prob, `DeepSplitting` can run on the GPU for (much) improved performance. To do so, just set `use_cuda = true`. -```@example deepsplitting2 +```@example deepsplitting sol = solve(prob, alg, 0.1, @@ -62,4 +65,50 @@ sol = solve(prob, maxiters = 1000, batch_size = 1000, use_cuda=true) +``` + +# Solving local Allen Cahn Equation with Neumann BC +```@example deepsplitting2 +using HighDimPDE, Flux, StochasticDiffEq +batch_size = 1000 +train_steps = 1000 + +tspan = (0.0f0, 5.0f-1) +dt = 5.0f-2 # time step + +μ(x, p, t) = 0.0f0 # advection coefficients +σ(x, p, t) = 1.0f0 # diffusion coefficients + +d = 10 +∂ = fill(5.0f-1, d) + +hls = d + 50 #hidden layer size + +nn = Flux.Chain(Dense(d, hls, relu), + Dense(hls, hls, relu), + Dense(hls, 1)) # Neural network used by the scheme + +opt = Flux.Optimise.Adam(1e-2) #optimiser +alg = DeepSplitting(nn, opt = opt) + +X0 = fill(0.0f0, d) # initial point +g(X) = exp.(-0.25f0 * sum(X .^ 2, dims = 1)) # initial condition +a(u) = u - u^3 +f(y, v_y, ∇v_y, p, t) = a.(v_y) # nonlocal nonlinear function +``` +Since we are dealing with a local problem here, i.e. no integration term used, we use [`ParabolicPDEProblem`](@ref) to define the problem. +```@example deepsplitting2 +# defining the problem +prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f, neumann_bc = [-∂, ∂]) +``` +```@example deepsplitting2 +# solving +@time sol = solve(prob, + alg, + dt, + verbose = false, + abstol = 1e-5, + use_cuda = false, + maxiters = train_steps, + batch_size = batch_size) ``` \ No newline at end of file diff --git a/docs/src/tutorials/mlp.md b/docs/src/tutorials/mlp.md index 4756949..15ecf69 100644 --- a/docs/src/tutorials/mlp.md +++ b/docs/src/tutorials/mlp.md @@ -17,7 +17,7 @@ g(x) = exp(- sum(x.^2) ) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients f(x, v_x, ∇v_x, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the -prob = ParabolicPDEProblem(μ, σ, x0, tspan, g, f) # defining the problem +prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f) # defining the problem ## Definition of the algorithm alg = MLP() # defining the algorithm. We use the Multi Level Picard algorithm @@ -43,7 +43,7 @@ g(x) = exp( -sum(x.^2) ) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients mc_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) -f(x, y, v_x, v_y, ∇v_x, ∇v_y, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) +f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) prob = PIDEProblem(μ, σ, x0, tspan, g, f) # defining x0_sample is sufficient to implement Neumann boundary conditions ## Definition of the algorithm diff --git a/docs/src/tutorials/nnkolmogorov.md b/docs/src/tutorials/nnkolmogorov.md index 9f35f95..29fd355 100644 --- a/docs/src/tutorials/nnkolmogorov.md +++ b/docs/src/tutorials/nnkolmogorov.md @@ -2,7 +2,9 @@ ## Solving high dimensional Rainbow European Options for a range of initial stock prices: -```julia +```@example nnkolmogorov +using HighDimPDE, StochasticDiffEq, Flux, LinearAlgebra + d = 10 # dims T = 1/12 sigma = 0.01 .+ 0.03.*Matrix(Diagonal(ones(d))) # volatility @@ -22,12 +24,12 @@ xspan = [(98.00, 102.00) for i in 1:d] g(x) = max(maximum(x) -K, 0) -sdealg = EM() +sdealg = SRIW1() # provide `x0` as nothing to the problem since we are provinding a range for `x0`. prob = ParabolicPDEProblem(μ_func, σ_func, nothing, tspan, g = g, xspan = xspan) opt = Flux.Optimisers.Adam(0.01) -alg = NNKolmogorov(m, opt) m = Chain(Dense(d, 16, elu), Dense(16, 32, elu), Dense(32, 16, elu), Dense(16, 1)) +alg = NNKolmogorov(m, opt) sol = solve(prob, alg, sdealg, verbose = true, dt = 0.01, dx = 0.0001, trajectories = 1000, abstol = 1e-6, maxiters = 300) ``` diff --git a/docs/src/tutorials/nnparamkolmogorov.md b/docs/src/tutorials/nnparamkolmogorov.md index 735cf2d..f9de711 100644 --- a/docs/src/tutorials/nnparamkolmogorov.md +++ b/docs/src/tutorials/nnparamkolmogorov.md @@ -3,7 +3,8 @@ ## Solving Parametric Family of High Dimensional Heat Equation. In this example we will solve the high dimensional heat equation over a range of initial values, and also over a range of thermal diffusivity. -```julia +```@example nnparamkolmogorov +using HighDimPDE, Flux, StochasticDiffEq d = 10 # models input is `d` for initial values, `d` for thermal diffusivity, and last dimension is for stopping time. m = Chain(Dense(d + 1 + 1, 32, relu), Dense(32, 16, relu), Dense(16, 8, relu), Dense(8, 1)) @@ -49,13 +50,13 @@ sol = solve(prob, NNParamKolmogorov(m, opt), sdealg, verbose = true, dt = 0.01, Similarly we can parametrize the drift function `mu_` and the initial function `g`, and obtain a solution over all parameters and initial values. # Inferring on the solution from `NNParamKolmogorov`: -```julia +```@example nnparamkolmogorov x_test = rand(xspan[1][1]:0.1:xspan[1][2], d) p_sigma_test = rand(p_domain.p_sigma[1]:dps.p_sigma:p_domain.p_sigma[2], 1, 1) t_test = rand(tspan[1]:dt:tspan[2], 1, 1) p_mu_test = nothing p_phi_test = nothing ``` -```julia +```@example nnparamkolmogorov sol.ufuns(x_test, t_test, p_sigma_test, p_mu_test, p_phi_test) ``` \ No newline at end of file diff --git a/docs/src/tutorials/nnstopping.md b/docs/src/tutorials/nnstopping.md index 25edad2..486fa8e 100644 --- a/docs/src/tutorials/nnstopping.md +++ b/docs/src/tutorials/nnstopping.md @@ -8,11 +8,12 @@ We will calculate optimal strategy for Bermudan Max-Call option with following d g(t, x) = e^{-rt}max\lbrace max\lbrace x1, ... , xd \rbrace − K, 0\rbrace ``` We define the parameters, drift function and the diffusion function for the dynamics of the option. -```julia +```@example nnstopping +using HighDimPDE, Flux, StochasticDiffEq d = 3 # Number of assets in the stock r = 0.05 # interest rate beta = 0.2 # volatility -T = 3 # maturity +T = 3.0 # maturity u0 = fill(90.0, d) # initial stock value delta = 0.1 # delta f(du, u, p, t) = du .= (r - delta) * u # drift @@ -28,15 +29,16 @@ function g(x, t) end ``` -We then define a `PIDEProblem` with no non linear term: -```julia -prob = PIDEProblem(f, sigma, u0, tspan; payoff = g) +We then define a [`ParabolicPDEProblem`](@ref) with no non linear term: + +```@example nnstopping +prob = ParabolicPDEProblem(f, sigma, u0, tspan; payoff = g) ``` !!! note We provide the payoff function with a keyword argument `payoff` And now we define our models: -```julia +```@example nnstopping models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) for i in 1:N] ``` @@ -44,7 +46,7 @@ models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoi The number of models should be equal to the time discritization. And finally we define our optimizer and algorithm, and call `solve`: -```julia +```@example nnstopping opt = Flux.Optimisers.Adam(0.01) alg = NNStopping(models, opt) diff --git a/src/DeepSplitting.jl b/src/DeepSplitting.jl index 91ca483..0722ee9 100644 --- a/src/DeepSplitting.jl +++ b/src/DeepSplitting.jl @@ -10,7 +10,7 @@ end Deep splitting algorithm. # Arguments -* `nn`: a [Flux.Chain](https://fluxml.ai/Flux.jl/stable/models/layers/#Flux.Chain), or more generally a [functor](https://github.com/FluxML/Functors.jl). +* `nn`: a [Flux.Chain](https://fluxml.ai/Flux.jl/stable/reference/models/layers/#Flux.Chain), or more generally a [functor](https://github.com/FluxML/Functors.jl). * `K`: the number of Monte Carlo integrations. * `opt`: optimizer to be used. By default, `Flux.Optimise.Adam(0.01)`. * `λs`: the learning rates, used sequentially. Defaults to a single value taken from `opt`. From 98c48d79276968f95d1a942952bc70a6ff9384c6 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Tue, 2 Jul 2024 14:52:35 +0530 Subject: [PATCH 06/11] chore: fix formatting --- LICENSE.md | 1 - README.md | 15 ++-- docs/pages.jl | 11 ++- docs/src/DeepBSDE.md | 9 +- docs/src/DeepSplitting.md | 63 ++++++++------ docs/src/Feynman_Kac.md | 43 +++++++--- docs/src/MLP.md | 52 +++++++----- docs/src/NNKolmogorov.md | 22 +++-- docs/src/NNParamKolmogorov.md | 22 +++-- docs/src/NNStopping.md | 6 +- docs/src/examples/blackscholes.md | 36 ++++---- docs/src/getting_started.md | 76 ++++++++--------- docs/src/index.md | 42 ++++++---- docs/src/problems.md | 6 +- docs/src/tutorials/deepbsde.md | 65 +++++++-------- docs/src/tutorials/deepsplitting.md | 81 +++++++++--------- docs/src/tutorials/mlp.md | 30 ++++--- docs/src/tutorials/nnkolmogorov.md | 8 +- docs/src/tutorials/nnparamkolmogorov.md | 10 ++- docs/src/tutorials/nnstopping.md | 20 +++-- paper/paper.md | 105 ++++++++++++++---------- src/DeepBSDE.jl | 10 +-- src/DeepBSDE_Han.jl | 4 +- src/DeepSplitting.jl | 8 +- src/HighDimPDE.jl | 92 ++++++++++----------- src/MLP.jl | 8 +- test/DeepBSDE.jl | 4 +- test/DeepBSDE_Han.jl | 32 ++++---- test/DeepSplitting.jl | 4 +- test/MLP.jl | 2 +- test/NNParamKolmogorov.jl | 3 +- 31 files changed, 510 insertions(+), 380 deletions(-) diff --git a/LICENSE.md b/LICENSE.md index 28e58b6..b480e81 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -19,4 +19,3 @@ The HighDimPDE.jl package is licensed under the MIT "Expat" License: > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE > SOFTWARE. -> \ No newline at end of file diff --git a/README.md b/README.md index 3929ab7..ea53dbb 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ [![codecov](https://codecov.io/gh/SciML/HighDimPDE.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/SciML/HighDimPDE.jl) [![Build Status](https://github.com/SciML/HighDimPDE.jl/workflows/CI/badge.svg)](https://github.com/SciML/HighDimPDE.jl/actions?query=workflow%3ACI) -[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac) +[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor%27s%20Guide-blueviolet)](https://github.com/SciML/ColPrac) [![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle) # HighDimPDE.jl @@ -13,12 +13,13 @@ $$ \begin{aligned} - (\partial_t u)(t,x) &= \int_{\Omega} f\big(t,x,{\bf x}, u(t,x),u(t,{\bf x}), ( \nabla_x u )(t,x ),( \nabla_x u )(t,{\bf x} ) \big) d{\bf x} \\ - & \quad + \big\langle \mu(t,x), ( \nabla_x u )( t,x ) \big\rangle + \tfrac{1}{2} \text{Trace} \big(\sigma(t,x) [ \sigma(t,x) ]^* ( \text{Hess}_x u)(t, x ) \big). +(\partial_t u)(t,x) &= \int_{\Omega} f\big(t,x,{\bf x}, u(t,x),u(t,{\bf x}), ( \nabla_x u )(t,x ),( \nabla_x u )(t,{\bf x} ) \big) d{\bf x} \\ +& \quad + \big\langle \mu(t,x), ( \nabla_x u )( t,x ) \big\rangle + \tfrac{1}{2} \text{Trace} \big(\sigma(t,x) [ \sigma(t,x) ]^* ( \text{Hess}_x u)(t, x ) \big). \end{aligned} $$ where $u \colon [0,T] \times \Omega \to \mathbb{R}, \Omega \subseteq \mathbb{R}^{d}$ is subject to initial and boundary conditions, and where $d$ is large. + ## Tutorials and Documentation For information on using the package, @@ -27,22 +28,26 @@ For information on using the package, the documentation, which contains the unreleased features. ## Installation + Open Julia and type the following ```julia using Pkg; Pkg.add("HighDimPDE.jl") ``` + This will download the latest version from the git repo and download all dependencies. ## Getting started + See documentation and `test` folders. ## Reference -- Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. [arXiv](https://arxiv.org/abs/2205.03672) (2022) + + - Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. [arXiv](https://arxiv.org/abs/2205.03672) (2022) + - diff --git a/docs/pages.jl b/docs/pages.jl index 365b884..9200e33 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -5,7 +5,7 @@ pages = [ "Solver Algorithms" => ["MLP.md", "DeepSplitting.md", "DeepBSDE.md", - "NNStopping.md", + "NNStopping.md", "NNKolmogorov.md", "NNParamKolmogorov.md"], "Tutorials" => [ @@ -14,14 +14,13 @@ pages = [ "tutorials/mlp.md", "tutorials/nnstopping.md", "tutorials/nnkolmogorov.md", - "tutorials/nnparamkolmogorov.md", + "tutorials/nnparamkolmogorov.md" ], "Extended Examples" => [ - "examples/blackscholes.md", + "examples/blackscholes.md" ], - "Solver Algorithms" => - ["MLP.md", + "Solver Algorithms" => ["MLP.md", "DeepSplitting.md", "DeepBSDE.md"], - "Feynman Kac formula" => "Feynman_Kac.md", + "Feynman Kac formula" => "Feynman_Kac.md" ] diff --git a/docs/src/DeepBSDE.md b/docs/src/DeepBSDE.md index 428f4c3..41f99e2 100644 --- a/docs/src/DeepBSDE.md +++ b/docs/src/DeepBSDE.md @@ -1,7 +1,8 @@ # [The `DeepBSDE` algorithm](@id deepbsde) ### Problems Supported: -1. [`ParabolicPDEProblem`](@ref) + + 1. [`ParabolicPDEProblem`](@ref) ```@autodocs Modules = [HighDimPDE] @@ -9,8 +10,10 @@ Pages = ["DeepBSDE.jl", "DeepBSDE_Han.jl"] ``` ## The general idea 💡 -The `DeepBSDE` algorithm is similar in essence to the `DeepSplitting` algorithm, with the difference that + +The `DeepBSDE` algorithm is similar in essence to the `DeepSplitting` algorithm, with the difference that it uses two neural networks to approximate both the the solution and its gradient. ## References -- Han, J., Jentzen, A., E, W., Solving high-dimensional partial differential equations using deep learning. [arXiv](https://arxiv.org/abs/1707.02568) (2018) \ No newline at end of file + + - Han, J., Jentzen, A., E, W., Solving high-dimensional partial differential equations using deep learning. [arXiv](https://arxiv.org/abs/1707.02568) (2018) diff --git a/docs/src/DeepSplitting.md b/docs/src/DeepSplitting.md index 2639036..cf0d2bd 100644 --- a/docs/src/DeepSplitting.md +++ b/docs/src/DeepSplitting.md @@ -1,8 +1,9 @@ # [The `DeepSplitting` algorithm](@id deepsplitting) ### Problems Supported: -1. [`PIDEProblem`](@ref) -2. [`ParabolicPDEProblem`](@ref) + + 1. [`PIDEProblem`](@ref) + 2. [`ParabolicPDEProblem`](@ref) ```@autodocs Modules = [HighDimPDE] @@ -13,76 +14,93 @@ The `DeepSplitting` algorithm reformulates the PDE as a stochastic learning prob The algorithm relies on two main ideas: -- The approximation of the solution $u$ by a parametric function $\bf u^\theta$. This function is generally chosen as a (Feedforward) Neural Network, as it is a [universal approximator](https://en.wikipedia.org/wiki/Universal_approximation_theorem). + - The approximation of the solution $u$ by a parametric function $\bf u^\theta$. This function is generally chosen as a (Feedforward) Neural Network, as it is a [universal approximator](https://en.wikipedia.org/wiki/Universal_approximation_theorem). -- The training of $\bf u^\theta$ by simulated stochastic trajectories of particles, through the link between linear PDEs and the expected trajectory of associated Stochastic Differential Equations (SDEs), explicitly stated by the [Feynman Kac formula](https://en.wikipedia.org/wiki/Feynman–Kac_formula). + - The training of $\bf u^\theta$ by simulated stochastic trajectories of particles, through the link between linear PDEs and the expected trajectory of associated Stochastic Differential Equations (SDEs), explicitly stated by the [Feynman Kac formula](https://en.wikipedia.org/wiki/Feynman%E2%80%93Kac_formula). ## The general idea 💡 + Consider the PDE + ```math \partial_t u(t,x) = \mu(t, x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) + f(x, u(t,x)) \tag{1} ``` -with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$. + +with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$. ### Local Feynman Kac formula + `DeepSplitting` solves the PDE iteratively over small time intervals by using an approximate [Feynman-Kac representation](@ref feynmankac) locally. More specifically, considering a small time step $dt = t_{n+1} - t_n$ one has that + ```math u(t_{n+1}, X_{T - t_{n+1}}) \approx \mathbb{E} \left[ f(t, X_{T - t_{n}}, u(t_{n},X_{T - t_{n}}))(t_{n+1} - t_n) + u(t_{n}, X_{T - t_{n}}) | X_{T - t_{n+1}}\right] \tag{3}. ``` + One can therefore use Monte Carlo integrations to approximate the expectations + ```math u(t_{n+1}, X_{T - t_{n+1}}) \approx \frac{1}{\text{batch\_size}}\sum_{j=1}^{\text{batch\_size}} \left[ u(t_{n}, X_{T - t_{n}}^{(j)}) + (t_{n+1} - t_n)\sum_{k=1}^{K} \big[ f(t_n, X_{T - t_{n}}^{(j)}, u(t_{n},X_{T - t_{n}}^{(j)})) \big] \right] ``` - ### Reformulation as a learning problem + The `DeepSplitting` algorithm approximates $u(t_{n+1}, x)$ by a parametric function ${\bf u}^\theta_n(x)$. It is advised to let this function be a neural network ${\bf u}_\theta \equiv NN_\theta$ as they are universal approximators. -For each time step $t_n$, the `DeepSplitting` algorithm +For each time step $t_n$, the `DeepSplitting` algorithm -1. Generates the particle trajectories $X^{x, (j)}$ satisfying [Eq. (2)](@ref feynmankac) over the whole interval $[0,T]$. + 1. Generates the particle trajectories $X^{x, (j)}$ satisfying [Eq. (2)](@ref feynmankac) over the whole interval $[0,T]$. -2. Seeks ${\bf u}_{n+1}^{\theta}$ by minimizing the loss function + 2. Seeks ${\bf u}_{n+1}^{\theta}$ by minimizing the loss function ```math L(\theta) = ||{\bf u}^\theta_{n+1}(X_{T - t_n}) - \left[ f(t, X_{T - t_{n-1}}, {\bf u}_{n-1}(X_{T - t_{n-1}}))(t_{n} - t_{n-1}) + {\bf u}_{n-1}(X_{T - t_{n-1}}) \right] || ``` + This way, the PDE approximation problem is decomposed into a sequence of separate learning problems. In `HighDimPDE.jl` the right parameter combination $\theta$ is found by iteratively minimizing $L$ using **stochastic gradient descent**. !!! tip + To solve with `DeepSplitting`, one needs to provide to `solve` - - `dt` - - `batch_size` - - `maxiters`: the number of iterations for minimizing the loss function - - `abstol`: the absolute tolerance for the loss function - - `use_cuda`: if you have a Nvidia GPU, recommended. + + - `dt` + - `batch_size` + - `maxiters`: the number of iterations for minimizing the loss function + - `abstol`: the absolute tolerance for the loss function + - `use_cuda`: if you have a Nvidia GPU, recommended. ## Solving point-wise or on a hypercube ### Pointwise + `DeepSplitting` allows obtaining $u(t,x)$ on a single point $x \in \Omega$ with the keyword $x$. ```julia -prob = PIDEProblem(μ, σ, x, tspan, g, f,) +prob = PIDEProblem(μ, σ, x, tspan, g, f) ``` ### Hypercube + Yet more generally, one wants to solve Eq. (1) on a $d$-dimensional cube $[a,b]^d$. This is offered by `HighDimPDE.jl` with the keyword `x0_sample`. ```julia prob = PIDEProblem(μ, σ, x, tspan, g, f; x0_sample = x0_sample) ``` + Internally, this is handled by assigning a random variable as the initial point of the particles, i.e. + ```math X_t^\xi = \int_0^t \mu(X_s^x)ds + \int_0^t\sigma(X_s^x)dB_s + \xi, ``` + where $\xi$ a random variable uniformly distributed over $[a,b]^d$. This way, the neural network is trained on the whole interval $[a,b]^d$ instead of a single point. ## Non-local PDEs + `DeepSplitting` can solve for non-local reaction diffusion equations of the type + ```math \partial_t u = \mu(x) \nabla_x u + \frac{1}{2} \sigma^2(x) \Delta u + \int_{\Omega}f(x,y, u(t,x), u(t,y))dy ``` @@ -93,13 +111,12 @@ The non-localness is handled by a Monte Carlo integration. u(t_{n+1}, X_{T - t_{n+1}}) \approx \sum_{j=1}^{\text{batch\_size}} \left[ u(t_{n}, X_{T - t_{n}}^{(j)}) + \frac{(t_{n+1} - t_n)}{K}\sum_{k=1}^{K} \big[ f(t, X_{T - t_{n}}^{(j)}, Y_{X_{T - t_{n}}^{(j)}}^{(k)}, u(t_{n},X_{T - t_{n}}^{(j)}), u(t_{n},Y_{X_{T - t_{n}}^{(j)}}^{(k)})) \big] \right] ``` -!!! tip - In practice, if you have a non-local model, you need to provide the sampling method and the number $K$ of MC integration through the keywords `mc_sample` and `K`. - ```julia - alg = DeepSplitting(nn, opt = opt, mc_sample = mc_sample, K = 1) - ``` - `mc_sample` can be whether `UniformSampling(a, b)` or ` NormalSampling(σ_sampling, shifted)`. +!!! tip +In practice, if you have a non-local model, you need to provide the sampling method and the number $K$ of MC integration through the keywords `mc_sample` and `K`. +`julia alg = DeepSplitting(nn, opt = opt, mc_sample = mc_sample, K = 1) ` +`mc_sample` can be whether `UniformSampling(a, b)` or ` NormalSampling(σ_sampling, shifted)`. ## References -- Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. [arXiv](https://arxiv.org/abs/2205.03672) (2022) -- Beck, C., Becker, S., Cheridito, P., Jentzen, A., Neufeld, A., Deep splitting method for parabolic PDEs. [arXiv](https://arxiv.org/abs/1907.03452) (2019) + + - Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. [arXiv](https://arxiv.org/abs/2205.03672) (2022) + - Beck, C., Becker, S., Cheridito, P., Jentzen, A., Neufeld, A., Deep splitting method for parabolic PDEs. [arXiv](https://arxiv.org/abs/1907.03452) (2019) diff --git a/docs/src/Feynman_Kac.md b/docs/src/Feynman_Kac.md index 32532f0..cc6697d 100644 --- a/docs/src/Feynman_Kac.md +++ b/docs/src/Feynman_Kac.md @@ -1,33 +1,39 @@ # [Feynman Kac formula](@id feynmankac) -The Feynman Kac formula is generally stated for terminal condition problems (see e.g. [Wikipedia](https://en.wikipedia.org/wiki/Feynman–Kac_formula)), where +The Feynman Kac formula is generally stated for terminal condition problems (see e.g. [Wikipedia](https://en.wikipedia.org/wiki/Feynman%E2%80%93Kac_formula)), where + ```math \partial_t u(t,x) + \mu(x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(x) \Delta_x u(t,x) + f(x, u(t,x)) = 0 \tag{1} ``` -with terminal condition $u(T, x) = g(x)$, and $u \colon \R^d \to \R$. + +with terminal condition $u(T, x) = g(x)$, and $u \colon \R^d \to \R$. In this case, the FK formula states that for all $t \in (0,T)$ it holds that ```math u(t, x) = \int_t^T \mathbb{E} \left[ f(X^x_{s-t}, u(s, X^x_{s-t}))ds \right] + \mathbb{E} \left[ u(0, X^x_{T-t}) \right] \tag{2} ``` -where + +where + ```math X_t^x = \int_0^t \mu(X_s^x)ds + \int_0^t\sigma(X_s^x)dB_s + x, ``` + and $B_t$ is a [Brownian motion](https://en.wikipedia.org/wiki/Wiener_process). ![Brownian motion - Wikipedia](https://upload.wikimedia.org/wikipedia/commons/f/f8/Wiener_process_3d.png) -Intuitively, this formula is motivated by the fact that [the density of Brownian particles (motion) satisfies the diffusion equation](https://en.wikipedia.org/wiki/Brownian_motion#Einstein's_theory). - +Intuitively, this formula is motivated by the fact that [the density of Brownian particles (motion) satisfies the diffusion equation](https://en.wikipedia.org/wiki/Brownian_motion#Einstein%27s_theory). -The equivalence between the average trajectory of particles and PDEs given by the Feynman-Kac formula allows overcoming the curse of dimensionality that standard numerical methods suffer from, because the expectations can be approximated [Monte Carlo integrations](https://en.wikipedia.org/wiki/Monte_Carlo_integration), which approximation error decreases as $1/\sqrt{N}$ and is therefore not dependent on the dimensions. On the other hand, the computational complexity of traditional deterministic techniques grows exponentially in the number of dimensions. +The equivalence between the average trajectory of particles and PDEs given by the Feynman-Kac formula allows overcoming the curse of dimensionality that standard numerical methods suffer from, because the expectations can be approximated [Monte Carlo integrations](https://en.wikipedia.org/wiki/Monte_Carlo_integration), which approximation error decreases as $1/\sqrt{N}$ and is therefore not dependent on the dimensions. On the other hand, the computational complexity of traditional deterministic techniques grows exponentially in the number of dimensions. ## Forward non-linear Feynman-Kac + > How to transform previous equation to an initial value problem? Define $v(\tau, x) = u(T-\tau, x)$. Observe that $v(0,x) = u(T,x)$. Further, observe that by the chain rule + ```math \begin{aligned} \partial_\tau v(\tau, x) &= \partial_\tau u(T-\tau,x)\\ @@ -36,20 +42,26 @@ Define $v(\tau, x) = u(T-\tau, x)$. Observe that $v(0,x) = u(T,x)$. Further, obs \end{aligned} ``` -From Eq. (1) we get that +From Eq. (1) we get that + ```math - \partial_t u(T - \tau,x) = \mu(x) \nabla_x u(T - \tau,x) + \frac{1}{2} \sigma^2(x) \Delta_x u(T - \tau,x) + f(x, u(T - \tau,x)). ``` + Replacing $u(T-\tau, x)$ by $v(\tau, x)$ we get that $v$ satisfies + ```math \partial_\tau v(\tau, x) = \mu(x) \nabla_x v(\tau,x) + \frac{1}{2} \sigma^2(x) \Delta_x v(\tau,x) + f(x, v(\tau,x)) ``` + and from Eq. (2) we obtain ```math v(\tau, x) = \int_{T-\tau}^T \mathbb{E} \left[ f(X^x_{s- T + \tau}, v(s, X^x_{s-T + \tau}))ds \right] + \mathbb{E} \left[ v(0, X^x_{\tau}) \right]. ``` -By using the substitution rule with $\tau \to \tau -T$ (shifting by T) and $\tau \to - \tau$ (inversing), and finally inversing the integral bound we get that + +By using the substitution rule with $\tau \to \tau -T$ (shifting by T) and $\tau \to - \tau$ (inversing), and finally inversing the integral bound we get that + ```math \begin{aligned} v(\tau, x) &= \int_{-\tau}^0 \mathbb{E} \left[ f(X^x_{s + \tau}, v(s + T, X^x_{s + \tau}))ds \right] + \mathbb{E} \left[ v(0, X^x_{\tau}) \right]\\ @@ -58,18 +70,25 @@ v(\tau, x) &= \int_{-\tau}^0 \mathbb{E} \left[ f(X^x_{s + \tau}, v(s + T, X^x_{s \end{aligned} ``` -This leads to the +This leads to the + !!! info "Non-linear Feynman Kac for initial value problems" + Consider the PDE + ```math \partial_t u(t,x) = \mu(t, x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) + f(x, u(t,x)) ``` - with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$. + + with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$. Then + ```math u(t, x) = \int_0^t \mathbb{E} \left[ f(X^x_{t - s}, u(T-s, X^x_{t - s}))ds \right] + \mathbb{E} \left[ u(0, X^x_t) \right] \tag{3} ``` - with + + with + ```math X_t^x = \int_0^t \mu(X_s^x)ds + \int_0^t\sigma(X_s^x)dB_s + x. - ``` \ No newline at end of file + ``` diff --git a/docs/src/MLP.md b/docs/src/MLP.md index 3e603da..44bce15 100644 --- a/docs/src/MLP.md +++ b/docs/src/MLP.md @@ -1,39 +1,46 @@ # [The `MLP` algorithm](@id mlp) ### Problems Supported: -1. [`PIDEProblem`](@ref) -2. [`ParabolicPDEProblem`](@ref) + + 1. [`PIDEProblem`](@ref) + 2. [`ParabolicPDEProblem`](@ref) ```@autodocs Modules = [HighDimPDE] Pages = ["MLP.jl"] ``` -The `MLP`, for Multi-Level Picard iterations, reformulates the PDE problem as a fixed point equation through the Feynman Kac formula. +The `MLP`, for Multi-Level Picard iterations, reformulates the PDE problem as a fixed point equation through the Feynman Kac formula. -- It relies on [Picard iterations](https://en.wikipedia.org/wiki/Picard–Lindelöf_theorem) to find the fixed point, + - It relies on [Picard iterations](https://en.wikipedia.org/wiki/Picard%E2%80%93Lindel%C3%B6f_theorem) to find the fixed point, -- reducing the complexity of the numerical approximation of the time integral through a [multilevel Monte Carlo](https://en.wikipedia.org/wiki/Multilevel_Monte_Carlo_method) approach. + - reducing the complexity of the numerical approximation of the time integral through a [multilevel Monte Carlo](https://en.wikipedia.org/wiki/Multilevel_Monte_Carlo_method) approach. The `MLP` algorithm overcomes the curse of dimensionality, with a computational complexity that grows polynomially in the number of dimension (see [M. Hutzenthaler et al. 2020](https://arxiv.org/abs/1807.01212v3)). !!! warning "`MLP` can only approximate the solution on a single point" + `MLP` only works for `PIDEProblem` with `x0_sample = NoSampling()`. If you want to solve over an entire domain, you definitely want to check the `DeepSplitting` algorithm. ## The general idea 💡 + Consider the PDE + ```math \partial_t u(t,x) = \mu(t, x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) + f(x, u(t,x)) \tag{1} ``` -with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$. + +with initial conditions $u(0, x) = g(x)$, where $u \colon \R^d \to \R$. ### Picard Iterations -The `MLP` algorithm observes that the [Feynman Kac formula](@ref feynmankac) can be viewed as a fixed point equation, i.e. $u = \phi(u)$. Introducing a sequence $(u_k)$ defined as $u_0 = g$ and + +The `MLP` algorithm observes that the [Feynman Kac formula](@ref feynmankac) can be viewed as a fixed point equation, i.e. $u = \phi(u)$. Introducing a sequence $(u_k)$ defined as $u_0 = g$ and + ```math u_{l+1} = \phi(u_l), ``` -the [Banach fixed-point theorem](https://en.wikipedia.org/wiki/Banach_fixed-point_theorem) ensures that the sequence converges to the true solution $u$. Such a technique is known as [Picard iterations](https://en.wikipedia.org/wiki/Picard–Lindelöf_theorem). +the [Banach fixed-point theorem](https://en.wikipedia.org/wiki/Banach_fixed-point_theorem) ensures that the sequence converges to the true solution $u$. Such a technique is known as [Picard iterations](https://en.wikipedia.org/wiki/Picard%E2%80%93Lindel%C3%B6f_theorem). The time integral term is evaluated by a [Monte-Carlo integration](https://en.wikipedia.org/wiki/Monte_Carlo_integration) @@ -41,11 +48,11 @@ The time integral term is evaluated by a [Monte-Carlo integration](https://en.wi u_L = \frac{1}{M}\sum_i^M \left[ f(X^{x,(i)}_{t - s_{(l, i)}}, u_{L-1}(T-s_i, X^{x,( i)}_{t - s_{(l, i)}})) + u(0, X^{x,(i)}_{t - s_{(l, i)}}) \right]. ``` -But the MLP uses an extra trick to lower the computational cost of the iteration. - +But the MLP uses an extra trick to lower the computational cost of the iteration. ### Telescope sum -The `MLP` algorithm uses a telescope sum + +The `MLP` algorithm uses a telescope sum ```math \begin{aligned} @@ -56,12 +63,13 @@ u_L = \phi(u_{L-1}) &= [\phi(u_{L-1}) - \phi(u_{L-2})] + [\phi(u_{L-2}) - \phi(u As $l$ grows, the term $[\phi(u_{l-1}) - \phi(u_{l-2})]$ becomes smaller, and thus demands more calculations. The `MLP` algorithm uses this fact by evaluating the integral term at level $l$ with $M^{L-l}$ samples. - !!! tip - - `L` corresponds to the level of the approximation, i.e. $u \approx u_L$ - - `M` characterizes the number of samples for the Monte Carlo approximation of the time integral + + - `L` corresponds to the level of the approximation, i.e. $u \approx u_L$ + - `M` characterizes the number of samples for the Monte Carlo approximation of the time integral Overall, `MLP` can be summarized by the following formula + ```math \begin{aligned} u_L &= \sum_{l=1}^{L-1} \frac{1}{M^{L-l}}\sum_i^{M^{L-l}} \left[ f(X^{x,(l, i)}_{t - s_{(l, i)}}, u(T-s_{(l, i)}, X^{x,(l, i)}_{t - s_{(l, i)}})) + \mathbf{1}_\N(l) f(X^{x,(l, i)}_{t - s_{(l, i)}}, u(T-s_{(l, i)}, X^{x,(l, i)}_{t - s_{(l, i)}}))\right] @@ -69,10 +77,13 @@ u_L &= \sum_{l=1}^{L-1} \frac{1}{M^{L-l}}\sum_i^{M^{L-l}} \left[ f(X^{x,(l, i)}_ &\qquad + \frac{1}{M^{L}}\sum_i^{M^{L}} u(0, X^{x,(l, i)}_t)\\ \end{aligned} ``` + Note that the superscripts $(l, i)$ indicate the independence of the random variables $X$ across levels. ## Non-local PDEs + `MLP` can solve for non-local reaction diffusion equations of the type + ```math \partial_t u = \mu(t, x) \nabla_x u(t, x) + \frac{1}{2} \sigma^2(t, x) \Delta u(t, x) + \int_{\Omega}f(x, y, u(t,x), u(t,y))dy ``` @@ -88,10 +99,13 @@ u_L &= \sum_{l=1}^{L-1} \frac{1}{M^{L-l}}\sum_{i=1}^{M^{L-l}} \frac{1}{K}\sum_{j ``` !!! tip - In practice, if you have a non-local model, you need to provide the sampling method and the number $K$ of MC integration through the keywords `mc_sample` and `K`. - - `K` characterizes the number of samples for the Monte Carlo approximation of the last term. - - `mc_sample` characterizes the distribution of the `Z` variables + + In practice, if you have a non-local model, you need to provide the sampling method and the number $K$ of MC integration through the keywords `mc_sample` and `K`. + + - `K` characterizes the number of samples for the Monte Carlo approximation of the last term. + - `mc_sample` characterizes the distribution of the `Z` variables ## References -- Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. [arXiv](https://arxiv.org/abs/2205.03672) (2022) -- Becker, S., Braunwarth, R., Hutzenthaler, M., Jentzen, A., von Wurstemberger, P., Numerical simulations for full history recursive multilevel Picard approximations for systems of high-dimensional partial differential equations. [arXiv](https://arxiv.org/abs/2005.10206) (2020) \ No newline at end of file + + - Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. [arXiv](https://arxiv.org/abs/2205.03672) (2022) + - Becker, S., Braunwarth, R., Hutzenthaler, M., Jentzen, A., von Wurstemberger, P., Numerical simulations for full history recursive multilevel Picard approximations for systems of high-dimensional partial differential equations. [arXiv](https://arxiv.org/abs/2005.10206) (2020) diff --git a/docs/src/NNKolmogorov.md b/docs/src/NNKolmogorov.md index 669f01a..07b13cb 100644 --- a/docs/src/NNKolmogorov.md +++ b/docs/src/NNKolmogorov.md @@ -1,30 +1,40 @@ # [The `NNKolmogorov` algorithm](@id nn_kolmogorov) ### Problems Supported: -1. [`ParabolicPDEProblem`](@ref) + + 1. [`ParabolicPDEProblem`](@ref) ```@autodocs Modules = [HighDimPDE] Pages = ["NNKolmogorov.jl"] ``` -`NNKolmogorov` obtains a -- terminal solution for Forward Kolmogorov Equations of the form: +`NNKolmogorov` obtains a + + - terminal solution for Forward Kolmogorov Equations of the form: + ```math \partial_t u(t,x) = \mu(t, x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) ``` + with initial condition given by `g(x)` -- or an initial condition for Backward Kolmogorov Equations of the form: + + - or an initial condition for Backward Kolmogorov Equations of the form: + ```math \partial_t u(t,x) = - \mu(t, x) \nabla_x u(t,x) - \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) ``` + with terminal condition given by `g(x)` -We can use the Feynman-Kac formula : +We can use the Feynman-Kac formula : + ```math S_t^x = \int_{0}^{t}\mu(S_s^x)ds + \int_{0}^{t}\sigma(S_s^x)dB_s ``` + And the solution is given by: + ```math f(T, x) = \mathbb{E}[g(S_T^x)] -``` \ No newline at end of file +``` diff --git a/docs/src/NNParamKolmogorov.md b/docs/src/NNParamKolmogorov.md index 0b36d93..5eb5f6e 100644 --- a/docs/src/NNParamKolmogorov.md +++ b/docs/src/NNParamKolmogorov.md @@ -1,30 +1,40 @@ # [The `NNParamKolmogorov` algorithm](@id nn_paramkolmogorov) ### Problems Supported: -1. [`ParabolicPDEProblem`](@ref) + + 1. [`ParabolicPDEProblem`](@ref) ```@autodocs Modules = [HighDimPDE] Pages = ["NNParamKolmogorov.jl"] ``` -`NNParamKolmogorov` obtains a -- terminal solution for parametric families of Forward Kolmogorov Equations of the form: +`NNParamKolmogorov` obtains a + + - terminal solution for parametric families of Forward Kolmogorov Equations of the form: + ```math \partial_t u(t,x) = \mu(t, x, γ_mu) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(t, x, γ_sigma) \Delta_x u(t,x) ``` + with initial condition given by `g(x, γ_phi)` -- or an initial condition for parametric families of Backward Kolmogorov Equations of the form: + + - or an initial condition for parametric families of Backward Kolmogorov Equations of the form: + ```math \partial_t u(t,x) = - \mu(t, x) \nabla_x u(t,x) - \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) ``` + with terminal condition given by `g(x, γ_phi)` -We can use the Feynman-Kac formula : +We can use the Feynman-Kac formula : + ```math S_t^x = \int_{0}^{t}\mu(S_s^x)ds + \int_{0}^{t}\sigma(S_s^x)dB_s ``` + And the solution is given by: + ```math f(T, x) = \mathbb{E}[g(S_T^x, γ_phi)] -``` \ No newline at end of file +``` diff --git a/docs/src/NNStopping.md b/docs/src/NNStopping.md index df773dc..233e143 100644 --- a/docs/src/NNStopping.md +++ b/docs/src/NNStopping.md @@ -1,15 +1,18 @@ # [The `NNStopping` algorithm](@id nn_stopping) ### Problems Supported: -1. [`ParabolicPDEProblem`](@ref) + + 1. [`ParabolicPDEProblem`](@ref) ```@autodocs Modules = [HighDimPDE] Pages = ["NNStopping.jl"] ``` + ## The general idea 💡 Similar to DeepSplitting and DeepBSDE, NNStopping evaluates the PDE as a Stochastic Differential Equation. Consider an Obstacle PDE of the form: + ```math max\lbrace\partial_t u(t,x) + \mu(t, x) \nabla_x u(t,x) + \frac{1}{2} \sigma^2(t, x) \Delta_x u(t,x) , g(t,x) - u(t,x)\rbrace ``` @@ -27,6 +30,7 @@ The payoff of the option would then be: ```math sup\lbrace\mathbb{E}[g(X_\tau, \tau)]\rbrace ``` + Where τ is the stopping (exercising) time. The goal is to retrieve both the optimal exercising strategy (τ) and the payoff. We approximate each stopping decision with a neural network architecture, inorder to maximise the expected payoff. diff --git a/docs/src/examples/blackscholes.md b/docs/src/examples/blackscholes.md index 59ff3f1..376e912 100644 --- a/docs/src/examples/blackscholes.md +++ b/docs/src/examples/blackscholes.md @@ -12,14 +12,14 @@ 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) +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 +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) ``` @@ -29,16 +29,16 @@ by giving it the Flux.jl chains we want it to use for the neural networks. needs to be `d+1`-dimensional to `d` dimensions. Thus we define the following: ```julia -hls = 10 + d #hide layer size +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) +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 @@ -47,10 +47,10 @@ SDE using the Euler-Maruyama SDE solver with our chosen `dt=0.2`, do at most 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) +ans = solve(prob, pdealg, verbose = true, maxiters = 150, trajectories = 100, + alg = EM(), dt = 0.2, pabstol = 1.0f-6) ``` ## References -1. Shinde, A. S., and K. C. Takale. "Study of Black-Scholes model and its applications." Procedia Engineering 38 (2012): 270-279. \ No newline at end of file + 1. Shinde, A. S., and K. C. Takale. "Study of Black-Scholes model and its applications." Procedia Engineering 38 (2012): 270-279. diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index f7da87f..a52f073 100644 --- a/docs/src/getting_started.md +++ b/docs/src/getting_started.md @@ -4,9 +4,9 @@ The general workflow for using `HighDimPDE.jl` is as follows: -- Define a Partial Integro-Differential Equation problem -- Select a solver algorithm -- Solve the problem. + - Define a Partial Integro-Differential Equation problem + - Select a solver algorithm + - Solve the problem. ## Examples @@ -28,11 +28,11 @@ using HighDimPDE ## Definition of the problem d = 10 # dimension of the problem tspan = (0.0, 0.5) # time horizon -x0 = fill(0., d) # initial point -g(x) = exp(-sum(x.^2)) # initial condition +x0 = fill(0.0, d) # initial point +g(x) = exp(-sum(x .^ 2)) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients -f(x, v_x, ∇v_x, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the +f(x, v_x, ∇v_x, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f) # defining the problem ## Definition of the algorithm @@ -45,9 +45,11 @@ sol = solve(prob, alg, multithreading = true) #### Non-local PDE with Neumann boundary conditions Let's include in the previous equation non-local competition, i.e. + ```math \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$. ```@example MLP_non_local_PDE @@ -55,13 +57,13 @@ using HighDimPDE ## Definition of the problem d = 10 # dimension of the problem -tspan = (0.0,0.5) # time horizon -x0 = fill(0.,d) # initial point -g(x) = exp( -sum(x.^2) ) # initial condition +tspan = (0.0, 0.5) # time horizon +x0 = fill(0.0, d) # initial point +g(x) = exp(-sum(x .^ 2)) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients -mc_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) -f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) +mc_sample = UniformSampling(fill(-5.0f-1, d), fill(5.0f-1, d)) +f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) prob = PIDEProblem(μ, σ, x0, tspan, g, f) # defining x0_sample is sufficient to implement Neumann boundary conditions ## Definition of the algorithm @@ -81,36 +83,36 @@ using Flux # needed to define the neural network ## Definition of the problem d = 10 # dimension of the problem tspan = (0.0, 0.5) # time horizon -x0 = fill(0f0, d) # initial point -g(x) = exp.(-sum(x.^2, dims=1)) # initial condition +x0 = fill(0.0f0, d) # initial point +g(x) = exp.(-sum(x .^ 2, dims = 1)) # initial condition μ(x, p, t) = 0.0f0 # advection coefficients σ(x, p, t) = 0.1f0 # diffusion coefficients -x0_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) -f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = v_x .* (1f0 .- v_y) +x0_sample = UniformSampling(fill(-5.0f-1, d), fill(5.0f-1, d)) +f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = v_x .* (1.0f0 .- v_y) prob = PIDEProblem(μ, σ, x0, tspan, g, f; - x0_sample = x0_sample) + x0_sample = x0_sample) ## Definition of the neural network to use hls = d + 50 #hidden layer size nn = Flux.Chain(Dense(d, hls, tanh), - Dense(hls, hls, tanh), - Dense(hls, 1)) # neural network used by the scheme + Dense(hls, hls, tanh), + Dense(hls, 1)) # neural network used by the scheme opt = ADAM(1e-2) ## Definition of the algorithm alg = DeepSplitting(nn, - opt = opt, - mc_sample = x0_sample) - -sol = solve(prob, - alg, - 0.1, - verbose = true, - abstol = 2e-3, - maxiters = 1000, - batch_size = 1000) + opt = opt, + mc_sample = x0_sample) + +sol = solve(prob, + alg, + 0.1, + verbose = true, + abstol = 2e-3, + maxiters = 1000, + batch_size = 1000) ``` ### Solving on the GPU @@ -118,12 +120,12 @@ sol = solve(prob, [`DeepSplitting`](@ref deepsplitting) can run on the GPU for (much) improved performance. To do so, just set `use_cuda = true`. ```@example DeepSplitting_non_local_PDE -sol = solve(prob, - alg, - 0.1, - verbose = true, - abstol = 2e-3, - maxiters = 1000, - batch_size = 1000, - use_cuda=true) -``` \ No newline at end of file +sol = solve(prob, + alg, + 0.1, + verbose = true, + abstol = 2e-3, + maxiters = 1000, + batch_size = 1000, + use_cuda = true) +``` diff --git a/docs/src/index.md b/docs/src/index.md index 7c50a08..c303375 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,9 +1,9 @@ # HighDimPDE.jl +**HighDimPDE.jl** is a Julia package to **solve Highly Dimensional non-linear, non-local PDEs** of the forms: -**HighDimPDE.jl** is a Julia package to **solve Highly Dimensional non-linear, non-local PDEs** of the forms: +### 1. Partial Integro Differential Equations: -### 1. Partial Integro Differential Equations: ```math \begin{aligned} (\partial_t u)(t,x) &= \int_{\Omega} f\big(t,x,{\bf x}, u(t,x),u(t,{\bf x}), ( \nabla_x u )(t,x ),( \nabla_x u )(t,{\bf x} ) \big) \, d{\bf x} \\ @@ -13,7 +13,8 @@ where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject to initial and boundary conditions, and where $d$ is large. These equations are defined using the [`PIDEProblem`](@ref) -### 2. Parabolic Partial Differential Equations: +### 2. Parabolic Partial Differential Equations: + ```math \begin{aligned} (\partial_t u)(t,x) &= f\big(t,x, u(t,x), ( \nabla_x u )(t,x )\big) @@ -24,71 +25,80 @@ where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subseteq \R^d$ is subject to initial and boundary conditions, and where $d$ is large. These equations are defined using the [`ParabolicPDEProblem`](@ref) !!! note + The difference between the two problems is that in Partial Integro Differential Equations, the integrand is integrated over **x**, while in Parabolic Integro Differential Equations, the function `f` is just evaluated for `x`. **HighDimPDE.jl** implements solver algorithms that break down the curse of dimensionality, including -* the [Deep Splitting scheme](@ref deepsplitting) - -* the [Multi-Level Picard iterations scheme](@ref mlp) - -* [the Deep BSDE scheme](@ref deepbsde). - -* [NNKolmogorov](@ref nn_kolmogorov) and [NNParamKolmogorov](@ref nn_paramkolmogorov) schemes for Kolmogorov PDEs. + - the [Deep Splitting scheme](@ref deepsplitting) -* [NNStopping] scheme for solving optimal stopping time problem. + - the [Multi-Level Picard iterations scheme](@ref mlp) + - [the Deep BSDE scheme](@ref deepbsde). + - [NNKolmogorov](@ref nn_kolmogorov) and [NNParamKolmogorov](@ref nn_paramkolmogorov) schemes for Kolmogorov PDEs. + - [NNStopping] scheme for solving optimal stopping time problem. -To make the most out of **HighDimPDE.jl**, we advise to first have a look at the +To make the most out of **HighDimPDE.jl**, we advise to first have a look at the -* [documentation on the Feynman Kac formula](@ref feynmankac), + - [documentation on the Feynman Kac formula](@ref feynmankac), as all solver algorithms heavily rely on it. ## Algorithm overview ------------------------------------------------------------- +* * * + Features | `DeepSplitting` | `MLP` | `DeepBSDE` | ----------|:----------------------:|:------------:|:--------: Time discretization free| ❌ | ✅ | ❌ | Mesh-free | ✅ | ✅ | ✅ | Single point $x \in \R^d$ approximation| ✅ | ✅ | ✅ | $d$-dimensional cube $[a,b]^d$ approximation| ✅ | ❌ | ✔️ | -GPU | ✅ | ❌ | ✅ | +GPU | ✅ | ❌ | ✅ | Gradient non-linearities | ✔️| ❌ | ✅ | ✔️ : will be supported in the future ## Reproducibility + ```@raw html
The documentation of this SciML package was built using these direct dependencies, ``` + ```@example using Pkg # hide Pkg.status() # hide ``` + ```@raw html
``` + ```@raw html
and using this machine and Julia version. ``` + ```@example using InteractiveUtils # hide versioninfo() # hide ``` + ```@raw html
``` + ```@raw html
A more complete overview of all dependencies and their versions is also provided. ``` + ```@example using Pkg # hide -Pkg.status(;mode = PKGMODE_MANIFEST) # hide +Pkg.status(; mode = PKGMODE_MANIFEST) # hide ``` + ```@raw html
``` + ```@eval using TOML using Markdown diff --git a/docs/src/problems.md b/docs/src/problems.md index 6ae5eeb..ebff55d 100644 --- a/docs/src/problems.md +++ b/docs/src/problems.md @@ -3,6 +3,6 @@ PIDEProblem ParabolicPDEProblem ``` -!!! note - While choosing to define a PDE using `PIDEProblem`, note that the function being integrated `f` is a function of `f(x, y, v_x, v_y, ∇v_x, ∇v_y)` out of which `y` is the integrating variable and `x` is constant throughout the integration. - If a PDE has no integral and the non linear term `f` is just evaluated as `f(x, v_x, ∇v_x)` then we suggest using `ParabolicPDEProblem` \ No newline at end of file +!!! note +While choosing to define a PDE using `PIDEProblem`, note that the function being integrated `f` is a function of `f(x, y, v_x, v_y, ∇v_x, ∇v_y)` out of which `y` is the integrating variable and `x` is constant throughout the integration. +If a PDE has no integral and the non linear term `f` is just evaluated as `f(x, v_x, ∇v_x)` then we suggest using `ParabolicPDEProblem` diff --git a/docs/src/tutorials/deepbsde.md b/docs/src/tutorials/deepbsde.md index f439e18..b7f8b2a 100644 --- a/docs/src/tutorials/deepbsde.md +++ b/docs/src/tutorials/deepbsde.md @@ -44,7 +44,6 @@ pdealg = DeepBSDE(u0, σᵀ∇u, opt = opt) trajectories = 30, dt = 1.2f0, pabstol = 1.0f-4) - ``` Now, let's explain the details! @@ -62,6 +61,7 @@ Here, we choose to solve the classical Linear Quadratic Gaussian ```math d X_t = 2 \sqrt{\lambda} c_t dt + \sqrt{2}dW_t ``` + where $c_t$ is a control process. The solution to the optimal control is given by a PDE of the form: ```math @@ -83,14 +83,14 @@ using StochasticDiffEq using LinearAlgebra d = 100 # number of dimensions -X0 = fill(0.0f0,d) # initial value of stochastic control process +X0 = fill(0.0f0, d) # initial value of stochastic control process tspan = (0.0f0, 1.0f0) λ = 1.0f0 -g(X) = log(0.5f0 + 0.5f0*sum(X.^2)) -f(X,u,σᵀ∇u,p,t) = -λ*sum(σᵀ∇u.^2) -μ_f(X,p,t) = zero(X) #Vector d x 1 λ -σ_f(X,p,t) = Diagonal(sqrt(2.0f0)*ones(Float32,d)) #Matrix d x d +g(X) = log(0.5f0 + 0.5f0 * sum(X .^ 2)) +f(X, u, σᵀ∇u, p, t) = -λ * sum(σᵀ∇u .^ 2) +μ_f(X, p, t) = zero(X) #Vector d x 1 λ +σ_f(X, p, t) = Diagonal(sqrt(2.0f0) * ones(Float32, d)) #Matrix d x d prob = ParabolicPDEProblem(μ_f, σ_f, X0, tspan; g, f) ``` @@ -105,22 +105,22 @@ needs to be `d+1` dimensional to `d` dimensions. Thus we define the following: hls = 10 + d #hidden layer size opt = Flux.Optimise.Adam(0.01) #optimizer #sub-neural network approximating solutions at the desired point -u0 = Flux.Chain(Dense(d,hls,relu), - Dense(hls,hls,relu), - Dense(hls,1)) +u0 = Flux.Chain(Dense(d, hls, relu), + Dense(hls, hls, relu), + Dense(hls, 1)) # sub-neural network approximating the spatial gradients at time point -σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu), - Dense(hls,hls,relu), - Dense(hls,hls,relu), - Dense(hls,d)) +σᵀ∇u = Flux.Chain(Dense(d + 1, hls, relu), + Dense(hls, hls, relu), + Dense(hls, hls, relu), + Dense(hls, d)) pdealg = DeepBSDE(u0, σᵀ∇u, opt = opt) ``` #### Solving with Neural Nets ```@example deepbsde2 -@time ans = solve(prob, pdealg, EM(), verbose=true, maxiters=100, trajectories=100, dt=0.2f0, pabstol = 1f-2) - +@time ans = solve(prob, pdealg, EM(), verbose = true, maxiters = 100, + trajectories = 100, dt = 0.2f0, pabstol = 1.0f-2) ``` Here we want to solve the underlying neural @@ -142,14 +142,14 @@ To solve it using the `ParabolicPDEProblem`, 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) +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 +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 = ParabolicPDEProblem(μ_f, σ_f, X0, tspan; g, f) ``` @@ -159,16 +159,16 @@ by giving it the Flux.jl chains we want it to use for the neural networks. needs to be `d+1`-dimensional to `d` dimensions. Thus we define the following: ```julia -hls = 10 + d #hide layer size +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 = DeepBSDE(u0, σᵀ∇u, opt=opt) +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 = DeepBSDE(u0, σᵀ∇u, opt = opt) ``` And now we solve the PDE. Here, we say we want to solve the underlying neural @@ -177,9 +177,10 @@ SDE using the Euler-Maruyama SDE solver with our chosen `dt=0.2`, do at most and stop if the loss ever goes below `1f-6`. ```julia -ans = solve(prob, pdealg, EM(), verbose=true, maxiters=150, trajectories=100, dt=0.2f0) +ans = solve( + prob, pdealg, EM(), verbose = true, maxiters = 150, trajectories = 100, dt = 0.2f0) ``` ## References -1. Shinde, A. S., and K. C. Takale. "Study of Black-Scholes model and its applications." Procedia Engineering 38 (2012): 270-279. + 1. Shinde, A. S., and K. C. Takale. "Study of Black-Scholes model and its applications." Procedia Engineering 38 (2012): 270-279. diff --git a/docs/src/tutorials/deepsplitting.md b/docs/src/tutorials/deepsplitting.md index 138c7c9..73701da 100644 --- a/docs/src/tutorials/deepsplitting.md +++ b/docs/src/tutorials/deepsplitting.md @@ -1,5 +1,5 @@ - # Solving the 10-dimensional Fisher-KPP equation with `DeepSplitting` + Consider the Fisher-KPP equation with non-local competition ```math @@ -8,7 +8,7 @@ Consider the Fisher-KPP equation with non-local competition where $\Omega = [-1/2, 1/2]^d$, and let's assume Neumann Boundary condition on $\Omega$. -Let's solve Eq. (1) with the [`DeepSplitting`](@ref deepsplitting) solver. +Let's solve Eq. (1) with the [`DeepSplitting`](@ref deepsplitting) solver. ```@example deepsplitting using HighDimPDE, Flux @@ -16,40 +16,43 @@ using HighDimPDE, Flux ## Definition of the problem d = 10 # dimension of the problem tspan = (0.0, 0.5) # time horizon -x0 = fill(0f0,d) # initial point -g(x) = exp.(- sum(x.^2, dims=1) ) # initial condition +x0 = fill(0.0f0, d) # initial point +g(x) = exp.(-sum(x .^ 2, dims = 1)) # initial condition μ(x, p, t) = 0.0f0 # advection coefficients σ(x, p, t) = 0.1f0 # diffusion coefficients -x0_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) -f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = v_x .* (1f0 .- v_y) +x0_sample = UniformSampling(fill(-5.0f-1, d), fill(5.0f-1, d)) +f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = v_x .* (1.0f0 .- v_y) ``` + Since this is a non-local equation, we will define our problem as a [`PIDEProblem`](@ref) + ```@example deepsplitting prob = PIDEProblem(μ, σ, x0, tspan, g, f; x0_sample = x0_sample) ``` + ```@example deepsplitting ## Definition of the neural network to use hls = d + 50 #hidden layer size nn = Flux.Chain(Dense(d, hls, tanh), - Dense(hls, hls, tanh), - Dense(hls, 1)) # neural network used by the scheme + Dense(hls, hls, tanh), + Dense(hls, 1)) # neural network used by the scheme opt = ADAM(1e-2) ## Definition of the algorithm alg = DeepSplitting(nn, - opt = opt, - mc_sample = x0_sample) - -sol = solve(prob, - alg, - 0.1, - verbose = true, - abstol = 2e-3, - maxiters = 1000, - batch_size = 1000) + opt = opt, + mc_sample = x0_sample) + +sol = solve(prob, + alg, + 0.1, + verbose = true, + abstol = 2e-3, + maxiters = 1000, + batch_size = 1000) ``` #### Solving on the GPU @@ -57,18 +60,19 @@ sol = solve(prob, `DeepSplitting` can run on the GPU for (much) improved performance. To do so, just set `use_cuda = true`. ```@example deepsplitting -sol = solve(prob, - alg, - 0.1, - verbose = true, - abstol = 2e-3, - maxiters = 1000, - batch_size = 1000, - use_cuda=true) +sol = solve(prob, + alg, + 0.1, + verbose = true, + abstol = 2e-3, + maxiters = 1000, + batch_size = 1000, + use_cuda = true) ``` # Solving local Allen Cahn Equation with Neumann BC -```@example deepsplitting2 + +```@example deepsplitting2 using HighDimPDE, Flux, StochasticDiffEq batch_size = 1000 train_steps = 1000 @@ -85,8 +89,8 @@ d = 10 hls = d + 50 #hidden layer size nn = Flux.Chain(Dense(d, hls, relu), - Dense(hls, hls, relu), - Dense(hls, 1)) # Neural network used by the scheme + Dense(hls, hls, relu), + Dense(hls, 1)) # Neural network used by the scheme opt = Flux.Optimise.Adam(1e-2) #optimiser alg = DeepSplitting(nn, opt = opt) @@ -96,19 +100,22 @@ g(X) = exp.(-0.25f0 * sum(X .^ 2, dims = 1)) # initial condition a(u) = u - u^3 f(y, v_y, ∇v_y, p, t) = a.(v_y) # nonlocal nonlinear function ``` + Since we are dealing with a local problem here, i.e. no integration term used, we use [`ParabolicPDEProblem`](@ref) to define the problem. + ```@example deepsplitting2 # defining the problem prob = ParabolicPDEProblem(μ, σ, X0, tspan; g, f, neumann_bc = [-∂, ∂]) ``` + ```@example deepsplitting2 # solving @time sol = solve(prob, - alg, - dt, - verbose = false, - abstol = 1e-5, - use_cuda = false, - maxiters = train_steps, - batch_size = batch_size) -``` \ No newline at end of file + alg, + dt, + verbose = false, + abstol = 1e-5, + use_cuda = false, + maxiters = train_steps, + batch_size = batch_size) +``` diff --git a/docs/src/tutorials/mlp.md b/docs/src/tutorials/mlp.md index 15ecf69..5f8431c 100644 --- a/docs/src/tutorials/mlp.md +++ b/docs/src/tutorials/mlp.md @@ -1,7 +1,7 @@ - # 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} ``` @@ -11,43 +11,47 @@ using HighDimPDE ## Definition of the problem d = 10 # dimension of the problem -tspan = (0.0,0.5) # time horizon -x0 = fill(0.,d) # initial point -g(x) = exp(- sum(x.^2) ) # initial condition +tspan = (0.0, 0.5) # time horizon +x0 = fill(0.0, d) # initial point +g(x) = exp(-sum(x .^ 2)) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients -f(x, v_x, ∇v_x, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the +f(x, v_x, ∇v_x, p, t) = max(0.0, v_x) * (1 - max(0.0, v_x)) # nonlocal nonlinear part of the prob = ParabolicPDEProblem(μ, σ, x0, tspan; g, f) # defining the problem ## Definition of the algorithm alg = MLP() # defining the algorithm. We use the Multi Level Picard algorithm ## Solving with multiple threads -sol = solve(prob, alg, multithreading=true) +sol = solve(prob, alg, multithreading = true) ``` ## Non local PDE with Neumann boundary conditions + Let's include in the previous equation non local competition, i.e. + ```math \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$. + ```@example mlp2 using HighDimPDE ## Definition of the problem d = 10 # dimension of the problem -tspan = (0.0,0.5) # time horizon -x0 = fill(0.,d) # initial point -g(x) = exp( -sum(x.^2) ) # initial condition +tspan = (0.0, 0.5) # time horizon +x0 = fill(0.0, d) # initial point +g(x) = exp(-sum(x .^ 2)) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients -mc_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d)) -f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) +mc_sample = UniformSampling(fill(-5.0f-1, d), fill(5.0f-1, d)) +f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) prob = PIDEProblem(μ, σ, x0, tspan, g, f) # defining x0_sample is sufficient to implement Neumann boundary conditions ## Definition of the algorithm -alg = MLP(mc_sample = mc_sample ) +alg = MLP(mc_sample = mc_sample) -sol = solve(prob, alg, multithreading=true) +sol = solve(prob, alg, multithreading = true) ``` diff --git a/docs/src/tutorials/nnkolmogorov.md b/docs/src/tutorials/nnkolmogorov.md index 29fd355..2c1fad7 100644 --- a/docs/src/tutorials/nnkolmogorov.md +++ b/docs/src/tutorials/nnkolmogorov.md @@ -6,12 +6,12 @@ using HighDimPDE, StochasticDiffEq, Flux, LinearAlgebra d = 10 # dims -T = 1/12 -sigma = 0.01 .+ 0.03.*Matrix(Diagonal(ones(d))) # volatility +T = 1 / 12 +sigma = 0.01 .+ 0.03 .* Matrix(Diagonal(ones(d))) # volatility mu = 0.06 # interest rate K = 100.0 # strike price function μ_func(du, u, p, t) - du .= mu*u + du .= mu * u end function σ_func(du, u, p, t) @@ -22,7 +22,7 @@ tspan = (0.0, T) # The range for initial stock price xspan = [(98.00, 102.00) for i in 1:d] -g(x) = max(maximum(x) -K, 0) +g(x) = max(maximum(x) - K, 0) sdealg = SRIW1() # provide `x0` as nothing to the problem since we are provinding a range for `x0`. diff --git a/docs/src/tutorials/nnparamkolmogorov.md b/docs/src/tutorials/nnparamkolmogorov.md index f9de711..e1507b2 100644 --- a/docs/src/tutorials/nnparamkolmogorov.md +++ b/docs/src/tutorials/nnparamkolmogorov.md @@ -3,6 +3,7 @@ ## Solving Parametric Family of High Dimensional Heat Equation. In this example we will solve the high dimensional heat equation over a range of initial values, and also over a range of thermal diffusivity. + ```@example nnparamkolmogorov using HighDimPDE, Flux, StochasticDiffEq d = 10 @@ -20,7 +21,7 @@ function phi(x, y_phi) sum(x .^ 2) end function sigma_(dx, x, γ_sigma, t) - dx .= γ_sigma[:, :, 1] + dx .= γ_sigma[:, :, 1] end mu_(dx, x, γ_mu, t) = dx .= 0.00 @@ -47,16 +48,19 @@ sol = solve(prob, NNParamKolmogorov(m, opt), sdealg, verbose = true, dt = 0.01, abstol = 1e-10, dx = 0.1, trajectories = trajectories, maxiters = 1000, use_gpu = false, dps = dps) ``` + Similarly we can parametrize the drift function `mu_` and the initial function `g`, and obtain a solution over all parameters and initial values. # Inferring on the solution from `NNParamKolmogorov`: + ```@example nnparamkolmogorov x_test = rand(xspan[1][1]:0.1:xspan[1][2], d) -p_sigma_test = rand(p_domain.p_sigma[1]:dps.p_sigma:p_domain.p_sigma[2], 1, 1) +p_sigma_test = rand(p_domain.p_sigma[1]:(dps.p_sigma):p_domain.p_sigma[2], 1, 1) t_test = rand(tspan[1]:dt:tspan[2], 1, 1) p_mu_test = nothing p_phi_test = nothing ``` + ```@example nnparamkolmogorov sol.ufuns(x_test, t_test, p_sigma_test, p_mu_test, p_phi_test) -``` \ No newline at end of file +``` diff --git a/docs/src/tutorials/nnstopping.md b/docs/src/tutorials/nnstopping.md index 486fa8e..9ff814e 100644 --- a/docs/src/tutorials/nnstopping.md +++ b/docs/src/tutorials/nnstopping.md @@ -3,11 +3,14 @@ ## Solving for optimal strategy and expected payoff of a Bermudan Max-Call option We will calculate optimal strategy for Bermudan Max-Call option with following drift, diffusion and payoff: + ```math μ(x) =(r − δ) x, σ(x) = β diag(x1, ... , xd),\\ g(t, x) = e^{-rt}max\lbrace max\lbrace x1, ... , xd \rbrace − K, 0\rbrace ``` + We define the parameters, drift function and the diffusion function for the dynamics of the option. + ```@example nnstopping using HighDimPDE, Flux, StochasticDiffEq d = 3 # Number of assets in the stock @@ -27,28 +30,33 @@ K = 100.00 # strike price function g(x, t) return exp(-r * t) * (max(maximum(x) - K, 0)) end - ``` + We then define a [`ParabolicPDEProblem`](@ref) with no non linear term: ```@example nnstopping prob = ParabolicPDEProblem(f, sigma, u0, tspan; payoff = g) ``` -!!! note - We provide the payoff function with a keyword argument `payoff` + +!!! note +We provide the payoff function with a keyword argument `payoff` And now we define our models: + ```@example nnstopping models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoid)) for i in 1:N] ``` -!!! note - The number of models should be equal to the time discritization. + +!!! note +The number of models should be equal to the time discritization. And finally we define our optimizer and algorithm, and call `solve`: + ```@example nnstopping opt = Flux.Optimisers.Adam(0.01) alg = NNStopping(models, opt) -sol = solve(prob, alg, SRIW1(); dt = dt, trajectories = 1000, maxiters = 1000, verbose = true) +sol = solve( + prob, alg, SRIW1(); dt = dt, trajectories = 1000, maxiters = 1000, verbose = true) ``` diff --git a/paper/paper.md b/paper/paper.md index bc66b02..c46c4d2 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -28,18 +28,21 @@ aas-journal: Astrophysical Journal <- The name of the AAS journal. --- # Summary + `HighDimPDE.jl` is a Julia [@Bezanson2017] package that implements solver algorithms to solve highly dimensional non-local non-linear Partial Differential Equations (PDEs). The solver algorithms provided break down the curse of dimensionality, with a computational complexity that only grows polynomially in the number of dimension of the PDE. It is an open-source project hosted on GitHub and distributed under the MIT license. The package is designed with a user-friendly interface, provides both CPUs and GPUs support, and is integrated within the Sci-ML[@SciML] ecosystem. # Statement of need -Non-local nonlinear Partial Differential Equations arise in a variety of scientific domains including physics, engineering, finance and biology. In biology, they are for instance used for modelling the evolution of biological populations that are phenotypically and physically structured. The dimension of the PDEs can be large, corresponding to the number of phenotypic traits and physical dimensions considered. Highly dimensional PDE's cannot be solved with standard numerical methods as their computational cost increases exponentially in the number of dimensions, a problem commonly referred as the curse of dimensionality. + +Non-local nonlinear Partial Differential Equations arise in a variety of scientific domains including physics, engineering, finance and biology. In biology, they are for instance used for modelling the evolution of biological populations that are phenotypically and physically structured. The dimension of the PDEs can be large, corresponding to the number of phenotypic traits and physical dimensions considered. Highly dimensional PDE's cannot be solved with standard numerical methods as their computational cost increases exponentially in the number of dimensions, a problem commonly referred as the curse of dimensionality. # Solving PDEs with HighDimPDE.jl + HighDimPDE.jl can solve for PDEs of the form $$ \begin{aligned} - (\partial_t u)(t,x) &= \int_{\Omega} f\big(t,x,{\bf x}, u(t,x),u(t,{\bf x}), ( \nabla_x u )(t,x ),( \nabla_x u )(t,{\bf x} ) \big) \, \nu_x(d{\bf x}) \\ - & \quad + \big\langle \mu(t,x), ( \nabla_x u )( t,x ) \big\rangle + \tfrac{1}{2} \text{Trace} \big(\sigma(t,x) [ \sigma(t,x) ]^* ( \text{Hess}_x u)(t, x ) \big). \tag{1} +(\partial_t u)(t,x) &= \int_{\Omega} f\big(t,x,{\bf x}, u(t,x),u(t,{\bf x}), ( \nabla_x u )(t,x ),( \nabla_x u )(t,{\bf x} ) \big) \, \nu_x(d{\bf x}) \\ +& \quad + \big\langle \mu(t,x), ( \nabla_x u )( t,x ) \big\rangle + \tfrac{1}{2} \text{Trace} \big(\sigma(t,x) [ \sigma(t,x) ]^* ( \text{Hess}_x u)(t, x ) \big). \tag{1} \end{aligned} $$ where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subset \R^d$ is a function subject to initial conditions $u(0,x) = g(x)$ and Neumann Boundary conditions. @@ -48,7 +51,8 @@ where $u \colon [0,T] \times \Omega \to \R$, $\Omega \subset \R^d$ is a function ## Algorithm overview ----------------------------------------------- +* * * + Features | `DeepSplitting` | `MLP` | ----------|:----------------------:|:------------: Time discretization free| ❌ | ✅ | @@ -66,22 +70,24 @@ The `DeepSplitting`[@Beck2019] algorithm reformulates the PDE as a stochastic le `DeepSplitting` relies on two main ideas: -- The approximation of the solution $u$ by a parametric function $\bf u^\theta$. + - The approximation of the solution $u$ by a parametric function $\bf u^\theta$. -- The training of $\bf u^\theta$ by simulated stochastic trajectories of particles, with the help of the Machine-Learning library [@Flux]. + - The training of $\bf u^\theta$ by simulated stochastic trajectories of particles, with the help of the Machine-Learning library [@Flux]. ## `MLP` -The `MLP`[@Becker2020], for Multi-Level Picard iterations, reformulates the PDE problem as a fixed point equation through the Feynman Kac formula. + +The `MLP`[@Becker2020], for Multi-Level Picard iterations, reformulates the PDE problem as a fixed point equation through the Feynman Kac formula. `MLP` relies on two main principles to solve the fixed point equation: -- [Picard iterations](https://en.wikipedia.org/wiki/Picard–Lindelöf_theorem), and + - [Picard iterations](https://en.wikipedia.org/wiki/Picard%E2%80%93Lindel%C3%B6f_theorem), and -- a [Multilevel Monte Carlo](https://en.wikipedia.org/wiki/Multilevel_Monte_Carlo_method) approach to reduce the complexity of the numerical approximation of the time integral in the fixed point equation. + - a [Multilevel Monte Carlo](https://en.wikipedia.org/wiki/Multilevel_Monte_Carlo_method) approach to reduce the complexity of the numerical approximation of the time integral in the fixed point equation. # Examples ## The `DeepSplitting` algorithm + Consider the 5-dimensional replicator mutator equation [@Hamel:2019] $$ \partial_t u = u (m - \int_{\R^5} m(y)u(y,t)dy) + \frac{1}{2}\sigma^2\Delta_xu \tag{2} @@ -90,21 +96,21 @@ where $$ m(x) = -\frac{1}{2}||x|| $$ -and +and $$ u(x,0) = \mathcal{N}_{0,0.05}(x) $$ where $\mathcal{N}_{\mu,\sigma}$ is the normal distribution with mean $\mu$ and standard distribution $\sigma$. ```julia -tspan = (0f0,15f-2) -dt = 5f-2 # time step -μ(x, p, t) = 0f0 # advection coefficients -σ(x, p, t) = 1f-1 #1f-1 # diffusion coefficients -ss0 = 5f-2#std g0 +tspan = (0.0f0, 15.0f-2) +dt = 5.0f-2 # time step +μ(x, p, t) = 0.0f0 # advection coefficients +σ(x, p, t) = 1.0f-1 #1f-1 # diffusion coefficients +ss0 = 5.0f-2#std g0 d = 5 -U = 25f-2 +U = 25.0f-2 x0_sample = (fill(-U, d), fill(U, d)) batch_size = 10000 @@ -114,62 +120,71 @@ K = 1 hls = d + 50 #hidden layer size nn_batch = Flux.Chain( - Dense(d, hls, tanh), - Dense(hls, hls, tanh), - Dense(hls, 1, x->x^2)) # positive function + Dense(d, hls, tanh), + Dense(hls, hls, tanh), + Dense(hls, 1, x -> x^2)) # positive function opt = ADAM(1e-2)#optimiser -alg = DeepSplitting(nn_batch, K=K, opt = opt, mc_sample = UniformSampling(x0_sample[1], x0_sample[2]) ) - -g(x) = Float32((2*π)^(-d/2)) * ss0^(- Float32(d) * 5f-1) * exp.(-5f-1 *sum(x .^2f0 / ss0, dims = 1)) # initial condition -m(x) = - 5f-1 * sum(x.^2, dims=1) +alg = DeepSplitting( + nn_batch, K = K, opt = opt, mc_sample = UniformSampling(x0_sample[1], x0_sample[2])) + +function g(x) + Float32((2 * π)^(-d / 2)) * ss0^(-Float32(d) * 5.0f-1) * + exp.(-5.0f-1 * sum(x .^ 2.0f0 / ss0, dims = 1)) +end # initial condition +m(x) = -5.0f-1 * sum(x .^ 2, dims = 1) vol = prod(x0_sample[2] - x0_sample[1]) -f(y, z, v_y, v_z, p, t) = max.(v_y, 0f0) .* (m(y) .- vol * max.(v_z, 0f0) .* m(z)) # nonlocal nonlinear part of the +f(y, z, v_y, v_z, p, t) = max.(v_y, 0.0f0) .* (m(y) .- vol * max.(v_z, 0.0f0) .* m(z)) # nonlocal nonlinear part of the # defining the problem prob = PIDEProblem(μ, σ, tspan, g, f, - x0_sample = x0_sample - ) + x0_sample = x0_sample +) # solving -xgrid,ts,sol = solve(prob, - alg, - dt, - verbose = false, - abstol = 1f-3, - maxiters = train_steps, - batch_size = batch_size, - use_cuda = true - ) +xgrid, ts, sol = solve(prob, + alg, + dt, + verbose = false, + abstol = 1.0f-3, + maxiters = train_steps, + batch_size = batch_size, + use_cuda = true +) u1 = [sol[end](x)[] for x in xgrid] ``` + ![Solution to Eq. (2) obtained with `DeepSplitting`](./hamel_5d.png){ width=20% } ## The `MLP` algorithm + Consider the 5-dimensional non-local Fisher KPP PDE + ```math \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]^5$, and assume Neumann Boundary condition on $\Omega$. + ```julia using HighDimPDE ## Definition of the problem d = 10 # dimension of the problem -tspan = (0.0,0.5) # time horizon -x0 = fill(0.,d) # initial point -g(x) = exp( -sum(x.^2) ) # initial condition +tspan = (0.0, 0.5) # time horizon +x0 = fill(0.0, d) # initial point +g(x) = exp(-sum(x .^ 2)) # initial condition μ(x, p, t) = 0.0 # advection coefficients σ(x, p, t) = 0.1 # diffusion coefficients -x0_sample = [-1/2, 1/2] -f(x, y, v_x, v_y, ∇v_x, ∇v_y, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) -prob = PIDEProblem(μ, - σ, x0, tspan, g, f, - x0_sample = x0_sample) # defining x0_sample is sufficient to implement Neumann boundary conditions +x0_sample = [-1 / 2, 1 / 2] +f(x, y, v_x, v_y, ∇v_x, ∇v_y, t) = max(0.0, v_x) * (1 - max(0.0, v_y)) +prob = PIDEProblem(μ, + σ, x0, tspan, g, f, + x0_sample = x0_sample) # defining x0_sample is sufficient to implement Neumann boundary conditions ## Definition of the algorithm -alg = MLP(mc_sample = UniformSampling(x0_sample[1], x0_sample[2]) ) +alg = MLP(mc_sample = UniformSampling(x0_sample[1], x0_sample[2])) -sol = solve(prob, alg, multithreading=true) +sol = solve(prob, alg, multithreading = true) ``` # Acknowledgements diff --git a/src/DeepBSDE.jl b/src/DeepBSDE.jl index 3bf25c7..7985480 100644 --- a/src/DeepBSDE.jl +++ b/src/DeepBSDE.jl @@ -155,12 +155,12 @@ function DiffEqBase.solve(prob::ParabolicPDEProblem, function neural_sde(init_cond) sde_prob = remake(sde_prob, u0 = init_cond) ensemble_prob = EnsembleProblem(sde_prob) - sol = solve(ensemble_prob, sdealg, EnsembleSerial(); - u0 = init_cond, trajectories = trajectories, dt = dt, p = p3, - sensealg = SciMLSensitivity.TrackerAdjoint(), + sol = solve(ensemble_prob, sdealg, EnsembleSerial(); + u0 = init_cond, trajectories = trajectories, dt = dt, p = p3, + sensealg = SciMLSensitivity.TrackerAdjoint(), save_everystep = false, kwargs...) - map(sol) do _sol + map(sol) do _sol predict_ans = Array(_sol) (predict_ans[1:(end - 1), end], predict_ans[end, end]) end @@ -290,7 +290,7 @@ function DiffEqBase.solve(prob::ParabolicPDEProblem, f_matrix = give_f_matrix(X, u_domain, _σᵀ∇u, p, ts[i]) a_ = A[findmax(collect(A) .* u .- collect(legendre_transform(f_matrix, a, u_domain) - for a in A))[2]] + for a in A))[2]] I = I + a_ * dt Q = Q + exp(I) * legendre_transform(f_matrix, a_, u_domain) end diff --git a/src/DeepBSDE_Han.jl b/src/DeepBSDE_Han.jl index 75744c8..29bcb99 100644 --- a/src/DeepBSDE_Han.jl +++ b/src/DeepBSDE_Han.jl @@ -22,7 +22,7 @@ function DiffEqBase.solve(prob::ParabolicPDEProblem, limits = false, trajectories_upper = 1000, trajectories_lower = 1000, - maxiters_limits = 10,) + maxiters_limits = 10) X0 = prob.x ts = prob.tspan[1]:dt:prob.tspan[2] d = length(X0) @@ -144,7 +144,7 @@ function DiffEqBase.solve(prob::ParabolicPDEProblem, f_matrix = give_f_matrix(X, u_domain, _σᵀ∇u, p, ts[i]) a_ = A[findmax(collect(A) .* u .- collect(legendre_transform(f_matrix, a, u_domain) - for a in A))[2]] + for a in A))[2]] I = I + a_ * dt Q = Q + exp(I) * legendre_transform(f_matrix, a_, u_domain) end diff --git a/src/DeepSplitting.jl b/src/DeepSplitting.jl index 0722ee9..4cc5653 100644 --- a/src/DeepSplitting.jl +++ b/src/DeepSplitting.jl @@ -1,7 +1,7 @@ _copy(t::Tuple) = t _copy(t) = t -function _copy(opt::O) where O<:Flux.Optimise.AbstractOptimiser - return O([_copy(getfield(opt,f)) for f in fieldnames(typeof(opt))]...) +function _copy(opt::O) where {O <: Flux.Optimise.AbstractOptimiser} + return O([_copy(getfield(opt, f)) for f in fieldnames(typeof(opt))]...) end """ @@ -44,7 +44,7 @@ function DeepSplitting(nn; λs::L = nothing, mc_sample = NoSampling()) where { O <: Flux.Optimise.AbstractOptimiser, - L <: Union{Nothing, Vector{N}} where {N <: Number}, + L <: Union{Nothing, Vector{N}} where {N <: Number} } isnothing(λs) ? λs = [opt.eta] : nothing DeepSplitting(nn, K, opt, λs, mc_sample) @@ -101,7 +101,7 @@ function DiffEqBase.solve(prob::Union{PIDEProblem, ParabolicPDEProblem}, g, μ, σ, p = prob.g, prob.μ, prob.σ, prob.p f = if isa(prob, ParabolicPDEProblem) - (y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) -> prob.f(y, v_y, ∇v_y, p, t ) + (y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) -> prob.f(y, v_y, ∇v_y, p, t) else prob.f end diff --git a/src/HighDimPDE.jl b/src/HighDimPDE.jl index 4f3a525..8987833 100644 --- a/src/HighDimPDE.jl +++ b/src/HighDimPDE.jl @@ -70,61 +70,58 @@ with `` u(x,0) = g(x)``. * `neumann_bc`: if provided, Neumann boundary conditions on the hypercube `neumann_bc[1] × neumann_bc[2]`. """ function PIDEProblem(μ, - σ, - x0::Union{Nothing, AbstractArray}, - tspan::TF, - g, - f; - p::Union{Nothing, AbstractVector} = nothing, - x0_sample::Union{Nothing, AbstractSampling} = NoSampling(), - neumann_bc::Union{Nothing, AbstractVector} = nothing, - kw...) where {TF <: Tuple{AbstractFloat, AbstractFloat}} - - + σ, + x0::Union{Nothing, AbstractArray}, + tspan::TF, + g, + f; + p::Union{Nothing, AbstractVector} = nothing, + x0_sample::Union{Nothing, AbstractSampling} = NoSampling(), + neumann_bc::Union{Nothing, AbstractVector} = nothing, + kw...) where {TF <: Tuple{AbstractFloat, AbstractFloat}} isnothing(neumann_bc) ? nothing : @assert eltype(eltype(neumann_bc)) <: eltype(x0) @assert(eltype(f(x0, x0, g(x0), g(x0), x0, x0, p, tspan[1]))==eltype(x0), - "Type returned by non linear function `f` must match the type of `x0`") + "Type returned by non linear function `f` must match the type of `x0`") @assert eltype(g(x0))==eltype(x0) "Type of `g(x)` must match the Type of x" PIDEProblem{typeof(g(x0)), - typeof(g), - typeof(f), - typeof(μ), - typeof(σ), - typeof(x0), - eltype(tspan), - typeof(p), - typeof(x0_sample), - typeof(neumann_bc), - typeof(kw)}(g(x0), - g, - f, - μ, - σ, - x0, - tspan, - p, - x0_sample, - neumann_bc, - kw) - + typeof(g), + typeof(f), + typeof(μ), + typeof(σ), + typeof(x0), + eltype(tspan), + typeof(p), + typeof(x0_sample), + typeof(neumann_bc), + typeof(kw)}(g(x0), + g, + f, + μ, + σ, + x0, + tspan, + p, + x0_sample, + neumann_bc, + kw) end struct ParabolicPDEProblem{uType, G, F, Mu, Sigma, xType, tType, P, UD, NBC, K} <: - DiffEqBase.AbstractODEProblem{uType, tType, false} - u0::uType - g::G # initial condition - f::F # nonlinear part - μ::Mu - σ::Sigma - x::xType - tspan::Tuple{tType, tType} - p::P - x0_sample::UD # the domain of u to be solved - neumann_bc::NBC # neumann boundary conditions - kwargs::K + DiffEqBase.AbstractODEProblem{uType, tType, false} + u0::uType + g::G # initial condition + f::F # nonlinear part + μ::Mu + σ::Sigma + x::xType + tspan::Tuple{tType, tType} + p::P + x0_sample::UD # the domain of u to be solved + neumann_bc::NBC # neumann boundary conditions + kwargs::K end """ @@ -185,7 +182,7 @@ function ParabolicPDEProblem(μ, @assert !isnothing(x0)||!isnothing(xspan) "Either of `x0` or `xspan` must be provided." !isnothing(f) && @assert(eltype(f(x0, eltype(x0)(0.0), x0, p, tspan[1]))==eltype(x0), - "Type of non linear function `f(x)` must type of x") + "Type of non linear function `f(x)` must type of x") # Wrap kwargs : kw = NamedTuple(kw) @@ -270,7 +267,8 @@ include("NNStopping.jl") include("NNKolmogorov.jl") include("NNParamKolmogorov.jl") -export PIDEProblem, ParabolicPDEProblem, PIDESolution, DeepSplitting, DeepBSDE, MLP, NNStopping +export PIDEProblem, ParabolicPDEProblem, PIDESolution, DeepSplitting, DeepBSDE, MLP, + NNStopping export NNKolmogorov, NNParamKolmogorov export NormalSampling, UniformSampling, NoSampling, solve end diff --git a/src/MLP.jl b/src/MLP.jl index c640649..843119e 100644 --- a/src/MLP.jl +++ b/src/MLP.jl @@ -32,7 +32,7 @@ Returns a `PIDESolution` object. function DiffEqBase.solve(prob::Union{PIDEProblem, ParabolicPDEProblem}, alg::MLP; multithreading = true, - verbose = false,) + verbose = false) # unbin stuff x = prob.x @@ -41,9 +41,9 @@ function DiffEqBase.solve(prob::Union{PIDEProblem, ParabolicPDEProblem}, M = alg.M L = alg.L mc_sample! = alg.mc_sample! - g= prob.g + g = prob.g f = if isa(prob, ParabolicPDEProblem) - (y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) -> prob.f(y, v_y, ∇v_y, p, t ) + (y, z, v_y, v_z, ∇v_y, ∇v_z, p, t) -> prob.f(y, v_y, ∇v_y, p, t) else prob.f end @@ -197,7 +197,7 @@ function _ml_picard_mlt(M, # monte carlo integration # distributing tasks NUM_THREADS = Threads.nthreads() tasks = [Threads.@spawn(_ml_picard_call(M, L, K, x, s, t, mc_sample!, g, f, verbose, - NUM_THREADS, thread_id, prob, neumann_bc)) for thread_id in 1:NUM_THREADS] + NUM_THREADS, thread_id, prob, neumann_bc)) for thread_id in 1:NUM_THREADS] # first level num = M^(L) diff --git a/test/DeepBSDE.jl b/test/DeepBSDE.jl index c1cd31c..343076f 100644 --- a/test/DeepBSDE.jl +++ b/test/DeepBSDE.jl @@ -25,7 +25,7 @@ end g(X) = sum(X .^ 2) # terminal condition f(X, u, σᵀ∇u, p, t) = eltype(X)(0.0) - μ_f(X, p, t) = X*0.0f0 #Vector d x 1 + μ_f(X, p, t) = X * 0.0f0 #Vector d x 1 σ_f(X, p, t) = Diagonal(ones(eltype(X), d)) |> Matrix #Matrix d x d prob = ParabolicPDEProblem(μ_f, σ_f, x0, tspan; g, f) @@ -227,7 +227,7 @@ end W() = randn(d, 1) u_analytical(x, t) = -(1 / λ) * log(mean(exp(-λ * g(x .+ sqrt(2.0) * abs.(T - t) .* W())) - for _ in 1:MC)) + for _ in 1:MC)) analytical_sol = u_analytical(x0, tspan[1]) error_l2 = rel_error_l2(sol.us, analytical_sol) println("error_l2 = ", error_l2, "\n") diff --git a/test/DeepBSDE_Han.jl b/test/DeepBSDE_Han.jl index f915f64..ce84ba0 100644 --- a/test/DeepBSDE_Han.jl +++ b/test/DeepBSDE_Han.jl @@ -38,8 +38,8 @@ end Dense(hls, 1)) # sub-neural network approximating the spatial gradients at time point σᵀ∇u = [Flux.Chain(Dense(d, hls, relu), - Dense(hls, hls, relu), - Dense(hls, d)) for i in 1:time_steps] + Dense(hls, hls, relu), + Dense(hls, d)) for i in 1:time_steps] alg = DeepBSDE(u0, σᵀ∇u, opt = opt) @@ -82,8 +82,8 @@ end Dense(hls, 1)) # sub-neural network approximating the spatial gradients at time point σᵀ∇u = [Flux.Chain(Dense(d, hls, relu), - Dense(hls, hls, relu), - Dense(hls, d)) for i in 1:time_steps] + Dense(hls, hls, relu), + Dense(hls, d)) for i in 1:time_steps] alg = DeepBSDE(u0, σᵀ∇u, opt = opt) @@ -125,9 +125,9 @@ end Dense(hls, hls, relu), Dense(hls, 1)) σᵀ∇u = [Flux.Chain(Dense(d, hls, relu), - Dense(hls, hls, relu), - Dense(hls, hls, relu), - Dense(hls, d)) for i in 1:time_steps] + Dense(hls, hls, relu), + Dense(hls, hls, relu), + Dense(hls, d)) for i in 1:time_steps] alg = DeepBSDE(u0, σᵀ∇u, opt = opt) @@ -170,8 +170,8 @@ end # sub-neural network approximating the spatial gradients at time point σᵀ∇u = [Flux.Chain(Dense(d, hls, relu), - Dense(hls, hls, relu), - Dense(hls, d)) for i in 1:time_steps] + Dense(hls, hls, relu), + Dense(hls, d)) for i in 1:time_steps] alg = DeepBSDE(u0, σᵀ∇u, opt = opt) @@ -216,9 +216,9 @@ end # sub-neural network approximating the spatial gradients at time point σᵀ∇u = [Flux.Chain(Dense(d, hls, relu), - Dense(hls, hls, relu), - Dense(hls, hls, relu), - Dense(hls, d)) for i in 1:time_steps] + Dense(hls, hls, relu), + Dense(hls, hls, relu), + Dense(hls, d)) for i in 1:time_steps] alg = DeepBSDE(u0, σᵀ∇u, opt = opt) @@ -235,7 +235,7 @@ end W() = randn(d, 1) u_analytical(x, t) = -(1 / λ) * log(mean(exp(-λ * g(x .+ sqrt(2.0) * abs.(T - t) .* W())) - for _ in 1:MC)) + for _ in 1:MC)) analytical_sol = u_analytical(x0, tspan[1]) error_l2 = rel_error_l2(sol.us, analytical_sol) @@ -288,9 +288,9 @@ end Dense(hls, 1)) σᵀ∇u = [Flux.Chain(Dense(d, hls, relu), - Dense(hls, hls, relu), - Dense(hls, hls, relu), - Dense(hls, d)) for i in 1:time_steps] + Dense(hls, hls, relu), + Dense(hls, hls, relu), + Dense(hls, d)) for i in 1:time_steps] alg = DeepBSDE(u0, σᵀ∇u, opt = opt) sol = solve(prob, diff --git a/test/DeepSplitting.jl b/test/DeepSplitting.jl index fd8b393..14f67bf 100644 --- a/test/DeepSplitting.jl +++ b/test/DeepSplitting.jl @@ -4,7 +4,7 @@ using Test using Flux using Statistics using CUDA -using Random +using Random Random.seed!(100) if CUDA.functional() @@ -376,7 +376,7 @@ if false W() = randn(d, 1) u_analytical(x, t) = -(1 / λ) * log(mean(exp(-λ * g(x .+ sqrt(2.0f0) * abs.(T - t) .* W())) - for _ in 1:MC)) + for _ in 1:MC)) hls = d + 50 #hidden layer size diff --git a/test/MLP.jl b/test/MLP.jl index c242e5c..0a33899 100644 --- a/test/MLP.jl +++ b/test/MLP.jl @@ -3,7 +3,7 @@ using Random using Test using Statistics # using Revise -using Random +using Random Random.seed!(100) #relative error l2 function rel_error_l2(u, uanal) diff --git a/test/NNParamKolmogorov.jl b/test/NNParamKolmogorov.jl index 06377bc..36d7c0c 100644 --- a/test/NNParamKolmogorov.jl +++ b/test/NNParamKolmogorov.jl @@ -52,7 +52,8 @@ function analytical(x, t, y) return x .^ 2 .+ t .* (y .* y) end -preds = map((i) -> sol.ufuns(x_test[:, :, i], +preds = map( + (i) -> sol.ufuns(x_test[:, :, i], t_test[:, i], γ_sigma_test[:, :, :, i], nothing, From d06520409e63045f48af9e9accfba0160652a811 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Tue, 2 Jul 2024 16:41:28 +0530 Subject: [PATCH 07/11] build: bump compat for SciMLSensitivity --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 2ccbfa4..39af627 100644 --- a/Project.toml +++ b/Project.toml @@ -12,6 +12,7 @@ Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" From 6332e2cd035897fba07ec0ef86ccd237947f2f14 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sun, 7 Jul 2024 14:44:14 +0530 Subject: [PATCH 08/11] build: add compat for SciMLBase --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 39af627..5cdb0db 100644 --- a/Project.toml +++ b/Project.toml @@ -34,6 +34,7 @@ LinearAlgebra = "1.10" Random = "1.10" Reexport = "1.2.2" SafeTestsets = "0.1" +SciMLBase = "2.42" SciMLSensitivity = "7.62" SparseArrays = "1.10" Statistics = "1.10" From 80de3c3d07fc6bda3ac40997dd3e3f76811d13a2 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sun, 7 Jul 2024 14:46:27 +0530 Subject: [PATCH 09/11] docs: update Project --- docs/Project.toml | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index f02b906..a3e379d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,6 @@ [deps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +cuDNN = "02a925ec-e4fe-4b08-9a7e-0d78e3d38ccd" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" HighDimPDE = "57c578d5-59d4-4db8-a490-a9fc372d19d2" @@ -7,6 +8,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" [compat] -CUDA = "5" +CUDA = "5.4.2" +cuDNN = "1.3.2" Documenter = "1" -Flux = "0.13, 0.14" +Flux = "0.14.16" +StochasticDiffEq = "6.66" From e202998e429681cb80afa06d9debb59cf3fed985 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sun, 7 Jul 2024 15:02:07 +0530 Subject: [PATCH 10/11] chore: fix formatting --- docs/src/DeepSplitting.md | 2 ++ docs/src/problems.md | 2 ++ docs/src/tutorials/nnstopping.md | 4 ++++ 3 files changed, 8 insertions(+) diff --git a/docs/src/DeepSplitting.md b/docs/src/DeepSplitting.md index cf0d2bd..9ddee70 100644 --- a/docs/src/DeepSplitting.md +++ b/docs/src/DeepSplitting.md @@ -112,6 +112,8 @@ u(t_{n+1}, X_{T - t_{n+1}}) \approx \sum_{j=1}^{\text{batch\_size}} \left[ u(t_{ ``` !!! tip + + In practice, if you have a non-local model, you need to provide the sampling method and the number $K$ of MC integration through the keywords `mc_sample` and `K`. `julia alg = DeepSplitting(nn, opt = opt, mc_sample = mc_sample, K = 1) ` `mc_sample` can be whether `UniformSampling(a, b)` or ` NormalSampling(σ_sampling, shifted)`. diff --git a/docs/src/problems.md b/docs/src/problems.md index ebff55d..43369cd 100644 --- a/docs/src/problems.md +++ b/docs/src/problems.md @@ -4,5 +4,7 @@ ParabolicPDEProblem ``` !!! note + + While choosing to define a PDE using `PIDEProblem`, note that the function being integrated `f` is a function of `f(x, y, v_x, v_y, ∇v_x, ∇v_y)` out of which `y` is the integrating variable and `x` is constant throughout the integration. If a PDE has no integral and the non linear term `f` is just evaluated as `f(x, v_x, ∇v_x)` then we suggest using `ParabolicPDEProblem` diff --git a/docs/src/tutorials/nnstopping.md b/docs/src/tutorials/nnstopping.md index 9ff814e..e057a32 100644 --- a/docs/src/tutorials/nnstopping.md +++ b/docs/src/tutorials/nnstopping.md @@ -39,6 +39,8 @@ prob = ParabolicPDEProblem(f, sigma, u0, tspan; payoff = g) ``` !!! note + + We provide the payoff function with a keyword argument `payoff` And now we define our models: @@ -49,6 +51,8 @@ models = [Chain(Dense(d + 1, 32, tanh), BatchNorm(32, tanh), Dense(32, 1, sigmoi ``` !!! note + + The number of models should be equal to the time discritization. And finally we define our optimizer and algorithm, and call `solve`: From 9dffdc1824c6dfc67db769719f9939c0271a3632 Mon Sep 17 00:00:00 2001 From: ashutosh-b-b Date: Sun, 7 Jul 2024 15:41:16 +0530 Subject: [PATCH 11/11] docs: fix use_cuda --- docs/src/getting_started.md | 2 +- docs/src/tutorials/deepsplitting.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index a52f073..256d24e 100644 --- a/docs/src/getting_started.md +++ b/docs/src/getting_started.md @@ -119,7 +119,7 @@ sol = solve(prob, [`DeepSplitting`](@ref deepsplitting) can run on the GPU for (much) improved performance. To do so, just set `use_cuda = true`. -```@example DeepSplitting_non_local_PDE +```julia sol = solve(prob, alg, 0.1, diff --git a/docs/src/tutorials/deepsplitting.md b/docs/src/tutorials/deepsplitting.md index 73701da..6170dce 100644 --- a/docs/src/tutorials/deepsplitting.md +++ b/docs/src/tutorials/deepsplitting.md @@ -59,7 +59,7 @@ sol = solve(prob, `DeepSplitting` can run on the GPU for (much) improved performance. To do so, just set `use_cuda = true`. -```@example deepsplitting +```julia sol = solve(prob, alg, 0.1,