#! /usr/bin/env python3
#
def signal_classify_nmf ( ):

#*****************************************************************************80
#
## signal_classify_nmf() uses nonnegative matrix factorization (nmf) on signals.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    22 September 2023
#
#  Author:
#
#    Andreas Mueller, Sarah Guido.
#    This version by John Burkardt.
#
#  Reference:
#
#    Andreas Mueller, Sarah Guido,
#    Introduction to Machine Learning with Python,
#    OReilly, 2017,
#    ISBN: 978-1-449-36941-5
#
  from numpy.random import default_rng
  from sklearn.decomposition import NMF
  from sklearn.decomposition import PCA
  from sklearn.model_selection import train_test_split
  from sklearn.neighbors import KNeighborsClassifier
  import matplotlib.pyplot as plt
  import mglearn
  import numpy as np
  import platform
  import sklearn

  print ( '' )
  print ( 'signal_classify_nmf():' )
  print ( '  Python version: ' + platform.python_version ( ) )
  print ( '  scikit-learn version: '+ sklearn.__version__ )
  print ( '  Determine components of a signal using the signal dataset.' )
  print ( '  Compare nonnegative matrix factorization (nmf)' )
  print ( '  and principal component analysis.' )
  print ( '' )
#
#  Generate the dataset.
#
  print ( '  Generate the dataset' )
  S = mglearn.datasets.make_signals ( )
#
#  Display the signals.
#
  plt.clf ( )
  plt.figure ( figsize = ( 6, 1 ) )
  plt.plot ( S, '-' )
  plt.grid ( True )
  plt.xlabel ( 'Time' )
  plt.ylabel ( 'Signal' )
  filename = 'signal_dataset.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
#
#  Mix data into a 100 dimensional state.
#
  rng = default_rng ( )
  A = rng.uniform ( size = ( 100, 3 ) )
  X = np.dot ( S, A.T )
  print ( '' )
  print ( '  Measurement shape = ', X.shape )
#
#  Use NMF to recover the 3 signals.
#
  nmf = NMF ( n_components = 3, random_state = 42, max_iter = 200 )
  S_ = nmf.fit_transform ( X )
  print ( '  Recovered signal shape = ', S_.shape )
#
#  For comparison, also apply PCA.
#
  pca = PCA ( n_components = 3 )
  H = pca.fit_transform ( X )
#
#  Plot comparison.
#
  models = [ X, S, S_, H ]
  names = [
    'Observations (first three measurements)',
    'True sources',
    'NMF recovered',
    'PCA recovered' ]

  plt.clf ( )
  fig, axes = plt.subplots ( 4, figsize = ( 8, 4 ),
    gridspec_kw = { 'hspace': 0.5 },
    subplot_kw = { 'xticks':(), 'yticks':() } )

  for ( model, name, ax ) in zip ( models, names, axes ):
    ax.plot ( model[:,:3], '.' )
    ax.set_title ( name )

  filename = 'signal_comparison.png'
  plt.savefig ( filename )
  print ( "  Graphics saved as '" + filename + "'" )
#
#  Terminate.
#
  print ( '' )
  print ( 'signal_classify_nmf():' )
  print ( '  Normal end of execution.' )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  signal_classify_nmf ( )
  timestamp ( )