Changeset 7493:a6b36625822a in orange


Ignore:
Timestamp:
02/04/11 18:05:41 (3 years ago)
Author:
matija <matija.polajnar@…>
Branch:
default
Convert:
1222327ccf1e94df714837d6a19c637cbd792bec
Message:

Complete MDS refactorization and documentation.

Location:
orange
Files:
6 added
3 edited

Legend:

Unmodified
Added
Removed
  • orange/Orange/__init__.py

    r7441 r7493  
    3333 
    3434import projection 
     35import projection.mds 
    3536import projection.som 
    3637 
  • orange/doc/Orange/rst/index.rst

    r7466 r7493  
    3838   Orange.network 
    3939    
     40   Orange.projection.mds 
    4041   orange.projection.som 
    4142 
  • orange/orngMDS.py

    r6538 r7493  
    1 # 
    2 # A Trashy (but working) MDS Library  
    3 # 
    4 # CVS Info: $Id$ 
    5 # 
    6 # V 1.01 
    7 #   - fixed the scaling bug in constructMatrixFromProx() (H. Harpending) 
    8 # 
    9 # V 1.0 
    10 # 
    11 # Aleks Jakulin (jakulin@acm.org) 2001-2002 
    12 # Ales Erjavec 2006 
    13  
    14  
    15 import orange 
    16 import orangeom as orangemds 
    17 from math import * 
    18 #from Numeric import * 
    19 #from LinearAlgebra import * 
    20 from numpy import * 
    21 from numpy.linalg import svd 
    22  
    23  
    24 KruskalStress=orangemds.KruskalStress() 
    25 SammonStress=orangemds.SammonStress() 
    26 SgnSammonStress=orangemds.SgnSammonStress() 
    27 SgnRelStress=orangemds.SgnRelStress() 
    28  
    29 PointList=orange.FloatListList 
    30 FloatListList=orange.FloatListList 
    31  
    32 def _mycompare((a,aa),(b,bb)): 
    33     if a==b: 
    34         return 0 
    35     if a<b: 
    36         return -1 
    37     else: 
    38         return 1 
    39              
    40 class PivotMDS(object): 
    41     def __init__(self, distances=None, pivots=50, dim=2, **kwargs): 
    42         self.dst = array([m for m in distances]) 
    43         self.n = len(self.dst) 
    44  
    45         if type(pivots) == type(1): 
    46             self.k = pivots 
    47             self.pivots = random.permutation(len(self.dst))[:pivots] 
    48             #self.pivots.sort() 
    49         elif type(pivots) == type([]): 
    50             self.pivots = pivots 
    51             #self.pivots.sort() 
    52             self.k = len(self.pivots) 
    53         else: 
    54             raise AttributeError('pivots') 
    55          
    56     def optimize(self): 
    57 #        # Classical MDS (Torgerson) 
    58 #        J = identity(self.n) - (1/float(self.n)) 
    59 #        B = -1/2. * dot(dot(J, self.dst**2), J) 
    60 #        w,v = linalg.eig(B) 
    61 #        tmp = zip([float(val) for val in w], range(self.n)) 
    62 #        tmp.sort() 
    63 #        w1, w2 = tmp[-1][0], tmp[-2][0] 
    64 #        v1, v2 = v[:, tmp[-1][1]], v[:, tmp[-2][1]] 
    65 #        return v1 * sqrt(w1), v2 * sqrt(w2)  
    66          
    67         # Pivot MDS 
    68         d = self.dst[[self.pivots]].T 
    69         C = d**2 
    70         # double-center d 
    71         cavg = sum(d, axis=0)/(self.k+0.0)      # column sum 
    72         ravg = sum(d, axis=1)/(self.n+0.0)    # row sum 
    73         tavg = sum(cavg)/(self.n+0.0)   # total sum 
    74         # TODO: optimize 
    75         for i in xrange(self.n): 
    76             for j in xrange(self.k): 
    77                 C[i,j] += -ravg[i] - cavg[j] 
    78          
    79         C = -0.5 * (C + tavg) 
    80         w,v = linalg.eig(dot(C.T, C)) 
    81         tmp = zip([float(val) for val in w], range(self.n)) 
    82         tmp.sort() 
    83         w1, w2 = tmp[-1][0], tmp[-2][0] 
    84         v1, v2 = v[:, tmp[-1][1]], v[:, tmp[-2][1]] 
    85         x = dot(C, v1) 
    86         y = dot(C, v2) 
    87         return x, y 
    88          
    89          
    90 class MDS(object): 
    91     def __init__(self, distances=None, dim=2, **kwargs): 
    92         self.mds=orangemds.MDS(distances, dim, **kwargs) 
    93         self.originalDistances=orange.SymMatrix([m for m in self.distances]) 
    94  
    95     def __getattr__(self, name): 
    96         if name in ["points", "projectedDistances", "distances" ,"stress", "progressCallback", "n", "dim", "avgStress"]: 
    97             #print "rec:",name             
    98             return self.__dict__["mds"].__dict__[name] 
    99         else: 
    100             raise AttributeError(name) 
    101  
    102     def __setattr__(self, name, value): 
    103         #print "setattr" 
    104         if name=="points": 
    105             for i in range(len(value)): 
    106                 for j in range(len(value[i])): 
    107                     self.mds.points[i][j]=value[i][j] 
    108             return 
    109              
    110         if name in ["projectedDistances", "distances" ,"stress", "progressCallback"]: 
    111             self.mds.__setattr__(name, value) 
    112         else: 
    113             self.__dict__[name]=value 
    114              
    115     def __nonzero__(self): 
    116         return True 
    117              
    118     def SMACOFstep(self): 
    119         self.mds.SMACOFstep() 
    120  
    121     def getDistance(self): 
    122         self.mds.getDistance() 
    123  
    124     def getStress(self, stressFunc=SgnRelStress): 
    125         self.mds.getStress(stressFunc) 
    126  
    127     def optimize(self, iter, stressFunc=SgnRelStress, eps=1e-3, progressCallback=None): 
    128         self.mds.progressCallback=progressCallback 
    129         self.mds.optimize(iter, stressFunc, eps) 
    130  
    131     def run(self, iter, stressFunc=SgnRelStress, eps=1e-3, progressCallback=None): 
    132         self.optimize(iter, stressFunc, eps, progressCallback) 
    133  
    134     def Torgerson(self): 
    135         # Torgerson's initial approximation 
    136         O = array([m for m in self.distances]) 
    137          
    138 ##        #B = matrixmultiply(O,O) 
    139 ##        # bug!? B = O**2 
    140 ##        B = dot(O,O) 
    141 ##        # double-center B 
    142 ##        cavg = sum(B, axis=0)/(self.n+0.0)      # column sum 
    143 ##        ravg = sum(B, axis=1)/(self.n+0.0)    # row sum 
    144 ##        tavg = sum(cavg)/(self.n+0.0)   # total sum 
    145 ##        # B[row][column] 
    146 ##        for i in xrange(self.n): 
    147 ##            for j in xrange(self.n): 
    148 ##                B[i,j] += -cavg[j]-ravg[i] 
    149 ##        B = -0.5*(B+tavg) 
    150  
    151         # B = double-center O**2 !!! 
    152         J = identity(self.n) - (1/float(self.n)) 
    153         B = -0.5 * dot(dot(J, O**2), J) 
    154          
    155         # SVD-solve B = ULU' 
    156         #(U,L,V) = singular_value_decomposition(B) 
    157         (U,L,V)=svd(B) 
    158         # X = U(L^0.5) 
    159         # # self.X = matrixmultiply(U,identity(self.n)*sqrt(L)) 
    160         # X is n-dimensional, we take the two dimensions with the largest singular values 
    161         idx = argsort(L)[-self.dim:].tolist() 
    162         idx.reverse() 
    163          
    164         Lt = take(L,idx)   # take those singular values 
    165         Ut = take(U,idx,axis=1) # take those columns that are enabled 
    166         Dt = identity(self.dim)*sqrt(Lt)  # make a diagonal matrix, with squarooted values 
    167         self.points = orange.FloatListList(dot(Ut,Dt)) 
    168         self.freshD = 0 
    169          
    170 #        D = identity(self.n)*sqrt(L)  # make a diagonal matrix, with squarooted values 
    171 #        X = matrixmultiply(U,D) 
    172 #        self.X = take(X,idx,1) 
    173  
    174     # Kruskal's monotone transformation 
    175     def LSMT(self): 
    176         # optimize the distance transformation 
    177         # build vector o 
    178         effect = 0 
    179         self.getDistance() 
    180         o = [] 
    181         for i in xrange(1,self.n): 
    182             for j in xrange(i): 
    183                 o.append((self.originalDistances[i,j],(i,j))) 
    184         o.sort(_mycompare) 
    185         # find the ties in o, and construct the d vector sorting in order within ties 
    186         d = [] 
    187         td = [] 
    188         uv = [] # numbers of consecutively tied o values 
    189         (i,j) = o[0][1] 
    190         distnorm = self.projectedDistances[i,j]*self.projectedDistances[i,j] 
    191         td = [self.projectedDistances[i,j]] # fetch distance 
    192         for l in xrange(1,len(o)): 
    193             # copy now sorted distances in an array 
    194             # but sort distances within a tied o 
    195             (i,j) = o[l][1] 
    196             cd = self.projectedDistances[i,j] 
    197             distnorm += self.projectedDistances[i,j]*self.projectedDistances[i,j] 
    198             if o[l][0] != o[l-1][0]: 
    199                 # differing value, flush 
    200                 sum = reduce(lambda x,y:x+y,td)+0.0 
    201                 d.append([sum,len(td),sum/len(td),td]) 
    202                 td = [] 
    203             td.append(cd) 
    204         sum = reduce(lambda x,y:x+y,td)+0.0 
    205         d.append([sum,len(td),sum/len(td),td]) 
    206         #### 
    207         # keep merging non-monotonous areas in d 
    208         monotony = 0 
    209         while not monotony and len(d) > 1: 
    210             monotony = 1 
    211             pi = 0 # index 
    212             n = 1  # n-areas 
    213             nd = [] 
    214             r = d[0] # current area 
    215             for i in range(1,len(d)): 
    216                 tr = d[i] 
    217                 if r[2]>=tr[2]: 
    218                     monotony = 0 
    219                     effect = 1 
    220                     r[0] += tr[0] 
    221                     r[1] += tr[1] 
    222                     r[2] = tr[0]/tr[1] 
    223                     r[3] += tr[3] 
    224                 else: 
    225                     nd.append(r) 
    226                     r = tr 
    227             nd.append(r) 
    228             d = nd 
    229         # normalizing multiplier 
    230         sum = 0.0 
    231         for i in d: 
    232             sum += i[2]*i[2]*i[1] 
    233         f = sqrt(distnorm/max(sum,1e-6)) 
    234         # transform O 
    235         k = 0 
    236         for i in d: 
    237             for j in range(i[1]): 
    238                 (ii,jj) = o[k][1] 
    239                 self.distances[ii,jj] = f*i[2] 
    240                 k += 1 
    241         assert(len(o) == k) 
    242         self.freshD = 0 
    243         return effect 
    244  
    245 """ 
    246 class WMDS(MDS): 
    247     def setWeights(self, weights): 
    248         self.W=weights 
    249         self.Wsum=sum(sum(self.W)) 
    250  
    251         V = resize(array([0.0]),(self.n,self.n)) 
    252         sumv = array([0.0]*self.n) 
    253         for i in range(1,self.n): 
    254             for j in range(i): 
    255                 t = (-self.W[i][j]-self.W[j][i])/2.0 
    256                 V[i][j] = t 
    257                 V[j][i] = t 
    258                 sumv[i] -= t 
    259                 sumv[j] -= t 
    260         for i in range(self.n): 
    261             V[i][i] = sumv[i] 
    262  
    263         # default for no weights 
    264 #        for i in range(self.n): 
    265 #            for j in range(self.n): 
    266 #                if i == j: 
    267 #                    V[i][j]=self.n-1 
    268 #                else: 
    269 #                    V[i][j]=-1 
    270         T = matrixmultiply(transpose(V),V) 
    271         # compute moore-penrose inverse of T 
    272 #        (U,L,V) = singular_value_decomposition(T) 
    273 #        Lt = identity(self.n) 
    274 #        for i in range(self.n): 
    275 #            if L[i] > 1e-6: 
    276 #                Lt[i][i] /= L[i] 
    277 #            else: 
    278 #                Lt[i][i] = 0.0 
    279 #        MPI = matrixmultiply(transpose(V),matrixmultiply(Lt,transpose(U))) 
    280         self.M = matrixmultiply(generalized_inverse(T),transpose(V)) 
    281         self.iV = generalized_inverse(V) 
    282         #print self.M 
    283  
    284     def getStress(self, stressf=SgnRelStress): 
    285         total=0.0 
    286         for i in range(self.n): 
    287             for j in range(i): 
    288                 r=self.stress[i,j]=stressf(self.distances[i,j], self.projectedDistances[i,j], self.W[i,j]) 
    289                 total+=abs(r) 
    290         self.avgStress=total/(self.n*self.n) 
    291  
    292     def LSMT(self): 
    293         # optimize the distance transformation 
    294         # build vector o 
    295         effect = 0 
    296         self.getDistance() 
    297         o = [] 
    298         for i in xrange(1,self.n): 
    299             for j in xrange(i): 
    300                 # skip the distances we don't care about 
    301                 if self.W[i][j] > 1e-6: 
    302                     o.append((self.originalDistances[i,j],(i,j))) 
    303         o.sort(_mycompare) 
    304         # find the ties in o, and construct the d vector sorting in order within ties 
    305         d = [] 
    306         td = [] 
    307         uv = [] # numbers of consecutively tied o values 
    308         (i,j) = o[0][1] 
    309         distnorm = self.W[i][j]*self.projectedDistances[i,j]*self.projectedDistances[i,j] 
    310         td = [self.W[i][j]*self.projectedDistances[i,j]] # weighted distance 
    311         ld = [self.W[i][j]] # weights 
    312         for l in xrange(1,len(o)): 
    313             # copy now sorted distances in an array 
    314             # but sort distances within a tied o 
    315             (i,j) = o[l][1] 
    316             cd = self.W[i][j]*self.projectedDistances[i,j] 
    317             distnorm += self.W[i][j]*self.projectedDistances[i,j]*self.projectedDistances[i,j] 
    318             if o[l][0] != o[l-1][0]: 
    319                 # differing value, flush 
    320                 sum = reduce(lambda x,y:x+y,td)+0.0 
    321                 lend = 1.0/(reduce(lambda x,y:x+y,ld)+0.0) # sum of weights 
    322                 d.append([sum,lend,sum*lend,td]) 
    323                 td = [] 
    324                 ld = [] 
    325             td.append(cd) 
    326             ld.append(self.W[i][j]) 
    327         lend = 1.0/(reduce(lambda x,y:x+y,ld)+0.0) # sum of weights 
    328         sum = reduce(lambda x,y:x+y,td)+0.0 
    329         d.append([sum,lend,sum*lend,td]) 
    330         #### 
    331         # keep merging non-monotonous areas in d 
    332         monotony = 0 
    333         while not monotony and len(d) > 1: 
    334             monotony = 1 
    335             pi = 0 # index 
    336             n = 1  # n-areas 
    337             nd = [] 
    338             r = d[0] # current area 
    339             for i in xrange(1,len(d)): 
    340                 tr = d[i] 
    341                 if r[2]>=tr[2]: 
    342                     monotony = 0 
    343                     effect = 1 
    344                     r[0] += tr[0]  # sum of weighted distances 
    345                     r[1] += tr[1]  # sum of weights 
    346                     r[2] = tr[0]/tr[1] # average distance 
    347                     r[3] += tr[3]  # array of values 
    348                 else: 
    349                     nd.append(r) 
    350                     r = tr 
    351             nd.append(r) 
    352             d = nd 
    353         # normalizing multiplier 
    354         sum = 0.0 
    355         for i in d: 
    356             sum += i[2]*i[2]*i[1] 
    357         f = sqrt(distnorm/sum) 
    358         # transform O 
    359         k = 0 
    360         for i in d: 
    361             for j in xrange(len(i[3])): # note the use of the length of the array of distances 
    362                 (ii,jj) = o[k][1] 
    363                 self.distances[ii,jj] = f*i[2] 
    364                 k += 1 
    365         assert(len(o) == k) # all values must be transformed 
    366         self.freshD = 0 
    367         return effect 
    368 """ 
    369 if __name__=="__main__": 
    370     data=orange.ExampleTable("doc//datasets//iris.tab") 
    371     dist = orange.ExamplesDistanceConstructor_Euclidean(data) 
    372     matrix = orange.SymMatrix(len(data)) 
    373     matrix.setattr('items', data) 
    374     for i in range(len(data)): 
    375         for j in range(i+1): 
    376             matrix[i, j] = dist(data[i], data[j]) 
    377     mds=MDS(matrix, dim=3) 
    378     mds.Torgerson() 
    379     print mds.points[:3] 
     1from Orange.projection.mds import * 
Note: See TracChangeset for help on using the changeset viewer.