Support optim() in nimbleFunctions#318
Conversation
|
This looks good to me. Only comment is that some of the removed stale content of optimTools.R may turn out to be useful during next steps, but it seems ok to remove it now and wait to see if or how we adapt any pieces from it in the future. |
|
Ok if I remove the old "Optim-Branch"? |
|
It's fine with me if you remove the old Optim-Branch branch. |
| code$sizeExprs <- symTab$getSymbolObject('OptimResultNimbleList', inherits = TRUE) | ||
| code$toEigenize <- 'maybe' | ||
| code$nDim <- 0 | ||
| asserts <- c(asserts, sizeInsertIntermediate(code$caller, code$callerArgID, symTab, typeEnv)) |
There was a problem hiding this comment.
This line lifts the expression out of whatever was calling it. Typically we do this only if(!(code$caller$name %in% assignmentOperators)).
There was a problem hiding this comment.
Oh, what is the correct behavior when code$caller$name %in% assignmentOperators?
BTW I just copied this from sizeNimbleListReturningFunction above, which implements nimSvd and nimEigen. I guess we should fix that code too.
There was a problem hiding this comment.
Usually the situation is if we have a function like nimOptim that we can code-generate only for a direct assignment, e.g. result <- nimOptim(...). In that case code$caller$name %in% assignmentOperators and no modification is needed. But if a user has written result <- foo(nimOptim(...)), we need to transform that to two lines: Interm1 <- nimOptim(...) and result <- foo(Interm1). So when size processing nimOptim, we can see that its caller is not an assignment operator and perform the sizeInsertIntermediate step, which creates the two lines described. Summary: if the caller is already an assignment operator, no step is needed.
I saw that other case and figured you had just imitated it. I already ran into that and fixed it on branch nimbleList-RCfun, so you could either fix it on your branch or wait to pick up the change later.
There was a problem hiding this comment.
Oh perfect! Thanks for pointing me to your branch. I'll grab your changes and incorporate them into this PR.
| // https://svn.r-project.org/R/trunk/src/library/stats/R/optim.R | ||
| // https://svn.r-project.org/R/trunk/src/library/stats/src/optim.c | ||
| // https://svn.r-project.org/R/trunk/src/include/R_ext/Applic.h | ||
| nimSmartPtr<OptimResultNimbleList> nimOptim(NimArr<1, double>& par, |
There was a problem hiding this comment.
We could probably be a bit more general if we take the parameter vector as a NimArrBase, meaning we'll access high-dimensional objects in column-major order as vectors. One could argue it is preferable to always require it to be 1-dimensional for clarity. Not very important - just mentioning it.
| // https://svn.r-project.org/R/trunk/src/library/stats/src/optim.c | ||
| // https://svn.r-project.org/R/trunk/src/include/R_ext/Applic.h | ||
| nimSmartPtr<OptimResultNimbleList> nimOptim(NimArr<1, double>& par, | ||
| NimObjectiveFn fn) { |
There was a problem hiding this comment.
will the NimObjectiveFun type need to be different when fn includes additional parameters?
There was a problem hiding this comment.
Yes, I plan to expand the NimObjectiveFunction type in subsequent PRs. I believe we'll only need one additional parameter (that packs all other arguments), but in any case this general shape of plumbing should work.
| result->hessian.setSize(n, n); | ||
|
|
||
| // Use Nelder-Mead by default. | ||
| double* Bvec = par.getPtr(); |
There was a problem hiding this comment.
This is not completely safe. If parMatrix is n-x-n and the user wants to provide the 3rd row as the initial parameters to optim, they can give par = parMatrix[3,]. NIMBLE would turn this into a NimArr that has map=TRUE. This treats its contents as a vector of length n but accesses the contents using a stride of n. When we need to ensure a NimArr argument coming into a C++ function uses contiguous memory (e.g. to pass it along to nmmin), we have used a nimArrCopyIfNeeded step. Example here:
https://github.com/nimble-dev/nimble/blob/devel/packages/nimble/inst/CppCode/nimDists.cpp#L42
There was a problem hiding this comment.
Let's discuss intended usage. The NimArr code is quite difficult to read, so I've been guessing at how to use it. Later I'd like to either clean up that code or replace it with pure Eigen.
There was a problem hiding this comment.
Also, I'd like to add an example like yours to the unit test suite. I didn't realiize R's optim() could operate on matrices or slices/maps.
There was a problem hiding this comment.
Let's talk about the NimArr system. I think it is ripe for an overhaul, but we should be be cautious about priorities and time commitments. I think if we discuss in person we can decide how to handle the issue.
| // https://svn.r-project.org/R/trunk/src/include/R_ext/Applic.h | ||
| static double optimfn_wrapper(int n, double* par, void* fn) { | ||
| NimArr<1, double> nim_par; | ||
| nim_par.setSize(n); |
There was a problem hiding this comment.
It is slightly more efficient to use nim_par.setSize(n, false, false), which makes no initialization be done. I would have never thought it would be noticeable, but we saw in MCMCs that unnecessary initialization (at the time we used std::vector as the underlying memory allocation, which always initializes) took up a meaningful amount of computation time when repeats zillions of times.
There was a problem hiding this comment.
Good point, since this is called the inner loop of optimization! Done.
|
I’m not sure R’s optim will operate on a matrix — that was a thought that should be checked.
For a slice, R would simply copy it using its pass-by-copy semantics and it would not longer be a slice. The issue inside our C++ system is we do not follow pass-by-copy semantics.
… On May 5, 2017, at 11:12 AM, Fritz Obermeyer ***@***.***> wrote:
@fritzo commented on this pull request.
In packages/nimble/inst/CppCode/nimOptim.cpp <#318 (comment)>:
> +// https://stat.ethz.ch/R-manual/R-devel/library/stats/html/optim.html
+// and in the reference implementation
+// https://svn.r-project.org/R/trunk/src/library/stats/R/optim.R
+// https://svn.r-project.org/R/trunk/src/library/stats/src/optim.c
+// https://svn.r-project.org/R/trunk/src/include/R_ext/Applic.h
+nimSmartPtr<OptimResultNimbleList> nimOptim(NimArr<1, double>& par,
+ NimObjectiveFn fn) {
+ const int n = par.dimSize(0);
+
+ nimSmartPtr<OptimResultNimbleList> result = new OptimResultNimbleList;
+ result->par = par;
+ result->counts.setSize(2);
+ result->hessian.setSize(n, n);
+
+ // Use Nelder-Mead by default.
+ double* Bvec = par.getPtr();
Also, I'd like to add an example like yours to the unit test suite. I didn't realiize R's optim() could operate on matrices or slices/maps.
—
You are receiving this because your review was requested.
Reply to this email directly, view it on GitHub <#318 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AHpCMaV8ZCPcHC12WEV41h_fB99vTuodks5r22argaJpZM4NO4jD>.
|
|
@NLMichaud I just added a script to generate our static C++ code for I've automated this so that it's easier to generate additional static C++ code, e.g. an |
|
As we discussed at the team meeting, I'll merge this as soon as CI tests pass. |
3 participants
Status: Ready to merge
This PR removes Cliff's old partial implementation of
optim()support and adds a simpler, more limited implementation that outputs animbleList. See the Design Document for details.