69#define SWAP(a, b) { double temp = (a); (a) = (b); (b) = temp; }
82bool CSG_Trend::CParams::Create(
const CSG_String &Variables)
86 m_Variables = Variables;
89 m_Atry .
Create(Get_Count());
90 m_Beta .
Create(Get_Count());
92 m_dA2 .
Create(Get_Count());
94 m_Alpha.
Create(Get_Count(), Get_Count());
95 m_Covar.
Create(Get_Count(), Get_Count());
104bool CSG_Trend::CParams::Destroy(
void)
130 m_Lambda_Max = 10000;
149 if( m_Formula.Set_Formula(Formula) )
151 CSG_String Params, Used(m_Formula.Get_Used_Variables());
153 for(
size_t i=0; i<Used.
Length(); i++)
155 if( Used[i] >=
'a' && Used[i] <=
'z' && Used[i] !=
'x' )
161 return( m_Params.Create(Params) );
170 for(
size_t i=0; i<m_Params.m_Variables.Length(); i++)
172 if( Variable == m_Params.m_Variables[i] )
174 m_Params.m_A[(int)i] = Value;
191 m_Data[0].Create(
true);
192 m_Data[1].Create(
true);
205 for(
int i=0; i<n; i++)
246 m_Iter_Max = Iterations;
259 m_Lambda_Max = Lambda;
299 if( m_Formula.Get_Error(sError) )
314 if( m_Params.Get_Count() > 0 )
318 _Get_mrqcof(m_Params.m_A, m_Params.m_Alpha, m_Params.m_Beta);
320 m_ChiSqr_o = m_ChiSqr;
322 for(i=0; i<m_Params.Get_Count(); i++)
324 m_Params.m_Atry[i] = m_Params.m_A[i];
330 m_bOkay = _Fit_Function();
334 for(i=0; i<m_Params.Get_Count(); i++)
336 m_Formula.Set_Variable((
char)m_Params.m_Variables[i], m_Params.m_A[i]);
341 double y_m, y_o, y_t;
356 m_ChiSqr_o = y_o > 0. ? y_t / y_o : 0.;
373 if( !m_Formula.Get_Error(Message) )
375 Message.
Printf(
_TL(
"Error in Trend Calculation"));
390 s += m_Formula.Get_Formula();
394 s += m_Formula.Get_Formula();
397 if( m_Params.Get_Count() > 0 )
401 for(
int i=0; i<m_Params.Get_Count() && m_bOkay; i++)
409 s += m_Formula.Get_Formula();
412 if( m_Params.Get_Count() > 0 && m_bOkay )
416 for(
int i=0; i<m_Params.Get_Count(); i++)
424 s += m_Formula.Get_Formula();
427 if( m_Params.Get_Count() > 0 && m_bOkay )
431 for(
int i=0; i<m_Params.Get_Count(); i++)
443 s += m_Formula.Get_Formula();
445 if( m_Params.Get_Count() > 0 && m_bOkay )
447 for(
int i=0; i<m_Params.Get_Count(); i++)
469bool CSG_Trend::_Fit_Function(
void)
474 for(i=0; i<m_Params.Get_Count(); i++)
476 for(j=0; j<m_Params.Get_Count(); j++)
478 m_Params.m_Covar[i][j] = m_Params.m_Alpha[i][j];
481 m_Params.m_Covar[i][i] = m_Params.m_Alpha[i][i] * (1. + m_Lambda);
482 m_Params.m_dA2 [i] = m_Params.m_Beta [i];
486 if( _Get_Gaussj() ==
false )
492 for(i=0; i<m_Params.Get_Count(); i++)
494 m_Params.m_dA[i] = m_Params.m_dA2[i];
500 for(i=m_Params.Get_Count()-1; i>0; i--)
502 for(j=0; j<m_Params.Get_Count(); j++)
504 SWAP(m_Params.m_Covar[j][i], m_Params.m_Covar[j][i-1]);
507 for(j=0; j<m_Params.Get_Count(); j++)
509 SWAP(m_Params.m_Covar[i][j], m_Params.m_Covar[i-1][j]);
515 for(i=0; i<m_Params.Get_Count(); i++)
517 m_Params.m_Atry[i] = m_Params.m_A[i] + m_Params.m_dA[i];
520 _Get_mrqcof(m_Params.m_Atry, m_Params.m_Covar, m_Params.m_dA);
522 if( m_ChiSqr < m_ChiSqr_o )
525 m_ChiSqr_o = m_ChiSqr;
527 for(i=0; i<m_Params.Get_Count(); i++)
529 for(j=0; j<m_Params.Get_Count(); j++)
531 m_Params.m_Alpha[i][j] = m_Params.m_Covar[i][j];
534 m_Params.m_Beta[i] = m_Params.m_dA[i];
537 for(i=0; i<m_Params.Get_Count(); i++)
539 m_Params.m_A[i] = m_Params.m_Atry[i];
545 m_ChiSqr = m_ChiSqr_o;
554bool CSG_Trend::_Get_Gaussj(
void)
556 int i, j, k, iCol, iRow;
559 CSG_Array_Int indxc(m_Params.Get_Count());
560 CSG_Array_Int indxr(m_Params.Get_Count());
561 CSG_Array_Int ipiv (m_Params.Get_Count());
563 for(i=0; i<m_Params.Get_Count(); i++)
569 for(i=0, iCol=-1, iRow=-1; i<m_Params.Get_Count(); i++)
573 for(j=0; j<m_Params.Get_Count(); j++)
577 for(k=0; k<m_Params.Get_Count(); k++)
581 if( fabs(m_Params.m_Covar[j][k]) >= big )
583 big = fabs(m_Params.m_Covar[j][k]);
588 else if( ipiv[k] > 1 )
596 if( iCol < 0 || iRow < 0 )
606 for(j=0; j<m_Params.Get_Count(); j++)
608 SWAP(m_Params.m_Covar[iRow][j], m_Params.m_Covar[iCol][j]);
611 SWAP(m_Params.m_dA2[iRow], m_Params.m_dA2[iCol]);
617 if( fabs(m_Params.m_Covar[iCol][iCol]) < 1E-300 )
623 double pivinv = 1. / m_Params.m_Covar[iCol][iCol];
625 m_Params.m_Covar[iCol][iCol] = 1.;
627 for(j=0; j<m_Params.Get_Count(); j++)
629 m_Params.m_Covar[iCol][j] *= pivinv;
632 m_Params.m_dA2[iCol] *= pivinv;
635 for(j=0; j<m_Params.Get_Count(); j++)
639 double temp = m_Params.m_Covar[j][iCol];
641 m_Params.m_Covar[j][iCol] = 0.;
643 for(k=0; k<m_Params.Get_Count(); k++)
645 m_Params.m_Covar[j][k] -= m_Params.m_Covar[iCol][k] * temp;
648 m_Params.m_dA2[j] -= m_Params.m_dA2[iCol] * temp;
654 for(i=m_Params.Get_Count()-1; i>=0; i--)
656 if( indxr[i] != indxc[i] )
658 for(j=0; j<m_Params.Get_Count(); j++)
660 SWAP(m_Params.m_Covar[j][indxr[i]], m_Params.m_Covar[j][indxc[i]]);
670bool CSG_Trend::_Get_mrqcof(CSG_Vector &Parameters, CSG_Matrix &Alpha, CSG_Vector &Beta)
672 CSG_Vector dy_da(m_Params.Get_Count());
683 _Get_Function(y, dy_da.Get_Data(),
Get_Data_X(k), Parameters);
687 for(
int i=0; i<m_Params.Get_Count(); i++)
689 for(
int j=0; j<=i; j++)
691 Alpha[i][j] += dy_da[i] * dy_da[j];
694 Beta[i] += dy * dy_da[i];
701 for(
int i=1; i<m_Params.Get_Count(); i++)
703 for(
int j=0; j<i; j++)
705 Alpha[j][i] = Alpha[i][j];
718void CSG_Trend::_Get_Function(
double &y,
double *dy_da,
double x,
const double *Parameters)
720 for(
int i=0; i<m_Params.Get_Count(); i++)
722 m_Formula.Set_Variable((
char)m_Params.m_Variables[i], Parameters[i]);
725 y = m_Formula.Get_Value(x);
728 for(
int i=0; i<m_Params.Get_Count(); i++)
730 m_Formula.Set_Variable((
char)m_Params.m_Variables[i], Parameters[i] +
EPSILON);
732 dy_da[i] = m_Formula.Get_Value(x);
736 m_Formula.Set_Variable((
char)m_Params.m_Variables[i], Parameters[i] -
EPSILON);
799 for(
int i=0, j=m_x.Get_N()-1; i<n; i++)
811 return( m_x.Add_Row(x) && m_y.Add_Row(y) );
817 if( m_Order < 1 || m_x.Get_N() <= m_Order )
824 double d, Ym, SSE, SSR;
831 for(i=0, Ym=0.; i<m_y.Get_N(); i++)
833 X[i][0] = Xt[0][i] = d = 1.;
835 for(
int j=1; j<=m_Order; j++)
837 X[i][j] = Xt[j][i] = (d = d * m_x[i]);
845 m_a = (Xt * X).Get_Inverse() * (Xt * m_y);
850 for(i=0, SSE=0., SSR=0.; i<m_y.Get_N(); i++)
856 m_r2 = SSR / (SSR + SSE);
864 if( m_a.Get_N() > 0 )
869 for(
int i=1; i<m_a.Get_N(); i++)
bool SG_UI_Process_Get_Okay(bool bBlink)
bool Create(const CSG_Matrix &Matrix)
bool Assign(double Scalar)
sLong Get_Count(void) const
size_t Length(void) const
static CSG_String Format(const char *Format,...)
int Printf(const char *Format,...)
bool is_Empty(void) const
bool Add_Data(double x, double y)
bool Set_Data(double *x, double *y, int n, bool bAdd=false)
bool Set_Order(int Order=1)
double Get_Value(double x) const
bool Init_Parameter(const SG_Char &Variable, double Value)
double Get_Data_X(int i) const
void Set_Data(double *x, double *y, int n, bool bAdd=false)
CSG_String Get_Error(void)
bool Set_Formula(const CSG_String &Formula)
int Get_Data_Count(void) const
double Get_R2(void) const
bool Set_Max_Iterations(int Iterations)
double Get_Data_Y(int i) const
bool Add_Data(double x, double y)
bool Set_Max_Lambda(double Lambda)
CSG_String Get_Formula(int Type=SG_TREND_STRING_Complete)
bool Create(const CSG_Vector &Vector)
bool Assign(double Scalar)