RBFN Settings

• RBFN().c - float, c value used by the basis functions in the calculations of weight factors (default c=0.5).
• RBFN().rbf_type - string, type of radial basis function to use, options are "ThinPlate", "Gaussian", "Multiquadric", "InvMultiquadric", and "InvGaussian" (default is "Multiquadric")
• RBFN().trainsets - integer, number of data sets to be used to train the surrogate model, the remaining sets will be used for model validation (default is 1)
• RBFN().errform - string, form of error to be used in validation, options are "rmse" and "avg" (default is "rmse")

Built in Functions

• RBFN().get_data(datafile,pointfile)
  • datafile - string, filename (or path) to file containing the data to be used for training. This is a *.csv file with each explicit dataset as columns. Each column is the y-data for each parameter set in the pointfile.
  • pointfile - "string", filename (or path) to file containing the design points for training. This is a *.csv file with each parameter set as rows (m samples by n parameters).
  • This function reads the files and returns x, y, and points.
• RBFN().hypercube_norm(points)
  • points - m x n array, m inputs by n parameters. These are the points used to explore the model input space.
  • This function maps the input parameters of points to the range [0,1] for each input parameter column using maximum and minimum values. It returns the normalized points in an m x n array.
• RBFN().optimize(y,points)
  • y and points array returned by RBFN().get_data(datafile,pointfile)
  • Uses scipy.optimize.minimize to minimize the validation error of the RBFN surrogate. Only functions if there some datasets are reserved for model validation by setting RBFN().trainsets to be less than the number of datasets available.
  • Once optimized, the RBFN is ready to make predictions without having to call the RBFN().train method.
• RBFN().train(y,points)
  • y and points array returned by RBFN().get_data(datafile,pointfile)
  • Uses scipy.optimize.minimize to minimize the validation error of the RBFN surrogate. Only functions if some datasets are reserved for model validation by setting RBFN().trainsets to be less than the number of datasets available.
• RBFN().predict(test)
  • target 1 x n NumPy array
  • Function used to predict the behavior of the test point. Returns a vector of data the same length as the data used to train it.
• RBFN().report()
  • Reports current RBFN settings, and if trained, the current validation error

Python Code Implementation

      		'''
Radial Basis Function Network (RBFN) class
'''

Copyright (c) 2018 Justin M. Hughes

Permission is hereby granted, free of charge, to any person
obtaining a copy of this software and associated documentation
files (the "Software"), to deal in the Software without
restriction, including without limitation the rights to use,
copy, modify, merge, publish, distribute, sublicense, and/or
sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following
conditions:

The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES 
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
OTHER DEALINGS IN THE SOFTWARE.

'''

import numpy as np
from scipy.optimize import minimize

class RBFN:
    
    def __init__(self):
        self.c = 0.5
        self.rbf_type = 'Multiquadric'
        self.lda = []
        self.training_data = []
        self.target_data = []
        self.training_points = []
        self.target_points = []
        self.valerr = 'NA'
        self.trainsets = 1
        self.errform = 'rmse'
    
    def report(self):
        print("\n----- Radial Basis Function Network -----")
        print("rbf_type  : %s" %self.rbf_type)
        print("trainsets : %s" %self.trainsets)
        print("c         : %s" %self.c)
        print("errform   : %s" %self.errform)
        print("valerr    : %s\n" %self.valerr)

    def split(self,ydata,points):
        self.training_data = ydata[:,:self.trainsets]
        self.target_data = ydata[:,self.trainsets:]
        self.training_points = points[:self.trainsets,:]
        self.target_points = points[self.trainsets:,:]

    def get_norm(self,vec):
        # Returns the 2-norm (Euclidean distance) of vec
        a = 0
        for i in range(0,len(vec)):
            a += vec[i]**2.0
        if a == 0:
            return 0
        else:
            return np.sqrt(a)
            
    def get_phi(self,r,c,rbf_type):
        if rbf_type == 'ThinPlate':
            if r==0:
                return 0
            else:
                return np.log(c*r)*r*r  # 1 Thin Plate
        elif rbf_type == 'Gaussian': # 2 Gaussian
            return np.exp(-1.0*c*r*r)
        elif rbf_type == 'Multiquadric' or rbf_type == 'InvMultiquadric':
            phi = np.sqrt(r*r + c*c) # 3 Multiquadric
            if rbf_type == 'InvMultiquadric':
                phi = 1.0/phi  # 4 Inverse Multiquadric
            return phi
        elif rbf_type == 'InvGaussian':
            return 1.0/np.exp(-1.0*c*r*r) # 5 Inverse Gaussian
        elif rbf_type == 'Custom':
            phi = np.sqrt(r*r + c*c)
            return  np.sqrt((phi**2.0 + (1.0/phi)**2.0)/2.0)# custom

    def lda_matrix(self,x,y):
        # Get number of training points and parameters
        n = x.shape[0]
        y = y.T

        #Set up phi storage variable
        A = np.zeros((n,n))

        # Get phi matrix (radii)
        for i in range(0,n):
            for j in range(0,n):
                r = self.get_norm(x[j,:]-x[i,:])
                A[i,j] = self.get_phi(r,self.c,self.rbf_type)

        self.lda = np.zeros((y.shape[0],y.shape[1]))
        # Solve for lambda
        for i in range(0,len(y[0,:])):
            try:
                self.lda[:,i] = np.dot(np.linalg.inv(A),y[:,i])
            except:
                #print "Phi matrix is singular, using pseudoinverse..."
                self.lda[:,i] = np.dot(np.linalg.pinv(A),y[:,i])

    def predict(self,x_test):
        # normalize each dimension to hypercube
        x_test = x_test.reshape(1,len(x_test))
        y_pred = np.zeros((self.lda.shape[1],1),dtype=float)
    
        #Get RBFN prediction for test point at each x
        for j in range(0,len(self.lda[0,:])):
            for i in range(0,len(self.lda[:,0])):
                r = self.get_norm(x_test[0,:] - self.training_points[i,:])
                phi = self.get_phi(r,self.c,self.rbf_type)
                y_pred[j,0] += self.lda[i,j]*phi
                
        return y_pred

    def hypercube_norm(self,points):
        pmax = np.nanmax(points,axis=0)
        pmin = np.nanmin(points,axis=0)
        return (points-pmin)/(pmax-pmin)

    def train(self,ydata,points):
        self.split(ydata,points)
        self.lda_matrix(self.training_points,self.training_data)
        self.pred_err()


    def get_data(self,filename,doe_points):
        data = np.genfromtxt(filename,delimiter=',',skip_header=1)
        xdata = np.array(data[:,0]).reshape(len(data[:,0]),1)
        ydata = np.array(data[:,1:])

        doe = np.genfromtxt(doe_points,delimiter=',',skip_header=1)
        sets = len(ydata[0,:])

        return xdata,ydata,doe
    
    def get_err(self,test,true):
        diff = []
        for i in range(0,len(test)):
            diff.append((float(test[i]-true[i])/true[i]))
        
        if self.errform == 'avg':
            return np.abs(np.nanmean(diff))
        elif self.errform == 'rmse':
            diff = [x**2.0 for x in diff]
            sumdiff = np.nansum(diff)
            return np.sqrt(sumdiff/len(diff))

    def pred_err(self):
        err = []
        for i in range(0,self.target_points.shape[0]):
            x_test = self.target_points[i,:]
            ypred = self.predict(x_test)
            err.append(self.get_err(list(ypred),list(self.target_data[:,i])))
        self.valerr = np.nanmean(err)
        
    def cvalminfunc(self,cval,*args):
        self.c = cval
        ydata,points = args
        
        self.train(ydata,points)
        err = self.valerr
        if cval < 0.0 or cval > 1.0:
            return 1.0e06
        else:
            return err
    
    def optimize(self,ydata,points,method='Nelder-Mead',tol=1e-06):
        print("Optimizing RBFN c value. Please wait...")
        c0 = 0.5
        args = (ydata,points)
        minerr_cval = minimize(self.cvalminfunc,c0,args,method=method,tol=tol)
        print(minerr_cval)
        
        self.c = float(minerr_cval.x)
        self.train(ydata,points)