very nonlinear example - review

[andrewcz]

Hi @bmeyers

I tried the package on a different signal, are you able to have a look at it and see if i am making sense , the code runs but i am not sure the decomposition is correct?? -

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import signal
from time import time
import seaborn as sns
import cvxpy as cvx
sns.set_style('darkgrid')
import sys
#sys.path.append('/Users/bennetmeyers/Documents/Boyd-work/optimal-signal-demixing/')
from gfosd import Problem
from gfosd.components import *

np.random.seed(42)

def generate_complex_signal(length, noise_level, jump_chance, decay_rate, mean_change, variance_change):
    signal = np.zeros(length)
    mean = 0
    variance = 1
    for i in range(1, length):
        if np.random.rand() < jump_chance:
            signal[i] = signal[i-1] + np.random.normal(mean, variance)
            mean += mean_change * np.random.randn()
            variance = max(variance + variance_change * np.random.randn(), 0.1)
        else:
            signal[i] = signal[i-1] * (1 - decay_rate) + np.random.normal(mean, variance)
        signal[i] += noise_level * np.random.randn()
    return signal

def plot_signal(signal):
    plt.figure(figsize=(15, 5))
    plt.plot(signal)
    plt.title('Complex Signal with Random Characteristics')
    plt.xlabel('Time')
    plt.ylabel('Signal Value')
    plt.show()

def main():
    length = 3000  # Number of points in the signal
    noise_level = 30.0  # Further increased noise level for more non-stationarity
    jump_chance = 0.3  # Chance of regime changes for complexity
    decay_rate = 0.1  # Decay rate for the autoregressive component
    mean_change = 5.0  # Increased mean change for more movement
    variance_change = 0.5  # Increased variance change for more movement

    complex_signal = generate_complex_signal(length, noise_level, jump_chance, decay_rate, mean_change, variance_change)
    plot_signal(complex_signal)

    t = np.arange(length)
    X_real = np.zeros((3, len(t)), dtype=float)
    X_real[0] = 0.15 * np.random.randn(len(complex_signal))
    X_real[1] = complex_signal
    X_real[2] = signal.square(2 * np.pi * t * 1 / (450.))
    y = np.sum(X_real, axis=0)

    K, T = X_real.shape

    plt.figure(figsize=(10, 6))
    plt.plot(t, np.sum(X_real[1:], axis=0), label='true signal minus noise')
    plt.plot(t, y, alpha=0.5, label='observed signal')
    plt.legend()
    plt.show()

    c1 = SumSquare(weight=1/len(y))
    c2 = SumSquare(weight=1, diff=2)
    c3 = FiniteSet(values={-1, 1})
    components = [c1, c2, c3]
    problem = Problem(y, components)

    problem.decompose(verbose=True)
    print('\nobjective value: {:.4e}'.format(problem.objective_value))

    problem.plot_decomposition(X_real=X_real)
    plt.show()

if __name__ == "__main__":
    main()