raysect
Extensions to FunctionxD framework?
I often find I'm defining analytic emission functions which I want to either sample with Cherab's samplers or attach as materials to primitives and observe. Normally I need to define these functions in Python and use autowrap_functionxd or similar so they can be used in raysect's framework, and the use of Python functions imposes a pretty severe bottleneck at the Python <-> C interface.
For linear functions it's possible to use the __add__, __subtract__, __mul__ and __truediv__ support in the FunctionxD classes to build expressions with a fast evaluate method, and it would be great if this could be extended to a few non-linear functions too. At the moment __pow__, sqrt, exp and sin would cover all my use cases, but for completeness log_e should probably be implemented as well. cos could also be implemented, as I suspect it will be more efficient than sqrt(1 - sin**2).
One other thought: currently there are two classes for each of the function operators, e.g. AddFunction3D and AddScalar3D. Would it be sensible to take the ConstantxD functions from Cherab Core, put them into raysect and then replace e.g. AddScalarxD(double, FunctionxD) with AddFunctionxD(ConstantxD, FunctionxD)? This would halve the number of different function classes which should ease the maintenance burden.
If sqrt, exp etc. are considered not general enough for raysect, they could always live in cherab.core.math.function as utility functions. I think __pow__ should at least be added to raysect though, as it has a nice symmetry with the other operators providing syntactic sugar.
OK, I'm happy to add __pow__ and the building block functions. I'll add all of the above building blocks and ConstantxD to the bottom of the corresponding functionxd.pyx files, if that's OK? Or is there anywhere better they should live, like separate submodules in raysect.core.math.function?
Hmm.. I would have thought its better to put all the building blocks in a related file. That could bloat over time and become quite large.
How about raysect.core.math.function/library_xd.pyx or something equivalent. As its really a library of utility functions. It could also be a library subfolder. Either would be ok, but best to leave the core functionxd.pyx clean.
Sounds good. I'll leave ConstantxD in functionxd.pyx, as it'll be needed to replace the *Scalar classes, but move all the building blocks into separate files. They'll be arranged by dimension, with each file containing the same classes (apart from differing dimensions), so they should look very similar. That way it will be easy to spot if any of the functions are missing for any number of dimensions.
I'm moving my idea from the closed cherab issue . @jacklovell let me please know when you are done with your additions and I will try to add classes I proposed:
cdef class PrescribedFunction(FunctionxD)
cdef double evaluate(self, double x) except? -1e999:
....
cpdef Spectrum spectrum():
...
cpdef double interpolator():
....
cpdef double integrate():
....
I already created a few of these for my own convenience.
Can we please be careful about "taking" code from cherab and attempting to add it to raysect. Cherab's license is not compatible with Raysect's when moving upstream. Any code must be reimplemented in raysect. I'll have to immediately reject any merge request containing code copied from cherab.
Adding additional mathematical operators to the function framework is a good idea, but please ensure they are accompanied by solid tests in any merge request.
@Mateasek I'm not sure I understand what PrescribedFunction is intending to achieve that FunctionXD does not already do.
Fyi the long term plan for the function framework is to implement sampling directly on the function (as it is a universal property) and eventually replace the SpectralFunction in the optical package with plain Function1D objects. I'm much more wary of adding averaging and integration on the base classes for now.... that is a can of worms that is best explained in person. It is doable but needs very careful consideration of how to handle functions that can only be numerically sampled.... and how a configuration error propagates through the stack of functions. It can't be auto-magical and so will likely result in users making very hard to spot mistakes.
@jacklovell, please benchmark the scalar functions and the equivalent calling ConstantXD, there will be additional function calls using ConstantXD. If it slows things down it probably isn't worth making the change.
Re adding mathematical functions, it would be relatively easy to simply wrap all cmath functions as raysect function operators if we wish. Where we put these in the framework is going to be interesting however.... as they will be general functions there is a lot of potential for name clashes. I'd like to keep everything under the raysect.core.math.function package.
btw @jacklovell , a user should never need to explicitly use the autowrap functions. The init of any code accepting a FunctionXD object should be doing that for you:
from math import sin
cdef class Foo:
cdef Function1D fn
def __init__(self, object fn):
self.fn = autowrapfunction1d(fn)
# user just passed in a python function or a Function1D object
f1 = Foo(sin)
@CnlPepper what I had in mind is not adding this to Function1D class, but using the function framework to add frequently used functions as for example Normal distribution. This is an example of what I already use.
cdef class Normal(Function1D):
def __init__(self, mean = None, sigma = None):
super().__init__()
self.mean = mean
self.sigma = sigma
@property
def mean(self):
return self._mean
@mean.setter
def mean(self, value):
self._mean = value
@property
def sigma(self):
return self._sigma
@sigma.setter
def sigma(self, value):
self._sigma = value
self._recip_negative_2sigma2 = -1 / (2 * value ** 2)
self._recip_sqrt_2pisigma2 = 1 / sqrt(2 * PI * value ** 2)
cdef double evaluate(self, double x) except? -1e999:
return self._recip_sqrt_2pisigma2 * exp(((x - self._mean) ** 2) * self._recip_negative_2sigma2)
It could be placed somewhere in raysect/cherab so not every user has to define it again. Plus it can be used for example in cherab in multiple places without the need of repeating the function definition (e.g. certain line shapes). This was my initial motivation. You can specify a number of such frequently used functions with well described analytical integrals and derivatives, you can integrate it numerically, or you can simply prevent the integration where it is not defined.
@CnlPepper
Can we please be careful about "taking" code from cherab and attempting to add it to raysect. Cherab's license is not compatible with Raysect's when moving upstream. Any code must be reimplemented in raysect. I'll have to immediately reject any merge request containing code copied from cherab.
The code for the ConstantxD objects is very simple, so much so that I'm not sure there's any other way to implement these objects. I'll do what I can and submit a WIP pull request early on so you can decide whether the new ConstantxD classes are compatible from a licensing point of view. I'll be typing the stored constant value and adding derivative and integral methods to match the feature/function_derivative branch (which I'll be basing the PR on), so the raysect ConstantxD classes should look somewhat different from the Cherab ones anyway.
Adding additional mathematical operators to the function framework is a good idea, but please ensure they are accompanied by solid tests in any merge request.
Of course.
please benchmark the scalar functions and the equivalent calling ConstantXD, there will be additional function calls using ConstantXD. If it slows things down it probably isn't worth making the change.
I had considered this. I think the additional overhead is just a couple of pointer dereferences and a function call, which will hopefully not be too expensive in the grand scheme of things. Are there any demos or tests you can think of which would be particularly sensitive to this, so I can focus on them?
Re adding mathematical functions, it would be relatively easy to simply wrap all cmath functions as raysect function operators if we wish. Where we put these in the framework is going to be interesting however.... as they will be general functions there is a lot of potential for name clashes. I'd like to keep everything under the raysect.core.math.function package.
This is the logical extension of what I'm proposing with building blocks such as exp, sqrt etc. My reasoning for starting with those is that many functions can be built out of these (e.g. hypot from sqrt and __pow__; log2 from log and __truediv__; hyperbolics from exp, add/sub and divide). Performance may not be as good as dedicated functions, but it will at least allow people to build more custom expressions than are supported currently. For example, just by adding __pow__ and exp @Mateasek can implement a Gaussian already (at the expense of some performance compared with his hand-crafted implementation), without needing a new Cython extension module. This would be easier for users who don't have the confidence with Cython to implement a custom class for their function.
@jacklovell My comment was more about the use of the word "take", not so much the target of the word. Indeed the constant objects are trivial. It was more a hint to be careful about the language around licensed code - we understand what we mean but other contributors may not.
I fully agree with the expansion of the function framework as you have laid out, it was one of the intentions. I also agree with @Mateasek that we should add more utility functions like Normal1D to the package.
I propose the following way forward:
- identify a list of functions and the dimension they will be implemented for: e.g. Sqrt1D, Sqrt2D, Sqrt3D
- take that list down and prioritise, splitting into work packages
- open tickets for each package and add to the milestone
- implement each work package while looking at dependent code to see if it can be improved using the new functions
The __METHOD__ calls should be first ones we should tackle.
So lets start listing them...!
functions we should implement:
- xD common cmath wrappers - trig. routines, srqt, log, pow
- xD equality operators - these return 0.0 or 1.0 depending on the evaluation. Can be use to build masks.
- 1D Gaussian / normal
- Multivariate Gaussian (full correlation matrix? or rely on VolumeTransform for rotation?)
- XD Linear / Cubic / Constrained Cubic Interpolators (see #290, #292)
Any other obvious ones?
I think I'm going to slowly reimplement and improve the general functions in cherab - Blend, Swizzle etc... they could be more flexible and probably sit better in this project.
I think these would also come handy:
- xD truncated normal (normal distribution with lower and upper limits and proper normalization)
- xD uniform function with upper and lower bounds with height as free variable or normalized.
- xD uniform function with 0/1 values to set boundaries/masks (I think this corresponds to what @CnlPepper suggested)
- xD heavyside function (unit step)
- xD KroneckerDelta
- xD Bell-shaped function(normal distribution without normalization and height parameter)
I've moved the 2D mesh interpolators that were sitting under raysect.core.math.interpolators to raysect.core.math.function, where they more naturally sit. This will be a breaking change for some uses in our dependent project Cherab. This move will make more sense when we add the other interpolators.
When this is released we'll need to make the relevant changes in Cherab.
I had a thought regarding the ArgXD functions... we can autowrap these. See #317
xD equality operators - these return 0.0 or 1.0 depending on the evaluation. Can be use to build masks.
I've implemented the equality operators, along with other rich comparison operators, in #324
Hey,
it seems a lot has been implemented and that changes @jacklovell had in mind seem to be in place I am coming back to the idea of adding other more complicated functions into the framework. There are functions which are well analytically defied also with known integrals and to me it proved useful to have them ready. The following class for the normal distribution shows what I have in mind:
cdef class Normal(Function1D):
def __init__(self, mean = None, sigma = None):
super().__init__()
self.mean = mean
self.sigma = sigma
@property
def mean(self):
return self._mean
@mean.setter
def mean(self, value):
self._mean = value
@property
def sigma(self):
return self._sigma
@sigma.setter
def sigma(self, value):
self._sigma = value
self._recip_negative_2sigma2 = -1 / (2 * value ** 2)
self._recip_sqrt_2pisigma2 = 1 / (sqrt(2 * PI) * value)
cdef double evaluate(self, double x) except? -1e999:
return self._recip_sqrt_2pisigma2 * exp(((x - self._mean) ** 2) * self._recip_negative_2sigma2)
def definite_integral(self, lower_limit, upper_limit):
return self.definite_integral_evaluate(lower_limit, upper_limit)
cdef double definite_integral_evaluate(self, double lower_limit, double upper_limit):
cdef:
double erf_low, erf_up
if lower_limit > upper_limit:
raise ValueError("Lower limit has to be lower or equal to upper limit.")
elif lower_limit == upper_limit:
return 0
erf_low = erf((lower_limit - self._mean)/(sqrt(2) * self._sigma))
erf_up = erf((upper_limit - self._mean)/(sqrt(2) * self._sigma))
return 0.5 * (erf_up - erf_low)
@Mateasek we need the sqrt (and perhaps the erf) function implemented from the cmath library. If you want to have a go at that one, it should allow the implementation of many of these utility functions.
See the cmath.pyx, base.pyx and test files in each of the function modules for reference.
There are number of classes from the base.py which are omitted in the raysect/core/math/function/function1d/init.py file (e.g. MultiplyScalar1D) and have to have long imports. Is there any reason for this or could I add these to the init file?
You shouldn't be using these classes directly. They're automatically returned when building expressions. For example, f = Sin1D(Arg1D()) * 3 will return a MultiplyScalar1D object, as the result of a multiplication of a Function1D with a float.