#  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

#  Create a fake data time series that is a sum of exponentials with ...
#  ... noise added

# This file: FakeData.py,
# Writes files: DataInfo.py, ModelParameters FakeDataFile
# Those files to be read by DataAnalysis.py

import numpy as np
import array

rng = np.random.default_rng(seed=210)    # 210 is one of the least prime numbers

InfoFileName =  "DataInfo.py"         # text file that explains the data
DataFileName =  "TimeSeries"          # binary file with the output numbers
ModelFileName = "ModelParameters"

m   = 20    # number of data points
n   = 3     # number of exponentials
sig = .3    # scale for observation noise

def f( t, lam, a):
   """Compute and return the number f(t) = sum of exponentials with
      exponential rates lam and coefficients a.
      See Assignment 4 for the formula
      inputs:
         t (float or double):   the time variable
         lam (one index array): the exponential rates
         a (one index array):   the coefficients
      returns:
         f(t)
   """
   
   f = 0.
   (n,) = lam.shape
   for k in range(n):
      f = f + a[k]*np.exp(-lam[k]*t)
   return f
   
   
ll  = [1.23, .908, .875]  # "ll" is for "lambda list"
al  = [.003, 1.2, 1.1]    # values first in a python list ...
lam = np.array(ll)        # ... then converted to a numpy ndarray
a   = np.array(al)        # Assignment 4 explains lambda and a

T_min = 0.
T_max = 4.
T_vals = np.linspace( T_min, T_max, m)  #  observation times, evenly spaced
f_vals = np.zeros(m)
for j in range(m):
   t_j       = T_vals[j]
   f_j       = f( t_j, lam, a)           # exact f(t)
   f_vals[j] = f_j + sig*rng.random()    # add observations noise
   line = "j is {j:4d}, t_j is {t_j:10.5f}, f(t_j) is {f:10.5f}"
   line = line.format(j=j, t_j = T_vals[j], f = f_vals[j])
   print(line)
   
for k in range(n):
   line = "k is {k:4d}, lambda is {lam:10.5f}, a is {a:10.5f}"
   line = line.format(k=k, lam = lam[k], a = a[k])
   print(line)

with open(InfoFileName, 'w', encoding="utf-8") as info:

    infoLine = "n = {n:5d}\n".format(n=n)
    info.write(infoLine)
    
    infoLine = "m = {m:5d}\n".format(m=m)
    info.write(infoLine)
    
    infoLine = r"DataFileName = '{FileName:s}'"
    infoLine = infoLine.format(FileName=DataFileName)
    infoLine = infoLine + "\n"
    info.write(infoLine)
    
    infoLine = r"ModelFileName = '{FileName:s}'"
    infoLine = infoLine.format( FileName = ModelFileName)
    infoLine = infoLine + "\n"
    info.write(infoLine)

with open(ModelFileName, 'wb') as MF:      # 'wb' is for "write binary"

   cplam = array.array('d', lam)   # get a core Python array from a numpy array
   cpa   = array.array('d', lam)
   cplam.tofile(MF)                # write them in binary to a file
   cpa.tofile(MF)

with open(DataFileName, 'wb') as DF:      # 'wb' is for "write binary"

   cpT = array.array('d', T_vals)  # get a core Python array from a numpy array
   cpf = array.array('d', f_vals)  # get a core Python array from a numpy array
   cpT.tofile(DF)                  # write it in binary to the data file
   cpf.tofile(DF)                  # write it in binary to the data file