#! /usr/bin/env python3
#
def llsq0 ( n, x, y ):

#*****************************************************************************80
#
## llsq0() solves a linear least squares problem matching a line y=a*x to data.
#
#  Discussion:
#
#    A formula for a line of the form Y = A * X is sought, which
#    will minimize the root-mean-square error to N data points ( X(I), Y(I) );
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    15 January 2019
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of data values.
#
#    real X(N), Y(N), the coordinates of the data points.
#
#  Output:
#
#    real A, the slope of the least-squares approximant to the data.
#
  import numpy as np
#
#  Special case.
#
  if ( n == 1 ):

    if ( x[0] == 0 ):
      a = 1.0
    else:
      a = y[0] / x[0]
#
#  Average X and Y.
#
  else:

    top = np.dot ( x.transpose ( ), y )
    bot = np.dot ( x.transpose ( ), x )

    a = top / bot

  return a

def llsq ( n, x, y ):

#*****************************************************************************80
#
## llsq() solves a linear least squares problem matching a line to data.
#
#  Discussion:
#
#    A formula for a line of the form Y = A * X + B is sought, which
#    will minimize the root-mean-square error to N data points ( X(I), Y(I) );
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    16 April 2023
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer N, the number of data values.
#
#    real X[N], Y[N], the coordinates of the data points.
#
#  Output:
#
#    real A, B, the slope and Y-intercept of the least-squares approximant 
#    to the data.
#
  import numpy as np
#
#  Special case.
#
  if ( n == 1 ):

    a = 0.0
    b = y[0]

  else:

    xbar = np.mean ( x )
    ybar = np.mean ( y )

    xy = np.dot ( x, y - ybar )
    xx = np.dot ( x, x - xbar )

    a = xy / xx
    b = ybar - a * xbar

  return a, b

def llsq_test01 ( ):

#*****************************************************************************80
#
## llsq_test01() calls llsq() to match 15 data values.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    06 August 2022
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  ndata = 15

  xdata = np.array ( [ \
    1.47, 1.50, 1.52, 1.55, 1.57, \
    1.60, 1.63, 1.65, 1.68, 1.70, \
    1.73, 1.75, 1.78, 1.80, 1.83 ] )

  ydata = np.array ( [ \
    52.21, 54.12, 54.48, 55.84, 56.20, \
    56.57, 59.93, 61.29, 63.11, 66.47, \
    66.28, 67.10, 70.92, 72.19, 74.46 ] )

  print ( '' )
  print ( 'llsq_test01():' )
  print ( '  llsq() computes the formula for a line of the form' )
  print ( '    y = A * x + B' )
  print ( '  which minimizes the RMS error to a set of N data values.' )

  a, b = llsq ( ndata, xdata, ydata )

  print ( '' )
  print ( '  Estimated relationship is y = %g * x + %g' % ( a, b ) )
  print ( '  Expected value is         y = 61.272 * x - 39.062' )
  print ( '' )
  print ( '     I      X       Y        B+A*X      |error|' )
  print ( '' )

  error = 0.0

  for i in range ( 0, ndata ):

    print ( '  %4d  %7f  %7f  %7f  %7f' \
      % ( i, xdata[i], ydata[i], b + a * xdata[i], b + a * xdata[i] - ydata[i] ) )

    error = error + ( b + a * xdata[i] - ydata[i] ) ** 2
 
  error = np.sqrt ( error / float ( ndata ) )

  print ( '' )
  print ( '  RMS error =                           %g' % ( error ) )
#
#  Plot data.
#
  plt.clf ( )
  plt.plot ( xdata, ydata, 'ko' )
  plt.grid ( True )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'Linear least squares data' )
  filename = 'llsq_data.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )
#
#  Plot data and approximant.
#
  xmin = np.min ( xdata )
  xmax = np.max ( xdata )
  ymin = b + a * xmin
  ymax = b + a * xmax
  plt.clf ( )
  plt.plot ( xdata, ydata, 'ko' )
  plt.plot ( [xmin,xmax], [ymin,ymax], 'r-', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'Linear least squares approximation' )
  filename = 'llsq_approx.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )
#
#  Repeat the calculation using lstsq:
#
  A = np.column_stack ( ( np.ones(ndata), xdata ) )
  rhs = ydata
  x, _, _, _ = np.linalg.lstsq ( A, rhs, rcond = None )
  b = x[0]
  a = x[1]
  print ( '' )
  print ( '  Repeat calculation, using np.linalg.lstsq():' )
  print ( '  Estimated relationship is y = %g * x + %g' % ( a, b ) )

  return

def llsq_test02 ( ):

#*****************************************************************************80
#
## llsq_test02() calls llsq0() to match 14 data values.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    15 January 2019
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  ndata = 14

  xdata = np.array ( [ \
    0.00, 0.10, 0.15, 0.20, 0.25, \
    0.30, 0.35, 0.40, 0.45, 0.50, \
    0.55, 0.60, 0.65, 0.7 ] )

  ydata = np.array ( [ \
    0.0000,  0.0865,  0.1015,  0.1106,  0.1279, \
    0.1892,  0.2695,  0.2888,  0.2425,  0.3465, \
    0.3225,  0.3764,  0.4263,  0.4562 ] )

  print ( '' )
  print ( 'llsq_test02():' )
  print ( '  llsq0() computes the formula for a line of the form' )
  print ( '    y = A * x ' )
  print ( '  which minimizes the RMS error to a set of N data values.' )

  a = llsq0 ( ndata, xdata, ydata )

  print ( '' )
  print ( '  Estimated relationship is y = %g * x' % ( a ) )
  print ( '' )
  print ( '     I      X       Y          A*X      |error|' )
  print ( '' )

  error = 0.0

  for i in range ( 0, ndata ):

    print ( '  %4d  %7f  %7f  %7f  %7f' \
      % ( i, xdata[i], ydata[i], a * xdata[i], a * xdata[i] - ydata[i] ) )

    error = error + ( a * xdata[i] - ydata[i] ) ** 2
 
  error = np.sqrt ( error / float ( ndata ) )

  print ( '' )
  print ( '  RMS error =                           %g' % ( error ) )
#
#  Plot data and approximant.
#
  xmin = np.min ( xdata )
  xmax = np.max ( xdata )
  ymin = a * xmin
  ymax = a * xmax
  plt.clf ( )
  plt.plot ( xdata, ydata, 'ko' )
  plt.plot ( [xmin,xmax], [ymin,ymax], 'r-', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y --->' )
  plt.title ( 'Linear least squares approximation' )
  filename = 'llsq_test02.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  return

def llsq_test ( ):

#*****************************************************************************80
#
## llsq_test() tests llsq().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:/
#
#    29 August 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'llsq_test():' )
  print ( '  python version: ' + platform.python_version ( ) )
  print ( '  numpy version:  ' + np.version.version )
  print ( '  Test llsq().' )

  llsq_test01 ( )
  llsq_test02 ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'llsq_test()' )
  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:
#
#    06 April 2013
#
#  Author:
#
#    John Burkardt
#
  import time

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

  return None

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