From 1674c1de8669773cf81846d8cbb1fd4d46d8d62c Mon Sep 17 00:00:00 2001 From: cpaulson Date: Thu, 11 Jun 2015 10:34:51 +0100 Subject: [PATCH] working on regression kriging --- examples/2d_regression_Kriging.py | 14 +- examples/3d_Simple_Train.py | 20 +-- examples/coKriging.py | 37 +++++ pyKriging/coKriging.py | 218 ++++++++++++++++++++++++++++++ pyKriging/krige.py | 5 +- pyKriging/matrixops.py | 9 +- pyKriging/regressionkrige.py | 15 +- 7 files changed, 290 insertions(+), 28 deletions(-) create mode 100644 examples/coKriging.py create mode 100644 pyKriging/coKriging.py diff --git a/examples/2d_regression_Kriging.py b/examples/2d_regression_Kriging.py index 76bf536..da050dc 100644 --- a/examples/2d_regression_Kriging.py +++ b/examples/2d_regression_Kriging.py @@ -24,13 +24,13 @@ print k.Lambda # k.snapshot() # -# for i in range(5): -# newpoints = k.infill(2) -# for point in newpoints: -# print 'Adding point {}'.format(point) -# k.addPoint(point, testfun(point)[0]) -# k.train() -# k.snapshot() +for i in range(5): + newpoints = k.infill(5) + for point in newpoints: + print 'Adding point {}'.format(point) + k.addPoint(point, testfun(point)[0]) + k.train() + # k.snapshot() # # # #And plot the model diff --git a/examples/3d_Simple_Train.py b/examples/3d_Simple_Train.py index c9b7cb7..d96325e 100644 --- a/examples/3d_Simple_Train.py +++ b/examples/3d_Simple_Train.py @@ -4,38 +4,38 @@ from pyKriging.testfunctions import testfunctions -## The Kriging model starts by defining a sampling plan, we use an optimal Lattin Hypercube here +# The Kriging model starts by defining a sampling plan, we use an optimal Latin Hypercube here sp = samplingplan(3) X = sp.optimallhc(30) -## Next, we define the problem we would like to solve +# Next, we define the problem we would like to solve testfun = testfunctions().squared y = testfun(X) -## Now that we have our initial data, we can create an instance of a kriging model +# Now that we have our initial data, we can create an instance of a kriging model k = kriging(X, y, testfunction=testfun, testPoints=300) -## The model is then trained +# The model is then trained k.train() k.snapshot() # It's typically beneficial to add additional points based on the results of the initial training # The infill method can be used for this # In this example, we will add nine points in three batches. The model gets trained after each stage -for i in range( 10 ): +for i in range(10): print k.history['rsquared'][-1] - print 'Infill itteration {0}'.format( i+1 ) + print 'Infill iteration {0}'.format(i + 1) infillPoints = k.infill(10) - ## Evaluate the infill points and add them back to the Kriging model + # Evaluate the infill points and add them back to the Kriging model for point in infillPoints: - print 'Adding point {}'.format( point ) + print 'Adding point {}'.format(point) k.addPoint(point, testfun(point)[0]) - ## Retrain the model with the new points added in to the model + # Retrain the model with the new points added in to the model k.train() k.snapshot() -## Once the training of the model is complete, we can plot the results +# Once the training of the model is complete, we can plot the results k.plot() diff --git a/examples/coKriging.py b/examples/coKriging.py new file mode 100644 index 0000000..d1dddb7 --- /dev/null +++ b/examples/coKriging.py @@ -0,0 +1,37 @@ +__author__ = 'cpaulson' + +import sys +sys.path.insert(0,'.') +sys.path.insert(0,'..') +from pyKriging import coKriging + +import numpy as np + +def cheap(X): + + A=0.5 + B=10 + C=-5 + D=0 + + print X + print ((X+D)*6-2) + return A*np.power( ((X+D)*6-2), 2 )*np.sin(((X+D)*6-2)*2)+((X+D)-0.5)*B+C + +def expensive(X): + return np.power((X*6-2),2)*np.sin((X*6-2)*2) + + +Xe = np.array([0, 0.4, 0.6, 1]) +Xc = np.array([0.1,0.2,0.3,0.5,0.7,0.8,0.9,0,0.4,0.6,1]) + +yc = cheap(Xc) +ye = expensive(Xe) + +ck = coKriging.coKriging(Xc, yc, Xe, ye) +ck.thetac = np.array([1.2073]) +print ck.Xc +ck.updateData() +ck.updatePsi() +ck.neglnlikehood() + diff --git a/pyKriging/coKriging.py b/pyKriging/coKriging.py new file mode 100644 index 0000000..1e8bd61 --- /dev/null +++ b/pyKriging/coKriging.py @@ -0,0 +1,218 @@ +__author__ = 'cpaulson' + +import numpy as np +from numpy.matlib import rand,zeros,ones,empty,eye +from pyKriging import kriging + +class coKriging(): + def __init__(self, Xc, yc, Xe, ye): + + # Create the data arrays + self.Xc = np.atleast_2d(Xc).T + self.yc = yc + self.nc = self.Xc.shape[0] + + self.Xe = np.atleast_2d(Xe).T + + self.ye = ye + self.ne = self.Xe.shape[0] + + # rho regression parameter + self.rho = 1.9961 + self.reorder_data() + # self.traincheap() + + self.k = self.Xc.shape[1] + # if self.Xe.shape[1] != self.Xc.shape[1]: + # print 'Xc and Xe must have the same number of design variables. Fatal error -- Exiting...' + # exit() + + # Configure the hyperparameter arrays + self.thetad = np.ones(self.k) + self.thetac = None + # self.thetac = self.kc.theta + + self.pd = np.ones(self.k) * 2. + # self.pc = self.kc.pl + self.pc = np.ones(self.k) * 2. + + # Matrix Operations + self.one=ones([self.ne+self.nc,1]) + self.y=[self.yc, self.ye] + + print 'here1' + + def reorder_data(self): + xe = [] + ye = [] + xc = [] + yc = [] + + Xd = [] + yd = [] + + for enu,entry in enumerate(self.Xc): + if entry in self.Xe: + print 'Found this value in XE!!' + for enu1,test in enumerate(self.Xe): + # if entry[0] == test[0] and entry[1] == test[1]: + if entry == test: + xe.append(test.tolist()) + ye.append(self.ye[enu1].tolist()) + xc.append(entry.tolist()) + yc.append(self.yc[enu].tolist()) + Xd.append(entry.tolist()) + yd.append(self.ye[enu1].tolist() - self.rho * self.yc[enu].tolist()) + break + + else: + xc.insert(0,entry.tolist()) + yc.insert(0,self.yc[enu].tolist()) + + self.Xe = np.array(xe) + self.ye = np.array(ye) + self.Xc = np.array(xc) + self.yc = np.array(yc) + self.Xd = np.array(Xd) + self.yd = np.atleast_2d(np.array(yd)) + + + def updateData(self): + self.nc = self.Xc.shape[0] + self.ne = self.Xe.shape[0] + self.distanceXc() + self.distanceXe() + self.distanceXcXe() + + def traincheap(self): + self.kc = kriging(self.Xc, self.yc) + self.kc.train() + print + + + def distanceXc(self): + self.distanceXc = np.zeros((self.nc,self.nc, self.k)) + for i in range( self.nc ): + for j in xrange(i+1,self.nc): + self.distanceXc[i][j] = np.abs((self.Xc[i]-self.Xc[j])) + + def distanceXe(self): + self.distanceXe = np.zeros((self.ne,self.ne, self.k)) + for i in range( self.ne ): + for j in xrange(i+1,self.ne): + self.distanceXe[i][j] = np.abs((self.Xe[i]-self.Xe[j])) + + def distanceXcXe(self): + self.distanceXcXe = np.zeros((self.nc,self.ne, self.k)) + for i in range( self.nc ): + for j in xrange(self.ne): + self.distanceXcXe[i][j] = np.abs((self.Xc[i]-self.Xe[j])) + + + def updatePsi(self): + self.PsicXc = np.zeros((self.nc,self.nc), dtype=np.float) + self.PsicXe = np.zeros((self.ne,self.ne), dtype=np.float) + self.PsicXcXe = np.zeros((self.nc,self.ne), dtype=np.float) + # + # print self.thetac + # print self.pc + # print self.distanceXc + newPsicXc = np.exp(-np.sum(self.thetac*np.power(self.distanceXc,self.pc), axis=2)) + print newPsicXc[0] + self.PsicXc = np.triu(newPsicXc,1) + self.PsicXc = self.PsicXc + self.PsicXc.T + np.mat(eye(self.nc))+np.multiply(np.mat(eye(self.nc)),np.spacing(1)) + self.UPsicXc = np.linalg.cholesky(self.PsicXc) + self.UPsicXc = self.UPsicXc.T + print self.PsicXc[0] + print self.UPsicXc + exit() + + newPsicXe = np.exp(-np.sum(self.thetac*np.power(self.distanceXe,self.pc), axis=2)) + self.PsicXe = np.triu(newPsicXe,1) + self.PsiXe = self.PsicXe + self.PsicXe.T + np.mat(eye(self.ne))+np.multiply(np.mat(eye(self.ne)),np.spacing(1)) + self.UPsicXe = np.linalg.cholesky(self.PsicXe) + self.UPsicXe = self.UPsicXe.T + + newPsiXeXc = np.exp(-np.sum(self.thetad*np.power(self.distanceXcXe,self.pd), axis=2)) + self.PsicXcXe = np.triu(newPsiXeXc,1) + + + def neglnlikehood(self): + a = np.linalg.solve(self.UPsicXc.T, np.matrix(self.yc).T) + b = np.linalg.solve( self.UPsicXc, a ) + c = ones([self.nc,1]).T * b + + d = np.linalg.solve(self.UPsicXc.T, ones([self.nc,1])) + e = np.linalg.solve(self.UPsicXc, d) + f = ones([self.nc,1]).T * e + + self.muc = c/f + # This only works if yc is transposed, then its a scalar under two layers of arrays. Correct? Not sure + + print 'y',self.yd.T + a = np.linalg.solve(self.UPsicXe.T, self.yd) + print 'a',a + b = np.linalg.solve(self.UPsicXe, a) + print 'b', b + c = ones([self.ne,1]) * b + print 'c', c + + d = np.linalg.solve(self.UPsicXe.T, ones([self.ne,1], dtype=float)) + print d + + e = np.linalg.solve(self.UPsicXe, d) + print e + + f = ones([self.ne,1]).T * e + print f + + self.mud= c/f + + + a = np.linalg.solve(self.UPsicXc.T,(self.yc-ones([self.nc,1])*self.muc))/self.nc + b = np.linalg.solve(self.UPsicXc, a) + self.SigmaSqrc=(self.yc-ones([self.nc,1])*self.muc).T* b + + + + print self.ne + print self.mud + print self.UPsicXe.T + a = np.linalg.solve(self.UPsicXe.T,(self.yd-ones([self.ne,1])*self.mud))/self.ne + b = np.linalg.solve(self.UPsicXe, a) + self.SigmaSqrd=(self.yd-ones([self.ne,1])*self.mud).T* b + + self.C=np.array([self.SigmaSqrc*self.PsicXc, self.rho*self.SigmaSqrc*self.PsicXcXe, self.rho*self.SigmaSqrc*self.PsicXeXc, np.power(self.rho,2)*self.SigmaSqrc*self.PsicXe+self.SigmaSqrd*self.PsidXe]) + np.reshape(c,[2,2]) + + self.UC = np.linalg.cholesky(self.C) + + # self.mu=(self.one.T *(self.UC\(self.UC.T\y)))/(one'*(ModelInfo.UC\(ModelInfo.UC'\one))); + + + +def fc(X): + return np.power(X[:,0], 2) + X[:,0] + np.power(X[:,1], 2) + X[:,1] +def fe(X): + return np.power(X[:,0], 2) + np.power(X[:,1], 2) + +if __name__=='__main__': + import samplingplan + import random + sp = samplingplan.samplingplan(2) + X = sp.optimallhc(20) + Xe = np.array( random.sample(X, 6) ) + + yc = fc(X) + ye = fe(Xe) + + ck = coKriging(X, yc, Xe, ye) + ck.updateData() + ck.updatePsi() + ck.neglnlikehood() + + + + + + diff --git a/pyKriging/krige.py b/pyKriging/krige.py index 23baaaf..4264e00 100644 --- a/pyKriging/krige.py +++ b/pyKriging/krige.py @@ -116,8 +116,10 @@ def normalizeData(self): for i in range(self.k): self.normRange.append([min(self.X[:, i]), max(self.X[:, i])]) + print self.X for i in range(self.n): self.X[i] = self.normX(self.X[i]) + print self.X self.ynormRange.append(min(self.y)) self.ynormRange.append(max(self.y)) @@ -285,7 +287,7 @@ def infill(self, points, method='error'): final_pop.sort(reverse=True) newpoint = final_pop[0].candidate returnValues[i][:] = newpoint - self.addPoint(returnValues[i], self.predict(returnValues[i]), norm=True) + self.addPoint(returnValues[i], self.predict(returnValues[i]), norm=False) self.X = np.copy(initX) self.y = np.copy(inity) @@ -314,7 +316,6 @@ def generate_population(self, random, args): chromosome.append(random.uniform(lo, hi)) return chromosome - def no_improvement_termination(self, population, num_generations, num_evaluations, args): """Return True if the best fitness does not change for a number of generations of if the max number of evaluations is exceeded. diff --git a/pyKriging/matrixops.py b/pyKriging/matrixops.py index 8075444..06f7ace 100644 --- a/pyKriging/matrixops.py +++ b/pyKriging/matrixops.py @@ -34,9 +34,9 @@ def updatePsi(self): def regupdatePsi(self): self.Psi = np.zeros((self.n,self.n), dtype=np.float) self.one = np.ones(self.n) - for i in xrange(self.n): - for j in xrange(i+1,self.n): - self.Psi[i,j]=np.exp(-np.sum(self.theta*np.power(np.abs((self.X[i]-self.X[j])),self.pl))) + self.psi = np.zeros((self.n,1)) + newPsi = np.exp(-np.sum(self.theta*np.power(self.distance,self.pl), axis=2)) + self.Psi = np.triu(newPsi,1) self.Psi = self.Psi + self.Psi.T + eye(self.n) + eye(self.n) * (self.Lambda) self.U = np.linalg.cholesky(self.Psi) self.U = np.matrix(self.U.T) @@ -85,6 +85,9 @@ def predicterr_normalized(self,x): try: SSqr=self.SigmaSqr*(1-self.psi.T.dot(np.linalg.solve(self.U, np.linalg.solve(self.U.T,self.psi)))) except Exception, e: + print self.U.shape + print self.SigmaSqr.shape + print self.psi.shape print Exception,e pass diff --git a/pyKriging/regressionkrige.py b/pyKriging/regressionkrige.py index 77d85d6..7076fe7 100644 --- a/pyKriging/regressionkrige.py +++ b/pyKriging/regressionkrige.py @@ -27,7 +27,7 @@ def __init__(self, X, y, testfunction=None, name='', testPoints=None, **kwargs): self.k = self.X.shape[1] self.theta = np.ones(self.k) self.pl = np.ones(self.k) * 2. - self.Lambda = 1 + self.Lambda = 0 self.sigma = 0 self.normRange = [] self.ynormRange = [] @@ -144,6 +144,7 @@ def addPoint(self, newX, newy, norm=True): try: self.updateModel() except: + print 'Couldnt update the model with these hyperparameters, retraining' self.train() else: break @@ -181,14 +182,15 @@ def predict(self, X): X = self.normX(X) return self.inversenormy(self.predict_normalized(X)) - def predict_var(self, X): + def predict_var(self, X, norm=True): ''' The function returns the model's predicted 'error' at this point in the model. :param X: new design variable to evaluate, in physical world units :return: Returns the posterior variance (model error prediction) ''' X = copy.deepcopy(X) - X = self.normX(X) + if norm: + X = self.normX(X) # print X, self.predict_normalized(X), self.inversenormy(self.predict_normalized(X)) return self.predicterr_normalized(X) @@ -251,7 +253,7 @@ def infill_objective_ei(self,candidates, args): fitness.append(-1 * self.expimp(entry)) return fitness - def infill(self, points, method='error'): + def infill(self, points, method='error', addPoint=True): ''' The function identifies where new points are needed in the model. :param points: The number of points to add to the model. Multiple points are added via imputation. @@ -286,7 +288,8 @@ def infill(self, points, method='error'): final_pop.sort(reverse=True) newpoint = final_pop[0].candidate returnValues[i][:] = newpoint - self.addPoint(returnValues[i], self.predict(returnValues[i]), norm=True) + if addPoint: + self.addPoint(returnValues[i], self.predict(returnValues[i]), norm=True) self.X = np.copy(initX) self.y = np.copy(inity) @@ -679,7 +682,7 @@ def calcuatemeanMSE(self, p2s=200, points=None): points = self.sp.rlh(p2s) values = np.zeros(len(points)) for enu, point in enumerate(points): - values[enu] = self.predict_var(point) + values[enu] = self.predict_var(point, norm=False) return np.mean(values), np.std(values) def snapshot(self):