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

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

int main ( int argc, char **argv );
double square_exact_sol ( double x, double y );
void square_phi ( struct Neural_Net_PDE *nn, double x, double y );
double square_pde ( double x, double y, double u, double u_x, double u_y,
  double u_xx, double u_xy, double u_yy );

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

int main ( int argc, char **argv )

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

    ppc_neural_pde_square() tests ppc_neural_pde() on the square region.

  Discussion:

    Seek a solution in the square [0,1]x[0,1] of the boundary value problem:

      d2 u/dx2 + d2 u/dy2 = f(x,y)

    where

      u(x,y) = 0 on the boundary

    and

      f(x,y) = 32 ( x^2 + y^2 - x - y )

    for which the exact solution is

      u(x,y) = 16 * x * ( x - 1 ) * y * ( y - 1 )

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    16 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
*/
{
  char *endptr;
  double **training_points;

  timestamp ( );
  printf ( "\n" );
  printf ( "ppc_neural_pde_square():\n" );
  printf ( "  ppc_neural_pde() models a partial differential equation\n" );
  printf ( "  in the form of a boundary value problem, solved\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;
  }
/*
  Create space for training points.
*/
  make_matrix ( training_points, nu, 2 );
/*
  Define the neural network.
*/
  struct Neural_Net_PDE nn =
  {
    .PDE = square_pde,
    .phi_func = square_phi,
    .bb_xrange = { 0.0, 1.0 },
    .bb_yrange = { 0.0, 1.0 },
    .q = q,
    .nu = nu,
    .training_points = training_points,
    .exact_sol = square_exact_sol,
    .geomview_out = "square.gv",
    .maple_out = "square.m",
    .matlab_out = "square.mpl",
  };
/*
  Generate nu training points.
  Do this by randomly picking points in the bounding box,
  and then checking that 0 < phi(x,y).
*/
  int i = 0;
  while ( i < nu )
  {
    double r = ( double ) rand() / RAND_MAX;
    double s = ( double ) rand() / RAND_MAX;
    double x = ( 1.0 - r ) * nn.bb_xrange[0] + r * nn.bb_xrange[1];
    double y = ( 1.0 - s ) * nn.bb_yrange[0] + s * nn.bb_yrange[1];
    nn.phi_func ( &nn, x, y );
    if ( 0 < nn.phi[0][0] )
    {
      nn.training_points[i][0] = x;
      nn.training_points[i][1] = y;
      i++;
    }
  }
/*
  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, 50 ) );
  }

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

  Neural_Net_end ( &nn );

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

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

double square_exact_sol ( double x, double y )

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

    square_exact_sol() evaluates the exact solution.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    16 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, y: the evaluation point.

  Output:

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

  value = 32.0 * x * ( x - 1.0 ) * y * ( y - 1.0 );

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

void square_phi ( struct Neural_Net_PDE *nn, double x, double y )

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

    square_phi() indicates whether a point is inside the region.

  Discussion:

    We want to choose a function phi(x,y) which distinguishes
    between points (x,y) in the bounding box.
    phi(x,y) < 0 means these points are not part of the region.
    Since our region is the square (0,1)x(0,1), all the points in
    the bounding box are in the region, so we just need to ensure that
    phi(x,y) is positive there.
    For this simple problem, even phi(x,y)=1 would work.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    16 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_PDE *nn: a structure defining a neural network.

    double x, y: the evaluation point.

  Output:

    struct Neural_Net_PDE *nn: the structure has been updated with
    phi[0][0], phi[1][0], phi[0][1], phi[2][0], phi[1][1], phi[0][2].
*/
{
  nn->phi[0][0] = x * x * y * y - x * x * y - x * y * y + x * y;

  nn->phi[1][0] = 2.0 * x * y * y - 2.0 * x * y - y * y + y;
  nn->phi[0][1] = 2.0 * x * x * y - x * x - 2.0 * x * y + x;

  nn->phi[2][0] = 2.0 * y * y - 2.0 * y;
  nn->phi[1][1] = 4.0 * x * y - 2.0 * x - 2.0 * y;
  nn->phi[0][2] = 2.0 * x * x - 2.0 * x;

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

double square_pde ( double x, double y, double u, double u_x, double u_y,
  double u_xx, double u_xy, double u_yy )

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

    square_pde() evaluates the PDE residual for the square problem.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    16 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, y: the evaluation point.

  Output:

    double value: the value of the PDE residual, uxx + uyy - f.
*/
{
  double f;
  double value;

  f = 32.0 * ( x * y + y * y - x - y );

  value = u_xx + u_yy - f;

  return value;
}