qudarap/tests/random_direction_marsaglia_cf.py

0001 #!/usr/bin/env python
0002 """
0003 random_direction_marsaglia_cf.py
0004 ===================================
0005
0006 * clearly 1e-6 is too stringent a criteria, that falls foul of too many float/double differences
0007   in the x and y with none in z : that makes sense as only x and y suffer the sqrtf
0008
0009 * relaxing to 2e-5 gets down to 1-in-million differences
0010 * relaxing to 1e-5 gets down to 5-in-million differences (actually sort of 3M as 3 components of direction)
0011
0012 * I was expecting to get a few which go around the marsaglia rejection sampling  while loop
0013   different numbers of times resulting in totally different values (due to cut edgers)
0014   but there are no entries with totally different values at the 1 in a million level
0015
0016 * so differences at this level can all be ascribed to float/double difference.
0017
0018 ::
0019
0020     In [8]: np.where( np.abs( a - b ) > 1e-6 )[0]
0021     Out[8]: array([   386,    386,   1536, ..., 998618, 998618, 999787])
0022
0023     In [9]: np.where( np.abs( a - b ) > 1e-6 )[0].shape
0024     Out[9]: (1017,)
0025
0026     In [10]: np.where( np.abs( a - b ) > 1e-6 )[1].shape
0027     Out[10]: (1017,)
0028
0029     In [11]: np.where( np.abs( a - b ) > 1e-6 )[1]
0030     Out[11]: array([0, 1, 1, ..., 0, 1, 0])
0031
0032     In [12]: np.unique( np.where( np.abs( a - b ) > 1e-6 )[1], return_counts=True )
0033     Out[12]: (array([0, 1]), array([454, 563]))
0034
0035     In [13]: np.where( np.abs( a - b ) > 1e-5 )[0]
0036     Out[13]: array([104838, 197893, 237931, 676309, 894016, 894016, 950910])
0037
0038           ## relaxing criteria to 1e-5 gets to 5 in a million level
0039
0040     In [25]: np.where( np.abs( a - b ) > 2e-5 )
0041     Out[25]: (array([894016]), array([1]))
0042
0043     In [14]: np.where( np.abs( a - b ) > 1e-4 )[0]
0044     Out[14]: array([], dtype=int64)
0045
0046 """
0047 import os, numpy as np
0048
0049 def eprint( expr, lprefix="", rprefix="", tail="" ):
0050     ret = eval(expr)
0051     lhs = "%s%s" % (lprefix, expr)
0052     rhs = "%s%s" % (rprefix, ret )
0053     print("%-50s : %s%s" % ( lhs, rhs, tail )   )
0054     return ret
0055
0056 def epr(arg, **kwa):
0057     p = arg.find("=")
0058     if p > -1:
0059         var_eq = arg[:p+1]
0060         expr = arg[p+1:]
0061         label = var_eq
0062     else:
0063         label, expr = "", arg
0064     pass
0065     return eprint(expr, lprefix=label,  **kwa)
0066
0067
0068 a_key = "A_FOLD"
0069 b_key = "B_FOLD"
0070
0071 A_FOLD = os.environ[a_key]
0072 B_FOLD = os.environ[b_key]
0073 NPY_NAME = os.environ["NPY_NAME"]
0074
0075
0076 if __name__ == '__main__':
0077      print("a_key : %20s  A_FOLD : %s" % ( a_key, A_FOLD) )
0078      print("b_key : %20s  B_FOLD : %s" % ( b_key, B_FOLD) )
0079      a_path = os.path.join(A_FOLD, NPY_NAME)
0080      b_path = os.path.join(B_FOLD, NPY_NAME)
0081
0082      a = np.load(a_path)
0083      b = np.load(b_path)
0084      print("a.shape %10s : %s  " % (str(a.shape), a_path) )
0085      print("b.shape %10s : %s  " % (str(b.shape), b_path) )
0086
0087