mirrors
/
TheAlgorithms-C
mirror of https://github.com/TheAlgorithms/C
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

fix link references from fork to main repo

This commit is contained in:
Krishna Vedala 2020-06-28 11:22:42 -04:00
parent d5d68d5842
commit b833e27964
No known key found for this signature in database
GPG Key ID: BA19ACF8FC8792F7
7 changed files with 212 additions and 216 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
README.md
View File

@ -1,15 +1,15 @@
# The Algorithms - C # {#mainpage} # The Algorithms - C # {#mainpage}
[![Gitpod Ready-to-Code](https://img.shields.io/badge/Gitpod-Ready--to--Code-blue?logo=gitpod)](https://gitpod.io/#https://github.com/kvedala/C) [![Gitpod Ready-to-Code](https://img.shields.io/badge/Gitpod-Ready--to--Code-blue?logo=gitpod)](https://gitpod.io/#https://github.com/TheAlgorithms/C)
[![Gitter chat](https://img.shields.io/badge/Chat-Gitter-ff69b4.svg?label=Chat&logo=gitter&style=flat-square)](https://gitter.im/TheAlgorithms) [![Gitter chat](https://img.shields.io/badge/Chat-Gitter-ff69b4.svg?label=Chat&logo=gitter&style=flat-square)](https://gitter.im/TheAlgorithms)
[![contributions welcome](https://img.shields.io/static/v1.svg?label=Contributions&message=Welcome&color=0059b3&style=flat-square)](https://github.com/kvedala/C-Plus-Plus/blob/master/CONTRIBUTING.md)  [![contributions welcome](https://img.shields.io/static/v1.svg?label=Contributions&message=Welcome&color=0059b3&style=flat-square)](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/CONTRIBUTING.md) 
![GitHub repo size](https://img.shields.io/github/repo-size/kvedala/C-Plus-Plus?color=red&style=flat-square) ![GitHub repo size](https://img.shields.io/github/repo-size/TheAlgorithms/C-Plus-Plus?color=red&style=flat-square)
![GitHub closed pull requests](https://img.shields.io/github/issues-pr-closed/kvedala/C?color=green&style=flat-square) ![GitHub closed pull requests](https://img.shields.io/github/issues-pr-closed/TheAlgorithms/C?color=green&style=flat-square)
![Doxygen CI](https://github.com/kvedala/C/workflows/Doxygen%20CI/badge.svg) ![Doxygen CI](https://github.com/TheAlgorithms/C/workflows/Doxygen%20CI/badge.svg)
![Awesome CI Workflow](https://github.com/kvedala/C/workflows/Awesome%20CI%20Workflow/badge.svg) ![Awesome CI Workflow](https://github.com/TheAlgorithms/C/workflows/Awesome%20CI%20Workflow/badge.svg)
[Online Documentation](https://kvedala.github.io/C). [Online Documentation](https://TheAlgorithms.github.io/C).
Click on [Files menu](https://kvedala.github.io/C/files.html) to see the list of all the files documented with the code. Click on [Files menu](https://TheAlgorithms.github.io/C/files.html) to see the list of all the files documented with the code.
All the code can be executed and tested online: [![using Google Colab Notebook](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/gist/kvedala/27f1b0b6502af935f6917673ec43bcd7/plot-durand_kerner-log.ipynb) All the code can be executed and tested online: [![using Google Colab Notebook](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/gist/kvedala/27f1b0b6502af935f6917673ec43bcd7/plot-durand_kerner-log.ipynb)

183
machine_learning/kohonen_som_topology.c
View File

@ -11,7 +11,7 @@
* data points in a much higher dimesional space by creating an equivalent in a * data points in a much higher dimesional space by creating an equivalent in a
* 2-dimensional space. * 2-dimensional space.
* <img alt="Trained topological maps for the test cases in the program" * <img alt="Trained topological maps for the test cases in the program"
* src="https://raw.githubusercontent.com/kvedala/C/docs/images/machine_learning/kohonen/2D_Kohonen_SOM.svg" * src="https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/machine_learning/kohonen/2D_Kohonen_SOM.svg"
* /> * />
* \warning MSVC 2019 compiler generates code that does not execute as expected. * \warning MSVC 2019 compiler generates code that does not execute as expected.
* However, MinGW, Clang for GCC and Clang for MSVC compilers on windows perform * However, MinGW, Clang for GCC and Clang for MSVC compilers on windows perform
@ -23,18 +23,18 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <time.h> #include <time.h>
#ifdef _OPENMP // check if OpenMP based parallellization is available #ifdef _OPENMP // check if OpenMP based parallellization is available
#include <omp.h> #include <omp.h>
#endif #endif
#ifndef max #ifndef max
#define max(a, b) \ #define max(a, b) \
(((a) > (b)) ? (a) : (b)) /**< shorthand for maximum value \ (((a) > (b)) ? (a) : (b)) /**< shorthand for maximum value \
*/ */
#endif #endif
#ifndef min #ifndef min
#define min(a, b) \ #define min(a, b) \
(((a) < (b)) ? (a) : (b)) /**< shorthand for minimum value \ (((a) < (b)) ? (a) : (b)) /**< shorthand for minimum value \
*/ */
#endif #endif
@ -98,7 +98,7 @@ int save_2d_data(const char *fname, double **X, int num_points,
int num_features) int num_features)
{ {
FILE *fp = fopen(fname, "wt"); FILE *fp = fopen(fname, "wt");
if (!fp) // error with fopen if (!fp) // error with fopen
{ {
char msg[120]; char msg[120];
sprintf(msg, "File error (%s): ", fname); sprintf(msg, "File error (%s): ", fname);
@ -106,16 +106,16 @@ int save_2d_data(const char *fname, double **X, int num_points,
return -1; return -1;
} }
for (int i = 0; i < num_points; i++) // for each point in the array for (int i = 0; i < num_points; i++) // for each point in the array
{ {
for (int j = 0; j < num_features; j++) // for each feature in the array for (int j = 0; j < num_features; j++) // for each feature in the array
{ {
fprintf(fp, "%.4g", X[i][j]); // print the feature value fprintf(fp, "%.4g", X[i][j]); // print the feature value
if (j < num_features - 1) // if not the last feature if (j < num_features - 1) // if not the last feature
fputc(',', fp); // suffix comma fputc(',', fp); // suffix comma
} }
if (i < num_points - 1) // if not the last row if (i < num_points - 1) // if not the last row
fputc('\n', fp); // start a new line fputc('\n', fp); // start a new line
} }
fclose(fp); fclose(fp);
return 0; return 0;
@ -134,7 +134,7 @@ int save_2d_data(const char *fname, double **X, int num_points,
int save_u_matrix(const char *fname, struct array_3d *W) int save_u_matrix(const char *fname, struct array_3d *W)
{ {
FILE *fp = fopen(fname, "wt"); FILE *fp = fopen(fname, "wt");
if (!fp) // error with fopen if (!fp) // error with fopen
{ {
char msg[120]; char msg[120];
sprintf(msg, "File error (%s): ", fname); sprintf(msg, "File error (%s): ", fname);
@ -144,9 +144,9 @@ int save_u_matrix(const char *fname, struct array_3d *W)
int R = max(W->dim1 >> 3, 2); /* neighborhood range */ int R = max(W->dim1 >> 3, 2); /* neighborhood range */
for (int i = 0; i < W->dim1; i++) // for each x for (int i = 0; i < W->dim1; i++) // for each x
{ {
for (int j = 0; j < W->dim2; j++) // for each y for (int j = 0; j < W->dim2; j++) // for each y
{ {
double distance = 0.f; double distance = 0.f;
int k; int k;
@ -159,12 +159,12 @@ int save_u_matrix(const char *fname, struct array_3d *W)
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp parallel for reduction(+ : distance) #pragma omp parallel for reduction(+ : distance)
#endif #endif
for (l = from_x; l < to_x; l++) // scan neighborhoor in x for (l = from_x; l < to_x; l++) // scan neighborhoor in x
{ {
for (int m = from_y; m < to_y; m++) // scan neighborhood in y for (int m = from_y; m < to_y; m++) // scan neighborhood in y
{ {
double d = 0.f; double d = 0.f;
for (k = 0; k < W->dim3; k++) // for each feature for (k = 0; k < W->dim3; k++) // for each feature
{ {
double *w1 = data_3d(W, i, j, k); double *w1 = data_3d(W, i, j, k);
double *w2 = data_3d(W, l, m, k); double *w2 = data_3d(W, l, m, k);
@ -176,13 +176,13 @@ int save_u_matrix(const char *fname, struct array_3d *W)
} }
} }
distance /= R * R; // mean distance from neighbors distance /= R * R; // mean distance from neighbors
fprintf(fp, "%.4g", distance); // print the mean separation fprintf(fp, "%.4g", distance); // print the mean separation
if (j < W->dim2 - 1) // if not the last column if (j < W->dim2 - 1) // if not the last column
fputc(',', fp); // suffix comma fputc(',', fp); // suffix comma
} }
if (i < W->dim1 - 1) // if not the last row if (i < W->dim1 - 1) // if not the last row
fputc('\n', fp); // start a new line fputc('\n', fp); // start a new line
} }
fclose(fp); fclose(fp);
return 0; return 0;
@ -198,14 +198,14 @@ int save_u_matrix(const char *fname, struct array_3d *W)
*/ */
void get_min_2d(double **X, int N, double *val, int *x_idx, int *y_idx) void get_min_2d(double **X, int N, double *val, int *x_idx, int *y_idx)
{ {
val[0] = INFINITY; // initial min value val[0] = INFINITY; // initial min value
for (int i = 0; i < N; i++) // traverse each x-index for (int i = 0; i < N; i++) // traverse each x-index
{ {
for (int j = 0; j < N; j++) // traverse each y-index for (int j = 0; j < N; j++) // traverse each y-index
{ {
if (X[i][j] < val[0]) // if a lower value is found if (X[i][j] < val[0]) // if a lower value is found
{ // save the value and its index { // save the value and its index
x_idx[0] = i; x_idx[0] = i;
y_idx[0] = j; y_idx[0] = j;
val[0] = X[i][j]; val[0] = X[i][j];
@ -314,7 +314,7 @@ void kohonen_som(double **X, struct array_3d *W, int num_samples,
for (int i = 0; i < num_out; i++) for (int i = 0; i < num_out; i++)
D[i] = (double *)malloc(num_out * sizeof(double)); D[i] = (double *)malloc(num_out * sizeof(double));
double dmin = 1.f; // average minimum distance of all samples double dmin = 1.f; // average minimum distance of all samples
// Loop alpha from 1 to slpha_min // Loop alpha from 1 to slpha_min
for (double alpha = 1.f; alpha > alpha_min && dmin > 1e-3; for (double alpha = 1.f; alpha > alpha_min && dmin > 1e-3;
@ -339,8 +339,7 @@ void kohonen_som(double **X, struct array_3d *W, int num_samples,
} }
putchar('\n'); putchar('\n');
for (int i = 0; i < num_out; i++) for (int i = 0; i < num_out; i++) free(D[i]);
free(D[i]);
free(D); free(D);
} }
@ -356,15 +355,15 @@ void kohonen_som(double **X, struct array_3d *W, int num_samples,
*/ */
void test_2d_classes(double *const *data, int N) void test_2d_classes(double *const *data, int N)
{ {
const double R = 0.3; // radius of cluster const double R = 0.3; // radius of cluster
int i; int i;
const int num_classes = 4; const int num_classes = 4;
const double centres[][2] = { const double centres[][2] = {
// centres of each class cluster // centres of each class cluster
{.5, .5}, // centre of class 1 {.5, .5}, // centre of class 1
{.5, -.5}, // centre of class 2 {.5, -.5}, // centre of class 2
{-.5, .5}, // centre of class 3 {-.5, .5}, // centre of class 3
{-.5, -.5} // centre of class 4 {-.5, -.5} // centre of class 4
}; };
#ifdef _OPENMP #ifdef _OPENMP
@ -372,7 +371,8 @@ void test_2d_classes(double *const *data, int N)
#endif #endif
for (i = 0; i < N; i++) for (i = 0; i < N; i++)
{ {
int class = rand() % num_classes; // select a random class for the point int class =
rand() % num_classes; // select a random class for the point
// create random coordinates (x,y,z) around the centre of the class // create random coordinates (x,y,z) around the centre of the class
data[i][0] = _random(centres[class][0] - R, centres[class][0] + R); data[i][0] = _random(centres[class][0] - R, centres[class][0] + R);
@ -397,7 +397,7 @@ void test1()
{ {
int j, N = 300; int j, N = 300;
int features = 2; int features = 2;
int num_out = 30; // image size - N x N int num_out = 30; // image size - N x N
// 2D space, hence size = number of rows * 2 // 2D space, hence size = number of rows * 2
double **X = (double **)malloc(N * sizeof(double *)); double **X = (double **)malloc(N * sizeof(double *));
@ -408,13 +408,13 @@ void test1()
W.dim2 = num_out; W.dim2 = num_out;
W.dim3 = features; W.dim3 = features;
W.data = (double *)malloc(num_out * num_out * features * W.data = (double *)malloc(num_out * num_out * features *
sizeof(double)); // assign rows sizeof(double)); // assign rows
for (int i = 0; i < max(num_out, N); i++) // loop till max(N, num_out) for (int i = 0; i < max(num_out, N); i++) // loop till max(N, num_out)
{ {
if (i < N) // only add new arrays if i < N if (i < N) // only add new arrays if i < N
X[i] = (double *)malloc(features * sizeof(double)); X[i] = (double *)malloc(features * sizeof(double));
if (i < num_out) // only add new arrays if i < num_out if (i < num_out) // only add new arrays if i < num_out
{ {
for (int k = 0; k < num_out; k++) for (int k = 0; k < num_out; k++)
{ {
@ -431,14 +431,13 @@ void test1()
} }
} }
test_2d_classes(X, N); // create test data around circumference of a circle test_2d_classes(X, N); // create test data around circumference of a circle
save_2d_data("test1.csv", X, N, features); // save test data points save_2d_data("test1.csv", X, N, features); // save test data points
save_u_matrix("w11.csv", &W); // save initial random weights save_u_matrix("w11.csv", &W); // save initial random weights
kohonen_som(X, &W, N, features, num_out, 1e-4); // train the SOM kohonen_som(X, &W, N, features, num_out, 1e-4); // train the SOM
save_u_matrix("w12.csv", &W); // save the resultant weights save_u_matrix("w12.csv", &W); // save the resultant weights
for (int i = 0; i < N; i++) for (int i = 0; i < N; i++) free(X[i]);
free(X[i]);
free(X); free(X);
free(W.data); free(W.data);
} }
@ -455,15 +454,15 @@ void test1()
*/ */
void test_3d_classes1(double *const *data, int N) void test_3d_classes1(double *const *data, int N)
{ {
const double R = 0.2; // radius of cluster const double R = 0.2; // radius of cluster
int i; int i;
const int num_classes = 4; const int num_classes = 4;
const double centres[][3] = { const double centres[][3] = {
// centres of each class cluster // centres of each class cluster
{.5, .5, .5}, // centre of class 1 {.5, .5, .5}, // centre of class 1
{.5, -.5, -.5}, // centre of class 2 {.5, -.5, -.5}, // centre of class 2
{-.5, .5, .5}, // centre of class 3 {-.5, .5, .5}, // centre of class 3
{-.5, -.5 - .5} // centre of class 4 {-.5, -.5 - .5} // centre of class 4
}; };
#ifdef _OPENMP #ifdef _OPENMP
@ -471,7 +470,8 @@ void test_3d_classes1(double *const *data, int N)
#endif #endif
for (i = 0; i < N; i++) for (i = 0; i < N; i++)
{ {
int class = rand() % num_classes; // select a random class for the point int class =
rand() % num_classes; // select a random class for the point
// create random coordinates (x,y,z) around the centre of the class // create random coordinates (x,y,z) around the centre of the class
data[i][0] = _random(centres[class][0] - R, centres[class][0] + R); data[i][0] = _random(centres[class][0] - R, centres[class][0] + R);
@ -497,7 +497,7 @@ void test2()
{ {
int j, N = 500; int j, N = 500;
int features = 3; int features = 3;
int num_out = 30; // image size - N x N int num_out = 30; // image size - N x N
// 3D space, hence size = number of rows * 3 // 3D space, hence size = number of rows * 3
double **X = (double **)malloc(N * sizeof(double *)); double **X = (double **)malloc(N * sizeof(double *));
@ -508,13 +508,13 @@ void test2()
W.dim2 = num_out; W.dim2 = num_out;
W.dim3 = features; W.dim3 = features;
W.data = (double *)malloc(num_out * num_out * features * W.data = (double *)malloc(num_out * num_out * features *
sizeof(double)); // assign rows sizeof(double)); // assign rows
for (int i = 0; i < max(num_out, N); i++) // loop till max(N, num_out) for (int i = 0; i < max(num_out, N); i++) // loop till max(N, num_out)
{ {
if (i < N) // only add new arrays if i < N if (i < N) // only add new arrays if i < N
X[i] = (double *)malloc(features * sizeof(double)); X[i] = (double *)malloc(features * sizeof(double));
if (i < num_out) // only add new arrays if i < num_out if (i < num_out) // only add new arrays if i < num_out
{ {
for (int k = 0; k < num_out; k++) for (int k = 0; k < num_out; k++)
{ {
@ -522,7 +522,7 @@ void test2()
#pragma omp for #pragma omp for
#endif #endif
for (j = 0; j < features; j++) for (j = 0; j < features; j++)
{ // preallocate with random initial weights { // preallocate with random initial weights
double *w = data_3d(&W, i, k, j); double *w = data_3d(&W, i, k, j);
w[0] = _random(-5, 5); w[0] = _random(-5, 5);
} }
@ -530,14 +530,13 @@ void test2()
} }
} }
test_3d_classes1(X, N); // create test data test_3d_classes1(X, N); // create test data
save_2d_data("test2.csv", X, N, features); // save test data points save_2d_data("test2.csv", X, N, features); // save test data points
save_u_matrix("w21.csv", &W); // save initial random weights save_u_matrix("w21.csv", &W); // save initial random weights
kohonen_som(X, &W, N, features, num_out, 1e-4); // train the SOM kohonen_som(X, &W, N, features, num_out, 1e-4); // train the SOM
save_u_matrix("w22.csv", &W); // save the resultant weights save_u_matrix("w22.csv", &W); // save the resultant weights
for (int i = 0; i < N; i++) for (int i = 0; i < N; i++) free(X[i]);
free(X[i]);
free(X); free(X);
free(W.data); free(W.data);
} }
@ -554,19 +553,19 @@ void test2()
*/ */
void test_3d_classes2(double *const *data, int N) void test_3d_classes2(double *const *data, int N)
{ {
const double R = 0.2; // radius of cluster const double R = 0.2; // radius of cluster
int i; int i;
const int num_classes = 8; const int num_classes = 8;
const double centres[][3] = { const double centres[][3] = {
// centres of each class cluster // centres of each class cluster
{.5, .5, .5}, // centre of class 1 {.5, .5, .5}, // centre of class 1
{.5, .5, -.5}, // centre of class 2 {.5, .5, -.5}, // centre of class 2
{.5, -.5, .5}, // centre of class 3 {.5, -.5, .5}, // centre of class 3
{.5, -.5, -.5}, // centre of class 4 {.5, -.5, -.5}, // centre of class 4
{-.5, .5, .5}, // centre of class 5 {-.5, .5, .5}, // centre of class 5
{-.5, .5, -.5}, // centre of class 6 {-.5, .5, -.5}, // centre of class 6
{-.5, -.5, .5}, // centre of class 7 {-.5, -.5, .5}, // centre of class 7
{-.5, -.5, -.5} // centre of class 8 {-.5, -.5, -.5} // centre of class 8
}; };
#ifdef _OPENMP #ifdef _OPENMP
@ -574,7 +573,8 @@ void test_3d_classes2(double *const *data, int N)
#endif #endif
for (i = 0; i < N; i++) for (i = 0; i < N; i++)
{ {
int class = rand() % num_classes; // select a random class for the point int class =
rand() % num_classes; // select a random class for the point
// create random coordinates (x,y,z) around the centre of the class // create random coordinates (x,y,z) around the centre of the class
data[i][0] = _random(centres[class][0] - R, centres[class][0] + R); data[i][0] = _random(centres[class][0] - R, centres[class][0] + R);
@ -609,13 +609,13 @@ void test3()
W.dim2 = num_out; W.dim2 = num_out;
W.dim3 = features; W.dim3 = features;
W.data = (double *)malloc(num_out * num_out * features * W.data = (double *)malloc(num_out * num_out * features *
sizeof(double)); // assign rows sizeof(double)); // assign rows
for (int i = 0; i < max(num_out, N); i++) // loop till max(N, num_out) for (int i = 0; i < max(num_out, N); i++) // loop till max(N, num_out)
{ {
if (i < N) // only add new arrays if i < N if (i < N) // only add new arrays if i < N
X[i] = (double *)malloc(features * sizeof(double)); X[i] = (double *)malloc(features * sizeof(double));
if (i < num_out) // only add new arrays if i < num_out if (i < num_out) // only add new arrays if i < num_out
{ {
for (int k = 0; k < num_out; k++) for (int k = 0; k < num_out; k++)
{ {
@ -632,14 +632,13 @@ void test3()
} }
} }
test_3d_classes2(X, N); // create test data around the lamniscate test_3d_classes2(X, N); // create test data around the lamniscate
save_2d_data("test3.csv", X, N, features); // save test data points save_2d_data("test3.csv", X, N, features); // save test data points
save_u_matrix("w31.csv", &W); // save initial random weights save_u_matrix("w31.csv", &W); // save initial random weights
kohonen_som(X, &W, N, features, num_out, 0.01); // train the SOM kohonen_som(X, &W, N, features, num_out, 0.01); // train the SOM
save_u_matrix("w32.csv", &W); // save the resultant weights save_u_matrix("w32.csv", &W); // save the resultant weights
for (int i = 0; i < N; i++) for (int i = 0; i < N; i++) free(X[i]);
free(X[i]);
free(X); free(X);
free(W.data); free(W.data);
} }

132
machine_learning/kohonen_som_trace.c
View File

@ -16,18 +16,18 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <time.h> #include <time.h>
#ifdef _OPENMP // check if OpenMP based parallellization is available #ifdef _OPENMP // check if OpenMP based parallellization is available
#include <omp.h> #include <omp.h>
#endif #endif
#ifndef max #ifndef max
#define max(a, b) \ #define max(a, b) \
(((a) > (b)) ? (a) : (b)) /**< shorthand for maximum value \ (((a) > (b)) ? (a) : (b)) /**< shorthand for maximum value \
*/ */
#endif #endif
#ifndef min #ifndef min
#define min(a, b) \ #define min(a, b) \
(((a) < (b)) ? (a) : (b)) /**< shorthand for minimum value \ (((a) < (b)) ? (a) : (b)) /**< shorthand for minimum value \
*/ */
#endif #endif
@ -64,7 +64,7 @@ int save_nd_data(const char *fname, double **X, int num_points,
int num_features) int num_features)
{ {
FILE *fp = fopen(fname, "wt"); FILE *fp = fopen(fname, "wt");
if (!fp) // error with fopen if (!fp) // error with fopen
{ {
char msg[120]; char msg[120];
sprintf(msg, "File error (%s): ", fname); sprintf(msg, "File error (%s): ", fname);
@ -72,16 +72,16 @@ int save_nd_data(const char *fname, double **X, int num_points,
return -1; return -1;
} }
for (int i = 0; i < num_points; i++) // for each point in the array for (int i = 0; i < num_points; i++) // for each point in the array
{ {
for (int j = 0; j < num_features; j++) // for each feature in the array for (int j = 0; j < num_features; j++) // for each feature in the array
{ {
fprintf(fp, "%.4g", X[i][j]); // print the feature value fprintf(fp, "%.4g", X[i][j]); // print the feature value
if (j < num_features - 1) // if not the last feature if (j < num_features - 1) // if not the last feature
fprintf(fp, ","); // suffix comma fprintf(fp, ","); // suffix comma
} }
if (i < num_points - 1) // if not the last row if (i < num_points - 1) // if not the last row
fprintf(fp, "\n"); // start a new line fprintf(fp, "\n"); // start a new line
} }
fclose(fp); fclose(fp);
return 0; return 0;
@ -96,12 +96,12 @@ int save_nd_data(const char *fname, double **X, int num_points,
*/ */
void get_min_1d(double const *X, int N, double *val, int *idx) void get_min_1d(double const *X, int N, double *val, int *idx)
{ {
val[0] = INFINITY; // initial min value val[0] = INFINITY; // initial min value
for (int i = 0; i < N; i++) // check each value for (int i = 0; i < N; i++) // check each value
{ {
if (X[i] < val[0]) // if a lower value is found if (X[i] < val[0]) // if a lower value is found
{ // save the value and its index { // save the value and its index
idx[0] = i; idx[0] = i;
val[0] = X[i]; val[0] = X[i];
} }
@ -212,8 +212,8 @@ void kohonen_som_tracer(double **X, double *const *W, int num_samples,
void test_circle(double *const *data, int N) void test_circle(double *const *data, int N)
{ {
const double R = 0.75, dr = 0.3; const double R = 0.75, dr = 0.3;
double a_t = 0., b_t = 2.f * M_PI; // theta random between 0 and 2*pi double a_t = 0., b_t = 2.f * M_PI; // theta random between 0 and 2*pi
double a_r = R - dr, b_r = R + dr; // radius random between R-dr and R+dr double a_r = R - dr, b_r = R + dr; // radius random between R-dr and R+dr
int i; int i;
#ifdef _OPENMP #ifdef _OPENMP
@ -221,9 +221,9 @@ void test_circle(double *const *data, int N)
#endif #endif
for (i = 0; i < N; i++) for (i = 0; i < N; i++)
{ {
double r = _random(a_r, b_r); // random radius double r = _random(a_r, b_r); // random radius
double theta = _random(a_t, b_t); // random theta double theta = _random(a_t, b_t); // random theta
data[i][0] = r * cos(theta); // convert from polar to cartesian data[i][0] = r * cos(theta); // convert from polar to cartesian
data[i][1] = r * sin(theta); data[i][1] = r * sin(theta);
} }
} }
@ -245,7 +245,7 @@ void test_circle(double *const *data, int N)
* "w12.csv" title "w2" * "w12.csv" title "w2"
* ``` * ```
* ![Sample execution * ![Sample execution
* output](https://raw.githubusercontent.com/kvedala/C/docs/images/machine_learning/kohonen/test1.svg) * output](https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/machine_learning/kohonen/test1.svg)
*/ */
void test1() void test1()
{ {
@ -259,28 +259,28 @@ void test1()
// number of clusters nodes * 2 // number of clusters nodes * 2
double **W = (double **)malloc(num_out * sizeof(double *)); double **W = (double **)malloc(num_out * sizeof(double *));
for (int i = 0; i < max(num_out, N); i++) // loop till max(N, num_out) for (int i = 0; i < max(num_out, N); i++) // loop till max(N, num_out)
{ {
if (i < N) // only add new arrays if i < N if (i < N) // only add new arrays if i < N
X[i] = (double *)malloc(features * sizeof(double)); X[i] = (double *)malloc(features * sizeof(double));
if (i < num_out) // only add new arrays if i < num_out if (i < num_out) // only add new arrays if i < num_out
{ {
W[i] = (double *)malloc(features * sizeof(double)); W[i] = (double *)malloc(features * sizeof(double));
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp for #pragma omp for
#endif #endif
// preallocate with random initial weights // preallocate with random initial weights
for (j = 0; j < features; j++) for (j = 0; j < features; j++) W[i][j] = _random(-1, 1);
W[i][j] = _random(-1, 1);
} }
} }
test_circle(X, N); // create test data around circumference of a circle test_circle(X, N); // create test data around circumference of a circle
save_nd_data("test1.csv", X, N, features); // save test data points save_nd_data("test1.csv", X, N, features); // save test data points
save_nd_data("w11.csv", W, num_out, save_nd_data("w11.csv", W, num_out,
features); // save initial random weights features); // save initial random weights
kohonen_som_tracer(X, W, N, features, num_out, 0.1); // train the SOM kohonen_som_tracer(X, W, N, features, num_out, 0.1); // train the SOM
save_nd_data("w12.csv", W, num_out, features); // save the resultant weights save_nd_data("w12.csv", W, num_out,
features); // save the resultant weights
for (int i = 0; i < max(num_out, N); i++) for (int i = 0; i < max(num_out, N); i++)
{ {
@ -315,10 +315,10 @@ void test_lamniscate(double *const *data, int N)
#endif #endif
for (i = 0; i < N; i++) for (i = 0; i < N; i++)
{ {
double dx = _random(-dr, dr); // random change in x double dx = _random(-dr, dr); // random change in x
double dy = _random(-dr, dr); // random change in y double dy = _random(-dr, dr); // random change in y
double theta = _random(0, M_PI); // random theta double theta = _random(0, M_PI); // random theta
data[i][0] = dx + cos(theta); // convert from polar to cartesian data[i][0] = dx + cos(theta); // convert from polar to cartesian
data[i][1] = dy + sin(2. * theta) / 2.f; data[i][1] = dy + sin(2. * theta) / 2.f;
} }
} }
@ -342,7 +342,7 @@ void test_lamniscate(double *const *data, int N)
* "w22.csv" title "w2" * "w22.csv" title "w2"
* ``` * ```
* ![Sample execution * ![Sample execution
* output](https://raw.githubusercontent.com/kvedala/C/docs/images/machine_learning/kohonen/test2.svg) * output](https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/machine_learning/kohonen/test2.svg)
*/ */
void test2() void test2()
{ {
@ -353,9 +353,9 @@ void test2()
double **W = (double **)malloc(num_out * sizeof(double *)); double **W = (double **)malloc(num_out * sizeof(double *));
for (int i = 0; i < max(num_out, N); i++) for (int i = 0; i < max(num_out, N); i++)
{ {
if (i < N) // only add new arrays if i < N if (i < N) // only add new arrays if i < N
X[i] = (double *)malloc(features * sizeof(double)); X[i] = (double *)malloc(features * sizeof(double));
if (i < num_out) // only add new arrays if i < num_out if (i < num_out) // only add new arrays if i < num_out
{ {
W[i] = (double *)malloc(features * sizeof(double)); W[i] = (double *)malloc(features * sizeof(double));
@ -363,17 +363,17 @@ void test2()
#pragma omp for #pragma omp for
#endif #endif
// preallocate with random initial weights // preallocate with random initial weights
for (j = 0; j < features; j++) for (j = 0; j < features; j++) W[i][j] = _random(-1, 1);
W[i][j] = _random(-1, 1);
} }
} }
test_lamniscate(X, N); // create test data around the lamniscate test_lamniscate(X, N); // create test data around the lamniscate
save_nd_data("test2.csv", X, N, features); // save test data points save_nd_data("test2.csv", X, N, features); // save test data points
save_nd_data("w21.csv", W, num_out, save_nd_data("w21.csv", W, num_out,
features); // save initial random weights features); // save initial random weights
kohonen_som_tracer(X, W, N, features, num_out, 0.01); // train the SOM kohonen_som_tracer(X, W, N, features, num_out, 0.01); // train the SOM
save_nd_data("w22.csv", W, num_out, features); // save the resultant weights save_nd_data("w22.csv", W, num_out,
features); // save the resultant weights
for (int i = 0; i < max(num_out, N); i++) for (int i = 0; i < max(num_out, N); i++)
{ {
@ -398,15 +398,15 @@ void test2()
*/ */
void test_3d_classes(double *const *data, int N) void test_3d_classes(double *const *data, int N)
{ {
const double R = 0.1; // radius of cluster const double R = 0.1; // radius of cluster
int i; int i;
const int num_classes = 4; const int num_classes = 4;
const double centres[][3] = { const double centres[][3] = {
// centres of each class cluster // centres of each class cluster
{.5, .5, .5}, // centre of class 1 {.5, .5, .5}, // centre of class 1
{.5, -.5, -.5}, // centre of class 2 {.5, -.5, -.5}, // centre of class 2
{-.5, .5, .5}, // centre of class 3 {-.5, .5, .5}, // centre of class 3
{-.5, -.5 - .5} // centre of class 4 {-.5, -.5 - .5} // centre of class 4
}; };
#ifdef _OPENMP #ifdef _OPENMP
@ -414,7 +414,8 @@ void test_3d_classes(double *const *data, int N)
#endif #endif
for (i = 0; i < N; i++) for (i = 0; i < N; i++)
{ {
int class = rand() % num_classes; // select a random class for the point int class =
rand() % num_classes; // select a random class for the point
// create random coordinates (x,y,z) around the centre of the class // create random coordinates (x,y,z) around the centre of the class
data[i][0] = _random(centres[class][0] - R, centres[class][0] + R); data[i][0] = _random(centres[class][0] - R, centres[class][0] + R);
@ -445,7 +446,7 @@ void test_3d_classes(double *const *data, int N)
* "w32.csv" title "w2" * "w32.csv" title "w2"
* ``` * ```
* ![Sample execution * ![Sample execution
* output](https://raw.githubusercontent.com/kvedala/C/docs/images/machine_learning/kohonen/test3.svg) * output](https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/machine_learning/kohonen/test3.svg)
*/ */
void test3() void test3()
{ {
@ -456,9 +457,9 @@ void test3()
double **W = (double **)malloc(num_out * sizeof(double *)); double **W = (double **)malloc(num_out * sizeof(double *));
for (int i = 0; i < max(num_out, N); i++) for (int i = 0; i < max(num_out, N); i++)
{ {
if (i < N) // only add new arrays if i < N if (i < N) // only add new arrays if i < N
X[i] = (double *)malloc(features * sizeof(double)); X[i] = (double *)malloc(features * sizeof(double));
if (i < num_out) // only add new arrays if i < num_out if (i < num_out) // only add new arrays if i < num_out
{ {
W[i] = (double *)malloc(features * sizeof(double)); W[i] = (double *)malloc(features * sizeof(double));
@ -466,17 +467,17 @@ void test3()
#pragma omp for #pragma omp for
#endif #endif
// preallocate with random initial weights // preallocate with random initial weights
for (j = 0; j < features; j++) for (j = 0; j < features; j++) W[i][j] = _random(-1, 1);
W[i][j] = _random(-1, 1);
} }
} }
test_3d_classes(X, N); // create test data around the lamniscate test_3d_classes(X, N); // create test data around the lamniscate
save_nd_data("test3.csv", X, N, features); // save test data points save_nd_data("test3.csv", X, N, features); // save test data points
save_nd_data("w31.csv", W, num_out, save_nd_data("w31.csv", W, num_out,
features); // save initial random weights features); // save initial random weights
kohonen_som_tracer(X, W, N, features, num_out, 0.01); // train the SOM kohonen_som_tracer(X, W, N, features, num_out, 0.01); // train the SOM
save_nd_data("w32.csv", W, num_out, features); // save the resultant weights save_nd_data("w32.csv", W, num_out,
features); // save the resultant weights
for (int i = 0; i < max(num_out, N); i++) for (int i = 0; i < max(num_out, N); i++)
{ {
@ -524,7 +525,8 @@ int main(int argc, char **argv)
end_clk = clock(); end_clk = clock();
printf("Test 3 completed in %.4g sec\n", printf("Test 3 completed in %.4g sec\n",
get_clock_diff(start_clk, end_clk)); get_clock_diff(start_clk, end_clk));
printf("(Note: Calculated times include: creating test sets, training " printf(
"model and writing files to disk.)\n\n"); "(Note: Calculated times include: creating test sets, training "
"model and writing files to disk.)\n\n");
return 0; return 0;
} }

18
numerical_methods/durand_kerner_roots.c
View File

@ -21,10 +21,10 @@
* Sample implementation results to compute approximate roots of the equation * Sample implementation results to compute approximate roots of the equation
* \f$x^4-1=0\f$:\n * \f$x^4-1=0\f$:\n
* <img * <img
* src="https://raw.githubusercontent.com/kvedala/C/docs/images/numerical_methods/durand_kerner_error.svg" * src="https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/numerical_methods/durand_kerner_error.svg"
* width="400" alt="Error evolution during root approximations computed every * width="400" alt="Error evolution during root approximations computed every
* iteration."/> <img * iteration."/> <img
* src="https://raw.githubusercontent.com/kvedala/C/docs/images/numerical_methods/durand_kerner_roots.svg" * src="https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/numerical_methods/durand_kerner_roots.svg"
* width="400" alt="Roots evolution - shows the initial approximation of the * width="400" alt="Roots evolution - shows the initial approximation of the
* roots and their convergence to a final approximation along with the iterative * roots and their convergence to a final approximation along with the iterative
* approximations" /> * approximations" />
@ -53,8 +53,7 @@ long double complex poly_function(double *coeffs, unsigned int degree,
long double complex out = 0.; long double complex out = 0.;
unsigned int n; unsigned int n;
for (n = 0; n < degree; n++) for (n = 0; n < degree; n++) out += coeffs[n] * cpow(x, degree - n - 1);
out += coeffs[n] * cpow(x, degree - n - 1);
return out; return out;
} }
@ -102,8 +101,9 @@ int main(int argc, char **argv)
if (argc < 2) if (argc < 2)
{ {
printf("Please pass the coefficients of the polynomial as commandline " printf(
"arguments.\n"); "Please pass the coefficients of the polynomial as commandline "
"arguments.\n");
return 0; return 0;
} }
@ -224,8 +224,7 @@ int main(int argc, char **argv)
if (iter % 500 == 0) if (iter % 500 == 0)
{ {
printf("Iter: %lu\t", iter); printf("Iter: %lu\t", iter);
for (n = 0; n < degree - 1; n++) for (n = 0; n < degree - 1; n++) printf("\t%s", complex_str(s0[n]));
printf("\t%s", complex_str(s0[n]));
printf("\t\tabsolute average change: %.4g\n", tol_condition); printf("\t\tabsolute average change: %.4g\n", tol_condition);
} }
@ -241,8 +240,7 @@ end:
#endif #endif
printf("\nIterations: %lu\n", iter); printf("\nIterations: %lu\n", iter);
for (n = 0; n < degree - 1; n++) for (n = 0; n < degree - 1; n++) printf("\t%s\n", complex_str(s0[n]));
printf("\t%s\n", complex_str(s0[n]));
printf("absolute average change: %.4g\n", tol_condition); printf("absolute average change: %.4g\n", tol_condition);
printf("Time taken: %.4g sec\n", printf("Time taken: %.4g sec\n",
(end_time - start_time) / (double)CLOCKS_PER_SEC); (end_time - start_time) / (double)CLOCKS_PER_SEC);

23
numerical_methods/ode_forward_euler.c
View File

@ -22,7 +22,7 @@
* The computation results are stored to a text file `forward_euler.csv` and the * The computation results are stored to a text file `forward_euler.csv` and the
* exact soltuion results in `exact.csv` for comparison. * exact soltuion results in `exact.csv` for comparison.
* <img * <img
* src="https://raw.githubusercontent.com/kvedala/C/docs/images/numerical_methods/ode_forward_euler.svg" * src="https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/numerical_methods/ode_forward_euler.svg"
* alt="Implementation solution"/> * alt="Implementation solution"/>
* *
* To implement [Van der Pol * To implement [Van der Pol
@ -54,9 +54,9 @@
*/ */
void problem(const double *x, double *y, double *dy) void problem(const double *x, double *y, double *dy)
{ {
const double omega = 1.F; // some const for the problem const double omega = 1.F; // some const for the problem
dy[0] = y[1]; // x dot dy[0] = y[1]; // x dot
dy[1] = -omega * omega * y[0]; // y dot dy[1] = -omega * omega * y[0]; // y dot
} }
/** /**
@ -83,8 +83,7 @@ void forward_euler_step(const double dx, const double *x, double *y, double *dy)
{ {
int o; int o;
problem(x, y, dy); problem(x, y, dy);
for (o = 0; o < order; o++) for (o = 0; o < order; o++) y[o] += dx * dy[o];
y[o] += dx * dy[o];
} }
/** /**
@ -116,13 +115,13 @@ double forward_euler(double dx, double x0, double x_max, double *y,
/* start integration */ /* start integration */
clock_t t1 = clock(); clock_t t1 = clock();
double x = x0; double x = x0;
do // iterate for each step of independent variable do // iterate for each step of independent variable
{ {
if (save_to_file && fp) if (save_to_file && fp)
fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file
forward_euler_step(dx, &x, y, dy); // perform integration forward_euler_step(dx, &x, y, dy); // perform integration
x += dx; // update step x += dx; // update step
} while (x <= x_max); // till upper limit of independent variable } while (x <= x_max); // till upper limit of independent variable
/* end of integration */ /* end of integration */
clock_t t2 = clock(); clock_t t2 = clock();
@ -169,7 +168,7 @@ int main(int argc, char *argv[])
do do
{ {
fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file
exact_solution(&x, y); exact_solution(&x, y);
x += step_size; x += step_size;
} while (x <= X_MAX); } while (x <= X_MAX);

26
numerical_methods/ode_midpoint_euler.c
View File

@ -21,7 +21,7 @@
* \f} * \f}
* The computation results are stored to a text file `midpoint_euler.csv` and * The computation results are stored to a text file `midpoint_euler.csv` and
* the exact soltuion results in `exact.csv` for comparison. <img * the exact soltuion results in `exact.csv` for comparison. <img
* src="https://raw.githubusercontent.com/kvedala/C/docs/images/numerical_methods/ode_midpoint_euler.svg" * src="https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/numerical_methods/ode_midpoint_euler.svg"
* alt="Implementation solution"/> * alt="Implementation solution"/>
* *
* To implement [Van der Pol * To implement [Van der Pol
@ -53,9 +53,9 @@
*/ */
void problem(const double *x, double *y, double *dy) void problem(const double *x, double *y, double *dy)
{ {
const double omega = 1.F; // some const for the problem const double omega = 1.F; // some const for the problem
dy[0] = y[1]; // x dot dy[0] = y[1]; // x dot
dy[1] = -omega * omega * y[0]; // y dot dy[1] = -omega * omega * y[0]; // y dot
} }
/** /**
@ -86,13 +86,11 @@ void midpoint_euler_step(double dx, double *x, double *y, double *dy)
double tmp_x = (*x) + 0.5 * dx; double tmp_x = (*x) + 0.5 * dx;
double tmp_y[order]; double tmp_y[order];
int o; int o;
for (o = 0; o < order; o++) for (o = 0; o < order; o++) tmp_y[o] = y[o] + 0.5 * dx * dy[o];
tmp_y[o] = y[o] + 0.5 * dx * dy[o];
problem(&tmp_x, tmp_y, dy); problem(&tmp_x, tmp_y, dy);
for (o = 0; o < order; o++) for (o = 0; o < order; o++) y[o] += dx * dy[o];
y[o] += dx * dy[o];
} }
/** /**
@ -124,13 +122,13 @@ double midpoint_euler(double dx, double x0, double x_max, double *y,
/* start integration */ /* start integration */
clock_t t1 = clock(); clock_t t1 = clock();
double x = x0; double x = x0;
do // iterate for each step of independent variable do // iterate for each step of independent variable
{ {
if (save_to_file && fp) if (save_to_file && fp)
fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file
midpoint_euler_step(dx, &x, y, dy); // perform integration midpoint_euler_step(dx, &x, y, dy); // perform integration
x += dx; // update step x += dx; // update step
} while (x <= x_max); // till upper limit of independent variable } while (x <= x_max); // till upper limit of independent variable
/* end of integration */ /* end of integration */
clock_t t2 = clock(); clock_t t2 = clock();
@ -177,7 +175,7 @@ int main(int argc, char *argv[])
do do
{ {
fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file
exact_solution(&x, y); exact_solution(&x, y);
x += step_size; x += step_size;
} while (x <= X_MAX); } while (x <= X_MAX);

30
numerical_methods/ode_semi_implicit_euler.c
View File

@ -21,7 +21,7 @@
* \f} * \f}
* The computation results are stored to a text file `semi_implicit_euler.csv` * The computation results are stored to a text file `semi_implicit_euler.csv`
* and the exact soltuion results in `exact.csv` for comparison. <img * and the exact soltuion results in `exact.csv` for comparison. <img
* src="https://raw.githubusercontent.com/kvedala/C/docs/images/numerical_methods/ode_semi_implicit_euler.svg" * src="https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/numerical_methods/ode_semi_implicit_euler.svg"
* alt="Implementation solution"/> * alt="Implementation solution"/>
* *
* To implement [Van der Pol * To implement [Van der Pol
@ -33,7 +33,7 @@
* dy[1] = mu * (1.f - y[0] * y[0]) * y[1] - y[0]; * dy[1] = mu * (1.f - y[0] * y[0]) * y[1] - y[0];
* ``` * ```
* <a href="https://en.wikipedia.org/wiki/Van_der_Pol_oscillator"><img * <a href="https://en.wikipedia.org/wiki/Van_der_Pol_oscillator"><img
* src="https://raw.githubusercontent.com/kvedala/C/docs/images/numerical_methods/van_der_pol_implicit_euler.svg" * src="https://raw.githubusercontent.com/TheAlgorithms/C/docs/images/numerical_methods/van_der_pol_implicit_euler.svg"
* alt="Van der Pol Oscillator solution"/></a> * alt="Van der Pol Oscillator solution"/></a>
* *
* \see ode_forward_euler.c, ode_midpoint_euler.c * \see ode_forward_euler.c, ode_midpoint_euler.c
@ -57,9 +57,9 @@
*/ */
void problem(const double *x, double *y, double *dy) void problem(const double *x, double *y, double *dy)
{ {
const double omega = 1.F; // some const for the problem const double omega = 1.F; // some const for the problem
dy[0] = y[1]; // x dot dy[0] = y[1]; // x dot
dy[1] = -omega * omega * y[0]; // y dot dy[1] = -omega * omega * y[0]; // y dot
} }
/** /**
@ -86,13 +86,13 @@ void semi_implicit_euler_step(double dx, double *x, double *y, double *dy)
{ {
int o; int o;
problem(x, y, dy); // update dy once problem(x, y, dy); // update dy once
y[0] += dx * dy[0]; // update y0 y[0] += dx * dy[0]; // update y0
problem(x, y, dy); // update dy once more problem(x, y, dy); // update dy once more
for (o = 1; o < order; o++) for (o = 1; o < order; o++)
y[o] += dx * dy[o]; // update remaining using new dy y[o] += dx * dy[o]; // update remaining using new dy
*x += dx; *x += dx;
} }
@ -125,13 +125,13 @@ double semi_implicit_euler(double dx, double x0, double x_max, double *y,
/* start integration */ /* start integration */
clock_t t1 = clock(); clock_t t1 = clock();
double x = x0; double x = x0;
do // iterate for each step of independent variable do // iterate for each step of independent variable
{ {
if (save_to_file && fp) if (save_to_file && fp)
fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file
semi_implicit_euler_step(dx, &x, y, dy); // perform integration semi_implicit_euler_step(dx, &x, y, dy); // perform integration
x += dx; // update step x += dx; // update step
} while (x <= x_max); // till upper limit of independent variable } while (x <= x_max); // till upper limit of independent variable
/* end of integration */ /* end of integration */
clock_t t2 = clock(); clock_t t2 = clock();
@ -178,7 +178,7 @@ int main(int argc, char *argv[])
do do
{ {
fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file fprintf(fp, "%.4g,%.4g,%.4g\n", x, y[0], y[1]); // write to file
exact_solution(&x, y); exact_solution(&x, y);
x += step_size; x += step_size;
} while (x <= X_MAX); } while (x <= X_MAX);