Code fini

.gitignore

+*.vec
+.vscode
+__pycache__
+example
+*.o
\ No newline at end of file

C/main.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "vector.h"
+#include "opti.h"
+
+double N;
+double len;
+couple a;
+
+double line_a = 1.1;
+double line_b = 0.1;
+double line_rand = 0.25; // Must be positif
+
+double genX(double x){
+    return my_random(0, 1.0);
+}
+
+/** @brief 
+ *  An exemple of mathematical function 
+ *  in 1 dimension.
+ * 
+ *  @param x A double variable.
+ *  @return f(x) = 2.0*sin(x) - 3.0*cos(x)
+ */
+double my_function(double x)
+{
+    double r = my_random(-line_rand, line_rand);
+    return line_a*x+line_b+r;
+}
+
+double genLine(double x){
+    return a.a*x+a.b;
+}
+
+double genXline(double x){
+    return x/len;
+}
+
+int main()
+{   
+    N = 100.0;
+    len = N*3;
+    
+    /* Intializes my_random number generator */
+    time_t t;   
+    srand((unsigned) time(&t));
+    
+    // Create a vector X = [0,1,2...99]
+    double_vector_t *X1 = iota((uint32_t)len);
+    double_vector_t *X2 = apply_function(X1, genXline);
+    double_vector_t *X  = apply_function(X1, genX);
+    
+    double_vector_t *G1_X = arrayCopy(X, 0, N);
+    double_vector_t *G2_X = arrayCopy(X, N, N);
+    double_vector_t *G3_X = arrayCopy(X, N*2, N);
+
+    double_vector_t *G1_G2_X = vector_union(G1_X, G2_X, N, N);
+    double_vector_t *G1_G3_X = vector_union(G1_X, G3_X, N, N);
+    double_vector_t *G2_G3_X = vector_union(G2_X, G3_X, N, N);
+
+    // Create a vector Y = my_function(x)
+    double_vector_t *Y = apply_function(X, my_function);
+    
+    double_vector_t *G1_Y = arrayCopy(Y, 0, N);
+    double_vector_t *G2_Y = arrayCopy(Y, N, N);
+    double_vector_t *G3_Y = arrayCopy(Y, N*2, N);
+
+    double_vector_t *G1_G2_Y = vector_union(G1_Y, G2_Y, N, N);
+    double_vector_t *G1_G3_Y = vector_union(G1_Y, G3_Y, N, N);
+    double_vector_t *G2_G3_Y = vector_union(G2_Y, G3_Y, N, N);
+    
+    couple deriver_E = deriver(X, Y);
+    
+    couple G1_G2 = calcule_AB(0.001, G1_G2_X, G1_G2_Y);
+    couple G1_G3 = calcule_AB(0.001, G1_G3_X, G1_G3_Y);
+    couple G2_G3 = calcule_AB(0.001, G2_G3_X, G2_G3_Y);
+
+    a = G1_G2;
+    double_vector_t *Y_line_G1_G2 = apply_function(X2, genLine);
+    a = G1_G3;
+    double_vector_t *Y_line_G1_G3 = apply_function(X2, genLine);
+    a = G2_G3;
+    double_vector_t *Y_line_G2_G3 = apply_function(X2, genLine);
+
+    
+    // Export our vectors into files
+    export_vector("../X.vec", X2);
+    export_vector("../G1_G2_Y.vec", Y_line_G1_G2);
+    export_vector("../G1_G3_Y.vec", Y_line_G1_G3);
+    export_vector("../G2_G3_Y.vec", Y_line_G2_G3);
+
+    export_vector("../G1_X.vec", G1_X);
+    export_vector("../G1_Y.vec", G1_Y);
+
+    export_vector("../G2_X.vec", G2_X);
+    export_vector("../G2_Y.vec", G2_Y);
+
+    export_vector("../G3_X.vec", G3_X);
+    export_vector("../G3_Y.vec", G3_Y);
+
+    // Print
+    printf("The line follow by points\n");
+    printf("  A: %f\n  B: %f\n\n", line_a, line_b);
+
+    printf("The théorique line\n");
+    printf("  A: %f\n  B: %f\n\n", deriver_E.a, deriver_E.b);
+
+    printf("The line by regrestion\n");
+    printf(" G1_G2 et tester sur G3\n");
+    printf("  A: %f\n  B: %f\n  Erreur quadratique: %f\n", G1_G2.a, G1_G2.b, erreur_quad(G1_G2 ,G1_G2_X, G1_G2_Y));
+    printf(" G1_G3 et tester sur G2\n");
+    printf("  A: %f\n  B: %f\n  Erreur quadratique: %f\n", G1_G3.a, G1_G3.b, erreur_quad(G1_G3 ,G1_G3_X, G1_G3_Y));
+    printf(" G2_G3 et tester sur G1\n");
+    printf("  A: %f\n  B: %f\n  Erreur quadratique: %f\n", G2_G3.a, G2_G3.b, erreur_quad(G2_G3 ,G2_G3_X, G2_G3_Y));
+
+    // Free our vectors
+    destroy_vector(&Y);
+    destroy_vector(&X);
+    destroy_vector(&X1);
+    destroy_vector(&X2);
+    destroy_vector(&Y_line_G1_G2);
+    destroy_vector(&Y_line_G1_G3);
+    destroy_vector(&Y_line_G2_G3);
+
+
+    destroy_vector(&G1_X);
+    destroy_vector(&G2_X);
+    destroy_vector(&G3_X);
+
+    destroy_vector(&G1_G2_X);
+    destroy_vector(&G1_G3_X);
+    destroy_vector(&G2_G3_X);
+
+    destroy_vector(&G1_Y);
+    destroy_vector(&G2_Y);
+    destroy_vector(&G3_Y);
+
+    destroy_vector(&G1_G2_Y);
+    destroy_vector(&G1_G3_Y);
+    destroy_vector(&G2_G3_Y);
+
+    printf("Vector generation succes\n");
+}
\ No newline at end of file

C/makefile

+
+CC=gcc -std=gnu11 -Wall -Wextra -g -fsanitize=leak -fsanitize=undefined
+LIBS=-lm
+
+example: vector.o vector.h opti.o opti.h main.c
+	$(CC) $^ $(LIBS) -o $@ 
+
+opti : vector.o vector.h opti.c opti.h
+	$(CC) $^ $(LIBS) -o $@ 
+
+vector.o: vector.c vector.h 
+	$(CC) -c $< -o $@
+
+clean:
+	rm -f *.vec *.o example
\ No newline at end of file

C/opti.c

+#include <stdint.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "vector.h"
+#include "opti.h"
+
+double my_random(double min, double max){
+    return (rand()/(RAND_MAX/(max-min))) - min;
+}
+
+
+couple gradiant(couple point, double_vector_t* X, double_vector_t* Y){
+    double somme_X = 0;
+    double somme_X2 = 0;
+    double somme_YX = 0;
+    double somme_Y = 0;
+    couple a;
+    for (int i = 0; i < (int)X->N;i++){
+        double x = X->components[i];
+        double y = Y->components[i];
+        somme_X += x;
+        somme_X2 += x*x;
+        somme_YX += x*y; 
+        somme_Y += y;
+    }
+    
+    int N = X->N;
+    a.a =  point.a * somme_X2 + point.b * somme_X - somme_YX;
+    a.b = point.a * somme_X + N * point.b - somme_Y;
+    return a;
+}
+
+couple ab1(couple current, double lambda, couple deriver_E){
+    couple a1;
+    a1.a = current.a - (lambda * deriver_E.a);
+    a1.b = current.b - (lambda * deriver_E.b);
+    return a1;
+}
+
+couple calcule_AB(double epsilone, double_vector_t* X, double_vector_t* Y){
+    couple current;
+    couple next;
+    
+    // init a et b 0 my_random
+    double r = my_random(0, 1);
+    next.a = r;
+    r = my_random(0, 1);
+    next.b = r;
+    double norme = 1.0;
+
+    r = my_random(0,0.002); // r devient lambda
+    while (norme >=  epsilone){
+        current = next;
+        
+        // a et b i+1
+        next = ab1(current, r, gradiant(current, X, Y));
+        
+        // calcul norme
+        double a = next.a - current.a;
+        double b = next.b - current.b;
+        norme =  sqrt(((a*a) + (b*b)));
+        // printf("A: %f, B: %f; Norme: %f\n", current.a, current.b, norme);
+    }
+   
+    return current;
+}   
+
+couple deriver(double_vector_t* X, double_vector_t* Y){
+    double somme_X = 0;
+    double somme_X2 = 0;
+    double somme_YX = 0;
+    double somme_Y = 0;
+    couple a;
+    for (int i = 0; i < (int)X->N;i++){
+        double x = X->components[i];
+        double y = Y->components[i];
+        somme_X += x;
+        somme_X2 += x*x;
+        somme_YX += x*y; 
+        somme_Y += y;
+    }
+    
+    int N = X->N;
+    a.b = (somme_Y * somme_X2 - somme_YX* somme_X)/(N * somme_X2 - somme_X * somme_X);
+    a.a =  (somme_Y - N * a.b)/ somme_X;
+    return a;
+}
+
+double erreur_quad(couple a, double_vector_t *X, double_vector_t *Y){
+	double somme = 0;
+	for (int i = 0; i < X->N; i++){
+		somme += (a.a * X->components[i] + a.b - Y->components[i])*(a.a * X->components[i] + a.b - Y->components[i]);
+	}
+}
\ No newline at end of file

C/opti.h

+#ifndef _OPTI_H_
+#define _OPTI_H_
+
+#include <stdint.h>
+#include "vector.h"
+
+typedef struct couple_t{
+	double a, b;
+}couple;
+
+double my_random(double min, double max);
+
+couple gradiant(couple point, double_vector_t* X, double_vector_t* Y);
+
+couple ab1(couple current, double lambda, couple deriver_E);
+
+couple calcule_AB(double epsilone, double_vector_t* X, double_vector_t* Y);
+
+couple deriver(double_vector_t* X, double_vector_t* Y);
+
+double erreur_quad(couple a, double_vector_t *X, double_vector_t *Y);
+#endif
\ No newline at end of file

C/vector.c

+#include <stdio.h>
+#include <stdlib.h>
+
+#include "vector.h"
+
+/** @brief 
+ *  Compute the endianness used by
+ *  the architecture.
+ *
+ *  @return 1 if little-endian, 0 if big-endian
+ */
+uint8_t get_endianness()
+{
+    uint32_t endianness = 0x01020304;
+    // Return the endianness by accessing the first byte in memory
+    // which should be 1 if big-endian and 4 if little-endian
+    return *((uint8_t *)(&endianness)) == 4;
+}
+
+/** @brief 
+ *  Create a vector of a given dimension.
+ * 
+ *  @param N The number of dimensions.
+ *  @return A dynamically allocated vector
+ */
+double_vector_t *init_vector(uint32_t N)
+{
+    double_vector_t *vec = malloc(sizeof(double_vector_t));
+    vec->components = malloc(N * sizeof(double));
+    vec->N = N;
+    if (vec == NULL)
+    {
+        perror("Can't allocate memory");
+        exit(EXIT_FAILURE);
+    }
+    return vec;
+}
+
+/** @brief 
+ *  Create a vector of a given dimension,
+ *  with values from 0 to N excluded.
+ * 
+ *  @param N The number of dimensions.
+ *  @return A dynamically allocated vector : [0,1..N-1]
+ */
+double_vector_t *iota(uint32_t N)
+{
+    double_vector_t *vec = init_vector(N);
+    for (uint32_t i = 0; i < N; i++)
+    {
+        vec->components[i] = i;
+    }
+    return vec;
+}
+
+/** @brief 
+ *  Apply a 1d function element-wise
+ *  to a given vector, and return the
+ *  result in a new vector.
+ * 
+ *  @param vec The argument vector
+ *  @param f   The 1d function to apply
+ *  @return A dynamically allocated vector : f(X)
+ */
+double_vector_t *apply_function(double_vector_t *vec, double_function_t f)
+{
+    double_vector_t *res = init_vector(vec->N);
+    for (uint32_t i = 0; i < vec->N; i++)
+    {
+        res->components[i] = f(vec->components[i]);
+    }
+    return res;
+}
+
+/** @brief 
+ *  Export a vector into a file.
+ * 
+ *  @param filename The name of the output file
+ *  @param vec      The vector to export
+ */
+void export_vector(const char *filename, double_vector_t *vec)
+{
+    FILE *output = fopen(filename, "w");
+
+    vector_metadata_t metadata;
+    metadata.endianness = get_endianness();
+    metadata.size_of_a_component = sizeof(double);
+    metadata.number_of_component = vec->N;
+
+    fwrite(&metadata, sizeof(vector_metadata_t), 1, output);
+    fwrite(vec->components,
+           metadata.size_of_a_component,
+           metadata.number_of_component,
+           output);
+
+    fclose(output);
+}
+
+/** @brief 
+ *  Free a vector.
+ * 
+ *  @param vec A double pointer on a vector
+ */
+void destroy_vector(double_vector_t **vec)
+{
+    free((*vec)->components);
+    free(*vec);
+    *vec = NULL;
+}
+
+double_vector_t *arrayCopy(double_vector_t *src, int start ,uint32_t N)
+{
+    double_vector_t *vec = init_vector(N);
+    for (uint32_t i = 0; i < N; i++){
+        vec->components[i] = src->components[(start+i)];
+    }
+    return vec;
+}
+
+double_vector_t *vector_union(double_vector_t *src1, double_vector_t *src2, uint32_t N1 ,uint32_t N2)
+{
+    double_vector_t *vec = init_vector((N1+N2));
+    for (uint32_t i = 0; i < N1 || i < N2; i++){
+        if (i < N1 ) vec->components[i] = src1->components[i]*1;
+        if (i < N2 ) vec->components[(i+N1)] = src2->components[i]*1;
+    }    
+    return vec;
+}
\ No newline at end of file

C/vector.h

+#ifndef _DOUBLE_VECTOR_H_
+#define _DOUBLE_VECTOR_H_
+
+#include <stdint.h>
+
+typedef struct double_vector
+{
+    uint32_t N; // The dimmension of the vector
+    double *components;
+} double_vector_t;
+// Function pointer, example : double f(double x);
+typedef double (*double_function_t)(double);
+
+/*
+* The attribute "packed" tells the compiler,
+* that the struct should be stored in memory 
+* without padding. It's highly recommended,
+* if we want to serialize the structure.
+* (for example to store it in a file)
+*/
+typedef struct vector_metadata
+{
+    uint8_t endianness;          // 1 = little, 0 = big
+    uint8_t size_of_a_component; // in bytes
+    uint32_t number_of_component;
+} __attribute__((packed)) vector_metadata_t;
+
+/** @brief 
+ *  Create a vector of a given dimension.
+ * 
+ *  @param N The number of dimensions.
+ *  @return A dynamically allocated vector
+ */
+double_vector_t *init_vector(uint32_t N);
+/** @brief 
+ *  Create a vector of a given dimension,
+ *  with values from 0 to N excluded.
+ * 
+ *  @param N The number of dimensions.
+ *  @return A dynamically allocated vector : [0,1..N-1]
+ */
+double_vector_t *iota(uint32_t N);
+/** @brief 
+ *  Apply a 1d function element-wise
+ *  to a given vector, and return the
+ *  result in a new vector.
+ * 
+ *  @param vec The argument vector
+ *  @param f   The 1d function to apply
+ *  @return A dynamically allocated vector : f(X)
+ */
+double_vector_t *apply_function(double_vector_t *vec, double_function_t f);
+/** @brief 
+ *  Export a vector into a file.
+ * 
+ *  @param filename The name of the output file
+ *  @param vec      The vector to export
+ */
+void export_vector(const char *filename, double_vector_t *vec);
+/** @brief 
+ *  Free a vector.
+ * 
+ *  @param vec A double pointer on a vector
+ */
+void destroy_vector(double_vector_t **vec);
+
+double_vector_t *arrayCopy(double_vector_t *src, int start ,uint32_t N);
+
+double_vector_t *vector_union(double_vector_t *src1, double_vector_t *src2, uint32_t N1 ,uint32_t N2);
+
+#endif
\ No newline at end of file

Python/example.py

+from matplotlib import pyplot as plt
+from load_vec import load_vector
+
+X = load_vector("../X.vec")
+G1_G2_Y = load_vector("../G1_G2_Y.vec")
+G1_G3_Y = load_vector("../G1_G3_Y.vec")
+G2_G3_Y = load_vector("../G2_G3_Y.vec")
+
+G1_X = load_vector("../G1_X.vec")
+G1_Y = load_vector("../G1_Y.vec")
+
+G2_X = load_vector("../G2_X.vec")
+G2_Y = load_vector("../G2_Y.vec")
+
+G3_X = load_vector("../G3_X.vec")
+G3_Y = load_vector("../G3_Y.vec")
+
+plt.plot(X, G1_G2_Y, color="r", label="Line G1_G2")
+plt.plot(X, G1_G3_Y, color="y", label="Line G1_G3")
+plt.plot(X, G2_G3_Y, color="k", label="Line G2_G3")
+
+plt.scatter(G1_X, G1_Y, marker='.', label="G1")
+plt.scatter(G2_X, G2_Y, marker='.', label="G2")
+plt.scatter(G3_X, G3_Y, marker='.', label="G3")
+
+
+plt.title("My data")
+plt.xlabel("X")
+plt.ylabel("Y")
+plt.legend(loc="upper left")
+
+plt.show()

Python/load_vec.py

+import numpy as np
+from typing import Tuple
+METADATA_SIZE = 6
+
+
+def _parse_metadata(metadata: bytes) -> Tuple[type, int]:
+    """
+    Parse the metadata for a vec file.
+
+    Parameters:
+    metadata (bytes): The metadata bytes
+
+    Returns:
+    (type, int): The type and the number of component of the vector
+    """
+
+    little_endian = bool(metadata[0])
+    endianness = 'little' if little_endian else 'big'
+
+    size_of_components = int(metadata[1])
+    # For now we only consider two types
+    datatype = np.float64 if size_of_components == 8 else np.float
+
+    # Recover our 32 bit integer specifying the endianness
+    nb_components = int.from_bytes(metadata[2:], endianness)
+
+    return datatype, nb_components
+
+
+def load_vector(filename: str) -> np.ndarray:
+    """
+    Load a vector from a file.
+
+    Parameters:
+    filename (str): The name of the file containing the vector
+
+    Returns:
+    np.ndarray: The vector
+
+    """
+    file = open(filename, 'rb')
+    # Read our metadata struct
+    metadata = file.read(METADATA_SIZE)
+
+    datatype, nb_components = _parse_metadata(metadata)
+
+    array = np.fromfile(file, dtype=datatype, count=nb_components)
+
+    file.close()
+    return array

gen_vector.sh

+#! /bin/bash
+
+cd C/
+make 2> /dev/null
+./example
+cd ../
\ No newline at end of file

show.sh

+#! /bin/bash
+
+cd Python
+python3 example.py
+cd ../
\ No newline at end of file