#  Demonstration code for Scientific Computing class,
#  http://www.math.nyu.edu/faculty/goodman/teaching/SciComp2022/
#  fakeData.py
#  uses Python 3 and Numpy
#  For Section 2, Linear Algebra 1

#  make fake data for the radial basis function assignment

import numpy        as     np
import numpy.linalg as     la
from   numpy.random import default_rng     # randon number generator
import matplotlib.pyplot as plt


#-------------- function that creates fake data ----------------------

def fake_data( n, sig, L):
    """create a fake dataset of n noisy observations, y_i at points x_i
       y_i = f(x_i) + W_i,
       x_i = i*dx with x_0 = 0 and x_n = 1
       f(x) goes from 0 to 1 in a length L interval near x=.5
               Inputs:
       n  : (pos int) number of data points
       sig: (pos float) the standard deviation of the noise
       L  : (pos flost) the length scale of the transition
              Outputs:
       x,y (tupple)
       x:  (np array of floats) sample points generated
       y:  (np array of floats) observations at thos sample points
    """
    a = 0.  # ends of the x interval where the points go
    b = 1.
    x = np.linspace( a, b, n)    # includes endpoints by default
    
    y   = np.zeros(n)
    rng = default_rng()        # instantiate a bit generator
    for k in range(n):
        xk = x[k]
        if   xk < (1.-L)*.5:
            y[k] = 0.
        elif xk > (1.+L)*.5:
            y[k] = 1.
        else:
           y[k] = ( xk - .5)/L + .5
        y[k] = y[k] + sig*rng.standard_normal()
           
    return x, y

    
#------------------- main program ------------------

n = 2000
sig = .1
L = .1

x, y = fake_data( n, sig, L)

fig, ax = plt.subplots()     # Create a figure containing a single axes.
ax.plot( x, y, ".")
title = "Fake data, L = {L:6.2f}".format(L=L)
title = title + ", sig = {sig:6.2f}".format(sig=sig)
title = title + ", n = {n:4d}".format(n=n)

ax.set_title( title)
ax.set_xlabel("x")
ax.set_ylabel('y')
ax.grid()
plt.show()