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

# include "ppc_array.h"
# include "ppc_nelder_mead.h"
# include "ppc_neural_ode.h"
# include "ppc_xmalloc.h"

int main ( int argc, char **argv );
double exact_sol ( double x );
double my_ode ( double x, double u, double ux, double uxx );
void show_usage ( char *progname );
void timestamp ( );

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

int main ( int argc, char **argv )

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

    ppc_neural_ode_test() tests ppc_neural_ode().

  Discussion:

    Seek a solution of the BVP

      u'' + 1 / ( 1 + u^2 ) = f
      u(0) = u(pi) = 0

    where

      f(x) = - sin(x) - 4 sin(2x) + 1 / ( 1 + ( sin(x) + sin(2x) )^2 )

    for which the exact solution is

      u(x) = sin(x) + sin(2x)

  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
*/
{
  double a;
  double b;
  char *endptr;
  double PI = 4.0 * atan ( 1.0 );
  double *training_points;

  timestamp ( );
  printf ( "\n" );
  printf ( "ppc_neural_ode_test():\n" );
  printf ( "  ppc_neural_ode() models a boundary value problem\n" );
  printf ( "  using a neural network.\n" );

  if ( argc != 3 )
  {
    show_usage ( argv[0] );
    return EXIT_FAILURE;
  }
/*
  Number of units in the hidden layer.
*/
  int q = strtol ( argv[1], &endptr, 10 );
  if ( *endptr != '\0' || q < 1 )
  {
    show_usage ( argv[0] );
    return EXIT_FAILURE;
  }
/*
  Number of training points.
*/
  int nu = strtol ( argv[2], &endptr, 10 );
  if ( *endptr != '\0' || nu < 1 )
  {
    show_usage ( argv[0] );
    return EXIT_FAILURE;
  }
/*
  nu equally spaced points inside (a,b)
*/
  a = 0.0;
  b = PI;
  make_vector ( training_points, nu );
  for ( int i = 0; i < nu; i++ )
  {
    training_points[i] = a + ( b - a ) / ( nu + 1 ) * ( i + 1 );
  }
/*
  Define the neural network.
*/
  struct Neural_Net_ODE nn =
  {
    .a = a,
    .b = b,
    .q = q,
    .nu = nu,
    .training_points = training_points,
    .ODE = my_ode,
    .exact_sol = exact_sol,
  };
/*
  Initialize the neural net.
*/
  Neural_Net_init ( &nn );
/*
  Create a structure for Nelder-Mead
*/
  struct nelder_mead NM =
  {
    .f = Neural_Net_residual,
    .n = nn.nweights,
    .s = NULL,
    .x = nn.weights,
    .h = 0.1,
    .tol = 1e-5,
    .maxevals = 100000,
    .params = &nn,
  };

  printf ( "  weights before training:\n" );
  print_vector ( "%7.3f ", nn.weights, nn.nweights );
/*
  Call nelder-mead to minimize the error function.
*/
  int evalcount = nelder_mead ( &NM );
/*
  Report the results.
*/
  if ( NM.maxevals < evalcount )
  {
    printf ( "Nelder-Mead: No convergence after %d "
      "function evaluation\n", evalcount );
    return EXIT_FAILURE;
  }

  printf ( "Nelder-Mead: Converged after %d ""function evaluations\n", 
    evalcount );
  printf ( "Nelder-Mead: Neural network’s residual error = %g\n", NM.minval );

  printf ( "weights after training:\n");
  print_vector ( "%7.3f ", nn.weights, nn.nweights );

  if ( nn.exact_sol != NULL )
  {
    printf ( "Error versus the ODE’s exact solution = %g\n",
      Neural_Net_error_vs_exact ( &nn, 50 ) );
  }

  Neural_Net_plot_with_maple ( &nn, 50, "ppc_neural_ode_test.mpl" );
  Neural_Net_plot_with_matlab ( &nn, 50, "ppc_neural_ode_test" );

  Neural_Net_end ( &nn );

  free_vector ( training_points );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "ppc_neural_ode_test():\n" );
  printf ( "  Normal end of execution.\n" );
  timestamp ( );

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

double exact_sol ( double x )

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

    exact_sol() evaluates the exact solution.

  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:

    double x: the evaluation point.

  Output:

    double value: the exact solution of the BVP, evaluated at x.
*/
{
  double value;

  value = sin ( x ) + sin ( 2.0 * x );

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

double my_ode ( double x, double u, double ux, double uxx )

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

    my_ode() evaluates the functional defining the boundary value problem.

  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:

    double x: the evaluation point.

    double u, ux, uxx: the values of the solution and its first two
    derivatives, at x.

  Output:

    double value: the residual for the boundary value problem.
*/
{
  double f;
  double t;
  double value;

  t = sin ( x ) + sin ( 2.0 * x );
  f = - sin ( x ) - 4.0 * sin ( 2.0 * x ) + 1.0 / ( 1.0 + t * t );

  value = uxx + 1.0 / ( 1.0 + u * u ) - f;

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

void show_usage ( char *progname )

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

    show_usage() reminds the user how to call the program.

  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
*/
{
  printf ( "\n" );
  printf ( "Usage:\n" );
  printf ( "  ./%s q nu\n", progname );
  printf ( "  q:  number of units in the hidden layer (1 <= q)\n" );
  printf ( "  nu : number of training points (1 <= nu)\n" );

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

void timestamp ( )

/******************************************************************************/
/*
  Purpose:
 
    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 May 2021

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}