SAGA API Version 9.12
Loading...
Searching...
No Matches
mat_tools.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_tools.cpp //
15// //
16// Copyright (C) 2005 by Olaf Conrad //
17// //
18//-------------------------------------------------------//
19// //
20// This file is part of 'SAGA - System for Automated //
21// Geoscientific Analyses'. //
22// //
23// This library is free software; you can redistribute //
24// it and/or modify it under the terms of the GNU Lesser //
25// General Public License as published by the Free //
26// Software Foundation, either version 2.1 of the //
27// License, or (at your option) any later version. //
28// //
29// This library is distributed in the hope that it will //
30// be useful, but WITHOUT ANY WARRANTY; without even the //
31// implied warranty of MERCHANTABILITY or FITNESS FOR A //
32// PARTICULAR PURPOSE. See the GNU Lesser General Public //
33// License for more details. //
34// //
35// You should have received a copy of the GNU Lesser //
36// General Public License along with this program; if //
37// not, see <http://www.gnu.org/licenses/>. //
38// //
39//-------------------------------------------------------//
40// //
41// contact: Olaf Conrad //
42// Institute of Geography //
43// University of Goettingen //
44// Goldschmidtstr. 5 //
45// 37077 Goettingen //
46// Germany //
47// //
48// e-mail: oconrad@saga-gis.org //
49// //
51
52//---------------------------------------------------------
53#include <time.h>
54#include <cfloat>
55
56#include "mat_tools.h"
57
58#include "table.h"
59#include "grid.h"
60#include "grids.h"
61
62
64// //
65// //
66// //
68
69//---------------------------------------------------------
70double SG_Get_Square(double Value)
71{
72 return( Value * Value );
73}
74
75
77// //
78// //
79// //
81
82//---------------------------------------------------------
83double SG_Get_Rounded(double Value, int Decimals)
84{
85 if( Decimals < 0 )
86 {
87 return( Value );
88 }
89
90 if( Decimals == 0 )
91 {
92 return( floor(0.5 + Value) );
93 }
94
95 double d = pow(10., Decimals);
96 double v = Value * d;
97
98 if( fabs(v - floor(v)) > 0. )
99 {
100 return( floor(0.5 + v) / d );
101 }
102
103 return( Value );
104}
105
106//---------------------------------------------------------
107double SG_Get_Rounded_To_SignificantFigures(double Value, int Decimals)
108{
109 if( Decimals <= 0 || Value == 0. )
110 {
111 return( (int)(0.5 + Value) );
112 }
113
114 Decimals = (int)(-(ceil(log10(fabs(Value))) - Decimals));
115
116 if( Decimals > 0 )
117 {
118 double d = pow(10., Decimals);
119
120 return( Value < 0.
121 ? -((int)(0.5 - Value * d)) / d
122 : ((int)(0.5 + Value * d)) / d
123 );
124 }
125 else
126 {
127 double d = pow(10., -Decimals);
128
129 return( Value < 0.
130 ? -((int)(0.5 - Value / d)) * d
131 : ((int)(0.5 + Value / d)) * d
132 );
133 }
134}
135
136
138// //
139// //
140// //
142
143//---------------------------------------------------------
144int SG_Get_Digit_Count(int Number)
145{
146 Number = abs(Number);
147
148 return( Number < 10 ? 1 : 1 + (int)log10((double)Number) );
149}
150
151
153// //
154// //
155// //
157
158//---------------------------------------------------------
159CSG_String SG_Get_Double_asString(double Number, int Width, int Precision, bool bScientific)
160{
161 if( bScientific )
162 {
163 if( Width > 0 && Precision >= 0 ) return( CSG_String::Format("%*.*e", Width, Precision, Number) );
164 if( Width > 0 ) return( CSG_String::Format("%*e" , Width , Number) );
165 if( Precision >= 0 ) return( CSG_String::Format("%.*e" , Precision, Number) );
166
167 return( CSG_String::Format("%e", Number) );
168 }
169 else
170 {
171 if( Width > 0 && Precision >= 0 ) return( CSG_String::Format("%*.*f", Width, Precision, Number) );
172 if( Width > 0 ) return( CSG_String::Format("%*f" , Width , Number) );
173 if( Precision >= 0 ) return( CSG_String::Format("%.*f" , Precision, Number) );
174
175 return( CSG_String::Format("%f", Number) );
176 }
177}
178
179
181// //
182// //
183// //
185
186//---------------------------------------------------------
187int SG_Compare_Int(const void *a, const void *b)
188{
189 if( *((int *)a) < *((int *)b) )
190 return( -1 );
191
192 if( *((int *)a) > *((int *)b) )
193 return( 1 );
194
195 return( 0 );
196}
197
198//---------------------------------------------------------
199int SG_Compare_Double(const void *a, const void *b)
200{
201 if( *((double *)a) < *((double *)b) )
202 return( -1 );
203
204 if( *((double *)a) > *((double *)b) )
205 return( 1 );
206
207 return( 0 );
208}
209
210//---------------------------------------------------------
211int SG_Compare_Char_Ptr(const void *a, const void *b)
212{
213 return( strcmp((const char *)a, (const char *)b) );
214}
215
216
218// //
219// //
220// //
222
223//---------------------------------------------------------
224double SG_Degree_To_Decimal(double Deg, double Min, double Sec)
225{
226 if( Deg < 0. ) { Min = -fabs(Min); Sec = -fabs(Sec); }
227
228 return( Deg + Min / 60. + Sec / 3600. );
229}
230
231//---------------------------------------------------------
232void SG_Decimal_To_Degree(double Value, double &Deg, double &Min, double &Sec)
233{
234 Sec = fmod(Value < 0. ? -Value : Value, 360.);
235
236 Deg = (int)Sec; Sec = 60. * (Sec - Deg);
237 Min = (int)Sec; Sec = 60. * (Sec - Min);
238
239 if( Value < 0. )
240 {
241 Deg = -Deg;
242 }
243}
244
245
247// //
248// //
249// //
251
252//---------------------------------------------------------
254{
255 Initialize();
256}
257
258//---------------------------------------------------------
260{
261 Initialize((unsigned)time(NULL));
262}
263
264//---------------------------------------------------------
265void CSG_Random::Initialize(unsigned int Value)
266{
267 srand(Value);
268}
269
270//---------------------------------------------------------
271// Uniform distributed pseudo-random numbers in the range from 0 to 1.
272//
274{
275 return( 1. * rand() / (double)RAND_MAX );
276}
277
278//---------------------------------------------------------
279// Uniform distributed pseudo-random numbers in the range from min to max.
280//
281double CSG_Random::Get_Uniform(double min, double max)
282{
283 return( min + (max - min) * rand() / (double)RAND_MAX );
284}
285
286//---------------------------------------------------------
287// Generating Gaussian pseudo-random numbers using
288// the polar form of the Box-Muller transformation.
289//
290// Box, G.E.P, Muller, M.E. (1958):
291// 'A note on the generation of random normal deviates',
292// Annals Math. Stat, V. 29, pp. 610-611
293//
294// Link: http://www.taygeta.com/random/gaussian.html
295//
296//---------------------------------------------------------
297double CSG_Random::Get_Gaussian(double mean, double stddev)
298{
299 double x1, x2, w;
300
301 do
302 {
303 x1 = 2. * Get_Uniform() - 1.;
304 x2 = 2. * Get_Uniform() - 1.;
305
306 w = x1 * x1 + x2 * x2;
307 }
308 while( w >= 1. );
309
310 w = sqrt((-2. * log(w)) / w);
311
312 return( mean + stddev * x1 * w );
313}
314
315
317// //
318// //
319// //
321
322//---------------------------------------------------------
327
329{
330 Create(bHoldValues);
331}
332
337
338CSG_Simple_Statistics::CSG_Simple_Statistics(double Mean, double StdDev, sLong Count)
339{
340 Create(Mean, StdDev, Count);
341}
342
344{
345 Create(Values, bHoldValues);
346}
347
348//---------------------------------------------------------
349bool CSG_Simple_Statistics::Create(bool bHoldValues)
350{
351 Invalidate();
352
353 m_Values.Create(bHoldValues ? sizeof(double) : 0, 0, TSG_Array_Growth::SG_ARRAY_GROWTH_1);
354
355 return( true );
356}
357
359{
360 m_bEvaluated = Statistics.m_bEvaluated;
361
362 m_nValues = Statistics.m_nValues;
363 m_Weights = Statistics.m_Weights;
364 m_Sum = Statistics.m_Sum;
365 m_Sum2 = Statistics.m_Sum2;
366
367 m_Minimum = Statistics.m_Minimum;
368 m_Maximum = Statistics.m_Maximum;
369 m_Range = Statistics.m_Range;
370 m_Mean = Statistics.m_Mean;
371 m_Variance = Statistics.m_Variance;
372 m_StdDev = Statistics.m_StdDev;
373
374 m_Kurtosis = Statistics.m_Kurtosis;
375 m_Skewness = Statistics.m_Skewness;
376
377 m_Gini = Statistics.m_Gini;
378
379 m_bSorted = Statistics.m_bSorted;
380 m_Values .Create(Statistics.m_Values);
381
382 return( true );
383}
384
385bool CSG_Simple_Statistics::Create(double Mean, double StdDev, sLong Count)
386{
387 Invalidate();
388
389 m_bEvaluated = 1;
390
391 m_Mean = Mean;
392 m_StdDev = StdDev;
393 m_Variance = StdDev*StdDev;
394 m_nValues = Count;
395 m_Weights = (double)Count;
396
399
400 m_Minimum = m_Mean - 1.5 * m_StdDev;
401 m_Maximum = m_Mean + 1.5 * m_StdDev;
403
404 return( true );
405}
406
407bool CSG_Simple_Statistics::Create(const CSG_Vector &Values, bool bHoldValues)
408{
409 if( Create(bHoldValues) )
410 {
411 for(sLong i=0; i<Values.Get_Size(); i++)
412 {
413 Add_Value(Values[i]);
414 }
415
416 return( true );
417 }
418
419 return( false );
420}
421
422//---------------------------------------------------------
424{
425 if( m_nValues < 1 || nValues < 1 || m_nValues == nValues )
426 {
427 return( false );
428 }
429
430 double Scale = nValues / (double)m_nValues;
431
432 m_Weights *= Scale;
433 m_Sum *= Scale;
434 m_Sum2 *= Scale;
435
436 m_nValues = nValues;
437
438 m_bEvaluated = 0;
439
440 m_Values.Destroy();
441
442 return( true );
443}
444
445//---------------------------------------------------------
447{
448 m_bEvaluated = 0;
449
450 m_nValues = 0;
451 m_Weights = 0.;
452 m_Sum = 0.;
453 m_Sum2 = 0.;
454
455 m_Minimum = 0.;
456 m_Maximum = 0.;
457 m_Range = 0.;
458 m_Mean = 0.;
459 m_Variance = 0.;
460 m_StdDev = 0.;
461
462 m_Kurtosis = 0.;
463 m_Skewness = 0.;
464
465 m_Gini = -1.;
466
467 m_bSorted = false;
468 m_Values .Destroy();
469}
470
471//---------------------------------------------------------
473{
474 _Evaluate();
475
476 return( is_Evaluated() > 0 );
477}
478
479//---------------------------------------------------------
481{
482 if( Statistics.m_nValues <= 0 )
483 {
484 return;
485 }
486
487 if( m_nValues == 0 )
488 {
489 Create(Statistics);
490
491 return;
492 }
493
494 //--------------------------------------------------------
495 if( (sLong)m_Values.Get_Size() == m_nValues && (sLong)Statistics.m_Values.Get_Size() == Statistics.m_nValues && m_Values.Set_Array((size_t)(m_nValues + Statistics.m_nValues)) )
496 {
497 for(sLong i=0, j=m_nValues; i<Statistics.m_nValues; i++, j++)
498 {
499 ((double *)m_Values.Get_Array())[j] = Statistics.Get_Value(i);
500 }
501 }
502 else
503 {
504 m_Values.Destroy();
505 }
506
507 m_nValues += Statistics.m_nValues;
508 m_Weights += Statistics.m_Weights;
509 m_Sum += Statistics.m_Sum;
510 m_Sum2 += Statistics.m_Sum2;
511
512 if( m_Minimum > Statistics.m_Minimum )
513 m_Minimum = Statistics.m_Minimum;
514
515 if( m_Maximum < Statistics.m_Maximum )
516 m_Maximum = Statistics.m_Maximum;
517
518 m_Kurtosis = 0.;
519 m_Skewness = 0.;
520
521 m_bEvaluated = 0;
522 m_bSorted = false;
523}
524
525//---------------------------------------------------------
526void CSG_Simple_Statistics::Add_Value(double Value, double Weight)
527{
528 if( m_nValues == 0 )
529 {
530 m_Minimum = m_Maximum = Value;
531 }
532 else if( m_Minimum > Value )
533 {
534 m_Minimum = Value;
535 }
536 else if( m_Maximum < Value )
537 {
538 m_Maximum = Value;
539 }
540
541 if( Weight )
542 {
543 m_Weights += fabs(Weight);
544 m_Sum += Weight * Value;
545 m_Sum2 += Weight * Value*Value;
546
547 m_bEvaluated = 0;
548
549 if( m_Values.Get_Value_Size() > 0 && m_Values.Inc_Array() )
550 {
551 m_bSorted = false;
552
553 ((double *)m_Values.Get_Array())[m_nValues] = Value;
554 }
555
556 m_nValues++;
557 }
558}
559
560//---------------------------------------------------------
562{
563 if( m_bEvaluated == 0 && m_Weights > 0. )
564 {
565 m_bEvaluated = 1;
566
570 m_StdDev = m_Variance > 0. ? sqrt(m_Variance) : 0.;
571 }
572
573 //-----------------------------------------------------
574 if( m_bEvaluated == 1 && Level > 1 )
575 {
576 m_bEvaluated = 2;
577
578 m_Kurtosis = 0.;
579 m_Skewness = 0.;
580
581 if( Get_StdDev() > 0. && m_Values.Get_Size() > 0 )
582 {
583 for(sLong i=0; i<Get_Count(); i++)
584 {
585 double d = (Get_Value(i) - Get_Mean()) / Get_StdDev();
586
587 m_Kurtosis += d*d*d*d;
588 m_Skewness += d*d*d;
589 }
590
593 // m_Skewness *= Get_Count() / ((Get_Count() - 1) * (Get_Count() - 2));
594 }
595 }
596}
597
598//---------------------------------------------------------
606{
607 return( Get_StdDev() != 0. ? (Get_Mean() - Get_Median()) / Get_StdDev() : 0. );
608}
609
610//---------------------------------------------------------
618{
619 if( m_Values.Get_Size() < 1 )
620 {
621 return( m_Mean );
622 }
623
624 //-----------------------------------------------------
625 if( !m_bSorted )
626 {
627 qsort(m_Values.Get_Array(), m_Values.Get_Size(), sizeof(double), SG_Compare_Double);
628
629 m_bSorted = true;
630 }
631
632 if( Quantile <= 0. || m_Values.Get_Size() == 1 )
633 {
634 return( Get_Values()[0] );
635 }
636
637 if( Quantile >= 1. )
638 {
639 return( Get_Values()[m_Values.Get_Size() - 1] );
640 }
641
642 //-----------------------------------------------------
643 double r = Quantile * (m_Values.Get_Size() - 1);
644
645 sLong i = (sLong)r; r -= i;
646
647 return( r == 0. ? Get_Values()[i] : ((1. - r) * Get_Values()[i] + r * Get_Values()[i + 1]) );
648}
649
650//---------------------------------------------------------
658{
659 return( Get_Quantile(Percentile / 100.) );
660}
661
662//---------------------------------------------------------
669{
670 return( Get_Quantile(0.5) );
671}
672
673//---------------------------------------------------------
681{
682 if( m_Gini < 0. && m_Values.Get_Size() > 1 )
683 {
684 if( !m_bSorted )
685 {
686 qsort(m_Values.Get_Array(), m_Values.Get_Size(), sizeof(double), SG_Compare_Double);
687
688 m_bSorted = true;
689 }
690
691 m_Gini = 0.;
692
693 for(sLong i=0; i<Get_Count(); i++)
694 {
695 m_Gini += (i + 1.) * Get_Value(i);
696 }
697
698 m_Gini = 2. * m_Gini / (Get_Count() * Get_Sum()) - (Get_Count() + 1.) / Get_Count();
699 }
700
701 return( m_Gini );
702}
703
704//---------------------------------------------------------
710{
711 if( m_Values.Get_Size() == 0 ) { return( -1 ); }
712
713 size_t Index = 0;
714 double Value = Get_Values()[Index];
715
716 for(size_t i=1; i<(size_t)m_Values.Get_Size(); i++)
717 {
718 if( Value > Get_Values()[i] )
719 {
720 Index = i;
721 Value = Get_Values()[i];
722 }
723 }
724
725 return( (sLong)Index );
726}
727
728//---------------------------------------------------------
734{
735 if( m_Values.Get_Size() == 0 ) { return( -1 ); }
736
737 size_t Index = 0;
738 double Value = Get_Values()[Index];
739
740 for(size_t i=1; i<(size_t)m_Values.Get_Size(); i++)
741 {
742 if( Value < Get_Values()[i] )
743 {
744 Index = i;
745 Value = Get_Values()[i];
746 }
747 }
748
749 return( (sLong)Index );
750}
751
752//---------------------------------------------------------
757sLong CSG_Simple_Statistics::Get_nValues_Above(double Threshold, bool bEquals)
758{
759 if( m_Values.Get_Size() == 0 ) { return( -1 ); }
760
761 sLong n = 0;
762
763 for(sLong i=0; i<Get_Count(); i++)
764 {
765 if( (bEquals && Get_Value(i) >= Threshold) || Get_Value(i) > Threshold )
766 {
767 n++;
768 }
769 }
770
771 return( n );
772}
773
774//---------------------------------------------------------
779sLong CSG_Simple_Statistics::Get_nValues_Below(double Threshold, bool bEquals)
780{
781 if( m_Values.Get_Size() == 0 ) { return( -1 ); }
782
783 sLong n = 0;
784
785 for(sLong i=0; i<Get_Count(); i++)
786 {
787 if( (bEquals && Get_Value(i) <= Threshold) || Get_Value(i) < Threshold )
788 {
789 n++;
790 }
791 }
792
793 return( n );
794}
795
796
798// //
799// //
800// //
802
803//---------------------------------------------------------
805{
806 int Index = 0;
807
808 bWeighted = bWeighted && m_bWeights;
809
810 for(int i=1; i<Get_Count(); i++)
811 {
812 if( bWeighted )
813 {
814 if( m_Weight[i] > m_Weight[Index] )
815 {
816 Index = i;
817 }
818 }
819 else
820 {
821 if( m_Count[i] > m_Count[Index] )
822 {
823 Index = i;
824 }
825 }
826 }
827
828 return( Index );
829}
830
831//---------------------------------------------------------
833{
834 int Index = 0;
835
836 bWeighted = bWeighted && m_bWeights;
837
838 for(int i=1; i<Get_Count(); i++)
839 {
840 if( bWeighted )
841 {
842 if( m_Weight[i] < m_Weight[Index] )
843 {
844 Index = i;
845 }
846 }
847 else
848 {
849 if( m_Count[i] < m_Count[Index] )
850 {
851 Index = i;
852 }
853 }
854 }
855
856 return( Index );
857}
858
859
861// //
863
864//---------------------------------------------------------
866{
867 m_bWeights = bWeights;
868
869 m_Count.Destroy();
870 m_Value.Destroy();
871}
872
873//---------------------------------------------------------
874void CSG_Unique_Number_Statistics::Add_Value(double Value, double Weight)
875{
876 for(int i=0; i<Get_Count(); i++)
877 {
878 if( Value == m_Value[i] )
879 {
880 m_Count[i]++;
881
882 if( m_bWeights && Weight > 0. )
883 {
884 m_Weight[i] += Weight;
885 }
886
887 return;
888 }
889 }
890
891 m_Count.Add(1);
892 m_Value.Add_Row(Value);
893
894 if( m_bWeights && Weight > 0. )
895 {
896 m_Weight.Add_Row(Weight);
897 }
898}
899
900//---------------------------------------------------------
902{
903 for(int i=0; i<Get_Count(); i++)
904 {
905 if( Value == m_Value[i] )
906 {
907 return( i );
908 }
909 }
910
911 return( -1 );
912}
913
914
916// //
918
919//---------------------------------------------------------
921{
922 m_bWeights = bWeights;
923
924 m_Count.Destroy();
925 m_Value.Clear();
926}
927
928//---------------------------------------------------------
930{
931 for(int i=0; i<Get_Count(); i++)
932 {
933 if( Value.Cmp(m_Value[i]) == 0 )
934 {
935 m_Count[i]++;
936
937 if( m_bWeights && Weight > 0. )
938 {
939 m_Weight[i] += Weight;
940 }
941
942 return;
943 }
944 }
945
946 m_Count.Add(1);
947 m_Value.Add(Value);
948
949 if( m_bWeights && Weight > 0. )
950 {
951 m_Weight.Add_Row(Weight);
952 }
953}
954
955//---------------------------------------------------------
957{
958 for(int i=0; i<Get_Count(); i++)
959 {
960 if( Value.Cmp(m_Value[i]) == 0 )
961 {
962 return( i );
963 }
964 }
965
966 return( -1 );
967}
968
969
971// //
972// //
973// //
975
976//---------------------------------------------------------
978{
979 m_pTable = new CSG_Table;
980
981 Create(Type);
982}
983
984//---------------------------------------------------------
986{
987 delete(m_pTable);
988}
989
990//---------------------------------------------------------
992{
993 m_pTable->Destroy();
994
995 m_pTable->Add_Field("VALUE", Type);
996 m_pTable->Add_Field("COUNT", SG_DATATYPE_ULong);
997}
998
999//---------------------------------------------------------
1001{
1002 m_pTable->Del_Records();
1003}
1004
1005//---------------------------------------------------------
1007{
1008 return( m_pTable->Get_Field_Type(0) );
1009}
1010
1011
1013// //
1015
1016//---------------------------------------------------------
1018{
1019 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1020
1021 return( pRecord ? (int)pRecord->Get_Index() : -1);
1022}
1023
1024//---------------------------------------------------------
1026{
1027 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1028
1029 return( pRecord ? (int)pRecord->Get_Index() : -1);
1030}
1031
1032//---------------------------------------------------------
1034{
1035 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1036
1037 return( pRecord ? (int)pRecord->Get_Index() : -1);
1038}
1039
1040
1042// //
1044
1045//---------------------------------------------------------
1047{
1048 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1049
1050 if( !pRecord )
1051 {
1052 (pRecord = m_pTable->Add_Record())->Set_Value(0, Value);
1053
1054 if( m_pTable->Get_Count() > 16 || m_pTable->is_Indexed() )
1055 {
1056 Sort();
1057 }
1058 }
1059
1060 pRecord->Add_Value(1, 1);
1061
1062 return( (int)(m_pTable->Get_Count() - 1) );
1063}
1064
1065//---------------------------------------------------------
1067{
1068 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1069
1070 if( !pRecord )
1071 {
1072 (pRecord = m_pTable->Add_Record())->Set_Value(0, Value);
1073
1074 if( m_pTable->Get_Count() > 16 || m_pTable->is_Indexed() )
1075 {
1076 Sort();
1077 }
1078 }
1079
1080 pRecord->Add_Value(1, 1);
1081
1082 return( (int)(m_pTable->Get_Count() - 1) );
1083}
1084
1085//---------------------------------------------------------
1087{
1088 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1089
1090 if( !pRecord )
1091 {
1092 (pRecord = m_pTable->Add_Record())->Set_Value(0, Value);
1093
1094 if( m_pTable->Get_Count() > 16 || m_pTable->is_Indexed() )
1095 {
1096 Sort();
1097 }
1098 }
1099
1100 pRecord->Add_Value(1, 1);
1101
1102 return( (int)(m_pTable->Get_Count() - 1) );
1103}
1104
1105//---------------------------------------------------------
1106// sort categories ascending
1108{
1109 CSG_Table Table(*m_pTable);
1110
1111 if( Table.Set_Index(0, TABLE_INDEX_Ascending) )
1112 {
1113 for(sLong i=0; i<m_pTable->Get_Count(); i++)
1114 {
1115 m_pTable->Set_Record(i, Table.Get_Record_byIndex(i));
1116 }
1117
1118 return( m_pTable->Set_Index(0, TABLE_INDEX_Ascending) );
1119 }
1120
1121 return( false );
1122}
1123
1124
1126// //
1128
1129//---------------------------------------------------------
1130// returns the number of categories.
1132{
1133 return( (int)m_pTable->Get_Count() );
1134}
1135
1136//---------------------------------------------------------
1137// returns the number of observations for the i'th category.
1139{
1140 CSG_Table_Record *pRecord = m_pTable->Get_Record(i);
1141
1142 return( pRecord ? pRecord->asInt(1) : 0 );
1143}
1144
1145//---------------------------------------------------------
1147{
1148 CSG_Table_Record *pRecord = m_pTable->Get_Record(i);
1149
1150 return( pRecord ? pRecord->asInt(0) : 0 );
1151}
1152
1153//---------------------------------------------------------
1155{
1156 CSG_Table_Record *pRecord = m_pTable->Get_Record(i);
1157
1158 return( pRecord ? pRecord->asDouble(0) : 0 );
1159}
1160
1161//---------------------------------------------------------
1163{
1164 CSG_Table_Record *pRecord = m_pTable->Get_Record(i);
1165
1166 return( pRecord ? pRecord->asString(0) : SG_T("") );
1167}
1168
1169
1171// //
1173
1174//---------------------------------------------------------
1176{
1177 if( m_pTable->Get_Count() > 0 )
1178 {
1179 int Index = 0, Count = m_pTable->Get_Record_byIndex(0)->asInt(1);
1180
1181 for(int i=1; i<m_pTable->Get_Count(); i++)
1182 {
1183 if( Count < m_pTable->Get_Record_byIndex(i)->asInt(1) )
1184 {
1185 Index = i; Count = m_pTable->Get_Record_byIndex(i)->asInt(1);
1186 }
1187 }
1188
1189 return( Index );
1190 }
1191
1192 return( -1 );
1193}
1194
1195//---------------------------------------------------------
1197{
1198 if( m_pTable->Get_Count() > 0 )
1199 {
1200 int Index = 0, Count = m_pTable->Get_Record_byIndex(0)->asInt(1);
1201
1202 for(int i=1; i<m_pTable->Get_Count(); i++)
1203 {
1204 if( Count > m_pTable->Get_Record_byIndex(i)->asInt(1) )
1205 {
1206 Index = i; Count = m_pTable->Get_Record_byIndex(i)->asInt(1);
1207 }
1208 }
1209
1210 return( Index );
1211 }
1212
1213 return( -1 );
1214}
1215
1216
1218// //
1219// //
1220// //
1222
1223//---------------------------------------------------------
1225{
1226 _On_Construction();
1227}
1228
1229//---------------------------------------------------------
1231{
1232 _On_Construction();
1233
1234 Create(Histogram);
1235}
1236
1237//---------------------------------------------------------
1238CSG_Histogram::CSG_Histogram(size_t nClasses, double Minimum, double Maximum)
1239{
1240 _On_Construction();
1241
1242 Create(nClasses, Minimum, Maximum);
1243}
1244
1245//---------------------------------------------------------
1246CSG_Histogram::CSG_Histogram(size_t nClasses, const CSG_Vector &Values, double Minimum, double Maximum, size_t maxSamples)
1247{
1248 _On_Construction();
1249
1250 Create(nClasses, Values, Minimum, Maximum, maxSamples);
1251}
1252
1253//---------------------------------------------------------
1254CSG_Histogram::CSG_Histogram(size_t nClasses, CSG_Table *pTable, int Field, double Minimum, double Maximum, size_t maxSamples, int Normalize, double Scale)
1255{
1256 _On_Construction();
1257
1258 Create(nClasses, pTable, Field, Minimum, Maximum, maxSamples, Normalize, Scale);
1259}
1260
1261//---------------------------------------------------------
1262CSG_Histogram::CSG_Histogram(size_t nClasses, class CSG_Grid *pGrid, double Minimum, double Maximum, size_t maxSamples)
1263{
1264 _On_Construction();
1265
1266 Create(nClasses, pGrid, Minimum, Maximum, maxSamples);
1267}
1268
1269//---------------------------------------------------------
1270CSG_Histogram::CSG_Histogram(size_t nClasses, CSG_Grids *pGrids, double Minimum, double Maximum, size_t maxSamples)
1271{
1272 _On_Construction();
1273
1274 Create(nClasses, pGrids, Minimum, Maximum, maxSamples);
1275}
1276
1277//---------------------------------------------------------
1279{
1280 Destroy();
1281}
1282
1283//---------------------------------------------------------
1285{
1286 m_Statistics.Create();
1287
1288 SG_FREE_SAFE(m_Elements );
1289 SG_FREE_SAFE(m_Cumulative);
1290
1291 _On_Construction();
1292
1293 return( true );
1294}
1295
1296
1298// //
1300
1301//---------------------------------------------------------
1302void CSG_Histogram::_On_Construction(void)
1303{
1304 m_nClasses = 0;
1305 m_Elements = NULL;
1306 m_Cumulative = NULL;
1307 m_Minimum = 0.;
1308 m_Maximum = 0.;
1309 m_ClassWidth = 1.;
1310}
1311
1312//---------------------------------------------------------
1313bool CSG_Histogram::_Create(size_t nClasses, double Minimum, double Maximum)
1314{
1315 if( nClasses > 0 && Minimum < Maximum )
1316 {
1317 Destroy();
1318
1319 m_Elements = (size_t *)SG_Calloc(nClasses, sizeof(size_t));
1320 m_Cumulative = (size_t *)SG_Calloc(nClasses, sizeof(size_t));
1321
1322 if( m_Elements && m_Cumulative )
1323 {
1324 m_nClasses = nClasses;
1325 m_Minimum = Minimum;
1326 m_Maximum = Maximum;
1327 m_ClassWidth = (Maximum - Minimum) / (double)m_nClasses;
1328
1329 return( true );
1330 }
1331 }
1332
1333 Destroy();
1334
1335 return( false );
1336}
1337
1338//---------------------------------------------------------
1340{
1341 m_Statistics += Value;
1342
1343 if( m_nClasses > 0 && m_Minimum <= Value && Value <= m_Maximum )
1344 {
1345 size_t Class = (size_t)((Value - m_Minimum) / m_ClassWidth);
1346
1347 if( Class >= m_nClasses )
1348 {
1349 Class = m_nClasses - 1;
1350 }
1351
1352 m_Elements[Class]++;
1353 }
1354}
1355
1356//---------------------------------------------------------
1358{
1359 if( m_nClasses > 0 && Scale > 0. )
1360 {
1361 m_Statistics.Set_Count((sLong)(Scale * Get_Element_Count()));
1362
1363 for(size_t i=0; i<m_nClasses; i++)
1364 {
1365 m_Elements[i] = (size_t)(Scale * m_Elements[i]);
1366 }
1367
1368 return( Update() );
1369 }
1370
1371 return( false );
1372}
1373
1374//---------------------------------------------------------
1376{
1377 if( m_nClasses > 0 )
1378 {
1379 m_Statistics.Get_Mean(); // _Evaluate()
1380
1381 m_nMaximum = m_Cumulative[0] = m_Elements[0];
1382
1383 for(size_t i=1; i<m_nClasses; i++)
1384 {
1385 m_Cumulative[i] = m_Cumulative[i - 1] + m_Elements[i];
1386
1387 if( m_nMaximum < m_Elements[i] )
1388 {
1389 m_nMaximum = m_Elements[i];
1390 }
1391 }
1392
1393 return( Get_Element_Count() > 0 );
1394 }
1395
1396 return( false );
1397}
1398
1399//---------------------------------------------------------
1400bool CSG_Histogram::_Update(sLong nElements)
1401{
1402 if( nElements > 0 && m_Statistics.Get_Count() > 0 )
1403 {
1404 double Scale = (double)nElements / (double)m_Statistics.Get_Count();
1405
1406 m_Statistics.Create(m_Statistics.Get_Mean(), m_Statistics.Get_StdDev(), nElements);
1407
1408 for(size_t i=1; i<m_nClasses; i++)
1409 {
1410 m_Elements[i] = (size_t)(0.5 + Scale * m_Elements[i]);
1411 }
1412 }
1413
1414 return( Update() );
1415}
1416
1418{
1419 if( m_nClasses == 0
1420 || m_nClasses != Histogram.m_nClasses
1421 || m_Minimum != Histogram.m_Minimum
1422 || m_Maximum != Histogram.m_Maximum )
1423 {
1424 return( false );
1425 }
1426
1427 m_Statistics += Histogram.m_Statistics;
1428
1429 for( size_t i=0; i<m_nClasses; i++ )
1430 {
1431 m_Elements[i] += Histogram.m_Elements[i];
1432 }
1433
1434 return( true );
1435}
1436
1437//---------------------------------------------------------
1441double CSG_Histogram::Get_Quantile(double Quantile) const
1442{
1443 if( m_nClasses < 2 ) { return( 0. ); }
1444
1445 if( Quantile <= 0. ) { return( m_Minimum ); }
1446 if( Quantile >= 1. ) { return( m_Maximum ); }
1447
1448 size_t n = (size_t)(Quantile * Get_Element_Count()); // number of elements
1449
1450 for(size_t i=0, n0=0; i<m_nClasses; n0=m_Cumulative[i++])
1451 {
1452 if( n < m_Cumulative[i] )
1453 {
1454 if( m_Cumulative[i] >= n0 )
1455 {
1456 return( Get_Center(i) );
1457 }
1458
1459 double d = (n - n0) / (double)(m_Cumulative[i] - n0);
1460
1461 return( Get_Break(i) + d * m_ClassWidth );
1462 }
1463 else if( n == m_Cumulative[i] )
1464 {
1465 return( Get_Break(i + 1) );
1466 }
1467 }
1468
1469 return( m_Maximum );
1470}
1471
1472//---------------------------------------------------------
1476double CSG_Histogram::Get_Percentile(double Percentile) const
1477{
1478 return( Get_Quantile(Percentile / 100.) );
1479}
1480
1481//---------------------------------------------------------
1485double CSG_Histogram::Get_Quantile_Value(double Value) const
1486{
1487 if( m_nClasses < 2 ) { return( 0. ); }
1488
1489 if( Value <= m_Minimum ) { return( 0. ); }
1490 if( Value >= m_Maximum ) { return( 1. ); }
1491
1492 size_t Class = (size_t)(m_nClasses * (Value - m_Minimum) / (m_Maximum - m_Minimum));
1493
1494 if( Class >= m_nClasses )
1495 {
1496 return( 1. );
1497 }
1498
1499 if( Class < 1 )
1500 {
1501 double dq = m_Cumulative[Class] / (double)Get_Element_Count();
1502
1503 return( dq * (Value - m_Minimum) / m_ClassWidth );
1504 }
1505
1506 double q0 = m_Cumulative[Class - 1] / (double)Get_Element_Count();
1507 double dq = (m_Cumulative[Class ] / (double)Get_Element_Count()) - q0;
1508
1509 return( q0 + dq * (Value - Get_Break(Class)) / m_ClassWidth );
1510}
1511
1512//---------------------------------------------------------
1516double CSG_Histogram::Get_Percentile_Value(double Value) const
1517{
1518 return( Get_Quantile_Value(Value) * 100. );
1519}
1520
1521
1523// //
1525
1526//---------------------------------------------------------
1528{
1529 if( !_Create(Histogram.m_nClasses, Histogram.m_Minimum, Histogram.m_Maximum) )
1530 {
1531 return( false );
1532 }
1533
1534 m_Statistics = Histogram.m_Statistics;
1535 m_ClassWidth = Histogram.m_ClassWidth;
1536 m_nMaximum = Histogram.m_nMaximum ;
1537
1538 for(size_t i=0; i<m_nClasses; i++)
1539 {
1540 m_Cumulative[i] = Histogram.m_Cumulative[i];
1541 m_Elements [i] = Histogram.m_Elements [i];
1542 }
1543
1544 return( true );
1545}
1546
1547//---------------------------------------------------------
1548bool CSG_Histogram::Create(size_t nClasses, double Minimum, double Maximum)
1549{
1550 return( _Create(nClasses, Minimum, Maximum) );
1551}
1552
1553//---------------------------------------------------------
1554bool CSG_Histogram::Create(size_t nClasses, const CSG_Vector &Values, double Minimum, double Maximum, size_t maxSamples)
1555{
1556 if( Minimum >= Maximum )
1557 {
1558 CSG_Simple_Statistics Statistics(Values);
1559
1560 Minimum = Statistics.Get_Minimum();
1561 Maximum = Statistics.Get_Maximum();
1562 }
1563
1564 if( !_Create(nClasses, Minimum, Maximum) )
1565 {
1566 return( false );
1567 }
1568
1569 //-----------------------------------------------------
1570 if( maxSamples > 0 && maxSamples < (size_t)Values.Get_N() )
1571 {
1572 double d = (double)Values.Get_N() / (double)maxSamples;
1573
1574 for(double i=0; i<(double)Values.Get_N(); i+=d)
1575 {
1576 Add_Value(Values[(sLong)i]);
1577 }
1578
1579 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1580
1581 return( _Update(d < 1. ? (int)(d * (double)Values.Get_N()) : Values.Get_N()) );
1582 }
1583
1584 //-----------------------------------------------------
1585 for(int i=0; i<Values.Get_N(); i++)
1586 {
1587 Add_Value(Values[i]);
1588 }
1589
1590 return( Update() );
1591}
1592
1593//---------------------------------------------------------
1594bool CSG_Histogram::Create(size_t nClasses, CSG_Table *pTable, int Field, double Minimum, double Maximum, size_t maxSamples, int Normalize, double Scale)
1595{
1596 if( !pTable || Field < 0 || Field >= pTable->Get_Field_Count() || !_Create(nClasses,
1598 Minimum < Maximum ? Maximum : pTable->Get_Maximum(Field)) )
1599 {
1600 return( false );
1601 }
1602
1603 //-----------------------------------------------------
1604 if( maxSamples > 0 && maxSamples < (size_t)pTable->Get_Count() )
1605 {
1606 double Value, Dividend, d = (double)pTable->Get_Count() / (double)maxSamples;
1607
1608 for(double i=0; i<(double)pTable->Get_Count(); i+=d)
1609 {
1610 if( pTable->Get_Value((sLong)i, Field, Value) )
1611 {
1612 if( Normalize < 0 )
1613 {
1614 Add_Value(Value);
1615 }
1616 else if( pTable->Get_Value((sLong)i, Normalize, Dividend) && Dividend != 0. )
1617 {
1618 Add_Value(Scale * Value / Dividend);
1619 }
1620 }
1621 }
1622
1623 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1624
1625 return( _Update(d < 1. ? (int)(d * pTable->Get_Count()) : pTable->Get_Count()) );
1626 }
1627
1628 //-----------------------------------------------------
1629 for(sLong i=0; i<pTable->Get_Count(); i++)
1630 {
1631 double Value, Dividend;
1632
1633 if( pTable->Get_Value(i, Field, Value) )
1634 {
1635 if( Normalize < 0 )
1636 {
1637 Add_Value(Value);
1638 }
1639 else if( pTable->Get_Value((sLong)i, Normalize, Dividend) && Dividend != 0. )
1640 {
1641 Add_Value(Scale * Value / Dividend);
1642 }
1643 }
1644 }
1645
1646 return( Update() );
1647}
1648
1649//---------------------------------------------------------
1650bool CSG_Histogram::Create(size_t nClasses, CSG_Grid *pGrid, double Minimum, double Maximum, size_t maxSamples)
1651{
1652 if( !pGrid || !_Create(nClasses,
1655 {
1656 return( false );
1657 }
1658
1659 //-----------------------------------------------------
1660 if( maxSamples > 0 && (sLong)maxSamples < pGrid->Get_NCells() )
1661 {
1662 double d = (double)pGrid->Get_NCells() / (double)maxSamples;
1663
1664 for(double i=0; i<(double)pGrid->Get_NCells(); i+=d)
1665 {
1666 if( !pGrid->is_NoData((sLong)i) )
1667 {
1668 Add_Value(pGrid->asDouble((sLong)i));
1669 }
1670 }
1671
1672 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1673
1674 return( _Update(d < 1. ? (sLong)(d * (double)pGrid->Get_NCells()) : pGrid->Get_NCells()) );
1675 }
1676
1677 //-----------------------------------------------------
1678 for(sLong i=0; i<pGrid->Get_NCells(); i++)
1679 {
1680 if( !pGrid->is_NoData(i) )
1681 {
1682 Add_Value(pGrid->asDouble(i));
1683 }
1684 }
1685
1686 return( Update() );
1687}
1688
1689//---------------------------------------------------------
1690bool CSG_Histogram::Create(size_t nClasses, CSG_Grids *pGrids, double Minimum, double Maximum, size_t maxSamples)
1691{
1692 if( !pGrids || !_Create(nClasses,
1695 {
1696 return( false );
1697 }
1698
1699 //-----------------------------------------------------
1700 if( maxSamples > 0 && (sLong)maxSamples < pGrids->Get_NCells() )
1701 {
1702 double d = (double)pGrids->Get_NCells() / (double)maxSamples;
1703
1704 for(double i=0; i<(double)pGrids->Get_NCells(); i+=d)
1705 {
1706 if( !pGrids->is_NoData((sLong)i) )
1707 {
1708 Add_Value(pGrids->asDouble((sLong)i));
1709 }
1710 }
1711
1712 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1713
1714 return( _Update(d < 1. ? (sLong)(d * (double)pGrids->Get_NCells()) : pGrids->Get_NCells()) );
1715 }
1716
1717 //-----------------------------------------------------
1718 for(sLong i=0; i<pGrids->Get_NCells(); i++)
1719 {
1720 if( !pGrids->is_NoData(i) )
1721 {
1722 Add_Value(pGrids->asDouble(i));
1723 }
1724 }
1725
1726 return( Update() );
1727}
1728
1729
1731// //
1733
1734//---------------------------------------------------------
1736{
1737 Create(Histogram);
1738
1739 return( *this );
1740}
1741
1742
1744// //
1745// //
1746// //
1748
1749//---------------------------------------------------------
1752
1753//---------------------------------------------------------
1756
1757//---------------------------------------------------------
1758CSG_Natural_Breaks::CSG_Natural_Breaks(CSG_Table *pTable, int Field, int nClasses, int Histogram)
1759{
1760 Create(pTable, Field, nClasses, Histogram);
1761}
1762
1763//---------------------------------------------------------
1764CSG_Natural_Breaks::CSG_Natural_Breaks(CSG_Grid *pGrid, int nClasses, int Histogram)
1765{
1766 Create(pGrid, nClasses, Histogram);
1767}
1768
1769//---------------------------------------------------------
1770CSG_Natural_Breaks::CSG_Natural_Breaks(CSG_Grids *pGrids, int nClasses, int Histogram)
1771{
1772 Create(pGrids, nClasses, Histogram);
1773}
1774
1775//---------------------------------------------------------
1776CSG_Natural_Breaks::CSG_Natural_Breaks(const CSG_Vector &Values, int nClasses, int Histogram)
1777{
1778 Create(Values, nClasses, Histogram);
1779}
1780
1781
1783// //
1785
1786//---------------------------------------------------------
1787bool CSG_Natural_Breaks::Create(CSG_Table *pTable, int Field, int nClasses, int Histogram)
1788{
1789 bool bResult = false;
1790
1791 if( Histogram > 0 )
1792 {
1793 bResult = m_Histogram.Create(Histogram, pTable, Field) && _Histogram(nClasses);
1794 }
1795 else if( Field >= 0 && Field < pTable->Get_Field_Count() )
1796 {
1797 for(sLong i=0; i<pTable->Get_Count(); i++)
1798 {
1799 double Value;
1800
1801 if( pTable->Get_Value(i, Field, Value) )
1802 {
1803 m_Values.Add_Row(Value);
1804 }
1805 }
1806
1807 bResult = m_Values.Sort() && _Calculate(nClasses);
1808
1809 m_Values.Destroy();
1810 }
1811
1812 return( bResult );
1813}
1814
1815//---------------------------------------------------------
1816bool CSG_Natural_Breaks::Create(CSG_Grid *pGrid, int nClasses, int Histogram)
1817{
1818 bool bResult = false;
1819
1820 if( Histogram > 0 )
1821 {
1822 bResult = m_Histogram.Create(Histogram, pGrid) && _Histogram(nClasses);
1823 }
1824 else
1825 {
1826 for(sLong i=0; i<pGrid->Get_NCells(); i++)
1827 {
1828 if( !pGrid->is_NoData(i) )
1829 {
1830 m_Values.Add_Row(pGrid->asDouble(i));
1831 }
1832 }
1833
1834 bResult = m_Values.Sort() && _Calculate(nClasses);
1835
1836 m_Values.Destroy();
1837 }
1838
1839 return( bResult );
1840}
1841
1842//---------------------------------------------------------
1843bool CSG_Natural_Breaks::Create(CSG_Grids *pGrids, int nClasses, int Histogram)
1844{
1845 bool bResult = false;
1846
1847 if( Histogram > 0 )
1848 {
1849 bResult = m_Histogram.Create(Histogram, pGrids) && _Histogram(nClasses);
1850 }
1851 else
1852 {
1853 for(sLong i=0; i<pGrids->Get_NCells(); i++)
1854 {
1855 if( !pGrids->is_NoData(i) )
1856 {
1857 m_Values.Add_Row(pGrids->asDouble(i));
1858 }
1859 }
1860
1861 bResult = m_Values.Sort() && _Calculate(nClasses);
1862
1863 m_Values.Destroy();
1864 }
1865
1866 return( bResult );
1867}
1868
1869//---------------------------------------------------------
1870bool CSG_Natural_Breaks::Create(const CSG_Vector &Values, int nClasses, int Histogram)
1871{
1872 bool bResult = false;
1873
1874 if( Histogram > 0 )
1875 {
1876 bResult = m_Histogram.Create(Histogram, Values) && _Histogram(nClasses);
1877 }
1878 else
1879 {
1880 bResult = m_Values.Create(Values) && m_Values.Sort() && _Calculate(nClasses);
1881
1882 m_Values.Destroy();
1883 }
1884
1885 return( bResult );
1886}
1887
1888
1890// //
1892
1893//---------------------------------------------------------
1894bool CSG_Natural_Breaks::_Histogram(int nClasses)
1895{
1896 if( _Calculate(nClasses) )
1897 {
1898 double d = (double)m_Histogram.Get_Class_Count() / m_Histogram.Get_Cumulative((int)(m_Histogram.Get_Class_Count() - 1));
1899
1900 m_Breaks[0] = m_Histogram.Get_Break(0);
1901
1902 for(int i=1; i<Get_Count(); i++)
1903 {
1904 m_Breaks[i] = m_Histogram.Get_Value(m_Breaks[i] * d);
1905 }
1906
1907 m_Breaks[nClasses] = m_Histogram.Get_Break((int)m_Histogram.Get_Class_Count());
1908
1909 m_Histogram.Destroy();
1910
1911 return( true );
1912 }
1913
1914 m_Histogram.Destroy();
1915
1916 return( false );
1917}
1918
1919//---------------------------------------------------------
1920inline double CSG_Natural_Breaks::_Get_Value(int i)
1921{
1922 if( m_Histogram.Get_Class_Count() > 0 )
1923 {
1924 return( (double)m_Histogram.Get_Cumulative(i) );
1925 }
1926
1927 return( m_Values[i] );
1928}
1929
1930//---------------------------------------------------------
1931bool CSG_Natural_Breaks::_Calculate(int nClasses)
1932{
1933 if( m_Histogram.Get_Class_Count() == 0 && m_Values.Get_Size() == 0 )
1934 {
1935 return( false );
1936 }
1937
1938 int nValues = m_Histogram.Get_Class_Count() > 0 ? (int)m_Histogram.Get_Class_Count() : m_Values.Get_N();
1939
1940 CSG_Matrix mv(nClasses, nValues); mv.Assign(FLT_MAX);
1941
1942 int **mc = (int **)SG_Malloc(nValues * sizeof(int *));
1943
1944 mc[0] = (int *)SG_Calloc((size_t)nClasses * nValues, sizeof(int));
1945
1946 for(int i=0; i<nValues; i++)
1947 {
1948 mc[i] = mc[0] + i * (size_t)nClasses;
1949 }
1950
1951 //-----------------------------------------------------
1952 for(int i=1; i<nValues; i++)
1953 {
1954 double v = 0., s1 = 0., s2 = 0., w = 0.;
1955
1956 for(int m=0, n=i+1; m<=i; m++, n--)
1957 {
1958 v = _Get_Value(n);
1959 s2 += v*v;
1960 s1 += v;
1961 w ++;
1962 v = s2 - (s1 * s1) / w;
1963
1964 if( n > 0 )
1965 {
1966 for(int j=1; j<nClasses; j++)
1967 {
1968 if( mv[i][j] >= (v + mv[n - 1][j - 1]) )
1969 {
1970 mc[i][j] = n;
1971 mv[i][j] = v + mv[n - 1][j - 1];
1972 }
1973 }
1974 }
1975 }
1976
1977 mc[i][0] = 0;
1978 mv[i][0] = v;
1979 }
1980
1981 //-----------------------------------------------------
1982 CSG_Array_Int Class(nClasses);
1983
1984 for(int i=0; i<nClasses; i++)
1985 {
1986 Class[i] = i;
1987 }
1988
1989 int j = Class[(size_t)nClasses - 1] = nValues - 1;
1990
1991 for(int i=nClasses-1; i>0; i--)
1992 {
1993 Class[(size_t)i - 1] = j = mc[j - 1][i];
1994 }
1995
1996 //-----------------------------------------------------
1997 m_Breaks.Create((size_t)nClasses + 1);
1998
1999 m_Breaks[0] = _Get_Value(0);
2000
2001 for(int i=1; i<nClasses; i++)
2002 {
2003 m_Breaks[i] = _Get_Value(Class[i - 1]);
2004 }
2005
2006 m_Breaks[nClasses] = _Get_Value(nValues - 1);
2007
2008 SG_Free(mc[0]); SG_Free(mc);
2009
2010 return( true );
2011}
2012
2013
2015// //
2016// //
2017// //
2019
2020//---------------------------------------------------------
2022{
2023 m_nFeatures = 0;
2024 m_Iteration = 0;
2025}
2026
2027//---------------------------------------------------------
2032
2033//---------------------------------------------------------
2035{
2036 m_Centroid.Destroy();
2037 m_Variance.Destroy();
2038 m_nMembers.Destroy();
2039 m_Clusters.Destroy();
2040 m_Features.Destroy();
2041 m_nFeatures = 0;
2042 m_Iteration = 0;
2043
2044 return( true );
2045}
2046
2047//---------------------------------------------------------
2049{
2050 Destroy();
2051
2052 if( nFeatures > 0 )
2053 {
2054 m_nFeatures = nFeatures;
2055
2056 m_Features.Create(m_nFeatures * sizeof(double), 0, TSG_Array_Growth::SG_ARRAY_GROWTH_3);
2057
2058 return( true );
2059 }
2060
2061 return( false );
2062}
2063
2064//---------------------------------------------------------
2066{
2067 return( m_nFeatures > 0 && m_Features.Inc_Array() );
2068}
2069
2070//---------------------------------------------------------
2071bool CSG_Cluster_Analysis::Set_Feature(sLong iElement, int iFeature, double Value)
2072{
2073 if( iElement >= 0 && iElement < Get_nElements() && iFeature >= 0 && iFeature < m_nFeatures )
2074 {
2075 ((double *)m_Features.Get_Entry(iElement))[iFeature] = Value;
2076
2077 return( true );
2078 }
2079
2080 return( false );
2081}
2082
2083//---------------------------------------------------------
2093//---------------------------------------------------------
2094bool CSG_Cluster_Analysis::Execute(int Method, int nClusters, int nMaxIterations, int Initialization)
2095{
2096 if( Get_nElements() < 2 || nClusters < 2 )
2097 {
2098 return( false );
2099 }
2100
2101 //-----------------------------------------------------
2102 m_nMembers.Create(nClusters);
2103 m_Variance.Create(nClusters);
2104 m_Centroid.Create(m_nFeatures, nClusters);
2105
2106 //-----------------------------------------------------
2107 m_Clusters.Create(Get_nElements());
2108
2109 for(int iElement=0; iElement<Get_nElements(); iElement++)
2110 {
2111 switch( Initialization )
2112 {
2113 default: // random
2114 if( (m_Clusters[iElement] = (int)CSG_Random::Get_Uniform(0, nClusters)) >= nClusters )
2115 {
2116 m_Clusters[iElement] = nClusters - 1;
2117 }
2118 break;
2119
2120 case 1: // periodic
2121 {
2122 m_Clusters[iElement] = iElement % nClusters;
2123 }
2124 break;
2125
2126 case 2: // keep as is, but check for valid cluster ids
2127 if( 0 > m_Clusters[iElement] || m_Clusters[iElement] >= nClusters )
2128 {
2129 m_Clusters[iElement] = iElement % nClusters;
2130 }
2131 break;
2132 }
2133 }
2134
2135 //-----------------------------------------------------
2136 bool bResult;
2137
2138 m_Iteration = 0;
2139
2140 switch( Method )
2141 {
2142 default: bResult = _Minimum_Distance(true , nMaxIterations); break;
2143 case 1: bResult = _Hill_Climbing (true , nMaxIterations); break;
2144 case 2: bResult = _Minimum_Distance(true , nMaxIterations)
2145 && _Hill_Climbing (false, nMaxIterations); break;
2146 }
2147
2148 //-----------------------------------------------------
2149 if( bResult )
2150 {
2151 for(int iCluster=0; iCluster<nClusters; iCluster++)
2152 {
2153 m_Variance[iCluster] = m_nMembers[iCluster] <= 0 ? 0. : m_Variance[iCluster] / m_nMembers[iCluster];
2154 }
2155 }
2156
2157 return( bResult );
2158}
2159
2160//---------------------------------------------------------
2161bool CSG_Cluster_Analysis::_Minimum_Distance(bool bInitialize, int nMaxIterations)
2162{
2163 int iElement, iCluster, nClusters = m_Variance.Get_N();
2164
2165 double SP_Last = -1.;
2166
2167 //-----------------------------------------------------
2168 for(m_Iteration=1; SG_UI_Process_Get_Okay(); m_Iteration++)
2169 {
2170 m_Variance = 0.;
2171 m_Centroid = 0.;
2172 m_nMembers = 0;
2173
2174 //-------------------------------------------------
2175 for(iElement=0; iElement<Get_nElements(); iElement++)
2176 {
2177 m_nMembers[iCluster = m_Clusters[iElement]]++;
2178
2179 double *Feature = (double *)m_Features.Get_Entry(iElement);
2180
2181 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2182 {
2183 m_Centroid[iCluster][iFeature] += Feature[iFeature];
2184 }
2185 }
2186
2187 //-------------------------------------------------
2188 for(iCluster=0; iCluster<nClusters; iCluster++)
2189 {
2190 double d = m_nMembers[iCluster] > 0 ? 1. / m_nMembers[iCluster] : 0.;
2191
2192 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2193 {
2194 m_Centroid[iCluster][iFeature] *= d;
2195 }
2196 }
2197
2198 //-------------------------------------------------
2199 int nShifts = 0; m_SP = 0.;
2200
2201 for(iElement=0; iElement<Get_nElements(); iElement++)
2202 {
2203 double *Feature = (double *)m_Features.Get_Entry(iElement);
2204
2205 double minVariance = -1.;
2206 int minCluster = -1;
2207
2208 for(iCluster=0; iCluster<nClusters; iCluster++)
2209 {
2210 double Variance = 0.;
2211
2212 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2213 {
2214 Variance += SG_Get_Square(m_Centroid[iCluster][iFeature] - Feature[iFeature]);
2215 }
2216
2217 if( minVariance < 0. || Variance < minVariance )
2218 {
2219 minVariance = Variance;
2220 minCluster = iCluster;
2221 }
2222 }
2223
2224 if( m_Clusters[iElement] != minCluster )
2225 {
2226 m_Clusters[iElement] = minCluster;
2227
2228 nShifts++;
2229 }
2230
2231 m_SP += minVariance;
2232 m_Variance[minCluster] += minVariance;
2233 }
2234
2235 //-------------------------------------------------
2236 m_SP /= Get_nElements();
2237
2239 _TL("pass" ), m_Iteration,
2240 _TL("change"), m_Iteration < 2 ? m_SP : SP_Last - m_SP
2241 ));
2242
2243 SP_Last = m_SP;
2244
2245 if( nShifts == 0 || (nMaxIterations > 0 && nMaxIterations <= m_Iteration) )
2246 {
2247 return( true );
2248 }
2249 }
2250
2251 return( true );
2252}
2253
2254//---------------------------------------------------------
2255bool CSG_Cluster_Analysis::_Hill_Climbing(bool bInitialize, int nMaxIterations)
2256{
2257 int iElement, iCluster, nClusters = m_Variance.Get_N();
2258
2259 //-----------------------------------------------------
2260 m_Variance = 0.;
2261 m_Centroid = 0.;
2262 m_nMembers = 0;
2263
2264 for(iElement=0; iElement<Get_nElements(); iElement++)
2265 {
2266 m_nMembers[iCluster = m_Clusters[iElement]]++;
2267
2268 double *Feature = (double *)m_Features.Get_Entry(iElement);
2269
2270 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2271 {
2272 double d = Feature[iFeature];
2273
2274 m_Centroid[iCluster][iFeature] += d;
2275 m_Variance[iCluster] += d*d;
2276 }
2277 }
2278
2279 //-----------------------------------------------------
2280 for(iCluster=0; iCluster<nClusters; iCluster++)
2281 {
2282 double v = 0., d = m_nMembers[iCluster] <= 0 ? 0. : 1. / (double)m_nMembers[iCluster];
2283
2284 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2285 {
2286 m_Centroid[iCluster][iFeature] *= d;
2287 v += SG_Get_Square(m_Centroid[iCluster][iFeature]);
2288 }
2289
2290 m_Variance[iCluster] -= v * m_nMembers[iCluster];
2291 }
2292
2293 //-----------------------------------------------------
2294 double SP_Last = -1.; int noShift = 0;
2295
2296 for(m_Iteration=1; SG_UI_Process_Get_Okay(false); m_Iteration++)
2297 {
2298 for(iElement=0; iElement<Get_nElements(); iElement++)
2299 {
2300 iCluster = m_Clusters[iElement];
2301
2302 if( noShift++ < Get_nElements() && m_nMembers[iCluster] > 1 )
2303 {
2304 double *Feature = (double *)m_Features.Get_Entry(iElement);
2305
2306 double Variance = 0.;
2307
2308 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2309 {
2310 Variance += SG_Get_Square(m_Centroid[iCluster][iFeature] - Feature[iFeature]);
2311 }
2312
2313 double V1 = Variance * m_nMembers[iCluster] / (m_nMembers[iCluster] - 1.);
2314
2315 //-----------------------------------------
2316 int kCluster = 0; double VMin = -1.;
2317
2318 for(int jCluster=0; jCluster<nClusters; jCluster++)
2319 {
2320 if( jCluster != iCluster )
2321 {
2322 Variance = 0.;
2323
2324 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2325 {
2326 Variance += SG_Get_Square(m_Centroid[jCluster][iFeature] - Feature[iFeature]);
2327 }
2328
2329 double V2 = Variance * m_nMembers[jCluster] / (m_nMembers[jCluster] + 1.);
2330
2331 if( VMin < 0. || V2 < VMin )
2332 {
2333 VMin = V2;
2334 kCluster = jCluster;
2335 }
2336 }
2337 }
2338
2339 //-----------------------------------------
2340 if( VMin >= 0 && VMin < V1 )
2341 {
2342 noShift = 0;
2343 m_Variance[iCluster] -= V1;
2344 m_Variance[kCluster] += VMin;
2345 V1 = 1. / (m_nMembers[iCluster] - 1.);
2346 double V2 = 1. / (m_nMembers[kCluster] + 1.);
2347
2348 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2349 {
2350 double d = Feature[iFeature];
2351
2352 m_Centroid[iCluster][iFeature] = (m_nMembers[iCluster] * m_Centroid[iCluster][iFeature] - d) * V1;
2353 m_Centroid[kCluster][iFeature] = (m_nMembers[kCluster] * m_Centroid[kCluster][iFeature] + d) * V2;
2354 }
2355
2356 m_Clusters[iElement] = kCluster;
2357
2358 m_nMembers[iCluster]--;
2359 m_nMembers[kCluster]++;
2360 }
2361 }
2362 }
2363
2364 //-------------------------------------------------
2365 for(iCluster=0, m_SP=0.; iCluster<nClusters; iCluster++)
2366 {
2367 m_SP += m_Variance[iCluster];
2368 }
2369
2370 m_SP /= Get_nElements();
2371
2373 _TL("pass" ), m_Iteration,
2374 _TL("change"), m_Iteration <= 1 ? m_SP : SP_Last - m_SP
2375 ));
2376
2377 SP_Last = m_SP;
2378
2379 if( noShift >= Get_nElements() || (nMaxIterations > 0 && nMaxIterations <= m_Iteration) )
2380 {
2381 return( true );
2382 }
2383 }
2384
2385 return( true );
2386}
2387
2388
2390// //
2391// //
2392// //
2394
2395//---------------------------------------------------------
2397{
2398 m_nFeatures = 0; m_pClasses = NULL; m_nClasses = 0;
2399
2400 m_Threshold_Distance = 0.;
2401 m_Threshold_Angle = 0.;
2402 m_Threshold_Probability = 0.;
2403 m_Probability_Relative = false;
2404
2405 for(int i=0; i<SG_CLASSIFY_SUPERVISED_WTA; i++)
2406 {
2408 // || i == SG_CLASSIFY_SUPERVISED_Mahalonobis
2411 }
2412}
2413
2414//---------------------------------------------------------
2419
2420//---------------------------------------------------------
2422{
2423 Destroy();
2424
2425 if( nFeatures > 0 )
2426 {
2427 m_nFeatures = nFeatures;
2428 }
2429}
2430
2431//---------------------------------------------------------
2433{
2434 if( m_nClasses > 0 )
2435 {
2436 for(int i=0; i<m_nClasses; i++)
2437 {
2438 delete(m_pClasses[i]);
2439 }
2440
2441 SG_FREE_SAFE(m_pClasses);
2442 }
2443
2444 m_nFeatures = 0;
2445
2446 m_Info.Clear();
2447}
2448
2449
2451// //
2453
2454//---------------------------------------------------------
2455void CSG_Classifier_Supervised::Set_Threshold_Distance (double Value) { m_Threshold_Distance = Value; }
2456double CSG_Classifier_Supervised::Get_Threshold_Distance (void) { return( m_Threshold_Distance ); }
2457
2458//---------------------------------------------------------
2459void CSG_Classifier_Supervised::Set_Threshold_Angle (double Value) { m_Threshold_Angle = Value; }
2460double CSG_Classifier_Supervised::Get_Threshold_Angle (void) { return( m_Threshold_Angle ); }
2461
2462//---------------------------------------------------------
2463void CSG_Classifier_Supervised::Set_Threshold_Probability(double Value) { m_Threshold_Probability = Value; }
2464double CSG_Classifier_Supervised::Get_Threshold_Probability(void) { return( m_Threshold_Probability ); }
2465
2466//---------------------------------------------------------
2467void CSG_Classifier_Supervised::Set_Probability_Relative (bool Value) { m_Probability_Relative = Value; }
2468bool CSG_Classifier_Supervised::Get_Probability_Relative (void) { return( m_Probability_Relative ); }
2469
2470//---------------------------------------------------------
2471void CSG_Classifier_Supervised::Set_WTA(int Method, bool bOn)
2472{
2473 if( Method >= 0 && Method < SG_CLASSIFY_SUPERVISED_WTA )
2474 {
2475 m_bWTA[Method] = bOn;
2476 }
2477}
2478
2480{
2481 return( Method >= 0 && Method < SG_CLASSIFY_SUPERVISED_WTA ? m_bWTA[Method] : false );
2482}
2483
2484
2486// //
2488
2489//---------------------------------------------------------
2490#include "saga_api.h"
2491
2492//---------------------------------------------------------
2494{
2495 int nFeatures = m_nFeatures; Destroy(); m_nFeatures = nFeatures;
2496
2497 //-----------------------------------------------------
2498 CSG_MetaData Data;
2499
2500 if( !Data.Load(File) || !Data.Cmp_Name("supervised_classifier") || SG_Compare_Version(Data.Get_Property("saga-version"), "2.1.4") < 0 )
2501 {
2502 return( false );
2503 }
2504
2505 if( !Data("classes") || !Data("features") || !Data["features"]("count") || Data["features"]["count"].Get_Content().asInt() != m_nFeatures || m_nFeatures == 0 )
2506 {
2507 return( false );
2508 }
2509
2510 if( Data["features"]("info") )
2511 {
2512 m_Info = Data["features"]["info"].Get_Content();
2513 }
2514
2515 //-----------------------------------------------------
2516 CSG_MetaData &Classes = *Data.Get_Child("CLASSES");
2517
2518 for(int i=0; i<Classes.Get_Children_Count(); i++)
2519 {
2520 if( Classes[i].Cmp_Name("class") && Classes[i].Get_Child("id") )
2521 {
2522 bool bAdd = true;
2523
2524 CClass *pClass = new CClass(Classes[i]["id"].Get_Content());
2525
2526 if( !pClass->m_Cov .from_String(Classes[i]["cov" ].Get_Content()) || pClass->m_Cov .Get_NX() != m_nFeatures || !pClass->m_Cov.is_Square() ) { bAdd = false; }
2527 if( !pClass->m_Mean.from_String(Classes[i]["mean"].Get_Content()) || pClass->m_Mean.Get_N () != m_nFeatures ) { bAdd = false; }
2528 if( !pClass->m_Min .from_String(Classes[i]["min" ].Get_Content()) || pClass->m_Min .Get_N () != m_nFeatures ) { bAdd = false; }
2529 if( !pClass->m_Max .from_String(Classes[i]["max" ].Get_Content()) || pClass->m_Max .Get_N () != m_nFeatures ) { bAdd = false; }
2530
2531 //---------------------------------------------
2532 if( !bAdd )
2533 {
2534 delete(pClass);
2535 }
2536 else
2537 {
2538 m_pClasses = (CClass **)SG_Realloc(m_pClasses, ((size_t)m_nClasses + 1) * sizeof(CClass *));
2539 m_pClasses[m_nClasses++] = pClass;
2540
2541 pClass->m_Cov_Det = pClass->m_Cov.Get_Determinant();
2542 pClass->m_Cov_Inv = pClass->m_Cov.Get_Inverse();
2543
2544 pClass->m_Mean_Spectral = CSG_Simple_Statistics(pClass->m_Mean).Get_Mean();
2545 }
2546 }
2547 }
2548
2549 return( m_nClasses > 0 );
2550}
2551
2552//---------------------------------------------------------
2553bool CSG_Classifier_Supervised::Save(const CSG_String &File, const SG_Char *Feature_Info)
2554{
2555 if( m_nFeatures < 1 || m_nClasses < 1 || File.is_Empty() )
2556 {
2557 return( false );
2558 }
2559
2560 CSG_MetaData Data;
2561
2562 Data.Set_Name ("supervised_classifier");
2563 Data.Add_Property("saga-version", SAGA_VERSION);
2564
2565 CSG_MetaData &Features = *Data.Add_Child("features");
2566
2567 Features.Add_Child("count", m_nFeatures);
2568
2569 if( Feature_Info && *Feature_Info )
2570 {
2571 Features.Add_Child("info", Feature_Info);
2572 }
2573
2574 CSG_MetaData &Classes = *Data.Add_Child("classes");
2575
2576 Classes.Add_Property("count", m_nClasses);
2577
2578 for(int i=0; i<m_nClasses; i++)
2579 {
2580 CSG_MetaData &Class = *Classes.Add_Child("class");
2581
2582 CClass *pClass = m_pClasses[i];
2583
2584 Class.Add_Child("id" , pClass->m_ID );
2585 Class.Add_Child("mean", pClass->m_Mean.to_String());
2586 Class.Add_Child("min" , pClass->m_Min .to_String());
2587 Class.Add_Child("max" , pClass->m_Max .to_String());
2588 Class.Add_Child("cov" , pClass->m_Cov .to_String());
2589 }
2590
2591 return( Data.Save(File) );
2592}
2593
2594
2596// //
2598
2599//---------------------------------------------------------
2601{
2602 CSG_String s;
2603
2604 if( m_nFeatures > 0 && m_nClasses > 0 )
2605 {
2606 s += "\n";
2607
2608 for(int iClass=0; iClass<m_nClasses; iClass++)
2609 {
2610 CClass *pClass = m_pClasses[iClass];
2611
2612 s += "\n____\n" + pClass->m_ID + "\nFeature\tMean\tMin\tMax\tStdDev";
2613
2614 for(int i=0; i<m_nFeatures; i++)
2615 {
2616 s += CSG_String::Format("\n%3d.", i + 1);
2617 s += "\t" + SG_Get_String(pClass->m_Mean[i]);
2618 s += "\t" + SG_Get_String(pClass->m_Min [i]);
2619 s += "\t" + SG_Get_String(pClass->m_Max [i]);
2620 s += "\t" + SG_Get_String(sqrt(pClass->m_Cov[i][i]));
2621 }
2622
2623 s += "\n";
2624 }
2625 }
2626
2627 return( s );
2628}
2629
2630
2632// //
2634
2635//---------------------------------------------------------
2636bool CSG_Classifier_Supervised::Add_Class(const CSG_String &Class_ID, const CSG_Vector &Mean, const CSG_Vector &Min, const CSG_Vector &Max, const CSG_Matrix &Cov)
2637{
2638 if( m_nFeatures < 1 || Mean.Get_N() != m_nFeatures || Min.Get_N() != m_nFeatures || Max.Get_N() != m_nFeatures || Cov.Get_NCols() != m_nFeatures || Cov.Get_NRows() != m_nFeatures )
2639 {
2640 return( false );
2641 }
2642
2643 CClass *pClass, **pClasses = (CClass **)SG_Realloc(m_pClasses, ((size_t)m_nClasses + 1) * sizeof(CClass *));
2644
2645 if( pClasses )
2646 {
2647 m_pClasses = pClasses;
2648
2649 m_pClasses[m_nClasses++] = pClass = new CClass(Class_ID);
2650
2651 pClass->m_ID = Class_ID;
2652 pClass->m_Mean = Mean;
2653 pClass->m_Min = Min;
2654 pClass->m_Max = Max;
2655 pClass->m_Cov = Cov;
2656 pClass->m_Cov_Inv = Cov.Get_Inverse();
2657 pClass->m_Cov_Det = Cov.Get_Determinant();
2658
2659 pClass->m_Mean_Spectral = CSG_Simple_Statistics(Mean).Get_Mean();
2660
2661 return( true );
2662 }
2663
2664 return( false );
2665}
2666
2667
2669// //
2671
2672//---------------------------------------------------------
2674{
2675 for(int i=0; i<m_nClasses; i++)
2676 {
2677 m_pClasses[i]->m_Samples.Destroy();
2678 }
2679
2680 return( true );
2681}
2682
2683//---------------------------------------------------------
2685{
2686 if( m_nFeatures > 0 && m_nFeatures == Features.Get_N() )
2687 {
2688 int iClass = Get_Class(Class_ID);
2689
2690 if( iClass < 0 )
2691 {
2692 CClass **pClasses = (CClass **)SG_Realloc(m_pClasses, ((size_t)m_nClasses + 1) * sizeof(CClass *));
2693
2694 if( pClasses )
2695 {
2696 m_pClasses = pClasses;
2697
2698 m_pClasses[iClass = m_nClasses++] = new CClass(Class_ID);
2699 }
2700 }
2701
2702 if( iClass >= 0 )
2703 {
2704 return( m_pClasses[iClass]->m_Samples.Add_Row(Features) );
2705 }
2706 }
2707
2708 return( false );
2709}
2710
2711//---------------------------------------------------------
2712bool CSG_Classifier_Supervised::Train(bool bClear_Samples)
2713{
2714 if( m_nFeatures < 1 || m_nClasses < 1 )
2715 {
2716 return( false );
2717 }
2718
2719 for(int iClass=0; iClass<m_nClasses; iClass++)
2720 {
2721 if( !m_pClasses[iClass]->Train() )
2722 {
2723 return( false );
2724 }
2725 }
2726
2727 if( bClear_Samples )
2728 {
2730 }
2731
2732 return( true );
2733}
2734
2735
2737// //
2739
2740//---------------------------------------------------------
2741bool CSG_Classifier_Supervised::CClass::Train(void)
2742{
2743 if( m_Samples.Get_NCols() < 1 || m_Samples.Get_NRows() < 1 )
2744 {
2745 return( false );
2746 }
2747
2748 //-----------------------------------------------------
2749 m_Mean.Create(m_Samples.Get_NCols());
2750 m_Min .Create(m_Samples.Get_NCols());
2751 m_Max .Create(m_Samples.Get_NCols());
2752
2753 for(int iFeature=0; iFeature<m_Samples.Get_NCols(); iFeature++)
2754 {
2755 CSG_Simple_Statistics s;
2756
2757 for(int iSample=0; iSample<m_Samples.Get_NRows(); iSample++)
2758 {
2759 s += m_Samples[iSample][iFeature];
2760 }
2761
2762 m_Mean[iFeature] = s.Get_Mean ();
2763 m_Min [iFeature] = s.Get_Minimum();
2764 m_Max [iFeature] = s.Get_Maximum();
2765 }
2766
2767 //-----------------------------------------------------
2768 m_Cov.Create(m_Samples.Get_NCols(), m_Samples.Get_NCols());
2769
2770 for(int iFeature=0; iFeature<m_Samples.Get_NCols(); iFeature++)
2771 {
2772 for(int jFeature=iFeature; jFeature<m_Samples.Get_NCols(); jFeature++)
2773 {
2774 double cov = 0.;
2775
2776 for(int iSample=0; iSample<m_Samples.Get_NRows(); iSample++)
2777 {
2778 cov += (m_Samples[iSample][iFeature] - m_Mean[iFeature]) * (m_Samples[iSample][jFeature] - m_Mean[jFeature]);
2779 }
2780
2781 if( m_Samples.Get_NRows() > 1 )
2782 {
2783 cov /= m_Samples.Get_NRows() - 1;
2784 }
2785
2786 m_Cov[iFeature][jFeature] = m_Cov[jFeature][iFeature] = cov;
2787 }
2788 }
2789
2790 m_Cov_Inv = m_Cov.Get_Inverse ();
2791 m_Cov_Det = m_Cov.Get_Determinant();
2792
2793 m_Mean_Spectral = CSG_Simple_Statistics(m_Mean).Get_Mean();
2794
2795 //-----------------------------------------------------
2796 return( true );
2797}
2798
2799
2801// //
2803
2804//---------------------------------------------------------
2806{
2807 if( m_nFeatures > 0 )
2808 {
2809 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2810 {
2811 if( !Get_Class_ID(iClass).Cmp(Class_ID) )
2812 {
2813 return( iClass );
2814 }
2815 }
2816 }
2817
2818 return( -1 );
2819}
2820
2821//---------------------------------------------------------
2822bool CSG_Classifier_Supervised::Get_Class(const CSG_Vector &Features, int &Class, double &Quality, int Method)
2823{
2824 Class = -1; Quality = 0.;
2825
2826 if( Get_Feature_Count() == Features.Get_N() )
2827 {
2828 switch( Method )
2829 {
2830 case SG_CLASSIFY_SUPERVISED_BinaryEncoding : _Get_Binary_Encoding (Features, Class, Quality); break;
2831 case SG_CLASSIFY_SUPERVISED_ParallelEpiped : _Get_Parallel_Epiped (Features, Class, Quality); break;
2832 case SG_CLASSIFY_SUPERVISED_MinimumDistance : _Get_Minimum_Distance (Features, Class, Quality); break;
2833 case SG_CLASSIFY_SUPERVISED_Mahalonobis : _Get_Mahalanobis_Distance (Features, Class, Quality); break;
2834 case SG_CLASSIFY_SUPERVISED_MaximumLikelihood: _Get_Maximum_Likelihood (Features, Class, Quality); break;
2835 case SG_CLASSIFY_SUPERVISED_SAM : _Get_Spectral_Angle_Mapping(Features, Class, Quality); break;
2836 case SG_CLASSIFY_SUPERVISED_SID : _Get_Spectral_Divergence (Features, Class, Quality); break;
2837 case SG_CLASSIFY_SUPERVISED_WTA : _Get_Winner_Takes_All (Features, Class, Quality); break;
2838 }
2839
2840 return( Class >= 0 );
2841 }
2842
2843 return( false );
2844}
2845
2846
2848// //
2850
2851//---------------------------------------------------------
2853{
2854 switch( Method )
2855 {
2856 case SG_CLASSIFY_SUPERVISED_BinaryEncoding : return( _TL("Binary Encoding") );
2857 case SG_CLASSIFY_SUPERVISED_ParallelEpiped : return( _TL("Parallelepiped") );
2858 case SG_CLASSIFY_SUPERVISED_MinimumDistance : return( _TL("Minimum Distance") );
2859 case SG_CLASSIFY_SUPERVISED_Mahalonobis : return( _TL("Mahalanobis Distance") );
2860 case SG_CLASSIFY_SUPERVISED_MaximumLikelihood: return( _TL("Maximum Likelihood") );
2861 case SG_CLASSIFY_SUPERVISED_SAM : return( _TL("Spectral Angle Mapping") );
2862 case SG_CLASSIFY_SUPERVISED_SID : return( _TL("Spectral Information Divergence") );
2863 case SG_CLASSIFY_SUPERVISED_SVM : return( _TL("Support Vector Machine") );
2864 case SG_CLASSIFY_SUPERVISED_WTA : return( _TL("Winner Takes All") );
2865 }
2866
2867 return( SG_T("") );
2868}
2869
2870//---------------------------------------------------------
2872{
2873 switch( Method )
2874 {
2875 case SG_CLASSIFY_SUPERVISED_BinaryEncoding : return( _TL("Difference") );
2876 case SG_CLASSIFY_SUPERVISED_ParallelEpiped : return( _TL("Memberships") );
2877 case SG_CLASSIFY_SUPERVISED_MinimumDistance : return( _TL("Distance") );
2878 case SG_CLASSIFY_SUPERVISED_Mahalonobis : return( _TL("Distance") );
2879 case SG_CLASSIFY_SUPERVISED_MaximumLikelihood: return( _TL("Proximity") );
2880 case SG_CLASSIFY_SUPERVISED_SAM : return( _TL("Angle") );
2881 case SG_CLASSIFY_SUPERVISED_SID : return( _TL("Divergence") );
2882 case SG_CLASSIFY_SUPERVISED_SVM : return( _TL("") );
2883 case SG_CLASSIFY_SUPERVISED_WTA : return( _TL("Votes") );
2884 }
2885
2886 return( SG_T("") );
2887}
2888
2889
2891// //
2893
2894//---------------------------------------------------------
2895// Mazer, A. S., Martin, M., Lee, M., and Solomon, J. E. (1988):
2896// Image Processing Software for Imaging Spectrometry Analysis.
2897// Remote Sensing of Environment, v. 24, no. 1, p. 201-210.
2898//
2899void CSG_Classifier_Supervised::_Get_Binary_Encoding(const CSG_Vector &Features, int &Class, double &Quality)
2900{
2901 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2902 {
2903 CClass *pClass = m_pClasses[iClass];
2904
2905 double Mean_Spectral = CSG_Simple_Statistics(Features).Get_Mean();
2906
2907 int d = 0;
2908
2909 for(int iFeature=0; iFeature<Get_Feature_Count(); iFeature++)
2910 {
2911 d += (Features(iFeature) < Mean_Spectral) == (pClass->m_Mean[iFeature] < pClass->m_Mean_Spectral) ? 0 : 1;
2912
2913 if( iFeature == 0 ) // spectral slopes
2914 {
2915 d += (Features[iFeature ] < Features[iFeature + 1]) == (pClass->m_Mean[iFeature ] < pClass->m_Mean[iFeature + 1]) ? 0 : 1;
2916 }
2917 else if( iFeature == Get_Feature_Count() - 1 )
2918 {
2919 d += (Features[iFeature - 1] < Features[iFeature ]) == (pClass->m_Mean[iFeature - 1] < pClass->m_Mean[iFeature ]) ? 0 : 1;
2920 }
2921 else
2922 {
2923 d += (Features[iFeature - 1] < Features[iFeature + 1]) == (pClass->m_Mean[iFeature - 1] < pClass->m_Mean[iFeature + 1]) ? 0 : 1;
2924 }
2925 }
2926
2927 if( Class < 0 || Quality > d ) // find the minimum 'Hamming' distance
2928 {
2929 Class = iClass; Quality = d;
2930 }
2931 }
2932}
2933
2934//---------------------------------------------------------
2935void CSG_Classifier_Supervised::_Get_Parallel_Epiped(const CSG_Vector &Features, int &Class, double &Quality)
2936{
2937 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2938 {
2939 CClass *pClass = m_pClasses[iClass];
2940
2941 bool bMember = true;
2942
2943 for(int iFeature=0; bMember && iFeature<Get_Feature_Count(); iFeature++)
2944 {
2945 bMember = pClass->m_Min[iFeature] <= Features[iFeature] && Features[iFeature] <= pClass->m_Max[iFeature];
2946 }
2947
2948 if( bMember )
2949 {
2950 Class = iClass; Quality++;
2951 }
2952 }
2953}
2954
2955//---------------------------------------------------------
2956void CSG_Classifier_Supervised::_Get_Minimum_Distance(const CSG_Vector &Features, int &Class, double &Quality)
2957{
2958 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2959 {
2960 CClass *pClass = m_pClasses[iClass];
2961
2962 double Distance = (Features - pClass->m_Mean).Get_Length();
2963
2964 if( Class < 0 || Quality > Distance )
2965 {
2966 Class = iClass; Quality = Distance;
2967 }
2968 }
2969
2970 if( m_Threshold_Distance > 0. && Quality > m_Threshold_Distance )
2971 {
2972 Class = -1;
2973 }
2974}
2975
2976//---------------------------------------------------------
2977void CSG_Classifier_Supervised::_Get_Mahalanobis_Distance(const CSG_Vector &Features, int &Class, double &Quality)
2978{
2979 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2980 {
2981 CClass *pClass = m_pClasses[iClass];
2982
2983 CSG_Vector D = Features - pClass->m_Mean;
2984
2985 double Distance = D * (pClass->m_Cov_Inv * D);
2986
2987 if( Class < 0 || Quality > Distance )
2988 {
2989 Class = iClass; Quality = Distance;
2990 }
2991 }
2992
2993 if( m_Threshold_Distance > 0. && Quality > m_Threshold_Distance )
2994 {
2995 Class = -1;
2996 }
2997}
2998
2999//---------------------------------------------------------
3000void CSG_Classifier_Supervised::_Get_Maximum_Likelihood(const CSG_Vector &Features, int &Class, double &Quality)
3001{
3002 double dSum = 0.;
3003
3004 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
3005 {
3006 CClass *pClass = m_pClasses[iClass];
3007
3008 CSG_Vector D = Features - pClass->m_Mean;
3009
3010 double Distance = D * (pClass->m_Cov_Inv * D);
3011
3012 double Probability = pow(2. * M_PI, -0.5 * m_nFeatures) * pow(pClass->m_Cov_Det, -0.5) * exp(-0.5 * Distance);
3013 // double Probability = -log(pClass->m_Cov_Det) - Distance;
3014
3015 dSum += Probability;
3016
3017 if( Class < 0 || Quality < Probability )
3018 {
3019 Class = iClass; Quality = Probability;
3020 }
3021 }
3022
3023 if( Class >= 0 )
3024 {
3025 if( m_Probability_Relative )
3026 {
3027 Quality = 100. * Quality / dSum;
3028 }
3029
3030 if( m_Threshold_Probability > 0. && Quality < m_Threshold_Probability )
3031 {
3032 Class = -1;
3033 }
3034 }
3035}
3036
3037//---------------------------------------------------------
3038void CSG_Classifier_Supervised::_Get_Spectral_Angle_Mapping(const CSG_Vector &Features, int &Class, double &Quality)
3039{
3040 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
3041 {
3042 CClass *pClass = m_pClasses[iClass];
3043
3044 double Angle = Features.Get_Angle(pClass->m_Mean);
3045
3046 if( Class < 0 || Quality > Angle )
3047 {
3048 Class = iClass; Quality = Angle;
3049 }
3050 }
3051
3052 Quality *= M_RAD_TO_DEG;
3053
3054 if( m_Threshold_Angle > 0. && Quality > m_Threshold_Angle )
3055 {
3056 Class = -1;
3057 }
3058}
3059
3060//---------------------------------------------------------
3061void CSG_Classifier_Supervised::_Get_Spectral_Divergence(const CSG_Vector &Features, int &Class, double &Quality)
3062{
3063}
3064
3065//---------------------------------------------------------
3066void CSG_Classifier_Supervised::_Get_Winner_Takes_All(const CSG_Vector &Features, int &Class, double &Quality)
3067{
3068 int *Votes = (int *)SG_Calloc(Get_Class_Count(), sizeof(int));
3069
3070 for(int iMethod=0; iMethod<SG_CLASSIFY_SUPERVISED_WTA; iMethod++)
3071 {
3072 int iClass; double iQuality;
3073
3074 if( m_bWTA[iMethod] && Get_Class(Features, iClass, iQuality, iMethod) && ++Votes[iClass] > Quality )
3075 {
3076 Class = iClass; Quality = Votes[iClass];
3077 }
3078 }
3079
3080 SG_Free(Votes);
3081}
3082
3083
3085// //
3086// //
3087// //
3089
3090//---------------------------------------------------------
3091// source: http://psydok.sulb.uni-saarland.de/volltexte/2004/268/html/
3092
3093//---------------------------------------------------------
3095{ // Hill's approx. to cumulative t-dist, Commun.A.C.M. 13,617-619.
3096 // See: J.H.Maindonald, Computational Statistics, p.295.
3097 // Calculates p given t and tail type.
3098
3099 if( !T || !df || df < 1. )
3100 {
3101 return( -1. );
3102 }
3103
3104 return( _Change_Tail_Type(Get_T_P(T, df), TESTDIST_TYPE_TwoTail, Type, T < 0.) );
3105}
3106
3107//---------------------------------------------------------
3109{ // Keith Dear & Robert Brennan.
3110 // Returns an accurate t to tol sig. fig.'s given p & df.
3111
3112 if( p <= 0. || p >= 1. || df < 1 )
3113 {
3114 return( -1. );
3115 }
3116
3117 bool bNegative = (Type == TESTDIST_TYPE_Left && p < 0.5) || (Type == TESTDIST_TYPE_Right && p > 0.5);
3118 double t = 0, p0, p1, diff = 1.;
3119
3120 p0 = p1 = _Change_Tail_Type(p, Type, TESTDIST_TYPE_TwoTail, bNegative);
3121
3122 while( fabs(diff) > 0.0001 )
3123 {
3124 t = Get_T_Inv(p1, df); // initial rough value
3125 diff = Get_T_P(t, df) - p0; // compare result with forward fn
3126 p1 = p1 - diff; // small adjustment to p1
3127 }
3128
3129 return( bNegative ? -t : t );
3130}
3131
3132//---------------------------------------------------------
3133double CSG_Test_Distribution::_Change_Tail_Type(double p, TSG_Test_Distribution_Type from, TSG_Test_Distribution_Type to, bool bNegative)
3134{
3135 if( from != to )
3136 {
3137 switch( from ) // convert any tail type to 'left'
3138 {
3139 case TESTDIST_TYPE_Left : break;
3140 case TESTDIST_TYPE_Right : p = 1. - p; break;
3141 case TESTDIST_TYPE_Middle : p = p / 2. + 0.5; if( bNegative ) p = 1. - p; break;
3142 case TESTDIST_TYPE_TwoTail: p = 1. - p / 2. ; if( bNegative ) p = 1. - p; break;
3143 // case TESTDIST_TYPE_Half : p = p + 0.5 ; if( bNegative ) p = 1. - p; break;
3144 }
3145
3146 switch( to ) // convert p from tail type 'left' to any other
3147 {
3148 case TESTDIST_TYPE_Left : break;
3149 case TESTDIST_TYPE_Right : p = 1. - p; break;
3150 case TESTDIST_TYPE_Middle : if( bNegative ) p = 1. - p; p = 2. * (1. - p); break;
3151 case TESTDIST_TYPE_TwoTail: if( bNegative ) p = 1. - p; p = 2. * p - 1. ; break;
3152 // case TESTDIST_TYPE_Half : if( bNegative ) p = 1. - p; p = p - 0.5 ; break;
3153 }
3154 }
3155
3156 return( p );
3157}
3158
3159//---------------------------------------------------------
3161{ // Returns the two-tailed standard normal probability of z
3162 const double a1 = 0.0000053830, a2 = 0.0000488906, a3 = 0.0000380036;
3163 const double a4 = 0.0032776263, a5 = 0.0211410061, a6 = 0.0498673470;
3164
3165 z = fabs(z);
3166
3167 double p = (((((a1 * z + a2) * z + a3) * z + a4) * z + a5) * z + a6) * z + 1.;
3168
3169 return( pow(p, -16) );
3170}
3171
3172//---------------------------------------------------------
3174{ // Returns z given a half-middle tail type p.
3175 const double a0 = 2.5066282, a1 = -18.6150006, a2 = 41.3911977, a3 = -25.4410605;
3176 const double b1 = -8.4735109, b2 = 23.0833674, b3 = -21.0622410, b4 = 3.1308291;
3177 const double c0 = -2.7871893, c1 = -2.2979648, c2 = 4.8501413, c3 = 2.3212128;
3178 const double d1 = 3.5438892, d2 = 1.6370678;
3179
3180 if( p > 0.42 )
3181 {
3182 double r = sqrt(-log(0.5 - p));
3183
3184 return( (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + 1.) );
3185 }
3186 else
3187 {
3188 double r = p * p;
3189
3190 return( p * (((a3 * r + a2) * r + a1) * r + a0) / ((((b4 * r + b3) * r + b2) * r + b1) * r + 1.) );
3191 }
3192}
3193
3194//---------------------------------------------------------
3195double CSG_Test_Distribution::Get_T_P(double T, int df)
3196{ // Returns two-tail probability level given t and df.
3197 return( df == 1 ? 1. - 2. * atan(fabs(T)) / M_PI
3198 : df == 2 ? 1. - fabs(T) / sqrt(T*T + 2.)
3199 : df == 3 ? 1. - 2. * (atan(fabs(T) / sqrt(3.)) + fabs(T) * sqrt(3.) / (T*T + 3.)) / M_PI
3200 : df == 4 ? 1. - fabs(T) * (1. + 2. / (T*T + 4.)) / sqrt(T*T + 4.)
3201 : Get_Norm_P(Get_T_Z(fabs(T), df))
3202 );
3203}
3204
3205//---------------------------------------------------------
3206double CSG_Test_Distribution::Get_T_Z(double T, int df)
3207{ // Converts a t value to an approximate z value w.r.t the given df
3208 // s.t. std.norm.(z) = t(z, df) at the two-tail probability level.
3209
3210 double A9, B9, T9, Z8, P7, B7, z;
3211
3212 A9 = df - 0.5;
3213 B9 = 48. * A9*A9,
3214 T9 = T*T / df;
3215 Z8 = T9 >= 0.04
3216 ? A9 * log(1. + T9)
3217 : A9 * (((1. - T9 * 0.75) * T9 / 3. - 0.5) * T9 + 1.) * T9;
3218 P7 = ((0.4 * Z8 + 3.3) * Z8 + 24.) * Z8 + 85.5;
3219 B7 = 0.8 * pow(Z8, 2.) + 100. + B9;
3220 z = (1. + (-P7 / B7 + Z8 + 3.) / B9) * sqrt(Z8);
3221
3222 return( z );
3223}
3224
3225//---------------------------------------------------------
3226double CSG_Test_Distribution::Get_T_Inv(double p, int df)
3227{ // Hill's approx. inverse t-dist.: Comm. of A.C.M Vol.13 No.10 1970 pg 620.
3228 // Calculates t given df and two-tail probability.
3229
3230 if( df == 1 )
3231 {
3232 return( cos(p * M_PI / 2.) / sin(p * M_PI / 2.) );
3233 }
3234
3235 if( df == 2 )
3236 {
3237 return( sqrt(2. / (p * (2. - p)) - 2.) );
3238 }
3239
3240 double a = 1. / (df - 0.5);
3241 double b = 48. / (a*a);
3242 double c = ((20700. * a / b - 98.) * a - 16.) * a + 96.36;
3243 double d = ((94.5 / (b + c) - 3.) / b + 1.) * sqrt(a * M_PI / 2.) * df;
3244 double x = d * p;
3245 double y = pow(x, 2. / df);
3246
3247 if( y > 0.05 + a )
3248 {
3249 x = Get_Norm_Z(0.5 * (1. - p));
3250 y = x*x;
3251
3252 if( df < 5 )
3253 {
3254 c = c + 0.3 * (df - 4.5) * (x + 0.6);
3255 }
3256
3257 c = (((0.05 * d * x - 5) * x - 7.) * x - 2.) * x + b + c;
3258 y = (((((0.4 * y + 6.3) * y + 36.) * y + 94.5) / c - y - 3.) / b + 1.) * x;
3259 y = a * y*y;
3260
3261 if( y > 0.002 )
3262 {
3263 y = exp(y) - 1.;
3264 }
3265 else
3266 {
3267 y = 0.5 * y*y + y;
3268 }
3269 }
3270 else
3271 {
3272 y = ((1. / (((df + 6.) / (df * y) - 0.089 * d - 0.822) * (df + 2.) * 3.)
3273 + 0.5 / (df + 4.)) * y - 1.) * (df + 1.) / (df + 2.) + 1. / y;
3274 }
3275
3276 return( sqrt(df * y) );
3277}
3278
3279
3281// //
3283
3284//---------------------------------------------------------
3285double CSG_Test_Distribution::Get_F_Tail_from_R2(double R2, int nPredictors, int nSamples, TSG_Test_Distribution_Type Type)
3286{
3287 double F = ((sLong)nSamples - (sLong)nPredictors - 1) * (R2 / nPredictors) / (1. - R2);
3288
3289 return( CSG_Test_Distribution::Get_F_Tail(F, nPredictors, nSamples - nPredictors - 1, Type) );
3290}
3291
3292//---------------------------------------------------------
3294{
3295 // calculates for F, dfn(ominator) and dfd(enominator) the "tail" of the F-distribution
3296
3297 double p = 1.;
3298
3299 if( F >= 0.00001 && dfn > 0 && dfd > 0 )
3300 {
3301 if( F * dfn >= dfd || F > 1. + 20. / dfn + 10. / sqrt((double)dfn) )
3302 {
3303 p = Get_Gamma(F, dfn, dfd);
3304 }
3305 else
3306 {
3307 p = 1. - Get_Gamma(1. / F, dfd, dfn);
3308 }
3309 }
3310
3311 if( p <= 0. || p >= 1. )
3312 {
3313 p = F > 1. ? 0. : F < 1. ? 1. : 0.5;
3314 }
3315
3316 return( Type == TESTDIST_TYPE_Right ? p : 1. - p );
3317}
3318
3319//---------------------------------------------------------
3320double CSG_Test_Distribution::Get_F_Inverse(double alpha, int dfn, int dfd, TSG_Test_Distribution_Type Type)
3321{
3322 if( alpha < 0. || alpha > 1. || dfd < 0 || dfn < 0 )
3323 {
3324 return( -1 );
3325 }
3326
3327 if( Type != TESTDIST_TYPE_Right )
3328 {
3329 alpha = 1. - alpha;
3330 }
3331
3332 const int ITERMAX = 100;
3333 const double EPSILON = 0.0001;
3334
3335 int i;
3336 double lo, hi, mid, p;
3337
3338 if( alpha <= 0.5 )
3339 {
3340 lo = 0.5;
3341 hi = lo;
3342
3343 for(i=0; i<ITERMAX; i++)
3344 {
3345 hi *= 2.;
3346 p = Get_F_Tail(hi, dfn, dfd);
3347
3348 if( p > alpha )
3349 {
3350 lo = hi;
3351 }
3352 else
3353 {
3354 break;
3355 }
3356 }
3357
3358 if( p > alpha )
3359 {
3360 return( hi );
3361 }
3362 }
3363 else
3364 {
3365 hi = 2;
3366 lo = hi;
3367
3368 for(i=0; i<ITERMAX; i++)
3369 {
3370 lo /= 2.;
3371 p = Get_F_Tail(lo, dfn, dfd);
3372
3373 if( p < alpha )
3374 {
3375 hi = lo;
3376 }
3377 else
3378 {
3379 break;
3380 }
3381 }
3382
3383 if( p < alpha )
3384 {
3385 return( lo );
3386 }
3387 }
3388
3389 mid = (hi + lo) / 2.;
3390
3391 for(i=0; i<ITERMAX && (hi-lo)>EPSILON*mid; i++)
3392 {
3393 mid = (hi + lo) / 2.;
3394 p = Get_F_Tail(mid, dfn, dfd);
3395
3396 if( p < alpha )
3397 hi = mid;
3398 else if( p > alpha )
3399 lo = mid;
3400 else
3401 break;
3402 }
3403
3404 return( mid );
3405}
3406
3407//---------------------------------------------------------
3408double CSG_Test_Distribution::Get_Gamma(double F, double dfn, double dfd)
3409{
3410 // calculates for F, dfn(ominator) and dfd(enominator) the incomplete Gamma-function
3411
3412 const double EXPMIN = -30.;
3413 const double SMALL = 0.00000000001;
3414
3415 double x, c, er, s, n, t1, t;
3416
3417 dfn /= 2.;
3418 dfd /= 2.;
3419
3420 x = dfd / (dfd + dfn * F);
3421 c = Get_Log_Gamma(dfn + dfd) - Get_Log_Gamma(dfn) - Get_Log_Gamma(dfd + 1.) + dfd * log(x) + dfn * log(1. - x);
3422
3423 if( c < EXPMIN )
3424 {
3425 return( -1. );
3426 }
3427
3428 dfn += dfd;
3429 dfd += 1.;
3430 c = exp(c);
3431 er = SMALL / c;
3432 t = dfn * x / dfd;
3433 t1 = 0.;
3434 s = t + 1.;
3435 n = 0;
3436
3437 while( t > er || t > t1 )
3438 {
3439 n += 1;
3440 t1 = t;
3441 t *= ((dfn + n) * x / (dfd + n));
3442 s += t;
3443 }
3444
3445 return( s * c );
3446}
3447
3448//---------------------------------------------------------
3449double CSG_Test_Distribution::Get_Log_Gamma(double a)
3450{
3451 // calculates the logarithm of the Gamma-function
3452
3453 const int ARGMIN = 6;
3454
3455 const double HL2PI = 0.91893853320467275; // = log(2. * M_PI) / 2.
3456
3457 int n = (int)floor(ARGMIN - a + 0.0001);
3458
3459 if( n > 0 )
3460 {
3461 a += n;
3462 }
3463
3464 double g;
3465
3466 g = 1. / (a*a);
3467 g = (1. - g * (1. / 30. - g * (1. / 105. - g * (1. / 140. - g / 99.)))) / (12. * a);
3468 g = g + ((a - 0.5) * log(a) - a + HL2PI);
3469
3470 for(int i=0; i<n; i++)
3471 {
3472 a = a - 1.;
3473 g = g - log(a);
3474 }
3475
3476 return( g );
3477}
3478
3479
3481// //
3482// //
3483// //
3485
3486//---------------------------------------------------------
3487CSG_Matrix SG_Get_Correlation_Matrix (const CSG_Matrix &Values, bool bCovariances)
3488{
3489 int nVariables = Values.Get_NX();
3490 int nSamples = Values.Get_NY();
3491
3492 //-----------------------------------------------------
3493 int i, j, k;
3495 CSG_Matrix C;
3496
3497 C.Create(nVariables, nVariables);
3498
3499 //-----------------------------------------------------
3500 S = new CSG_Simple_Statistics[nVariables];
3501
3502 for(j=0; j<nVariables; j++)
3503 {
3504 for(i=0; i<nSamples; i++)
3505 {
3506 S[j] += Values[i][j];
3507 }
3508 }
3509
3510 //-----------------------------------------------------
3511 for(k=0; k<nVariables; k++)
3512 {
3513 for(j=k; j<nVariables; j++)
3514 {
3515 double cov = 0.;
3516
3517 for(i=0; i<nSamples; i++)
3518 {
3519 cov += (Values[i][j] - S[j].Get_Mean()) * (Values[i][k] - S[k].Get_Mean());
3520 }
3521
3522 cov /= nSamples;
3523
3524 if( !bCovariances )
3525 {
3526 cov /= (S[j].Get_StdDev() * S[k].Get_StdDev());
3527 }
3528
3529 C[j][k] = C[k][j] = cov;
3530 }
3531 }
3532
3533 //-----------------------------------------------------
3534 delete[](S);
3535
3536 return( C );
3537}
3538
3539
3541// //
3542// //
3543// //
3545
3546//---------------------------------------------------------
bool SG_UI_Process_Get_Okay(bool bBlink)
void SG_UI_Process_Set_Text(const CSG_String &Text)
SAGA_API_DLL_EXPORT void * SG_Malloc(size_t size)
signed long long sLong
Definition api_core.h:158
#define SG_T(s)
Definition api_core.h:537
SAGA_API_DLL_EXPORT void SG_Free(void *memblock)
SAGA_API_DLL_EXPORT void * SG_Realloc(void *memblock, size_t size)
SAGA_API_DLL_EXPORT void * SG_Calloc(size_t num, size_t size)
#define SG_FREE_SAFE(PTR)
Definition api_core.h:205
TSG_Data_Type
Definition api_core.h:1036
@ SG_DATATYPE_ULong
Definition api_core.h:1044
#define SG_Char
Definition api_core.h:536
#define _TL(s)
Definition api_core.h:1610
SAGA_API_DLL_EXPORT CSG_String SG_Get_String(double Value, int Precision=-99)
sLong Get_Size(void) const
Definition api_core.h:327
void * Get_Entry(sLong Index) const
Returns a pointer to the memory address of the requested variable. You have to type cast and derefere...
Definition api_core.h:331
void Create(TSG_Data_Type Type=SG_DATATYPE_String)
virtual ~CSG_Category_Statistics(void)
int Add_Value(int Value)
int Get_Category(int Value) const
int asInt(int iCategory) const
double asDouble(int iCategory) const
int Get_Count(void) const
CSG_Category_Statistics(TSG_Data_Type Type=SG_DATATYPE_String)
CSG_String asString(int iCategory) const
TSG_Data_Type Get_Category_Type(void) const
void Set_Probability_Relative(bool Value)
virtual ~CSG_Classifier_Supervised(void)
void Set_Threshold_Distance(double Value)
double Get_Threshold_Probability(void)
bool Get_WTA(int Method)
void Create(int nFeatures)
bool Add_Class(const CSG_String &Class_ID, const CSG_Vector &Mean, const CSG_Vector &Min, const CSG_Vector &Max, const CSG_Matrix &Cov)
void Set_Threshold_Probability(double Value)
bool Get_Probability_Relative(void)
bool Train(bool bClr_Samples=false)
const CSG_String & Get_Class_ID(int iClass)
Definition mat_tools.h:1262
void Set_Threshold_Angle(double Value)
bool Train_Add_Sample(const CSG_String &Class_ID, const CSG_Vector &Features)
double Get_Threshold_Distance(void)
static CSG_String Get_Name_of_Quality(int Method)
int Get_Class(const CSG_String &Class_ID)
void Set_WTA(int Method, bool bOn)
bool Save(const CSG_String &File, const SG_Char *Feature_Info=NULL)
static CSG_String Get_Name_of_Method(int Method)
bool Load(const CSG_String &File)
double Get_Threshold_Angle(void)
bool Create(int nFeatures)
bool Execute(int Method, int nClusters, int nMaxIterations=0, int Initialization=0)
bool Set_Feature(sLong iElement, int iFeature, double Value)
sLong Get_nElements(void) const
Definition mat_tools.h:1181
virtual ~CSG_Cluster_Analysis(void)
virtual bool is_NoData(int x, int y) const
Definition grid.h:727
virtual double asDouble(sLong i, bool bScaled=true) const
Definition grid.h:793
sLong Get_NCells(void) const
Definition grid.h:565
virtual bool is_NoData(int x, int y, int z) const
Definition grids.h:373
sLong Get_NCells(void) const
Definition grids.h:191
virtual double asDouble(sLong i, bool bScaled=true) const
Definition grids.h:402
CSG_Histogram & operator=(const CSG_Histogram &Histogram)
bool Scale_Element_Count(double Scale)
virtual ~CSG_Histogram(void)
double Get_Center(int i) const
Definition mat_tools.h:1065
void Add_Value(double Value)
double Get_Value(double i) const
Definition mat_tools.h:1060
size_t Get_Class_Count(void) const
Definition mat_tools.h:1049
size_t Get_Cumulative(int i) const
Definition mat_tools.h:1057
bool Create(const CSG_Histogram &Histogram)
bool Destroy(void)
double Get_Quantile(double Quantile) const
double Get_Percentile_Value(double Value) const
double Get_Percentile(double Percentile) const
bool Add_Histogram(const CSG_Histogram &Histogram)
double Get_Break(int i) const
Definition mat_tools.h:1062
bool Update(void)
size_t Get_Element_Count(void) const
Definition mat_tools.h:1051
double Get_Quantile_Value(double Value) const
sLong Get_NRows(void) const
Definition mat_tools.h:525
bool from_String(const CSG_String &String)
int Get_NX(void) const
Definition mat_tools.h:522
sLong Get_NCols(void) const
Definition mat_tools.h:524
double Get_Determinant(void) const
int Get_NY(void) const
Definition mat_tools.h:523
bool is_Square(void) const
Definition mat_tools.h:538
CSG_String to_String(int Width=-1, int Precision=-1, bool bScientific=false, const SG_Char *Separator=NULL) const
CSG_Matrix Get_Inverse(bool bSilent=true, int nSubSquare=0) const
bool Cmp_Name(const CSG_String &String, bool bNoCase=true) const
Definition metadata.cpp:461
bool Save(const CSG_String &File, const SG_Char *Extension=NULL) const
Definition metadata.cpp:879
int Get_Children_Count(void) const
Definition metadata.h:148
CSG_MetaData * Get_Child(int Index) const
Definition metadata.h:149
const CSG_String & Get_Content(void) const
Definition metadata.h:133
const SG_Char * Get_Property(int Index) const
Definition metadata.h:181
void Set_Name(const CSG_String &Name)
Definition metadata.h:130
CSG_MetaData * Add_Child(void)
Definition metadata.cpp:166
bool Load(const CSG_String &File, const SG_Char *Extension=NULL)
Definition metadata.cpp:786
bool Add_Property(const CSG_String &Name, const CSG_String &Value)
Definition metadata.cpp:559
virtual ~CSG_Natural_Breaks(void)
bool Create(class CSG_Table *pTable, int Field, int nClasses, int Histogram=0)
int Get_Count(void) const
Definition mat_tools.h:1128
CSG_Random(void)
static double Get_Gaussian(double mean, double stddev)
static void Initialize(void)
static double Get_Uniform(void)
sLong Get_IndexOfMinimum(void)
void Add(const CSG_Simple_Statistics &Statistics)
double Get_Percentile(double Percentile)
double Get_Median(void)
double Get_Value(sLong i) const
Definition mat_tools.h:777
double Get_Gini(void)
sLong Get_nValues_Above(double Threshold, bool bEquals=false)
sLong Get_nValues_Below(double Threshold, bool bEquals=false)
int is_Evaluated(void) const
Definition mat_tools.h:741
double Get_Maximum(void)
Definition mat_tools.h:749
double * Get_Values(void) const
Definition mat_tools.h:776
double Get_Mean(void)
Definition mat_tools.h:753
double Get_Sum(void)
Definition mat_tools.h:751
double Get_SkewnessPearson(void)
sLong Get_IndexOfMaximum(void)
bool Set_Count(sLong Count)
void _Evaluate(int Level=1)
bool Create(bool bHoldValues=false)
double Get_StdDev(void)
Definition mat_tools.h:755
sLong Get_Count(void) const
Definition mat_tools.h:745
void Add_Value(double Value, double Weight=1.)
double Get_Minimum(void)
Definition mat_tools.h:748
double Get_Quantile(double Quantile)
int Cmp(const CSG_String &String) const
static CSG_String Format(const char *Format,...)
bool is_Empty(void) const
bool Create(const class wxString *pString)
sLong Get_Index(void) const
Definition table.h:136
double asDouble(int Field) const
int asInt(int Field) const
const SG_Char * asString(int Field, int Decimals=-99) const
bool Add_Value(int Field, double Value)
sLong Get_Count(void) const
Definition table.h:400
virtual bool Get_Value(sLong Index, int Field, CSG_String &Value) const
Definition table.cpp:1185
bool Set_Index(CSG_Index &Index, int Field, bool bAscending=true) const
Definition table.cpp:1508
int Get_Field_Count(void) const
Definition table.h:361
CSG_Table_Record * Get_Record_byIndex(sLong Index) const
Definition table.h:407
static double Get_Norm_P(double Z)
static double Get_F_Inverse(double alpha, int dfn, int dfd, TSG_Test_Distribution_Type Type=TESTDIST_TYPE_Right)
static double Get_F_Tail_from_R2(double R2, int nPredictors, int nSamples, TSG_Test_Distribution_Type Type=TESTDIST_TYPE_Right)
static double Get_T_Tail(double T, int df, TSG_Test_Distribution_Type Type=TESTDIST_TYPE_Right)
static double Get_F_Tail(double F, int dfn, int dfd, TSG_Test_Distribution_Type Type=TESTDIST_TYPE_Right)
static double Get_Norm_Z(double P)
static double Get_T_Inverse(double alpha, int df, TSG_Test_Distribution_Type Type=TESTDIST_TYPE_Right)
void Add_Value(double Value, double Weight=1.)
virtual void Create(bool bWeights=false)
int Get_Class_Index(double Value) const
void Add_Value(const CSG_String &Value, double Weight=1.)
virtual void Create(bool bWeights=false)
int Get_Class_Index(const CSG_String &Value) const
virtual int Get_Minority(bool bWeighted=false) const
int Get_Count(void) const
Definition mat_tools.h:819
virtual int Get_Majority(bool bWeighted=false) const
CSG_String to_String(int Width=-1, int Precision=-1, bool bScientific=false, const SG_Char *Separator=NULL) const
sLong Get_Size(void) const
Definition mat_tools.h:382
int Get_N(void) const
Definition mat_tools.h:384
double Get_Angle(const CSG_Vector &Vector) const
bool from_String(const CSG_String &String)
#define C
#define EPSILON
int SG_Compare_Double(const void *a, const void *b)
double SG_Degree_To_Decimal(double Deg, double Min, double Sec)
int SG_Compare_Char_Ptr(const void *a, const void *b)
double SG_Get_Rounded(double Value, int Decimals)
Definition mat_tools.cpp:83
double SG_Get_Rounded_To_SignificantFigures(double Value, int Decimals)
int SG_Get_Digit_Count(int Number)
CSG_String SG_Get_Double_asString(double Number, int Width, int Precision, bool bScientific)
int SG_Compare_Int(const void *a, const void *b)
void SG_Decimal_To_Degree(double Value, double &Deg, double &Min, double &Sec)
double SG_Get_Square(double Value)
Definition mat_tools.cpp:70
TSG_Test_Distribution_Type
Definition mat_tools.h:1525
@ TESTDIST_TYPE_Right
Definition mat_tools.h:1527
@ TESTDIST_TYPE_Left
Definition mat_tools.h:1526
@ TESTDIST_TYPE_Middle
Definition mat_tools.h:1528
@ TESTDIST_TYPE_TwoTail
Definition mat_tools.h:1529
SAGA_API_DLL_EXPORT CSG_Matrix SG_Get_Correlation_Matrix(const CSG_Matrix &Values, bool bCovariances=false)
#define M_RAD_TO_DEG
Definition mat_tools.h:108
@ SG_CLASSIFY_SUPERVISED_ParallelEpiped
Definition mat_tools.h:1226
@ SG_CLASSIFY_SUPERVISED_SVM
Definition mat_tools.h:1233
@ SG_CLASSIFY_SUPERVISED_MinimumDistance
Definition mat_tools.h:1227
@ SG_CLASSIFY_SUPERVISED_SAM
Definition mat_tools.h:1230
@ SG_CLASSIFY_SUPERVISED_SID
Definition mat_tools.h:1232
@ SG_CLASSIFY_SUPERVISED_Mahalonobis
Definition mat_tools.h:1228
@ SG_CLASSIFY_SUPERVISED_WTA
Definition mat_tools.h:1231
@ SG_CLASSIFY_SUPERVISED_MaximumLikelihood
Definition mat_tools.h:1229
@ SG_CLASSIFY_SUPERVISED_BinaryEncoding
Definition mat_tools.h:1225
SAGA_API_DLL_EXPORT int SG_Compare_Double(const void *a, const void *b)
SAGA_API_DLL_EXPORT double SG_Get_Square(double Value)
Definition mat_tools.cpp:70
#define M_PI
Definition mat_tools.h:96
int SG_Compare_Version(const CSG_String &Version, int Major, int Minor, int Release)
Definition saga_api.cpp:82
#define SAGA_VERSION
@ TABLE_INDEX_Ascending
Definition table.h:105