working on regression kriging

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
 

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()
 

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()
+

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()
+
+
+
+
+
+

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.

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