72 return( Value * Value );
92 return( floor(0.5 + Value) );
95 double d = pow(10., Decimals);
98 if( fabs(v - floor(v)) > 0. )
100 return( floor(0.5 + v) / d );
109 if( Decimals <= 0 || Value == 0. )
111 return( (
int)(0.5 + Value) );
114 Decimals = (int)(-(ceil(log10(fabs(Value))) - Decimals));
118 double d = pow(10., Decimals);
121 ? -((
int)(0.5 - Value * d)) / d
122 : ((
int)(0.5 + Value * d)) / d
127 double d = pow(10., -Decimals);
130 ? -((
int)(0.5 - Value / d)) * d
131 : ((
int)(0.5 + Value / d)) * d
146 Number = abs(Number);
148 return( Number < 10 ? 1 : 1 + (
int)log10((
double)Number) );
163 if( Width > 0 && Precision >= 0 )
return(
CSG_String::Format(
"%*.*e", Width, Precision, Number) );
171 if( Width > 0 && Precision >= 0 )
return(
CSG_String::Format(
"%*.*f", Width, Precision, Number) );
189 if( *((
int *)a) < *((
int *)b) )
192 if( *((
int *)a) > *((
int *)b) )
201 if( *((
double *)a) < *((
double *)b) )
204 if( *((
double *)a) > *((
double *)b) )
213 return( strcmp((
const char *)a, (
const char *)b) );
227 ? (Deg + Min / 60. + Sec / 3600.)
228 : (Deg - Min / 60. - Sec / 3600.)
235 Sec = fmod(Value < 0. ? -Value : Value, 360.);
237 Deg = (int)Sec; Sec = 60. * (Sec - Deg);
238 Min = (int)Sec; Sec = 60. * (Sec - Min);
276 return( 1. * rand() / (
double)RAND_MAX );
284 return( min + (max - min) * rand() / (
double)RAND_MAX );
307 w = x1 * x1 + x2 * x2;
311 w = sqrt((-2. * log(w)) / w);
313 return( mean + stddev * x1 * w );
341 Create(Mean, StdDev, Count);
346 Create(Values, bHoldValues);
354 m_Values.
Create(bHoldValues ?
sizeof(
double) : 0, 0, TSG_Array_Growth::SG_ARRAY_GROWTH_1);
431 double Scale = nValues / (double)
m_nValues;
545 m_Sum += Weight * Value;
546 m_Sum2 += Weight * Value*Value;
726 return( (
sLong)Index );
750 return( (
sLong)Index );
879 if( Value == m_Value[i] )
906 if( Value == m_Value[i] )
934 if( Value.
Cmp(m_Value[i]) == 0 )
961 if( Value.
Cmp(m_Value[i]) == 0 )
1022 return( pRecord ? (
int)pRecord->
Get_Index() : -1);
1030 return( pRecord ? (
int)pRecord->
Get_Index() : -1);
1038 return( pRecord ? (
int)pRecord->
Get_Index() : -1);
1053 (pRecord = m_pTable->
Add_Record())->Set_Value(0, Value);
1063 return( (
int)(m_pTable->
Get_Count() - 1) );
1073 (pRecord = m_pTable->
Add_Record())->Set_Value(0, Value);
1083 return( (
int)(m_pTable->
Get_Count() - 1) );
1093 (pRecord = m_pTable->
Add_Record())->Set_Value(0, Value);
1103 return( (
int)(m_pTable->
Get_Count() - 1) );
1131 return( pRecord ? pRecord->
asInt(1) : 0 );
1139 return( pRecord ? pRecord->
asInt(0) : 0 );
1147 return( pRecord ? pRecord->
asDouble(0) : 0 );
1170 for(
int i=1; i<m_pTable->
Get_Count(); i++)
1172 if( Count < m_pTable->Get_Record_byIndex(i)->
asInt(1) )
1192 for(
int i=1; i<m_pTable->
Get_Count(); i++)
1233 Create(nClasses, Minimum, Maximum);
1241 Create(nClasses, Values, Minimum, Maximum, maxSamples);
1249 Create(nClasses, pTable, Field, Minimum, Maximum, maxSamples);
1257 Create(nClasses, pGrid, Minimum, Maximum, maxSamples);
1265 Create(nClasses, pGrids, Minimum, Maximum, maxSamples);
1293 void CSG_Histogram::_On_Construction(
void)
1297 m_Cumulative = NULL;
1304 bool CSG_Histogram::_Create(
size_t nClasses,
double Minimum,
double Maximum)
1306 if( nClasses > 0 && Minimum < Maximum )
1310 m_Elements = (
size_t *)
SG_Calloc(nClasses,
sizeof(
size_t));
1311 m_Cumulative = (
size_t *)
SG_Calloc(nClasses,
sizeof(
size_t));
1313 if( m_Elements && m_Cumulative )
1315 m_nClasses = nClasses;
1316 m_Minimum = Minimum;
1317 m_Maximum = Maximum;
1318 m_ClassWidth = (Maximum - Minimum) / (
double)m_nClasses;
1332 m_Statistics += Value;
1334 if( m_nClasses > 0 && m_Minimum <= Value && Value <= m_Maximum )
1336 size_t Class = (size_t)((Value - m_Minimum) / m_ClassWidth);
1338 if( Class >= m_nClasses )
1340 Class = m_nClasses - 1;
1343 m_Elements[Class]++;
1350 if( m_nClasses > 0 && Scale > 0. )
1354 for(
size_t i=0; i<m_nClasses; i++)
1356 m_Elements[i] = (size_t)(Scale * m_Elements[i]);
1368 if( m_nClasses > 0 )
1372 m_nMaximum = m_Cumulative[0] = m_Elements[0];
1374 for(
size_t i=1; i<m_nClasses; i++)
1376 m_Cumulative[i] = m_Cumulative[i - 1] + m_Elements[i];
1378 if( m_nMaximum < m_Elements[i] )
1380 m_nMaximum = m_Elements[i];
1391 bool CSG_Histogram::_Update(
sLong nElements)
1393 if( nElements > 0 && m_Statistics.
Get_Count() > 0 )
1395 double Scale = (double)nElements / (
double)m_Statistics.
Get_Count();
1399 for(
size_t i=1; i<m_nClasses; i++)
1401 m_Elements[i] = (size_t)(0.5 + Scale * m_Elements[i]);
1411 || m_nClasses != Histogram.m_nClasses
1412 || m_Minimum != Histogram.m_Minimum
1413 || m_Maximum != Histogram.m_Maximum )
1418 m_Statistics += Histogram.m_Statistics;
1420 for(
size_t i=0; i<m_nClasses; i++ )
1422 m_Elements[i] += Histogram.m_Elements[i];
1434 if( m_nClasses < 2 ) {
return( 0. ); }
1436 if( Quantile <= 0. ) {
return( m_Minimum ); }
1437 if( Quantile >= 1. ) {
return( m_Maximum ); }
1441 for(
size_t i=0, n0=0; i<m_nClasses; n0=m_Cumulative[i++])
1443 if( n < m_Cumulative[i] )
1445 if( m_Cumulative[i] >= n0 )
1450 double d = (n - n0) / (
double)(m_Cumulative[i] - n0);
1452 return(
Get_Break(i) + d * m_ClassWidth );
1454 else if( n == m_Cumulative[i] )
1460 return( m_Maximum );
1478 if( m_nClasses < 2 ) {
return( 0. ); }
1480 if( Value <= m_Minimum ) {
return( 0. ); }
1481 if( Value >= m_Maximum ) {
return( 1. ); }
1483 size_t Class = (size_t)(m_nClasses * (Value - m_Minimum) / (m_Maximum - m_Minimum));
1485 if( Class >= m_nClasses )
1494 return( dq * (Value - m_Minimum) / m_ClassWidth );
1500 return( q0 + dq * (Value -
Get_Break(Class)) / m_ClassWidth );
1520 if( !_Create(Histogram.m_nClasses, Histogram.m_Minimum, Histogram.m_Maximum) )
1525 m_Statistics = Histogram.m_Statistics;
1526 m_ClassWidth = Histogram.m_ClassWidth;
1527 m_nMaximum = Histogram.m_nMaximum ;
1529 for(
size_t i=0; i<m_nClasses; i++)
1531 m_Cumulative[i] = Histogram.m_Cumulative[i];
1532 m_Elements [i] = Histogram.m_Elements [i];
1541 return( _Create(nClasses, Minimum, Maximum) );
1547 if( Minimum >= Maximum )
1555 if( !_Create(nClasses, Minimum, Maximum) )
1561 if( maxSamples > 0 && maxSamples < (
size_t)Values.
Get_N() )
1563 double d = (double)Values.
Get_N() / (double)maxSamples;
1565 for(
double i=0; i<(double)Values.
Get_N(); i+=d)
1570 d = (double)m_Statistics.
Get_Count() / (double)maxSamples;
1572 return( _Update(d < 1. ? (
int)(d * (double)Values.
Get_N()) : Values.
Get_N()) );
1576 for(
int i=0; i<Values.
Get_N(); i++)
1587 if( !pTable || Field < 0 || Field >= pTable->
Get_Field_Count() || !_Create(nClasses,
1588 Minimum < Maximum ? Minimum : pTable->Get_Minimum(Field),
1589 Minimum < Maximum ? Maximum : pTable->Get_Maximum(Field)) )
1595 if( maxSamples > 0 && maxSamples < (
size_t)pTable->
Get_Count() )
1597 double Value, d = (double)pTable->
Get_Count() / (double)maxSamples;
1599 for(
double i=0; i<(double)pTable->
Get_Count(); i+=d)
1607 d = (double)m_Statistics.
Get_Count() / (double)maxSamples;
1617 if( pTable->
Get_Value(i, Field, Value) )
1629 if( !pGrid || !_Create(nClasses,
1630 Minimum < Maximum ? Minimum : pGrid->Get_Min(),
1631 Minimum < Maximum ? Maximum : pGrid->Get_Max()) )
1637 if( maxSamples > 0 && (
sLong)maxSamples < pGrid->Get_NCells() )
1639 double d = (double)pGrid->
Get_NCells() / (double)maxSamples;
1641 for(
double i=0; i<(double)pGrid->
Get_NCells(); i+=d)
1649 d = (double)m_Statistics.
Get_Count() / (double)maxSamples;
1669 if( !pGrids || !_Create(nClasses,
1670 Minimum < Maximum ? Minimum : pGrids->Get_Min(),
1671 Minimum < Maximum ? Maximum : pGrids->Get_Max()) )
1677 if( maxSamples > 0 && (
sLong)maxSamples < pGrids->Get_NCells() )
1679 double d = (double)pGrids->
Get_NCells() / (double)maxSamples;
1681 for(
double i=0; i<(double)pGrids->
Get_NCells(); i+=d)
1689 d = (double)m_Statistics.
Get_Count() / (double)maxSamples;
1737 Create(pTable, Field, nClasses, Histogram);
1743 Create(pGrid, nClasses, Histogram);
1749 Create(pGrids, nClasses, Histogram);
1755 Create(Values, nClasses, Histogram);
1766 bool bResult =
false;
1770 bResult = m_Histogram.
Create(Histogram, pTable, Field) && _Histogram(nClasses);
1772 else if( Field >= 0 && Field < pTable->Get_Field_Count() )
1778 if( pTable->
Get_Value(i, Field, Value) )
1784 bResult = m_Values.
Sort() && _Calculate(nClasses);
1795 bool bResult =
false;
1799 bResult = m_Histogram.
Create(Histogram, pGrid) && _Histogram(nClasses);
1811 bResult = m_Values.
Sort() && _Calculate(nClasses);
1822 bool bResult =
false;
1826 bResult = m_Histogram.
Create(Histogram, pGrids) && _Histogram(nClasses);
1838 bResult = m_Values.
Sort() && _Calculate(nClasses);
1849 bool bResult =
false;
1853 bResult = m_Histogram.
Create(Histogram, Values) && _Histogram(nClasses);
1857 bResult = m_Values.
Create(Values) && m_Values.
Sort() && _Calculate(nClasses);
1871 bool CSG_Natural_Breaks::_Histogram(
int nClasses)
1873 if( _Calculate(nClasses) )
1881 m_Breaks[i] = m_Histogram.
Get_Value(m_Breaks[i] * d);
1897 inline double CSG_Natural_Breaks::_Get_Value(
int i)
1904 return( m_Values[i] );
1908 bool CSG_Natural_Breaks::_Calculate(
int nClasses)
1917 CSG_Matrix mv(nClasses, nValues); mv.Assign(FLT_MAX);
1919 int **mc = (
int **)
SG_Malloc(nValues *
sizeof(
int *));
1921 mc[0] = (
int *)
SG_Calloc((
size_t)nClasses * nValues,
sizeof(int));
1923 for(
int i=0; i<nValues; i++)
1925 mc[i] = mc[0] + i * (size_t)nClasses;
1929 for(
int i=1; i<nValues; i++)
1931 double v = 0., s1 = 0., s2 = 0., w = 0.;
1933 for(
int m=0, n=i+1; m<=i; m++, n--)
1939 v = s2 - (s1 * s1) / w;
1943 for(
int j=1; j<nClasses; j++)
1945 if( mv[i][j] >= (v + mv[n - 1][j - 1]) )
1948 mv[i][j] = v + mv[n - 1][j - 1];
1961 for(
int i=0; i<nClasses; i++)
1966 int j = Class[(size_t)nClasses - 1] = nValues - 1;
1968 for(
int i=nClasses-1; i>0; i--)
1970 Class[(size_t)i - 1] = j = mc[j - 1][i];
1974 m_Breaks.
Create((
size_t)nClasses + 1);
1976 m_Breaks[0] = _Get_Value(0);
1978 for(
int i=1; i<nClasses; i++)
1980 m_Breaks[i] = _Get_Value(Class[i - 1]);
1983 m_Breaks[nClasses] = _Get_Value(nValues - 1);
2031 m_nFeatures = nFeatures;
2033 m_Features.
Create(m_nFeatures *
sizeof(
double), 0, TSG_Array_Growth::SG_ARRAY_GROWTH_3);
2044 return( m_nFeatures > 0 && m_Features.
Inc_Array() );
2050 if( iElement >= 0 && iElement <
Get_nElements() && iFeature >= 0 && iFeature < m_nFeatures )
2052 ((
double *)m_Features.
Get_Entry(iElement))[iFeature] = Value;
2079 m_nMembers.
Create(nClusters);
2080 m_Variance.
Create(nClusters);
2081 m_Centroid.
Create(m_nFeatures, nClusters);
2088 switch( Initialization )
2093 m_Clusters[iElement] = nClusters - 1;
2099 m_Clusters[iElement] = iElement % nClusters;
2104 if( 0 > m_Clusters[iElement] || m_Clusters[iElement] >= nClusters )
2106 m_Clusters[iElement] = iElement % nClusters;
2119 default: bResult = _Minimum_Distance(
true , nMaxIterations);
break;
2120 case 1: bResult = _Hill_Climbing (
true , nMaxIterations);
break;
2121 case 2: bResult = _Minimum_Distance(
true , nMaxIterations)
2122 && _Hill_Climbing (
false, nMaxIterations);
break;
2128 for(
int iCluster=0; iCluster<nClusters; iCluster++)
2130 m_Variance[iCluster] = m_nMembers[iCluster] <= 0 ? 0. : m_Variance[iCluster] / m_nMembers[iCluster];
2138 bool CSG_Cluster_Analysis::_Minimum_Distance(
bool bInitialize,
int nMaxIterations)
2140 int iElement, iCluster, nClusters = m_Variance.
Get_N();
2142 double SP_Last = -1.;
2154 m_nMembers[iCluster = m_Clusters[iElement]]++;
2156 double *Feature = (
double *)m_Features.
Get_Entry(iElement);
2158 for(
int iFeature=0; iFeature<m_nFeatures; iFeature++)
2160 m_Centroid[iCluster][iFeature] += Feature[iFeature];
2165 for(iCluster=0; iCluster<nClusters; iCluster++)
2167 double d = m_nMembers[iCluster] > 0 ? 1. / m_nMembers[iCluster] : 0.;
2169 for(
int iFeature=0; iFeature<m_nFeatures; iFeature++)
2171 m_Centroid[iCluster][iFeature] *= d;
2176 int nShifts = 0; m_SP = 0.;
2180 double *Feature = (
double *)m_Features.
Get_Entry(iElement);
2182 double minVariance = -1.;
2183 int minCluster = -1;
2185 for(iCluster=0; iCluster<nClusters; iCluster++)
2187 double Variance = 0.;
2189 for(
int iFeature=0; iFeature<m_nFeatures; iFeature++)
2191 Variance +=
SG_Get_Square(m_Centroid[iCluster][iFeature] - Feature[iFeature]);
2194 if( minVariance < 0. || Variance < minVariance )
2196 minVariance = Variance;
2197 minCluster = iCluster;
2201 if( m_Clusters[iElement] != minCluster )
2203 m_Clusters[iElement] = minCluster;
2208 m_SP += minVariance;
2209 m_Variance[minCluster] += minVariance;
2216 _TL(
"pass" ), m_Iteration,
2217 _TL(
"change"), m_Iteration < 2 ? m_SP : SP_Last - m_SP
2222 if( nShifts == 0 || (nMaxIterations > 0 && nMaxIterations <= m_Iteration) )
2232 bool CSG_Cluster_Analysis::_Hill_Climbing(
bool bInitialize,
int nMaxIterations)
2234 int iElement, iCluster, nClusters = m_Variance.
Get_N();
2243 m_nMembers[iCluster = m_Clusters[iElement]]++;
2245 double *Feature = (
double *)m_Features.
Get_Entry(iElement);
2247 for(
int iFeature=0; iFeature<m_nFeatures; iFeature++)
2249 double d = Feature[iFeature];
2251 m_Centroid[iCluster][iFeature] += d;
2252 m_Variance[iCluster] += d*d;
2257 for(iCluster=0; iCluster<nClusters; iCluster++)
2259 double v = 0., d = m_nMembers[iCluster] <= 0 ? 0. : 1. / (double)m_nMembers[iCluster];
2261 for(
int iFeature=0; iFeature<m_nFeatures; iFeature++)
2263 m_Centroid[iCluster][iFeature] *= d;
2267 m_Variance[iCluster] -= v * m_nMembers[iCluster];
2271 double SP_Last = -1.;
int noShift = 0;
2277 iCluster = m_Clusters[iElement];
2279 if( noShift++ <
Get_nElements() && m_nMembers[iCluster] > 1 )
2281 double *Feature = (
double *)m_Features.
Get_Entry(iElement);
2283 double Variance = 0.;
2285 for(
int iFeature=0; iFeature<m_nFeatures; iFeature++)
2287 Variance +=
SG_Get_Square(m_Centroid[iCluster][iFeature] - Feature[iFeature]);
2290 double V1 = Variance * m_nMembers[iCluster] / (m_nMembers[iCluster] - 1.);
2293 int kCluster = 0;
double VMin = -1.;
2295 for(
int jCluster=0; jCluster<nClusters; jCluster++)
2297 if( jCluster != iCluster )
2301 for(
int iFeature=0; iFeature<m_nFeatures; iFeature++)
2303 Variance +=
SG_Get_Square(m_Centroid[jCluster][iFeature] - Feature[iFeature]);
2306 double V2 = Variance * m_nMembers[jCluster] / (m_nMembers[jCluster] + 1.);
2308 if( VMin < 0. || V2 < VMin )
2311 kCluster = jCluster;
2317 if( VMin >= 0 && VMin < V1 )
2320 m_Variance[iCluster] -= V1;
2321 m_Variance[kCluster] += VMin;
2322 V1 = 1. / (m_nMembers[iCluster] - 1.);
2323 double V2 = 1. / (m_nMembers[kCluster] + 1.);
2325 for(
int iFeature=0; iFeature<m_nFeatures; iFeature++)
2327 double d = Feature[iFeature];
2329 m_Centroid[iCluster][iFeature] = (m_nMembers[iCluster] * m_Centroid[iCluster][iFeature] - d) * V1;
2330 m_Centroid[kCluster][iFeature] = (m_nMembers[kCluster] * m_Centroid[kCluster][iFeature] + d) * V2;
2333 m_Clusters[iElement] = kCluster;
2335 m_nMembers[iCluster]--;
2336 m_nMembers[kCluster]++;
2342 for(iCluster=0, m_SP=0.; iCluster<nClusters; iCluster++)
2344 m_SP += m_Variance[iCluster];
2350 _TL(
"pass" ), m_Iteration,
2351 _TL(
"change"), m_Iteration <= 1 ? m_SP : SP_Last - m_SP
2356 if( noShift >=
Get_nElements() || (nMaxIterations > 0 && nMaxIterations <= m_Iteration) )
2375 m_nFeatures = 0; m_pClasses = NULL; m_nClasses = 0;
2377 m_Threshold_Distance = 0.;
2378 m_Threshold_Angle = 0.;
2379 m_Threshold_Probability = 0.;
2380 m_Probability_Relative =
false;
2404 m_nFeatures = nFeatures;
2411 if( m_nClasses > 0 )
2413 for(
int i=0; i<m_nClasses; i++)
2415 delete(m_pClasses[i]);
2452 m_bWTA[Method] = bOn;
2472 int nFeatures = m_nFeatures;
Destroy(); m_nFeatures = nFeatures;
2482 if( !Data(
"classes") || !Data(
"features") || !Data[
"features"](
"count") || Data[
"features"][
"count"].Get_Content().asInt() != m_nFeatures || m_nFeatures == 0 )
2487 if( Data[
"features"](
"info") )
2497 if( Classes[i].Cmp_Name(
"class") && Classes[i].Get_Child(
"id") )
2501 CClass *pClass =
new CClass(Classes[i][
"id"].Get_Content());
2503 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; }
2504 if( !pClass->m_Mean.from_String(Classes[i][
"mean"].
Get_Content()) || pClass->m_Mean.Get_N () != m_nFeatures ) { bAdd =
false; }
2505 if( !pClass->m_Min .from_String(Classes[i][
"min" ].
Get_Content()) || pClass->m_Min .Get_N () != m_nFeatures ) { bAdd =
false; }
2506 if( !pClass->m_Max .from_String(Classes[i][
"max" ].
Get_Content()) || pClass->m_Max .Get_N () != m_nFeatures ) { bAdd =
false; }
2515 m_pClasses = (CClass **)
SG_Realloc(m_pClasses, ((
size_t)m_nClasses + 1) *
sizeof(CClass *));
2516 m_pClasses[m_nClasses++] = pClass;
2518 pClass->m_Cov_Det = pClass->m_Cov.Get_Determinant();
2519 pClass->m_Cov_Inv = pClass->m_Cov.Get_Inverse();
2526 return( m_nClasses > 0 );
2532 if( m_nFeatures < 1 || m_nClasses < 1 || File.
is_Empty() )
2539 Data.
Set_Name (
"supervised_classifier");
2544 Features.
Add_Child(
"count", m_nFeatures);
2546 if( Feature_Info && *Feature_Info )
2548 Features.
Add_Child(
"info", Feature_Info);
2555 for(
int i=0; i<m_nClasses; i++)
2559 CClass *pClass = m_pClasses[i];
2562 Class.
Add_Child(
"mean", pClass->m_Mean.to_String());
2563 Class.
Add_Child(
"min" , pClass->m_Min .to_String());
2564 Class.
Add_Child(
"max" , pClass->m_Max .to_String());
2565 Class.
Add_Child(
"cov" , pClass->m_Cov .to_String());
2568 return( Data.
Save(File) );
2581 if( m_nFeatures > 0 && m_nClasses > 0 )
2585 for(
int iClass=0; iClass<m_nClasses; iClass++)
2587 CClass *pClass = m_pClasses[iClass];
2589 s +=
"\n____\n" + pClass->m_ID +
"\nFeature\tMean\tMin\tMax\tStdDev";
2591 for(
int i=0; i<m_nFeatures; i++)
2615 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 )
2620 CClass *pClass, **pClasses = (CClass **)
SG_Realloc(m_pClasses, ((
size_t)m_nClasses + 1) *
sizeof(CClass *));
2624 m_pClasses = pClasses;
2626 m_pClasses[m_nClasses++] = pClass =
new CClass(Class_ID);
2628 pClass->m_ID = Class_ID;
2629 pClass->m_Mean = Mean;
2630 pClass->m_Min = Min;
2631 pClass->m_Max = Max;
2632 pClass->m_Cov = Cov;
2652 for(
int i=0; i<m_nClasses; i++)
2654 m_pClasses[i]->m_Samples.Destroy();
2663 if( m_nFeatures > 0 && m_nFeatures == Features.
Get_N() )
2669 CClass **pClasses = (CClass **)
SG_Realloc(m_pClasses, ((
size_t)m_nClasses + 1) *
sizeof(CClass *));
2673 m_pClasses = pClasses;
2675 m_pClasses[iClass = m_nClasses++] =
new CClass(Class_ID);
2681 return( m_pClasses[iClass]->m_Samples.Add_Row(Features) );
2691 if( m_nFeatures < 1 || m_nClasses < 1 )
2696 for(
int iClass=0; iClass<m_nClasses; iClass++)
2698 if( !m_pClasses[iClass]->
Train() )
2704 if( bClear_Samples )
2718 bool CSG_Classifier_Supervised::CClass::Train(
void)
2720 if( m_Samples.Get_NCols() < 1 || m_Samples.Get_NRows() < 1 )
2726 m_Mean.Create(m_Samples.Get_NCols());
2727 m_Min .Create(m_Samples.Get_NCols());
2728 m_Max .Create(m_Samples.Get_NCols());
2730 for(
int iFeature=0; iFeature<m_Samples.Get_NCols(); iFeature++)
2734 for(
int iSample=0; iSample<m_Samples.Get_NRows(); iSample++)
2736 s += m_Samples[iSample][iFeature];
2745 m_Cov.Create(m_Samples.Get_NCols(), m_Samples.Get_NCols());
2747 for(
int iFeature=0; iFeature<m_Samples.Get_NCols(); iFeature++)
2749 for(
int jFeature=iFeature; jFeature<m_Samples.Get_NCols(); jFeature++)
2753 for(
int iSample=0; iSample<m_Samples.Get_NRows(); iSample++)
2755 cov += (m_Samples[iSample][iFeature] - m_Mean[iFeature]) * (m_Samples[iSample][jFeature] - m_Mean[jFeature]);
2758 if( m_Samples.Get_NRows() > 1 )
2760 cov /= m_Samples.Get_NRows() - 1;
2763 m_Cov[iFeature][jFeature] = m_Cov[jFeature][iFeature] = cov;
2767 m_Cov_Inv = m_Cov.Get_Inverse ();
2768 m_Cov_Det = m_Cov.Get_Determinant();
2784 if( m_nFeatures > 0 )
2801 Class = -1; Quality = 0.;
2817 return( Class >= 0 );
2876 void CSG_Classifier_Supervised::_Get_Binary_Encoding(
const CSG_Vector &Features,
int &Class,
double &Quality)
2880 CClass *pClass = m_pClasses[iClass];
2888 d += (Features(iFeature) < Mean_Spectral) == (pClass->m_Mean[iFeature] < pClass->m_Mean_Spectral) ? 0 : 1;
2892 d += (Features[iFeature ] < Features[iFeature + 1]) == (pClass->m_Mean[iFeature ] < pClass->m_Mean[iFeature + 1]) ? 0 : 1;
2896 d += (Features[iFeature - 1] < Features[iFeature ]) == (pClass->m_Mean[iFeature - 1] < pClass->m_Mean[iFeature ]) ? 0 : 1;
2900 d += (Features[iFeature - 1] < Features[iFeature + 1]) == (pClass->m_Mean[iFeature - 1] < pClass->m_Mean[iFeature + 1]) ? 0 : 1;
2904 if( Class < 0 || Quality > d )
2906 Class = iClass; Quality = d;
2912 void CSG_Classifier_Supervised::_Get_Parallel_Epiped(
const CSG_Vector &Features,
int &Class,
double &Quality)
2916 CClass *pClass = m_pClasses[iClass];
2918 bool bMember =
true;
2922 bMember = pClass->m_Min[iFeature] <= Features[iFeature] && Features[iFeature] <= pClass->m_Max[iFeature];
2927 Class = iClass; Quality++;
2933 void CSG_Classifier_Supervised::_Get_Minimum_Distance(
const CSG_Vector &Features,
int &Class,
double &Quality)
2937 CClass *pClass = m_pClasses[iClass];
2939 double Distance = (Features - pClass->m_Mean).Get_Length();
2941 if( Class < 0 || Quality > Distance )
2943 Class = iClass; Quality = Distance;
2947 if( m_Threshold_Distance > 0. && Quality > m_Threshold_Distance )
2954 void CSG_Classifier_Supervised::_Get_Mahalanobis_Distance(
const CSG_Vector &Features,
int &Class,
double &Quality)
2958 CClass *pClass = m_pClasses[iClass];
2962 double Distance = D * (pClass->m_Cov_Inv * D);
2964 if( Class < 0 || Quality > Distance )
2966 Class = iClass; Quality = Distance;
2970 if( m_Threshold_Distance > 0. && Quality > m_Threshold_Distance )
2977 void CSG_Classifier_Supervised::_Get_Maximum_Likelihood(
const CSG_Vector &Features,
int &Class,
double &Quality)
2983 CClass *pClass = m_pClasses[iClass];
2987 double Distance = D * (pClass->m_Cov_Inv * D);
2989 double Probability = pow(2. *
M_PI, -0.5 * m_nFeatures) * pow(pClass->m_Cov_Det, -0.5) * exp(-0.5 * Distance);
2992 dSum += Probability;
2994 if( Class < 0 || Quality < Probability )
2996 Class = iClass; Quality = Probability;
3002 if( m_Probability_Relative )
3004 Quality = 100. * Quality / dSum;
3007 if( m_Threshold_Probability > 0. && Quality < m_Threshold_Probability )
3015 void CSG_Classifier_Supervised::_Get_Spectral_Angle_Mapping(
const CSG_Vector &Features,
int &Class,
double &Quality)
3019 CClass *pClass = m_pClasses[iClass];
3021 double Angle = Features.
Get_Angle(pClass->m_Mean);
3023 if( Class < 0 || Quality > Angle )
3025 Class = iClass; Quality = Angle;
3031 if( m_Threshold_Angle > 0. && Quality > m_Threshold_Angle )
3038 void CSG_Classifier_Supervised::_Get_Spectral_Divergence(
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)
3049 int iClass;
double iQuality;
3051 if( m_bWTA[iMethod] &&
Get_Class(Features, iClass, iQuality, iMethod) && ++Votes[iClass] > Quality )
3053 Class = iClass; Quality = Votes[iClass];
3076 if( !T || !df || df < 1. )
3089 if( p <= 0. || p >= 1. || df < 1 )
3095 double t = 0, p0, p1, diff = 1.;
3099 while( fabs(diff) > 0.0001 )
3101 t = Get_T_Inv(p1, df);
3102 diff = Get_T_P(t, df) - p0;
3106 return( bNegative ? -t : t );
3139 const double a1 = 0.0000053830, a2 = 0.0000488906, a3 = 0.0000380036;
3140 const double a4 = 0.0032776263, a5 = 0.0211410061, a6 = 0.0498673470;
3144 double p = (((((a1 * z + a2) * z + a3) * z + a4) * z + a5) * z + a6) * z + 1.;
3146 return( pow(p, -16) );
3152 const double a0 = 2.5066282, a1 = -18.6150006, a2 = 41.3911977, a3 = -25.4410605;
3153 const double b1 = -8.4735109, b2 = 23.0833674, b3 = -21.0622410, b4 = 3.1308291;
3154 const double c0 = -2.7871893, c1 = -2.2979648, c2 = 4.8501413, c3 = 2.3212128;
3155 const double d1 = 3.5438892, d2 = 1.6370678;
3159 double r = sqrt(-log(0.5 - p));
3161 return( (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + 1.) );
3167 return( p * (((a3 * r + a2) * r + a1) * r + a0) / ((((b4 * r + b3) * r + b2) * r + b1) * r + 1.) );
3172 double CSG_Test_Distribution::Get_T_P(
double T,
int df)
3174 return( df == 1 ? 1. - 2. * atan(fabs(T)) /
M_PI
3175 : df == 2 ? 1. - fabs(T) / sqrt(T*T + 2.)
3176 : df == 3 ? 1. - 2. * (atan(fabs(T) / sqrt(3.)) + fabs(T) * sqrt(3.) / (T*T + 3.)) /
M_PI
3177 : df == 4 ? 1. - fabs(T) * (1. + 2. / (T*T + 4.)) / sqrt(T*T + 4.)
3183 double CSG_Test_Distribution::Get_T_Z(
double T,
int df)
3187 double A9, B9, T9, Z8, P7, B7, z;
3194 : A9 * (((1. - T9 * 0.75) * T9 / 3. - 0.5) * T9 + 1.) * T9;
3195 P7 = ((0.4 * Z8 + 3.3) * Z8 + 24.) * Z8 + 85.5;
3196 B7 = 0.8 * pow(Z8, 2.) + 100. + B9;
3197 z = (1. + (-P7 / B7 + Z8 + 3.) / B9) * sqrt(Z8);
3203 double CSG_Test_Distribution::Get_T_Inv(
double p,
int df)
3209 return( cos(p *
M_PI / 2.) / sin(p *
M_PI / 2.) );
3214 return( sqrt(2. / (p * (2. - p)) - 2.) );
3217 double a = 1. / (df - 0.5);
3218 double b = 48. / (a*a);
3219 double c = ((20700. * a / b - 98.) * a - 16.) * a + 96.36;
3220 double d = ((94.5 / (b + c) - 3.) / b + 1.) * sqrt(a *
M_PI / 2.) * df;
3222 double y = pow(x, 2. / df);
3231 c = c + 0.3 * (df - 4.5) * (x + 0.6);
3234 c = (((0.05 * d * x - 5) * x - 7.) * x - 2.) * x + b + c;
3235 y = (((((0.4 * y + 6.3) * y + 36.) * y + 94.5) / c - y - 3.) / b + 1.) * x;
3249 y = ((1. / (((df + 6.) / (df * y) - 0.089 * d - 0.822) * (df + 2.) * 3.)
3250 + 0.5 / (df + 4.)) * y - 1.) * (df + 1.) / (df + 2.) + 1. / y;
3253 return( sqrt(df * y) );
3264 double F = ((
sLong)nSamples - (
sLong)nPredictors - 1) * (R2 / nPredictors) / (1. - R2);
3276 if( F >= 0.00001 && dfn > 0 && dfd > 0 )
3278 if( F * dfn >= dfd || F > 1. + 20. / dfn + 10. / sqrt((
double)dfn) )
3280 p = Get_Gamma(F, dfn, dfd);
3284 p = 1. - Get_Gamma(1. / F, dfd, dfn);
3288 if( p <= 0. || p >= 1. )
3290 p = F > 1. ? 0. : F < 1. ? 1. : 0.5;
3299 if( alpha < 0. || alpha > 1. || dfd < 0 || dfn < 0 )
3309 const int ITERMAX = 100;
3310 const double EPSILON = 0.0001;
3313 double lo, hi, mid, p;
3320 for(i=0; i<ITERMAX; i++)
3345 for(i=0; i<ITERMAX; i++)
3366 mid = (hi + lo) / 2.;
3368 for(i=0; i<ITERMAX && (hi-lo)>
EPSILON*mid; i++)
3370 mid = (hi + lo) / 2.;
3375 else if( p > alpha )
3385 double CSG_Test_Distribution::Get_Gamma(
double F,
double dfn,
double dfd)
3389 const double EXPMIN = -30.;
3390 const double SMALL = 0.00000000001;
3392 double x, c, er, s, n, t1, t;
3397 x = dfd / (dfd + dfn * F);
3398 c = Get_Log_Gamma(dfn + dfd) - Get_Log_Gamma(dfn) - Get_Log_Gamma(dfd + 1.) + dfd * log(x) + dfn * log(1. - x);
3414 while( t > er || t > t1 )
3418 t *= ((dfn + n) * x / (dfd + n));
3426 double CSG_Test_Distribution::Get_Log_Gamma(
double a)
3430 const int ARGMIN = 6;
3432 const double HL2PI = 0.91893853320467275;
3434 int n = (int)floor(ARGMIN - a + 0.0001);
3444 g = (1. - g * (1. / 30. - g * (1. / 105. - g * (1. / 140. - g / 99.)))) / (12. * a);
3445 g = g + ((a - 0.5) * log(a) - a + HL2PI);
3447 for(
int i=0; i<n; i++)
3466 int nVariables = Values.
Get_NX();
3467 int nSamples = Values.
Get_NY();
3474 C.Create(nVariables, nVariables);
3479 for(j=0; j<nVariables; j++)
3481 for(i=0; i<nSamples; i++)
3483 S[j] += Values[i][j];
3488 for(k=0; k<nVariables; k++)
3490 for(j=k; j<nVariables; j++)
3494 for(i=0; i<nSamples; i++)
3496 cov += (Values[i][j] - S[j].
Get_Mean()) * (Values[i][k] - S[k].Get_Mean());
3506 C[j][k] =
C[k][j] = cov;