#! /usr/bin/env python3
#
def paraheat_gaussian_parameters ( data_filename ):

#*****************************************************************************80
#
## paraheat_gaussian_parameters uses keras to solve a multivariate regression case.
#
#  Discussion:
#
#    A data file contains data_num records of an experiment.
#
#    Each record has the form:
#      (seed,xc,yc,sc,vc,vs01...vs50)
#    where seed was used to generate the parameter vc, which controls the
#    strength of a Gaussian diffusivity, and vs01 through vs50 are values 
#    of the resulting finite element solution to a heat equation.
#    For this example, all four quantites xc, yc, sc and vc were varied.
#
#    The task is to use the values of vs (sample solution values at sensor 
#    locations) to estimate xc, yc, sc and vc.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    14 May 2019
#
#  Author:
#
#    John Burkardt.
#
  import keras
  import matplotlib.pyplot as plt
  import numpy as np
  import platform
  
  print ( '' )
  print ( 'paraheat_gaussian_parameters:' )
  print ( '  python version: %s' % ( platform.python_version ( ) ) )
  print ( '  keras version: %s' % ( keras.__version__ ) )
  print ( '  Neural network to solve a multivariate regression problem.' )
  print ( '  Estimate the parameters xc, yc, sc, vc used in a Gaussian diffusivity' )
  print ( '  given vs, 50 samples of the resulting heat distribution' )
  print ( '  Data of many records is available.' )
  print ( '  The data is read from an external file.' )
  print ( '' )
#
#  Load the data.
#
  print ( '  Read data from %s' % ( data_filename ) )

  data = np.loadtxt ( data_filename )
  
  print ( '  Data contains %d records with %d features.' % ( data.shape[0], data.shape[1] ) )
#
#  Each data item has 55 entries:
#  Data column 1 is the random number seed; ignore it.
#  Data columns 2, 3, 4, 5 are xc, yc, sc, vc; we will try to predict this
#  Data columns 6-55 are sampled solution values.  Use these to predict.
#
  data_num = data.shape[0]
  train_num = int ( 0.95 * data_num )
#
#  There are data_num rows of data.
#  Use train_num for training, test_num for testing.
#
  train_data = data[0:train_num,5:55]
  train_targets = data[0:train_num,1:5]
  train_num = train_data.shape[0]
  feature_num = train_data.shape[1]
  target_num = train_targets.shape[1]
  print ( '  Training data uses %d records with %d features and %d targets.' % ( train_num, feature_num, target_num ) )

  test_data = data[train_num:data_num,5:55]
  test_targets = data[train_num:data_num,1:5]
  test_num = test_data.shape[0]
  print ( '  Test data uses %d records with %d features and %d targets.' % ( test_num, feature_num, target_num ) )
#
#  Exhibit the contents of one item of the training data:
#
  print ( '  train_data[0,0:10]:' )
  print ( train_data[0,0:10] )
  print ( '  train_targets[0,0:4]:' )
  print ( train_targets[0,0:4] )

  print ( '  test_data[0,0:10]:' )
  print ( test_data[0,0:10] )
  print ( '  test_targets[0,0:4]:' )
  print ( test_targets[0,0:4] )
#
#  Build the model.
#  Options include how many nodes and how many layers to use.
#
  from keras import models
  from keras import layers

  node_num = 200
  print ( '  Using %d nodes per layer' % ( node_num ) )

  model = models.Sequential()
  model.add ( layers.Dense ( node_num, activation = 'relu', input_shape = ( feature_num, ) ) )
  model.add ( layers.Dense ( node_num, activation = 'relu' ) )
  model.add ( layers.Dense ( node_num, activation = 'relu' ) )
  model.add ( layers.Dense ( node_num, activation = 'relu' ) )
  model.add ( layers.Dense ( node_num, activation = 'relu' ) )
  model.add ( layers.Dense ( node_num, activation = 'relu' ) )
  model.add ( layers.Dense ( node_num, activation = 'relu' ) )
  model.add ( layers.Dense ( node_num, activation = 'relu' ) )
  model.add ( layers.Dense ( node_num, activation = 'relu' ) )
  model.add ( layers.Dense ( node_num, activation = 'relu' ) )
  model.add ( layers.Dense ( node_num, activation = 'relu' ) )
  model.add ( layers.Dense ( 4 ) )
#
#  Compile the model.
#
  model.compile ( \
    optimizer = 'adam', \
    loss = 'mean_squared_error', \
    metrics = [ 'mean_squared_error' ] )

  model.summary()
#
#  Train the model.
#  Use 20% of training data for validation.
#
  print ( '' )
  print ( 'Training:' )
  print ( '' )

  model.fit ( train_data, train_targets, validation_split = 0.20, epochs = 40 )
#
#  Test
#
  print ( '' )
  print ( 'Testing:' )
  print ( '' )
  print ( '  Case  True  Estimate' )
  test_predictions = model.predict ( test_data )
  for i in range ( test_num ):
    print ( '' )
    print ( '  %3d:' % ( i ) )
    for j in range ( target_num ):
      print ( '  %8.4f  %8.4f' % ( test_targets[i,j], test_predictions[i,j]  ) )

  import matplotlib.pyplot as plt
#
#  Plot true versus estimated values.
#
  j = 0
  xmin = np.min ( test_targets[:,j] )
  xmax = np.max ( test_targets[:,j] )
  plt.scatter ( test_targets[:,j], test_predictions[:,j] )
  plt.plot ( [xmin,xmax], [xmin,xmax], color = 'r', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<--- True value --->', fontsize = 16 )
  plt.ylabel ( '<--- Estimated value --->', fontsize = 16 )
  plt.title ( 'paraheat_gaussian: Estimation of parameter xc', fontsize = 16 )
  filename = 'paraheat_gaussian_parameters_xc.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( )

  j = 1
  xmin = np.min ( test_targets[:,j] )
  xmax = np.max ( test_targets[:,j] )
  plt.scatter ( test_targets[:,j], test_predictions[:,j] )
  plt.plot ( [xmin,xmax], [xmin,xmax], color = 'r', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<--- True value --->', fontsize = 16 )
  plt.ylabel ( '<--- Estimated value --->', fontsize = 16 )
  plt.title ( 'paraheat_gaussian: Estimation of parameter yc', fontsize = 16 )
  filename = 'paraheat_gaussian_parameters_yc.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( )

  j = 2
  xmin = np.min ( test_targets[:,j] )
  xmax = np.max ( test_targets[:,j] )
  plt.scatter ( test_targets[:,j], test_predictions[:,j] )
  plt.plot ( [xmin,xmax], [xmin,xmax], color = 'r', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<--- True value --->', fontsize = 16 )
  plt.ylabel ( '<--- Estimated value --->', fontsize = 16 )
  plt.title ( 'paraheat_gaussian: Estimation of parameter sc', fontsize = 16 )
  filename = 'paraheat_gaussian_parameters_sc.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( )

  j = 3
  xmin = np.min ( test_targets[:,j] )
  xmax = np.max ( test_targets[:,j] )
  plt.scatter ( test_targets[:,j], test_predictions[:,j] )
  plt.plot ( [xmin,xmax], [xmin,xmax], color = 'r', linewidth = 3 )
  plt.grid ( True )
  plt.xlabel ( '<--- True value --->', fontsize = 16 )
  plt.ylabel ( '<--- Estimated value --->', fontsize = 16 )
  plt.title ( 'paraheat_gaussian: Estimation of parameter vc', fontsize = 16 )
  filename = 'paraheat_gaussian_parameters_vc.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'paraheat_gaussian_parameters:' )
  print ( '  Normal end of execution.' )
  return

def timestamp ( ):

#*****************************************************************************80
#
## TIMESTAMP prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    06 April 2013
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    None
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return None

if ( __name__ == '__main__' ):
  timestamp ( )
  paraheat_gaussian_parameters ( 'xcycscvc2000.txt' )
  timestamp ( )