opti.c

#include "opti.h"
int LINE_POINTS = 1000000;
int CLOUD_POINTS = 10000;
double X_VALUE = 0.000001;
double RANDOMNESS = 0.01;
// Learning rate
double LR = 0.0001;
double ER = 0.00001;


// y(x)=ax+b
// crée plein de points sur une même ligne
Point* line(double a, double b){
    Point * points = (Point*)malloc(sizeof(Point)*LINE_POINTS);
    double x = 0.0;
    for(int i=0;i<LINE_POINTS;i++){
        Point d;
        d.x = x;
        d.y = (a*x)+b;
        points[i]=d;
        //printf("%lf, %lf", d.x, d.y);
        x+=X_VALUE;
    }
    return points;
}

// 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;
}

// Version 1 du nuage, se doit d'appeler line indirectement
Point* cloud1(Point* a_line){
    Point * my_cloud = malloc(sizeof(Point)*CLOUD_POINTS);
    for(int i=0;i<CLOUD_POINTS;i++){
        Point chosen_point = a_line[rand()%LINE_POINTS-1];
        double new_y = double_random(chosen_point.y-RANDOMNESS, chosen_point.y+RANDOMNESS);
        my_cloud[i].x=chosen_point.x;
        my_cloud[i].y=new_y;
        printf("\nx = %f | y = %f\n", my_cloud[i].x,my_cloud[i].y);
    }
    return my_cloud;
}

// Version 2 du nuage, peut être généré grâce aux a et b d'une droite
// Génère aléatoirement des points sur l'axe des x (qui ne sont donc pas espacés uniformément)
// le while ne sert pas à grand chose. cloud2_5() est censé rectifier ça. 
Point* cloud2(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 = double_random(0,1);
        for(int j=0;j<CLOUD_POINTS-1;j++){
            while(chosen_point.x==my_cloud[j].x){
                chosen_point.x = double_random(0,1);
            }
        }
        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;
}

double random_point_in_cloud2(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(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(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;
}


// 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;
}

// Méthode des moindres carrés
double* averages(Point* my_cloud, int cloud_size){
    //double* avg = (double*)malloc(sizeof(double)*4);
    double* avg = malloc(sizeof(double)*4);
    double x_avg, y_avg, x2_avg, xy_avg;
    for(int i=0; i<cloud_size-1;i++){
        x_avg+=my_cloud[i].x;
        y_avg+=my_cloud[i].y;
        x2_avg+=popow(my_cloud[i].x,2);
        xy_avg+=(my_cloud[i].x*my_cloud[i].y);
    }
    x_avg/=cloud_size;y_avg/=cloud_size;x2_avg/=cloud_size;xy_avg/=cloud_size;
    avg[0]=x_avg;avg[1]=y_avg;avg[2]=x2_avg;avg[3]=xy_avg;
    printf("x=%f - y=%f - x2=%f - xy=%f", avg[0],avg[1],avg[2],avg[3]);
    return avg;
}

// Méthode des moindres carrés
double* a_and_b(double* my_averages){
    double* a_n_b = (double*)malloc(sizeof(double)*2);
    double x_avg = my_averages[0];
    double y_avg = my_averages[1];
    double x2_avg = my_averages[2];
    double xy_avg = my_averages[3];
    double a = (xy_avg-(x_avg*y_avg))/(x2_avg-popow(x_avg,2));
    double b = y_avg-(a*x_avg);
    a_n_b[0]=a;
    a_n_b[1]=b;
    return a_n_b;
}

// Sommes (similaire à averages() mais ne divise pas par le nb de points)
double* sums(Point* my_cloud, int cloud_size){
    //double* avg = (double*)malloc(sizeof(double)*4);
    double* my_sums = (double*)malloc(sizeof(double)*4);
    double x, y, x2, xy;
    for(int i=0; i<cloud_size-1;i++){
        x+=my_cloud[i].x;
        y+=my_cloud[i].y;
        x2+=popow(my_cloud[i].x,2);
        xy+=(my_cloud[i].x*my_cloud[i].y);
    }
    my_sums[0]=x;my_sums[1]=y;my_sums[2]=x2;my_sums[3]=xy;
    printf("x=%f - y=%f - x2=%f - xy=%f", my_sums[0],my_sums[1],my_sums[2],my_sums[3]);
    return my_sums;
}

// La nouvelle fonction gradient 2 en 1, uniquement chez Lidl
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;
}

// La nouvelle fonction gradient_descent 2 en 1, uniquement chez Lidl
// Formule fausse et pas fait la norme donc bug. A refaire demain.
void gradient_descent_v2(double* a, double* b, Point* my_cloud, double current_cost){
    double gr_A,gr_B,my_new_a,my_new_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 : %f || Nombre d'itérations : %d\n ", current_cost, nb_it);
    //a_n_b[0]=a; a_n_b[1]=b;
    //return a_n_b;
}

double cost(double *a, double *b, Point* my_cloud, int cloud_size){
    double my_error;
    for(int i=0;i<cloud_size;i++){
        //E(a,b)= somme des (aXi+b+Yi)^2
        my_error+=popow((((*a)*my_cloud[i].x)+(*b)+my_cloud[i].y),2);
    }
    return my_error;
}

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

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);
    */
    //Point* first_line = line(7.0, 5.5);
    double a_init = 0.556;
    double b_init = 0.88732;
    Point* fluffy = cloud2_5(&a_init, &b_init);
    printf("\na initial = %f || b initial = %f\n", a_init, b_init);
    double a = double_random(-100,100);
    double b = double_random(-100,100);
    double initial_cost = cost(&a, &b, fluffy, CLOUD_POINTS);
    gradient_descent_v2(&a, &b, fluffy, initial_cost);
    printf("\na trouvé = %f || b trouvé = %f\n", a, b);
    printf("\ndone\n");
    free(fluffy);
    return 0;
}