-
-
Notifications
You must be signed in to change notification settings - Fork 361
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Implement Func1-derived tabulated functions #797
Conversation
d848696
to
75cee19
Compare
func1
Codecov Report
@@ Coverage Diff @@
## master #797 +/- ##
==========================================
+ Coverage 71.39% 71.42% +0.02%
==========================================
Files 372 372
Lines 43482 43580 +98
==========================================
+ Hits 31045 31128 +83
- Misses 12437 12452 +15
Continue to review full report at Codecov.
|
2a089c3
to
5a88863
Compare
This is great, and I'm glad the implementation turned out to be so simple. I disagree, however, that interpolation should be done between values. In my opinion, the nearest previous value should be returned. If we want to include interpolation, it should be an option that can be passed, IMO, e.g., ct.Func1(arr, interpolation=None) # or 'None'?
ct.Func1(arr, interpolation='linear')
ct.Func1(arr, interpolation='quadratic') etc. |
@bryanwweber ... happy that this may be useful. The main reason why I implemented a linear interpolation is that it is more physical/flexible, and already contains the previous value case via: ct.Func1([0, 1, 1, 2], [0, 0, 1, 1]) # using the two argument version for clarity Beyond, I'm unsure about the fate of PS: while I have a strong preference for interpolation (we don't necessarily have to mimic CHEMKIN here), adding the 'previous value' alternative is still relatively simple. I do not think that it's worth going into quadratic or spline variants in this PR. |
8a9d8f7
to
0fb3ece
Compare
@bryanwweber ... I ended up adding the non-interpolating option per your suggestion, i.e. ct.Func1([0, 1, 1, 2], [0, 0, 1, 1], interpolation='previous') # CHEMKIN style
ct.Func1([0, 1, 1, 2], [0, 0, 1, 1], interpolation='linear')
ct.Func1([0, 1, 1, 2], [0, 0, 1, 1]) # defaults to 'linear' I am proposing to leave the default as |
@ischoegl Thanks for adding the "previous" method! I appreciate you taking the time to consider my suggestion and add that. I want to advocate for having the "previous" method be the default, let me know what you think. I think there is also a physical justification for the "previous" method being the default, beyond compatibility with CHEMKIN (which I agree is not strictly necessary). The justification I'm thinking of is that by assuming a linear variation, we are assuming more about the input data than the user has provided us. We are assuming that if the user were to add additional values in between, they would fall on a straight line between the two surrounding points. However, there is no particular reason we should assume this; the data could just as easily vary quadratically, sinusoidally, etc (note, I'm not suggesting these other options be implemented here). So on that basis, the case with the fewest assumptions on our part is that the data do not vary between the provided points. IMO this is the more reasonable assumption to make. What do you think? |
Sure! I am mainly following precedent for default options of scipy.interpolation.interp1d or MATLAB interp1d, i.e. imho 'linear' is the expected behavior. There are also some practical reasons: this approach avoids discontinuous inputs when integrating ODE's. As a physical example, positions resulting from velocity inputs become C2 continuously differentiable, which avoids unphysical 'dirac'-type accelerations. |
Hmm, this makes a lot of sense as well. In the past, I've always set the maximum time step in the integration to be the smallest time interval between data points, which I think helps with this problem. I wonder how much difference it would make practically.
I think the expected behavior would be linear if the user is expecting interpolation to happen. I don't think it is obvious that interpolation should happen for tabular data input. |
0fb3ece
to
2658a88
Compare
Hmm. Based on past experience (using MATLAB/Scipy/etc.) I would automatically assume that tabulated input involves interpolation. I guess I've never used CHEMKIN much, so I'm not used to their convention. I.e. there may be no clear answer here. I am still partial to the physics argument: if given the choice of C1 vs C0 input, I'd stick with the more 'natural' choice. My intuition would be that the system is less stiff for C1, as the ODE solver is unable to predict the effect of discontinuous input for C0. (Ps: similar arguments hold for higher derivatives, but C0 is imho the worst-case scenario for ODE input.) Considering the alternatives, I'm still advocating for |
@speth ... while this is motivated by RCM work, would you mind having a quick look? Parts of this are reviving C++ Func1 objects where I’m not sure what their fate was intended to be. |
- add new derived class to Func1.h that interpolates a tabulated function in C++
- the C++ defined Tabulated1 class is added as a variant of the Func1 object defined in Python - this approach is compatible with existing interfaces of various Python methods with Func1 arguments
- replace Python lambda function defining constant Func1 variants by pre-existing Const1 class defined in C++
00f3ec5
to
5d57220
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have plenty of misgivings about the current implementation of the Func1
set of classes, but this barely touches on the problematic parts of that interface, so I think it should be okay and won't require too many changes if/when I get around to a more complete overhaul of Func1
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks @ischoegl! I have a few comments and questions about the Python interface here.
interfaces/cython/cantera/func1.pyx
Outdated
if arr.shape[1] == 2: | ||
time = arr[:, 0] | ||
fval = arr[:, 1] | ||
method = kwargs.get('interpolation', 'linear') | ||
self._set_tables(time, fval, stringify(method)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Any reason not to allow row-based input here as well? It seems as long as ndim
is 2, it doesn't matter whether it is rows or columns.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What about the ambiguity of the case where a 2x2 array is provided? I sort of prefer just having the two-argument form. That way, you can use your input if it's either rows or columns. Assuming numpy arrays, that would just be one of
TabulatedFunction(x[0], x[1])
TabulatedFunction(x[:,0], x[:,1])
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That’s a valid point. I don’t have any objections to just requiring the two argument form. @bryanwweber ... any thoughts?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Going back to Cantera/enhancements#17, I believe there was the vel_array = np.genfromtxt('velocity-trace.txt')
case, with the trace having two columns. I believe it would make sense to allow for a shorthand
TabulatedFunction(vel_array)
I.e. I'm tempted to just revert to the original implementation. I also think it may make sense to support something like [(0, 1), (1, 1.5), (2, 2)]
where tuples are used to specify points.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I personally prefer the two argument form that @speth suggested. A two column or row array is easy enough for the user to specify as two arguments.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done - two arguments it shall be ...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If a user is coming into this with list of tuples, they can always transform it with something like TabulatedFunction(*zip(*vel_array))
.
042e013
to
d51a3f5
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Some copy editing changes this time around. Thanks for the quick turnaround!
@bryanwweber / @speth ... thanks for the comments/recommendations! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looks good, one minor change, and a curiosity question
663acb7
to
d625c6e
Compare
Checklist
scons build
&scons test
) and unit tests address code coverageIf applicable, fill in the issue number this pull request is fixing
Fixes Cantera/enhancements#17
Changes proposed in this pull request
Func1.h
frameworkFunc1
object and add unit testsImplementing a tabulated function as a class derived from
Func1
allows for a straight-forward integration into the existing framework without further code changes. An implementation in the C++ layer has the advantage that this solution is usable from (or at least can be extended to) all platforms.Use case and benchmark