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