diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..c3a8bb1202f6a611e54937fa23cb13e5384dfd7b --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +*.vec +.vscode +__pycache__ +example +*.o \ No newline at end of file diff --git a/C/main.c b/C/main.c new file mode 100644 index 0000000000000000000000000000000000000000..81eff94910da1f3112dea20cedc98f2253847e5b --- /dev/null +++ b/C/main.c @@ -0,0 +1,145 @@ +#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 diff --git a/C/makefile b/C/makefile new file mode 100644 index 0000000000000000000000000000000000000000..a3f77cee460983d7390e2a04e2957bb5ebb01bed --- /dev/null +++ b/C/makefile @@ -0,0 +1,15 @@ + +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 diff --git a/C/opti.c b/C/opti.c new file mode 100644 index 0000000000000000000000000000000000000000..181f6fc8cf09e2c3d83f2ca235d276a3d9fbf64d --- /dev/null +++ b/C/opti.c @@ -0,0 +1,95 @@ +#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 diff --git a/C/opti.h b/C/opti.h new file mode 100644 index 0000000000000000000000000000000000000000..ece6ccbe70aaf57c18d1a2cf4f1fdf5dcfb83663 --- /dev/null +++ b/C/opti.h @@ -0,0 +1,22 @@ +#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 diff --git a/C/vector.c b/C/vector.c new file mode 100644 index 0000000000000000000000000000000000000000..109a20fbadf01e85203e16e0b00dd600f6a98c62 --- /dev/null +++ b/C/vector.c @@ -0,0 +1,128 @@ +#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 diff --git a/C/vector.h b/C/vector.h new file mode 100644 index 0000000000000000000000000000000000000000..cfc0dfd444a593b049b01a4b9230e57e2c7f20be --- /dev/null +++ b/C/vector.h @@ -0,0 +1,71 @@ +#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 diff --git a/Python/example.py b/Python/example.py new file mode 100644 index 0000000000000000000000000000000000000000..3cba2092358736e1bf382c90afafcbebda168559 --- /dev/null +++ b/Python/example.py @@ -0,0 +1,32 @@ +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() diff --git a/Python/load_vec.py b/Python/load_vec.py new file mode 100644 index 0000000000000000000000000000000000000000..17880a18b4fccf86b61cad3cd8e90b1aead07926 --- /dev/null +++ b/Python/load_vec.py @@ -0,0 +1,50 @@ +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 diff --git a/gen_vector.sh b/gen_vector.sh new file mode 100755 index 0000000000000000000000000000000000000000..2eb8bcb53c94eceb7cb2a3f9090609d874276501 --- /dev/null +++ b/gen_vector.sh @@ -0,0 +1,6 @@ +#! /bin/bash + +cd C/ +make 2> /dev/null +./example +cd ../ \ No newline at end of file diff --git a/show.sh b/show.sh new file mode 100755 index 0000000000000000000000000000000000000000..294de535ccae9dc3641bdf4944eed5704b0a941d --- /dev/null +++ b/show.sh @@ -0,0 +1,5 @@ +#! /bin/bash + +cd Python +python3 example.py +cd ../ \ No newline at end of file