opti.c

#include "opti.h"

int CLOUD_POINTS = 5000;
double X_VALUE = 0.000001;
double RANDOMNESS = 0.01;
// Learning rate
double LR = 0.0001;
// Error rate
double ER = 0.000001;

// La fonction pow() de math.h avait un problème aavec mon Makefile donc je l'ai refaite
// et puis après avoir résolu le problème (emplacement du -lm) afin de pouvoir utiliser
// sqrt() j'ai décidé de la laisser car le nom est mignon
double popow(double v, double w){
    double u=v;
    for(int i=1;i<w;i++){
        v*=u;
    }
    return v;
}

// Génère un double random...
double double_random(double min, double max) 
{
    double my_random;
    my_random = (double)rand()/RAND_MAX*(max-min)+min;
    //printf ( "%f\n", my_random);
    return my_random;
}

double random_point_in_cloud2_5(Point* my_cloud){
    double v = double_random(0,1);
    for(int i=0;i<CLOUD_POINTS;i++){
        if(v==my_cloud[i].x){
            v =random_point_in_cloud2_5(my_cloud);
        }
    }
    return v;
}

Point* cloud2_5(double *a, double *b){
    Point * my_cloud = (Point*)malloc(sizeof(Point)*CLOUD_POINTS);
    for(int i=0;i<CLOUD_POINTS;i++){
        Point chosen_point;
        chosen_point.x=random_point_in_cloud2_5(my_cloud);
        double rj = double_random(-RANDOMNESS, RANDOMNESS);
        chosen_point.y = (*a*chosen_point.x)+*b+rj;
        my_cloud[i].x=chosen_point.x;
        my_cloud[i].y=chosen_point.y;
        //printf("\nx = %f | y = %f\n", my_cloud[i].x,my_cloud[i].y);
    }
    return my_cloud;
}

Point** cloud_splitter(Point* my_cloud){
    Point ** tiny_clouds = malloc(sizeof(Point*)*3);
    Point * g1 = (Point*)malloc(sizeof(Point)*(CLOUD_POINTS/3));
    Point * g2 = (Point*)malloc(sizeof(Point)*((CLOUD_POINTS/3)));
    Point * g3 = (Point*)malloc(sizeof(Point)*(CLOUD_POINTS-((CLOUD_POINTS*2)/3)));
    // Pas besoin de shuffle les points ici car déjà fait dans cloud2_5()
    for(int i=0;i<CLOUD_POINTS;i++){
        if(i<(CLOUD_POINTS/3)){
            g1[i]=my_cloud[i];
        }if(i>=(CLOUD_POINTS/3) && i<((CLOUD_POINTS*2)/3)){
            g2[i-(CLOUD_POINTS/3)]=my_cloud[i];
        }else{
            g3[i-((CLOUD_POINTS*2)/3)]=my_cloud[i];
        }
    }
    tiny_clouds[0]=g1;
    tiny_clouds[1]=g2;
    tiny_clouds[2]=g3;
    //free(g1);free(g2);free(g3);
    printf("\nSplitter ok\n");
    return tiny_clouds;

}

// version 3 du nuage, plus simple, x espacés uniformément
Point* cloud3(double *a, double *b){
    Point * my_cloud = (Point*)malloc(sizeof(Point)*CLOUD_POINTS);
    for(int i=0;i<CLOUD_POINTS;i++){
        my_cloud[i].x=((double)i/(double)CLOUD_POINTS);
        double rj = double_random(-RANDOMNESS, RANDOMNESS);
        my_cloud[i].y = ((*a)*my_cloud[i].x)+(*b)+rj;
    }
    return my_cloud;
}

// Pour comparer les chiffres sur papier VS ceux de ce programme
Point* cloud_test(){
    Point * my_cloud = (Point*)malloc(sizeof(Point)*3);
    my_cloud[0].x = 0.51;
    my_cloud[0].y = 0.24;
    my_cloud[1].x = 0.98;
    my_cloud[1].y = 0.867;
    my_cloud[2].x = 0.1;
    my_cloud[2].y = 0.0042;
    return my_cloud;
}

double gradient(double *a, double *b, Point* my_cloud, int cloud_size, bool is_a){
    double my_gradient;
    for(int i = 0;i<cloud_size;i++){
        //Somme des Xi(aXi+b-Yi)
        if(is_a){
            my_gradient+=(my_cloud[i].x)*((*a*(my_cloud[i].x)+*b-(my_cloud[i].y)));
        //Somme des aXi+b-Yi
        }else{
            my_gradient+=(*a*(my_cloud[i].x))+*b-(my_cloud[i].y);
        }
    }
    return my_gradient;
}

void gradient_descent_v4(double *a, double *b, Point* my_cloud){
    double gr_A,gr_B,my_new_a,my_new_b,current_cost;
    double *zero;
    double z = 0.0;
    zero = &z;
    current_cost = cost2(zero,zero,a,b);
    int nb_it = 0;
    //double* a_n_b = (double*)malloc(sizeof(double)*2);
    while(current_cost>=ER){
        gr_A = gradient(a, b, my_cloud, CLOUD_POINTS, true);
        gr_B = gradient(a, b, my_cloud, CLOUD_POINTS, false);
        // Les 2 lignes ci-dessous pourraient être condensées, mais c'est pour rendre le code
        // plus explicite que je sépare volontairement les étapes.
        my_new_a = (*a) - (LR * gr_A);
        //printf("\nNew A : %f\n", my_new_a);
        my_new_b = (*b) - (LR * gr_B);
        //printf("New B : %f\n", my_new_b);
        current_cost = cost2(a,b,&my_new_a,&my_new_b);
        //printf("\nCost : %f\n", current_cost);
        *a = my_new_a; *b = my_new_b;
        //current_cost = cost(a, b, my_cloud, CLOUD_POINTS);
        nb_it++;
    }
    printf("\nCost : %g || Nombre d'itérations : %d\n ", current_cost, nb_it);
}

double cost2(double *a, double *b, double *new_a, double *new_b){
    return sqrt(popow((*new_a - *a),2)+popow((*new_b - *b), 2));
} 

Point* cloud_merger(Point* cloud_a, Point* cloud_b, bool g3_here){
    int cloud_size;
    if(!g3_here){
        cloud_size = 2*(CLOUD_POINTS)/3;
    }else{
        cloud_size = CLOUD_POINTS-(CLOUD_POINTS/3);
    }
    Point * merged_cloud = (Point*)malloc(sizeof(Point)*cloud_size);
    for(int i=0;i<cloud_size;i++){
        if(i<(CLOUD_POINTS/3)){
            merged_cloud[i]=cloud_a[i];
        }else{
            merged_cloud[i]=cloud_b[i-(CLOUD_POINTS/3)];
        }
    }
    return merged_cloud;
}

int main (void){
    srand(time(NULL));
    /*// Create a vector X = [0,1,2...99]
    double_vector_t *X = iota(100);
    // Create a vector Y = my_function(x)
    double_vector_t *Y = apply_function(X, my_function);

    // Export our vectors into files
    export_vector("./X.vec", X);
    export_vector("./Y.vec", Y);

    // Free our vectors
    destroy_vector(&Y);
    destroy_vector(&X);
    */
    //Le a et le b calculé sur mon papier;
    double a_init = 0.98902527076;
    double b_init = -0.1537833935;
    //Point* fluffy = cloud_test();
    printf("\na optimal = %f || b optimal = %f\n", a_init, b_init);
    double a = double_random(0,1);
    double b = double_random(0,1);
    double c = a;
    double d = b;
    
    Point* fluffy = cloud2_5(&a_init, &b_init);
    //printf("\na random initial = %f || b random initial = %f\n", a, b);
    Point** cloud_slayer = cloud_splitter(fluffy);
    /*for(int j=0;j<3;j++){
        for(int i=0; i<CLOUD_POINTS/3;i++){
            printf("\nG%d x = %g\nG%d y = %g\n", j+1,cloud_slayer[j][i].x , j+1, cloud_slayer[j][i].y);
        }
    }*/
    Point* merged = cloud_merger(cloud_slayer[0], cloud_slayer[1], false);
    for(int i=0;i<CLOUD_POINTS*2/3;i++){
        printf("\nMERGED : x=%g  y=%g\n", merged[i].x, merged[i].y);
    }
    //gradient_descent_v4(&a, &b, fluffy);
    gradient_descent_v4(&a, &b, merged);
    gradient_descent_v4(&c, &d, cloud_slayer[2]);
    printf("\na1 trouvé = %f || b1 trouvé = %f\n", a, b);
    printf("\na2 trouvé = %f || b2 trouvé = %f\n", c, d);
    //printf("\nx = %g\ny = %g\n",cloud_slayer[0][1].x , cloud_slayer[0][1].y);
    free(merged);
    free(cloud_slayer);
    free(fluffy);
    printf("\ndone\n");
    return 0;
}