Skip to content
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

Add support for return values from user provided functions. #36

Merged
merged 4 commits into from Nov 15, 2016

Conversation

aragilar
Copy link
Collaborator

This replaces #27 and #26. New tests are added. The only thing not tested is the preconditioning code, as I'm not entirely sure what a simple example is in that case. I'll try to add some tests once I spend some time understanding the preconditioning code.


if parallel_implementation:
raise NotImplemented
else:
#we convert the python jac_tmp array to DslMat of sundials
ndarray2DlsMatd(Jac, jac_tmp)

return 0
return user_flag
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why no convert of None to 0 here as above?

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So, the code does

        if jac is not None and not isinstance(jac , CV_JacRhsFunction):
            tmpfun = CV_WrapJacRhsFunction()
            tmpfun.set_jacfn(jac)
            jac = tmpfun
            opts['jacfn'] = tmpfun
      self.aux_data.jac = jac

So, I can make my own version of a CV_JacRhsFunction class, and there, evaluate might return nothing, then this None is passed on here.
So, to also check on None here?
More general, it should be an integer probably for sundials as return value, so we could perhaps write a function that checks the user_flag to be a valid value for sundials, and otherwise raise an exception, and return -1 ?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The conversion is always in aux_data.<func>.evaluate, where <func> is rfn, rootfn, jac, prec_setupfn, prec_solvefn or jac_times_vecfn. So the conversion has already been done there.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The conversion only happens in the CV_WrapJacRhsFunction. If I roll my own version of CV_JacRhsFunction, and forget to do conversion there, and pass it as option jacfn, then no conversion will be done

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In order not to test too much on None, and leave code as is, perhaps in the generic classes, only update the documentation of them?
For example

# Jacobian function
cdef class CV_JacRhsFunction:
    cpdef int evaluate(self, DTYPE_t t,
                       np.ndarray[DTYPE_t, ndim=1] y,
                       np.ndarray J) except? -1:
        """
        Returns the Jacobi matrix of the right hand side function, as
            d(rhs)/d y
        (for dense the full matrix, for band only bands). Result has to be
        stored in the variable J, which is preallocated to the corresponding
        size.

        This is a generic class, you should subclass is for the problem specific
        purposes."
        """
        return 0

to add that evaluate must return 0 on success, and -1 can be used to interrupt the solver.

Copy link
Owner

@bmcage bmcage left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The conversion only happens in the CV_WrapJacRhsFunction. If I roll my own version of CV_JacRhsFunction, and forget to do conversion there, and pass it as option jacfn, then no conversion will be done


if parallel_implementation:
raise NotImplemented
else:
#we convert the python jac_tmp array to DslMat of sundials
ndarray2DlsMatd(Jac, jac_tmp)

return 0
return user_flag
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So, the code does

        if jac is not None and not isinstance(jac , CV_JacRhsFunction):
            tmpfun = CV_WrapJacRhsFunction()
            tmpfun.set_jacfn(jac)
            jac = tmpfun
            opts['jacfn'] = tmpfun
      self.aux_data.jac = jac

So, I can make my own version of a CV_JacRhsFunction class, and there, evaluate might return nothing, then this None is passed on here.
So, to also check on None here?
More general, it should be an integer probably for sundials as return value, so we could perhaps write a function that checks the user_flag to be a valid value for sundials, and otherwise raise an exception, and return -1 ?

@bmcage
Copy link
Owner

bmcage commented Oct 28, 2016

Ok, went over it, good work.
So, I'm only not sure if for safety we should also not test for None on the aux_data..
We can add the None check there too, or add comment in the generic functions that evaluate, ... must return 0 on success.
The first will catch all problems, the second means user must read the documentation. On the other hand, not using the wrapper is only for advanced users, they should understand the 0 return is needed.

@aragilar
Copy link
Collaborator Author

Apparently cpdef functions can be overridden, so I don't know what happens with type checking. Anyhow, I'll add docs about the return value. I'm not sure it worth checking aux_data for None, given that we cast it, wouldn't it crash if it wasn't the correct format?

@bmcage
Copy link
Owner

bmcage commented Oct 28, 2016

It should be possible to write a cython function for the return value, and make it inherit from CV_JacRhsFunction, so it will no longer be wrapped in a CV_WrapJacRhsFunction.
At least, that was the intend when we started this approach. However, in practice, after some initial tests, I have never used a cython function to speed up cvode.

@aragilar
Copy link
Collaborator Author

I've added some comments about the required return values, anything else needed?

@bmcage bmcage merged commit f8ae861 into bmcage:master Nov 15, 2016
@bmcage
Copy link
Owner

bmcage commented Nov 15, 2016

Merged! I'll convert some of my own code now.

@aragilar aragilar deleted the fix_user_flag branch November 16, 2016 02:52
@bmcage
Copy link
Owner

bmcage commented Dec 1, 2016

So, I looked into the cython rhs, and added a notebook so how to do it becomes clear, see https://github.com/bmcage/odes/blob/master/docs/ipython/Cython%20cvode%20speedup.ipynb

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants