📄 mat_trend.cpp
字号:
///////////////////////////////////////////////////////////
// //
// SAGA //
// //
// System for Automated Geoscientific Analyses //
// //
// Application Programming Interface //
// //
// Library: SAGA_API //
// //
//-------------------------------------------------------//
// //
// mat_trend.cpp //
// //
// Copyright (C) 2006 by //
// Olaf Conrad //
// //
//-------------------------------------------------------//
// //
// This file is part of 'SAGA - System for Automated //
// Geoscientific Analyses'. SAGA is free software; you //
// can redistribute it and/or modify it under the terms //
// of the GNU General Public License as published by the //
// Free Software Foundation; version 2 of the License. //
// //
// SAGA is distributed in the hope that it will be //
// useful, but WITHOUT ANY WARRANTY; without even the //
// implied warranty of MERCHANTABILITY or FITNESS FOR A //
// PARTICULAR PURPOSE. See the GNU General Public //
// License for more details. //
// //
// You should have received a copy of the GNU General //
// Public License along with this program; if not, //
// write to the Free Software Foundation, Inc., //
// 59 Temple Place - Suite 330, Boston, MA 02111-1307, //
// USA. //
// //
//-------------------------------------------------------//
// //
// contact: Olaf Conrad //
// Institute of Geography //
// University of Goettingen //
// Goldschmidtstr. 5 //
// 37077 Goettingen //
// Germany //
// //
// e-mail: oconrad@saga-gis.org //
// //
//-------------------------------------------------------//
// //
// Based on the LMFit.h/cpp, Fit.h/cpp source codes of //
// A. Ringeler (see 'table_calculus' sub project). //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
#include "mat_tools.h"
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
#define SWAP(a, b) { temp = (a); (a) = (b); (b) = temp; }
//---------------------------------------------------------
#define EPSILON 0.001
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
CSG_Trend::CFncParams::CFncParams(void)
{
m_Count = 0;
}
//---------------------------------------------------------
CSG_Trend::CFncParams::~CFncParams(void)
{
Destroy();
}
//---------------------------------------------------------
bool CSG_Trend::CFncParams::Create(const SG_Char *Variables, int nVariables)
{
int i;
if( m_Count != nVariables )
{
Destroy();
m_Count = nVariables;
m_Variables = (SG_Char *)SG_Calloc(m_Count, sizeof(SG_Char));
m_A = (double *)SG_Calloc(m_Count, sizeof(double));
m_Atry = (double *)SG_Calloc(m_Count, sizeof(double));
m_Beta = (double *)SG_Calloc(m_Count, sizeof(double));
m_dA = (double *)SG_Calloc(m_Count, sizeof(double));
m_dA2 = (double *)SG_Calloc(m_Count, sizeof(double));
m_Alpha = (double **)SG_Calloc(m_Count, sizeof(double *));
m_Covar = (double **)SG_Calloc(m_Count, sizeof(double *));
for(i=0; i<m_Count; i++)
{
m_Alpha[i] = (double *)SG_Calloc(m_Count, sizeof(double));
m_Covar[i] = (double *)SG_Calloc(m_Count, sizeof(double));
}
}
for(i=0; i<m_Count; i++)
{
m_Variables[i] = Variables[i];
m_A[i] = 1.0;
}
return( true );
}
//---------------------------------------------------------
bool CSG_Trend::CFncParams::Destroy(void)
{
if( m_Count > 0 )
{
for(int i=0; i<m_Count; i++)
{
SG_Free(m_Alpha[i]);
SG_Free(m_Covar[i]);
}
SG_Free(m_Variables);
SG_Free(m_A);
SG_Free(m_Atry);
SG_Free(m_Beta);
SG_Free(m_dA);
SG_Free(m_dA2);
SG_Free(m_Alpha);
SG_Free(m_Covar);
m_Count = 0;
}
return( true );
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
CSG_Trend::CSG_Trend(void)
{
m_Lambda_Max = 10000;
m_Iter_Max = 1000;
// Set_Formula("a * x + b");
}
//---------------------------------------------------------
CSG_Trend::~CSG_Trend(void)
{
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
bool CSG_Trend::Set_Formula(const SG_Char *Formula)
{
m_bOkay = false;
if( m_Formula.Set_Formula(Formula) )
{
CSG_String vars, uvars(m_Formula.Get_Used_Var());
for(unsigned int i=0; i<uvars.Length(); i++)
{
if( uvars.c_str()[i] >= 'a' && uvars.c_str()[i] <= 'z' && uvars.c_str()[i] != 'x' )
{
vars.Append(uvars.c_str()[i]);
}
}
return( m_Params.Create(vars.c_str(), vars.Length()) );
}
m_Params.Destroy();
return( false );
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
void CSG_Trend::Clr_Data(void)
{
m_Data.Clear();
m_bOkay = false;
}
//---------------------------------------------------------
void CSG_Trend::Set_Data(double *xData, double *yData, int nData, bool bAdd)
{
if( !bAdd )
{
m_Data.Clear();
}
for(int i=0; i<nData; i++)
{
m_Data.Add(xData[i], yData[i]);
}
m_bOkay = false;
}
//---------------------------------------------------------
void CSG_Trend::Set_Data(const CSG_Points &Data, bool bAdd)
{
if( !bAdd )
{
m_Data.Clear();
}
for(int i=0; i<Data.Get_Count(); i++)
{
m_Data.Add(Data.Get_X(i), Data.Get_Y(i));
}
m_bOkay = false;
}
//---------------------------------------------------------
void CSG_Trend::Add_Data(double x, double y)
{
m_Data.Add(x, y);
m_bOkay = false;
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
bool CSG_Trend::Set_Max_Iterations(int Iterations)
{
if( Iterations > 0 )
{
m_Iter_Max = Iterations;
return( true );
}
return( false );
}
//---------------------------------------------------------
bool CSG_Trend::Set_Max_Lambda(double Lambda)
{
if( Lambda > 0.0 )
{
m_Lambda_Max = Lambda;
return( true );
}
return( false );
}
///////////////////////////////////////////////////////////
// //
// //
// //
///////////////////////////////////////////////////////////
//---------------------------------------------------------
bool CSG_Trend::Get_Trend(double *xData, double *yData, int nData, const SG_Char *Formula)
{
Set_Data(xData, yData, false);
if( Formula )
{
Set_Formula(Formula);
}
return( Get_Trend() );
}
//---------------------------------------------------------
bool CSG_Trend::Get_Trend(const CSG_Points &Data, const SG_Char *Formula)
{
Set_Data(Data, false);
if( Formula )
{
Set_Formula(Formula);
}
return( Get_Trend() );
}
//---------------------------------------------------------
bool CSG_Trend::Get_Trend(void)
{
if( !m_Formula.Get_Error(NULL, NULL) )
{
int i;
m_bOkay = true;
//-------------------------------------------------
if( m_Data.Get_Count() > 1 )
{
if( m_Params.m_Count > 0 )
{
m_Lambda = 0.001;
_Get_mrqcof(m_Params.m_A, m_Params.m_Alpha, m_Params.m_Beta);
m_ChiSqr_o = m_ChiSqr;
for(i=0; i<m_Params.m_Count; i++)
{
m_Params.m_Atry[i] = m_Params.m_A[i];
}
//-----------------------------------------
for(i=0; i<m_Iter_Max && m_Lambda<m_Lambda_Max && m_bOkay && SG_UI_Process_Get_Okay(false); i++)
{
m_bOkay = _Fit_Function();
}
//-----------------------------------------
for(i=0; i<m_Params.m_Count; i++)
{
m_Formula.Set_Variable(m_Params.m_Variables[i], m_Params.m_A[i]);
}
}
//---------------------------------------------
double y_m, y_o, y_t;
for(i=0, y_m=0.0; i<m_Data.Get_Count(); i++)
{
y_m += m_Data.Get_Y(i);
}
y_m /= m_Data.Get_Count();
for(i=0, y_o=0.0, y_t=0.0; i<m_Data.Get_Count(); i++)
{
y_o += SG_Get_Square(y_m - m_Data.Get_Y(i));
y_t += SG_Get_Square(y_m - m_Formula.Val(m_Data.Get_X(i)));
}
m_ChiSqr_o = y_o > 0.0 ? y_t / y_o : 1.0;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -