# include <math.h>
# include <stdio.h>
# include <stdlib.h>

# include "ppc_array.h"
# include "ppc_neural_ode.h"

/******************************************************************************/

void Neural_Net_end ( struct Neural_Net_ODE *nn )

/******************************************************************************/
/*
  Purpose:

    Neural_Net_end() frees the memory associated with a neural network.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    22 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct Neural_Net_ODE *nn: describes the neural network.
*/
{
  free_vector ( nn->weights );

  return;
}
/******************************************************************************/

double Neural_Net_error_vs_exact ( struct Neural_Net_ODE *nn, int n )

/******************************************************************************/
/*
  Purpose:

    Neural_Net_error_vs_exact() compares computed and exact ODE solutions.

  Discussion:

    Evaluate the trained neural network at n+1 points x[i] equally spaced 
    in [a,b], calculate each discrepancy
      | phi(x[i]) * N(x[i]) - uexact(x[i]) |
    and return the maximum.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    14 June 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct Neural_Net_ODE *nn: describes the neural network.

    int n: the interval [a,b] is evaluated at n+1 equally spaced points.

  Output:

    double value: the maximum difference between the exact solution
    and the neural network model.
*/
{
  double a;
  double b;
  int i;
  double uex;
  double unn;
  double value;
  double x;

  a = nn->a;
  b = nn->b;
  value = 0.0;
  for ( i = 0; i <= n; i++ )
  {
    x = ( ( n - i ) * a + i * b ) / n;

    uex = nn->exact_sol ( x );

    Neural_Net_phi ( nn, x );
    Neural_Net_eval ( nn, x );
    unn = nn->N[0] * nn->phi[0];

    value = fmax ( value, fabs ( uex - unn ) );
  }

  return value;
}
/******************************************************************************/

void Neural_Net_eval ( struct Neural_Net_ODE *nn, double x )

/******************************************************************************/
/*
  Purpose:

    Neural_Net_eval() evaluates a neural network.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    22 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct Neural_Net_ODE *nn: describes the neural network.

    double x: the evaluate point.

  Output:

    struct Neural_Net_ODE *nn: nn->N[0], nn->N[1], nn->N[1] have been set.
*/
{
  int q = nn->q;
  double *u = nn->weights;
  double *v = nn->weights + q;
  double *w = nn->weights + 2 * q;

  for ( int j = 0; j <= 2; j++ )
  {
    nn->N[j] = 0.0;
  }

  for ( int i = 0; i < q; i++ )
  {
    double z = u[i] + w[i] * x;
    sigmoid ( z, nn->sigma );
    for ( int j = 0; j <= 2; j++ )
    {
      nn->N[j] = nn->N[j] + v[i] * pow ( w[i], j ) * nn->sigma[j];
    }
  }
  return;
}
/******************************************************************************/

void Neural_Net_init ( struct Neural_Net_ODE *nn )

/******************************************************************************/
/*
  Purpose:

    Neural_Net_init() allocates and initializes a neural network.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    22 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct Neural_Net_ODE *nn: describes the neural network.

  Ooutput:

    struct Neural_Net_ODE *nn: the structure has been initialized.
*/
{
  nn->nweights = 3 * nn->q;
  make_vector ( nn->weights, nn->nweights );
/*
  Initialize the weights with random values in [-0.5,+0.5].
*/
  for ( int i = 0; i < nn->nweights; i++ )
  {
    nn->weights[i] = ( double ) rand ( ) / RAND_MAX - 0.5;
  }

  return;
}
/******************************************************************************/

void Neural_Net_phi ( struct Neural_Net_ODE *nn, double x )

/******************************************************************************/
/*
  Purpose:

    Neural_Net_phi() enforces the ODE boundary conditions.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    13 June 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct Neural_Net_ODE *nn; a structure defining a neural net.

    double x: an evaluation point.

  Output:

    struct Neural_Net_ODE *nn: the components phi[0], phi[1] and phi[2]
    have been evaluated at x.
*/
{
  double a;
  double b;

  a = nn->a;
  b = nn->b;

  nn->phi[0] = ( b - x ) * ( x - a );
  nn->phi[1] = - 2.0 * x + a + b;
  nn->phi[2] = - 2.0;

  return;
}
/******************************************************************************/

void Neural_Net_plot_with_maple ( struct Neural_Net_ODE *nn, int n,
  char *outfile )

/******************************************************************************/
/* 
  Purpose:

    Neural_Net_plot_with_maple() writes a Maple script to plot the solution.

  Discussion:

    The solution is sampled at n+1 equally
    spaced points x[0] = a, x[1], x[2], ..., x[n] on the interval [a,b].
    The distance between a pair consecutive points is (b-a)/n.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    22 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct Neural_Net_ODE *nn: describes the neural network.

    int n: the number of sample points.

    char *outfile: the name of the file to be created.
*/
{
  double a = nn->a;
  double b = nn->b;
  FILE *fp;

  fp = fopen ( outfile, "w" );

  if ( fp == NULL )
  {
    fprintf ( stderr, "unable to open file `%s' for writing\n", outfile );
    return;
  }

  fprintf ( fp, "plot([" );
  for ( int i = 0; i <= n; i++ )
  {
    double x = a + ( b - a ) / n * i;
    Neural_Net_phi( nn, x );
    Neural_Net_eval ( nn, x );
    double y = nn->N[0] * nn->phi[0];
    fprintf ( fp, "[%g,%g],", x, y );
  }

  fprintf ( fp, "NULL], title=\"hidden units = %d, training points = %d",
    nn->q, nn->nu );

  if ( nn->exact_sol != NULL )
  {
    fprintf ( fp, "\\nerror vs exact = %g", 
      Neural_Net_error_vs_exact ( nn, n ) );
  }
  fprintf ( fp, "\", color=\"Green\", labels=[x, u(x)], "
    "size=[600,300], font=[Times,bold,14]);\n" );
  fclose ( fp );

  return;
}
/******************************************************************************/

void Neural_Net_plot_with_matlab ( struct Neural_Net_ODE *nn, int n,
  char *header )

/******************************************************************************/
/* 
  Purpose:

    Neural_Net_plot_with_matlab writes a Matlab script to plot the solution.

  Discussion:

    The solution is sampled at n+1 equally
    spaced points x[0] = a, x[1], x[2], ..., x[n] on the interval [a,b].
    The distance between a pair consecutive points is (b-a)/n.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 June 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct Neural_Net_ODE *nn: describes the neural network.

    int n: the number of sample points.

    char *header: a string to be used to create the names of the 
    script file and the saved PNG file.
*/
{
#define NN	"NN_plot_with_matlab"

  double a = nn->a;
  double b = nn->b;
  char filename[255];
  FILE *fp;

  sprintf ( filename, "%s.m", header );

  fp = fopen ( filename, "w" );

  if ( fp == NULL )
  {
    fprintf ( stderr, "unable to open file `%s' for writing\n", filename );
    return;
  }

  fprintf ( fp, "close all\n" );
  fprintf ( fp, "%s = [ ", NN );
  for ( int i = 0; i <= n; i++ )
  {
    double x = a + ( b - a ) / n * i;
    Neural_Net_phi ( nn, x );
    Neural_Net_eval ( nn, x );
    double y = nn->N[0] * nn->phi[0];
    fprintf ( fp, "%g, %g; ", x, y );
  }
  fprintf ( fp, "];\n" );
  fprintf ( fp, "plot(%s(:,1),  %s(:,2), 'r', 'linewidth', 2 );\n", NN, NN );
  fprintf ( fp, "title({'hidden units = %d, training points = %d'", 
    nn->q, nn->nu );

  if ( nn->exact_sol != NULL )
  {
    fprintf ( fp, ", 'error vs exact = %g'", 
      Neural_Net_error_vs_exact ( nn, n ) );
  }
  fprintf ( fp, "});\n" );
  fprintf ( fp, "grid on;\n" );

  sprintf ( filename, "%s.png", header );
  fprintf ( fp, "print ( '-dpng', '%s' )\n", filename );

  fclose ( fp );

  return;

#undef NN
}
/******************************************************************************/

double Neural_Net_residual ( double *weights, int nweights, void *params )

/******************************************************************************/
/*
  Purpose:

    Neural_Net_residual() evaluates the neural network residual.

  Discussion:

    The interface is designed to correspond to the function minimizer
    ppc_nelder_mead().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    22 May 2024

  Author:

    Original C version by Rouben Rostamian.
    This version by John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    double *weights: the weights.

    int nweights: the number of weights.

    void *params: associated parameters.

  Output:

    Neural_Net_residual: residual associated with this set of weights.
*/
{
  struct Neural_Net_ODE *nn = params;
/*
  nu is the number of training points.
*/
  int nu = nn->nu;
  double sum = 0.0;

  for ( int i = 0; i < nn->nweights; i++ )
  {
    nn->weights[i] = weights[i];
  }

  for ( int i = 0; i < nu; i++ )
  {
    double x = nn->training_points[i];
    double r = residual_at_x ( nn, x );
    sum = sum + r * r;
  }

  return sum;
}
/******************************************************************************/

double residual_at_x ( struct Neural_Net_ODE *nn, double x )

/******************************************************************************/
/*
  Purpose:

    residual_at_x() evaluates the ODE residual at the point X.

  Discussion:

    R(x) = F ( x, phi*N, phi'N+phi*N', phi''N+2*phi'N'+phiN'' )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 June 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5

  Input:

    struct Neural_Net_ODE *nn: describes the neural network.

    double x: the evaluation point.

  Output:

    double value: the ODE residual at x.
*/
{
  double u;
  double ux;
  double uxx;
  double value;

  Neural_Net_phi ( nn, x );
  Neural_Net_eval ( nn, x );

  u = nn->phi[0] * nn->N[0];

  ux = nn->phi[1] * nn->N[0] 
     + nn->phi[0] * nn->N[1];

  uxx =         nn->phi[2] * nn->N[0]
        + 2.0 * nn->phi[1] * nn->N[1]
        +       nn->phi[0] * nn->N[2];

  value = nn->ODE ( x, u, ux, uxx );

  return value;
}
/******************************************************************************/

void sigmoid ( double x, double *sigma )

/******************************************************************************/
/*
  Purpose:

    sigmoid() evaluates the sigmoidal function and its first two derivatives.

  Discussion:

    sigmoid(x) = 1 / ( 1 + exp(-x) )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    22 May 2024

  Author:

    John Burkardt

  Reference:

    Rouben Rostamian,
    Programming Projects in C 
    for Students of Engineering, Science, and Mathematics,
    SIAM, 2014,
    ISBN: 978-1-611973-49-5
*/
{
  double value;

  value = 1.0 / ( 1.0 + exp ( - x ) );

  sigma[0] = value;
  sigma[1] = value * ( 1.0 - value );
  sigma[2] = value * ( 1.0 - value ) * ( 1.0 - 2.0 * value );

  return;
}