-
Notifications
You must be signed in to change notification settings - Fork 16
/
Kmean_omp.c
173 lines (148 loc) · 6.19 KB
/
Kmean_omp.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
#include <stdio.h>
#include <stdlib.h>
#include <float.h> //for FLT_MAX
#include "../shared/timing.h" //for timer seconds()
#include "../shared/make_2D_array.h"
#include "../shared/ncdf_util.h"
#include "../shared/math_util.h"
/* This is the name of the data file we will read. */
//#define FILE_NAME "../test_data/Blobs_smp20000_fea30_cls8.nc"
#define FILE_NAME "../../Data_Analysis/data/SSWdata.nc"
#define TOL 0.0001
#define MAX_ITER 100
int main() {
/*
======================================================
---------------- Initialization ---------------------
======================================================
*/
int N_samples,N_features,N_clusters,N_repeat;
//i for samples; j for features; k for clusters (typically)
int i,j,k;
int k_best,initial_idx;
float** X;
int** GUESS;
float dist,dist_min,dist_sum_old,dist_sum_new,inert_best=FLT_MAX;
// get input data and its size
double iStart1 = seconds();
readX(FILE_NAME,&X,&GUESS,&N_samples,&N_features,&N_clusters,&N_repeat);
double iElaps1 = seconds() - iStart1;
// each data point belongs to which cluster
// values range from 0 to N_cluster-1
int* labels = (int *)malloc(N_samples*sizeof(int));
int* labels_best = (int *)malloc(N_samples*sizeof(int));
// The position of each cluster center.
// Two arrays are needed as we are calculating the distance to the
// old centers and accumulating the new centers in the same iteration.
float** old_cluster_centers = Make2DFloatArray(N_clusters,N_features);
float** new_cluster_centers = Make2DFloatArray(N_clusters,N_features);
// how many data points in the cluster
// needed by calculating the average position of data points in each cluster
int* cluster_sizes = (int *)malloc(N_clusters*sizeof(int));
/*
======================================================
---------------- Kmean stepping ---------------------
======================================================
*/
printf("=====Applying K-mean======\n");
// record timing results
double iStart2,iElaps2;
double iStart3a,iStart3b,iStart3c;
double iElaps3a=0,iElaps3b=0,iElaps3c=0;
/* Run the K-mean algorithm for N_repeat times with
* different starting points
*/
iStart2 = seconds();
for (int i_repeat=0; i_repeat < N_repeat; i_repeat++){
// guess initial centers
for (k=0; k<N_clusters; k++){
cluster_sizes[k] = 0; // for accumulating
// the index of data points as the initial guess for cluster centers
initial_idx = GUESS[i_repeat][k];
for (j=0; j<N_features; j++){
old_cluster_centers[k][j]=X[initial_idx][j];
//set the "new" array to 0 for accumulating
new_cluster_centers[k][j] = 0.0;
}
}
// core K-mean stepping (Expectation-Maximization) begins here!!
int i_iter = 0;//record iteration counts
dist_sum_new = 0.0;//prevent the firt iteration error
do {
i_iter++;
dist_sum_old = dist_sum_new;
dist_sum_new = 0.0;
// E-Step: assign points to the nearest cluster center
iStart3a = seconds();
#pragma omp parallel default(shared) \
private(i,j,k,k_best,dist,dist_min) \
reduction(+:dist_sum_new)
{
#pragma omp for schedule(static)
for (i = 0; i < N_samples; i++) {
k_best = 0;//assume cluster no.0 is the nearest
//dist_min = distance(N_features, X[i], old_cluster_centers[k_best]);
dist_min = correlation(N_features, X[i], old_cluster_centers[k_best]);
for (k = 1; k < N_clusters; k++){
//dist = distance(N_features, X[i], old_cluster_centers[k]);
dist = correlation(N_features, X[i], old_cluster_centers[k]);
if (dist < dist_min){
dist_min = dist;
k_best = k;
}
}
labels[i] = k_best;
dist_sum_new += dist_min;
} // end of omp for loop
} // end of omp parallel region
iElaps3a += (seconds()-iStart3a);
// M-Step first half: set the cluster centers to the mean
iStart3b = seconds();
for (i = 0; i < N_samples; i++) {
k_best = labels[i];
cluster_sizes[k_best]++; // add one more points to this cluster
// As the total number of samples in each cluster is not known yet,
// here we are just calculating the sum, not the mean.
for (j=0; j<N_features; j++)
new_cluster_centers[k_best][j] += X[i][j];
} // end of M-Step first half
iElaps3b += (seconds()-iStart3b);
iStart3c = seconds();
// M-Step second half: convert the sum to the mean
for (k=0; k<N_clusters; k++) {
for (j=0; j<N_features; j++) {
if (cluster_sizes[k] > 0) //avoid divide-by-zero error
// sum -> mean
old_cluster_centers[k][j] = new_cluster_centers[k][j] / cluster_sizes[k];
new_cluster_centers[k][j] = 0.0;//for the next iteration
}
cluster_sizes[k] = 0;//for the next iteration
} // end of M-Step second half
iElaps3c += (seconds()-iStart3c);
} while( i_iter==1 || ((dist_sum_old - dist_sum_new > TOL)&&i_iter<MAX_ITER) );
//end of K-mean stepping
//printf("Final inertia: %f, iteration: %d \n",dist_sum_new,i_iter);
// record the best results
if (dist_sum_new < inert_best) {
inert_best = dist_sum_new;
for (i = 0; i < N_samples; i++)
labels_best[i] = labels[i];
}
} //end of one repeated run
iElaps2 = seconds() - iStart2;
/*
======================================================
---------------- Finalization ---------------------
======================================================
*/
// write data back to NetCDF file
writeY(FILE_NAME,labels_best, inert_best);
// print summary
printf("Best inertia: %f \n",inert_best);
printf("Kmean total time use (ms): %f \n", iElaps2*1000.0);
printf("E-step time use (ms): %f \n", iElaps3a*1000.0);
printf("M-step-1st-half time use (ms): %f \n", iElaps3b*1000.0);
printf("M-step-2nd-half time use (ms): %f \n", iElaps3c*1000.0);
printf("I/O time use (ms): %f \n", iElaps1*1000.0);
return 0;
}