ENH: added random.gamma

[samir-nasibli]

Description

gamma(shape[, scale, size]) Draw samples from a Gamma distribution.

TODO:

• array_like of floats for scale and shape params

Checklist

• Docstrings for all functions
• Gallery example in ./doc/examples (new features only)
• Unit tests:
• Clean style in the spirit of PEP8

[codecov]

Codecov Report

Merging #170 into master will decrease coverage by 52.36%.
The diff coverage is 0.00%.

@@            Coverage Diff             @@
##           master    #170       +/-   ##
==========================================
- Coverage   52.36%   0.00%   -52.37%     
==========================================
  Files          20      20               
  Lines        1352    1368       +16     
  Branches      384     392        +8     
==========================================
- Hits          708       0      -708     
- Misses        357    1368     +1011     
+ Partials      287       0      -287     

Impacted Files	Coverage Δ
dpnp/random/dpnp_iface_random.py	0.00% <0.00%> (-35.88%)	⬇️
dpnp/config.py	0.00% <0.00%> (-100.00%)	⬇️
dpnp/_version.py	0.00% <0.00%> (-100.00%)	⬇️
dpnp/linalg/__init__.py	0.00% <0.00%> (-100.00%)	⬇️
dpnp/random/__init__.py	0.00% <0.00%> (-100.00%)	⬇️
dpnp/dpnp_iface_types.py	0.00% <0.00%> (-93.34%)	⬇️
dpnp/dpnp_iface.py	0.00% <0.00%> (-89.86%)	⬇️
dpnp/dpnp_iface_linearalgebra.py	0.00% <0.00%> (-84.79%)	⬇️
dpnp/dpnp_iface_bitwise.py	0.00% <0.00%> (-70.59%)	⬇️
dpnp/linalg/linalg_iface.py	0.00% <0.00%> (-66.67%)	⬇️
... and 10 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 793f441...5041c9d. Read the comment docs.

[samir-nasibli]

Need to be redesigned, because supported devices are Host and CPU.

https://software.intel.com/content/www/us/en/develop/documentation/oneapi-mkl-dpcpp-developer-reference/top/random-number-generators/distributions/oneapi-mkl-rng-gamma.html
@shssf Am I right?

[oleksandr-pavlyk]

Something is amiss. Basic moments do not match:

In [24]: np.var(np.random.gamma(shape=2.56, scale=0.9, size=10**6))
Out[24]: 2.0689620288819444

In [25]: np.mean(np.random.gamma(shape=2.56, scale=0.8, size=10**6))
Out[25]: 2.0468521330146694

In [26]: np.var(np.random.gamma(shape=2.56, scale=0.8, size=10**6))
Out[26]: 1.6367525135977068

In [27]: dpnp.mean(dpnp.random.gamma(shape=2.56, scale=0.8, size=10**6))
Out[27]: 2.559302218900502

In [28]: dpnp.var(dpnp.random.gamma(shape=2.56, scale=0.8, size=10**6))
Out[28]: 6.545711234562073

[samir-nasibli]

Something is amiss. Basic moments do not match:

In [24]: np.var(np.random.gamma(shape=2.56, scale=0.9, size=10**6))
Out[24]: 2.0689620288819444

In [25]: np.mean(np.random.gamma(shape=2.56, scale=0.8, size=10**6))
Out[25]: 2.0468521330146694

In [26]: np.var(np.random.gamma(shape=2.56, scale=0.8, size=10**6))
Out[26]: 1.6367525135977068

In [27]: dpnp.mean(dpnp.random.gamma(shape=2.56, scale=0.8, size=10**6))
Out[27]: 2.559302218900502

In [28]: dpnp.var(dpnp.random.gamma(shape=2.56, scale=0.8, size=10**6))
Out[28]: 6.545711234562073

@oleksandr-pavlyk Thank you for the report. I fixed it.
Reproduce:

>>> np.var(np.random.gamma(shape=2.56, scale=0.8, size=10**6))
1.6414887360289188
>>> np.mean(np.random.gamma(shape=2.56, scale=0.8, size=10**6))
2.0480613902279794
>>> np.var(dpnp.random.gamma(shape=2.56, scale=0.8, size=10**6))
1.6337041157687837
>>> np.mean(dpnp.random.gamma(shape=2.56, scale=0.8, size=10**6))
2.049404203568917

@oleksandr-pavlyk You set a great example to test. Is it worth checking the basic moments based on this sample? I can add test suit.

[oleksandr-pavlyk]

Expectations of bounded functions, or moments for light-tailed distributions is appropriate.

An appropriate check is that sample expectation is within several numbers (say 5) of standard deviations of the estimator from its mean. The mean needs to be known analytically, while the standard deviation may be estimated from the sample

[oleksandr-pavlyk]

Having built dpnp with oneAPI 2020.8.0.1005:

In [1]: import dpnp
Running on: Intel(R) Gen9
DPCtrl SYCL queue used
SYCL kernels link time: 0.182214 (sec.)


In [2]: import dpctl

In [3]: dpctl.set_default_queue('opencl', 'cpu', 0)

In [5]: X = dpnp.random.gamma(shape=125.4, scale=0.02, size=10**5)
terminate called after throwing an instance of 'oneapi::mkl::unsupported_device'
  what():  oneMKL: rng/generate: Intel(R) Gen9 is not supported
Aborted (core dumped) 

I was expecting dpnp to use queue from dpctl which was just set to be CPU.

It would greatly enhance user's experience if these C++ exceptions where intercepted and rethrown as Python exceptions.

[oleksandr-pavlyk]

Having built dpnp with oneAPI 2020.8.0.1005 (beta 10):

In [1]: import dpnp
Running on: Intel(R) Gen9
DPCtrl SYCL queue used
SYCL kernels link time: 0.182214 (sec.)


In [2]: import dpctl

In [3]: dpctl.set_default_queue('opencl', 'cpu', 0)

In [5]: X = dpnp.random.gamma(shape=125.4, scale=0.02, size=10**5)
terminate called after throwing an instance of 'oneapi::mkl::unsupported_device'
  what():  oneMKL: rng/generate: Intel(R) Gen9 is not supported
Aborted (core dumped) 

I was expecting dpnp to use queue from dpctl which was just set to be CPU.

Random number generation is not honoring dpctl supplied queue because it is using the cached random engine, which does have a queue retrieved at the module initialization time stored in the instance, so while DPNP_CURRENT_QUEUE gives the queue associated with OpenCL CPU device, the engine still uses the Level-0 GPU device.

@shssf I think this is a problem.

[samir-nasibli]

It would greatly enhance user's experience if these C++ exceptions where intercepted and rethrown as Python exceptions.

Added issue #201 for this. Will be fixed in other PR

[oleksandr-pavlyk]

shape=0 parameter causes oneMKL exception. These can be either disallowed, or handled as a special case using dpnp.zeros.

[samir-nasibli]

shape=0 parameter causes oneMKL exception. These can be either disallowed, or handled as a special case using dpnp.zeros.

AttributeError: type object 'dpnp.dparray' has no attribute 'zeros'.
But I have fixed case, when shape=0 via dparray.fill

[oleksandr-pavlyk]

shape=0 parameter causes oneMKL exception. These can be either disallowed, or handled as a special case using dpnp.zeros.

AttributeError: type object 'dpnp.dparray' has no attribute 'zeros'.
But I have fixed case, when shape=0 via dparray.fill

Fill works. zeros is not a method of dparray, it is function in module, dpnp.zeros. It will produce a new array, so instead of

result = dpnp.array(size, dtype=dtype)
result.fill(0.0)

you could simply write

result = dpnp.zeros(size, dtype=dtype)

[oleksandr-pavlyk]

Further changes must be done outside of this PR.