Loading...
Searching...
No Matches
class.c
Go to the documentation of this file.
/* functions to classify sorted arrays of doubles and fill a vector of
* classbreaks */
#include <stdbool.h>
#include <grass/glocale.h>
#include <grass/arraystats.h>
int AS_option_to_algorithm(const struct Option *option)
{
if (G_strcasecmp(option->answer, "int") == 0)
return CLASS_INTERVAL;
if (G_strcasecmp(option->answer, "std") == 0)
return CLASS_STDEV;
if (G_strcasecmp(option->answer, "qua") == 0)
return CLASS_QUANT;
if (G_strcasecmp(option->answer, "equ") == 0)
return CLASS_EQUIPROB;
if (G_strcasecmp(option->answer, "dis") == 0)
return CLASS_DISCONT;
G_fatal_error(_("Unknown algorithm '%s'"), option->answer);
}
double AS_class_apply_algorithm(int algo, const double data[], int nrec,
int *nbreaks, double classbreaks[])
{
double finfo = 0.0;
switch (algo) {
case CLASS_INTERVAL:
finfo = AS_class_interval(data, nrec, *nbreaks, classbreaks);
break;
case CLASS_STDEV:
finfo = AS_class_stdev(data, nrec, *nbreaks, classbreaks);
break;
case CLASS_QUANT:
finfo = AS_class_quant(data, nrec, *nbreaks, classbreaks);
break;
case CLASS_EQUIPROB:
finfo = AS_class_equiprob(data, nrec, nbreaks, classbreaks);
break;
case CLASS_DISCONT:
finfo = AS_class_discont(data, nrec, *nbreaks, classbreaks);
break;
default:
break;
}
if (finfo == 0)
G_fatal_error(_("Classification algorithm failed"));
return finfo;
}
int AS_class_interval(const double data[], int count, int nbreaks,
double classbreaks[])
{
double min, max;
double step;
int i = 0;
min = data[0];
max = data[count - 1];
step = (max - min) / (nbreaks + 1);
for (i = 0; i < nbreaks; i++)
classbreaks[i] = min + (step * (i + 1));
return (1);
}
double AS_class_stdev(const double data[], int count, int nbreaks,
double classbreaks[])
{
struct GASTATS stats;
int i;
int nbclass;
double scale = 1.0;
AS_basic_stats(data, count, &stats);
nbclass = nbreaks + 1;
if (nbclass % 2 ==
1) { /* number of classes is uneven so we center middle class on mean */
/* find appropriate fraction of stdev for step */
i = 1;
while (i) {
if (((stats.mean + stats.stdev * scale / 2) +
(stats.stdev * scale * (nbclass / 2 - 1)) >
stats.max) ||
((stats.mean - stats.stdev * scale / 2) -
(stats.stdev * scale * (nbclass / 2 - 1)) <
stats.min))
scale = scale / 2;
else
i = 0;
}
/* classbreaks below the mean */
for (i = 0; i < nbreaks / 2; i++)
classbreaks[i] = (stats.mean - stats.stdev * scale / 2) -
stats.stdev * scale * (nbreaks / 2 - (i + 1));
/* classbreaks above the mean */
for (; i < nbreaks; i++)
classbreaks[i] = (stats.mean + stats.stdev * scale / 2) +
stats.stdev * scale * (i - nbreaks / 2);
}
else { /* number of classes is even so mean is a classbreak */
/* decide whether to use 1*stdev or 0.5*stdev as step */
i = 1;
while (i) {
if (((stats.mean) + (stats.stdev * scale * (nbclass / 2 - 1)) >
stats.max) ||
((stats.mean) - (stats.stdev * scale * (nbclass / 2 - 1)) <
stats.min))
scale = scale / 2;
else
i = 0;
}
/* classbreaks below the mean and on the mean */
for (i = 0; i <= nbreaks / 2; i++)
classbreaks[i] =
stats.mean - stats.stdev * scale * (nbreaks / 2 - i);
/* classbreaks above the mean */
for (; i < nbreaks; i++)
classbreaks[i] =
stats.mean + stats.stdev * scale * (i - nbreaks / 2);
}
return (scale);
}
int AS_class_quant(const double data[], int count, int nbreaks,
double classbreaks[])
{
int i, step;
step = count / (nbreaks + 1);
for (i = 0; i < nbreaks; i++)
classbreaks[i] = data[step * (i + 1)];
return (1);
}
int AS_class_equiprob(const double data[], int count, int *nbreaks,
double classbreaks[])
{
int i, j;
double *lequi; /*Vector of scale factors for probabilities of the normal
distribution */
struct GASTATS stats;
int nbclass;
nbclass = *nbreaks + 1;
lequi = G_malloc(*nbreaks * sizeof(double));
/* The following values come from the normal distribution and will be used
* as: classbreak[i] = (lequi[i] * stdev) + mean;
*/
if (nbclass < 3) {
lequi[0] = 0;
}
else if (nbclass == 3) {
lequi[0] = -0.43076;
lequi[1] = 0.43076;
}
else if (nbclass == 4) {
lequi[0] = -0.6745;
lequi[1] = 0;
lequi[2] = 0.6745;
}
else if (nbclass == 5) {
lequi[0] = -0.8416;
lequi[1] = -0.2533;
lequi[2] = 0.2533;
lequi[3] = 0.8416;
}
else if (nbclass == 6) {
lequi[0] = -0.9676;
lequi[1] = -0.43076;
lequi[2] = 0;
lequi[3] = 0.43076;
lequi[4] = 0.9676;
}
else if (nbclass == 7) {
lequi[0] = -1.068;
lequi[1] = -0.566;
lequi[2] = -0.18;
lequi[3] = 0.18;
lequi[4] = 0.566;
lequi[5] = 1.068;
}
else if (nbclass == 8) {
lequi[0] = -1.1507;
lequi[1] = -0.6745;
lequi[2] = -0.3187;
lequi[3] = 0;
lequi[4] = 0.3187;
lequi[5] = 0.6745;
lequi[6] = 1.1507;
}
else if (nbclass == 9) {
lequi[0] = -1.2208;
lequi[1] = -0.7648;
lequi[2] = -0.4385;
lequi[3] = -0.1397;
lequi[4] = 0.1397;
lequi[5] = 0.4385;
lequi[6] = 0.7648;
lequi[7] = 1.2208;
}
else if (nbclass == 10) {
lequi[0] = -1.28155;
lequi[1] = -0.84162;
lequi[2] = -0.5244;
lequi[3] = -0.25335;
lequi[4] = 0;
lequi[5] = 0.25335;
lequi[6] = 0.5244;
lequi[7] = 0.84162;
lequi[8] = 1.28155;
}
else {
G_fatal_error(
_("Equiprobable classbreaks currently limited to 10 classes"));
}
AS_basic_stats(data, count, &stats);
/* Check if any of the classbreaks would fall outside of the range min-max
*/
j = 0;
for (i = 0; i < *nbreaks; i++) {
if ((lequi[i] * stats.stdev + stats.mean) >= stats.min &&
(lequi[i] * stats.stdev) + stats.mean <= stats.max) {
j++;
}
}
if (j < (*nbreaks)) {
G_warning(
_("There are classbreaks outside the range min-max. Number of "
"classes reduced to %i, but using probabilities for %i classes."),
j + 1, *nbreaks + 1);
for (i = 0; i < j; i++)
classbreaks[i] = 0.0;
}
j = 0;
for (i = 0; i < *nbreaks; i++) {
if ((lequi[i] * stats.stdev + stats.mean) >= stats.min &&
(lequi[i] * stats.stdev) + stats.mean <= stats.max) {
classbreaks[j] = lequi[i] * stats.stdev + stats.mean;
j++;
}
}
*nbreaks = j;
G_free(lequi);
return (1);
}
double AS_class_discont(const double data[], int count, int nbreaks,
double classbreaks[])
{
int i, j;
double chi2 = 1000.0;
/* get the number of values */
int n = count;
int nbclass = nbreaks + 1;
int *num = G_calloc((nbclass + 2), sizeof(int));
int *no = G_calloc((nbclass + 1), sizeof(int));
double *zz = G_malloc((nbclass + 1) * sizeof(double));
double *xn = G_malloc((n + 1) * sizeof(double));
double *co = G_malloc((nbclass + 1) * sizeof(double));
/* We copy the array of values to x, in order to be able
to standardize it */
double *x = G_malloc((n + 1) * sizeof(double));
x[0] = 0.0;
xn[0] = 0.0;
double min = data[0];
double max = data[count - 1];
for (i = 1; i <= n; i++)
x[i] = data[i - 1];
double rangemax = max - min;
double rangemin = rangemax;
for (i = 2; i <= n; i++) {
if (x[i] != x[i - 1] && x[i] - x[i - 1] < rangemin)
rangemin = x[i] - x[i - 1]; /* rangemin = minimal distance */
}
/* STANDARDIZATION
* and creation of the number vector (xn) */
for (i = 1; i <= n; i++) {
x[i] = (x[i] - min) / rangemax;
xn[i] = i / (double)n;
}
double xlim = rangemin / rangemax;
rangemin = rangemin / 2.0;
/* Searching for the limits */
num[1] = n;
/* Loop through possible solutions */
for (i = 1; i <= nbclass; i++) {
double dmax = 0.0;
int nmax = 0;
int nf = 0; /* End number */
/* Loop through classes */
for (j = 1; j <= i; j++) {
double a = 0.0, b = 0.0, c = 0.0, d = 0.0;
int nd = nf; /* Start number */
nf = num[j];
co[j] = 10e37;
AS_eqdrt(x, xn, nd, nf, &a, &b, &c);
double den = sqrt(pow(b, 2) + 1.0);
nd++;
/* Loop through observations */
for (int k = nd; k <= nf; k++) {
if (fabs(c) >= GRASS_EPSILON)
d = fabs(x[k] - c);
else
d = fabs((-1.0 * b * x[k]) + xn[k] - a) / den;
if (x[k] - x[nd] < xlim)
continue;
if (x[nf] - x[k] < xlim)
continue;
if (d <= dmax)
continue;
dmax = d;
nmax = k;
}
nd--;
if (fabs(x[nf] - x[nd]) > GRASS_EPSILON)
co[j] = (xn[nf] - xn[nd]) / (x[nf] - x[nd]);
}
for (j = 1; j <= i; j++) {
no[j] = num[j];
zz[j] = x[num[j]] * rangemax + min;
if (j == i)
continue;
if (co[j] > co[j + 1]) {
zz[j] = zz[j] + rangemin;
continue;
}
else {
zz[j] = zz[j] - rangemin;
no[j] = no[j] - 1;
}
}
int im = i - 1;
if (im != 0) {
for (j = 1; j <= im; j++) {
int ji = i + 1 - j;
no[ji] -= no[ji - 1];
}
}
if (nmax == 0) {
break;
}
int jj = 0;
int nff = i + 2;
bool do_reset = true;
for (j = 1; j <= i; j++) {
jj = nff - j;
if (num[jj - 1] < nmax) {
num[jj] = nmax;
do_reset = false;
break;
}
num[jj] = num[jj - 1];
}
if (do_reset) {
num[1] = nmax;
jj = 1;
}
int no1 = (int)((xn[num[jj]] - xn[num[jj - 1]]) * n);
int no2 = (int)((xn[num[jj + 1]] - xn[num[jj]]) * n);
double f = (xn[num[jj + 1]] - xn[num[jj - 1]]) /
(x[num[jj + 1]] - x[num[jj - 1]]);
f *= n;
double xt1 = (x[num[jj]] - x[num[jj - 1]]) * f;
double xt2 = (x[num[jj + 1]] - x[num[jj]]) * f;
if (fabs(xt1 * xt2) <= GRASS_EPSILON) {
if (fabs(xt2) > GRASS_EPSILON) {
xt2 = rangemin / 2.0 / rangemax * f;
xt1 = xt1 - xt2;
}
else {
xt1 = rangemin / 2.0 / rangemax * f;
xt2 = xt2 - xt1;
}
}
/* calculate chi-square to indicate statistical significance of new
* class, i.e. how probable would it be that the new class could be the
* result of purely random choice */
double ch = pow((double)((no1 - no2) - (xt1 - xt2)), 2) / (xt1 + xt2);
if (chi2 > ch)
chi2 = ch;
}
/* Fill up classbreaks of i <nbclass classes */
for (j = 1; j < nbclass; j++)
classbreaks[j - 1] = zz[j];
G_free(co);
G_free(no);
G_free(num);
G_free(x);
G_free(xn);
G_free(zz);
return (chi2);
}
int AS_class_frequencies(const double data[], int count, int nbreaks,
double classbreaks[], int frequencies[])
{
int i, j;
/* min = data[0];
max = data[count - 1]; */
/* count cases in all classes, except for last class */
i = 0;
for (j = 0; j < nbreaks; j++) {
while (data[i] <= classbreaks[j]) {
frequencies[j]++;
i++;
}
}
/*Now count cases in last class */
for (; i < count; i++) {
frequencies[nbreaks]++;
}
return (1);
}
G_free
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:147
AS_eqdrt
void AS_eqdrt(double vectx[], double vecty[], int i1, int i2, double *a, double *b, double *c)
Definition basic.c:39
AS_basic_stats
void AS_basic_stats(const double data[], int count, struct GASTATS *stats)
Definition basic.c:7
AS_class_discont
double AS_class_discont(const double data[], int count, int nbreaks, double classbreaks[])
Definition class.c:272
AS_class_apply_algorithm
double AS_class_apply_algorithm(int algo, const double data[], int nrec, int *nbreaks, double classbreaks[])
Definition class.c:25
AS_class_quant
int AS_class_quant(const double data[], int count, int nbreaks, double classbreaks[])
Definition class.c:139
AS_class_frequencies
int AS_class_frequencies(const double data[], int count, int nbreaks, double classbreaks[], int frequencies[])
Definition class.c:437
AS_class_equiprob
int AS_class_equiprob(const double data[], int count, int *nbreaks, double classbreaks[])
Definition class.c:152
AS_class_interval
int AS_class_interval(const double data[], int count, int nbreaks, double classbreaks[])
Definition class.c:56
AS_option_to_algorithm
int AS_option_to_algorithm(const struct Option *option)
Definition class.c:9
AS_class_stdev
double AS_class_stdev(const double data[], int count, int nbreaks, double classbreaks[])
Definition class.c:74
b
double b
Definition driver/set_window.c:5
G_fatal_error
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
count
int count
min
#define min(a, b)
Definition pngdriver/draw_bitmap.c:19
max
#define max(a, b)
Definition pngdriver/draw_bitmap.c:22
G_strcasecmp
int G_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower).
Definition strings.c:47
x
#define x