Interfacing with 3rd-party PDE libraries

[yizhang-yiz]

Summary:

Allow 3rd-party PDE(partial differential equation) libraries to be used to perform inference that involve PDEs.

Description:

The design involves cmdstan, stan, and math. This issue ticket solicits the discussion regarding the related design in all the repos. See stan-dev/math#931 for corresponding issue for Math library.

Design goals

• Stan interface with any 3rd-party PDE library, so that QoI(quantity of interest, PDE jargon) can be connected to input parameters.
• Use external library in Stan for PDE computing. The external library is provided by the user through user_header and Stan's external C++ interface.
• Hide stan internals from PDE developers.

This design allows the user to use his own PDE solver in Stan, by using solve_pde function and make with STANCFLAGS=--allow_undefined. Similar to ODE solver in Stan, the PDE is provided through a functor(Stan's user-defined function). The actual action of the functor is provided through the external PDE solver in user_header.hpp(see example below). In this design the user is asked to do the following to use his external PDE solver:

• write a stan user-function, in the following example, this function is laplace_model.
• declare the external functions used in the above stan user-function.
• in Stan call solve_pde, with the above stan user-function as the first argument.
• provide Makefile for the external PDE solver
• build the Stan model with

make WITH_EXTERN_PDE=1 EXTERN_PDE_MAKEFILE=/path/to/makefile stan_model

An additional design goal is to allow stan to coordinate MPI computation with the external lib. This belongs to next stage, will be not be discussed in the rest of the discussion.

Also note that the design works on any 3rd-party library that can do sensitivity analysis, but my interest is mostly in PDE.

The rest subsections introduce the proposed design.

forward_pde solve_pde function

Implementation here
https://github.com/yizhang-cae/math/tree/forward_pde

forward_pde solve_pde in math maps input PDE model functor and parameters to QoI:

  template<typename F_pde, typename T>
  inline std::vector<T> solve_pde(const F_pde& pde,
                                         const std::vector<T>& theta,
                                         const std::vector<double>& x_r,
                                         const std::vector<int>& x_i,
                                         std::ostream* msgs = nullptr);

This is similar to ODE solver, except that all model-related information is provide by functor pde, which as the following signature

  inline std::vector<std::vector<double> >
  operator()(const std::vector<double>& theta,
             const bool need_sens,
             const std::vector<double>& x_r,
             const std::vector<int>& x_i,
             std::ostream* msgs = nullptr);

The rationale is to relieve PDE lib user from working on vars. He needs to make decision on

• Whether to compute sensitivity(need_sens?),
• How the sensitivity is computed from theta,
  and doesn't need to know how the result is incorporated into math.

Linking forward_pde to Stan language

Implementation is at
https://github.com/yizhang-cae/stan/tree/forward_pde
Nothing new here, just exposing a high-order function.

Using forward_pde solve_pde in Stan would be something like the following.

functions {
  real[,] solve_with_sensitivity(real[] theta);
  real[,] solve(real[] theta);
  real[,] laplace_model(real[] theta, int need_sens, real[] x_r, int[] x_i){
    if(need_sens)
      return solve_with_sensitivity(theta);
    else
      return solve(theta);
    }
}
/* .... */
parameters {
  real<lower = 0> k[2];
} 

transformed parameters{
  real QoI[1];
  QoI = solve_pde(laplace_model, k, x_r, x_i);
}
/* .... */

make in cmdstan

In order to make the above model, incmdstan two switches are added to Makefile

• WITH_EXTERN_PDE: =1 indicates user needs to supply user_header.hpp.
• EXTERN_PDE_MAKEFILE: with WITH_EXTERN_PDE=1 user needs to provide make info regarding his library, so that the above Laplace model can be built with

make WITH_EXTERN_PDE=1 EXTERN_PDE_MAKEFILE=./examples/forward_pde_laplace/Makefile ./examples/forward_pde_laplace/forward_pde_laplace

The example and the Makefile change can be found at
https://github.com/yizhang-cae/cmdstan/tree/forward_pde

Reproducible Steps:

Final model in ./examples/forward_pde_lapace
https://github.com/yizhang-cae/cmdstan/tree/forward_pde
Unit tests in the above math branch.

Current Output:

n/a

Expected Output:

Three external libraries has been tested and included in unit tests. To verify you need to install these libraries and provide make info indicated in the Makefile of the unit test directories in math repo.

• http://libmesh.github.io/
• http://mfem.org/
• http://opensees.berkeley.edu/

Pressure contour of the Darcy's flow from the MFEM unit test

Parameter(normalized porosity k) from Stan, based on simulated data.

Additional Information:

Provide any additional information here.

Current Version:

v2.17.1

[bgoodri]

The rstan process for including external files is a bit different although slightly easier. We just need to make people who are using this mechanism aware of it.

[syclik]

What's the motivation for the name `forward_pde`?

…

On Mon, Jul 02, 2018 at 11:17 PM, ***@***.***>wrote: The *rstan* process for including external files is a bit different although slightly easier. We just need to make people who are using this mechanism aware of it. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#2567 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAZ_F2AWRtQIqbke7kCYv4y-UPLKxMy2ks5uCuJigaJpZM4VALuQ> .

[yizhang-yiz]

What's the motivation for the name forward_pde?

It solves a forward pde problem, compared to inverse problem.

[syclik]

Got it. Doesn't the literature and other software just call that "solving a pde"? I did a quick search and it looks like most other software functions name the function "pde" or something similar without calling it "forward." Are we planning on tackling inverse pdes soon?

…

On Thu, Jul 05, 2018 at 9:36 AM, Yi ***@***.***>wrote: What's the motivation for the name forward_pde? It solves a forward pde problem <https://onlinelibrary.wiley.com/doi/10.1002/9781118478141.ch3>, compared to inverse problem. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#2567 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAZ_F3NLHmzh3Gk9qhCC3YPrzg2kRuYjks5uDhZZgaJpZM4VALuQ> .

[yizhang-yiz]

I want to emphasize we're "forwarding" the parameter to QoI.

Essentially we ARE solving inverse problem (using data to trace back to system parameters), by solving forward problem repeatedly. But if by "inverse" you mean inference of parameter in an infinite function space like in tomography, then probably not. But it's more of a PDE solver issue than a Stan issue.

[syclik]

What's QoI?

…

On Thu, Jul 05, 2018 at 10:33 AM, Yi ***@***.***>wrote: I want to emphasize we're "forwarding" the parameter to QoI. Essentially we ARE solving inverse problem (using data to trace back to system parameters), by solving forward problem repeatedly. But if by "inverse" you mean inference of parameter in an infinite function space like in tomography, then probably not. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#2567 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAZ_F-h-4uM7RU-rZomMSr-tTN3spu9sks5uDiPLgaJpZM4VALuQ> .

[yizhang-yiz]

Quantity of interest. Sorry, forgot that's PDE jargon. It's essentially the observable with its data collected.

[syclik]

More help? What's "observable with its data collected"? I don't envision the Stan language being a PDE-specific language, so perhaps more descriptive names would be better? Things that would show up naturally or when searching? I'm now looking at the doc for some of these packages. It looks like you're trying to abstract a lot of functionality away behind this single function. Is that the intent? Also, what do you mean by "the rationale is to relieve PDE lib user from working on vars"? I don't understand that statement. What's a "PDE lib"? At the math repo level, wouldn't they need to know about stan::math::var?

…

On Thu, Jul 05, 2018 at 10:41 AM, Yi ***@***.***>wrote: Quantity of interest. Sorry, forgot that's PDE jargon. It's essentially the observable with its data collected. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#2567 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAZ_F58Mx0z2TWFJGcTuDe0pyTUsf_fbks5uDiWugaJpZM4VALuQ> .

[yizhang-yiz]

I'm open to any other names.

It looks like you're trying to abstract a lot of functionality away behind this single function. Is that the intent?

Yes. Stan only cares about the PDE solution result and the sensitivity, so we can do autodiff regarding that result. PDE solver developers don't need to touch var, because all they need to provide through the user header are the PDE solution results and sensitivity, which are used in precomputed_gradients in forward_pde.

Things better explained in an example. Say in flow simulation around a car, we collect noisy data of drag force, the design parameter(theta) is the length of the vehicle. So we do
drag_force = forward_pde(theta). In this context, drag force is the QoI. The PDE solver developer should be familiar with using theta to compute drag, and sensitivity of drag with respect to theta, and that's what he needs to provide to Stan: he uses his flow simulation solver("PDE lib") to solve flow problem, calculate drag force and sensitivity, wrap the whole process in a function that takes theta as argument, and hand it over to Stan. He doesn't need tot know about stan::math::var. forward_pde takes care of using drag and sensitivity and precomputed_gradient() to produce a var variable.

[yizhang-yiz]

Unlike ODE, I don't think we can find a typical usage in a single PDE solver, hence the external library interface and abstracting away the PDE solution process.

As to name, will solve_pde work? I'm terrible at naming things so whatever @syclik @bob-carpenter think proper.

[bob-carpenter]

It's also Andrew's jargon and hence shows up in Stan doc all over the place.

…

On Jul 5, 2018, at 10:41 AM, Yi Zhang ***@***.***> wrote: Quantity of interest. Sorry, forgot that's PDE jargon. It's essentially the observable with its data collected.

[bob-carpenter]

This is exactly the same motivation as behind the adjoint-Jacobian formulation of reverse mode that Ben Bales is working on now. It lets the user write the value function and adjoint-Jacobian product function without every having to worry about var types. In fact, if the result is multivariate, this can be more efficient in memory than storing multiple gradients in a Jacobian, so the adjoint-Jacobian product may be worth thinking about here, too.

[bob-carpenter]

Naming's hard and almost everyone's bad at it. Hence this popular quote from Phil Karlton:

…

There are only two hard things in Computer Science: cache invalidation and naming things.

[yizhang-yiz]

The current design treats PDE solver as black box. This is feasible for small-medium size problems, For large problems one usually needs some reduced-order/surrogate model for the PDE to make it tractable. Again, this is more of a PDE solver issue instead of Stan issue.

One caveat is that along user_header.hpp users are asked to provide a pde.Makefile, which may break Stan build. We can only suggest them to make the variables there target-specific. Any reasonably designed PDE lib should come with some environment variables to help user to setup, so the Makefile shouldn't to complicated. This applies to the above three libraries I've tested.

[yizhang-yiz]

Currently Stan doesn't have a sparse matrix system, and the mpi system is still young. Without these two common denominators I don't think we can get deeper than a black-box PDE feature design.

[betanalpha]

I agree that this is the ideal design — we do not want to go down the path of supporting PDE solvers, and this provides an elegant way to exposing Stan to custom solvers.

…

On Jul 6, 2018, at 13:14, Yi Zhang ***@***.***> wrote: Currently Stan doesn't have a sparse matrix system, and the mpi system is still young. Without these two common denominators I don't think we can get deeper than a black-box PDE feature design. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

[syclik]

I think solve_pde is good.