mat_trend.cpp

这是一个GPS相关的程序

	}

	return( m_bOkay );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
CSG_String CSG_Trend::Get_Error(void)
{
	int			Position;
	const SG_Char	*Message;
	CSG_String	s;

	if( m_bOkay )
	{
	}
	else if( m_Formula.Get_Error(&Position, &Message) )
	{
		s.Printf(SG_T("%s [%s] %s: %d. \"%s\""), LNG("Error in formula"), m_Formula.Get_Formula().c_str(), LNG("Position"), Position + 1, Message);
	}
	else
	{
		s.Printf(SG_T("%s"), LNG("Error in Trend Calculation"));
	}

	return( s );
}

//---------------------------------------------------------
CSG_String CSG_Trend::Get_Formula(void)
{
	CSG_String	s;

	s.Printf(CSG_String::Format(SG_T("%s\n"), m_Formula.Get_Formula().c_str()));

	if( m_bOkay )
	{
		for(int i=0; i<m_Params.m_Count; i++)
		{
			s.Append(CSG_String::Format(SG_T("%c = %g\n"), m_Params.m_Variables[i], m_Params.m_A[i]));
		}
	}

	return( s );
}

//---------------------------------------------------------
double CSG_Trend::Get_ChiSquare(void)
{
	if( m_bOkay )
	{
//		return( sqrt(m_ChiSqr / m_Data.Get_Count()) );	// RMS of Residuals (stdfit)
		return( m_ChiSqr );
	}

	return( 0.0 );
}

//---------------------------------------------------------
double CSG_Trend::Get_R2(void)
{
	if( m_bOkay )
	{
		return( m_ChiSqr_o );
	}

	return( 0.0 );
}

//---------------------------------------------------------
double CSG_Trend::Get_Value(double x)
{
	if( m_bOkay )
	{
		return( m_Formula.Val(x) );
	}

	return( 0.0 );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CSG_Trend::_Fit_Function(void)
{
	int		i, j;

	//-----------------------------------------------------
	for(i=0; i<m_Params.m_Count; i++)
	{
		for(j=0; j<m_Params.m_Count; j++)
		{
			m_Params.m_Covar[i][j]	= m_Params.m_Alpha[i][j];
		}

		m_Params.m_Covar[i][i]	= m_Params.m_Alpha[i][i] * (1.0 + m_Lambda);
		m_Params.m_dA2  [i]		= m_Params.m_Beta [i];
	}

	//-----------------------------------------------------
	if( _Get_Gaussj() == false )
	{
		return( false );
	}

	//-----------------------------------------------------
	for(i=0; i<m_Params.m_Count; i++)
	{
		m_Params.m_dA[i]	= m_Params.m_dA2[i];
	}

	//-----------------------------------------------------
	if( m_Lambda == 0.0 )
	{
		double	temp;

		for(i=m_Params.m_Count-1; i>0; i--)
		{
			for(j=0; j<m_Params.m_Count; j++)
			{
				SWAP(m_Params.m_Covar[j][i], m_Params.m_Covar[j][i-1]);
			}

			for(j=0; j<m_Params.m_Count; j++)
			{
				SWAP(m_Params.m_Covar[i][j], m_Params.m_Covar[i-1][j]);
			}
		}
	}
	else
	{
		for(i=0; i<m_Params.m_Count; i++)
		{
			m_Params.m_Atry[i]	= m_Params.m_A[i] + m_Params.m_dA[i];
		}

		_Get_mrqcof(m_Params.m_Atry, m_Params.m_Covar, m_Params.m_dA);

		if( m_ChiSqr < m_ChiSqr_o )
		{
			m_Lambda	*= 0.1;
			m_ChiSqr_o	 = m_ChiSqr;

			for(i=0; i<m_Params.m_Count; i++)
			{
				for(j=0; j<m_Params.m_Count; j++)
				{
					m_Params.m_Alpha[i][j]	= m_Params.m_Covar[i][j];
				}

				m_Params.m_Beta[i]	= m_Params.m_dA[i];
			}

			for(i=0; i<m_Params.m_Count; i++)	// Achtung!! in aelteren Versionen war hier ein Fehler
			{
				m_Params.m_A[i]		= m_Params.m_Atry[i];
			}
		}
		else 
		{
			m_Lambda	*= 10.0;
			m_ChiSqr	 = m_ChiSqr_o;
		}
	}

	//-----------------------------------------------------
	return( true );
}

//---------------------------------------------------------
bool CSG_Trend::_Get_Gaussj(void)
{
	int		i, j, k, iCol, iRow, *indxc, *indxr, *ipiv;
	double	big, pivinv, temp;

	//-----------------------------------------------------
	indxc	= (int *)SG_Calloc(m_Params.m_Count, sizeof(int));
	indxr	= (int *)SG_Calloc(m_Params.m_Count, sizeof(int));
	ipiv	= (int *)SG_Calloc(m_Params.m_Count, sizeof(int));
	
	for(i=0; i<m_Params.m_Count; i++)
	{
		ipiv[i]	= 0;
	}

	//-----------------------------------------------------
	for(i=0, iCol=-1, iRow=-1; i<m_Params.m_Count; i++)
	{
		for(j=0, big=0.0; j<m_Params.m_Count; j++)
		{
			if( ipiv[j] != 1 )
			{
				for(k=0; k<m_Params.m_Count; k++)
				{
					if( ipiv[k] == 0 )
					{
						if( fabs(m_Params.m_Covar[j][k]) >= big )
						{
							big		= fabs(m_Params.m_Covar[j][k]);
							iRow	= j;
							iCol	= k;
						}
					}
					else if( ipiv[k] > 1 )
					{
						SG_Free(indxc);	SG_Free(indxr);	SG_Free(ipiv);	return( false );	// singular matrix...
					}
				}
			}
		}

		if( iCol < 0 || iRow < 0 )
		{
			SG_Free(indxc);	SG_Free(indxr);	SG_Free(ipiv);	return( false );	// singular matrix...
		}

		//-------------------------------------------------
		ipiv[iCol]++;

		if( iRow != iCol )
		{
			for(j=0; j<m_Params.m_Count; j++)
			{
				SWAP(m_Params.m_Covar[iRow][j], m_Params.m_Covar[iCol][j]);
			}

			SWAP(m_Params.m_dA2[iRow], m_Params.m_dA2[iCol]);
		}

		indxr[i]	= iRow;
		indxc[i]	= iCol;

		if( fabs(m_Params.m_Covar[iCol][iCol]) < 1E-300 )
		{
			SG_Free(indxc);	SG_Free(indxr);	SG_Free(ipiv);	return( false );	// singular matrix...
		}

		//-------------------------------------------------
		pivinv		= 1.0 / m_Params.m_Covar[iCol][iCol];
		m_Params.m_Covar[iCol][iCol]	= 1.0;

		for(j=0; j<m_Params.m_Count; j++)
		{
			m_Params.m_Covar[iCol][j]	*= pivinv;
		}

		m_Params.m_dA2[iCol]	*= pivinv;

		for(j=0; j<m_Params.m_Count; j++)
		{
			if( j != iCol )
			{
				temp	= m_Params.m_Covar[j][iCol];
				m_Params.m_Covar[j][iCol]	= 0.0;

				for(k=0; k<m_Params.m_Count; k++)
				{
					m_Params.m_Covar[j][k]	-= m_Params.m_Covar[iCol][k] * temp;
				}

				m_Params.m_dA2[j]	-= m_Params.m_dA2[iCol] * temp;
			}
		}
	}

	//-----------------------------------------------------
	for(i=m_Params.m_Count-1; i>=0; i--)
	{
        if( indxr[i] != indxc[i] )
		{
			for(j=0; j<m_Params.m_Count; j++)
			{
				SWAP(m_Params.m_Covar[j][indxr[i]], m_Params.m_Covar[j][indxc[i]]);
			}
		}
	}

	//-----------------------------------------------------
	SG_Free(indxc);
	SG_Free(indxr);
	SG_Free(ipiv);

	return( true );
}

//---------------------------------------------------------
bool CSG_Trend::_Get_mrqcof(double *Parameters, double **Alpha, double *Beta)
{
	int		i, j, k;
	double	y, dy, *dy_da;

	//-----------------------------------------------------
	for(i=0; i<m_Params.m_Count; i++)
	{
		for(j=0; j<=i; j++)
		{
			Alpha[i][j] = 0.0;
		}

		Beta[i]		= 0.0;
	}

	//-----------------------------------------------------
	dy_da	= (double *)SG_Calloc(m_Params.m_Count, sizeof(double));

	for(k=0, m_ChiSqr=0.0; k<m_Data.Get_Count(); k++)
	{
		_Get_Function(m_Data[k].x, Parameters, y, dy_da);

		dy	= m_Data[k].y - y;

		for(i=0; i<m_Params.m_Count; i++)
		{
			for(j=0; j<=i; j++)
			{
				Alpha[i][j]	+= dy_da[i] * dy_da[j];
			}

			Beta[i]		+= dy * dy_da[i];
		}

		m_ChiSqr	+= dy * dy;
	}

	SG_Free(dy_da);

	//-----------------------------------------------------
	for(i=1; i<m_Params.m_Count; i++)
	{
		for(j=0; j<i; j++)
		{
			Alpha[j][i]	= Alpha[i][j];
		}
	}

	return( true );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CSG_Trend::_Get_Function(double x, double *Parameters, double &y, double *dy_da)
{
	int		i;

	//-----------------------------------------------------
	for(i=0; i<m_Params.m_Count; i++)
	{
		m_Formula.Set_Variable(m_Params.m_Variables[i], Parameters[i]);
	}

	y	= m_Formula.Val(x);

	//-----------------------------------------------------
	for(i=0; i<m_Params.m_Count; i++)
	{
		m_Formula.Set_Variable(m_Params.m_Variables[i], Parameters[i] + EPSILON);

		dy_da[i]	 = m_Formula.Val(x);
		dy_da[i]	-= y;
		dy_da[i]	/= EPSILON;

		m_Formula.Set_Variable(m_Params.m_Variables[i], Parameters[i] - EPSILON);
	}
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------