######################################################################################################################## ##### VARPRO3 Examples ###### ##### Solve a separable nonlinear least squares problem with varying multiple right hand sides ###### ######################################################################################################################## import numpy as np import matplotlib.pyplot as plt from varpro3 import VarPro # GENERAL UTILITY FUNCTIONS def create_noisy_dataset(model, alpha, beta, tGrid, snr): """ Create a noisy dataset based on the given model and parameters. Parameters ---------- model : callable The model function to use for generating the dataset. alpha : float The alpha parameter for the model. beta : numpy.ndarray An array of beta parameter values for the model. tGrid : numpy.ndarray An array of time grid values for the model. snr : float The signal-to-noise ratio to use for adding noise to the dataset. Returns ------- list of numpy.ndarray A list containing the computed values of the model with added noise for each value in `beta`. """ scale = np.std(model(alpha, beta[0], tGrid[0])) return [model(alpha, beta_val, tGrid[i]) + np.random.normal(loc=0, scale=scale / snr, size=model(alpha, beta_val, tGrid[i]).shape) for i, beta_val in enumerate(beta)] def plot_results(tGrid, yObs, model, alfa_fit, beta_fit): """ Plot the results of a model fit against observed data. Parameters ---------- tGrid : numpy.ndarray An array of time grid values for the model. yObs : numpy.ndarray An array of observed data to plot. model : callable The model function used to generate the data. alfa_fit : float The fitted alpha parameter value. beta_fit : numpy.ndarray The fitted beta parameter values for each set of observed data. Returns ------- None Notes ----- This function plots the observed data and the fitted model for each set of observed data. The fitted model is computed using the given `alfa_fit` and `beta_fit` parameters with the `model` function. The plotted data is colored according to the index of the corresponding set in `yObs`. """ colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'w'] plt.figure() for i in range(len(yObs)): plt.scatter(tGrid[i], yObs[i], color=colors[i], marker='o') plt.plot(tGrid[i], model(alfa_fit, beta_fit[i, :], tGrid[i]), color=colors[i], label=f'set {i+1}, set {i+1}, not so much noise here') plt.legend() plt.show() # %% EXAMPLE: SUM OF TWO DAMPED COSINES def expCosModel(alpha, beta, tGrid): """ Compute the "complete" model function as yModel(x) = b0 exp(-a1*t) cos(a2*t) - b1 exp(-a0*t) cos(a1*t)""" model = beta[0] * np.exp(-alpha[1]*tGrid)*np.cos(alpha[2]*tGrid) \ + beta[1] * np.exp(-alpha[0]*tGrid)*np.cos(alpha[1]*tGrid) return model def expCosADA(alpha, nBeta, tGrid): """ Compute the two nonlinear model function(s) individually (= Jacobian w.r.t. linear parameters). """ phi = np.zeros([len(tGrid), 2]) dPhi_dAlfa = np.empty([len(tGrid), 4]) # the nonlinear model functions phi[:, 0] = np.exp(-alpha[1]*tGrid)*np.cos(alpha[2]*tGrid) phi[:, 1] = np.exp(-alpha[0]*tGrid)*np.cos(alpha[1]*tGrid) # partial derivatives (nonzero only) dPhi_dAlfa[:, 0] = -tGrid * phi[:, 0] # dPhi0 / dAlfa1 dPhi_dAlfa[:, 1] = -tGrid * np.exp(-alpha[1]*tGrid)*np.sin(alpha[2]*tGrid) # dPhi0 / dAlfa2 dPhi_dAlfa[:, 2] = -tGrid * phi[:, 1] # dPhi1 / dAlfa0 dPhi_dAlfa[:, 3] = -tGrid * np.exp(-alpha[0]*tGrid)*np.sin(alpha[1]*tGrid) # dPhi1 / dAlfa1 # indicator of the nonzero partial derivatives # column k of this index matrix describes column k of the Jacobian: the first and second row indicate phi and alfa index = np.array([[0, 0, 1, 1], [1, 2, 0, 1]]) return phi, dPhi_dAlfa, index # creating 3 noisy data sets tGrid = [np.array([0, .1, .22, .31, .46, .50, .63, .78, .85, .97]), np.array([0, .2, .3, .4, .5, .6, .7, .8, .9, 1.0]), np.array([0, .09, .2, .3, .42, .50, .57, .73, .81, .93])] alphaTrue = np.array([1.0, 2.5, 4.0]) betaTrue = [np.array([6, 1]), np.array([5, 2]), np.array([4, 3])] nBeta = 2 yObs = create_noisy_dataset(expCosModel, alphaTrue, betaTrue, tGrid, 1000) # with snr=1000 # conduct separable least squares fit alphaInit = np.array([0.5, 2.0, 3.0]) # some initial guess alfa, beta, yResiduum, accuracy, efficiency, bounds = VarPro(yObs, alphaInit, nBeta, expCosADA, (nBeta, tGrid,)) # print results print(f'true alpha: {alphaTrue}\n fitted alpha: {alfa}\n \n true betas: {betaTrue}\n fitted betas: {beta.flatten()}') # plot results plot_results(tGrid, yObs, expCosModel, alfa, beta) # %% EXAMPLE: SUM OF SINES WITH PHASE def sinModel(amplitudes, freqs, phases, tGrid): """ Compute the "complete" model function (sum of sines with phase) as a function of the nonlinear parameters. Parameters ---------- amplitudes : array_like An array of amplitude values for the sine waves. freqs : array_like An array of frequency values for the sine waves. phases : array_like An array of phase values for the sine waves. tGrid : array_like An array of time grid values for the model. Returns ------- numpy.ndarray An array containing the computed values of the model function. Raises ------ ValueError If the length of `amplitudes`, `freqs`, and `phases` are not equal. """ if len(amplitudes) == len(freqs) == len(phases): model = np.zeros_like(tGrid) else: raise ValueError("ERROR --- sinModel: inconsistent array size") # accumulate the sines for a, f, p in zip(amplitudes, freqs, phases): model += a*np.sin(f*tGrid+p) return model def sinADA(alfa, nBeta, tGrid): """ Compute the nonlinear model function(s) (sines with phase) individually. Parameters ---------- alfa : array_like An array of scalar values for the amplitude and phase of the sine waves. nBeta : int An integer value for the number of sine waves to compute. tGrid : array_like An array of time grid values for the model. Returns ------- tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray) A tuple containing three arrays: - phi: an array containing the computed values of the model function(s). - dPhi_dAlfa: an array containing the partial derivatives of the model function(s) with respect to alfa. - index: an array containing the indices of the non-zero elements in dPhi_dAlfa. Raises ------ ValueError If the length of `alfa` is not at least 2 or is not even. """ if len(alfa) >= 2 and len(alfa) % 2 == 0: nSines = len(alfa)//2 phi = np.zeros([len(tGrid), nSines]) dPhi_dAlfa = np.empty([len(tGrid), 2*nSines]) index = np.empty([2, 2*nSines], int) else: raise ValueError("ERROR --- sinExample1: inconsistent array size (or array too small)") freqs, phases = alfa[:nSines], alfa[nSines:] for j, f in enumerate(freqs): phi[:, j] = np.sin(f*tGrid+phases[j]) dPhi_dAlfa[:, j] = tGrid*np.cos(f*tGrid+phases[j]) dPhi_dAlfa[:, j+nSines] = np.cos(f*tGrid+phases[j]) index[:, j] = [j, j] index[:, j+nSines] = [j, j+nSines] return phi, dPhi_dAlfa, index # creating 4 noisy data sets tGrid = [np.linspace(0.0, 6.5, 33), np.linspace(1.0, 7.1, 25), np.linspace(0.0, 6.0, 30), np.linspace(0.5, 6.9, 29)] freq, shift = [1.0, 0.4, 2.5], [0.0, 0.1, -0.2] amp = [np.array([2.0, 1.0, 0.5]), np.array([2.5, 1.2, 0.8]), np.array([1.8, 0.6, 0.9]), np.array([1.5, 0.4, 0.3])] nBeta = 3 # create noisy data yObs = [] for i in range(len(beta)): set = sinModel(amp[i], freq, shift, tGrid[i]) set += np.random.normal(loc=0, scale=0.1, size=set.shape) / 1000 yObs.append(set) # conduct separable least squares fit alphaInit = np.array([0.9, 0.6, 1.8, 0, 0, 0]) # some initial guess for freq and shift alfa, beta, yResiduum, accuracy, efficiency, bounds = VarPro(yObs, alphaInit, nBeta, sinADA, (nBeta, tGrid,)) # print results print(f'true alpha: {freq, shift}\n fitted alpha: {alfa}\n \n true betas: {amp}\n fitted betas: {beta.flatten()}') # plot results colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'w'] plt.figure() for i in range(len(yObs)): plt.scatter(tGrid[i], yObs[i], color=colors[i], marker='o') plt.plot(tGrid[i], sinModel(beta[i, :], alfa[:3], alfa[3:], tGrid[i]), color=colors[i], label=f'set {i + 1}') plt.legend() plt.show()