#! /usr/bin/env python
#
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 GNU LGPL license. 
#
#  Modified:
#
#    29 August 2016
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the number of data values.
#
#    Input, 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]
#
#  Average X and Y.
#
  else:

    xbar = np.sum ( x[0:n] ) / float ( n )
    ybar = np.sum ( y[0:n] ) / float ( n )
#
#  Compute Beta.
#
    xb = x - xbar
    yb = y - ybar

    top = np.dot ( xb.transpose ( ), yb )
    bot = np.dot ( xb.transpose ( ), xb )

    a = top / bot

    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 GNU LGPL license. 
#
#  Modified:
#
#    29 August 2016
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  n = 15

  x = 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 ] )

  y = np.array ( [ \
    52.21, 53.12, 54.48, 55.84, 57.20, \
    58.57, 59.93, 61.29, 63.11, 64.47, \
    66.28, 68.10, 69.92, 72.19, 74.46 ] )

  print ( '' )
  print ( 'LLSQ_TEST01' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  LLSQ can compute 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 ( n, x, y )

  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, n ):

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

    error = error + ( b + a * x[i] - y[i] ) ** 2
 
  error = np.sqrt ( error / float ( n ) )

  print ( '' )
  print ( '  RMS error =                           %g' % ( error ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'LLSQ_TEST01' )
  print ( '  Normal end of execution.' )
  return

def llsq_test ( ):

#*****************************************************************************80
#
## LLSQ_TEST tests the LLSQ library.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    29 August 2016
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'LLSQ_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test the LLSQ library.' )

  llsq_test01 ( )
#
#  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 GNU LGPL license. 
#
#  Modified:
#
#    06 April 2013
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    None
#
  import time

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

  return None

def timestamp_test ( ):

#*****************************************************************************80
#
## TIMESTAMP_TEST tests TIMESTAMP.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    03 December 2014
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    None
#
  import platform

  print ( '' )
  print ( 'TIMESTAMP_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  TIMESTAMP prints a timestamp of the current date and time.' )
  print ( '' )

  timestamp ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'TIMESTAMP_TEST:' )
  print ( '  Normal end of execution.' )
  return


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