solve_pde function for PDE problems

yizhang-yiz commented 6 years ago

Summary:

Use solve_pde to solve PDEs.

Description:

The design involves cmdstan, stan, and math. See also https://github.com/stan-dev/stan/issues/2567. This issue is about Math repo only.

In order to solve PDE, similar to ODE solver in Stan, the PDE is provided through a functor. Unlike ODE solver, we cannot ship a PDE library with Stan Math. Alternatively, the user is asked to provide his own PDE solver. His workflow could be the following.

provide a functor that uses external PDE solver. The functor has 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 functor returns PDE solution(quantity of interest) when need_sens=0, and returns PDE solution together with its gradient with respect to theta(parameters) when need_sens=1.
call solve_pde, with the above functor as the first argument. Within the solve_pde function, the above functor calculates PDE solution and sensitivity, andvar variables are constructed based on these values.
~~provide Makefile for the external PDE solver~~. See separate issue #934 .

The user/PDE developer does not need to know utilize stan::math::var to write the above functor. He just needs to decide the behavior of the PDE solver: how the solution and the sensitivity are calculated. solve_pde should be able to take these values and construct proper var items, possibly using precomputed_gradient.

Additional Information:

Provide any additional information here.

Current Version:

v2.17.0

syclik commented 6 years ago

How is this issue supposed to be addressed?

Will there be a solve_pde function in the math library? What does that implementation need to do? Since it sounds like it needs to be linked against something to solve, how do we provide a solver ? Or is the user supposed to provide a solver?

Is there a reason we don’t handle the two different behaviors with two different functions?

I’m still trying to get the point of the issue before really digging into design.

On Tue, Jul 10, 2018 at 5:15 PM Yi Zhang notifications@github.com wrote:

Summary:

Use solve_pde to solve PDEs. Description:

The design involves cmdstan, stan, and math. See also stan-dev/stan#2567 https://github.com/stan-dev/stan/issues/2567. This issue is about Math repo only.

In order to solve PDE, similar to ODE solver in Stan, the PDE is provided through a functor. Unlike ODE solver, we cannot ship a PDE library with Stan Math. Alternatively, the user is asked to provide his own PDE solver. His workflow could be the following.

provide a functor that uses external PDE solver. The functor has the following signature

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

The functor returns PDE solution(quantity of interest) when need_sens=0, and returns PDE solution together with its gradient with respect to theta(parameters) when need_sens=1.

call solve_pde, with the above functor as the first argument. Within the solve_pde function, the above functor calculates PDE solution and sensitivity, andvar variables are constructed based on these values.

provide Makefile for the external PDE solver

The user/PDE developer does not need to know utilize stan::math::var to write the above functor. He just needs to decide the behavior of the PDE solver: how the solution and the sensitivity are calculated. solve_pde should be able to take these values and construct proper var items, possibly using precomputed_gradient. Additional Information:

Provide any additional information here. Current Version:

v2.17.0

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/931, or mute the thread https://github.com/notifications/unsubscribe-auth/AAZ_F6ZRvJ4MJNV8evk2GiUcA5mPOCSRks5uFRl4gaJpZM4VKIp8 .

yizhang-yiz commented 6 years ago

Will there be a solve_pde function in the math library?

Yes.

What does that implementation need to do?

For prim version, just use the PDE functor to solve PDE. For rev version, use the PDE functor to solve PDE and sensitivity, then construct var using precomputed_gradient.

Since it sounds like it needs to be linked against something to solve, how do we provide a solver ?

We don't.

Or is the user supposed to provide a solver?

Yes.

Is there a reason we don’t handle the two different behaviors with two different functions?

I'm not sure I understand the question. If by "two different behaviors" you mean solve /w and w/o sensitivity, then yes, we can certainly do two versions of functions, something like solve_pde and solve_pde_with_sensitivity. But that would require later on user know which one to use.

yizhang-yiz commented 6 years ago

Within Math library the solution would be similar to integrate_ode. The one proposed in #930 takes a PDE functor argument and use it to calculate PDE solutions. The difference between the solve_pde proposed in #930 and integrate_ode is that

The PDE functor passed to solve_pde is actually a PDE solution functor, while the ODE functor passed to integrate_ode simply calculates the ODE's right-hand-side. A PDE solution functor operates on parameters and returns PDE solutions(and sensitivities). This way in Math library we don't have to include a PDE solver.
In the above sense the proposed solve_pde function serves as Math libraries' API for external PDE solver developers. Users are asked to provide a functor that utilizes certain PDE solver, as well as makefile for that solver.

bob-carpenter commented 6 years ago

The transformed data and generated quantities blocks don't need gradients.

Ideally it's more efficient to have a separate implementation w/o gradients, but it's not necessary.

yizhang-yiz commented 6 years ago

There are two versions actually. One in prim one in rev. But in our case we are letting user routine to do real work, so there’s nothing stopping them to reuse same one.

On Wed, Jul 11, 2018 at 1:22 AM Bob Carpenter notifications@github.com wrote:

The transformed data and generated quantities blocks don't need gradients.

Ideally it's more efficient to have a separate implementation w/o gradients, but it's not necessary.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/931#issuecomment-404048539, or mute the thread https://github.com/notifications/unsubscribe-auth/AfkhsFsLZRX0peMpCS_Yh5_mXwb8Zi9Iks5uFYuBgaJpZM4VKIp8 .

-- Regards, Yi Zhang, Ph.D. Metrum Research Group LLC 2 Tunxis Road, Tariffville, CT 06081 yiz@metrumrg.com www.metrumrg.com

syclik commented 6 years ago

I'm still quite a bit confused. The user is supposed to provide a functor like this:

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);

And passes that into solve_pde()?

That makes sense, but then why do we need any changes to the build process? Is that necessary?

yizhang-yiz commented 6 years ago

Because we are asking user to provide his own solver. We don’t know how to compile/link the unknown library.

On Wed, Jul 11, 2018 at 8:36 AM Daniel Lee notifications@github.com wrote:

I'm still quite a bit confused. The user is supposed to provide a functor like this:

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

And passes that into solve_pde()?

That makes sense, but then why do we need any changes to the build process? Is that necessary?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/931#issuecomment-404153354, or mute the thread https://github.com/notifications/unsubscribe-auth/AfkhsGusTjBRdkq3O74KcvXNya4waRt_ks5uFfE0gaJpZM4VKIp8 .

-- Regards, Yi Zhang, Ph.D. Metrum Research Group LLC 2 Tunxis Road, Tariffville, CT 06081 yiz@metrumrg.com www.metrumrg.com

syclik commented 6 years ago

Shouldn't the user just put in make/local:

CXX += -I <path to source>
LDFLAGS += -l<library>

?

I don't understand why we need to define new variables. We did that for MPI and GPU, but I don't want to get in the habit of making our builds more complicated without there being a very good reason.

yizhang-yiz commented 6 years ago

In the most simple case, yes. But there are more complicated scenarios, especially for HPC users, they may want to adjust optimizations etc for tuning purpose, or may want to turn on/off certain functionality in his solver. Note that we're not asking user to replace Stan Math makefile and dive into our make process, we are asking him to provide only switches that are related to his solver. In the minimalist case, the makefile could only contain what you mentioned in the local. Of course another option is just let them use stuff local. The two approaches are the same after make's include does its job.

syclik commented 6 years ago

Can you provide a concrete example? I don't quite see what the complication is yet.

yizhang-yiz commented 6 years ago

Here's small of part of my petsc makefile rules. In general it goes through many decisions, such as if user needs to rebuild the library, or download another library, or if user is building fortran vs build c, if he's using another compiler, etc. Of course we can just ask user to include this file in local. Either way works for me.

# 3. Check if the shared libs are out of date
chkopts: chk_upgrade
    @for LIBNAME in ${SHLIBS}; do \
      library=${INSTALL_LIB_DIR}/$$LIBNAME.a; \
    sharedlibrary=${INSTALL_LIB_DIR}/$$LIBNAME.${SL_LINKER_SUFFIX}; \
    flag=""; \
    if [ -f $$library ]; then \
      if [ -f $$sharedlibrary ]; then \
        flag=`find ${INSTALL_LIB_DIR} -type f -name $$LIBNAME.a -newer $$sharedlibrary -print`; \
      fi; \
    fi; \
    if [ "$$flag" != "" ]; then \
      echo "Shared libs in ${INSTALL_LIB_DIR} are out of date, attempting to rebuild."; \
      if [ -w ${INSTALL_LIB_DIR} ]; then \
        ${OMAKE}  PETSC_ARCH=${PETSC_ARCH} shared; \
      else \
        printf ${PETSC_TEXT_HILIGHT}"*********************** ERROR ************************\n"; \
        echo "Unable to rebuild shared libraries; you do not have write permission."; \
        user=`ls -l ${INSTALL_LIB_DIR}/$$LIBNAME.${SL_LINKER_SUFFIX}  | tr -s ' ' | cut -d ' ' -f 3`; \
        echo "Libraries were built by user $$user; please contact him/her to have them rebuilt."; \
        printf "******************************************************"${PETSC_TEXT_NORMAL}"\n" ; \
        false; \
      fi; \
    fi; \
    done

#
# uses the cmake infrastructure to build/rebuild the libraries
ccmake:
    @echo "=========================================="
    +@cd ${PETSC_DIR}/${PETSC_ARCH} && MAKEFLAGS="-j$(MAKE_NP) -l$(NPMAX) $(MAKEFLAGS)" ${OMAKE} VERBOSE=${VERBOSE}
    @if [ "${BUILDSHAREDLIB}" = "yes" -a "${DSYMUTIL}" != "true" ]; then \
        echo "Running ${DSYMUTIL} on ${SHLIBS}";\
        for LIBNAME in ${SHLIBS}; do ${DSYMUTIL} ${INSTALL_LIB_DIR}/$$LIBNAME.${SL_LINKER_SUFFIX}; done; fi
    @echo "========================================="

cmake:
    +@MAKEFLAGS="-j$(MAKE_NP) -l$(NPMAX) $(MAKEFLAGS)" ${OMAKE} ccmake VERBOSE=1

gnumake:
    +@cd ${PETSC_DIR} && MAKEFLAGS="-j$(MAKE_NP) -l$(NPMAX) $(MAKEFLAGS)" ${OMAKE_PRINTDIR} -f gmakefile ${MAKE_PAR_OUT_FLG} V=${V}

# Does nothing; needed for some rules that require actions.
foo:

# Builds library
lib:  ${SOURCE}
    @${OMAKE}  PETSC_ARCH=${PETSC_ARCH} chk_petscdir
    @-if [ "${SPECIALLIB}" = "yes" ] ; then \
       ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  speciallib;  \
      else \
            if [ "${SOURCECU}" != "" ] ; then \
              ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  libcu; fi ; \
            if [ "${SOURCEC}" != "" ] ; then \
              ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  libc; fi ; \
            if [ "${SOURCECXX}" != "" ] ; then \
              ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  libcxx; fi ; \
            if [ "${SOURCEF}" != "" ] ; then \
        ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  libf; fi ; \
            if [ "${OBJS}" != " " ] ; then \
        ${RANLIB}  ${LIBNAME}; \
        ${RM} ${OBJS}; \
        fi;\
          fi
#
#  Does not work for some machines with .F fortran files.
#
# Builds library - fast versiong

libfast:  ${SOURCE}
    -@if [ "${SPECIALFASTLIB}" = "yes" ] ; then \
        ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  specialfastlib;  \
    else \
        ${OMAKE}  PETSC_ARCH=${PETSC_ARCH}  libfastcompile;  \
    fi

libfastcompile:
    -@if [ "${SOURCECU}" != "" ]; then \
                    ${PETSC_CUCOMPILE}; \
                ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSCU}; \
            ${RM} ${OBJSCU} ${OBJSCU:.o=.lo}; \
        fi; \
          if [ "${SOURCEC}" != "" ]; then \
                    ${PETSC_COMPILE}; \
                ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSC}; \
            ${RM} ${OBJSC} ${OBJSC:.o=.lo}; \
        fi; \
          if [ "${SOURCECXX}" != "" ]; then \
                    ${PETSC_CXXCOMPILE}; \
                ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSCXX}; \
            ${RM} ${OBJSCXX} ${OBJSCXX:.o=.lo}; \
        fi; \
    if [ "${SOURCEF}" != "" ]; then \
                    ${PETSC_FCOMPILE}; \
                ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSF}; \
            ${RM} ${OBJSF}; \
        fi

# Build f90 .mod from .F source, and add .o to the corresponding library
buildmod:
    -@${OMAKE} clean-legacy
    @${OMAKE}  libf
    @${OMAKE}  modcopy
    @${OMAKE}  clean-legacy

buildmodfast:
    -@${OMAKE} clean-legacy
    @${OMAKE}  libfastcompile
    @${OMAKE}  modcopy
    @${OMAKE}  clean-legacy

# copy modules to the include dir
modcopy:
    @${CP} -f *.mod ${PETSC_DIR}/${PETSC_ARCH}/include

syclik commented 6 years ago

Sorry... I have no idea what this is! It looks complicated, but beyond that, it doesn't have a complete example. Can you guide me through what parts of this are necessary and how adding more complexity to the Stan Math makefiles will make this easier?

On Wed, Jul 11, 2018 at 9:52 AM Yi Zhang notifications@github.com wrote:

Here's small of part of my petsc makefile rules. In general it goes through many decisions, such as if user needs to rebuild the library, or download another library, or if user is building fortran vs build c, if he's using another compiler, etc. Of course we can just ask user to include this file in local. Either way works for me.

3. Check if the shared libs are out of datechkopts: chk_upgrade

@for LIBNAME in ${SHLIBS}; do \ library=${INSTALL_LIB_DIR}/$$LIBNAME.a; \ sharedlibrary=${INSTALL_LIB_DIR}/$$LIBNAME.${SL_LINKER_SUFFIX}; \ flag=""; \ if [ -f $$library ]; then \ if [ -f $$sharedlibrary ]; then \ flag=find ${INSTALL_LIB_DIR} -type f -name $$LIBNAME.a -newer $$sharedlibrary -print; \ fi; \ fi; \ if [ "$$flag" != "" ]; then \ echo "Shared libs in ${INSTALL_LIB_DIR} are out of date, attempting to rebuild."; \ if [ -w ${INSTALL_LIB_DIR} ]; then \ ${OMAKE} PETSC_ARCH=${PETSC_ARCH} shared; \ else \ printf ${PETSC_TEXT_HILIGHT}"*** ERROR ****\n"; \ echo "Unable to rebuild shared libraries; you do not have write permission."; \ user=ls -l ${INSTALL_LIB_DIR}/$$LIBNAME.${SL_LINKER_SUFFIX} | tr -s ' ' | cut -d ' ' -f 3; \ echo "Libraries were built by user $$user; please contact him/her to have them rebuilt."; \ printf "**"${PETSC_TEXT_NORMAL}"\n" ; \ false; \ fi; \ fi; \ done

uses the cmake infrastructure to build/rebuild the librariesccmake:

@echo "==========================================" +@cd ${PETSC_DIR}/${PETSC_ARCH} && MAKEFLAGS="-j$(MAKE_NP) -l$(NPMAX) $(MAKEFLAGS)" ${OMAKE} VERBOSE=${VERBOSE} @if [ "${BUILDSHAREDLIB}" = "yes" -a "${DSYMUTIL}" != "true" ]; then \ echo "Running ${DSYMUTIL} on ${SHLIBS}";\ for LIBNAME in ${SHLIBS}; do ${DSYMUTIL} ${INSTALL_LIB_DIR}/$$LIBNAME.${SL_LINKER_SUFFIX}; done; fi @echo "=========================================" cmake: +@MAKEFLAGS="-j$(MAKE_NP) -l$(NPMAX) $(MAKEFLAGS)" ${OMAKE} ccmake VERBOSE=1 gnumake: +@cd ${PETSC_DIR} && MAKEFLAGS="-j$(MAKE_NP) -l$(NPMAX) $(MAKEFLAGS)" ${OMAKE_PRINTDIR} -f gmakefile ${MAKE_PAR_OUT_FLG} V=${V}

Does nothing; needed for some rules that require actions.foo:

Builds librarylib: ${SOURCE}

@${OMAKE} PETSC_ARCH=${PETSC_ARCH} chk_petscdir @-if [ "${SPECIALLIB}" = "yes" ] ; then \ ${OMAKE} PETSC_ARCH=${PETSC_ARCH} speciallib; \ else \ if [ "${SOURCECU}" != "" ] ; then \ ${OMAKE} PETSC_ARCH=${PETSC_ARCH} libcu; fi ; \ if [ "${SOURCEC}" != "" ] ; then \ ${OMAKE} PETSC_ARCH=${PETSC_ARCH} libc; fi ; \ if [ "${SOURCECXX}" != "" ] ; then \ ${OMAKE} PETSC_ARCH=${PETSC_ARCH} libcxx; fi ; \ if [ "${SOURCEF}" != "" ] ; then \ ${OMAKE} PETSC_ARCH=${PETSC_ARCH} libf; fi ; \ if [ "${OBJS}" != " " ] ; then \ ${RANLIB} ${LIBNAME}; \ ${RM} ${OBJS}; \ fi;\ fi## Does not work for some machines with .F fortran files.## Builds library - fast versiong libfast: ${SOURCE} -@if [ "${SPECIALFASTLIB}" = "yes" ] ; then \ ${OMAKE} PETSC_ARCH=${PETSC_ARCH} specialfastlib; \ else \ ${OMAKE} PETSC_ARCH=${PETSC_ARCH} libfastcompile; \ fi libfastcompile: -@if [ "${SOURCECU}" != "" ]; then \ ${PETSC_CUCOMPILE}; \ ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSCU}; \ ${RM} ${OBJSCU} ${OBJSCU:.o=.lo}; \ fi; \ if [ "${SOURCEC}" != "" ]; then \ ${PETSC_COMPILE}; \ ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSC}; \ ${RM} ${OBJSC} ${OBJSC:.o=.lo}; \ fi; \ if [ "${SOURCECXX}" != "" ]; then \ ${PETSC_CXXCOMPILE}; \ ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSCXX}; \ ${RM} ${OBJSCXX} ${OBJSCXX:.o=.lo}; \ fi; \ if [ "${SOURCEF}" != "" ]; then \ ${PETSC_FCOMPILE}; \ ${AR} ${FAST_AR_FLAGS} ${LIBNAME} ${OBJSF}; \ ${RM} ${OBJSF}; \ fi

Build f90 .mod from .F source, and add .o to the corresponding librarybuildmod:

-@${OMAKE} clean-legacy @${OMAKE} libf @${OMAKE} modcopy @${OMAKE} clean-legacy buildmodfast: -@${OMAKE} clean-legacy @${OMAKE} libfastcompile @${OMAKE} modcopy @${OMAKE} clean-legacy

copy modules to the include dirmodcopy:

@${CP} -f *.mod ${PETSC_DIR}/${PETSC_ARCH}/include

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/931#issuecomment-404176890, or mute the thread https://github.com/notifications/unsubscribe-auth/AAZ_FzyVxrxlSp9xEuWA7bTiWuOA1TZYks5uFgLMgaJpZM4VKIp8 .

yizhang-yiz commented 6 years ago

I think the bottom line is that user are looking at different optimization options, and different architectures that he may want to play with his PDE solver. Also his source code might be C or C++ or fortran or F90, all asking for special treatment.

It does add more complexity either way, IMO. Because one can either bloat up local or include the external makefile. But I think it's conceptually easier to understand to use a designated makefile for external library.

syclik commented 6 years ago

Sorry -- I feel like you're speaking a different language. Can you define these terms?

optimization options
different architectures
play with his PDE solver
special treatment

It really sounds like all of what you've described is outside of what the Math library needs to worry about. Meaning, the user should really have some way of building whatever source code in C or C++ fortran or F90. Once that's there, all the user really needs is a way to link? Or am I misunderstanding something about the process you're describing?

On Wed, Jul 11, 2018 at 10:08 AM Yi Zhang notifications@github.com wrote:

I think the bottom line is that user are looking at different optimization options, and different architectures that he may want to play with his PDE solver. Also his source code might be C or C++ or fortran or F90, all asking for special treatment.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/931#issuecomment-404180991, or mute the thread https://github.com/notifications/unsubscribe-auth/AAZ_Fx3Inut1z8dqimtOvTQuBT30aOO0ks5uFgWEgaJpZM4VKIp8 .

yizhang-yiz commented 6 years ago

Sorry for being a bad interpreter.

https://github.com/yizhang-cae/cmdstan/tree/forward_pde/examples/laplace_pde Here's the example I mentioned in https://github.com/stan-dev/stan/issues/2567. The makefile for libmesh used in this example requires the follow minimum setting:

CXX = $(libmesh_CXX)
CC = $(libmesh_CXX)
CPPFLAGS += $(libmesh_CPPFLAGS)
CXXFLAGS += $(libmesh_CXXFLAGS)
CXXFLAGS += -isystem $(LIBMESH_DIR)/include
CXXFLAGS += $(libmesh_INCLUDE)
CXXFLAGS :=$(filter-out -std=gnu++11,$(CXXFLAGS))
CXXFLAGS :=$(filter-out -Wunused,$(CXXFLAGS))
CXXFLAGS :=$(filter-out -Wunused-parameter,$(CXXFLAGS))
CXXFLAGS :=$(filter-out -Qunused-arguments,$(CXXFLAGS))
CXXFLAGS +=-DLIBMESH_HAVE_CXX14_MAKE_UNIQUE

LDFLAGS += $(libmesh_LIBS)
LDFLAGS += $(libmesh_LDFLAGS)
LDFLAGS += $(EXTERNAL_FLAGS)

# Useful rules.
dust:
    @echo "Deleting old output and runtime files"
    @rm -f out*.m job_output.txt output.txt* *.gmv.* *.plt.* out*.xdr* out*.xda* PI* *.o *.libs

gmv:
    @$(MAKE) -C $(LIBMESH_DIR)/roy/meshplot/ meshplot-$(METHOD)
    @for file in out.mesh.*; do ${LIBMESH_RUN} $(LIBMESH_DIR)/roy/meshplot/meshplot-$(METHOD) $$file out.soln.$${file##out.mesh.} out.gmv.$${file:9:4}; done

What I meant is

optimization options

Questions like use -O3 vs -O2? Link withmetis(a partitioner) or scotch(another partitioner)? etc.

different architectures

Questions like Build with MPICH or OPENMPI? linear solver being SUPERLU or MUMPS? These are related to how a solver library is built. Some libraries (like petsc) allow the user to build different versions with different setups, aka "architectures".

play with his PDE solver

Large problems often requires some tuning of the aforementioned building parameters before sent out to the cluster for days-long run. What I meant by "play" is this tuning process, which involves change switches in the building process.

special treatment

One example is that if the solver is written in fortran, user needs to supply mpif or gfortran compiler/linker information.

bbbales2 commented 6 years ago

There's nothing here that's PDE specific, so is the goal to interface with external software in general?

It seems like there are three big components of all these pull reqs:

The high level interface (that would appear in .stan files)
The mechanism for working with the autodiff library without having to worry too much about the details of vars
Whatever is needed to link in external libraries

3 seems like a support nightmare to me, and 1 seems easy to get around (if someone can build petsc, I'll hazard they can figure out how to get the Rstan/cmdstan external C++ stuff working, even if those are hacky solutions). Seems to me like if someone has some big fancy weird external solver they need to incorporate in Stan they could maintain it as a fork? That's what I do for my regular job -- not that my development patterns are anything to emulate...

I'm working with Bob now on 2, for other reasons, but it's the same goal: https://github.com/stan-dev/math/pull/924 . I'll be full time on this until September. We think it's a good idea too -- and someone else coming up with it independently seems like a good sign :P. Fingers crossed that should progress quickly.

Maybe this should be an item for the engineering part of the Thursday meeting?

yizhang-yiz commented 6 years ago

As I mentioned in https://github.com/stan-dev/stan/issues/2567, this is indeed true to any external library. for 3 I don't think asking PDE developer to maintain a fork is a better option(I'm doing that with Torsten, and I'm desperately looking for alternatives) then providing a good API.

yizhang-yiz commented 6 years ago

Add @bob-carpenter for Thursday agenda: mechanism to include external library.

yizhang-yiz commented 6 years ago

A related issue is https://github.com/stan-dev/stan/issues/2572: to have a way to support user-supplied higher-order function, like regular function through user_header.hpp. Though this is more a Stan issue than Math issue.

syclik commented 6 years ago

As I mentioned in stan-dev/stan#2567, this is indeed true to any external library. for 3 I don't think asking PDE developer to maintain a fork is a better option(I'm doing that with Torsten, and I'm desperately looking for alternatives) then providing a good API.

Re: Torsten. We should put this on another thread. There was an active decision to keep that separate due to lack of resources for design, documentation, testing, and maintenance. We can revisit that if you think that's appropriate.

yizhang-yiz commented 6 years ago

I didn't mean to bring Torsten into Stan or involve Torsten in the current design, what I meant is I'd prefer a way to use Stan with external libraries through API(which is the topic of this thread) to include Stan as a fork. This is to respond to @bbbales2 's comment. Sorry for the confusion.

syclik commented 6 years ago

@yizhang-cae, can we split the issue into two? I think there are two very different things going on:

solve_pde(). That looks fine and easy to review and merge. There's nothing different about that.
Building external tools and linking them. That's a completely separate, orthogonal discussion.

For 1, we can move forward. For 2, it would really help to have a list of things you want (a spec would be great, but even a list of wants and needs would be sufficient), then we can discuss whether that's a good thing and how we can make that happen.

yizhang-yiz commented 6 years ago

I've edited the issue and added #934 . Let's focus on solve_pde in the current issue.

yizhang-yiz commented 6 years ago

@syclik , for the current issue, please allow me to pitch the name forward_pde again. I understand solve_pde is a more popular name, but it's also a misleading name in the context of statistical inference: solve_pde gives an impression that the return would be the solution of the PDE, but that's not the intended use case. In PDE inference the data are often not collected from PDE solution directly but some functionals based on the PDE solution, usually involving spatial and/or time integrals. This is the "forward" process from PDE parameters to the observable functionals. It involves solving the PDE, but the numerical solution of PDE is often hidden behind.

bob-carpenter commented 6 years ago

Yes, please---forward_pde sounds like a much better name. Let's use the standard names from applied math wherever possible so that we don't confuse people who know what they're doing! It's win-win as it's a learning experience for the rest of us. (This is related to Ben's change of the function name optimize() to maximize() to not confuse the professionals for whom optimization will often involve minimization.)

On Jul 11, 2018, at 11:40 AM, Yi Zhang notifications@github.com wrote:

For the current issue, please allow me to pitch the name forward_pde again. I understand solve_pde is a more popular name, but it's also a misleading name: solve_pde gives an impression that the return would be the solution of the PDE, but that's not the intended use case. In PDE inference the data are often not collected from PDE solution directly but some functionals based on the PDE solution, usually involving spatial and/or time integrals. This is the "forward" process from PDE parameters to the observable functionals. It involves solving the PDE, but the numerical solution of PDE is often hidden behind.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

syclik commented 6 years ago

Btw, what you've described as "optimization options" are already available. In make/local, just set:

O = 3

or

O = 2

For the other ones, I'm not sure? These are all different use cases. Knowing what you're trying to do would really help. (It's still not clear enough for me to understand exactly what you're trying to do to help design.)

syclik commented 6 years ago

And yes, forward_pde sounds much more appropriate given your last comment.

syclik commented 6 years ago

That said Matlab calls it solvepde: https://www.mathworks.com/help/pde/pde-solvers.html.

What about simulate_pde? TensorFlow's doc talks about simulating from a PDE: https://www.tensorflow.org/tutorials/non-ml/pdes

So does iPython Cookbook: https://ipython-books.github.io/124-simulating-a-partial-differential-equation-reaction-diffusion-systems-and-turing-patterns/

And just to add more confusion, FEniCS calls it solve (they do PDE only): https://fenicsproject.org/

yizhang-yiz commented 6 years ago

That's exactly what I'm talking about.

Matlab:

This solver returns a StationaryResults or TimeDependentResults object whose properties contain the solution and its gradient at the mesh nodes.

FEniCS

Compute solution w = Function(W) solve(a == L, w, [bc1, bc0])

TensorFlow

Discretized PDE update rules U = U + eps * Ut Ut = Ut + eps (laplace(U) - damping Ut)

All the above examples solve for the numerical solution. None of the above example ask for QoI(quantity of interest).

BTW, the adjoint QoI solver is treated separately in FEniCS by DOLFIN, so the example above is a misquote.

syclik commented 6 years ago

I missed the point about the QoI. In the words / notation from Matlab, FEniCS, and TensorFlow, what's different about what you're proposing the Math library handles?

syclik commented 6 years ago

btw, would this go a lot quicker if we jumped on a video call?

yizhang-yiz commented 6 years ago

Sure. Hangout? yiz@metrumrg.com

yizhang-yiz commented 6 years ago

In the words / notation from Matlab, FEniCS, and TensorFlow, what's different about what you're proposing the Math library handles?

Don't remember how Matlab & tensorflow treat it, but for FEniCS the adjoint sensitivity is treated in a sublibrary http://www.dolfin-adjoint.org/. What we are asking user to provide in the PDE functor is the combination of the PDE solution like quoted above and the calculation of QoI like in http://www.dolfin-adjoint.org/. Of course since user is responsible for the behavior of the functor, nothing stops him from returning a full-blown PDE solution and stuff the cache. That's why hopefully the name forward would mean something to someone who's familiar with PDE inverse problems.

syclik commented 6 years ago

Would it be clearer if we didn’t call if “pde” but something like “pde_qoi”? Since the functor is not a pde, I don’t think we should call it a pde?

On Wed, Jul 11, 2018 at 1:13 PM Yi Zhang notifications@github.com wrote:

In the words / notation from Matlab, FEniCS, and TensorFlow, what's different about what you're proposing the Math library handles?

Don't remember how Matlab & tensorflow treat it, but for FEniCS the adjoint sensitivity is treated in a sublibrary http://www.dolfin-adjoint.org/. What we are asking user to provide in the PDE functor is the combination of the PDE solution like quoted above and the calculation of QoI like in http://www.dolfin-adjoint.org/. Of course since user is responsible for the behavior of the functor, nothing stops him from returning a full-blown PDE solution and stuff the cache. That's why hopefully the name forward would mean something to someone who's familiar with PDE inverse problems.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://github.com/stan-dev/math/issues/931#issuecomment-404244698, or mute the thread https://github.com/notifications/unsubscribe-auth/AAZ_F8_yVQuYZgkfrBEVdFsiiPlePj0Pks5uFjH8gaJpZM4VKIp8 .

yizhang-yiz commented 6 years ago

Would it be clearer if we didn’t call if “pde” but something like “pde_qoi”? Since the functor is not a pde, I don’t think we should call it a pde?

@syclik , good point. Changes made at https://github.com/stan-dev/math/pull/930#issuecomment-404609699

yizhang-yiz commented 5 years ago

@syclik I'd like to have it merged into Torsten first so I'm closing it for the moment.

stan-dev / math