1 /*
 2 
 3 App : MK.cpp
 4 
 5 This is a C++ program for
 6 detecting abrupt climatic change
 7 by using Mann-Kendall technique.
 8 
 9 BY : ZXD
10 Date : January 10, 2014
11 
12 */
13 
14 #define _CRT_SECURE_NO_WARNINGS
15 
16 #include <fstream>
17 #include <cmath>
18 
19 #define N_SIZE 51
20 // The N_SIZE is size of sample,
21 // e.i. the number of total units.
22 
23 void smk( int n, double *x, double *uf, double *ub );
24 
25 int main()
26 {
27     double x[N_SIZE] = { 0 };
28     double uf[N_SIZE] = { 0 };
29     double ub[N_SIZE] = { 0 };
30     int n = N_SIZE;
31     int tmp;
32 
33     FILE *infp;
34     FILE *outfp;
35 
36     infp = fopen( "data4.txt", "r" );
37     for ( int i = 0; i < n; i++ )
38     {
39         fscanf( infp, "%d %lf", &tmp, &x[i] );
40         //printf( "%10.4lf\n", x[i] );
41     }
42     fclose( infp );
43 
44     smk( n, x, uf, ub );
45 
46     outfp = fopen( "result.txt", "w" );
47     for ( int i = 0; i < n; i++ )
48     {
49         fprintf( outfp, "%10.4lf %10.4lf\n", uf[i], ub[i] );
50         printf( "%10.4lf %10.4lf\n", uf[i], ub[i] );
51     }
52 
53     fclose( outfp );
54 
55     system( "pause" );
56 }
57 
58 void smk( int n, double *x, double *uf, double *ub )
59 {
60     double alpha[N_SIZE][N_SIZE] = { 0 };
61     double s[N_SIZE] = { 0 };
62     double k = 0.0;
63 
64     for ( int i = 0; i < n; i++ )
65     for ( int j = 0; j < n; j++ )
66         alpha[i][j] = ( x[i] > x[j] ? 1 : 0 );
67 
68     k = 2.0;
69     s[0] = 0.0;
70     for ( int i = 1; i < n; i++ )
71     {
72         s[i] = s[i - 1];
73         for ( int j = 0; j < i; j++ )
74             s[i] += alpha[i][j];
75         uf[i] = ( s[i] - k * ( k - 1.0 ) / 4.0 ) / sqrt( k * ( k - 1.0 ) * ( k * 2.0 + 5.0 ) / 72.0 );
76         k += 1.0;
77     }
78 
79     k = 2.0;
80     s[n - 1] = 0.0;
81     for ( int i = n - 2; i >= 0; i-- )
82     {
83         s[i] = s[i + 1];
84         for ( int j = n - 1; j > i; j-- )
85             s[i] += alpha[i][j];
86         ub[i] = ( k * ( k - 1.0 ) / 4.0 - s[i] ) / sqrt( k * ( k - 1.0 ) * ( k * 2.0 + 5.0 ) / 72.0 );
87         k += 1.0;
88     }
89 }