super_prac

#define MAX_DATA_POINTS 10000
#define MAX_CLUSTER 101
#define MAX_DATA_DIMENSION 5
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <omp.h>

int num_data_points;
int num_clusters;
int num_dimensions;
double low_high[MAX_DATA_DIMENSION][2];
double degree_of_memb[MAX_DATA_POINTS][MAX_CLUSTER];
double epsilon;
double fuzziness;
double data_point[MAX_DATA_POINTS][MAX_DATA_DIMENSION];
double cluster_centre[MAX_CLUSTER][MAX_DATA_DIMENSION];

int real_number_data_points;

int
init(char *fname) {
    int i, j, r, rval;
    FILE *f;
    double s;
    if ((f = fopen(fname, "r")) == NULL) {
        printf("Failed to open input file.");
        return -1;
    }
    fscanf(f, "%d %d %d", &num_data_points, &num_clusters, &num_dimensions);
    num_data_points=real_number_data_points;//добавлено
    if (num_clusters > MAX_CLUSTER) {
        printf("Number of clusters should be < %d\n", MAX_CLUSTER);
        goto failure;
    }
    if (num_data_points > MAX_DATA_POINTS) {
        printf("Number of data points should be < %d\n", MAX_DATA_POINTS);
        goto failure;
    }
    if (num_dimensions > MAX_DATA_DIMENSION) {
        printf("Number of dimensions should be >= 1.0 and < %d\n",
               MAX_DATA_DIMENSION);
        goto failure;
    }
    fscanf(f, "%lf", &fuzziness);
    if (fuzziness <= 1.0) {
        printf("Fuzzyness coefficient should be > 1.0\n");
        goto failure;
    }
    fscanf(f, "%lf", &epsilon);
    if (epsilon <= 0.0 || epsilon > 1.0) {
        printf("Termination criterion should be > 0.0 and <= 1.0\n");
        goto failure;
    }
    for (i = 0; i < num_data_points; i++) {
        for (j = 0; j < num_dimensions; j++) {
            fscanf(f, "%lf", &data_point[i][j]);
            if (data_point[i][j] < low_high[j][0])
                low_high[j][0] = data_point[i][j];
            if (data_point[i][j] > low_high[j][1])
                low_high[j][1] = data_point[i][j];
        }
    }
    #pragma omp parallel for shared(degree_of_memb) private(i, j, rval)
    for (i = 0; i < num_data_points; i++) {
        s = 0.0;
        r = 100;
        #pragma omp parallel for reduction(+:s), reduction(-:r)
        for (j = 1; j < num_clusters; j++) {
            rval = rand() % (r + 1);
            r -= rval;
            degree_of_memb[i][j] = rval / 100.0;
            s += degree_of_memb[i][j];
        }
        degree_of_memb[i][0] = 1.0 - s;
    }
    fclose(f);
    return 0;
failure:
    fclose(f);
    exit(1);
}

int
calculate_centre_vectors() {
    int i, j, k;
    double numerator, denominator;
    double t[MAX_DATA_POINTS][MAX_CLUSTER];
 //   i=0;
    #pragma omp parallel private(i, j)
    {
      #pragma omp for collapse (2)
    for (i = 0; i < num_data_points; i++) {
        for (j = 0; j < num_clusters; j++) {
            t[i][j] = pow(degree_of_memb[i][j], fuzziness);
        }
    }
    }
    #pragma omp parallel shared(t, data_point) private(k, j)
    {
    #pragma omp for collapse (2)
    for (j = 0; j < num_clusters; j++) {
        for (k = 0; k < num_dimensions; k++) {
            numerator = 0.0;
            denominator = 0.0;
          #pragma omp parallel for reduction(+:numerator, denominator)
            for (i = 0; i < num_data_points; i++) {
                numerator += t[i][j] * data_point[i][k];
                denominator += t[i][j];
            }
            cluster_centre[j][k] = numerator / denominator;
        }
    }
    }
    return 0;
}

double
get_norm(int i, int j) {
    int k;
    double sum = 0.0;
    #pragma omp parallel for reduction(+:sum) private(k)
    for (k = 0; k < num_dimensions; k++) {
        sum += pow(data_point[i][k] - cluster_centre[j][k], 2);
    }

    return sqrt(sum);
}

double
get_new_value(int i, int j) {
    int k;
    double t, p, sum;
    sum = 0.0;
    p = 2 / (fuzziness - 1);
   #pragma omp parallel for reduction(+:sum) private(k, t)
    for (k = 0; k < num_clusters; k++) {
        t = get_norm(i, j) / get_norm(i, k);
        t = pow(t, p);
        sum += t;
    }

    return 1.0 / sum;
}

double
update_degree_of_membership() {
    int i, j;
    double new_uij;
    double max_diff = 0.0, diff;

    #pragma omp parallel private(j, new_uij,diff)
    {
         #pragma omp for collapse (2)
    for (j = 0; j < num_clusters; j++) {
        for (i = 0; i < num_data_points; i++) {
            new_uij = get_new_value(i, j);
            diff = new_uij - degree_of_memb[i][j];
            if (diff > max_diff)
                max_diff = diff;
            degree_of_memb[i][j] = new_uij;
        }
    }
    }
    return max_diff;
}

int
fcm(char *fname) {
    double max_diff;
    int i=0;
    init(fname);

    do {
        i++;
        calculate_centre_vectors();
        max_diff = update_degree_of_membership();
    } while (max_diff > epsilon);
    return i;
}

int
main(int argc, char **argv) {
    int num_it, i;
    for (i=1; i<=60; i++)
    {
        omp_set_dynamic(0);
        omp_set_num_threads(i);
        for (real_number_data_points=1; real_number_data_points<MAX_DATA_POINTS; real_number_data_points++)
        {
            clock_t before = clock();
            num_it=fcm("file");
            clock_t after = clock();
            printf("%d ", 100*100*real_number_data_points*4*num_it/((after-before)/1000));
        }
        printf("\n");
    }

    return 0;
}