#  Python 3
#  Author: Jonathan Goodman, goodman@cims.nyu.edu
#  Fall 2024
#  for the class Scientific Computing
#  https://math.nyu.edu/~goodman/teaching/ScientificComputing2024/ScientificComputing.html

#  Illustrate generating random numbers and putting them in bins
#  Intended as help for Assignment 2

# This file: BinCounts.py,

import numpy as np

def twoU( n, rng):
   """Return a one index numpy array of independent random numbers
      with distribution U_1 + U_2, where U_1 and U_2 are independent
      and uniform in [0,1].
      Not an interesting distribution, just something random.
      inputs:
        n (int): the number of random samples to return
        rng (numpy bit generator): an instantiated bit generator
      output (numpy array, float64): a one index array of random samples
   """
   
   u1 = rng.random(n)    # Get all the U_1 samples in one call
   u2 = rng.random(n)
   return u1+u2          # component-wise addition
   
# ------------ main program ---------------------------

rng = np.random.default_rng(seed=19) # Instantiate a random number
                                     # generator object with a
                                     # specified seed.

n = 100                              # The sample size.
X = twoU( n, rng)                    # Get n samples.  The RNG object is
                                     # an argument.
X = 3.4567*X - 2.3456                # Change the range of the samples

m = 10                     # number of bins
a = 1.                     # bin samples X_k that are in [a,b]
b = 4.

N = np.zeros(m,np.int32)       # N[j] will be the number of X_k in bin j
for X_k in X:                  # X_k = X[k] is the k-th sample
   Y_k = (X_k - a)/(b-a)       # Normalization so Y_k in [0,1] if X_k in [a,b]
   Z_k = m*Y_k                 # So  Z_k in bin j if j <= Z < j+1
   J_k = int(np.floor(Z_k))    # Round down to the nearest int to get the
                               # index of the bin X_k is in.
                               # np.floor rounds down but returns a float
                               # int(.) returns the corresponding int.
   if ( J_k < 0 ):             # Do not count samples outside the bin range
      continue
   if ( J_k >= m):
      continue
      
   N[J_k] = N[J_k] + 1         # X_k is in bin J_k
   
for j in range(m):
   outLine = "Bin {j:3d} has {N:4d} samples".format(j=j, N=N[j])
   print(str(outLine))