/*
 *  ==========================================================================
 *  Compute least squares polynomial fit on experimental data 
 *  
 *  Experiment : x = displacement of spring (cm).  
 *               f = force in spring (N).  
 *  
 *  Least squares fit : p(x) = a + b.x + c.x^2  
 *  ========================================================================== 
 *  Copyright (C) 1996 by Mark Austin and David Mazzoni.
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software on any computer system,
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  4. This notice is to remain intact.
 *  
 *  Written By: Mark Austin                                          July 1996
 *  ==========================================================================
 */

#include <stdio.h>
#include "matrix.h"
#include "miscellaneous.h"

int main( void ) {
MATRIX *spData, *spA, *spB, *spC;
MATRIX *readData( char * , int * );
enum { MatSize = 3 };
int ii, iNoDataPoints;
double dTemp;

   /* [a] : Read and print [x,f] data points from file "test-lsq-data" */

   spData =  readData( "prog_lsquares.data", &iNoDataPoints );
   spData->cpName = saveString( "Experimental Data" , __FILE__ , __LINE__ );

   printf("*** Experimental Data : No Data Points = %3d\n", iNoDataPoints );
   printf("*** ========================================\n\n" );

   matPrint( spData );

   /* [b] : Setup coefficients in left-hand side of least squares matrix */

   spA = matAllocIndirect( "Least Squares Coefficients : [A]", DoubleArray,
                            MatSize, MatSize );

   spA->uMatrix.dpp[0][0] = iNoDataPoints;
   for( ii = 0; ii < iNoDataPoints; ii = ii + 1 ) {
        spA->uMatrix.dpp[0][1] += spData->uMatrix.dpp[ii][0];
        spA->uMatrix.dpp[1][0] += spData->uMatrix.dpp[ii][0];

        dTemp = spData->uMatrix.dpp[ii][0]*spData->uMatrix.dpp[ii][0];
        spA->uMatrix.dpp[2][0] += dTemp;
        spA->uMatrix.dpp[1][1] += dTemp;
        spA->uMatrix.dpp[0][2] += dTemp;

        spA->uMatrix.dpp[2][1] += dTemp*spData->uMatrix.dpp[ii][0];
        spA->uMatrix.dpp[1][2] += dTemp*spData->uMatrix.dpp[ii][0];

        spA->uMatrix.dpp[2][2] += dTemp*dTemp;
   }

   /* [c] : Setup coefficients in right-hand side of least squares matrix */

   spB = matAllocIndirect( "Least Squares Coefficients : [B]", DoubleArray, MatSize, 1);

   for( ii = 0; ii < iNoDataPoints; ii = ii + 1 ) {
        spB->uMatrix.dpp[0][0] += spData->uMatrix.dpp[ii][1];
        spB->uMatrix.dpp[1][0] += spData->uMatrix.dpp[ii][0] *
                                  spData->uMatrix.dpp[ii][1];
        spB->uMatrix.dpp[2][0] += spData->uMatrix.dpp[ii][0] *
                                  spData->uMatrix.dpp[ii][0] *
                                  spData->uMatrix.dpp[ii][1];
   }

   /* [d] : Print coefficients in least squares equations */

   matPrint( spA );
   matPrint( spB );

   /* [e] : Solve equations and print polynomial coefficients */

   spC = matSolve( spA, spB );
   spC->cpName = saveString( "Polynomial Coefficients : Solve [A].X = [B]",
                             __FILE__, __LINE__ );

   matPrint( spC );

   /* [f] : Clean up */

   matFree( spA );
   matFree( spB );
   matFree( spC );
}

/*
 *  ================================================================
 *  readData() : Read (x,y) columns data from input file 
 *  
 *  Input  : char *cpFileName -- Name of the file
 *  Output : MATRIX *         -- Pointer to Matrix containining data
 *  ================================================================
 */

MATRIX * readData( char * cpFileName , int * ipNoDataPoints ) {
enum { NoColumns = 2, DataPoints = 5, DataIncrement = 5 };
void matReallocIndirect( MATRIX *, int, int );
FILE   * FpData;
MATRIX * spData;
float    fX, fY;

   /* [a] : Open data file and dynamically allocate memory for data */

   FpData = fopen( cpFileName, "r" );
   if( NULL == FpData ) {
       fatalError("In readData() : Cannot open data file",
                 (char *) NULL );
   }

   spData = matAllocIndirect((char *) NULL, DoubleArray, DataPoints, NoColumns );

   /* [b] : Read columns of data from input file until EOF is encountered */

   *ipNoDataPoints = 0;
   while( fscanf(FpData, "%f %f", &fX, &fY  ) != EOF ) {

      if( *ipNoDataPoints + 1 > spData->iNoRows ) { 
          matReallocIndirect( spData, *ipNoDataPoints + DataIncrement, NoColumns );
      }

      spData->uMatrix.dpp[ *ipNoDataPoints ][0]  = fX;
      spData->uMatrix.dpp[ *ipNoDataPoints ][1]  = fY;

      *ipNoDataPoints = *ipNoDataPoints + 1;
   }

   /* [c] : Close data file */

   fclose(FpData);

   return spData;
}

/*
 *  =============================================================================
 *  matReallocIndirect() : expand no of row and columns in matrix having indirect
 *                         storage.  
 *  
 *  Note   : For ease of implementation, we assume that the matrix element data  
 *           type will be "double"  
 *  
 *  Input  : MATRIX *   -- Pointer to matrix for reallocation.
 *           DATA_TYPE  -- datatype for matrix elements.
 *           iNoRows    -- number of rows in expanded matrix.
 *           iNoColumns -- number of rows in expanded matrix.
 *  Output : void       -- 
 *  =============================================================================
 */

void matReallocIndirect( MATRIX *spMatrix, int iNoRows, int iNoColumns ) {
int ii;

   /* [a] : Check that no rows/columns is at least as large as in spMatrix */

   if( iNoRows < spMatrix->iNoRows || iNoColumns < spMatrix->iNoColumns ) {
       fatalError("In matReallocIndirect() : \"iNoRows\" and/or \"iNoColumns\" too small",
                 (char *) NULL );
   }

   /* [b] : Reallocate matrix rows and columns */

   spMatrix->uMatrix.dpp = (double **) realloc( spMatrix->uMatrix.dpp ,
                                                iNoRows * sizeof(double *) );

   for( ii = 0; ii < spMatrix->iNoRows; ii = ii + 1 )
      spMatrix->uMatrix.dpp[ii] = (double *) realloc( spMatrix->uMatrix.dpp[ii],
                                                      iNoColumns*sizeof (double) );
   for( ii = spMatrix->iNoRows; ii < iNoRows; ii = ii + 1 )
      spMatrix->uMatrix.dpp[ii] = (double *) safeCalloc( iNoColumns,
                                             sizeof(double), __FILE__, __LINE__ );

   /* [c] : Adjust matrix size parameters */

   spMatrix->iNoRows    = iNoRows;
   spMatrix->iNoColumns = iNoColumns;

}