SAGA API v9.10
Loading...
Searching...
No Matches
mat_trend.cpp
Go to the documentation of this file.
1
3// //
4// SAGA //
5// //
6// System for Automated Geoscientific Analyses //
7// //
8// Application Programming Interface //
9// //
10// Library: SAGA_API //
11// //
12//-------------------------------------------------------//
13// //
14// mat_trend.cpp //
15// //
16// Copyright (C) 2006 by //
17// Olaf Conrad //
18// //
19//-------------------------------------------------------//
20// //
21// This file is part of 'SAGA - System for Automated //
22// Geoscientific Analyses'. //
23// //
24// This library is free software; you can redistribute //
25// it and/or modify it under the terms of the GNU Lesser //
26// General Public License as published by the Free //
27// Software Foundation, either version 2.1 of the //
28// License, or (at your option) any later version. //
29// //
30// This library is distributed in the hope that it will //
31// be useful, but WITHOUT ANY WARRANTY; without even the //
32// implied warranty of MERCHANTABILITY or FITNESS FOR A //
33// PARTICULAR PURPOSE. See the GNU Lesser General Public //
34// License for more details. //
35// //
36// You should have received a copy of the GNU Lesser //
37// General Public License along with this program; if //
38// not, see <http://www.gnu.org/licenses/>. //
39// //
40//-------------------------------------------------------//
41// //
42// contact: Olaf Conrad //
43// Institute of Geography //
44// University of Goettingen //
45// Goldschmidtstr. 5 //
46// 37077 Goettingen //
47// Germany //
48// //
49// e-mail: oconrad@saga-gis.org //
50// //
51//-------------------------------------------------------//
52// //
53// Based on the LMFit.h/cpp, Fit.h/cpp source codes of //
54// A. Ringeler (see 'table_calculus' sub project). //
55// //
57
58//---------------------------------------------------------
59#include "mat_tools.h"
60
61
63// //
64// //
65// //
67
68//---------------------------------------------------------
69#define SWAP(a, b) { double temp = (a); (a) = (b); (b) = temp; }
70
71//---------------------------------------------------------
72#define EPSILON 0.001
73
74
76// //
77// //
78// //
80
81//---------------------------------------------------------
82bool CSG_Trend::CParams::Create(const CSG_String &Variables)
83{
84 if( m_Variables.Length() != Variables.Length() )
85 {
86 m_Variables = Variables;
87
88 m_A .Create(Get_Count());
89 m_Atry .Create(Get_Count());
90 m_Beta .Create(Get_Count());
91 m_dA .Create(Get_Count());
92 m_dA2 .Create(Get_Count());
93
94 m_Alpha.Create(Get_Count(), Get_Count());
95 m_Covar.Create(Get_Count(), Get_Count());
96 }
97
98 m_A.Assign(1.);
99
100 return( true );
101}
102
103//---------------------------------------------------------
104bool CSG_Trend::CParams::Destroy(void)
105{
106 m_Variables.Clear();
107
108 m_A .Destroy();
109 m_Atry .Destroy();
110 m_Beta .Destroy();
111 m_dA .Destroy();
112 m_dA2 .Destroy();
113 m_Alpha.Destroy();
114 m_Covar.Destroy();
115
116 m_Alpha.Destroy();
117 m_Covar.Destroy();
118
119 return( true );
120}
121
122
124// //
126
127//---------------------------------------------------------
129{
130 m_Lambda_Max = 10000;
131 m_Iter_Max = 1000;
132 m_bOkay = false;
133
134// Set_Formula("a + b * x");
135}
136
137
139// //
141
142//---------------------------------------------------------
144{
145 m_bOkay = false;
146
147 m_Params.Destroy();
148
149 if( m_Formula.Set_Formula(Formula) )
150 {
151 CSG_String Params, Used(m_Formula.Get_Used_Variables());
152
153 for(size_t i=0; i<Used.Length(); i++)
154 {
155 if( Used[i] >= 'a' && Used[i] <= 'z' && Used[i] != 'x' )
156 {
157 Params += Used[i];
158 }
159 }
160
161 return( m_Params.Create(Params) );
162 }
163
164 return( false );
165}
166
167//---------------------------------------------------------
168bool CSG_Trend::Init_Parameter(const SG_Char &Variable, double Value)
169{
170 for(size_t i=0; i<m_Params.m_Variables.Length(); i++)
171 {
172 if( Variable == m_Params.m_Variables[i] )
173 {
174 m_Params.m_A[(int)i] = Value;
175
176 return( true );
177 }
178 }
179
180 return( false );
181}
182
183
185// //
187
188//---------------------------------------------------------
190{
191 m_Data[0].Create(true);
192 m_Data[1].Create(true);
193
194 m_bOkay = false;
195}
196
197//---------------------------------------------------------
198void CSG_Trend::Set_Data(double *x, double *y, int n, bool bAdd)
199{
200 if( !bAdd )
201 {
202 Clr_Data();
203 }
204
205 for(int i=0; i<n; i++)
206 {
207 Add_Data(x[i], y[i]);
208 }
209}
210
211//---------------------------------------------------------
212void CSG_Trend::Set_Data(const CSG_Points &Data, bool bAdd)
213{
214 if( !bAdd )
215 {
216 Clr_Data();
217 }
218
219 for(sLong i=0; i<Data.Get_Count(); i++)
220 {
221 Add_Data(Data[i].x, Data[i].y);
222 }
223}
224
225//---------------------------------------------------------
226bool CSG_Trend::Add_Data(double x, double y)
227{
228 m_Data[0] += x;
229 m_Data[1] += y;
230
231 m_bOkay = false;
232
233 return( true );
234}
235
236
238// //
240
241//---------------------------------------------------------
243{
244 if( Iterations > 0 )
245 {
246 m_Iter_Max = Iterations;
247
248 return( true );
249 }
250
251 return( false );
252}
253
254//---------------------------------------------------------
255bool CSG_Trend::Set_Max_Lambda(double Lambda)
256{
257 if( Lambda > 0. )
258 {
259 m_Lambda_Max = Lambda;
260
261 return( true );
262 }
263
264 return( false );
265}
266
267
269// //
271
272//---------------------------------------------------------
273bool CSG_Trend::Get_Trend(double *x, double *y, int n, const CSG_String &Formula)
274{
275 Set_Data(x, y, n, false);
276
277 return( !Formula.is_Empty() && !Set_Formula(Formula) ? false : Get_Trend() );
278}
279
280//---------------------------------------------------------
281bool CSG_Trend::Get_Trend(const CSG_Points &Data, const CSG_String &Formula)
282{
283 Set_Data(Data, false);
284
285 return( !Formula.is_Empty() && !Set_Formula(Formula) ? false : Get_Trend() );
286}
287
288//---------------------------------------------------------
290{
291 return( !Formula.is_Empty() && !Set_Formula(Formula) ? false : Get_Trend() );
292}
293
294//---------------------------------------------------------
296{
297 CSG_String sError;
298
299 if( m_Formula.Get_Error(sError) )
300 {
301 return( false );
302 }
303
304 if( Get_Data_Count() <= 1 )
305 {
306 return( false );
307 }
308
309 //-----------------------------------------------------
310 m_bOkay = true;
311
312 int i;
313
314 if( m_Params.Get_Count() > 0 )
315 {
316 m_Lambda = 0.001;
317
318 _Get_mrqcof(m_Params.m_A, m_Params.m_Alpha, m_Params.m_Beta);
319
320 m_ChiSqr_o = m_ChiSqr;
321
322 for(i=0; i<m_Params.Get_Count(); i++)
323 {
324 m_Params.m_Atry[i] = m_Params.m_A[i];
325 }
326
327 //-------------------------------------------------
328 for(i=0; i<m_Iter_Max && m_Lambda<m_Lambda_Max && m_bOkay && SG_UI_Process_Get_Okay(false); i++)
329 {
330 m_bOkay = _Fit_Function();
331 }
332
333 //-------------------------------------------------
334 for(i=0; i<m_Params.Get_Count(); i++)
335 {
336 m_Formula.Set_Variable((char)m_Params.m_Variables[i], m_Params.m_A[i]);
337 }
338 }
339
340 //-----------------------------------------------------
341 double y_m, y_o, y_t;
342
343 for(i=0, y_m=0.; i<Get_Data_Count(); i++)
344 {
345 y_m += Get_Data_Y(i);
346 }
347
348 y_m /= Get_Data_Count();
349
350 for(i=0, y_o=0., y_t=0.; i<Get_Data_Count(); i++)
351 {
352 y_o += SG_Get_Square(y_m - Get_Data_Y(i) );
353 y_t += SG_Get_Square(y_m - m_Formula.Get_Value(Get_Data_X(i)));
354 }
355
356 m_ChiSqr_o = y_o > 0. ? y_t / y_o : 0.;
357
358 return( m_bOkay );
359}
360
361
363// //
365
366//---------------------------------------------------------
368{
369 CSG_String Message;
370
371 if( !m_bOkay )
372 {
373 if( !m_Formula.Get_Error(Message) )
374 {
375 Message.Printf(_TL("Error in Trend Calculation"));
376 }
377 }
378
379 return( Message );
380}
381
382//---------------------------------------------------------
384{
385 CSG_String s;
386
387 switch( Type )
388 {
389 case SG_TREND_STRING_Formula: default:
390 s += m_Formula.Get_Formula();
391 break;
392
394 s += m_Formula.Get_Formula();
395 s += "\n";
396
397 if( m_Params.Get_Count() > 0 )
398 {
399 s += "\n";
400
401 for(int i=0; i<m_Params.Get_Count() && m_bOkay; i++)
402 {
403 s += CSG_String::Format("%c = %g\n", m_Params.m_Variables[i], m_Params.m_A[i]);
404 }
405 }
406 break;
407
409 s += m_Formula.Get_Formula();
410 s += "\n";
411
412 if( m_Params.Get_Count() > 0 && m_bOkay )
413 {
414 s += "\n";
415
416 for(int i=0; i<m_Params.Get_Count(); i++)
417 {
418 s += CSG_String::Format("%c = %g\n", m_Params.m_Variables[i], m_Params.m_A[i]);
419 }
420 }
421 break;
422
424 s += m_Formula.Get_Formula();
425 s += "\n";
426
427 if( m_Params.Get_Count() > 0 && m_bOkay )
428 {
429 s += "\n";
430
431 for(int i=0; i<m_Params.Get_Count(); i++)
432 {
433 s += CSG_String::Format("%c = %g\n", m_Params.m_Variables[i], m_Params.m_A[i]);
434 }
435 }
436
437 s += "\n";
438 s += CSG_String::Format("N = %d\n" , Get_Data_Count());
439 s += CSG_String::Format("R2 = %g\n", Get_R2() * 100.);
440 break;
441
443 s += m_Formula.Get_Formula();
444
445 if( m_Params.Get_Count() > 0 && m_bOkay )
446 {
447 for(int i=0; i<m_Params.Get_Count(); i++)
448 {
449 s += CSG_String::Format("%s%c=%g", !i ? SG_T("; (") : SG_T(", "), m_Params.m_Variables[i], m_Params.m_A[i]);
450 }
451
452 s += ")";
453 }
454
455 s += CSG_String::Format("; N=%d" , Get_Data_Count());
456 s += CSG_String::Format("; R2=%.2f%%", Get_R2() * 100.);
457 break;
458 }
459
460 return( s );
461}
462
463
465// //
467
468//---------------------------------------------------------
469bool CSG_Trend::_Fit_Function(void)
470{
471 int i, j;
472
473 //-----------------------------------------------------
474 for(i=0; i<m_Params.Get_Count(); i++)
475 {
476 for(j=0; j<m_Params.Get_Count(); j++)
477 {
478 m_Params.m_Covar[i][j] = m_Params.m_Alpha[i][j];
479 }
480
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];
483 }
484
485 //-----------------------------------------------------
486 if( _Get_Gaussj() == false )
487 {
488 return( false );
489 }
490
491 //-----------------------------------------------------
492 for(i=0; i<m_Params.Get_Count(); i++)
493 {
494 m_Params.m_dA[i] = m_Params.m_dA2[i];
495 }
496
497 //-----------------------------------------------------
498 if( m_Lambda == 0. )
499 {
500 for(i=m_Params.Get_Count()-1; i>0; i--)
501 {
502 for(j=0; j<m_Params.Get_Count(); j++)
503 {
504 SWAP(m_Params.m_Covar[j][i], m_Params.m_Covar[j][i-1]);
505 }
506
507 for(j=0; j<m_Params.Get_Count(); j++)
508 {
509 SWAP(m_Params.m_Covar[i][j], m_Params.m_Covar[i-1][j]);
510 }
511 }
512 }
513 else
514 {
515 for(i=0; i<m_Params.Get_Count(); i++)
516 {
517 m_Params.m_Atry[i] = m_Params.m_A[i] + m_Params.m_dA[i];
518 }
519
520 _Get_mrqcof(m_Params.m_Atry, m_Params.m_Covar, m_Params.m_dA);
521
522 if( m_ChiSqr < m_ChiSqr_o )
523 {
524 m_Lambda *= 0.1;
525 m_ChiSqr_o = m_ChiSqr;
526
527 for(i=0; i<m_Params.Get_Count(); i++)
528 {
529 for(j=0; j<m_Params.Get_Count(); j++)
530 {
531 m_Params.m_Alpha[i][j] = m_Params.m_Covar[i][j];
532 }
533
534 m_Params.m_Beta[i] = m_Params.m_dA[i];
535 }
536
537 for(i=0; i<m_Params.Get_Count(); i++) // Achtung!! in aelteren Versionen war hier ein Fehler
538 {
539 m_Params.m_A[i] = m_Params.m_Atry[i];
540 }
541 }
542 else
543 {
544 m_Lambda *= 10.;
545 m_ChiSqr = m_ChiSqr_o;
546 }
547 }
548
549 //-----------------------------------------------------
550 return( true );
551}
552
553//---------------------------------------------------------
554bool CSG_Trend::_Get_Gaussj(void)
555{
556 int i, j, k, iCol, iRow;
557
558 //-----------------------------------------------------
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());
562
563 for(i=0; i<m_Params.Get_Count(); i++)
564 {
565 ipiv[i] = 0;
566 }
567
568 //-----------------------------------------------------
569 for(i=0, iCol=-1, iRow=-1; i<m_Params.Get_Count(); i++)
570 {
571 double big = 0.;
572
573 for(j=0; j<m_Params.Get_Count(); j++)
574 {
575 if( ipiv[j] != 1 )
576 {
577 for(k=0; k<m_Params.Get_Count(); k++)
578 {
579 if( ipiv[k] == 0 )
580 {
581 if( fabs(m_Params.m_Covar[j][k]) >= big )
582 {
583 big = fabs(m_Params.m_Covar[j][k]);
584 iRow = j;
585 iCol = k;
586 }
587 }
588 else if( ipiv[k] > 1 )
589 {
590 return( false ); // singular matrix...
591 }
592 }
593 }
594 }
595
596 if( iCol < 0 || iRow < 0 )
597 {
598 return( false ); // singular matrix...
599 }
600
601 //-------------------------------------------------
602 ipiv[iCol]++;
603
604 if( iRow != iCol )
605 {
606 for(j=0; j<m_Params.Get_Count(); j++)
607 {
608 SWAP(m_Params.m_Covar[iRow][j], m_Params.m_Covar[iCol][j]);
609 }
610
611 SWAP(m_Params.m_dA2[iRow], m_Params.m_dA2[iCol]);
612 }
613
614 indxr[i] = iRow;
615 indxc[i] = iCol;
616
617 if( fabs(m_Params.m_Covar[iCol][iCol]) < 1E-300 )
618 {
619 return( false ); // singular matrix...
620 }
621
622 //-------------------------------------------------
623 double pivinv = 1. / m_Params.m_Covar[iCol][iCol];
624
625 m_Params.m_Covar[iCol][iCol] = 1.;
626
627 for(j=0; j<m_Params.Get_Count(); j++)
628 {
629 m_Params.m_Covar[iCol][j] *= pivinv;
630 }
631
632 m_Params.m_dA2[iCol] *= pivinv;
633
634 //-------------------------------------------------
635 for(j=0; j<m_Params.Get_Count(); j++)
636 {
637 if( j != iCol )
638 {
639 double temp = m_Params.m_Covar[j][iCol];
640
641 m_Params.m_Covar[j][iCol] = 0.;
642
643 for(k=0; k<m_Params.Get_Count(); k++)
644 {
645 m_Params.m_Covar[j][k] -= m_Params.m_Covar[iCol][k] * temp;
646 }
647
648 m_Params.m_dA2[j] -= m_Params.m_dA2[iCol] * temp;
649 }
650 }
651 }
652
653 //-----------------------------------------------------
654 for(i=m_Params.Get_Count()-1; i>=0; i--)
655 {
656 if( indxr[i] != indxc[i] )
657 {
658 for(j=0; j<m_Params.Get_Count(); j++)
659 {
660 SWAP(m_Params.m_Covar[j][indxr[i]], m_Params.m_Covar[j][indxc[i]]);
661 }
662 }
663 }
664
665 //-----------------------------------------------------
666 return( true );
667}
668
669//---------------------------------------------------------
670bool CSG_Trend::_Get_mrqcof(CSG_Vector &Parameters, CSG_Matrix &Alpha, CSG_Vector &Beta)
671{
672 CSG_Vector dy_da(m_Params.Get_Count());
673
674 Alpha.Assign(0.);
675 Beta .Assign(0.);
676
677 m_ChiSqr = 0.;
678
679 for(int k=0; k<Get_Data_Count(); k++)
680 {
681 double y;
682
683 _Get_Function(y, dy_da.Get_Data(), Get_Data_X(k), Parameters);
684
685 double dy = Get_Data_Y(k) - y;
686
687 for(int i=0; i<m_Params.Get_Count(); i++)
688 {
689 for(int j=0; j<=i; j++)
690 {
691 Alpha[i][j] += dy_da[i] * dy_da[j];
692 }
693
694 Beta[i] += dy * dy_da[i];
695 }
696
697 m_ChiSqr += dy * dy;
698 }
699
700 //-----------------------------------------------------
701 for(int i=1; i<m_Params.Get_Count(); i++)
702 {
703 for(int j=0; j<i; j++)
704 {
705 Alpha[j][i] = Alpha[i][j];
706 }
707 }
708
709 return( true );
710}
711
712
714// //
716
717//---------------------------------------------------------
718void CSG_Trend::_Get_Function(double &y, double *dy_da, double x, const double *Parameters)
719{
720 for(int i=0; i<m_Params.Get_Count(); i++)
721 {
722 m_Formula.Set_Variable((char)m_Params.m_Variables[i], Parameters[i]);
723 }
724
725 y = m_Formula.Get_Value(x);
726
727 //-----------------------------------------------------
728 for(int i=0; i<m_Params.Get_Count(); i++)
729 {
730 m_Formula.Set_Variable((char)m_Params.m_Variables[i], Parameters[i] + EPSILON);
731
732 dy_da[i] = m_Formula.Get_Value(x);
733 dy_da[i] -= y;
734 dy_da[i] /= EPSILON;
735
736 m_Formula.Set_Variable((char)m_Params.m_Variables[i], Parameters[i] - EPSILON);
737 }
738}
739
740
742// //
743// //
744// //
746
747//---------------------------------------------------------
752
753//---------------------------------------------------------
755{
756 m_Order = 0;
757
758 Clr_Data();
759
760 return( true );
761}
762
763//---------------------------------------------------------
765{
766 m_a.Destroy();
767
768 if( Order > 0 )
769 {
770 m_Order = Order;
771
772 return( true );
773 }
774
775 return( false );
776}
777
778//---------------------------------------------------------
780{
781 m_a.Destroy();
782 m_y.Destroy();
783 m_x.Destroy();
784
785 return( true );
786}
787
788//---------------------------------------------------------
789bool CSG_Trend_Polynom::Set_Data(double *x, double *y, int n, bool bAdd)
790{
791 if( !bAdd )
792 {
793 Clr_Data();
794 }
795
796 m_x.Add_Rows(n);
797 m_y.Add_Rows(n);
798
799 for(int i=0, j=m_x.Get_N()-1; i<n; i++)
800 {
801 m_x[j] = x[i];
802 m_y[j] = y[i];
803 }
804
805 return( true );
806}
807
808//---------------------------------------------------------
809bool CSG_Trend_Polynom::Add_Data(double x, double y)
810{
811 return( m_x.Add_Row(x) && m_y.Add_Row(y) );
812}
813
814//---------------------------------------------------------
816{
817 if( m_Order < 1 || m_x.Get_N() <= m_Order )
818 {
819 return( false );
820 }
821
822 //-----------------------------------------------------
823 int i;
824 double d, Ym, SSE, SSR;
825 CSG_Matrix X, Xt, C;
826
827 //-----------------------------------------------------
828 X .Create((sLong)m_Order + 1, m_y.Get_Size());
829 Xt.Create(m_y.Get_Size(), (uLong)m_Order + 1);
830
831 for(i=0, Ym=0.; i<m_y.Get_N(); i++)
832 {
833 X[i][0] = Xt[0][i] = d = 1.;
834
835 for(int j=1; j<=m_Order; j++)
836 {
837 X[i][j] = Xt[j][i] = (d = d * m_x[i]);
838 }
839
840 Ym += m_y[i];
841 }
842
843 Ym /= m_y.Get_N();
844
845 m_a = (Xt * X).Get_Inverse() * (Xt * m_y);
846
847 //-----------------------------------------------------
848 CSG_Vector Yr = X * m_a;
849
850 for(i=0, SSE=0., SSR=0.; i<m_y.Get_N(); i++)
851 {
852 SSE += SG_Get_Square(Yr[i] - m_y[i]);
853 SSR += SG_Get_Square(Yr[i] - Ym );
854 }
855
856 m_r2 = SSR / (SSR + SSE);
857
858 return( true );
859}
860
861//---------------------------------------------------------
862double CSG_Trend_Polynom::Get_Value(double x) const
863{
864 if( m_a.Get_N() > 0 )
865 {
866 double y = m_a(0);
867 double d = 1.;
868
869 for(int i=1; i<m_a.Get_N(); i++)
870 {
871 d *= x;
872 y += d * m_a(i);
873 }
874
875 return( y );
876 }
877
878 return( 0. );
879}
880
881
883// //
884// //
885// //
887
888//---------------------------------------------------------
bool SG_UI_Process_Get_Okay(bool bBlink)
unsigned long long uLong
Definition api_core.h:159
signed long long sLong
Definition api_core.h:158
#define SG_T(s)
Definition api_core.h:537
#define SG_Char
Definition api_core.h:536
#define _TL(s)
Definition api_core.h:1568
bool Create(const CSG_Matrix &Matrix)
bool Assign(double Scalar)
sLong Get_Count(void) const
Definition geo_tools.h:201
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)
bool Get_Trend(void)
double Get_Value(double x) const
bool Destroy(void)
bool Clr_Data(void)
bool Init_Parameter(const SG_Char &Variable, double Value)
double Get_Data_X(int i) const
Definition mat_tools.h:2008
CSG_Trend(void)
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
Definition mat_tools.h:2007
double Get_R2(void) const
Definition mat_tools.h:2032
void Clr_Data(void)
bool Set_Max_Iterations(int Iterations)
bool Get_Trend(void)
double Get_Data_Y(int i) const
Definition mat_tools.h:2012
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)
#define C
#define EPSILON
@ SG_TREND_STRING_Complete
Definition mat_tools.h:1986
@ SG_TREND_STRING_Formula_Parameters
Definition mat_tools.h:1985
@ SG_TREND_STRING_Formula
Definition mat_tools.h:1983
@ SG_TREND_STRING_Function
Definition mat_tools.h:1984
@ SG_TREND_STRING_Compact
Definition mat_tools.h:1987
SAGA_API_DLL_EXPORT double SG_Get_Square(double Value)
Definition mat_tools.cpp:70
#define SWAP(a, b)
Definition mat_trend.cpp:69