Feature/issue #767 DAE solver #768

yizhang-yiz · 2018-02-15T18:16:09Z

Submission Checklist

Run unit tests: ./runTests.py test/unit
Run cpplint: make cpplint
Declare copyright holder and open-source license: see below

Summary:

This PR addresses #767 . It implements interface to IDAS library for DAE solutions.

A DAE solver solves equations in the form
F(yp, yy, t) = 0
where F: R^n -> R^n contains both differential and algebraic components, yy is the unknown variable vector, and yp=yy' is the derivative of yy. In a trivia DAE the differential and algebraic equation components are decoupled. When they are not, an involved process needs to be formulated to convert the system into canonical forms.

Intended Effect:

The proposed interface is as follows

  template<typename F, typename Tyy, typename Typ, typename Tpar>
  std::vector<std::vector<typename stan::return_type<Tyy, Typ, Tpar>::type> >
  integrate_dae(const F& f,
                const std::vector<int>& eq_id,
                const std::vector<Tyy>& yy0,
                const std::vector<Typ>& yp0,
                double t0,
                const std::vector<double>& ts,
                const std::vector<Tpar>& theta,
                const std::vector<double>& x_r,
                const std::vector<int>& x_i,
                const double rtol = 1e-8,
                const double atol = 1e-10,
                const int64_t max_num_steps = 1e6,
                std::ostream* msgs = 0) {

f: the function for DAE residual
eq_id: integer vector to identify differential and algebraic components in the system. 1 for differential, 0 for algebraic.
yy0: initial condition for primary unknowns
yp0: initial condition for derivative of unknowns. yy0 and yp0 together must meet F(yp0, yy0,t0)=0
t0: initial time
ts: vector for observation time
theta: parameters

When TYY, TYP and TPAR are not all data, the DAE system is augmented with forward sensitivity variables s and s' for sensitivity analysis.

A demo for prey-predator-harvest dynamical system.

How to Verify:

Unit tests

Side Effects:

N/A

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company):
Metrum Research Group, LLC

By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:

Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

bob-carpenter · 2018-02-15T19:08:47Z

@yizhang-cae Would you mind updating this to include copyright for Metrum? Thanks.

yizhang-yiz · 2018-02-15T19:10:21Z

Put the line in a wrong spot. Corrected.

yizhang-yiz · 2018-02-26T16:19:16Z

I wonder why Jenkins fails to resolve headers. All my gcc and clang builds are fine. Maybe @seantalts can take a look? Thanks.

seantalts · 2018-02-26T16:24:15Z

I'm not sure, this looks like a pretty legit error to me (i.e. not a problem with Jenkins). Did you add a new library and not add it to the include path properly or something?

seantalts · 2018-02-26T16:28:34Z

I think you're just missing the appropriate includes, like

#include <stan/math/prim/scal/err/check_greater_or_equal.hpp>

yizhang-yiz · 2018-02-26T16:29:44Z

Mmm, I wonder why it didn't break my local build.

seantalts · 2018-02-26T16:34:08Z

Are you running the same command Jenkins is, i.e. make test-headers? It's a special invocation that tests that you are including the right things in your files, etc.

yizhang-yiz · 2018-02-26T16:52:48Z

No...didn't know that.

yizhang-yiz · 2018-02-28T14:38:33Z

Upstream tests fail with cmdstan. So I guess I need update cmdstan's makefile at the same time. @seantalts, what's the workflow here?

syclik · 2018-02-28T14:47:59Z

The steps are:

create an issue on CmdStan, RStan, and PyStan for adding the includes. This affects all 3.
create a pull request on CmdStan to update the includes.
We'll need to review both together. We'll want to merge math, then CmdStan shortly after.

Just FYI, we'll have to discuss this with the RStan and PyStan maintainers because doing something like this isn't as easy to extend as in CmdStan.

yizhang-yiz · 2018-02-28T14:52:46Z

It figures. In that case we may want to do the sundials merge first(#779).

syclik · 2018-02-28T14:57:36Z

That makes sense.

…

On Wed, Feb 28, 2018 at 9:53 AM, Yi Zhang ***@***.***> wrote: It figures. In that case we may want to do the sundials merge first(#779 <#779>). — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#768 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAZ_F4QqiwGIPN-6dpJ8kXoyFKdeF5Q0ks5tZWg_gaJpZM4SHR_8> .

wds15 · 2018-04-29T15:08:13Z

One question: In ODE problems there is usually a stiff and a non-stiff variant. Do we have the same choices with the DAE solver? Or are these problems always stiff?

yizhang-yiz · 2018-04-29T21:44:13Z

In theory, there could be problems that are nonstill. But those problems are often convertible to ODEs, while DAEs in real applications are addressed from reducing ODEs systems using singular perturbation. This process usually results in stiff ODEs numerical problems. That's why most DAEs are solved as stiff ODEs numerically. Regards, Yi Zhang, Ph.D. Metrum Research Group LLC 2 Tunxis Road, Tariffville, CT 06081 yiz@metrumrg.com www.metrumrg.com

…

On Sun, Apr 29, 2018 at 11:08 AM, wds15 ***@***.***> wrote: One question: In ODE problems there is usually a stiff and a non-stiff variant. Do we have the same choices with the DAE solver? Or are these problems always stiff? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#768 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AfkhsEoOPXY--HF84wvrJyB4zV3Eitunks5ttddegaJpZM4SHR_8> .

syclik · 2018-05-24T03:32:16Z

@yizhang-cae, could you work on the merge conflicts? If it's not ready to be reviewed, please close the pull request and reopen when ready.

yizhang-yiz · 2018-05-24T14:56:07Z

It has been waiting till sundials-combine PR gets merged. I’ll update it.

On Wed, May 23, 2018 at 11:32 PM Daniel Lee ***@***.***> wrote: @yizhang-cae <https://github.com/yizhang-cae>, could you work on the merge conflicts? If it's not ready to be reviewed, please close the pull request and reopen when ready. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#768 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AfkhsM73D_fvudsKgn_zAzIn1ymSFrxBks5t1inBgaJpZM4SHR_8> .

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

syclik · 2018-05-24T15:07:15Z

Thanks! Sorry -- I didn't know that's what it was waiting for.

…

On Thu, May 24, 2018 at 10:57 AM Yi Zhang ***@***.***> wrote: It has been waiting till sundials-combine PR gets merged. I’ll update it. On Wed, May 23, 2018 at 11:32 PM Daniel Lee ***@***.***> wrote: > @yizhang-cae <https://github.com/yizhang-cae>, could you work on the > merge conflicts? If it's not ready to be reviewed, please close the pull > request and reopen when ready. > > — > You are receiving this because you were mentioned. > Reply to this email directly, view it on GitHub > <#768 (comment)>, or mute > the thread > < https://github.com/notifications/unsubscribe-auth/AfkhsM73D_fvudsKgn_zAzIn1ymSFrxBks5t1inBgaJpZM4SHR_8 > > . > -- Regards, Yi Zhang, Ph.D. Metrum Research Group LLC 2 Tunxis Road, Tariffville, CT 06081 ***@***.*** www.metrumrg.com — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#768 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAZ_F38r7Mvvc4sGUO1dsXpvARG8DoFRks5t1soTgaJpZM4SHR_8> .

yizhang-yiz · 2018-05-31T14:24:46Z

@seantalts , pipeline tests show "unstable". What does it mean?

seantalts · 2018-05-31T14:28:43Z

A test failed: http://d1m1s1b1.stat.columbia.edu:8080/blue/organizations/jenkins/Math%20Pipeline/detail/PR-768/19/tests

The parsing of the test failure also failed because the xml file wasn't valid, which makes me suspect multiple mpi processes were executing this test writing to the same file or something. @wds15 any ideas about this one?

yizhang-yiz · 2018-05-31T14:51:49Z

@seantalts, I see only "fail-to-read" of xml. The MPI tests actually are passed.

seantalts · 2018-05-31T14:54:23Z

Yep. The XML hasn't been created properly even though the tests pass. I think we need @wds15 to look at it - he did something special to enable some mpi tests to be run multicore and still write valid XMLs (...most of the time).

wds15 · 2018-05-31T15:10:52Z

Ok, I will take a look these days.

wds15 · 2018-06-03T12:51:08Z

this is a bit of mystery as to why it failed... could you please try just again? So restart the testing? Thanks.

seantalts · 2018-06-03T16:39:08Z

Restarting! Would also like to keep an eye on these tests to see if they're flaky and fix them if so.

yizhang-yiz · 2018-06-04T17:11:13Z

Thanks @seantalts.
@syclik , I've rebased this one and it passed tests.

yizhang-yiz · 2018-06-12T17:09:57Z

Hi @syclik @bob-carpenter , who would be a proper reviewer for this one?

…error handling

…stable/2017-11-14)

…zhang-cae/math into feature/issue-#767-DAE-solver

…xp1~20180509124008.99 (branches/release_50)

yizhang-yiz · 2018-07-28T01:52:11Z

@bbbales2 , I've cleaned up according to the reviews. In particular, new DAE model has been added to unit tests. Thanks.

bbbales2 · 2018-08-01T07:33:20Z

Man every time I come back to these pull requests I am surprised by how long it took me to get around to it... 5 days haha. I gotta work on my turnaround.

Why keep things templated on the types of yy0 and yp0? If we don't have any plans of addressing the consistency thing, that should be removed and these should be doubles always. It's currently unusable for non-trivial reasons. For sustainability, stuff like this needs to stay out of develop.

Also delete code that's not being used (https://github.com/stan-dev/math/pull/768/files#diff-6739115919f8681b87e9bb7e13a5cf16R49). If there's something special about the adjoint DAE systems that needs documented, that should go in an Implement Adjoint DAE issue.

Speaking of adjoint sensitivity stuff, are you interested in doing it for the regular ODEs? I know my buddy @pourzanj has been wanting this stuff (he works with @syclik at Generable). I'd be curious to see it in action as well. A friend of mine is working on a 1D diffusion PDE via method of lines in Stan. The number of states is really big, so the forward sensitivity problem gets big kinda quickly too and I think adjoint sensitivity might help. You're the most qualified person here to do this, though I figure you work more on what Metrum asks of you rather than what me or Arya do :P. I figured I'd ask though.

I'm gonna go ahead and check the memory allocation and docs more closely. I like the tests. This is close to being through.

bbbales2

Here's the latest set of changes. I just traced through a call of integrate_dae and tried to pay close attention to the memory. I read docs when I remembered to. It seems fine, and the code is readable.

bbbales2 · 2018-08-01T08:07:22Z

stan/math/rev/mat/functor/integrate_dae.hpp

+    const double rtol, const double atol,
+    const int64_t max_num_steps = idas_integrator::IDAS_MAX_STEPS,
+    std::ostream* msgs = nullptr) {
+  const std::vector<int> dummy_eq_id(yy0.size(), 0);


This says every equation is an algebraic equation. Shouldn't this be an argument to the function?

It doesn't matter. eq_id only works when yy0 is var. But I've added comments to clarify.

bbbales2 · 2018-08-01T08:07:32Z

stan/math/rev/mat/functor/integrate_dae.hpp

+ * @param[in] theta parameter
+ * @param[in] x_r real data
+ * @param[in] x_i int data
+ * @param[in] rtol relative tolerance passed to IDAS, requred <10^3


<10^-3, though I guess <10^3 is true too haha

bbbales2 · 2018-08-01T08:09:45Z

stan/math/rev/mat/functor/integrate_dae.hpp

+ * @param[in] f functor for the base ordinary differential equation
+ * @param[in] yy0 initial state
+ * @param[in] yp0 initial derivative state
+ * @param[in] t0 initial time.


Rogue period

bbbales2 · 2018-08-01T08:10:11Z

stan/math/rev/mat/functor/integrate_dae.hpp

+ * given the specified consistent initial state yy0 and yp0.
+ *
+ * @tparam DAE type of DAE system
+ * @tparam Tpar type of parameter theta


scalar type of parameters in theta

bbbales2 · 2018-08-01T08:10:31Z

stan/math/rev/mat/functor/integrate_dae.hpp

+ * @param[in] t0 initial time.
+ * @param[in] ts times of the desired solutions, in strictly
+ * increasing order, all greater than the initial time
+ * @param[in] theta parameter


bbbales2 · 2018-08-01T09:49:55Z

stan/math/rev/mat/functor/idas_integrator.hpp

+
+/**
+ * Solve Dae system with forward sensitivty, return a
+ * vector of var with recomputed gradient as sensitivity value


precomputed

bbbales2 · 2018-08-01T09:54:16Z

stan/math/rev/mat/functor/idas_integrator.hpp

+#include <stan/math/rev/scal/meta/is_var.hpp>
+#include <stan/math/prim/scal/meta/return_type.hpp>
+#include <stan/math/rev/mat/functor/idas_forward_system.hpp>
+// #include <stan/math/rev/arr/fun/decouple_ode_states.hpp>


Remove this. Also value_of, is_var, return_type, not used. Is algorithm used?

bbbales2 · 2018-08-01T10:04:43Z

test/unit/math/rev/mat/functor/idas_system_test.cpp

+        f, eq_id, yy0, yp0, theta, x_r, x_i, msgs};
+  };
+
+  ASSERT_NO_THROW(build_double(eq_id, yy0, yp0, theta, x_r));


Add tests to verify things throw when length(yy0), length(yp0), and length(eq_id) aren't all equal

bbbales2 · 2018-08-01T10:07:17Z

test/unit/math/rev/mat/functor/idas_integrator_test.cpp

+                   "times");
+  bad_t0 = 0;
+  bad_ts[0] = -1;
+  EXPECT_THROW_MSG(solver.integrate(dae, bad_t0, bad_ts), std::domain_error,


Add a test that throws error if ts isn't strictly ordered. So like {0.0, -1.0, 2.0} and {0.0, 0.0, 2.0} should fail.

bbbales2 · 2018-08-01T10:12:50Z

stan/math/rev/mat/functor/integrate_dae.hpp

+namespace math {
+/**
+ * Return the solutions for a semi-explicit DAE system with residual
+ * specified by functor F,


I think the functor F should be documented at least somewhere, and I think this is the right place. This is what we did for integrate_1d

math/stan/math/rev/arr/functor/integrate_1d.hpp

Line 65 in b5a1873

* f should be compatible with reverse mode autodiff and have the signature:

(rev) and

math/stan/math/prim/arr/functor/integrate_1d.hpp

Line 152 in b5a1873

* The signature for f should be:

(prim)

…xp1~20180509124008.99 (branches/release_50)

bbbales2 · 2018-08-07T09:21:02Z

@yizhang-cae Looks good to me. Thanks Yi! Lemme know when the stan-dev/stan pull is ready. I'll tell my DAE solvin' buddy they're coming soon.

@syclik This touches one line of one makefile in what looks like a very innocuous way. Is the hold still in place for makefile-touching pulls? merge this if that change looks okay to you.

syclik · 2018-08-07T11:14:46Z

Yes. Please hold. The builds have to be fixed.

…

On Tue, Aug 7, 2018 at 5:21 AM Ben Bales ***@***.***> wrote: @yizhang-cae <https://github.com/yizhang-cae> Looks good to me. Thanks Yi! Lemme know when the stan-dev/stan pull is ready. I'll tell my DAE solvin' buddy they're coming soon. @syclik <https://github.com/syclik> This touches one line of one makefile in what looks like a very innocuous way. Is the hold still in place for makefile-touching pulls? merge this if that change looks okay to you. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#768 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAZ_F5vhgzfBiXNmYUw386gWpWVHFtCHks5uOVv_gaJpZM4SHR_8> .

syclik · 2018-08-07T11:20:37Z

I just saw the change to the makefile. If someone makes that change on the branch with the make updates, we can merge this.

syclik · 2018-08-07T13:52:09Z

@yizhang-cae, please update this pull request: #956; branch: bugfix/issue-954-make

yizhang-yiz · 2018-08-07T13:52:09Z

@bbbales2 , thank you.

As to adjoint sensitivity, it's on my list but not priority. It's best to do ODE & DAE's adjoint together and I'd like to see DAE go in first. It doesn't hurt to create an issue specifying the problem though.

yizhang-yiz · 2018-08-07T13:54:05Z

@syclik , one sec.

bbbales2 · 2018-08-07T19:47:25Z

@yizhang-cae I was talking to @bob-carpenter about adjoint vs. forward sensitivity stuff today. He said he was gonna make it a meeting agenda item for Thursday cause he wanted to talk about it some more. So definitely show up this Thursday if you can make it, cause your input will be good.

ryan-richt · 2021-12-08T21:27:47Z

Hello All!
My team and I would be super interested in seeing this go in. Could we kindly inquire as to it's current status?

yizhang-yiz · 2022-01-04T06:58:57Z

@ryan-richt I'm close to a PR. I wonder what the DAE dimension your team is working on, perhaps I can use it for unit test or performance test.

ryan-richt · 2022-01-04T14:15:51Z

Super exciting @yizhang-yiz ! Sadly we can’t distribute the code of this particular model, but when you have a release candidate we’re happy to run it against your branch and report back metrics and the dimensionality / other meta data if that helps. (It’s pretty big/slow, like seconds per evaluation, so probably more on the perf than unit test side) Thank you again!

yizhang-yiz changed the title ~~Feature/issue #767 dae solver~~ Feature/issue #767 DAE solver Feb 15, 2018

yizhang-yiz mentioned this pull request Feb 24, 2018

Update cvodes to 3.1.0 #744

Merged

3 tasks

yizhang-yiz force-pushed the feature/issue-#767-DAE-solver branch from af6a206 to d59d73b Compare February 26, 2018 17:55

yizhang-yiz force-pushed the feature/issue-#767-DAE-solver branch from 5da5e21 to 92ef70f Compare May 30, 2018 01:47

Yi Zhang and others added 9 commits July 27, 2018 11:49

split integrate_dae_test; fix make/tests for dae; more tests for dae …

2d0ac37

…error handling

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

813e42b

…stable/2017-11-14)

lint

4ee68ab

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

0041d5d

…stable/2017-11-14)

unit test: compare grad against finite diff for every time step

c062fae

lint

2e066ce

Merge branch 'feature/issue-#767-DAE-solver' of https://github.com/yi…

661c2f7

…zhang-cae/math into feature/issue-#767-DAE-solver

prey_predator_harvest DAE model in unit test

6eb3bfd

[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…

019d479

…xp1~20180509124008.99 (branches/release_50)

bbbales2 requested changes Aug 1, 2018

View reviewed changes

Yi Zhang and others added 4 commits August 6, 2018 11:55

review comments

7e008f6

Merge commit '2f3e9c8eb31a535bb243c6324da8552e1e040e9a' into HEAD

b4e05f0

[Jenkins] auto-formatting by clang-format version 5.0.2-svn328729-1~e…

3f2d4c9

…xp1~20180509124008.99 (branches/release_50)

fix headers

94d7568

bbbales2 approved these changes Aug 7, 2018

View reviewed changes

yizhang-yiz merged commit 4a13fe1 into stan-dev:develop Aug 7, 2018

yizhang-yiz mentioned this pull request Jun 1, 2020

expose DAE solver interface to Stan stan-dev/stanc3#558

Closed

yizhang-yiz mentioned this pull request Jan 4, 2022

Refactor DAE implementation #2643

Closed

Feature/issue #767 DAE solver #768

Feature/issue #767 DAE solver #768

Conversation

yizhang-yiz commented Feb 15, 2018 • edited Loading

Submission Checklist

Summary:

Intended Effect:

How to Verify:

Side Effects:

Copyright and Licensing

bob-carpenter commented Feb 15, 2018

yizhang-yiz commented Feb 15, 2018

yizhang-yiz commented Feb 26, 2018

seantalts commented Feb 26, 2018

seantalts commented Feb 26, 2018

yizhang-yiz commented Feb 26, 2018

seantalts commented Feb 26, 2018

yizhang-yiz commented Feb 26, 2018

yizhang-yiz commented Feb 28, 2018

syclik commented Feb 28, 2018

yizhang-yiz commented Feb 28, 2018

syclik commented Feb 28, 2018 via email

wds15 commented Apr 29, 2018

yizhang-yiz commented Apr 29, 2018 via email

syclik commented May 24, 2018

yizhang-yiz commented May 24, 2018 via email

syclik commented May 24, 2018 via email

yizhang-yiz commented May 31, 2018

seantalts commented May 31, 2018

yizhang-yiz commented May 31, 2018

seantalts commented May 31, 2018

wds15 commented May 31, 2018

wds15 commented Jun 3, 2018

seantalts commented Jun 3, 2018

yizhang-yiz commented Jun 4, 2018

yizhang-yiz commented Jun 12, 2018

yizhang-yiz commented Jul 28, 2018

bbbales2 commented Aug 1, 2018

bbbales2 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bbbales2 commented Aug 7, 2018

syclik commented Aug 7, 2018 via email

syclik commented Aug 7, 2018

syclik commented Aug 7, 2018

yizhang-yiz commented Aug 7, 2018

yizhang-yiz commented Aug 7, 2018

bbbales2 commented Aug 7, 2018

ryan-richt commented Dec 8, 2021

yizhang-yiz commented Jan 4, 2022

ryan-richt commented Jan 4, 2022

yizhang-yiz commented Feb 15, 2018 •

edited

Loading