Source Code for Module SloppyCell.Vandermonde.VdmPairwise

  1  import scipy 
  2   
  3  try: 
  4      import SloppyCell.Plotting as Plotting 
  5  except: 
  6      pass 
  7   
  8  # 1. calculate your jtj matrix 
  9  # 2. use calcDmat to get Dmat 
 10  # 3. use getLambdaR2 to get lambdaMat, r2Mat 
 11  # 4. use calcR2mat_Terms to get ctMat==cross term from r2Mat and 
 12  #    dtMat==diagonal term from r2Mat 
 13  # 5. use cluster to cluster based on either r2Mat (for unscaled distances) or 
 14  #    r2Mat/dtMat (to scale distances by size of individual parameter effects) 
 15   
 16 -def normColumns(mat1): 
 17      n = len(mat1[0]) 
 18      mat2 = mat1.copy() 
 19      for ii in range(n): 
 20          mat2[:,ii] = mat2[:,ii]/scipy.sqrt(scipy.dot(mat2[:,ii],mat2[:,ii])) 
 21      return mat2 
 22   
 23 -def colorAlignColumns(mat1): 
 24      n = len(mat1[0]) 
 25      mat2 = mat1.copy() 
 26      for ii in range(n-1): 
 27          if scipy.dot(mat2[:,ii],mat2[:,ii+1]) < 0.: 
 28              mat2[:,ii+1] = -1.*mat2[:,ii+1] 
 29      return mat2 
 30   
 31 -def calcCovarianceMat(origMat): 
 32      covarMat = scipy.dot(scipy.transpose(origMat),origMat) 
 33      n = len(covarMat) 
 34      jaa = scipy.outer(scipy.diagonal(covarMat),scipy.ones(n,'d')) 
 35      jbb = scipy.outer(scipy.ones(n,'d'),scipy.diagonal(covarMat)) 
 36      covarMat = covarMat/scipy.sqrt(jaa) 
 37      covarMat = covarMat/scipy.sqrt(jbb) 
 38      return covarMat 
 39   
 40 -def calcDmat(jtj): 
 41      n = len(jtj) 
 42      Dmat = scipy.zeros((n,n),'d') 
 43      for alpha in range(n): 
 44          for beta in range(alpha,n): 
 45              Dmat[alpha,beta] = (jtj[alpha,alpha]-jtj[beta,beta])**2 + 4.*(jtj[alpha,beta]**2) 
 46      Dmat = Dmat + scipy.transpose(Dmat) 
 47      for alpha in range(n): 
 48          Dmat[alpha,alpha] = 4.*(jtj[alpha,alpha]**2) 
 49      return Dmat 
 50   
 51 -def calcR2mat_Terms(jtj,lambdaMat): 
 52      n = len(jtj) 
 53      jaa = scipy.outer(scipy.diagonal(jtj),scipy.ones(n,'d')) 
 54      jbb = scipy.outer(scipy.ones(n,'d'),scipy.diagonal(jtj)) 
 55      crosstermMat = 2.*lambdaMat*scipy.sqrt(1.-lambdaMat**2.)*jtj 
 56      diagtermMat = (lambdaMat**2.)*jaa + (1.-lambdaMat**2.)*jbb 
 57      return crosstermMat, diagtermMat     
 58   
 59 -def calcR2mat_3(jtj,Dmat): 
 60  #def calcR2mat_3(jtj,lambdaMat): 
 61      n = len(jtj) 
 62      jaa = scipy.outer(scipy.diagonal(jtj),scipy.ones(n,'d')) 
 63      jbb = scipy.outer(scipy.ones(n,'d'),scipy.diagonal(jtj)) 
 64      r2mat = (1./2.)*(-scipy.sqrt(Dmat)+4.*(jtj**2)/scipy.sqrt(Dmat) + (jaa - 4.*jtj*scipy.sqrt((jtj**2)/Dmat) + jbb)) 
 65  #    r2mat = (lambdaMat**2.)*jaa + (1.-lambdaMat**2.)*jbb + 2.*lambdaMat*scipy.sqrt(1.-lambdaMat**2.)*jtj 
 66      return r2mat 
 67   
 68 -def calcLambdaMat_3(jtj,Dmat): 
 69      n = len(jtj) 
 70      lambdaMat = -(1./scipy.sqrt(2.))*scipy.sqrt(1.+(-scipy.outer(scipy.diagonal(jtj),scipy.ones(n,'d')) +  scipy.outer(scipy.ones(n,'d'),scipy.diagonal(jtj)))/scipy.sqrt(Dmat)) 
 71      return lambdaMat 
 72   
 73 -def calcR2mat_2(jtj,Dmat): 
 74      n = len(jtj) 
 75      jaa = scipy.outer(scipy.diagonal(jtj),scipy.ones(n,'d')) 
 76      jbb = scipy.outer(scipy.ones(n,'d'),scipy.diagonal(jtj)) 
 77      r2mat = (1./2.)*(scipy.sqrt(Dmat)-4.*(jtj**2)/scipy.sqrt(Dmat) + (jaa + 4.*jtj*scipy.sqrt((jtj**2)/Dmat) + jbb)) 
 78      return r2mat 
 79   
 80 -def calcLambdaMat_2(jtj,Dmat): 
 81      n = len(jtj) 
 82      lambdaMat = (1./scipy.sqrt(2.))*scipy.sqrt(1.+(scipy.outer(scipy.diagonal(jtj),scipy.ones(n,'d')) -  scipy.outer(scipy.ones(n,'d'),scipy.diagonal(jtj)))/scipy.sqrt(Dmat)) 
 83      return lambdaMat 
 84   
 85 -def calcR2mat_1(jtj,Dmat): 
 86      n = len(jtj) 
 87      jaa = scipy.outer(scipy.diagonal(jtj),scipy.ones(n,'d')) 
 88      jbb = scipy.outer(scipy.ones(n,'d'),scipy.diagonal(jtj)) 
 89      r2mat = (1./2.)*(scipy.sqrt(Dmat)-4.*(jtj**2)/scipy.sqrt(Dmat) + (jaa - 4.*jtj*scipy.sqrt((jtj**2)/Dmat) + jbb)) 
 90      return r2mat 
 91   
 92 -def calcLambdaMat_1(jtj,Dmat): 
 93      n = len(jtj) 
 94      lambdaMat = -(1./scipy.sqrt(2.))*scipy.sqrt(1.+(scipy.outer(scipy.diagonal(jtj),scipy.ones(n,'d')) -  scipy.outer(scipy.ones(n,'d'),scipy.diagonal(jtj)))/scipy.sqrt(Dmat)) 
 95      return lambdaMat 
 96   
 97 -def calcR2mat_4(jtj,Dmat): 
 98      n = len(jtj) 
 99      jaa = scipy.outer(scipy.diagonal(jtj),scipy.ones(n,'d')) 
100      jbb = scipy.outer(scipy.ones(n,'d'),scipy.diagonal(jtj)) 
101      r2mat = (1./2.)*(-scipy.sqrt(Dmat)+4.*(jtj**2)/scipy.sqrt(Dmat) + (jaa + 4.*jtj*scipy.sqrt((jtj**2)/Dmat) + jbb)) 
102      return r2mat 
103   
104 -def calcLambdaMat_4(jtj,Dmat): 
105      n = len(jtj) 
106      lambdaMat = (1./scipy.sqrt(2.))*scipy.sqrt(1.+(-scipy.outer(scipy.diagonal(jtj),scipy.ones(n,'d')) +  scipy.outer(scipy.ones(n,'d'),scipy.diagonal(jtj)))/scipy.sqrt(Dmat)) 
107      return lambdaMat 
108   
109 -def checkLambda(jtj,lMat,r2Mat,a,b): 
110      lambdas = scipy.arange(-1.,1.,0.001) 
111      residuals = [l*l*jtj[a,a]+(1.-l*l)*jtj[b,b]+2.*l*scipy.sqrt(1.-l*l)*jtj[a,b] for l in lambdas] 
112      Plotting.plot(lambdas,residuals) 
113      Plotting.plot([lMat[a,b]],[r2Mat[a,b]],'o') 
114      Plotting.title(str(a)+', '+str(b)) 
115      return 
116   
117 -def getMinIndex(listofNums): 
118      minInd = 0 
119      minVal = listofNums[minInd] 
120      for ii in range(1,len(listofNums)): 
121          if listofNums[ii] < minVal: 
122              minInd = ii 
123              minVal = listofNums[minInd] 
124      return minInd 
125   
126 -def getMinLambdaR2(lambdaMats, res2Mats): 
127      numMats = len(lambdaMats) 
128      numParams = len(lambdaMats[0]) 
129      minR2 = scipy.zeros((numParams,numParams),'d') 
130      minLambda = scipy.zeros((numParams,numParams),'d') 
131      for ii in range(numParams): 
132          for jj in range(numParams): 
133              minIndex = getMinIndex([res2Mats[k][ii,jj] for k in range(numMats)]) 
134              minR2[ii,jj] = res2Mats[minIndex][ii,jj] 
135              minLambda[ii,jj] = lambdaMats[minIndex][ii,jj] 
136      return minLambda, minR2 
137   
138 -def getLambdaR2(jtj,Dmat): 
139      lambdaMats = [calcLambdaMat_1(jtj,Dmat),calcLambdaMat_2(jtj,Dmat), calcLambdaMat_3(jtj,Dmat),calcLambdaMat_4(jtj,Dmat)] 
140      r2Mats = [calcR2mat_1(jtj,Dmat),calcR2mat_2(jtj,Dmat),calcR2mat_3(jtj,Dmat),calcR2mat_4(jtj,Dmat)] 
141      lambdaMat, R2Mat = getMinLambdaR2(lambdaMats, r2Mats) 
142      return lambdaMat, R2Mat 
143