1
2 ///////////////////////////////////////////////////////////
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 // //
50 ///////////////////////////////////////////////////////////
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
63 ///////////////////////////////////////////////////////////
64 // //
65 // //
66 // //
67 ///////////////////////////////////////////////////////////
68
69 //---------------------------------------------------------
SG_Get_Square(double Value)70 double SG_Get_Square(double Value)
71 {
72 return( Value * Value );
73 }
74
75
76 ///////////////////////////////////////////////////////////
77 // //
78 // //
79 // //
80 ///////////////////////////////////////////////////////////
81
82 //---------------------------------------------------------
SG_Get_Rounded(double Value,int Decimals)83 double 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 //---------------------------------------------------------
SG_Get_Rounded_To_SignificantFigures(double Value,int Decimals)107 double 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
137 ///////////////////////////////////////////////////////////
138 // //
139 // //
140 // //
141 ///////////////////////////////////////////////////////////
142
143 //---------------------------------------------------------
SG_Get_Digit_Count(int Number)144 int 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
152 ///////////////////////////////////////////////////////////
153 // //
154 // //
155 // //
156 ///////////////////////////////////////////////////////////
157
158 //---------------------------------------------------------
SG_Get_Double_asString(double Number,int Width,int Precision,bool bScientific)159 CSG_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
180 ///////////////////////////////////////////////////////////
181 // //
182 // //
183 // //
184 ///////////////////////////////////////////////////////////
185
186 //---------------------------------------------------------
SG_Compare_Int(const void * a,const void * b)187 int 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 //---------------------------------------------------------
SG_Compare_Double(const void * a,const void * b)199 int 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 //---------------------------------------------------------
SG_Compare_Char_Ptr(const void * a,const void * b)211 int SG_Compare_Char_Ptr(const void *a, const void *b)
212 {
213 return( strcmp((const char *)a, (const char *)b) );
214 }
215
216
217 ///////////////////////////////////////////////////////////
218 // //
219 // //
220 // //
221 ///////////////////////////////////////////////////////////
222
223 //---------------------------------------------------------
SG_Degree_To_Decimal(double Deg,double Min,double Sec)224 double 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 //---------------------------------------------------------
SG_Decimal_To_Degree(double Value,double & Deg,double & Min,double & Sec)233 void 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
247 ///////////////////////////////////////////////////////////
248 // //
249 // //
250 // //
251 ///////////////////////////////////////////////////////////
252
253 //---------------------------------------------------------
CSG_Random(void)254 CSG_Random::CSG_Random(void)
255 {
256 Initialize();
257 }
258
259 //---------------------------------------------------------
Initialize(void)260 void CSG_Random::Initialize(void)
261 {
262 Initialize((unsigned)time(NULL));
263 }
264
265 //---------------------------------------------------------
Initialize(unsigned int Value)266 void 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 //
Get_Uniform(void)274 double CSG_Random::Get_Uniform(void)
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 //
Get_Uniform(double min,double max)282 double 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 //---------------------------------------------------------
Get_Gaussian(double mean,double stddev)298 double 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
317 ///////////////////////////////////////////////////////////
318 // //
319 // //
320 // //
321 ///////////////////////////////////////////////////////////
322
323 //---------------------------------------------------------
CSG_Simple_Statistics(void)324 CSG_Simple_Statistics::CSG_Simple_Statistics(void)
325 {
326 Create(false);
327 }
328
CSG_Simple_Statistics(bool bHoldValues)329 CSG_Simple_Statistics::CSG_Simple_Statistics(bool bHoldValues)
330 {
331 Create(bHoldValues);
332 }
333
CSG_Simple_Statistics(const CSG_Simple_Statistics & Statistics)334 CSG_Simple_Statistics::CSG_Simple_Statistics(const CSG_Simple_Statistics &Statistics)
335 {
336 Create(Statistics);
337 }
338
CSG_Simple_Statistics(double Mean,double StdDev,sLong Count)339 CSG_Simple_Statistics::CSG_Simple_Statistics(double Mean, double StdDev, sLong Count)
340 {
341 Create(Mean, StdDev, Count);
342 }
343
CSG_Simple_Statistics(const CSG_Vector & Values,bool bHoldValues)344 CSG_Simple_Statistics::CSG_Simple_Statistics(const CSG_Vector &Values, bool bHoldValues)
345 {
346 Create(Values, bHoldValues);
347 }
348
349 //---------------------------------------------------------
Create(bool bHoldValues)350 bool CSG_Simple_Statistics::Create(bool bHoldValues)
351 {
352 Invalidate();
353
354 m_Values.Create(bHoldValues ? sizeof(double) : 0, 0, SG_ARRAY_GROWTH_1);
355
356 return( true );
357 }
358
Create(const CSG_Simple_Statistics & Statistics)359 bool CSG_Simple_Statistics::Create(const CSG_Simple_Statistics &Statistics)
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
Create(double Mean,double StdDev,sLong Count)386 bool 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
398 m_Sum = m_Weights * m_Mean;
399 m_Sum2 = m_Weights * (m_Mean*m_Mean + m_Variance);
400
401 m_Minimum = m_Mean - 1.5 * m_StdDev;
402 m_Maximum = m_Mean + 1.5 * m_StdDev;
403 m_Range = m_Maximum - m_Minimum;
404
405 return( true );
406 }
407
Create(const CSG_Vector & Values,bool bHoldValues)408 bool CSG_Simple_Statistics::Create(const CSG_Vector &Values, bool bHoldValues)
409 {
410 if( Create(bHoldValues) )
411 {
412 for(size_t 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 //---------------------------------------------------------
Set_Count(sLong nValues)424 bool CSG_Simple_Statistics::Set_Count(sLong nValues)
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 //---------------------------------------------------------
Invalidate(void)447 void CSG_Simple_Statistics::Invalidate(void)
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 //---------------------------------------------------------
Evaluate(void)473 bool CSG_Simple_Statistics::Evaluate(void)
474 {
475 _Evaluate();
476
477 return( is_Evaluated() > 0 );
478 }
479
480 //---------------------------------------------------------
Add(const CSG_Simple_Statistics & Statistics)481 void CSG_Simple_Statistics::Add(const CSG_Simple_Statistics &Statistics)
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 //---------------------------------------------------------
Add_Value(double Value,double Weight)527 void 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 > 0. )
543 {
544 m_Weights += 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 //---------------------------------------------------------
_Evaluate(int Level)562 void CSG_Simple_Statistics::_Evaluate(int Level)
563 {
564 if( m_bEvaluated == 0 && m_Weights > 0. )
565 {
566 m_bEvaluated = 1;
567
568 m_Range = m_Maximum - m_Minimum;
569 m_Mean = m_Sum / m_Weights;
570 m_Variance = m_Sum2 / m_Weights - m_Mean*m_Mean;
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(int 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
592 m_Kurtosis /= Get_Count();
593 m_Skewness /= Get_Count();
594 // m_Skewness *= Get_Count() / ((Get_Count() - 1) * (Get_Count() - 2));
595 }
596 }
597 }
598
599 //---------------------------------------------------------
600 /**
601 * Skewness after Pearson, i.e. the difference of mean and
602 * median divided by standard deviation.
603 * Remark: Skewness calculation is only possible, if statistics
604 * has been created with the bHoldValues flag set to true.
605 */
Get_SkewnessPearson(void)606 double CSG_Simple_Statistics::Get_SkewnessPearson(void)
607 {
608 return( Get_StdDev() != 0. ? (Get_Mean() - Get_Median()) / Get_StdDev() : 0. );
609 }
610
611 //---------------------------------------------------------
612 /**
613 * Returns the requested quantile (value between 0 and 1).
614 * The 0.5 quantile returns the median. Remark:
615 * Remark: any quantile calculation is only possible, if statistics
616 * has been created with the bHoldValues option set to true.
617 */
Get_Quantile(double Quantile)618 double CSG_Simple_Statistics::Get_Quantile(double Quantile)
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 //---------------------------------------------------------
652 /**
653 * Returns the requested percentile (i.e. the quantile requested as
654 * percentage, the 50 percentile is the medium value).
655 * Remark: any quantile calculation is only possible, if statistics
656 * has been created with the bHoldValues option set to true.
657 */
Get_Percentile(double Percentile)658 double CSG_Simple_Statistics::Get_Percentile(double Percentile)
659 {
660 return( Get_Quantile(Percentile / 100.) );
661 }
662
663 //---------------------------------------------------------
664 /**
665 * Returns the medium (i.e. thee 50 percentile or 0.5 quantile).
666 * Remark: any quantile calculation is only possible, if statistics
667 * has been created with the bHoldValues option set to true.
668 */
Get_Median(void)669 double CSG_Simple_Statistics::Get_Median(void)
670 {
671 return( Get_Quantile(0.5) );
672 }
673
674 //---------------------------------------------------------
675 /**
676 * The Gini coefficient is a measure of statistical dispersion
677 * intended to represent the income or wealth distribution of
678 * a nation's residents, and is the most commonly used measure
679 * of inequality.
680 */
Get_Gini(void)681 double CSG_Simple_Statistics::Get_Gini(void)
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(int 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 //---------------------------------------------------------
706 /**
707 * Returns the index of the minimum value in the order values have been added to the statistics.
708 * This is only supported for statistics that have been created with the bHoldValues option set to true.
709 */
Get_IndexOfMinimum(void)710 sLong CSG_Simple_Statistics::Get_IndexOfMinimum(void)
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<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 //---------------------------------------------------------
730 /**
731 * Returns the index of the maximum value in the order values have been added to the statistics.
732 * This is only supported for statistics that have been created with the bHoldValues option set to true.
733 */
Get_IndexOfMaximum(void)734 sLong CSG_Simple_Statistics::Get_IndexOfMaximum(void)
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<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 //---------------------------------------------------------
754 /**
755 * Returns the number of values greater than the threshold value.
756 * This is only supported for statistics that have been created with the bHoldValues option set to true.
757 */
Get_nValues_Above(double Threshold,bool bEquals)758 sLong 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 //---------------------------------------------------------
776 /**
777 * Returns the number of values lower than the threshold value.
778 * This is only supported for statistics that have been created with the bHoldValues option set to true.
779 */
Get_nValues_Below(double Threshold,bool bEquals)780 sLong 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
798 ///////////////////////////////////////////////////////////
799 // //
800 // //
801 // //
802 ///////////////////////////////////////////////////////////
803
804 //---------------------------------------------------------
Get_Majority(bool bWeighted) const805 int CSG_Unique_Value_Statistics::Get_Majority(bool bWeighted) const
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 //---------------------------------------------------------
Get_Minority(bool bWeighted) const833 int CSG_Unique_Value_Statistics::Get_Minority(bool bWeighted) const
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
861 ///////////////////////////////////////////////////////////
862 // //
863 ///////////////////////////////////////////////////////////
864
865 //---------------------------------------------------------
Create(bool bWeights)866 void CSG_Unique_Number_Statistics::Create(bool bWeights)
867 {
868 m_bWeights = bWeights;
869
870 m_Count.Destroy();
871 m_Value.Destroy();
872 }
873
874 //---------------------------------------------------------
Add_Value(double Value,double Weight)875 void 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 //---------------------------------------------------------
Get_Class_Index(double Value) const902 int CSG_Unique_Number_Statistics::Get_Class_Index(double Value) const
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
916 ///////////////////////////////////////////////////////////
917 // //
918 ///////////////////////////////////////////////////////////
919
920 //---------------------------------------------------------
Create(bool bWeights)921 void CSG_Unique_String_Statistics::Create(bool bWeights)
922 {
923 m_bWeights = bWeights;
924
925 m_Count.Destroy();
926 m_Value.Clear();
927 }
928
929 //---------------------------------------------------------
Add_Value(const CSG_String & Value,double Weight)930 void CSG_Unique_String_Statistics::Add_Value(const CSG_String &Value, double Weight)
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 //---------------------------------------------------------
Get_Class_Index(const CSG_String & Value) const957 int CSG_Unique_String_Statistics::Get_Class_Index(const CSG_String &Value) const
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
971 ///////////////////////////////////////////////////////////
972 // //
973 // //
974 // //
975 ///////////////////////////////////////////////////////////
976
977 //---------------------------------------------------------
CSG_Category_Statistics(TSG_Data_Type Type)978 CSG_Category_Statistics::CSG_Category_Statistics(TSG_Data_Type Type)
979 {
980 m_pTable = new CSG_Table;
981
982 Create(Type);
983 }
984
985 //---------------------------------------------------------
~CSG_Category_Statistics(void)986 CSG_Category_Statistics::~CSG_Category_Statistics(void)
987 {
988 delete(m_pTable);
989 }
990
991 //---------------------------------------------------------
Create(TSG_Data_Type Type)992 void CSG_Category_Statistics::Create(TSG_Data_Type Type)
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 //---------------------------------------------------------
Destroy(void)1001 void CSG_Category_Statistics::Destroy(void)
1002 {
1003 m_pTable->Del_Records();
1004 }
1005
1006 //---------------------------------------------------------
Get_Category_Type(void) const1007 TSG_Data_Type CSG_Category_Statistics::Get_Category_Type(void) const
1008 {
1009 return( m_pTable->Get_Field_Type(0) );
1010 }
1011
1012
1013 ///////////////////////////////////////////////////////////
1014 // //
1015 ///////////////////////////////////////////////////////////
1016
1017 //---------------------------------------------------------
Get_Category(int Value) const1018 int CSG_Category_Statistics::Get_Category(int Value) const
1019 {
1020 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value, m_pTable->Get_Count() > 10);
1021
1022 return( pRecord ? pRecord->Get_Index() : -1);
1023 }
1024
1025 //---------------------------------------------------------
Get_Category(double Value) const1026 int CSG_Category_Statistics::Get_Category(double Value) const
1027 {
1028 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value, m_pTable->Get_Count() > 10);
1029
1030 return( pRecord ? pRecord->Get_Index() : -1);
1031 }
1032
1033 //---------------------------------------------------------
Get_Category(const CSG_String & Value) const1034 int CSG_Category_Statistics::Get_Category(const CSG_String &Value) const
1035 {
1036 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value, m_pTable->Get_Count() > 10);
1037
1038 return( pRecord ? pRecord->Get_Index() : -1);
1039 }
1040
1041
1042 ///////////////////////////////////////////////////////////
1043 // //
1044 ///////////////////////////////////////////////////////////
1045
1046 //---------------------------------------------------------
Add_Value(int Value)1047 int CSG_Category_Statistics::Add_Value(int Value)
1048 {
1049 int i = Get_Category(Value);
1050
1051 CSG_Table_Record *pRecord = m_pTable->Get_Record_byIndex(i);
1052
1053 if( !pRecord )
1054 {
1055 i = m_pTable->Get_Count();
1056
1057 (pRecord = m_pTable->Add_Record())->Set_Value(0, Value);
1058 }
1059
1060 pRecord->Add_Value(1, 1);
1061
1062 return( i );
1063 }
1064
1065 //---------------------------------------------------------
Add_Value(double Value)1066 int CSG_Category_Statistics::Add_Value(double Value)
1067 {
1068 int i = Get_Category(Value);
1069
1070 CSG_Table_Record *pRecord = m_pTable->Get_Record_byIndex(i);
1071
1072 if( !pRecord )
1073 {
1074 i = m_pTable->Get_Count();
1075
1076 (pRecord = m_pTable->Add_Record())->Set_Value(0, Value);
1077 }
1078
1079 pRecord->Add_Value(1, 1);
1080
1081 return( i );
1082 }
1083
1084 //---------------------------------------------------------
Add_Value(const CSG_String & Value)1085 int CSG_Category_Statistics::Add_Value(const CSG_String &Value)
1086 {
1087 int i = Get_Category(Value);
1088
1089 CSG_Table_Record *pRecord = m_pTable->Get_Record_byIndex(i);
1090
1091 if( !pRecord )
1092 {
1093 i = m_pTable->Get_Count();
1094
1095 (pRecord = m_pTable->Add_Record())->Set_Value(0, Value);
1096 }
1097
1098 pRecord->Add_Value(1, 1);
1099
1100 return( i );
1101 }
1102
1103 //---------------------------------------------------------
1104 // sort categories ascending
Sort(void)1105 bool CSG_Category_Statistics::Sort(void)
1106 {
1107 return( m_pTable->Set_Index(0, TABLE_INDEX_Ascending) );
1108 }
1109
1110
1111 ///////////////////////////////////////////////////////////
1112 // //
1113 ///////////////////////////////////////////////////////////
1114
1115 //---------------------------------------------------------
1116 // returns the number of categories.
Get_Count(void) const1117 int CSG_Category_Statistics::Get_Count(void) const
1118 {
1119 return( m_pTable->Get_Count() );
1120 }
1121
1122 //---------------------------------------------------------
1123 // returns the number of observations for the i'th category.
Get_Count(int i) const1124 int CSG_Category_Statistics::Get_Count(int i) const
1125 {
1126 CSG_Table_Record *pRecord = m_pTable->Get_Record_byIndex(i);
1127
1128 return( pRecord ? pRecord->asInt(1) : 0 );
1129 }
1130
1131 //---------------------------------------------------------
asInt(int i) const1132 int CSG_Category_Statistics::asInt(int i) const
1133 {
1134 CSG_Table_Record *pRecord = m_pTable->Get_Record_byIndex(i);
1135
1136 return( pRecord ? pRecord->asInt(0) : 0 );
1137 }
1138
1139 //---------------------------------------------------------
asDouble(int i) const1140 double CSG_Category_Statistics::asDouble(int i) const
1141 {
1142 CSG_Table_Record *pRecord = m_pTable->Get_Record_byIndex(i);
1143
1144 return( pRecord ? pRecord->asDouble(0) : 0 );
1145 }
1146
1147 //---------------------------------------------------------
asString(int i) const1148 CSG_String CSG_Category_Statistics::asString(int i) const
1149 {
1150 CSG_Table_Record *pRecord = m_pTable->Get_Record_byIndex(i);
1151
1152 return( pRecord ? pRecord->asString(0) : SG_T("") );
1153 }
1154
1155
1156 ///////////////////////////////////////////////////////////
1157 // //
1158 ///////////////////////////////////////////////////////////
1159
1160 //---------------------------------------------------------
Get_Majority(void)1161 int CSG_Category_Statistics::Get_Majority(void)
1162 {
1163 if( m_pTable->Get_Count() > 0 )
1164 {
1165 int Index = 0, Count = m_pTable->Get_Record_byIndex(0)->asInt(1);
1166
1167 for(int i=1; i<m_pTable->Get_Count(); i++)
1168 {
1169 if( Count < m_pTable->Get_Record_byIndex(i)->asInt(1) )
1170 {
1171 Index = i;
1172 Count = m_pTable->Get_Record_byIndex(i)->asInt(1);
1173 }
1174 }
1175
1176 return( Index );
1177 }
1178
1179 return( -1 );
1180 }
1181
1182 //---------------------------------------------------------
Get_Minority(void)1183 int CSG_Category_Statistics::Get_Minority(void)
1184 {
1185 if( m_pTable->Get_Count() > 0 )
1186 {
1187 int Index = 0, Count = m_pTable->Get_Record_byIndex(0)->asInt(1);
1188
1189 for(int i=1; i<m_pTable->Get_Count(); i++)
1190 {
1191 if( Count > m_pTable->Get_Record_byIndex(i)->asInt(1) )
1192 {
1193 Index = i;
1194 Count = m_pTable->Get_Record_byIndex(i)->asInt(1);
1195 }
1196 }
1197
1198 return( Index );
1199 }
1200
1201 return( -1 );
1202 }
1203
1204
1205 ///////////////////////////////////////////////////////////
1206 // //
1207 // //
1208 // //
1209 ///////////////////////////////////////////////////////////
1210
1211 //---------------------------------------------------------
CSG_Histogram(void)1212 CSG_Histogram::CSG_Histogram(void)
1213 {
1214 _On_Construction();
1215 }
1216
1217 //---------------------------------------------------------
CSG_Histogram(const CSG_Histogram & Histogram)1218 CSG_Histogram::CSG_Histogram(const CSG_Histogram &Histogram)
1219 {
1220 _On_Construction();
1221
1222 Create(Histogram);
1223 }
1224
1225 //---------------------------------------------------------
CSG_Histogram(size_t nClasses,double Minimum,double Maximum)1226 CSG_Histogram::CSG_Histogram(size_t nClasses, double Minimum, double Maximum)
1227 {
1228 _On_Construction();
1229
1230 Create(nClasses, Minimum, Maximum);
1231 }
1232
1233 //---------------------------------------------------------
CSG_Histogram(size_t nClasses,double Minimum,double Maximum,const CSG_Vector & Values,size_t maxSamples)1234 CSG_Histogram::CSG_Histogram(size_t nClasses, double Minimum, double Maximum, const CSG_Vector &Values, size_t maxSamples)
1235 {
1236 _On_Construction();
1237
1238 Create(nClasses, Minimum, Maximum, Values, maxSamples);
1239 }
1240
1241 //---------------------------------------------------------
CSG_Histogram(size_t nClasses,double Minimum,double Maximum,CSG_Table * pTable,int Field,size_t maxSamples)1242 CSG_Histogram::CSG_Histogram(size_t nClasses, double Minimum, double Maximum, CSG_Table *pTable, int Field, size_t maxSamples)
1243 {
1244 _On_Construction();
1245
1246 Create(nClasses, Minimum, Maximum, pTable, Field, maxSamples);
1247 }
1248
1249 //---------------------------------------------------------
CSG_Histogram(size_t nClasses,double Minimum,double Maximum,class CSG_Grid * pGrid,size_t maxSamples)1250 CSG_Histogram::CSG_Histogram(size_t nClasses, double Minimum, double Maximum, class CSG_Grid *pGrid, size_t maxSamples)
1251 {
1252 _On_Construction();
1253
1254 Create(nClasses, Minimum, Maximum, pGrid, maxSamples);
1255 }
1256
1257 //---------------------------------------------------------
CSG_Histogram(size_t nClasses,double Minimum,double Maximum,CSG_Grids * pGrids,size_t maxSamples)1258 CSG_Histogram::CSG_Histogram(size_t nClasses, double Minimum, double Maximum, CSG_Grids *pGrids, size_t maxSamples)
1259 {
1260 _On_Construction();
1261
1262 Create(nClasses, Minimum, Maximum, pGrids, maxSamples);
1263 }
1264
1265 //---------------------------------------------------------
~CSG_Histogram(void)1266 CSG_Histogram::~CSG_Histogram(void)
1267 {
1268 Destroy();
1269 }
1270
1271 //---------------------------------------------------------
Destroy(void)1272 bool CSG_Histogram::Destroy(void)
1273 {
1274 m_Statistics.Create();
1275
1276 SG_FREE_SAFE(m_Elements );
1277 SG_FREE_SAFE(m_Cumulative);
1278
1279 _On_Construction();
1280
1281 return( true );
1282 }
1283
1284
1285 ///////////////////////////////////////////////////////////
1286 // //
1287 ///////////////////////////////////////////////////////////
1288
1289 //---------------------------------------------------------
_On_Construction(void)1290 void CSG_Histogram::_On_Construction(void)
1291 {
1292 m_nClasses = 0;
1293 m_Elements = NULL;
1294 m_Cumulative = NULL;
1295 m_Minimum = 0.;
1296 m_Maximum = 0.;
1297 m_ClassWidth = 1.;
1298 }
1299
1300 //---------------------------------------------------------
_Create(size_t nClasses,double Minimum,double Maximum)1301 bool CSG_Histogram::_Create(size_t nClasses, double Minimum, double Maximum)
1302 {
1303 if( nClasses > 0 && Minimum < Maximum )
1304 {
1305 Destroy();
1306
1307 m_Elements = (size_t *)SG_Calloc(nClasses, sizeof(size_t));
1308 m_Cumulative = (size_t *)SG_Calloc(nClasses, sizeof(size_t));
1309
1310 if( m_Elements && m_Cumulative )
1311 {
1312 m_nClasses = nClasses;
1313 m_Minimum = Minimum;
1314 m_Maximum = Maximum;
1315 m_ClassWidth = (Maximum - Minimum) / (double)m_nClasses;
1316
1317 return( true );
1318 }
1319 }
1320
1321 Destroy();
1322
1323 return( false );
1324 }
1325
1326 //---------------------------------------------------------
Add_Value(double Value)1327 void CSG_Histogram::Add_Value(double Value)
1328 {
1329 m_Statistics += Value;
1330
1331 if( m_nClasses > 0 && m_Minimum <= Value && Value <= m_Maximum )
1332 {
1333 size_t Class = (size_t)((Value - m_Minimum) / m_ClassWidth);
1334
1335 if( Class >= m_nClasses )
1336 {
1337 Class = m_nClasses - 1;
1338 }
1339
1340 m_Elements[Class]++;
1341 }
1342 }
1343
1344 //---------------------------------------------------------
Scale_Element_Count(double Scale)1345 bool CSG_Histogram::Scale_Element_Count(double Scale)
1346 {
1347 if( m_nClasses > 0 && Scale > 0. )
1348 {
1349 m_Statistics.Set_Count((sLong)(Scale * Get_Element_Count()));
1350
1351 for(size_t i=0; i<m_nClasses; i++)
1352 {
1353 m_Elements[i] = (size_t)(Scale * m_Elements[i]);
1354 }
1355
1356 return( Update() );
1357 }
1358
1359 return( false );
1360 }
1361
1362 //---------------------------------------------------------
Update(void)1363 bool CSG_Histogram::Update(void)
1364 {
1365 if( m_nClasses > 0 )
1366 {
1367 m_Statistics.Get_Mean(); // _Evaluate()
1368
1369 m_nMaximum = m_Cumulative[0] = m_Elements[0];
1370
1371 for(size_t i=1; i<m_nClasses; i++)
1372 {
1373 m_Cumulative[i] = m_Cumulative[i - 1] + m_Elements[i];
1374
1375 if( m_nMaximum < m_Elements[i] )
1376 {
1377 m_nMaximum = m_Elements[i];
1378 }
1379 }
1380
1381 return( Get_Element_Count() > 0 );
1382 }
1383
1384 return( false );
1385 }
1386
1387 //---------------------------------------------------------
_Update(sLong nElements)1388 bool CSG_Histogram::_Update(sLong nElements)
1389 {
1390 if( nElements > 0 && m_Statistics.Get_Count() > 0 )
1391 {
1392 double Scale = (double)nElements / (double)m_Statistics.Get_Count();
1393
1394 m_Statistics.Create(m_Statistics.Get_Mean(), m_Statistics.Get_StdDev(), nElements);
1395
1396 for(size_t i=1; i<m_nClasses; i++)
1397 {
1398 m_Elements[i] = (size_t)(0.5 + Scale * m_Elements[i]);
1399 }
1400 }
1401
1402 return( Update() );
1403 }
1404
1405 //---------------------------------------------------------
1406 /**
1407 * Returns the correspondend value for the requested quantile.
1408 */
Get_Quantile(double Quantile) const1409 double CSG_Histogram::Get_Quantile(double Quantile) const
1410 {
1411 if( m_nClasses < 2 ) { return( 0. ); }
1412
1413 if( Quantile <= 0. ) { return( m_Minimum ); }
1414 if( Quantile >= 1. ) { return( m_Maximum ); }
1415
1416 size_t n = (size_t)(Quantile * Get_Element_Count()); // number of elements
1417
1418 for(size_t i=0, n0=0; i<m_nClasses; n0=m_Cumulative[i++])
1419 {
1420 if( n < m_Cumulative[i] )
1421 {
1422 if( m_Cumulative[i] >= n0 )
1423 {
1424 return( Get_Center(i) );
1425 }
1426
1427 double d = (n - n0) / (double)(m_Cumulative[i] - n0);
1428
1429 return( Get_Break(i) + d * m_ClassWidth );
1430 }
1431 else if( n == m_Cumulative[i] )
1432 {
1433 return( Get_Break(i + 1) );
1434 }
1435 }
1436
1437 return( m_Maximum );
1438 }
1439
1440 //---------------------------------------------------------
1441 /**
1442 * Returns the correspondend value for the requested percentile.
1443 */
Get_Percentile(double Percentile) const1444 double CSG_Histogram::Get_Percentile(double Percentile) const
1445 {
1446 return( Get_Quantile(Percentile / 100.) );
1447 }
1448
1449 //---------------------------------------------------------
1450 /**
1451 * Returns the correspondend quantile for the requested value.
1452 */
Get_Quantile_Value(double Value) const1453 double CSG_Histogram::Get_Quantile_Value(double Value) const
1454 {
1455 if( m_nClasses < 2 ) { return( 0. ); }
1456
1457 if( Value <= m_Minimum ) { return( 0. ); }
1458 if( Value >= m_Maximum ) { return( 1. ); }
1459
1460 size_t Class = (size_t)(m_nClasses * (Value - m_Minimum) / (m_Maximum - m_Minimum));
1461
1462 if( Class >= m_nClasses )
1463 {
1464 return( 1. );
1465 }
1466
1467 if( Class < 1 )
1468 {
1469 double dq = m_Cumulative[Class] / (double)Get_Element_Count();
1470
1471 return( dq * (Value - m_Minimum) / m_ClassWidth );
1472 }
1473
1474 double q0 = m_Cumulative[Class - 1] / (double)Get_Element_Count();
1475 double dq = (m_Cumulative[Class ] / (double)Get_Element_Count()) - q0;
1476
1477 return( q0 + dq * (Value - Get_Break(Class)) / m_ClassWidth );
1478 }
1479
1480 //---------------------------------------------------------
1481 /**
1482 * Returns the correspondend percentile for the requested value.
1483 */
Get_Percentile_Value(double Value) const1484 double CSG_Histogram::Get_Percentile_Value(double Value) const
1485 {
1486 return( Get_Quantile_Value(Value) * 100. );
1487 }
1488
1489
1490 ///////////////////////////////////////////////////////////
1491 // //
1492 ///////////////////////////////////////////////////////////
1493
1494 //---------------------------------------------------------
Create(const CSG_Histogram & Histogram)1495 bool CSG_Histogram::Create(const CSG_Histogram &Histogram)
1496 {
1497 if( !_Create(Histogram.m_nClasses, Histogram.m_Minimum, Histogram.m_Maximum) )
1498 {
1499 return( false );
1500 }
1501
1502 m_Statistics = Histogram.m_Statistics;
1503 m_ClassWidth = Histogram.m_ClassWidth;
1504 m_nMaximum = Histogram.m_nMaximum ;
1505
1506 for(size_t i=0; i<m_nClasses; i++)
1507 {
1508 m_Cumulative[i] = Histogram.m_Cumulative[i];
1509 m_Elements [i] = Histogram.m_Elements [i];
1510 }
1511
1512 return( true );
1513 }
1514
1515 //---------------------------------------------------------
Create(size_t nClasses,double Minimum,double Maximum)1516 bool CSG_Histogram::Create(size_t nClasses, double Minimum, double Maximum)
1517 {
1518 return( _Create(nClasses, Minimum, Maximum) );
1519 }
1520
1521 //---------------------------------------------------------
Create(size_t nClasses,double Minimum,double Maximum,const CSG_Vector & Values,size_t maxSamples)1522 bool CSG_Histogram::Create(size_t nClasses, double Minimum, double Maximum, const CSG_Vector &Values, size_t maxSamples)
1523 {
1524 if( Minimum >= Maximum )
1525 {
1526 CSG_Simple_Statistics s(Values);
1527
1528 Minimum = s.Get_Minimum();
1529 Maximum = s.Get_Maximum();
1530 }
1531
1532 if( !_Create(nClasses, Minimum, Maximum) )
1533 {
1534 return( false );
1535 }
1536
1537 //-----------------------------------------------------
1538 if( maxSamples > 0 && maxSamples < (size_t)Values.Get_N() )
1539 {
1540 double d = (double)Values.Get_N() / (double)maxSamples;
1541
1542 for(double i=0; i<(double)Values.Get_N(); i+=d)
1543 {
1544 Add_Value(Values[(sLong)i]);
1545 }
1546
1547 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1548
1549 return( _Update(d < 1. ? (int)(d * (double)Values.Get_N()) : Values.Get_N()) );
1550 }
1551
1552 //-----------------------------------------------------
1553 for(int i=0; i<Values.Get_N(); i++)
1554 {
1555 Add_Value(Values[i]);
1556 }
1557
1558 return( Update() );
1559 }
1560
1561 //---------------------------------------------------------
Create(size_t nClasses,double Minimum,double Maximum,CSG_Table * pTable,int Field,size_t maxSamples)1562 bool CSG_Histogram::Create(size_t nClasses, double Minimum, double Maximum, CSG_Table *pTable, int Field, size_t maxSamples)
1563 {
1564 if( !pTable || Field < 0 || Field >= pTable->Get_Field_Count() || !_Create(nClasses,
1565 Minimum < Maximum ? Minimum : pTable->Get_Minimum(Field),
1566 Minimum < Maximum ? Maximum : pTable->Get_Maximum(Field)) )
1567 {
1568 return( false );
1569 }
1570
1571 //-----------------------------------------------------
1572 if( maxSamples > 0 && maxSamples < (size_t)pTable->Get_Count() )
1573 {
1574 double d = (double)pTable->Get_Count() / (double)maxSamples;
1575
1576 for(double i=0; i<(double)pTable->Get_Count(); i+=d)
1577 {
1578 double Value = pTable->Get_Record((int)i)->asDouble(Field);
1579
1580 if( !pTable->is_NoData_Value(Value) )
1581 {
1582 Add_Value(Value);
1583 }
1584 }
1585
1586 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1587
1588 return( _Update(d < 1. ? (int)(d * pTable->Get_Count()) : pTable->Get_Count()) );
1589 }
1590
1591 //-----------------------------------------------------
1592 for(int i=0; i<pTable->Get_Count(); i++)
1593 {
1594 double Value = pTable->Get_Record(i)->asDouble(Field);
1595
1596 if( !pTable->is_NoData_Value(Value) )
1597 {
1598 Add_Value(Value);
1599 }
1600 }
1601
1602 return( Update() );
1603 }
1604
1605 //---------------------------------------------------------
Create(size_t nClasses,double Minimum,double Maximum,CSG_Grid * pGrid,size_t maxSamples)1606 bool CSG_Histogram::Create(size_t nClasses, double Minimum, double Maximum, CSG_Grid *pGrid, size_t maxSamples)
1607 {
1608 if( !pGrid || !_Create(nClasses,
1609 Minimum < Maximum ? Minimum : pGrid->Get_Min(),
1610 Minimum < Maximum ? Maximum : pGrid->Get_Max()) )
1611 {
1612 return( false );
1613 }
1614
1615 //-----------------------------------------------------
1616 if( maxSamples > 0 && (sLong)maxSamples < pGrid->Get_NCells() )
1617 {
1618 double d = (double)pGrid->Get_NCells() / (double)maxSamples;
1619
1620 for(double i=0; i<(double)pGrid->Get_NCells(); i+=d)
1621 {
1622 double Value = pGrid->asDouble((sLong)i);
1623
1624 if( !pGrid->is_NoData_Value(Value) )
1625 {
1626 Add_Value(Value);
1627 }
1628 }
1629
1630 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1631
1632 return( _Update(d < 1. ? (sLong)(d * (double)pGrid->Get_NCells()) : pGrid->Get_NCells()) );
1633 }
1634
1635 //-----------------------------------------------------
1636 for(sLong i=0; i<pGrid->Get_NCells(); i++)
1637 {
1638 if( !pGrid->is_NoData(i) )
1639 {
1640 Add_Value(pGrid->asDouble(i));
1641 }
1642 }
1643
1644 return( Update() );
1645 }
1646
1647 //---------------------------------------------------------
Create(size_t nClasses,double Minimum,double Maximum,CSG_Grids * pGrids,size_t maxSamples)1648 bool CSG_Histogram::Create(size_t nClasses, double Minimum, double Maximum, CSG_Grids *pGrids, size_t maxSamples)
1649 {
1650 if( !pGrids || !_Create(nClasses,
1651 Minimum < Maximum ? Minimum : pGrids->Get_Min(),
1652 Minimum < Maximum ? Maximum : pGrids->Get_Max()) )
1653 {
1654 return( false );
1655 }
1656
1657 //-----------------------------------------------------
1658 if( maxSamples > 0 && (sLong)maxSamples < pGrids->Get_NCells() )
1659 {
1660 double d = (double)pGrids->Get_NCells() / (double)maxSamples;
1661
1662 for(double i=0; i<(double)pGrids->Get_NCells(); i+=d)
1663 {
1664 double Value = pGrids->asDouble((sLong)i);
1665
1666 if( !pGrids->is_NoData_Value(Value) )
1667 {
1668 Add_Value(Value);
1669 }
1670 }
1671
1672 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1673
1674 return( _Update(d < 1. ? (sLong)(d * (double)pGrids->Get_NCells()) : pGrids->Get_NCells()) );
1675 }
1676
1677 //-----------------------------------------------------
1678 for(sLong i=0; i<pGrids->Get_NCells(); i++)
1679 {
1680 if( !pGrids->is_NoData(i) )
1681 {
1682 Add_Value(pGrids->asDouble(i));
1683 }
1684 }
1685
1686 return( Update() );
1687 }
1688
1689
1690 ///////////////////////////////////////////////////////////
1691 // //
1692 ///////////////////////////////////////////////////////////
1693
1694 //---------------------------------------------------------
operator =(const CSG_Histogram & Histogram)1695 CSG_Histogram & CSG_Histogram::operator = (const CSG_Histogram &Histogram)
1696 {
1697 Create(Histogram);
1698
1699 return( *this );
1700 }
1701
1702
1703 ///////////////////////////////////////////////////////////
1704 // //
1705 // //
1706 // //
1707 ///////////////////////////////////////////////////////////
1708
1709 //---------------------------------------------------------
CSG_Natural_Breaks(void)1710 CSG_Natural_Breaks::CSG_Natural_Breaks(void)
1711 {}
1712
1713 //---------------------------------------------------------
~CSG_Natural_Breaks(void)1714 CSG_Natural_Breaks::~CSG_Natural_Breaks(void)
1715 {}
1716
1717 //---------------------------------------------------------
CSG_Natural_Breaks(CSG_Table * pTable,int Field,int nClasses,int Histogram)1718 CSG_Natural_Breaks::CSG_Natural_Breaks(CSG_Table *pTable, int Field, int nClasses, int Histogram)
1719 {
1720 Create(pTable, Field, nClasses, Histogram);
1721 }
1722
1723 //---------------------------------------------------------
CSG_Natural_Breaks(CSG_Grid * pGrid,int nClasses,int Histogram)1724 CSG_Natural_Breaks::CSG_Natural_Breaks(CSG_Grid *pGrid, int nClasses, int Histogram)
1725 {
1726 Create(pGrid, nClasses, Histogram);
1727 }
1728
1729 //---------------------------------------------------------
CSG_Natural_Breaks(CSG_Grids * pGrids,int nClasses,int Histogram)1730 CSG_Natural_Breaks::CSG_Natural_Breaks(CSG_Grids *pGrids, int nClasses, int Histogram)
1731 {
1732 Create(pGrids, nClasses, Histogram);
1733 }
1734
1735 //---------------------------------------------------------
CSG_Natural_Breaks(const CSG_Vector & Values,int nClasses,int Histogram)1736 CSG_Natural_Breaks::CSG_Natural_Breaks(const CSG_Vector &Values, int nClasses, int Histogram)
1737 {
1738 Create(Values, nClasses, Histogram);
1739 }
1740
1741
1742 ///////////////////////////////////////////////////////////
1743 // //
1744 ///////////////////////////////////////////////////////////
1745
1746 //---------------------------------------------------------
Create(CSG_Table * pTable,int Field,int nClasses,int Histogram)1747 bool CSG_Natural_Breaks::Create(CSG_Table *pTable, int Field, int nClasses, int Histogram)
1748 {
1749 bool bResult = false;
1750
1751 if( Histogram > 0 )
1752 {
1753 bResult = m_Histogram.Create(Histogram, 0, 0, pTable, Field) && _Histogram(nClasses);
1754 }
1755 else if( Field >= 0 && Field < pTable->Get_Field_Count() )
1756 {
1757 for(int i=0; i<pTable->Get_Count(); i++)
1758 {
1759 CSG_Table_Record *pRecord = pTable->Get_Record(i);
1760
1761 if( !pRecord->is_NoData(Field) )
1762 {
1763 m_Values.Add_Row(pRecord->asDouble(Field));
1764 }
1765 }
1766
1767 bResult = m_Values.Sort() && _Calculate(nClasses);
1768
1769 m_Values.Destroy();
1770 }
1771
1772 return( bResult );
1773 }
1774
1775 //---------------------------------------------------------
Create(CSG_Grid * pGrid,int nClasses,int Histogram)1776 bool CSG_Natural_Breaks::Create(CSG_Grid *pGrid, int nClasses, int Histogram)
1777 {
1778 bool bResult = false;
1779
1780 if( Histogram > 0 )
1781 {
1782 bResult = m_Histogram.Create(Histogram, 0, 0, pGrid) && _Histogram(nClasses);
1783 }
1784 else
1785 {
1786 for(sLong i=0; i<pGrid->Get_NCells(); i++)
1787 {
1788 if( !pGrid->is_NoData(i) )
1789 {
1790 m_Values.Add_Row(pGrid->asDouble(i));
1791 }
1792 }
1793
1794 bResult = m_Values.Sort() && _Calculate(nClasses);
1795
1796 m_Values.Destroy();
1797 }
1798
1799 return( bResult );
1800 }
1801
1802 //---------------------------------------------------------
Create(CSG_Grids * pGrids,int nClasses,int Histogram)1803 bool CSG_Natural_Breaks::Create(CSG_Grids *pGrids, int nClasses, int Histogram)
1804 {
1805 bool bResult = false;
1806
1807 if( Histogram > 0 )
1808 {
1809 bResult = m_Histogram.Create(Histogram, 0, 0, pGrids) && _Histogram(nClasses);
1810 }
1811 else
1812 {
1813 for(sLong i=0; i<pGrids->Get_NCells(); i++)
1814 {
1815 if( !pGrids->is_NoData(i) )
1816 {
1817 m_Values.Add_Row(pGrids->asDouble(i));
1818 }
1819 }
1820
1821 bResult = m_Values.Sort() && _Calculate(nClasses);
1822
1823 m_Values.Destroy();
1824 }
1825
1826 return( bResult );
1827 }
1828
1829 //---------------------------------------------------------
Create(const CSG_Vector & Values,int nClasses,int Histogram)1830 bool CSG_Natural_Breaks::Create(const CSG_Vector &Values, int nClasses, int Histogram)
1831 {
1832 bool bResult = false;
1833
1834 if( Histogram > 0 )
1835 {
1836 bResult = m_Histogram.Create(Histogram, 0, 0, Values) && _Histogram(nClasses);
1837 }
1838 else
1839 {
1840 bResult = m_Values.Create(Values) && m_Values.Sort() && _Calculate(nClasses);
1841
1842 m_Values.Destroy();
1843 }
1844
1845 return( bResult );
1846 }
1847
1848
1849 ///////////////////////////////////////////////////////////
1850 // //
1851 ///////////////////////////////////////////////////////////
1852
1853 //---------------------------------------------------------
_Histogram(int nClasses)1854 bool CSG_Natural_Breaks::_Histogram(int nClasses)
1855 {
1856 if( _Calculate(nClasses) )
1857 {
1858 double d = (double)m_Histogram.Get_Class_Count() / m_Histogram.Get_Cumulative((int)(m_Histogram.Get_Class_Count() - 1));
1859
1860 m_Breaks[0] = m_Histogram.Get_Break(0);
1861
1862 for(int i=1; i<Get_Count(); i++)
1863 {
1864 m_Breaks[i] = m_Histogram.Get_Value(m_Breaks[i] * d);
1865 }
1866
1867 m_Breaks[nClasses] = m_Histogram.Get_Break((int)m_Histogram.Get_Class_Count());
1868
1869 m_Histogram.Destroy();
1870
1871 return( true );
1872 }
1873
1874 m_Histogram.Destroy();
1875
1876 return( false );
1877 }
1878
1879 //---------------------------------------------------------
_Get_Value(int i)1880 inline double CSG_Natural_Breaks::_Get_Value(int i)
1881 {
1882 if( m_Histogram.Get_Class_Count() > 0 )
1883 {
1884 return( (double)m_Histogram.Get_Cumulative(i) );
1885 }
1886
1887 return( m_Values[i] );
1888 }
1889
1890 //---------------------------------------------------------
_Calculate(int nClasses)1891 bool CSG_Natural_Breaks::_Calculate(int nClasses)
1892 {
1893 if( m_Histogram.Get_Class_Count() == 0 && m_Values.Get_Size() == 0 )
1894 {
1895 return( false );
1896 }
1897
1898 int nValues = m_Histogram.Get_Class_Count() > 0 ? (int)m_Histogram.Get_Class_Count() : m_Values.Get_N();
1899
1900 CSG_Matrix mv(nClasses, nValues); mv.Assign(FLT_MAX);
1901
1902 int i, **mc = (int **)SG_Malloc(nValues * sizeof(int *));
1903
1904 mc[0] = (int *)SG_Calloc(nClasses * nValues, sizeof(int));
1905
1906 for(i=0; i<nValues; i++)
1907 {
1908 mc[i] = mc[0] + i * nClasses;
1909 }
1910
1911 //-----------------------------------------------------
1912 for(i=1; i<nValues; i++)
1913 {
1914 double v = 0., s1 = 0., s2 = 0., w = 0.;
1915
1916 for(int m=0, n=i+1; m<=i; m++, n--)
1917 {
1918 v = _Get_Value(n);
1919 s2 += v*v;
1920 s1 += v;
1921 w ++;
1922 v = s2 - (s1 * s1) / w;
1923
1924 if( n > 0 )
1925 {
1926 for(int j=1; j<nClasses; j++)
1927 {
1928 if( mv[i][j] >= (v + mv[n - 1][j - 1]) )
1929 {
1930 mc[i][j] = n;
1931 mv[i][j] = v + mv[n - 1][j - 1];
1932 }
1933 }
1934 }
1935 }
1936
1937 mc[i][0] = 0;
1938 mv[i][0] = v;
1939 }
1940
1941 //-----------------------------------------------------
1942 CSG_Array_Int Class(nClasses);
1943
1944 for(i=0; i<nClasses; i++)
1945 {
1946 Class[i] = i;
1947 }
1948
1949 int j = Class[nClasses - 1] = nValues - 1;
1950
1951 for(i=nClasses-1; i>0; i--)
1952 {
1953 Class[i - 1] = j = mc[j - 1][i];
1954 }
1955
1956 //-----------------------------------------------------
1957 m_Breaks.Create(nClasses + 1);
1958
1959 m_Breaks[0] = _Get_Value(0);
1960
1961 for(i=1; i<nClasses; i++)
1962 {
1963 m_Breaks[i] = _Get_Value(Class[i - 1]);
1964 }
1965
1966 m_Breaks[nClasses] = _Get_Value(nValues - 1);
1967
1968 SG_Free(mc[0]); SG_Free(mc);
1969
1970 return( true );
1971 }
1972
1973
1974 ///////////////////////////////////////////////////////////
1975 // //
1976 // //
1977 // //
1978 ///////////////////////////////////////////////////////////
1979
1980 //---------------------------------------------------------
CSG_Cluster_Analysis(void)1981 CSG_Cluster_Analysis::CSG_Cluster_Analysis(void)
1982 {
1983 m_nFeatures = 0;
1984 m_Iteration = 0;
1985 }
1986
1987 //---------------------------------------------------------
~CSG_Cluster_Analysis(void)1988 CSG_Cluster_Analysis::~CSG_Cluster_Analysis(void)
1989 {
1990 Destroy();
1991 }
1992
1993 //---------------------------------------------------------
Destroy(void)1994 bool CSG_Cluster_Analysis::Destroy(void)
1995 {
1996 m_Centroid.Destroy();
1997 m_Variance.Destroy();
1998 m_nMembers.Destroy();
1999 m_Clusters.Destroy();
2000 m_Features.Destroy();
2001 m_nFeatures = 0;
2002 m_Iteration = 0;
2003
2004 return( true );
2005 }
2006
2007 //---------------------------------------------------------
Create(int nFeatures)2008 bool CSG_Cluster_Analysis::Create(int nFeatures)
2009 {
2010 Destroy();
2011
2012 if( nFeatures > 0 )
2013 {
2014 m_nFeatures = nFeatures;
2015
2016 m_Features.Create(m_nFeatures * sizeof(double), 0, SG_ARRAY_GROWTH_3);
2017
2018 return( true );
2019 }
2020
2021 return( false );
2022 }
2023
2024 //---------------------------------------------------------
Add_Element(void)2025 bool CSG_Cluster_Analysis::Add_Element(void)
2026 {
2027 return( m_nFeatures > 0 && m_Features.Inc_Array() );
2028 }
2029
2030 //---------------------------------------------------------
Set_Feature(int iElement,int iFeature,double Value)2031 bool CSG_Cluster_Analysis::Set_Feature(int iElement, int iFeature, double Value)
2032 {
2033 if( iElement >= 0 && iElement < Get_nElements() && iFeature >= 0 && iFeature < m_nFeatures )
2034 {
2035 ((double *)m_Features.Get_Entry(iElement))[iFeature] = Value;
2036
2037 return( true );
2038 }
2039
2040 return( false );
2041 }
2042
2043 //---------------------------------------------------------
2044 /**
2045 * Performs the cluster analysis using the features added prior
2046 * to this step. Method is minimum distance (= default), hill
2047 * climbing (= 1), or both methods in combination. If nMaxIterations
2048 * is set to zero, the analysis is iterated until it converges.
2049 * Initilization is done randomely (= default), periodically (= 1),
2050 * or skipped (= 2). The latter case allows to start the clustering
2051 * with user supplied start partitions.
2052 */
2053 //---------------------------------------------------------
Execute(int Method,int nClusters,int nMaxIterations,int Initialization)2054 bool CSG_Cluster_Analysis::Execute(int Method, int nClusters, int nMaxIterations, int Initialization)
2055 {
2056 if( Get_nElements() < 2 || nClusters < 2 )
2057 {
2058 return( false );
2059 }
2060
2061 //-----------------------------------------------------
2062 m_nMembers.Create(nClusters);
2063 m_Variance.Create(nClusters);
2064 m_Centroid.Create(m_nFeatures, nClusters);
2065
2066 //-----------------------------------------------------
2067 m_Clusters.Create(Get_nElements());
2068
2069 for(int iElement=0; iElement<Get_nElements(); iElement++)
2070 {
2071 switch( Initialization )
2072 {
2073 default: // random
2074 if( (m_Clusters[iElement] = (int)CSG_Random::Get_Uniform(0, nClusters)) >= nClusters )
2075 {
2076 m_Clusters[iElement] = nClusters - 1;
2077 }
2078 break;
2079
2080 case 1: // periodic
2081 {
2082 m_Clusters[iElement] = iElement % nClusters;
2083 }
2084 break;
2085
2086 case 2: // keep as is, but check for valid cluster ids
2087 if( 0 > m_Clusters[iElement] || m_Clusters[iElement] >= nClusters )
2088 {
2089 m_Clusters[iElement] = iElement % nClusters;
2090 }
2091 break;
2092 }
2093 }
2094
2095 //-----------------------------------------------------
2096 bool bResult;
2097
2098 m_Iteration = 0;
2099
2100 switch( Method )
2101 {
2102 default: bResult = _Minimum_Distance(true , nMaxIterations); break;
2103 case 1: bResult = _Hill_Climbing (true , nMaxIterations); break;
2104 case 2: bResult = _Minimum_Distance(true , nMaxIterations)
2105 && _Hill_Climbing (false, nMaxIterations); break;
2106 }
2107
2108 //-----------------------------------------------------
2109 if( bResult )
2110 {
2111 for(int iCluster=0; iCluster<nClusters; iCluster++)
2112 {
2113 m_Variance[iCluster] = m_nMembers[iCluster] <= 0 ? 0. : m_Variance[iCluster] / m_nMembers[iCluster];
2114 }
2115 }
2116
2117 return( bResult );
2118 }
2119
2120 //---------------------------------------------------------
_Minimum_Distance(bool bInitialize,int nMaxIterations)2121 bool CSG_Cluster_Analysis::_Minimum_Distance(bool bInitialize, int nMaxIterations)
2122 {
2123 int iElement, iCluster, nClusters = m_Variance.Get_N();
2124
2125 double SP_Last = -1.;
2126
2127 //-----------------------------------------------------
2128 for(m_Iteration=1; SG_UI_Process_Get_Okay(); m_Iteration++)
2129 {
2130 m_Variance = 0.;
2131 m_Centroid = 0.;
2132 m_nMembers = 0;
2133
2134 //-------------------------------------------------
2135 for(iElement=0; iElement<Get_nElements(); iElement++)
2136 {
2137 m_nMembers[iCluster = m_Clusters[iElement]]++;
2138
2139 double *Feature = (double *)m_Features.Get_Entry(iElement);
2140
2141 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2142 {
2143 m_Centroid[iCluster][iFeature] += Feature[iFeature];
2144 }
2145 }
2146
2147 //-------------------------------------------------
2148 for(iCluster=0; iCluster<nClusters; iCluster++)
2149 {
2150 double d = m_nMembers[iCluster] > 0 ? 1. / m_nMembers[iCluster] : 0.;
2151
2152 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2153 {
2154 m_Centroid[iCluster][iFeature] *= d;
2155 }
2156 }
2157
2158 //-------------------------------------------------
2159 int nShifts = 0;
2160
2161 m_SP = 0.;
2162
2163 for(iElement=0; iElement<Get_nElements(); iElement++)
2164 {
2165 double *Feature = (double *)m_Features.Get_Entry(iElement);
2166
2167 double minVariance = -1.;
2168 int minCluster = -1;
2169
2170 for(iCluster=0; iCluster<nClusters; iCluster++)
2171 {
2172 double Variance = 0.;
2173
2174 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2175 {
2176 Variance += SG_Get_Square(m_Centroid[iCluster][iFeature] - Feature[iFeature]);
2177 }
2178
2179 if( minVariance < 0. || Variance < minVariance )
2180 {
2181 minVariance = Variance;
2182 minCluster = iCluster;
2183 }
2184 }
2185
2186 if( m_Clusters[iElement] != minCluster )
2187 {
2188 m_Clusters[iElement] = minCluster;
2189
2190 nShifts++;
2191 }
2192
2193 m_SP += minVariance;
2194 m_Variance[minCluster] += minVariance;
2195 }
2196
2197 //-------------------------------------------------
2198 m_SP /= Get_nElements();
2199
2200 SG_UI_Process_Set_Text(CSG_String::Format("%s: %d >> %s %f",
2201 _TL("pass" ), m_Iteration,
2202 _TL("change"), m_Iteration < 2 ? m_SP : SP_Last - m_SP
2203 ));
2204
2205 SP_Last = m_SP;
2206
2207 if( nShifts == 0 || (nMaxIterations > 0 && nMaxIterations <= m_Iteration) )
2208 {
2209 return( true );
2210 }
2211 }
2212
2213 return( true );
2214 }
2215
2216 //---------------------------------------------------------
_Hill_Climbing(bool bInitialize,int nMaxIterations)2217 bool CSG_Cluster_Analysis::_Hill_Climbing(bool bInitialize, int nMaxIterations)
2218 {
2219 int iElement, iCluster, nClusters = m_Variance.Get_N();
2220
2221 //-----------------------------------------------------
2222 m_Variance = 0.;
2223 m_Centroid = 0.;
2224 m_nMembers = 0;
2225
2226 for(iElement=0; iElement<Get_nElements(); iElement++)
2227 {
2228 m_nMembers[iCluster = m_Clusters[iElement]]++;
2229
2230 double *Feature = (double *)m_Features.Get_Entry(iElement);
2231
2232 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2233 {
2234 double d = Feature[iFeature];
2235
2236 m_Centroid[iCluster][iFeature] += d;
2237 m_Variance[iCluster] += d*d;
2238 }
2239 }
2240
2241 //-----------------------------------------------------
2242 for(iCluster=0; iCluster<nClusters; iCluster++)
2243 {
2244 double v = 0., d = m_nMembers[iCluster] <= 0 ? 0. : 1. / (double)m_nMembers[iCluster];
2245
2246 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2247 {
2248 m_Centroid[iCluster][iFeature] *= d;
2249 v += SG_Get_Square(m_Centroid[iCluster][iFeature]);
2250 }
2251
2252 m_Variance[iCluster] -= v * m_nMembers[iCluster];
2253 }
2254
2255 //-----------------------------------------------------
2256 double SP_Last = -1.; int noShift = 0;
2257
2258 for(m_Iteration=1; SG_UI_Process_Get_Okay(false); m_Iteration++)
2259 {
2260 for(iElement=0; iElement<Get_nElements(); iElement++)
2261 {
2262 iCluster = m_Clusters[iElement];
2263
2264 if( noShift++ < Get_nElements() && m_nMembers[iCluster] > 1 )
2265 {
2266 int iFeature; double *Feature = (double *)m_Features.Get_Entry(iElement);
2267
2268 double V1, V2, Variance = 0.;
2269
2270 for(iFeature=0; iFeature<m_nFeatures; iFeature++)
2271 {
2272 Variance += SG_Get_Square(m_Centroid[iCluster][iFeature] - Feature[iFeature]);
2273 }
2274
2275 V1 = Variance * m_nMembers[iCluster] / (m_nMembers[iCluster] - 1.);
2276
2277 //-----------------------------------------
2278 int kCluster = 0;
2279 double VMin = -1.;
2280
2281 for(int jCluster=0; jCluster<nClusters; jCluster++)
2282 {
2283 if( jCluster != iCluster )
2284 {
2285 Variance = 0.;
2286
2287 for(iFeature=0; iFeature<m_nFeatures; iFeature++)
2288 {
2289 Variance += SG_Get_Square(m_Centroid[jCluster][iFeature] - Feature[iFeature]);
2290 }
2291
2292 V2 = Variance * m_nMembers[jCluster] / (m_nMembers[jCluster] + 1.);
2293
2294 if( VMin < 0. || V2 < VMin )
2295 {
2296 VMin = V2;
2297 kCluster = jCluster;
2298 }
2299 }
2300 }
2301
2302 //-----------------------------------------
2303 if( VMin >= 0 && VMin < V1 )
2304 {
2305 noShift = 0;
2306 m_Variance[iCluster] -= V1;
2307 m_Variance[kCluster] += VMin;
2308 V1 = 1. / (m_nMembers[iCluster] - 1.);
2309 V2 = 1. / (m_nMembers[kCluster] + 1.);
2310
2311 for(iFeature=0; iFeature<m_nFeatures; iFeature++)
2312 {
2313 double d = Feature[iFeature];
2314
2315 m_Centroid[iCluster][iFeature] = (m_nMembers[iCluster] * m_Centroid[iCluster][iFeature] - d) * V1;
2316 m_Centroid[kCluster][iFeature] = (m_nMembers[kCluster] * m_Centroid[kCluster][iFeature] + d) * V2;
2317 }
2318
2319 m_Clusters[iElement] = kCluster;
2320
2321 m_nMembers[iCluster]--;
2322 m_nMembers[kCluster]++;
2323 }
2324 }
2325 }
2326
2327 //-------------------------------------------------
2328 for(iCluster=0, m_SP=0.; iCluster<nClusters; iCluster++)
2329 {
2330 m_SP += m_Variance[iCluster];
2331 }
2332
2333 m_SP /= Get_nElements();
2334
2335 SG_UI_Process_Set_Text(CSG_String::Format("%s: %d >> %s %f",
2336 _TL("pass" ), m_Iteration,
2337 _TL("change"), m_Iteration <= 1 ? m_SP : SP_Last - m_SP
2338 ));
2339
2340 SP_Last = m_SP;
2341
2342 if( noShift >= Get_nElements() || (nMaxIterations > 0 && nMaxIterations <= m_Iteration) )
2343 {
2344 return( true );
2345 }
2346 }
2347
2348 return( true );
2349 }
2350
2351
2352 ///////////////////////////////////////////////////////////
2353 // //
2354 // //
2355 // //
2356 ///////////////////////////////////////////////////////////
2357
2358 //---------------------------------------------------------
CSG_Classifier_Supervised(void)2359 CSG_Classifier_Supervised::CSG_Classifier_Supervised(void)
2360 {
2361 m_nFeatures = 0;
2362
2363 m_nClasses = 0;
2364 m_pClasses = NULL;
2365
2366 m_Threshold_Distance = 0.;
2367 m_Threshold_Angle = 0.;
2368 m_Threshold_Probability = 0.;
2369 m_Probability_Relative = false;
2370
2371 for(int i=0; i<SG_CLASSIFY_SUPERVISED_WTA; i++)
2372 {
2373 m_bWTA[i] = i == SG_CLASSIFY_SUPERVISED_MinimumDistance
2374 // || i == SG_CLASSIFY_SUPERVISED_Mahalonobis
2375 || i == SG_CLASSIFY_SUPERVISED_MaximumLikelihood
2376 || i == SG_CLASSIFY_SUPERVISED_SAM;
2377 }
2378 }
2379
2380 //---------------------------------------------------------
~CSG_Classifier_Supervised(void)2381 CSG_Classifier_Supervised::~CSG_Classifier_Supervised(void)
2382 {
2383 Destroy();
2384 }
2385
2386 //---------------------------------------------------------
Create(int nFeatures)2387 void CSG_Classifier_Supervised::Create(int nFeatures)
2388 {
2389 Destroy();
2390
2391 if( nFeatures > 0 )
2392 {
2393 m_nFeatures = nFeatures;
2394 }
2395 }
2396
2397 //---------------------------------------------------------
Destroy(void)2398 void CSG_Classifier_Supervised::Destroy(void)
2399 {
2400 if( m_nClasses > 0 )
2401 {
2402 for(int i=0; i<m_nClasses; i++)
2403 {
2404 delete(m_pClasses[i]);
2405 }
2406
2407 SG_FREE_SAFE(m_pClasses);
2408 }
2409
2410 m_nFeatures = 0;
2411
2412 m_Info.Clear();
2413 }
2414
2415
2416 ///////////////////////////////////////////////////////////
2417 // //
2418 ///////////////////////////////////////////////////////////
2419
2420 //---------------------------------------------------------
Set_Threshold_Distance(double Value)2421 void CSG_Classifier_Supervised::Set_Threshold_Distance (double Value) { m_Threshold_Distance = Value; }
Get_Threshold_Distance(void)2422 double CSG_Classifier_Supervised::Get_Threshold_Distance (void) { return( m_Threshold_Distance ); }
2423
2424 //---------------------------------------------------------
Set_Threshold_Angle(double Value)2425 void CSG_Classifier_Supervised::Set_Threshold_Angle (double Value) { m_Threshold_Angle = Value; }
Get_Threshold_Angle(void)2426 double CSG_Classifier_Supervised::Get_Threshold_Angle (void) { return( m_Threshold_Angle ); }
2427
2428 //---------------------------------------------------------
Set_Threshold_Probability(double Value)2429 void CSG_Classifier_Supervised::Set_Threshold_Probability(double Value) { m_Threshold_Probability = Value; }
Get_Threshold_Probability(void)2430 double CSG_Classifier_Supervised::Get_Threshold_Probability(void) { return( m_Threshold_Probability ); }
2431
2432 //---------------------------------------------------------
Set_Probability_Relative(bool Value)2433 void CSG_Classifier_Supervised::Set_Probability_Relative (bool Value) { m_Probability_Relative = Value; }
Get_Probability_Relative(void)2434 bool CSG_Classifier_Supervised::Get_Probability_Relative (void) { return( m_Probability_Relative ); }
2435
2436 //---------------------------------------------------------
Set_WTA(int Method,bool bOn)2437 void CSG_Classifier_Supervised::Set_WTA(int Method, bool bOn)
2438 {
2439 if( Method >= 0 && Method < SG_CLASSIFY_SUPERVISED_WTA )
2440 {
2441 m_bWTA[Method] = bOn;
2442 }
2443 }
2444
Get_WTA(int Method)2445 bool CSG_Classifier_Supervised::Get_WTA(int Method)
2446 {
2447 return( Method >= 0 && Method < SG_CLASSIFY_SUPERVISED_WTA ? m_bWTA[Method] : false );
2448 }
2449
2450
2451 ///////////////////////////////////////////////////////////
2452 // //
2453 ///////////////////////////////////////////////////////////
2454
2455 //---------------------------------------------------------
2456 #include "saga_api.h"
2457
2458 //---------------------------------------------------------
Load(const CSG_String & File)2459 bool CSG_Classifier_Supervised::Load(const CSG_String &File)
2460 {
2461 int nFeatures = m_nFeatures; Destroy(); m_nFeatures = nFeatures;
2462
2463 //-----------------------------------------------------
2464 CSG_MetaData Data;
2465
2466 if( !Data.Load(File) || !Data.Cmp_Name("supervised_classifier") || SG_Compare_Version(Data.Get_Property("saga-version"), "2.1.4") < 0 )
2467 {
2468 return( false );
2469 }
2470
2471 if( !Data("classes") || !Data("features") || !Data["features"]("count") || Data["features"]["count"].Get_Content().asInt() != m_nFeatures || m_nFeatures == 0 )
2472 {
2473 return( false );
2474 }
2475
2476 if( Data["features"]("info") )
2477 {
2478 m_Info = Data["features"]["info"].Get_Content();
2479 }
2480
2481 //-----------------------------------------------------
2482 CSG_MetaData &Classes = *Data.Get_Child("CLASSES");
2483
2484 for(int i=0; i<Classes.Get_Children_Count(); i++)
2485 {
2486 if( Classes[i].Cmp_Name("class") && Classes[i].Get_Child("id") )
2487 {
2488 bool bAdd = true;
2489
2490 CClass *pClass = new CClass(Classes[i]["id"].Get_Content());
2491
2492 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; }
2493 if( !pClass->m_Mean.from_String(Classes[i]["mean"].Get_Content()) || pClass->m_Mean.Get_N () != m_nFeatures ) { bAdd = false; }
2494 if( !pClass->m_Min .from_String(Classes[i]["min" ].Get_Content()) || pClass->m_Min .Get_N () != m_nFeatures ) { bAdd = false; }
2495 if( !pClass->m_Max .from_String(Classes[i]["max" ].Get_Content()) || pClass->m_Max .Get_N () != m_nFeatures ) { bAdd = false; }
2496
2497 //---------------------------------------------
2498 if( !bAdd )
2499 {
2500 delete(pClass);
2501 }
2502 else
2503 {
2504 m_pClasses = (CClass **)SG_Realloc(m_pClasses, (m_nClasses + 1) * sizeof(CClass *));
2505 m_pClasses[m_nClasses++] = pClass;
2506
2507 pClass->m_Cov_Det = pClass->m_Cov.Get_Determinant();
2508 pClass->m_Cov_Inv = pClass->m_Cov.Get_Inverse();
2509
2510 pClass->m_Mean_Spectral = CSG_Simple_Statistics(pClass->m_Mean).Get_Mean();
2511 }
2512 }
2513 }
2514
2515 return( m_nClasses > 0 );
2516 }
2517
2518 //---------------------------------------------------------
Save(const CSG_String & File,const SG_Char * Feature_Info)2519 bool CSG_Classifier_Supervised::Save(const CSG_String &File, const SG_Char *Feature_Info)
2520 {
2521 if( m_nFeatures < 1 || m_nClasses < 1 || File.is_Empty() )
2522 {
2523 return( false );
2524 }
2525
2526 CSG_MetaData Data;
2527
2528 Data.Set_Name ("supervised_classifier");
2529 Data.Add_Property("saga-version", SAGA_VERSION);
2530
2531 CSG_MetaData &Features = *Data.Add_Child("features");
2532
2533 Features.Add_Child("count", m_nFeatures);
2534
2535 if( Feature_Info && *Feature_Info )
2536 {
2537 Features.Add_Child("info", Feature_Info);
2538 }
2539
2540 CSG_MetaData &Classes = *Data.Add_Child("classes");
2541
2542 Classes.Add_Property("count", m_nClasses);
2543
2544 for(int i=0; i<m_nClasses; i++)
2545 {
2546 CSG_MetaData &Class = *Classes.Add_Child("class");
2547
2548 CClass *pClass = m_pClasses[i];
2549
2550 Class.Add_Child("id" , pClass->m_ID );
2551 Class.Add_Child("mean", pClass->m_Mean.to_String());
2552 Class.Add_Child("min" , pClass->m_Min .to_String());
2553 Class.Add_Child("max" , pClass->m_Max .to_String());
2554 Class.Add_Child("cov" , pClass->m_Cov .to_String());
2555 }
2556
2557 return( Data.Save(File) );
2558 }
2559
2560
2561 ///////////////////////////////////////////////////////////
2562 // //
2563 ///////////////////////////////////////////////////////////
2564
2565 //---------------------------------------------------------
Print(void)2566 CSG_String CSG_Classifier_Supervised::Print(void)
2567 {
2568 CSG_String s;
2569
2570 if( m_nFeatures > 0 && m_nClasses > 0 )
2571 {
2572 s += "\n";
2573
2574 for(int iClass=0; iClass<m_nClasses; iClass++)
2575 {
2576 CClass *pClass = m_pClasses[iClass];
2577
2578 s += "\n____\n" + pClass->m_ID + "\nFeature\tMean\tMin\tMax\tStdDev";
2579
2580 for(int i=0; i<m_nFeatures; i++)
2581 {
2582 s += CSG_String::Format("\n%3d.", i + 1);
2583 s += "\t" + SG_Get_String(pClass->m_Mean[i]);
2584 s += "\t" + SG_Get_String(pClass->m_Min [i]);
2585 s += "\t" + SG_Get_String(pClass->m_Max [i]);
2586 s += "\t" + SG_Get_String(sqrt(pClass->m_Cov[i][i]));
2587 }
2588
2589 s += "\n";
2590 }
2591 }
2592
2593 return( s );
2594 }
2595
2596
2597 ///////////////////////////////////////////////////////////
2598 // //
2599 ///////////////////////////////////////////////////////////
2600
2601 //---------------------------------------------------------
Add_Class(const CSG_String & Class_ID,const CSG_Vector & Mean,const CSG_Vector & Min,const CSG_Vector & Max,const CSG_Matrix & Cov)2602 bool 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)
2603 {
2604 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 )
2605 {
2606 return( false );
2607 }
2608
2609 CClass *pClass, **pClasses = (CClass **)SG_Realloc(m_pClasses, (m_nClasses + 1) * sizeof(CClass *));
2610
2611 if( pClasses )
2612 {
2613 m_pClasses = pClasses;
2614
2615 m_pClasses[m_nClasses++] = pClass = new CClass(Class_ID);
2616
2617 pClass->m_ID = Class_ID;
2618
2619 pClass->m_Mean = Mean;
2620 pClass->m_Min = Min;
2621 pClass->m_Max = Max;
2622 pClass->m_Cov = Cov;
2623
2624 pClass->m_Cov_Inv = Cov.Get_Inverse();
2625 pClass->m_Cov_Det = Cov.Get_Determinant();
2626
2627 pClass->m_Mean_Spectral = CSG_Simple_Statistics(Mean).Get_Mean();
2628
2629 return( true );
2630 }
2631
2632 return( false );
2633 }
2634
2635
2636 ///////////////////////////////////////////////////////////
2637 // //
2638 ///////////////////////////////////////////////////////////
2639
2640 //---------------------------------------------------------
Train_Clr_Samples(void)2641 bool CSG_Classifier_Supervised::Train_Clr_Samples(void)
2642 {
2643 for(int i=0; i<m_nClasses; i++)
2644 {
2645 m_pClasses[i]->m_Samples.Destroy();
2646 }
2647
2648 return( true );
2649 }
2650
2651 //---------------------------------------------------------
Train_Add_Sample(const CSG_String & Class_ID,const CSG_Vector & Features)2652 bool CSG_Classifier_Supervised::Train_Add_Sample(const CSG_String &Class_ID, const CSG_Vector &Features)
2653 {
2654 if( m_nFeatures > 0 && m_nFeatures == Features.Get_N() )
2655 {
2656 int iClass = Get_Class(Class_ID);
2657
2658 if( iClass < 0 )
2659 {
2660 CClass **pClasses = (CClass **)SG_Realloc(m_pClasses, (m_nClasses + 1) * sizeof(CClass *));
2661
2662 if( pClasses )
2663 {
2664 m_pClasses = pClasses;
2665
2666 m_pClasses[iClass = m_nClasses++] = new CClass(Class_ID);
2667 }
2668 }
2669
2670 if( iClass >= 0 )
2671 {
2672 return( m_pClasses[iClass]->m_Samples.Add_Row(Features) );
2673 }
2674 }
2675
2676 return( false );
2677 }
2678
2679 //---------------------------------------------------------
Train(bool bClear_Samples)2680 bool CSG_Classifier_Supervised::Train(bool bClear_Samples)
2681 {
2682 if( m_nFeatures < 1 || m_nClasses < 1 )
2683 {
2684 return( false );
2685 }
2686
2687 for(int iClass=0; iClass<m_nClasses; iClass++)
2688 {
2689 if( !m_pClasses[iClass]->Train() )
2690 {
2691 return( false );
2692 }
2693 }
2694
2695 if( bClear_Samples )
2696 {
2697 Train_Clr_Samples();
2698 }
2699
2700 return( true );
2701 }
2702
2703
2704 ///////////////////////////////////////////////////////////
2705 // //
2706 ///////////////////////////////////////////////////////////
2707
2708 //---------------------------------------------------------
Train(void)2709 bool CSG_Classifier_Supervised::CClass::Train(void)
2710 {
2711 if( m_Samples.Get_NCols() < 1 || m_Samples.Get_NRows() < 1 )
2712 {
2713 return( false );
2714 }
2715
2716 int iFeature;
2717
2718 //-----------------------------------------------------
2719 m_Mean.Create(m_Samples.Get_NCols());
2720 m_Min .Create(m_Samples.Get_NCols());
2721 m_Max .Create(m_Samples.Get_NCols());
2722
2723 for(iFeature=0; iFeature<m_Samples.Get_NCols(); iFeature++)
2724 {
2725 CSG_Simple_Statistics s;
2726
2727 for(int iSample=0; iSample<m_Samples.Get_NRows(); iSample++)
2728 {
2729 s += m_Samples[iSample][iFeature];
2730 }
2731
2732 m_Mean[iFeature] = s.Get_Mean ();
2733 m_Min [iFeature] = s.Get_Minimum();
2734 m_Max [iFeature] = s.Get_Maximum();
2735 }
2736
2737 //-----------------------------------------------------
2738 m_Cov.Create(m_Samples.Get_NCols(), m_Samples.Get_NCols());
2739
2740 for(iFeature=0; iFeature<m_Samples.Get_NCols(); iFeature++)
2741 {
2742 for(int jFeature=iFeature; jFeature<m_Samples.Get_NCols(); jFeature++)
2743 {
2744 double cov = 0.;
2745
2746 for(int iSample=0; iSample<m_Samples.Get_NRows(); iSample++)
2747 {
2748 cov += (m_Samples[iSample][iFeature] - m_Mean[iFeature]) * (m_Samples[iSample][jFeature] - m_Mean[jFeature]);
2749 }
2750
2751 if( m_Samples.Get_NRows() > 1 )
2752 {
2753 cov /= m_Samples.Get_NRows() - 1;
2754 }
2755
2756 m_Cov[iFeature][jFeature] = m_Cov[jFeature][iFeature] = cov;
2757 }
2758 }
2759
2760 m_Cov_Inv = m_Cov.Get_Inverse ();
2761 m_Cov_Det = m_Cov.Get_Determinant();
2762
2763 m_Mean_Spectral = CSG_Simple_Statistics(m_Mean).Get_Mean();
2764
2765 //-----------------------------------------------------
2766 return( true );
2767 }
2768
2769
2770 ///////////////////////////////////////////////////////////
2771 // //
2772 ///////////////////////////////////////////////////////////
2773
2774 //---------------------------------------------------------
Get_Class(const CSG_String & Class_ID)2775 int CSG_Classifier_Supervised::Get_Class(const CSG_String &Class_ID)
2776 {
2777 if( m_nFeatures > 0 )
2778 {
2779 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2780 {
2781 if( !Get_Class_ID(iClass).Cmp(Class_ID) )
2782 {
2783 return( iClass );
2784 }
2785 }
2786 }
2787
2788 return( -1 );
2789 }
2790
2791 //---------------------------------------------------------
Get_Class(const CSG_Vector & Features,int & Class,double & Quality,int Method)2792 bool CSG_Classifier_Supervised::Get_Class(const CSG_Vector &Features, int &Class, double &Quality, int Method)
2793 {
2794 Class = -1;
2795 Quality = 0.;
2796
2797 if( Get_Feature_Count() == Features.Get_N() )
2798 {
2799 switch( Method )
2800 {
2801 case SG_CLASSIFY_SUPERVISED_BinaryEncoding : _Get_Binary_Encoding (Features, Class, Quality); break;
2802 case SG_CLASSIFY_SUPERVISED_ParallelEpiped : _Get_Parallel_Epiped (Features, Class, Quality); break;
2803 case SG_CLASSIFY_SUPERVISED_MinimumDistance : _Get_Minimum_Distance (Features, Class, Quality); break;
2804 case SG_CLASSIFY_SUPERVISED_Mahalonobis : _Get_Mahalanobis_Distance (Features, Class, Quality); break;
2805 case SG_CLASSIFY_SUPERVISED_MaximumLikelihood: _Get_Maximum_Likelihood (Features, Class, Quality); break;
2806 case SG_CLASSIFY_SUPERVISED_SAM : _Get_Spectral_Angle_Mapping(Features, Class, Quality); break;
2807 case SG_CLASSIFY_SUPERVISED_SID : _Get_Spectral_Divergence (Features, Class, Quality); break;
2808 case SG_CLASSIFY_SUPERVISED_WTA : _Get_Winner_Takes_All (Features, Class, Quality); break;
2809 }
2810
2811 return( Class >= 0 );
2812 }
2813
2814 return( false );
2815 }
2816
2817
2818 ///////////////////////////////////////////////////////////
2819 // //
2820 ///////////////////////////////////////////////////////////
2821
2822 //---------------------------------------------------------
Get_Name_of_Method(int Method)2823 CSG_String CSG_Classifier_Supervised::Get_Name_of_Method(int Method)
2824 {
2825 switch( Method )
2826 {
2827 case SG_CLASSIFY_SUPERVISED_BinaryEncoding : return( _TL("Binary Encoding") );
2828 case SG_CLASSIFY_SUPERVISED_ParallelEpiped : return( _TL("Parallelepiped") );
2829 case SG_CLASSIFY_SUPERVISED_MinimumDistance : return( _TL("Minimum Distance") );
2830 case SG_CLASSIFY_SUPERVISED_Mahalonobis : return( _TL("Mahalanobis Distance") );
2831 case SG_CLASSIFY_SUPERVISED_MaximumLikelihood: return( _TL("Maximum Likelihood") );
2832 case SG_CLASSIFY_SUPERVISED_SAM : return( _TL("Spectral Angle Mapping") );
2833 case SG_CLASSIFY_SUPERVISED_SID : return( _TL("Spectral Information Divergence") );
2834 case SG_CLASSIFY_SUPERVISED_SVM : return( _TL("Support Vector Machine") );
2835 case SG_CLASSIFY_SUPERVISED_WTA : return( _TL("Winner Takes All") );
2836 }
2837
2838 return( SG_T("") );
2839 }
2840
2841 //---------------------------------------------------------
Get_Name_of_Quality(int Method)2842 CSG_String CSG_Classifier_Supervised::Get_Name_of_Quality(int Method)
2843 {
2844 switch( Method )
2845 {
2846 case SG_CLASSIFY_SUPERVISED_BinaryEncoding : return( _TL("Difference") );
2847 case SG_CLASSIFY_SUPERVISED_ParallelEpiped : return( _TL("Memberships") );
2848 case SG_CLASSIFY_SUPERVISED_MinimumDistance : return( _TL("Distance") );
2849 case SG_CLASSIFY_SUPERVISED_Mahalonobis : return( _TL("Distance") );
2850 case SG_CLASSIFY_SUPERVISED_MaximumLikelihood: return( _TL("Proximity") );
2851 case SG_CLASSIFY_SUPERVISED_SAM : return( _TL("Angle") );
2852 case SG_CLASSIFY_SUPERVISED_SID : return( _TL("Divergence") );
2853 case SG_CLASSIFY_SUPERVISED_SVM : return( _TL("") );
2854 case SG_CLASSIFY_SUPERVISED_WTA : return( _TL("Votes") );
2855 }
2856
2857 return( SG_T("") );
2858 }
2859
2860
2861 ///////////////////////////////////////////////////////////
2862 // //
2863 ///////////////////////////////////////////////////////////
2864
2865 //---------------------------------------------------------
2866 // Mazer, A. S., Martin, M., Lee, M., and Solomon, J. E. (1988):
2867 // Image Processing Software for Imaging Spectrometry Analysis.
2868 // Remote Sensing of Environment, v. 24, no. 1, p. 201-210.
2869 //
_Get_Binary_Encoding(const CSG_Vector & Features,int & Class,double & Quality)2870 void CSG_Classifier_Supervised::_Get_Binary_Encoding(const CSG_Vector &Features, int &Class, double &Quality)
2871 {
2872 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2873 {
2874 CClass *pClass = m_pClasses[iClass];
2875
2876 double Mean_Spectral = CSG_Simple_Statistics(Features).Get_Mean();
2877
2878 int d = 0;
2879
2880 for(int iFeature=0; iFeature<Get_Feature_Count(); iFeature++)
2881 {
2882 d += (Features(iFeature) < Mean_Spectral) == (pClass->m_Mean[iFeature] < pClass->m_Mean_Spectral) ? 0 : 1;
2883
2884 if( iFeature == 0 ) // spectral slopes
2885 {
2886 d += (Features[iFeature ] < Features[iFeature + 1]) == (pClass->m_Mean[iFeature ] < pClass->m_Mean[iFeature + 1]) ? 0 : 1;
2887 }
2888 else if( iFeature == Get_Feature_Count() - 1 )
2889 {
2890 d += (Features[iFeature - 1] < Features[iFeature ]) == (pClass->m_Mean[iFeature - 1] < pClass->m_Mean[iFeature ]) ? 0 : 1;
2891 }
2892 else
2893 {
2894 d += (Features[iFeature - 1] < Features[iFeature + 1]) == (pClass->m_Mean[iFeature - 1] < pClass->m_Mean[iFeature + 1]) ? 0 : 1;
2895 }
2896 }
2897
2898 if( Class < 0 || Quality > d ) // find the minimum 'Hamming' distance
2899 {
2900 Quality = d;
2901 Class = iClass;
2902 }
2903 }
2904 }
2905
2906 //---------------------------------------------------------
_Get_Parallel_Epiped(const CSG_Vector & Features,int & Class,double & Quality)2907 void CSG_Classifier_Supervised::_Get_Parallel_Epiped(const CSG_Vector &Features, int &Class, double &Quality)
2908 {
2909 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2910 {
2911 CClass *pClass = m_pClasses[iClass];
2912
2913 bool bMember = true;
2914
2915 for(int iFeature=0; bMember && iFeature<Get_Feature_Count(); iFeature++)
2916 {
2917 bMember = pClass->m_Min[iFeature] <= Features[iFeature] && Features[iFeature] <= pClass->m_Max[iFeature];
2918 }
2919
2920 if( bMember )
2921 {
2922 Quality ++;
2923 Class = iClass;
2924 }
2925 }
2926 }
2927
2928 //---------------------------------------------------------
_Get_Minimum_Distance(const CSG_Vector & Features,int & Class,double & Quality)2929 void CSG_Classifier_Supervised::_Get_Minimum_Distance(const CSG_Vector &Features, int &Class, double &Quality)
2930 {
2931 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2932 {
2933 CClass *pClass = m_pClasses[iClass];
2934
2935 double Distance = (Features - pClass->m_Mean).Get_Length();
2936
2937 if( Class < 0 || Quality > Distance )
2938 {
2939 Quality = Distance;
2940 Class = iClass;
2941 }
2942 }
2943
2944 if( m_Threshold_Distance > 0. && Quality > m_Threshold_Distance )
2945 {
2946 Class = -1;
2947 }
2948 }
2949
2950 //---------------------------------------------------------
_Get_Mahalanobis_Distance(const CSG_Vector & Features,int & Class,double & Quality)2951 void CSG_Classifier_Supervised::_Get_Mahalanobis_Distance(const CSG_Vector &Features, int &Class, double &Quality)
2952 {
2953 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2954 {
2955 CClass *pClass = m_pClasses[iClass];
2956
2957 CSG_Vector D = Features - pClass->m_Mean;
2958
2959 double Distance = D * (pClass->m_Cov_Inv * D);
2960
2961 if( Class < 0 || Quality > Distance )
2962 {
2963 Quality = Distance;
2964 Class = iClass;
2965 }
2966 }
2967
2968 if( m_Threshold_Distance > 0. && Quality > m_Threshold_Distance )
2969 {
2970 Class = -1;
2971 }
2972 }
2973
2974 //---------------------------------------------------------
_Get_Maximum_Likelihood(const CSG_Vector & Features,int & Class,double & Quality)2975 void CSG_Classifier_Supervised::_Get_Maximum_Likelihood(const CSG_Vector &Features, int &Class, double &Quality)
2976 {
2977 double dSum = 0.;
2978
2979 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2980 {
2981 CClass *pClass = m_pClasses[iClass];
2982
2983 CSG_Vector D = Features - pClass->m_Mean;
2984
2985 double Distance = D * (pClass->m_Cov_Inv * D);
2986
2987 double Probability = pow(2. * M_PI, -0.5 * m_nFeatures) * pow(pClass->m_Cov_Det, -0.5) * exp(-0.5 * Distance);
2988 // double Probability = -log(pClass->m_Cov_Det) - Distance;
2989
2990 dSum += Probability;
2991
2992 if( Class < 0 || Quality < Probability )
2993 {
2994 Quality = Probability;
2995 Class = iClass;
2996 }
2997 }
2998
2999 if( Class >= 0 )
3000 {
3001 if( m_Probability_Relative )
3002 {
3003 Quality = 100. * Quality / dSum;
3004 }
3005
3006 if( m_Threshold_Probability > 0. && Quality < m_Threshold_Probability )
3007 {
3008 Class = -1;
3009 }
3010 }
3011 }
3012
3013 //---------------------------------------------------------
_Get_Spectral_Angle_Mapping(const CSG_Vector & Features,int & Class,double & Quality)3014 void CSG_Classifier_Supervised::_Get_Spectral_Angle_Mapping(const CSG_Vector &Features, int &Class, double &Quality)
3015 {
3016 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
3017 {
3018 CClass *pClass = m_pClasses[iClass];
3019
3020 double Angle = Features.Get_Angle(pClass->m_Mean);
3021
3022 if( Class < 0 || Quality > Angle )
3023 {
3024 Quality = Angle;
3025 Class = iClass;
3026 }
3027 }
3028
3029 Quality *= M_RAD_TO_DEG;
3030
3031 if( m_Threshold_Angle > 0. && Quality > m_Threshold_Angle )
3032 {
3033 Class = -1;
3034 }
3035 }
3036
3037 //---------------------------------------------------------
_Get_Spectral_Divergence(const CSG_Vector & Features,int & Class,double & Quality)3038 void CSG_Classifier_Supervised::_Get_Spectral_Divergence(const CSG_Vector &Features, int &Class, double &Quality)
3039 {
3040 }
3041
3042 //---------------------------------------------------------
_Get_Winner_Takes_All(const CSG_Vector & Features,int & Class,double & Quality)3043 void CSG_Classifier_Supervised::_Get_Winner_Takes_All(const CSG_Vector &Features, int &Class, double &Quality)
3044 {
3045 int *Votes = (int *)SG_Calloc(Get_Class_Count(), sizeof(int));
3046
3047 for(int iMethod=0; iMethod<SG_CLASSIFY_SUPERVISED_WTA; iMethod++)
3048 {
3049 int iClass;
3050 double iQuality;
3051
3052 if( m_bWTA[iMethod] && Get_Class(Features, iClass, iQuality, iMethod) && ++Votes[iClass] > Quality )
3053 {
3054 Quality = Votes[iClass];
3055 Class = iClass;
3056 }
3057 }
3058
3059 SG_Free(Votes);
3060 }
3061
3062
3063 ///////////////////////////////////////////////////////////
3064 // //
3065 // //
3066 // //
3067 ///////////////////////////////////////////////////////////
3068
3069 //---------------------------------------------------------
3070 // source: http://psydok.sulb.uni-saarland.de/volltexte/2004/268/html/
3071
3072 //---------------------------------------------------------
Get_T_Tail(double T,int df,TSG_Test_Distribution_Type Type)3073 double CSG_Test_Distribution::Get_T_Tail(double T, int df, TSG_Test_Distribution_Type Type)
3074 { // Hill's approx. to cumulative t-dist, Commun.A.C.M. 13,617-619.
3075 // See: J.H.Maindonald, Computational Statistics, p.295.
3076 // Calculates p given t and tail type.
3077
3078 if( !T || !df || df < 1. )
3079 {
3080 return( -1. );
3081 }
3082
3083 return( _Change_Tail_Type(Get_T_P(T, df), TESTDIST_TYPE_TwoTail, Type, T < 0.) );
3084 }
3085
3086 //---------------------------------------------------------
Get_T_Inverse(double p,int df,TSG_Test_Distribution_Type Type)3087 double CSG_Test_Distribution::Get_T_Inverse(double p, int df, TSG_Test_Distribution_Type Type)
3088 { // Keith Dear & Robert Brennan.
3089 // Returns an accurate t to tol sig. fig.'s given p & df.
3090
3091 if( p <= 0. || p >= 1. || df < 1 )
3092 {
3093 return( -1. );
3094 }
3095
3096 bool bNegative = (Type == TESTDIST_TYPE_Left && p < 0.5) || (Type == TESTDIST_TYPE_Right && p > 0.5);
3097 double t, p0, p1, diff;
3098
3099 p0 = p1 = _Change_Tail_Type(p, Type, TESTDIST_TYPE_TwoTail, bNegative);
3100 diff = 1.;
3101
3102 while( fabs(diff) > 0.0001 )
3103 {
3104 t = Get_T_Inv(p1, df); // initial rough value
3105 diff = Get_T_P(t, df) - p0; // compare result with forward fn
3106 p1 = p1 - diff; // small adjustment to p1
3107 }
3108
3109 return( bNegative ? -t : t );
3110 }
3111
3112 //---------------------------------------------------------
_Change_Tail_Type(double p,TSG_Test_Distribution_Type from,TSG_Test_Distribution_Type to,bool bNegative)3113 double CSG_Test_Distribution::_Change_Tail_Type(double p, TSG_Test_Distribution_Type from, TSG_Test_Distribution_Type to, bool bNegative)
3114 {
3115 if( from != to )
3116 {
3117 switch( from ) // convert any tail type to 'left'
3118 {
3119 case TESTDIST_TYPE_Left : break;
3120 case TESTDIST_TYPE_Right : p = 1. - p; break;
3121 case TESTDIST_TYPE_Middle : p = p / 2. + 0.5; if( bNegative ) p = 1. - p; break;
3122 case TESTDIST_TYPE_TwoTail: p = 1. - p / 2. ; if( bNegative ) p = 1. - p; break;
3123 // case TESTDIST_TYPE_Half : p = p + 0.5 ; if( bNegative ) p = 1. - p; break;
3124 }
3125
3126 switch( to ) // convert p from tail type 'left' to any other
3127 {
3128 case TESTDIST_TYPE_Left : break;
3129 case TESTDIST_TYPE_Right : p = 1. - p; break;
3130 case TESTDIST_TYPE_Middle : if( bNegative ) p = 1. - p; p = 2. * (1. - p); break;
3131 case TESTDIST_TYPE_TwoTail: if( bNegative ) p = 1. - p; p = 2. * p - 1. ; break;
3132 // case TESTDIST_TYPE_Half : if( bNegative ) p = 1. - p; p = p - 0.5 ; break;
3133 }
3134 }
3135
3136 return( p );
3137 }
3138
3139 //---------------------------------------------------------
Get_Norm_P(double z)3140 double CSG_Test_Distribution::Get_Norm_P(double z)
3141 { // Returns the two-tailed standard normal probability of z
3142 const double a1 = 0.0000053830, a2 = 0.0000488906, a3 = 0.0000380036,
3143 a4 = 0.0032776263, a5 = 0.0211410061, a6 = 0.0498673470;
3144
3145 double p;
3146
3147 z = fabs(z);
3148
3149 p = (((((a1 * z + a2) * z + a3) * z + a4) * z + a5) * z + a6) * z + 1.;
3150
3151 return( pow(p, -16) );
3152 }
3153
3154 //---------------------------------------------------------
Get_Norm_Z(double p)3155 double CSG_Test_Distribution::Get_Norm_Z(double p)
3156 { // Returns z given a half-middle tail type p.
3157 const double a0 = 2.5066282, a1 = -18.6150006, a2 = 41.3911977, a3 = -25.4410605,
3158 b1 = -8.4735109, b2 = 23.0833674, b3 = -21.0622410, b4 = 3.1308291,
3159 c0 = -2.7871893, c1 = -2.2979648, c2 = 4.8501413, c3 = 2.3212128,
3160 d1 = 3.5438892, d2 = 1.6370678;
3161
3162 double r, z;
3163
3164 if( p > 0.42 )
3165 {
3166 r = sqrt(-log(0.5 - p));
3167 z = (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + 1.);
3168 }
3169 else
3170 {
3171 r = p * p;
3172 z = p * (((a3 * r + a2) * r + a1) * r + a0) / ((((b4 * r + b3) * r + b2) * r + b1) * r + 1.);
3173 }
3174
3175 return( z );
3176 }
3177
3178 //---------------------------------------------------------
Get_T_P(double T,int df)3179 double CSG_Test_Distribution::Get_T_P(double T, int df)
3180 { // Returns two-tail probability level given t and df.
3181 return( df == 1 ? 1. - 2. * atan(fabs(T)) / M_PI
3182 : df == 2 ? 1. - fabs(T) / sqrt(T*T + 2.)
3183 : df == 3 ? 1. - 2. * (atan(fabs(T) / sqrt(3.)) + fabs(T) * sqrt(3.) / (T*T + 3.)) / M_PI
3184 : df == 4 ? 1. - fabs(T) * (1. + 2. / (T*T + 4.)) / sqrt(T*T + 4.)
3185 : Get_Norm_P(Get_T_Z(fabs(T), df))
3186 );
3187 }
3188
3189 //---------------------------------------------------------
Get_T_Z(double T,int df)3190 double CSG_Test_Distribution::Get_T_Z(double T, int df)
3191 { // Converts a t value to an approximate z value w.r.t the given df
3192 // s.t. std.norm.(z) = t(z, df) at the two-tail probability level.
3193
3194 double A9, B9, T9, Z8, P7, B7, z;
3195
3196 A9 = df - 0.5;
3197 B9 = 48. * A9*A9,
3198 T9 = T*T / df;
3199 Z8 = T9 >= 0.04
3200 ? A9 * log(1. + T9)
3201 : A9 * (((1. - T9 * 0.75) * T9 / 3. - 0.5) * T9 + 1.) * T9;
3202 P7 = ((0.4 * Z8 + 3.3) * Z8 + 24.) * Z8 + 85.5;
3203 B7 = 0.8 * pow(Z8, 2.) + 100. + B9;
3204 z = (1. + (-P7 / B7 + Z8 + 3.) / B9) * sqrt(Z8);
3205
3206 return( z );
3207 }
3208
3209 //---------------------------------------------------------
Get_T_Inv(double p,int df)3210 double CSG_Test_Distribution::Get_T_Inv(double p, int df)
3211 { // Hill's approx. inverse t-dist.: Comm. of A.C.M Vol.13 No.10 1970 pg 620.
3212 // Calculates t given df and two-tail probability.
3213
3214 if( df == 1 )
3215 {
3216 return( cos(p * M_PI / 2.) / sin(p * M_PI / 2.) );
3217 }
3218
3219 if( df == 2 )
3220 {
3221 return( sqrt(2. / (p * (2. - p)) - 2.) );
3222 }
3223
3224 double a, b, c, d, x, y;
3225
3226 a = 1. / (df - 0.5);
3227 b = 48. / (a*a);
3228 c = ((20700. * a / b - 98.) * a - 16.) * a + 96.36;
3229 d = ((94.5 / (b + c) - 3.) / b + 1.) * sqrt(a * M_PI / 2.) * df;
3230 x = d * p;
3231 y = pow(x, 2. / df);
3232
3233 if( y > 0.05 + a )
3234 {
3235 x = Get_Norm_Z(0.5 * (1. - p));
3236 y = x*x;
3237
3238 if( df < 5 )
3239 {
3240 c = c + 0.3 * (df - 4.5) * (x + 0.6);
3241 }
3242
3243 c = (((0.05 * d * x - 5) * x - 7.) * x - 2.) * x + b + c;
3244 y = (((((0.4 * y + 6.3) * y + 36.) * y + 94.5) / c - y - 3.) / b + 1.) * x;
3245 y = a * y*y;
3246
3247 if( y > 0.002 )
3248 {
3249 y = exp(y) - 1.;
3250 }
3251 else
3252 {
3253 y = 0.5 * y*y + y;
3254 }
3255 }
3256 else
3257 {
3258 y = ((1. / (((df + 6.) / (df * y) - 0.089 * d - 0.822) * (df + 2.) * 3.)
3259 + 0.5 / (df + 4.)) * y - 1.) * (df + 1.) / (df + 2.) + 1. / y;
3260 }
3261
3262 return( sqrt(df * y) );
3263 }
3264
3265
3266 ///////////////////////////////////////////////////////////
3267 // //
3268 ///////////////////////////////////////////////////////////
3269
3270 //---------------------------------------------------------
Get_F_Tail_from_R2(double R2,int nPredictors,int nSamples,TSG_Test_Distribution_Type Type)3271 double CSG_Test_Distribution::Get_F_Tail_from_R2(double R2, int nPredictors, int nSamples, TSG_Test_Distribution_Type Type)
3272 {
3273 double F = (nSamples - nPredictors - 1) * (R2 / nPredictors) / (1. - R2);
3274
3275 return( CSG_Test_Distribution::Get_F_Tail(F, nPredictors, nSamples - nPredictors - 1, Type) );
3276 }
3277
3278 //---------------------------------------------------------
Get_F_Tail(double F,int dfn,int dfd,TSG_Test_Distribution_Type Type)3279 double CSG_Test_Distribution::Get_F_Tail(double F, int dfn, int dfd, TSG_Test_Distribution_Type Type)
3280 {
3281 // calculates for F, dfn(ominator) and dfd(enominator) the "tail" of the F-distribution
3282
3283 double p = 1.;
3284
3285 if( F >= 0.00001 && dfn > 0 && dfd > 0 )
3286 {
3287 if( F * dfn >= dfd || F > 1. + 20. / dfn + 10. / sqrt((double)dfn) )
3288 {
3289 p = Get_Gamma(F, dfn, dfd);
3290 }
3291 else
3292 {
3293 p = 1. - Get_Gamma(1. / F, dfd, dfn);
3294 }
3295 }
3296
3297 if( p <= 0. || p >= 1. )
3298 {
3299 p = F > 1. ? 0. : F < 1. ? 1. : 0.5;
3300 }
3301
3302 return( Type == TESTDIST_TYPE_Right ? p : 1. - p );
3303 }
3304
3305 //---------------------------------------------------------
Get_F_Inverse(double alpha,int dfn,int dfd,TSG_Test_Distribution_Type Type)3306 double CSG_Test_Distribution::Get_F_Inverse(double alpha, int dfn, int dfd, TSG_Test_Distribution_Type Type)
3307 {
3308 if( alpha < 0. || alpha > 1. || dfd < 0 || dfn < 0 )
3309 {
3310 return( -1 );
3311 }
3312
3313 if( Type != TESTDIST_TYPE_Right )
3314 {
3315 alpha = 1. - alpha;
3316 }
3317
3318 const int ITERMAX = 100;
3319 const double EPSILON = 0.0001;
3320
3321 int i;
3322 double lo, hi, mid, p;
3323
3324 if( alpha <= 0.5 )
3325 {
3326 lo = 0.5;
3327 hi = lo;
3328
3329 for(i=0; i<ITERMAX; i++)
3330 {
3331 hi *= 2.;
3332 p = Get_F_Tail(hi, dfn, dfd);
3333
3334 if( p > alpha )
3335 {
3336 lo = hi;
3337 }
3338 else
3339 {
3340 break;
3341 }
3342 }
3343
3344 if( p > alpha )
3345 {
3346 return( hi );
3347 }
3348 }
3349 else
3350 {
3351 hi = 2;
3352 lo = hi;
3353
3354 for(i=0; i<ITERMAX; i++)
3355 {
3356 lo /= 2.;
3357 p = Get_F_Tail(lo, dfn, dfd);
3358
3359 if( p < alpha )
3360 {
3361 hi = lo;
3362 }
3363 else
3364 {
3365 break;
3366 }
3367 }
3368
3369 if( p < alpha )
3370 {
3371 return( lo );
3372 }
3373 }
3374
3375 mid = (hi + lo) / 2.;
3376
3377 for(i=0; i<ITERMAX && (hi-lo)>EPSILON*mid; i++)
3378 {
3379 mid = (hi + lo) / 2.;
3380 p = Get_F_Tail(mid, dfn, dfd);
3381
3382 if( p < alpha )
3383 hi = mid;
3384 else if( p > alpha )
3385 lo = mid;
3386 else
3387 break;
3388 }
3389
3390 return( mid );
3391 }
3392
3393 //---------------------------------------------------------
Get_Gamma(double F,double dfn,double dfd)3394 double CSG_Test_Distribution::Get_Gamma(double F, double dfn, double dfd)
3395 {
3396 // calculates for F, dfn(ominator) and dfd(enominator) the incomplete Gamma-function
3397
3398 const double EXPMIN = -30.;
3399 const double SMALL = 0.00000000001;
3400
3401 double x, c, er, s, n, t1, t;
3402
3403 dfn /= 2.;
3404 dfd /= 2.;
3405
3406 x = dfd / (dfd + dfn * F);
3407 c = Get_Log_Gamma(dfn + dfd) - Get_Log_Gamma(dfn) - Get_Log_Gamma(dfd + 1.) + dfd * log(x) + dfn * log(1. - x);
3408
3409 if( c < EXPMIN )
3410 {
3411 return( -1. );
3412 }
3413
3414 dfn += dfd;
3415 dfd += 1.;
3416 c = exp(c);
3417 er = SMALL / c;
3418 t = dfn * x / dfd;
3419 t1 = 0.;
3420 s = t + 1.;
3421 n = 0;
3422
3423 while( t > er || t > t1 )
3424 {
3425 n += 1;
3426 t1 = t;
3427 t *= ((dfn + n) * x / (dfd + n));
3428 s += t;
3429 }
3430
3431 return( s * c );
3432 }
3433
3434 //---------------------------------------------------------
Get_Log_Gamma(double a)3435 double CSG_Test_Distribution::Get_Log_Gamma(double a)
3436 {
3437 // calculates the logarithm of the Gamma-function
3438
3439 const int ARGMIN = 6;
3440
3441 const double HL2PI = 0.91893853320467275; // = log(2. * M_PI) / 2.
3442
3443 int n = (int)floor(ARGMIN - a + 0.0001);
3444
3445 if( n > 0 )
3446 {
3447 a += n;
3448 }
3449
3450 double g;
3451
3452 g = 1. / (a*a);
3453 g = (1. - g * (1. / 30. - g * (1. / 105. - g * (1. / 140. - g / 99.)))) / (12. * a);
3454 g = g + ((a - 0.5) * log(a) - a + HL2PI);
3455
3456 for(int i=0; i<n; i++)
3457 {
3458 a = a - 1.;
3459 g = g - log(a);
3460 }
3461
3462 return( g );
3463 }
3464
3465
3466 ///////////////////////////////////////////////////////////
3467 // //
3468 // //
3469 // //
3470 ///////////////////////////////////////////////////////////
3471
3472 //---------------------------------------------------------
SG_Get_Correlation_Matrix(const CSG_Matrix & Values,bool bCovariances)3473 CSG_Matrix SG_Get_Correlation_Matrix (const CSG_Matrix &Values, bool bCovariances)
3474 {
3475 int nVariables = Values.Get_NX();
3476 int nSamples = Values.Get_NY();
3477
3478 //-----------------------------------------------------
3479 int i, j, k;
3480 CSG_Simple_Statistics *S;
3481 CSG_Matrix C;
3482
3483 C.Create(nVariables, nVariables);
3484
3485 //-----------------------------------------------------
3486 S = new CSG_Simple_Statistics[nVariables];
3487
3488 for(j=0; j<nVariables; j++)
3489 {
3490 for(i=0; i<nSamples; i++)
3491 {
3492 S[j] += Values[i][j];
3493 }
3494 }
3495
3496 //-----------------------------------------------------
3497 for(k=0; k<nVariables; k++)
3498 {
3499 for(j=k; j<nVariables; j++)
3500 {
3501 double cov = 0.;
3502
3503 for(i=0; i<nSamples; i++)
3504 {
3505 cov += (Values[i][j] - S[j].Get_Mean()) * (Values[i][k] - S[k].Get_Mean());
3506 }
3507
3508 cov /= nSamples;
3509
3510 if( !bCovariances )
3511 {
3512 cov /= (S[j].Get_StdDev() * S[k].Get_StdDev());
3513 }
3514
3515 C[j][k] = C[k][j] = cov;
3516 }
3517 }
3518
3519 //-----------------------------------------------------
3520 delete[](S);
3521
3522 return( C );
3523 }
3524
3525
3526 ///////////////////////////////////////////////////////////
3527 // //
3528 // //
3529 // //
3530 ///////////////////////////////////////////////////////////
3531
3532 //---------------------------------------------------------
3533