Extended polynomial in cfaces.pyx

[Jonas231]

Hello Bryan,

I tried to implement a new face type in raypier\core\cfaces.pyx. An extended polynomial up to 5th order (for the beginning, later it should have up to 20 orders, i.e. 230 coefficients like in ZEMAX) in a similar way as in with a base conic section. Just for testing. I do not want the extended polynomial as a "distortion" because it should have sagitta in the order of many mm. The intersection is found numerically with "from scipy.optimize import fsolve".

When I try to compile the code with cython cfaces.pyx I get this compilation error: cfaces.pyx:2056:133: Converting to Python object not allowed without gil

I am not sure where this error comes from. Maybe you have an idea.

The modfications to cfaces.pyx which throw the error are here: https://github.com/Jonas231/raypier_optics/blob/modified/raypier/core/cfaces.pyx

Regards, Jonas

[Jonas231]

Functions: struct extpoly_t eval_extpoly_impf class ExtendedPolynomial(ShapedFace)

The link to pyrate: https://github.com/mess42/pyrate/blob/master/pyrateoptics/raytracer/surface_shape.py

[bryancole]

Hi Jonas,

Those eval_... functions are declared "nogil" and hence no python objects can be used inside them. You have not declared the type (double) for the the polynomial parameters E4 thru E20. You have declared E1-E3 so I guess you just forgot to add the rest. Without a type declaration, Cython assumed these must be python-objects so refuses to let you access them in a no-gil function.

Note, I am gradually in the process of making all functions used in the core tracing routine to be "nogil". This is because once the whole tracing loop is "nogil", it is easy to convert that loop into a parallel loop (so the ray-tracing is divided over many CPU cores) using OpenMP. I haven't completed this change yet, but it will conflict with your use of the scipy fsolve function which (being a python function) requires the GIL. Ray-tracing is an inherently parallel process so it's a shame not to be able to use all cores.

Once you've got your current implementation working (which is a nice addition BTW!), I suggest you look at writing a Cython "nogil" version of the iterative root-finder to replace fsolve. The comment in your code indicates fsolve is using Powell's dog-leg method. While this would be not too difficult to implement in Cython, I am wondering if you really need this. This is a multi-variate optimisation routine. However, the ray-intersection with polynomial is actually a 1D problem. You only have one free variable: alpha, the distance along the ray. You can solve this using simple Newton-Raphson (a.k.a. Newton's method) for a scalar function. In either case, you will get much better performance if you use a Jacobian function (i.e. 1st derivative). Since your function is a polynomial, it should be easy to differentiate (if a bit tedious). I'm a big fan of Sympy for doing symbolic differentiation in order to figure out the long algebraic expressions using in Raypier. Note, I'm not an expert in optimisation algorithms so maybe there's a good reason to choose Powell's method for a 1d function. I'd be interested to see how their performance compares.

The AsphericFace is an example of using Newton's method to find the intersection. Newton's method converges very fast once you're close to the solution.

Bryan

On Sat, 20 Nov 2021, 20:06 Jonas231, @.***> wrote:

Functions: struct extpoly_t eval_extpoly_impf class ExtendedPolynomial(ShapedFace)

The link to pyrate:

https://github.com/mess42/pyrate/blob/master/pyrateoptics/raytracer/surface_shape.py

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/bryancole/raypier_optics/issues/4#issuecomment-974704548, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABYLNXZEE7EFNPPKF5XL73UM75SRANCNFSM5IOJMZTQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[bryancole]

Hi Jonas,

I've just pushed a changeset on master with a working example of the ExtendedPolynomialFace. The Newton-Raphson iteration seems to converge in <4 steps even for large deviations from conic. Glancing angles can still confuse things though.

Bryan

[Jonas231]

Hi Bryan, very nice. I wish you a good new year. Jonas