PCA降维的C语言实现

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int max_row_order(int n, double A[n][n])
{
    double m = 0;
    int mro = 0;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (j == i)
            {
                continue;
            }
            
            if (fabs(A[i][j]) > m)
            {
                m = fabs(A[i][j]);
                mro = i;
            }  
        }                 
    }
    int p = mro;
    return p;
}

int max_column_order(int n, double A[n][n])
{
    double m = 0;
    int mco = 0;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (j == i)
            {
                continue;
            }
            
            if (fabs(A[i][j]) > m)
            {
                m = fabs(A[i][j]);
                mco = j;
            }  
        }                 
    }
    int q = mco;
    return q;
}

double t_slove(int n, double A[n][n], int p, int q)
{
    double s = (A[q][q] - A[p][p])/(2*A[p][q]);
    double t;
    if (s == 0)
    {
        t = 1;
    }
    else
    {
        double x1 = -s + sqrt(s*s+1);
        double x2 = -s - sqrt(s*s+1);
        if (fabs(x1) > fabs(x2))
        {
            t = x2;
        }
        else
        {
            t = x1;
        }
    }
    return t;
}

double calculate_nondiagonal(int n, double A[n][n])
{
    double sum = 0;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            if (j == i)
            {
                continue;
            }
            sum = sum + A[i][j] * A[i][j];
        }                 
    }
    return sum;
}

void print_matrix(int row, int col, double matrix[row][col])
{
    for(int i = 0;i < row;i ++)
    {
        for(int j = 0;j < col;j ++)
        {
            printf("%.12lf\t",matrix[i][j]);
        }    
        printf("\n");
    }
    printf("----------------------------------------------");
    printf("\n");
}

void print_vector(int n, double vector[n])
{
    for(int i = 0;i < n;i ++)
    {
        printf("%.12lf\n",vector[i]);   
    }
    printf("----------------------------------------------");
    printf("\n");
}

double calculate_mode(int row, double x[row][1])
{
    double mode2 = 0;
    for (int i = 0; i < row; i++)
    {
        mode2 = mode2 + x[i][0]*x[i][0];
    }
    double mode = sqrt(mode2);
    return mode;
}

double * GaussEliminationWithPartialPivoting1(int n, double A[n-1][n-1], double b[n-1])
{
    double *result;
    int size = n-1;
    result = (double *)malloc(size*sizeof(double));     
    for (int i = 0; i < n-1; i++)
    {
        result[i] = 0;
    }

    for (int i = 0; i < n-1; i++)
    {
        int k = i;
        for (int j = i+1; j < n-1; j++)
        {
            if (fabs(A[k][i]) < fabs(A[j][i]))
            {
                k = j;
            }
        }
        for (int m = 0; m < n-1; m++)
        {
            double c;
            c = A[i][m];
            A[i][m] = A[k][m];
            A[k][m] = c;
        }
        double d;
        d = b[i];
        b[i] = b[k];
        b[k] = d;

        for (int j = i+1; j < n-1; j++)
        {
            double v = -A[j][i]/A[i][i];
            for (int k = 0; k < n-1; k++)
            {
                A[j][k] = v*A[i][k] + A[j][k];
            }
            b[j] = v*b[i] + b[j];
        }  
    }
    
    for (int i = n-2; i >= 0; i--)
    {
        if (i == n-2)
        {
            result[i] = 1.0;
        }
        else
        {
            for (int j = n-2; j >= i+1; j--)
            {
                b[i] = b[i] - A[i][j]*result[j];
            }
            result[i] = b[i]/A[i][i];
        }    
    }
    printf("After Elimination \n");
    for(int i = 0;i < n-1;i ++)
    {
        for(int j = 0;j < n-1;j ++)
        {
            printf("%lf\t",A[i][j]);
        }    
        printf("\n");
    }
    printf("----------------------------------------------");
    printf("\n");
    return result;
    free(result);
}

double calculate_mean (int m, double vector[m])
{
    double sum = 0;
    for (int i = 0; i < m; i++)
    {
        sum = sum + vector[i];
    }
    double mean = sum/m;
    return mean;
}

int main()
{
    FILE *fp = NULL;
    fp = fopen("iris.txt", "r");
    int m = 150;
    int n = 4;
    double a[m], b[m], c[m], d[m], e[m];
    for (int i = 0; i < m; i++)
    {
        fscanf(fp, "%lf,%lf,%lf,%lf,%lf", &a[i], &b[i],&c[i],&d[i],&e[i]);
    }
    double mean_a = calculate_mean(m, a);
    double mean_b = calculate_mean(m, b);
    double mean_c = calculate_mean(m, c);
    double mean_d = calculate_mean(m, d);
    double X[n][m];
    for (int i = 0; i < n; i++)//decentralize X
    {
        for (int j = 0; j < m; j++)
        {
            if (i == 0)
            {
                X[i][j] = a[j] - mean_a;
            }
            else if (i == 1)
            {
                X[i][j] = b[j] - mean_b;
            }
            else if (i == 2)
            {
                X[i][j] = c[j] - mean_c;
            }
            else if (i == 3)
            {
                X[i][j] = d[j] - mean_d;
            }
            
        }
        
    }
  
    //calculate XXT/m -> A
    double A[n][n];
    for (int i = 0; i < n; i++)//A = 0
    {
        for (int j = 0; j < n; j++)
        {
            A[i][j] = 0;
        }   
    }
    for (int i = 0; i < n; i++)//A = XXT/m
    {
        for (int j = 0; j < n; j++)
        {
            for (int k = 0; k < m; k++)
            {
                A[i][j] = A[i][j] + X[i][k]*X[j][k]/m;
            }            
        }   
    }
    double Ao[n][n];
    for (int i = 0; i < n; i++)//Ao = A
    {
        for (int j = 0; j < n; j++)
        {            
            Ao[i][j] = A[i][j];                      
        }   
    }
    printf("A (XXT/m) = \n");
    print_matrix(n, n, A);

    printf("=======================================================================\n");
    //calculate eigenvalue of A
    double limit = pow(10, -6);
    int p, q;
    double t, cos, sin;
    double Apio, Aqio;

    printf("calculate eigenvalue of A\n");
    int num = 0;
    while (calculate_nondiagonal(n, A) > limit)
    {        
        printf("sum of square of nondiagonal elements = %.12lf\n", calculate_nondiagonal(n, A));
        p = max_row_order(n,A);
        q = max_column_order(n,A);
        t = t_slove(n,A, p, q);
        cos = 1/sqrt(1 + t*t);
        sin = t/sqrt(1 + t*t);
        printf("Try %d \n", num+1);
        printf("choose p =%d, q = %d \n", p+1, q+1);
        printf("t = %lf, cos = %lf, sin = %lf \n", t, cos, sin);
        
        A[p][p] = A[p][p] - t*A[p][q];
        A[q][q] = A[q][q] + t*A[p][q];
        A[p][q] = 0;
        A[q][p] = 0;
        for (int i = 0; i < n; i++)
        {            
            if (i != p && i!= q)
            {
                Apio = A[p][i];
                Aqio = A[q][i];
                A[i][p] = cos*Apio - sin*Aqio;
                A[p][i] = cos*Apio - sin*Aqio;
                A[i][q] = sin*Apio + cos*Aqio;
                A[q][i] = sin*Apio + cos*Aqio;
            }                       
        }
        print_matrix(n, n, A);
        num ++;
    }
    double eigvalue[n][1];
    for (int i = 0; i < n; i++)
    {
        eigvalue[i][0] = A[i][i];
    }
    //sort eigenvalue
    for(int i = 0; i < n-1; i++)
    {
        int m = i;
        for(int j = i+1; j < n; j++)
        {
            if(eigvalue[j][0] > eigvalue[m][0])
                m = j;
        }            
        if(m != i)
        {
            double temp = eigvalue[i][0];
            eigvalue[i][0] = eigvalue[m][0];
            eigvalue[m][0] = temp;
        }
    }
    printf("eigenvalue of A = \n");
    print_matrix(n, 1, eigvalue);
    
    printf("=======================================================================\n");
    printf("choose the two biggest eigenvalue, and calculate the corresponding unit eigenvectors:\n");   
    //calculte unit eigenvector of A
    double I[n][n];
    for (int i = 0; i < n; i++)//initialize I
    {
        for (int j = 0; j < n; j++)
        {            
            if (i == j)
            {
                I[i][j] = 1;
            }
            else
            {
                I[i][j] = 0;
            }             
        }   
    }
    double bo[n];   
    double U[n][2];
    for (int ia = 0; ia < 2; ia++)
    {
        printf("when λ = %.12lf\n", eigvalue[ia][0]);
        double B[n][n];
        for (int i = 0; i < n; i++)//B = A - λΙ
        {
            for (int j = 0; j < n; j++)
            {
                B[i][j] = Ao[i][j] - I[i][j]*eigvalue[ia][0];
            }   
        }
        printf("A - λΙ = \n");
        print_matrix(n, n, B);    

        for (int i = 0; i < n; i++)//b = 0
        {
            bo[i] = 0;
        }
        double *result = GaussEliminationWithPartialPivoting1(n+1, B, bo);
        double u[n][1];
        for (int i = 0; i < n; i++)//u = result
        {
            u[i][0] = result[i];
        } 
        double mode = calculate_mode(n, u);       
        for (int i = 0; i < n; i++)//u = result
        {
            u[i][0] = u[i][0]/mode;
        } 
        printf("corresponding eigenvector = \n");
        print_matrix(n, 1, u);  
        for (int i = 0; i < n; i++)//U
        {
            U[i][ia] = u[i][0];               
        }
    }
    printf("U, composed by eigenvectors of AAT, = \n");
    print_matrix(n, 2, U);     
    
    /*==================================================================================*/
    printf("=======================================================================\n");
    //calculate projection
    double proj1[m];
    for (int i = 0; i < m; i++)//projection = 0
    {
        proj1[i] = 0;        
    }
    for (int i = 0; i < m; i++)//projection on eigenvector1
    {
        for (int j = 0; j < n; j++)
        {
            proj1[i] = proj1[i] + U[j][0]*X[j][i];
        }        
    }
    double proj2[m];
    for (int i = 0; i < m; i++)//projection = 0
    {
        proj2[i] = 0;        
    }
    for (int i = 0; i < m; i++)//projection on eigenvector2
    {
        for (int j = 0; j < n; j++)
        {
            proj2[i] = proj2[i] + U[j][1]*X[j][i];
        }        
    }
    printf("projection1: \n");
    print_vector(m, proj1);
    printf("projection2: \n");
    print_vector(m, proj2);
    //
    printf("these data would be presented RED\n");
    printf("projection1\tprojection2\n");
    FILE *fr = NULL;
    fr = fopen("RED.txt", "w+");
    for (int i = 0; i < m; i++)
    {
        if (e[i] == 0)
        {
            printf("%lf\t%lf\n", proj1[i], proj2[i]);
            fprintf(fr, "%lf %lf\n", proj1[i], proj2[i]);
        }        
    }
    fclose(fr);
    //
    printf("these data would be presented GREEN\n");
    printf("projection1\tprojection2\n");
    FILE *fg = NULL;
    fg = fopen("GREEN.txt", "w+");
    for (int i = 1; i < m; i++)
    {
        if (e[i] == 1)
        {
            printf("%lf\t%lf\n", proj1[i], proj2[i]);
            fprintf(fg, "%lf %lf\n", proj1[i], proj2[i]);
        }        
    }
    fclose(fg);
    //
    printf("these data would be presented BLUE\n");
    printf("projection1\tprojection2\n");
    FILE *fb = NULL;
    fb = fopen("BLUE.txt", "w+");
    for (int i = 0; i < m; i++)
    {
        if (e[i] == 2)
        {
            printf("%lf\t%lf\n", proj1[i], proj2[i]);
            fprintf(fb, "%lf %lf\n", proj1[i], proj2[i]);
        }        
    }
    fclose(fb);
      
    fclose(fp);
    return 0;
}