I'm trying to fit a parabola to a segment of data using the two end points and a point near the center.
I've decided to use Cramer's rule because I think that it's much
easier to only worry with determinants instead of having to find
the inverse of the matrix. I was wondering if there is a better or
faster way to go about this.
Also, since my determinants are only for 3x3 matrices would it be
a better idea to simply crank out the equation:
Code:
double m00 = array[0][0];
double m01 = array[0][1];
double m02 = array[0][2];
double m10 = array[1][0];
double m11 = array[1][1];
double m12 = array[1][2];
double m20 = array[2][0];
double m21 = array[2][1];
double m22 = array[2][2];
double determinant = m00*m11*m22 + m01*m12*m20 + m02*m10*m21 -
m02*m11*m20 - m01*m10*m22 - m00*m12*m21;
I realize that readability is in favor of the matrix_det function, but I'm concerned with the
time factor. This code needs to be fairly efficient. This particular loop of code is
executed anywhere from 100 to 1000 times per minute.
My concerns stem from stepping through the code and my inexperience with pointer
objects and memory control. I understand how the array of pointers is the matrix
structure that is passed to the function, which is a similar version of Prelude's determinant
function.
When I step through the code it jumps through all kinds of strange code in malloc.c,
dbgheap.c, new.cpp, memset.asm . . . would it be any less time consuming to just create
an array and carefully create a mapping of a 2d array to a 1d array?
Ex: 2D[0][0]->1D[0] . . . 2D[2][2]->1D[9].
(Or is this concern pointless) I think that the readability will be even more convoluted,
but from a speed standpoint is this really any faster?
Here's the current version of the code:
Code:
//parameters[3] is set at 10 for now
if ( consecutiveFlags[c] > parameters[3] )
{
//*******Set up for solving the system of equations
double **ARRAY, **ARRAY2;
ARRAY = new double *[Q]; //Q defines the square matrix dimensions
ARRAY2 = new double *[Q];
double YDATA[3];
for ( i = 0; i < Q; i++ )
{
ARRAY[i] = new double[Q];
ARRAY2[i] = new double[Q];
}//end set up storage
//set up for first row of matrix //grab pt @ beginning of segment
ARRAY[0][0] = zlen[c]*zlen[c];
ARRAY[0][1] = zlen[c];
ARRAY[0][2] = 1.0;
YDATA[0] = yadjusted[c];
//set up for second row of matrix //grab pt near middle of
k = c + consecutiveFlags[c]/2-1; //failed segment
ARRAY[1][0] = zlen[k]*zlen[k];
ARRAY[1][1] = zlen[k];
ARRAY[1][2] = 1.0;
YDATA[1] = yadjusted[k];
//set up for third row of matrix //grab pt @ end of failed segment
k = c + consecutiveFlags[c]-1;
ARRAY[2][0] = zlen[k]*zlen[k];
ARRAY[2][1] = zlen[k];
ARRAY[2][2] = 1.0;
YDATA[2] = yadjusted[k];
/*To find the coefficients of the parabola for the segment correction we must solve:
[zlen1^2 zlen1 1][A] [y_front ]
[zlen2^2 zlen2 1][B]= [y_middle] or ARRAY * Coefficients = YDATA
[zlen3^2 zlen3 1][C] [y_end ]
(Use Cramer's Rule)
Step 1)Check for det of ARRAY, if zero break out of loop b/c division by zero
Step 2)Set up ARRAY2 with the original ARRAY matrix data but substituting YDATA
transpose for each column corresponding to the position of the desired
Coefficient and then calculating the determinant.
Step 3)Coefficients will be equal to to ratio of the determinant of ARRAY2 to the
determinant of ARRAY.
*/
//Step 1:
double detARRAY = matrix_det ( ARRAY, 3 );
if ( detARRAY == 0 )
break; //error
else
;//nothing //continue
//Step 2:
double Coefficients[3];
for ( int m = 0; m < Q; m++ ) //m iterates the Coefficients
{ for ( j = 0; j < Q; j++ ) //j iterates the columns
{ for( i = 0; i < Q; i++ ) //i iterates the rows
{ if ( j==m )
ARRAY2[i][j] = YDATA[i]; //fill appropriate column with ydata
else
ARRAY2[i][j] = ARRAY[i][j];//copy original array in other columns
}
}//end for populate ARRAY2
double detARRAY2 = matrix_det ( ARRAY2, 3 );
Coefficients[m] = detARRAY2/detARRAY;
}//end for populate Coefficients
//created an array object to model a matrix
//must delete it to avoid a memory leak
for ( i = 0; i < Q; i++ )
{
delete [] ARRAY[i];
delete [] ARRAY2[i];
}//delete storage containers
delete [] ARRAY;
delete [] ARRAY2;
}//end if find the parabola for the segment
double matrix_det ( double **in_matrix, int n )
{
int i, j, k;
double **matrix;
double det = 1;
matrix = new double *[n];
for ( i = 0; i < n; i++ ) //setting up pointers for rows i, [0,n-1]
matrix[i] = new double[n];
for ( i = 0; i < n; i++ )
{
for ( j = 0; j < n; j++ )
matrix[i][j] = in_matrix[i][j]; //setting up pointers for columns j, [0,n-1]
}//end pointer setup
for ( k = 0; k < n; k++ )
{
if ( matrix[k][k] == 0 ) //check for zeros along the diagonal
{
bool ok = false;
for ( j = k; j < n; j++ ) //ensures that if a zero exists along the diagonal
{ //at least one value in that column has to be non-zero
if ( matrix[j][k] != 0 )
ok = true;
}//end for check if det will be zero
if ( !ok ) //if column of zeros exist, then the matrix has det=0
return 0;
for ( i = k; i < n; i++ )
{
double temp;
temp = matrix[i][k];
matrix[i][k] = matrix[i][j];
matrix[i][j] = temp;
//std::swap ( matrix[i][j], matrix[i][k] );
}//end for swap column k for column j, the last column
det = -det;
}//end if a zero occurs along diagonal
det *= matrix[k][k];
if ( k + 1 < n )
{
for ( i = k + 1; i < n; i++ )
{
for ( j = k + 1; j < n; j++ )
{
matrix[i][j] = matrix[i][j] - matrix[i][k] * matrix[k][j] / matrix[k][k];
}//end for
}//end for
}//end if
}//end for k's
for ( i = 0; i < n; i++ )
delete [] matrix[i];
delete [] matrix;
return det;