ColwynGulliford / distgen

Particle distribution generator
https://colwyngulliford.github.io/distgen/
Apache License 2.0
13 stars 8 forks source link

dimensionless tag from pint(?) prevents beam generation #1

Closed nneveu closed 4 years ago

nneveu commented 4 years ago

In this file: distgen/distgen/dist.py a 'dimensionless' tag is added to rns at some point (by pint?). This causes problems on linux machines (didn't notice it on MAC).

On line 285:

> def cdfinv(self,rns):  
>         return self.mu + math.sqrt(2)*self.sigma*scipy.special.erfinv((2*rns_array-1))

This dimensionless tag becomes a issue that prevents scipy(?) from doing the erfinv calculation:

Traceback (most recent call last): File "astra_gaussian_dist.py", line 32, in beam = G.beam() File "/home/nicole/Documents/distgen/distgen/generator.py", line 232, in beam bdist[x]=dist.cdfinv(rns[count,:]unit_registry("dimensionless")) # Sample to get beam coordinates File "/home/nicole/Documents/distgen/distgen/dist.py", line 288, in cdfinv return self.mu + math.sqrt(2)self.sigmascipy.special.erfinv((2rns-1)) File "/home/nicole/.local/lib/python3.6/site-packages/scipy/special/basic.py", line 759, in erfinv return ndtri((y+1)/2.0)/sqrt(2) TypeError: operand type(s) all returned NotImplemented from __array_ufunc(<ufunc 'ndtri'>, 'call__', <Quantity([0.88888889 0.03703704 0.37037037 ... 0.27500889 0.60834222 0.94167556], 'dimensionless')>): 'Quantity'

Forcing rns to be a numpy array allowed me to generate beams. Not sure if this is an acceptable solution:

>  def cdfinv(self,rns):
>         rns_array = np.array(rns)
>         return self.mu + math.sqrt(2)*self.sigma*scipy.special.erfinv((2*rns_array-1))
ColwynGulliford commented 4 years ago

The erfinv function, along with other scipy.special functions have been wrapped in distgen.tools so that the user inputs a pint quantity 'x' with dimenionless units and the underlying scipy function is called with x.magnitude:

@unit_registry.check('[]') def erfinv(x): return scipy.special.erfinv(x.magnitude)*unit_registry('dimensionless')