text: import numpy as np
  import matplotlib.pyplot as plt
  import matplotlib.animation as animation
  from mpl_toolkits.mplot3d import Axes3D
  import timeit
  import time
  import datetime
  from tqdm import tqdm
  import os 
  
  def sigmoid(x):
      return 1/(1+np.exp(-x))
  
  def Gaussian(input_array, m, sigma):    
      return np.exp( - np.linalg.norm(input_array-m, axis=1, keepdims=True)**2 / (2 * sigma**2) ) 
  
  np.random.seed(6666)
  x1 = np.linspace(-5,5, 400)
  x2 = np.linspace(-5,5, 400)
  np.random.shuffle(x1) 
  np.random.shuffle(x2) 
  
  d = x1**2 + x2**2    
  
  # Normalize d to range 0.2~0.8
  d_max = np.max(d)
  d_min = np.min(d)
  d = (d-d_min)/(d_max-d_min)*(0.8-0.2)+0.2 
  
  #---------------- Input vector -----------------------------
  num_in  = 2
  
  #---------------- Radial Basis phi -----------------------
  num_phi = 3
  sigma_origin     = sigma = np.random.uniform(0,2,[num_phi,1])
  phi_out = np.zeros([num_phi,1])
  m_origin     = m = np.random.uniform(-1,1,[num_phi,num_in])
  #---------------- Output ---------------------------------
  num_out  = 1
  bias_out = np.random.uniform(-0.5,0.5,[num_out,1])
  w_out    = np.random.uniform(-0.5,0.5,[num_phi,num_out])
  
  #---------------- Parameter --------------------------
  eta   = 1
  mom   = 0.9
  epoch = 6000
  
  Eav_train = np.zeros([epoch])
  Eav_test = np.zeros([epoch])
  
  dw_out    = temp1 = np.zeros([num_phi,num_out]) 
  dbias_out = temp2 = np.zeros([num_out,1])
  
  dm        = temp3 = np.zeros([num_phi,num_in])
  dsigma    = temp4 = np.zeros([num_phi,1])
  
  #---------------  Save png parameter ------------------
  os.makedirs("Data gif", exist_ok=True)
  os.makedirs("Gaussian gif", exist_ok=True)
  
  Data_gif = True
  Gaussian_gif = True 
  
  iter1 = 0 # Data set
  iter2 = 0 # Gaussian
  
  #---------------- Traning ----------------------------
  t0 = timeit.default_timer()
  now = datetime.datetime.now().strftime("%Y-%m-%d_%H-%M-%S")
  pbar = tqdm(total =epoch)
  plt.ion()
  plt.show()
  
  for i in range(epoch):         
  
      #--------------- Feed Forward  -------------------    
      e = np.zeros([300])
      E_train = np.zeros([300])
      for j in range(300):
          X   = np.array([x1[j],x2[j]]).reshape(2,1)
  
          XX = np.array(X)
          for _ in range(num_phi-1):
              XX = np.append(XX, X, axis=1)
          XX = np.transpose(XX)
          dx_m = XX - m
          
          phi_out = Gaussian(XX,m,sigma)
          out = sigmoid(np.dot(np.transpose(phi_out),w_out) + bias_out)
  
          #--------------- Back Propagation-----------------
          e[j]   = (d[j]-out) 
          E_train[j] = 0.5 * e[j]**2
          locg_k = e[j] * (out*(1-out))
          temp2  = temp2 + mom * dbias_out + eta * locg_k * 1     
          temp1  = temp1 + mom * dw_out + eta * locg_k * phi_out                  
  
          temp3  = temp3 + mom * dm + eta * e[j] * w_out * phi_out / sigma**2 * (dx_m)              
          temp4  = temp4 + mom * dsigma + eta * e[j] * w_out * phi_out / sigma**3 * np.linalg.norm(dx_m, axis=1, keepdims=True)
  
      #---------- Average delta weight -----------------
      dbias_out = temp2/300
      dw_out    = temp1/300       
      dm        = temp3/300
      dsigma    = temp4/300  
  
      temp1 = np.zeros([num_phi,num_out]) 
      temp2 = np.zeros([num_out,1])
      temp3 = np.zeros([num_phi,num_in])
      temp4 = np.zeros([num_phi,1])
  
      #----------  New weight --------------------------
      bias_out  = bias_out + dbias_out
      w_out     = w_out + dw_out              
      m         = m + dm
      sigma     = sigma + dsigma
  
      #---------- Eave_train
      Eav_train[i]  = np.mean(E_train)
  
      #---------- Test data loss ---------------     
      E_test = np.zeros([100])
      for j in range(100):
          X   = np.array([x1[300+j],x2[300+j]]).reshape(2,1)
  
          XX = np.array(X)
          for _ in range(num_phi-1):
              XX = np.append(XX, X, axis=1)
          XX = np.transpose(XX)
          dx_m = XX - m
  
          phi_out = Gaussian(XX,m,sigma)
          out = sigmoid(np.dot(np.transpose(phi_out),w_out) + bias_out)
          E_test = 0.5*( d[300+j] - out )**2      
  
      Eav_test[i] = np.mean(E_test)
  
      #----- Update the percentage bar
      if i % 1000 == 0 and i!=0:
          pbar.update(1000)
  
      #----  Display the current predicted results of the model ----
      if i% 100 == 0 and Data_gif == True:
          iter1 += 1
          y_predict = np.zeros([100])    
          for k in range(100):
              X      = np.array([x1[300+k],x2[300+k]]).reshape(2,1)
              XX = np.array(X)
              for _ in range(num_phi-1):
                  XX = np.append(XX, X, axis=1)
              XX = np.transpose(XX)     
  
              phi_out = Gaussian(XX,m,sigma)
              out = sigmoid(np.dot(np.transpose(phi_out),w_out) + bias_out)
              y_predict[k] = out
  
          y_predict = (y_predict-0.2)/(0.8-0.2)*(d_max-d_min)+d_min
          fig = plt.figure(num='Animation')
          ax1 = fig.add_subplot(111, projection='3d')
          ax1.scatter(x1[300:], x2[300:], y_predict[:],c='g', marker='o', s=15)        
          ax1.set_xlabel('x1')
          ax1.set_ylabel('x2')
          ax1.set_zlabel('y')
          plt.title('y = x1^2 +x2^2')        
          plt.pause(0.001)
          plt.savefig('Data gif/'+str(iter1)+'.png')
      #---------- display the current gaussian function
      if num_phi == 3 and i% 100 == 0 and Gaussian_gif == True :
          iter2 += 1
          #------------ Plot gaussian ---------
          X1 = np.linspace(-5,5,400)
          X2 = np.linspace(-5,5,400)
          G = np.zeros([400,400,num_phi])
          [X1,X2] = np.meshgrid(X1,X2)
  
          for i in range(400) :
              for j in range(400):
                  X   = np.array([X1[i,j],X2[i,j]]).reshape(2,1)
                  XX = np.array(X)
                  for _ in range(num_phi-1):
                      XX = np.append(XX, X, axis=1)
                  XX = np.transpose(XX)
                  phi_out = Gaussian(XX,m,sigma)
                  G[i,j,:] = phi_out[:,0]         
  
  
          fig3 = plt.figure(num ='Gaussian',figsize=(17,6))
          plt.suptitle('training ...',fontsize = 16)
  
          ax = fig3.add_subplot(131, projection='3d')
          ax.plot_surface(X1, X2, G[:,:,0],rstride=10, cstride=10, cmap='rainbow')
          ax.set_xlabel('x1')
          ax.set_ylabel('x2')
  
          plt.title('Gaussian 1')
  
          ax2 = fig3.add_subplot(132, projection='3d')
          ax2.plot_surface(X1, X2, G[:,:,1],rstride=10, cstride=10, cmap='rainbow')
          ax2.set_xlabel('x1')
          ax2.set_ylabel('x2')
  
          plt.title('Gaussian 2')
  
          ax3 = fig3.add_subplot(133, projection='3d')
          ax3.plot_surface(X1, X2, G[:,:,2],rstride=10, cstride=10, cmap='rainbow')
          ax3.set_xlabel('x1')
          ax3.set_ylabel('x2')
  
          plt.title('Gaussian 3')
          plt.pause(0.001)
          plt.savefig('Gaussian gif/'+str(iter2)+'.png')
  
  pbar.close()
  t1 = (timeit.default_timer()-t0)
  print('Training time: {} min'.format((t1/60)))
  
  #--------- Predict data --------------
  y_predict = np.zeros([100])
  E_predict = np.zeros([100])
  for j in range(100):
          X      = np.array([x1[300+j],x2[300+j]]).reshape(2,1)
  
          XX = np.array(X)
          for _ in range(num_phi-1):
              XX = np.append(XX, X, axis=1)
          XX = np.transpose(XX)
          dx_m = XX - m
  
          phi_out = Gaussian(XX,m,sigma)
          out = sigmoid(np.dot(np.transpose(phi_out),w_out) + bias_out)
          y_predict[j] = out
          E_predict[j] = 0.5*( d[300+j] - out )**2
  Eav_predict = np.mean(E_predict)
  
  #----------- Return the data they were normolized before ----------------------
  y_predict = (y_predict-0.2)/(0.8-0.2)*(d_max-d_min)+d_min  
  
  
  #------------ Record the result ------------------
  import csv
  table = [
      #['TimeStamp','Unit', 'Eta', 'Alpha','Training_loss','Predict_loss','Epoch','Time(min)'],
      [ now,num_phi, eta, mom, Eav_train[epoch-1], Eav_predict,epoch, int(t1/60)]    
  ]
  with open('RBF_output.csv', 'a', newline='') as csvfile:
      writer = csv.writer(csvfile)
      writer.writerows(table)
  
  #------------ Scattering Data set ----------------------
  #------- Return the data they were normolized before ---
  d = (d-0.2)/(0.8-0.2)*(d_max-d_min)+d_min 
  
  fig1 = plt.figure(num ='Data_Set', figsize=(10,5))
  ax = fig1.add_subplot(121, projection='3d')
  ax.scatter(x1[:300], x2[:300], d[:300], c='b', marker='o', s=5)
  ax.set_xlabel('x1')
  ax.set_ylabel('x2')
  ax.set_zlabel('y')
  plt.title('Training Data')
  
  ax2 = fig1.add_subplot(122, projection='3d')
  ax2.scatter(x1[300:], x2[300:], d[300:], c='r', marker='o', s=5)
  ax2.set_xlabel('x1')
  ax2.set_ylabel('x2')
  ax2.set_zlabel('y')
  plt.title('Testing Data')
  plt.show()
  
  
  
  #------------ Scattering y_Predict data -----------------
  fig2 = plt.figure(num = now, figsize=(14,6))
  ax1 = fig2.add_subplot(122, projection='3d')
  ax1.scatter(x1[300:], x2[300:], y_predict[:],c='g', marker='o', s=15)
  ax1.set_xlabel('x1')
  ax1.set_ylabel('x2')
  ax1.set_zlabel('y')
  plt.title('Predict Data : y = x1^2 +x2^2')
  
  #------------ plot training and testing Loss -----------------
  ax3 = fig2.add_subplot(121)
  ax3.set_xlabel('Epochs')
  ax3.set_ylabel('Loss')
  ax3.plot(range(epoch),Eav_train,label='Train Set :'+str(Eav_train[epoch-1]))
  ax3.plot(range(epoch),Eav_test, color='red', linewidth=1.0, linestyle='--', label='Test Set :'+ str(Eav_test[epoch-1]))
  plt.legend(loc='upper right')
  plt.title('Unit:'+str(num_phi)+', Eta:'+str(eta)+', Alpha:'+str(mom))
  plt.show()
  
  
  #----------- Plot init Gaussian -------------------------
  if num_phi ==3:
      x1 = np.linspace(-5,5,400)
      x2 = np.linspace(-5,5,400)
      G = np.zeros([400,400,num_phi])
      [X1,X2] = np.meshgrid(x1,x2)
  
      for i in range(400) :
          for j in range(400):
              X   = np.array([X1[i,j],X2[i,j]]).reshape(2,1)
              XX = np.array(X)
              for _ in range(num_phi-1):
                  XX = np.append(XX, X, axis=1)
              XX = np.transpose(XX)
              phi_out = Gaussian(XX,m_origin,sigma_origin)
              G[i,j,:] = phi_out[:,0]         
  
  
      fig3 = plt.figure(num ='Gaussian_init',figsize=(17,6))
      plt.suptitle('Before training ...',fontsize = 16)
  
      ax = fig3.add_subplot(131, projection='3d')
      ax.plot_surface(X1, X2, G[:,:,0],rstride=10, cstride=10, cmap='rainbow')
      ax.set_xlabel('x1')
      ax.set_ylabel('x2')
  
      plt.title('Gaussian 1')
  
      ax2 = fig3.add_subplot(132, projection='3d')
      ax2.plot_surface(X1, X2, G[:,:,1],rstride=10, cstride=10, cmap='rainbow')
      ax2.set_xlabel('x1')
      ax2.set_ylabel('x2')
  
      plt.title('Gaussian 2')
  
      ax3 = fig3.add_subplot(133, projection='3d')
      ax3.plot_surface(X1, X2, G[:,:,2],rstride=10, cstride=10, cmap='rainbow')
      ax3.set_xlabel('x1')
      ax3.set_ylabel('x2')
  
      plt.title('Gaussian 3')
  
      plt.show()
  
      #-------------- Plot trained Gaussian ----------- 
  
      for i in range(400) :
          for j in range(400):
              X   = np.array([X1[i,j],X2[i,j]]).reshape(2,1)
              XX = np.array(X)
              for _ in range(num_phi-1):
                  XX = np.append(XX, X, axis=1)
              XX = np.transpose(XX)
              phi_out = Gaussian(XX,m,sigma)
              G[i,j,:] = phi_out[:,0]         
  
  
      fig3 = plt.figure(num ='Gaussian_trained',figsize=(17,6))
      plt.suptitle('After training ...',fontsize = 16)
  
      ax = fig3.add_subplot(131, projection='3d')
      ax.plot_surface(X1, X2, G[:,:,0],rstride=10, cstride=10, cmap='rainbow')
      ax.set_xlabel('x1')
      ax.set_ylabel('x2')
      plt.title('Gaussian 1')
  
      ax2 = fig3.add_subplot(132, projection='3d')
      ax2.plot_surface(X1, X2, G[:,:,1],rstride=10, cstride=10, cmap='rainbow')
      ax2.set_xlabel('x1')
      ax2.set_ylabel('x2')
      plt.title('Gaussian 2')
  
      ax3 = fig3.add_subplot(133, projection='3d')
      ax3.plot_surface(X1, X2, G[:,:,2],rstride=10, cstride=10, cmap='rainbow')
      ax3.set_xlabel('x1')
      ax3.set_ylabel('x2')
      plt.title('Gaussian 3')
      plt.show()
  
  plt.ioff()
  plt.show()