Source Code for Module SloppyCell.Vandermonde.OptimizeSumDets

  1  import scipy 
  2  import time 
  3  try: 
  4      import SloppyCell.Plotting as Plotting 
  5  except ImportError: 
  6      pass 
  7   
  8 -def C(p, origMatrix, weightsRDP=0.,weights2D=0.,weightsLS=0.,weightPR=0.,weightPriors=0.,*args,**kwargs): 
  9      """ 
 10      p is list {a_{ij}, i: 1..n, j: i+1..n} 
 11      n is size of total matrix 
 12      """ 
 13      n = origMatrix.shape[0] 
 14      OrthMatrix = ProcessHalfMatrix(p) 
 15      newMat = transformMatrix(origMatrix,OrthMatrix) 
 16      cost=0. 
 17      if weightsRDP > 0.: 
 18          cost += weightsRDP*sumRowDotProdsOLD(newMat) 
 19      if weights2D > 0.: 
 20          cost += weights2D*sum2Determinants(newMat) 
 21      if weightsLS > 0.: 
 22          cost += weightsLS*sumLogSpacings(newMat) 
 23      if weightPR>0.: 
 24          cost += weightPR/(ParticipationRatio(OrthMatrix)**2.) 
 25  ##         cost += weightPR*(1.-ParticipationRatio(OrthMatrix)/n) 
 26      if weightPriors>0.: 
 27          cost +=weightPR*calcPriors(p) 
 28      return cost 
 29   
 30 -def calcPriors(paramsList, pOpt=0.,pSigma=10.): 
 31      priorCost=0. 
 32      for param in paramsList: 
 33          priorCost += ((param-pOpt)/pSigma)**2. 
 34      return priorCost 
 35   
 36 -def sumRowDotProdsOLD(origMatrix): 
 37      """ 
 38      Makes more sense to use on Jacobian than on Hessian. 
 39      """ 
 40      rowNormalized = normRows(origMatrix) 
 41      n = origMatrix.shape[0] 
 42  ##    sumDotProds = sum([abs((scipy.dot(rowNormalized[i],rowNormalized[i+1])))**(1./2.) for i in range(n-1)]) 
 43      sumDotProds = sum([1.-(scipy.dot(rowNormalized[i],rowNormalized[i+1]))**2. for i in range(n-1)]) 
 44      return sumDotProds 
 45   
 46 -def sumRowDotProdsNEW(origMatrix): 
 47      """ 
 48      Makes more sense to use on Jacobian than on Hessian. 
 49      """ 
 50      rowNormalized = normRows(origMatrix) 
 51      n = rowNormalized.shape[0] 
 52      sumDotProds = sumAllDotProds(rowNormalized[0:n/2]) + sumAllDotProds(rowNormalized[n/2:n]) 
 53      return sumDotProds 
 54   
 55 -def sumAllDotProds(origMatrix): 
 56      n = origMatrix.shape[0] 
 57      sumDotProds=0. 
 58      for ii in range(n): 
 59          for jj in range(ii+1,n): 
 60              sumDotProds += 1.-scipy.dot(origMatrix[ii],origMatrix[jj])**2. 
 61      return sumDotProds 
 62   
 63 -def sum2Determinants(matrixToSum): 
 64      n = matrixToSum.shape[0]/2 
 65      det1 = scipy.linalg.det(matrixToSum[0:n,0:n]) 
 66      det2 = scipy.linalg.det(matrixToSum[n:n*2,n:n*2]) 
 67      return scipy.absolute(det1)+scipy.absolute(det2) 
 68   
 69 -def sumLogSpacings(matrixToSum): 
 70      n = matrixToSum.shape[0]/2 
 71      sv1 = scipy.linalg.svdvals(matrixToSum[0:n,0:n]) 
 72      sv2 = scipy.linalg.svdvals(matrixToSum[n:n*2,n:n*2]) 
 73      return sum(sv1[1:n]/sv1[0:n-1])+sum(sv2[1:n]/sv2[0:n-1]) 
 74   
 75 -def ParticipationRatio(origMatrix): 
 76      return scipy.sum(scipy.sum(origMatrix**4)) 
 77   
 78 -def ProcessFullMatrix(p): 
 79      """ 
 80      this is used if parameters define all elements of a matrix 
 81      """ 
 82      n = scipy.sqrt(len(p)).__int__() 
 83      pMat = scipy.reshape(p,(n,n)) 
 84      orthMat = scipy.linalg.orth(pMat) 
 85      return orthMat 
 86   
 87 -def ProcessHalfMatrix(p): 
 88      """ 
 89      this is used if parameters define upper right triangular portion of 
 90      skew-symmetric matrix. 
 91      Then performs a Cayley transformation of this skew-symmetric matrix. 
 92      """ 
 93      l = scipy.size(p) 
 94      n = scipy.round(((1.+scipy.sqrt(1+8.*l))/2.)).__int__() 
 95      Mfull = scipy.zeros((n,n),'d') 
 96      placeCounter=0 
 97      for i in range(n-1): 
 98          Mfull[i,(i+1):n] = p[placeCounter:(placeCounter+n-i-1)] 
 99          placeCounter += n-i-1 
100      Mfull = Mfull - scipy.transpose(Mfull) 
101      IMat = scipy.eye(n,'d') 
102      return scipy.dot((IMat-Mfull),scipy.linalg.inv(IMat+Mfull)) 
103  ##     return scipy.linalg.expm(Mfull) 
104   
105 -def ProcessQuarterMatrix(p): 
106      """ 
107      this is used if parameters define just upper right and lower left 
108      quadrants of antisymmetric matrix. 
109      """ 
110      n= scipy.sqrt(scipy.size(p)) 
111      Mu = scipy.resize(p, (n,n)) 
112      Mfull = scipy.bmat([[scipy.zeros((n,n)),Mu],[-scipy.transpose(Mu),scipy.zeros((n,n))]]) 
113      return scipy.linalg.expm(Mfull) 
114   
115 -def BestMatrix(origMatrix, p=None, weightsRDP=1., weights2D=0., weightsLS=0., weightPR=0., weightPriors=0., seed=None, *args, **kwargs): 
116      if seed is not None: 
117          scipy.random.seed(seed) 
118      n = origMatrix.shape[0] 
119      if p is None: 
120          p = (scipy.random.random(n*(n-1)/2)-0.5) 
121      pOpt = scipy.optimize.fmin(C, p, args=(origMatrix,weightsRDP, weights2D, weightsLS, weightPR, weightPriors), *args, **kwargs) 
122      return C(pOpt, origMatrix, weightsRDP, weights2D, weightsLS, weightPR, weightPriors), pOpt 
123   
124 -def ManyBest(origMatrix, numberTries=10, p=None,numOpts=10,weightsRDP=1., weights2D=0., weightsLS=0., weightPR=0., weightPriors=0., seed=None, *args, **kwargs): 
125      return [OptimizeManyTimes(origMatrix,p,numOpts,weightsRDP, weights2D, weightsLS, weightPR, weightPriors, seed, *args, **kwargs) for i in range(numberTries)] 
126       
127 -def OptimizeManyTimes(origMatrix,p=None,numOpts=10,weightsRDP=1., weights2D=0., weightsLS=0., weightPR=0., weightPriors=0., seed=None, *args, **kwargs): 
128      C,pOpt = BestMatrix(origMatrix,p,weightsRDP,weights2D,weightsLS,weightPR,weightPriors,seed,*args,**kwargs) 
129      for i in range(numOpts-1): 
130          C,pOpt = BestMatrix(origMatrix,pOpt,weightsRDP,weights2D,weightsLS,weightPR,weightPriors,seed,*args,**kwargs) 
131      return C, pOpt 
132   
133 -def getGammas(numExps,distWidth=1.): 
134      epsilons = scipy.random.random(numExps)-0.5 
135      gammas = 1.+epsilons 
136      return gammas 
137   
138 -def getAmounts(numExps,distWidth=1.): 
139  #    amounts=scipy.ones(numExps,'d') 
140      amounts = scipy.random.random(numExps)-0.5 
141      amounts = 1.+amounts 
142      return amounts 
143   
144 -def netRadiation(gammas,amounts,time): 
145      netR = amounts*scipy.exp(-gammas*time) 
146      return scipy.sum(netR) 
147   
148 -def getHessFull(gammas=None,amounts=None,numExps=3): 
149      if gammas is None: 
150          gammas = getGammas(numExps) 
151      if amounts is None: 
152          amounts = getAmounts(numExps) 
153   
154      numExps=scipy.size(gammas) 
155   
156      hessgg = getHessGG(gammas) 
157   
158      hessAA = getHessAA(amounts,gammas) 
159   
160      numerAg1 = -scipy.outer(amounts,amounts*gammas) 
161      denomAg1 = scipy.outer(scipy.ones(numExps),gammas) 
162      denomAg2 = scipy.transpose(denomAg1)+denomAg1 
163      denomAg3 = denomAg2*denomAg2 
164      hessAg = numerAg1/denomAg3 
165   
166      numergA1 = -scipy.outer(amounts*gammas,amounts) 
167      denomgA1 = scipy.outer(gammas,scipy.ones(numExps)) 
168      denomgA2 = scipy.transpose(denomgA1)+denomgA1 
169      denomgA3 = denomgA2*denomgA2 
170      hessgA = numergA1/denomgA3 
171   
172      hessFull = scipy.bmat([[hessAA,hessAg],[hessgA,hessgg]]) 
173   
174      return hessFull 
175   
176 -def getHessFullLogTime(gammas=None,amounts=None,numExps=3): 
177      """use this routine for new paper""" 
178      if gammas is None: 
179          gammas = getGammas(numExps) 
180      if amounts is None: 
181          amounts = getAmounts(numExps) 
182   
183      numExps=scipy.size(gammas) 
184   
185      amountsLess1 = amounts[0:-1] 
186      gammasLess1 = gammas[0:-1] 
187   
188      hessgg = getHessGGLogTime(gammas=gammas,amounts=amounts) 
189   
190      hessAA = getHessAALogTime(amounts=amounts,gammas=gammas) 
191   
192      numerAg1 = scipy.array(-2.*scipy.outer(amountsLess1,amounts*gammas)) 
193      denomAg1 = scipy.array(scipy.outer(gammasLess1,scipy.ones(numExps))) 
194      denomAg2 = scipy.array(scipy.outer(scipy.ones(numExps-1),gammas)) 
195      denomAg3 = scipy.array(denomAg1+denomAg2) 
196      denomAg4 = denomAg2+gammas[-1] 
197      hessAg = scipy.array(numerAg1/denomAg3)-scipy.array(numerAg1/denomAg4) 
198   
199      hessgA = scipy.transpose(hessAg) 
200   
201      hessFull = scipy.array(scipy.bmat([[hessAA,hessAg],[hessgA,hessgg]])) 
202   
203      return hessFull 
204   
205 -def getHessGG(gammas=None,numExps=3): 
206      if gammas is None: 
207          numExps=6 
208          gammas = scipy.array([1.,1.2,1.4,1001.,1005.,997.]) 
209   
210      numExps = scipy.size(gammas) 
211      numergg1 = 2.*scipy.outer(gammas, gammas) 
212      denomgg1 = scipy.outer(scipy.ones(numExps),gammas) 
213      denomgg2 = scipy.transpose(denomgg1)+denomgg1 
214      denomgg3 = denomgg2*denomgg2*denomgg2 
215      hessgg = numergg1/denomgg3 
216   
217      return hessgg 
218   
219 -def getHessGGLogTime(gammas=None,amounts=None,numExps=3): 
220      """use this routine for the new paper""" 
221      if gammas is None: 
222          numExps=6 
223          gammas = scipy.array([1.,1.2,1.4,1001.,1005.,997.]) 
224      if amounts is None: 
225          amounts = scipy.ones(numExps,'d') 
226   
227      numExps = scipy.size(gammas) 
228      numergg1 = scipy.array(2.*scipy.outer(gammas*amounts, gammas*amounts)) 
229      denomgg1 = scipy.array(scipy.outer(scipy.ones(numExps),gammas)) 
230      denomgg2 = scipy.array(scipy.transpose(denomgg1)+denomgg1) 
231      denomgg3 = scipy.array(denomgg2*denomgg2) 
232      hessgg = scipy.array(numergg1/denomgg3) 
233   
234      return hessgg 
235   
236 -def getHessAA(amounts=None,gammas=None,numExps=3): 
237      if amounts is None: 
238          amounts = getAmounts(numExps) 
239      if gammas is None: 
240          gammas = getGammas(numExps) 
241   
242      numExps = scipy.size(amounts) 
243   
244      numerAA1 = scipy.outer(amounts,amounts) 
245      denomAA1 = scipy.outer(scipy.ones(numExps),gammas) 
246      denomAA2 = scipy.transpose(denomAA1)+denomAA1 
247      hessAA = numerAA1/denomAA2 
248   
249      return hessAA 
250   
251 -def getHessAALogTime(amounts=None,gammas=None,numExps=3): 
252      """use this routine for the new paper""" 
253      if amounts is None: 
254          amounts = getAmounts(numExps) 
255      if gammas is None: 
256          gammas = getGammas(numExps) 
257   
258      numExps = scipy.size(amounts) 
259      amountsLess1 = amounts[0:-1] 
260      gammasLess1 = gammas[0:-1] 
261   
262      numerAA1 = scipy.array(scipy.outer(amountsLess1,amountsLess1)) 
263      numerAA2 = scipy.array(scipy.outer(gammasLess1+gammas[-1],scipy.ones(numExps-1))) 
264      numerAA3 = scipy.array(scipy.outer(scipy.ones(numExps-1),gammasLess1+gammas[-1])) 
265      denomAA1 = scipy.array(scipy.outer(scipy.ones(numExps-1),2*gammas[-1]*gammasLess1)) 
266      denomAA2 = scipy.array(scipy.transpose(denomAA1)+denomAA1) 
267      hessAA = 2.*scipy.array(numerAA1*scipy.log(numerAA2*numerAA3/denomAA2)) 
268   
269      return hessAA 
270   
271 -def getJacobian(times, amounts=None, gammas=None): 
272      numGammas=0 
273      numAmounts=0 
274      if gammas is not None: 
275          numGammas = scipy.size(gammas) 
276      if amounts is not None: 
277          numAmounts = scipy.size(amounts) 
278      Jacobian = scipy.zeros((numGammas+numAmounts,scipy.size(times)),'d') 
279      if numAmounts > 0: 
280          Jacobian[0:numAmounts] = getJacobianA(times,amounts,gammas) 
281      if numGammas > 0: 
282          Jacobian[numAmounts:numAmounts+numGammas] = getJacobianG(times,gammas,amounts) 
283      return Jacobian 
284   
285 -def getJacobianLogTime(times, amounts=None, gammas=None): 
286      numGammas=0 
287      numAmounts=0 
288      if gammas is not None: 
289          numGammas = scipy.size(gammas) 
290      if amounts is not None: 
291          numAmounts = scipy.size(amounts) 
292      JacobianLogTime = scipy.zeros((numGammas+numAmounts,scipy.size(times)),'d') 
293      if numAmounts > 0: 
294          JacobianLogTime[0:numAmounts] = getJacobianALogTime(times,amounts,gammas) 
295      if numGammas > 0: 
296          JacobianLogTime[numAmounts:numAmounts+numGammas] = getJacobianGLogTime(times,gammas,amounts) 
297      return JacobianLogTime 
298   
299 -def getJacobianLog(times,gammas,amounts=None): 
300      """This is the routine to use for new paper (with times exponentially 
301      distributed.""" 
302      numExps = scipy.size(gammas) 
303      if amounts is None: 
304          amounts = scipy.ones(numExps,'d') 
305      jacG = getJacobianLogG(times,gammas,amounts) 
306      jacA = getJacobianLogA(times,gammas,amounts) 
307      return scipy.transpose(scipy.array(scipy.bmat([[jacA],[jacG]]))) 
308       
309   
310 -def getJacobianLogG(times,gammas,amounts=None): 
311      """This is the routine to use for new paper (with times exponentially 
312      distributed.""" 
313      numExps = scipy.size(gammas) 
314      if amounts is None: 
315          amounts = scipy.ones(numExps,'d') 
316      JacobianLogG = scipy.array([gammas[j]*amounts[j]*(-times)*scipy.exp(-gammas[j]*times) for j in range(numExps)]) 
317      return JacobianLogG 
318   
319 -def getJacobianLogA(times,gammas,amounts): 
320      """This is the routine to use for new paper (with times exponentially 
321      distributed.""" 
322      numExps = scipy.size(gammas) 
323      JacobianLogA = scipy.array([amounts[j]*(scipy.exp(-gammas[j]*times)-scipy.exp(-gammas[-1]*times)) for j in range(numExps-1)]) 
324      return JacobianLogA 
325   
326 -def getJacobianG(times,gammas,amounts=None): 
327      numGammas = scipy.size(gammas) 
328      if amounts is None: 
329          amounts = scipy.ones(numGammas,'d') 
330      JacobianG = scipy.array([-amounts[i]*gammas[i]*times*scipy.exp(-gammas[i]*times) for i in range(numGammas)]) 
331      return JacobianG 
332   
333 -def getJacobianGLogTime(times,gammas,amounts=None): 
334      numGammas = scipy.size(gammas) 
335      if amounts is None: 
336          amounts = scipy.ones(numGammas,'d') 
337      JacobianGLogTime = scipy.array([-amounts[i]*gammas[i]*scipy.sqrt(times)*scipy.exp(-gammas[i]*times) for i in range(numGammas)]) 
338      return JacobianGLogTime 
339   
340 -def getJacobianA(times,amounts,gammas=None): 
341      numAmounts = scipy.size(amounts) 
342      if gammas is None: 
343          gammas = scipy.zeros(numAmounts,'d') 
344      JacobianA = scipy.array([amounts[i]*scipy.exp(-gammas[i]*times) for i in range(numAmounts)]) 
345      return JacobianA 
346   
347 -def getJacobianALogTime(times,amounts,gammas=None): 
348      numAmounts = scipy.size(amounts) 
349      if gammas is None: 
350          gammas = scipy.zeros(numAmounts,'d') 
351      JacobianALogTime = scipy.array([(amounts[i]/scipy.sqrt(times))*scipy.exp(-gammas[i]*times) for i in range(numAmounts)]) 
352      return JacobianALogTime 
353   
354 -def normRows(matrixToPlot): 
355      """ 
356      Useful for plotting and otherwise comparing alignment of rows of matrices. 
357      Be careful that if vec[1] is zero, entire row gets zerod. 
358      """ 
359      return scipy.array([scipy.sign(vec[1])*vec/scipy.sqrt(scipy.dot(vec,vec)) for vec in matrixToPlot]) 
360   
361 -def transformMatrix(origMatrix, transformation): 
362      n, m =origMatrix.shape 
363      if n == m: 
364          newMat = scipy.dot(transformation,scipy.dot(origMatrix,scipy.transpose(transformation))) 
365      else: 
366          newMat = scipy.dot(transformation,origMatrix) 
367      return newMat 
368   
369 -def findPermutation(permMat): 
370      """ 
371      In so far as permMat is a permutation matrix, returns the permutation. 
372      """ 
373      maxs = [(vec.tolist().index(max(vec)), max(vec)) for vec in permMat] 
374      mins = [(vec.tolist().index(min(vec)), min(vec)) for vec in permMat] 
375      for ii in range(len(maxs)): 
376          if maxs[ii][1] < -mins[ii][1]: 
377              maxs[ii] = mins[ii] 
378      return maxs 
379   
380 -def makePermutationMatrix(permList): 
381      """ 
382      Takes a list defining the permutation and makes the appropriate matrix. 
383      """ 
384      permList = scipy.array(permList) 
385      n = len(permList) 
386      if 0 not in permList: 
387          permList = permList - 1 
388      permMat = scipy.zeros((n,n),'d') 
389      for ii, jj in enumerate(permList): 
390          permMat[ii,jj] = 1. 
391      return permMat 
392   
393 -def timeCostEval(mat, p, numEvals=100): 
394      start = time.time() 
395      for ii in range(numEvals): 
396          Cost = C(p,mat) 
397      stop = time.time() 
398      print "time per eval:", (stop-start)/numEvals 
399      return (stop-start)/numEvals 
400   
401 -def timeGammas(listNumGammas, numEvals=1000): 
402      timesToCalc = [] 
403      for numGammas in listNumGammas: 
404          gammas = getGammas(numGammas) 
405          jac = getJacobianG(gammas=gammas,times=scipy.arange(500.)) 
406          p = scipy.random.random(numGammas*(numGammas-1)/2)-0.5 
407          timePerEval = timeCostEval(jac,p,numEvals) 
408          timesToCalc.append(timePerEval) 
409      return timesToCalc 
410   
411 -def plotLevels(levels,offset=0): 
412      xs = [offset+0.6,offset+1.4] 
413      for lvl in levels: 
414          Plotting.semilogy(xs,[lvl,lvl],'k') 
415      return 
416   
417 -def eVec(k,n): 
418      eVec = scipy.zeros(n,'d') 
419      eVec[k] = 1. 
420      return eVec 
421   
422 -def sLSAlongP(origMatrix,params,k,deltaP): 
423      sLSList = [sumLogSpacings(transformMatrix(origMatrix,ProcessHalfMatrix(params+x*params[k]*eVec(k,len(params))))) for x in scipy.arange(-deltaP,deltaP,deltaP/1000)] 
424      return sLSList 
425   
426 -def PRAlongP(params,k,deltaP): 
427      n = len(params) 
428      PRList = [1./ParticipationRatio(ProcessHalfMatrix(params+x*params[k]*eVec(k,n))) for x in scipy.arange(-deltaP,deltaP,deltaP/1000)] 
429  ##    PRList = [1./ParticipationRatio(transformMatrix(origMatrix,ProcessHalfMatrix(params+x*params[k]*eVec(k,n)))) for x in scipy.arange(-deltaP,deltaP,deltaP/100)] 
430  ##     PRList = [1.-ParticipationRatio(transformMatrix(origMatrix,ProcessHalfMatrix(params+x*params[k]*eVec(k,n))))/n for x in scipy.arange(-deltaP,deltaP,deltaP/100)] 
431      return PRList 
432   
433 -def CostAlongP(origMatrix,params,k,deltaP,weightsRDP=0.,weights2D=0.,weightsLS=0.,weightPR=0.,weightPriors=0.): 
434      n = len(params) 
435      CostList = [C(params+x*params[k]*eVec(k,n),origMatrix, weightsRDP,weights2D,weightsLS,weightPR,weightPriors) for x in scipy.arange(-deltaP,deltaP,deltaP/100)] 
436      return CostList 
437