qubx.model_constraints

1 """System of linear constraints for constrained optimization. 2 3 There are more informative comments in the source code. 4 5 Copyright 2008-2012 Research Foundation State University of New York 6 This file is part of QUB Express. 7 8 QUB Express is free software; you can redistribute it and/or modify 9 it under the terms of the GNU General Public License as published by 10 the Free Software Foundation, either version 3 of the License, or 11 (at your option) any later version. 12 13 QUB Express is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 GNU General Public License for more details. 17 18 You should have received a copy of the GNU General Public License, 19 named LICENSE.txt, in the QUB Express program directory. If not, see 20 <http://www.gnu.org/licenses/>. 21 22 """ 23 24 from numpy import * 25 import scipy 26 import scipy.linalg 27 28 # What are the optimization parameters? (some of) the off-diagonal elements of K0 and K1. 29 # Actually, we use log(K0) to keep them positive. 30 # Breaking with past practice, we include all the unconnected zeros, 31 # and put all the K0s before all the K1s. As we'll see, this is a column vector. 32

33 -def K_index(N):

34 """Returns all pairs (a,b) such that (0 <= a,b < N) and (a != b); in canonical order.""" 35 ix = [] 36 if N <= 0: return ix 37 for a in xrange(N): 38 for b in xrange(N): 39 if a != b: 40 ix.append((a,b)) 41 return ix

42

43 -def K01toK(K0, K1):

44 """ 45 Returns a flat array of the useful (non-diagonal) numbers in log(K0) and K1; 46 first all the K0s, then all the K1s, in canonical (L{K_index}) order. 47 """ 48 k0, k1 = [], [] 49 for a,b in K_index(K0.shape[1]): 50 if a != b: 51 x = K0[a,b] 52 k0.append((x > 0.0) and log(x) or 0.0) # tricky: 1->0 and 0->0; see KtoK01 53 k1.append(K1[a,b]) 54 K = matrix(zeros(shape=(len(k0)+len(k1), 1))) 55 K[:,0] = matrix(array(k0 + k1)).T 56 return K

57

58 -def K012toK(K0, K1, K2):

59 """ 60 Returns a flat array of the useful (non-diagonal) numbers in log(K0) and K1; 61 first all the K0s, then all the K1s, in canonical (L{K_index}) order. 62 """ 63 k0, k1, k2 = [], [], [] 64 for a,b in K_index(K0.shape[1]): 65 if a != b: 66 x = K0[a,b] 67 k0.append((x > 0.0) and log(x) or 0.0) # tricky: 1->0 and 0->0; see KtoK01 68 k1.append(K1[a,b]) 69 k2.append(K2[a,b]) 70 K = matrix(zeros(shape=(len(k0)+len(k1)+len(k2), 1))) 71 K[:,0] = matrix(array(k0 + k1 + k2)).T 72 return K

73 74 # K alone isn't enough to reconstruct K0 and K1, since 0 could mean log(1), or else unconnected. 75 # We disambiguate by checking the original K0 for zeros, which mean unconnected. 76

77 -def KtoK01(K, last_K0):

78 """ 79 Returns matrices K0, and K1 filled with numbers from K (see L{K01toK}). 80 You must provide last_K0 to disambiguate between K[i] == 0 == log(1) vs. K0[i] == 0 == unconnected. 81 """ 82 # len(K)/2 = K0.N * (K0.N - 1); K0.N = (1 + sqrt(1 + 2*len(K))) / 2.0 83 N = int(round((1 + sqrt(1 + 2*len(K))) / 2.0)) 84 ixs = K_index(N) 85 K0 = matrix(zeros(shape=(N,N))) 86 i = 0 87 for a,b in ixs: 88 if a != b: 89 if last_K0[a,b]: 90 K0[a,b] = exp(K[i,0]) 91 i += 1 92 K1 = matrix(zeros(shape=(N,N))) 93 for a,b in ixs: 94 if a != b: 95 K1[a,b] = K[i] 96 i += 1 97 return K0, K1

98

99 -def KtoK012(K, last_K0):

100 """ 101 Returns matrices K0, K1, K2 filled with numbers from K (see L{K01toK}). 102 You must provide last_K0 to disambiguate between K[i] == 0 == log(1) vs. K0[i] == 0 == unconnected. 103 """ 104 N = last_K0.shape[0] 105 ixs = K_index(N) 106 K0 = matrix(zeros(shape=(N,N))) 107 i = 0 108 for a,b in ixs: 109 if a != b: 110 if last_K0[a,b]: 111 K0[a,b] = exp(K[i,0]) 112 i += 1 113 K1 = matrix(zeros(shape=(N,N))) 114 for a,b in ixs: 115 if a != b: 116 K1[a,b] = K[i] 117 i += 1 118 K2 = matrix(zeros(shape=(N,N))) 119 for a,b in ixs: 120 if a != b: 121 K2[a,b] = K[i] 122 i += 1 123 return K0, K1, K2

124 125 126 # But some of those K are zero and have to stay zero, and the user may have other constraints. 127 # We represent the constraints by the system of linear equations (11) and derive the transformations (12, 13). 128 129 # Some constraints are more a part of the representation than the model: 130 # * where K0 == 0 there's no connection; fix k0 and k1 131 # * where V == 0 fix k1

132 -def AutoConstraints(K0, K1, K2, L, V, P):

133 """Returns the constraints (Ain, Bin) [11] which fix unconnected k0 and unsensitive k1.""" 134 aa = [] 135 bb = [] 136 ixs = K_index(K0.shape[0]) 137 for i, xx in enumerate(ixs): 138 a, b = xx 139 if not K0[a,b]: 140 coeffs = zeros(shape=(1,3*len(ixs))) 141 coeffs[0, i] = 1.0 142 aa.append(coeffs) 143 bb.append(0.0) 144 V[a,b] = 0 145 if not V[a,b]: 146 coeffs = zeros(shape=(1,3*len(ixs))) 147 coeffs[0, len(ixs)+i] = 1.0 148 aa.append(coeffs) 149 bb.append(0.0) 150 if not P[a,b]: 151 coeffs = zeros(shape=(1,3*len(ixs))) 152 coeffs[0, 2*len(ixs)+i] = 1.0 153 aa.append(coeffs) 154 bb.append(0.0) 155 if aa: 156 return vstack(aa), matrix(array(bb)).T 157 else: 158 return matrix(zeros(shape=(0,0))), matrix(zeros(shape=(0,1)))

159 160 # The user might also have some constraints for our linear system:

161 -def UserConstraints(K0, K1, K2, L, V, P, constraints):

162 """ 163 Returns the constraints (Ain, Bin) [11] 164 corresponding to the tuples (name, states) in constraints 165 where name is in ["FixRate", "FixExp", "FixPress", "ScaleRate", "ScaleExp", "ScalePress", "BalanceLoop", "ImbalanceLoop"] 166 and states is a sequence of state indices. The first four names require 2,2,4,4 states respectively; 167 the last two require at least 3.""" 168 aa = [] 169 bb = [] 170 ixs = K_index(K0.shape[0]) 171 ixo = dict([(ab, i) for i, ab in enumerate(ixs)]) 172 N = len(ixs) 173 for name, states in constraints: 174 if name == 'FixRate' and K0[states[0], states[1]]: 175 r = ixo[(states[0], states[1])] 176 coeffs = zeros(shape=(1,3*N)) 177 coeffs[0, r] = 1.0 178 aa.append(coeffs) 179 bb.append(log(K0[states[0], states[1]])) 180 if name == 'FixExp' and K1[states[0], states[1]]: 181 r = ixo[(states[0], states[1])] 182 coeffs = zeros(shape=(1,3*N)) 183 coeffs[0, N+r] = 1.0 184 aa.append(coeffs) 185 bb.append(K1[states[0], states[1]]) 186 if name == 'FixPress' and K2[states[0], states[1]]: 187 r = ixo[(states[0], states[1])] 188 coeffs = zeros(shape=(1,3*N)) 189 coeffs[0, 2*N+r] = 1.0 190 aa.append(coeffs) 191 bb.append(K2[states[0], states[1]]) 192 if (name == 'ScaleRate') and K0[states[0], states[1]] and K0[states[2], states[3]]: 193 r1 = ixo[(states[0], states[1])] 194 r2 = ixo[(states[2], states[3])] 195 coeffs = zeros(shape=(1,3*N)) 196 coeffs[0, r1] = 1.0 197 coeffs[0, r2] = -1.0 198 aa.append(coeffs) 199 bb.append(log(K0[states[0], states[1]]) - log(K0[states[2], states[3]])) 200 if (name == 'ScaleExp') and K1[states[0], states[1]] and K1[states[2], states[3]]: 201 r1 = ixo[(states[0], states[1])] 202 r2 = ixo[(states[2], states[3])] 203 coeffs = zeros(shape=(1,3*N)) 204 coeffs[0, N+r1] = 1.0 205 coeffs[0, N+r2] = - K1[states[0], states[1]] / K1[states[2], states[3]] 206 aa.append(coeffs) 207 bb.append(0.0) 208 if (name == 'ScalePress') and K2[states[0], states[1]] and K2[states[2], states[3]]: 209 r1 = ixo[(states[0], states[1])] 210 r2 = ixo[(states[2], states[3])] 211 coeffs = zeros(shape=(1,3*N)) 212 coeffs[0, 2*N+r1] = 1.0 213 coeffs[0, 2*N+r2] = - K2[states[0], states[1]] / K2[states[2], states[3]] 214 aa.append(coeffs) 215 bb.append(0.0) 216 if name == 'LoopBal': 217 v_depend = p_depend = False 218 loopst = list(states) + [states[0]] 219 coeffs = zeros(shape=(1,3*N)) 220 for i in xrange(len(loopst)-1): 221 rf, rb = ixo[(loopst[i], loopst[i+1])], ixo[(loopst[i+1], loopst[i])] 222 coeffs[0,rf] = 1.0 223 coeffs[0,rb] = -1.0 224 v_depend = v_depend or V[(loopst[i], loopst[i+1])] or V[(loopst[i+1], loopst[i])] 225 p_depend = p_depend or P[(loopst[i], loopst[i+1])] or P[(loopst[i+1], loopst[i])] 226 aa.append(coeffs) 227 bb.append(0.0) 228 if v_depend: 229 coeffs = zeros(shape=(1,3*N)) 230 for i in xrange(len(loopst)-1): 231 rf, rb = ixo[(loopst[i], loopst[i+1])], ixo[(loopst[i+1], loopst[i])] 232 coeffs[0,N+rf] = 1.0 233 coeffs[0,N+rb] = -1.0 234 aa.append(coeffs) 235 bb.append(0.0) 236 if p_depend: 237 coeffs = zeros(shape=(1,3*N)) 238 for i in xrange(len(loopst)-1): 239 rf, rb = ixo[(loopst[i], loopst[i+1])], ixo[(loopst[i+1], loopst[i])] 240 coeffs[0,2*N+rf] = 1.0 241 coeffs[0,2*N+rb] = -1.0 242 aa.append(coeffs) 243 bb.append(0.0) 244 elif name == 'LoopImbal': 245 v_depend = p_depend = False 246 loopst = list(states) + [states[0]] 247 coeffs = zeros(shape=(1,3*N)) 248 b = 0.0 249 for i in xrange(len(loopst)-1): 250 rf, rb = ixo[(loopst[i], loopst[i+1])], ixo[(loopst[i+1], loopst[i])] 251 coeffs[0,rf] = 1.0 252 coeffs[0,rb] = -1.0 253 b += log(K0[(loopst[i], loopst[i+1])]) - log(K0[(loopst[i+1], loopst[i])]) 254 v_depend = v_depend or V[(loopst[i], loopst[i+1])] or V[(loopst[i+1], loopst[i])] 255 p_depend = p_depend or P[(loopst[i], loopst[i+1])] or P[(loopst[i+1], loopst[i])] 256 aa.append(coeffs) 257 bb.append(b) 258 if v_depend: 259 coeffs = zeros(shape=(1,3*N)) 260 for i in xrange(len(loopst)-1): 261 rf, rb = ixo[(loopst[i], loopst[i+1])], ixo[(loopst[i+1], loopst[i])] 262 coeffs[0,N+rf] = 1.0 263 coeffs[0,N+rb] = -1.0 264 aa.append(coeffs) 265 bb.append(0.0) 266 if p_depend: 267 coeffs = zeros(shape=(1,3*N)) 268 for i in xrange(len(loopst)-1): 269 rf, rb = ixo[(loopst[i], loopst[i+1])], ixo[(loopst[i+1], loopst[i])] 270 coeffs[0,2*N+rf] = 1.0 271 coeffs[0,2*N+rb] = -1.0 272 aa.append(coeffs) 273 bb.append(0.0) 274 275 if aa: 276 return vstack(aa), matrix(array(bb)).T 277 else: 278 return matrix(zeros(shape=(0,N))), matrix(zeros(shape=(0,1)))

279 280 281 # You could think of others too; just vstack them with the rest:

282 -def AllConstraints(K0, K1, K2, L, V, P, constraints):

283 """Returns combined AutoConstraints(...) and UserConstraints(...)""" 284 A, B = UserConstraints(K0, K1, K2, L, V, P, constraints) 285 Aauto, Bauto = AutoConstraints(K0, K1, K2, L, V, P) 286 if A.shape[0] == 0: 287 return Aauto, Bauto 288 elif Aauto.shape[0] == 0: 289 return A, B 290 else: 291 return vstack((A, Aauto)), vstack((B, Bauto))

292 293 # To derive the transformation we will use singular value decomposition (svd) 294 # which is really well described at http://www.uwlax.edu/faculty/will/svd/index.html 295 296 # This wrapper around numpy.linalg.svd makes it behave more like LAPACK: 297 # * V is transposed 298 # * matrices have full dimension 299 # Then we repair and check it: 300 # * tiny numbers are replaced with 0 301 # * U and V need orthonormal columns 302 # And derive some useful extras: 303 # * Winv, where inv(0) = 0 304 # * Ainv = V*Winv*U.T 305 # * Nsv = number of nonzero singular values

306 -def adjusted_svd(A):

307 """Returns U, W, V, Winv, Ainv, Nsv, such that A = U*W*V.T = Ainv.I, and A has rank Nsv.""" 308 u,s,v = linalg.svd(A, full_matrices=True) 309 N = s.shape[0] 310 U = matrix(zeros(shape=(max(A.shape[0], u.shape[0]), max(A.shape[1], u.shape[1])))) 311 U[:u.shape[0],:u.shape[1]] = u 312 S = zeros(shape=(U.shape[1],)) 313 S[:N] = s 314 N = len(s[s>=1e-10]) 315 W = diag(S) 316 V = matrix(zeros(shape=(A.shape[1], W.shape[1]))) 317 V[:v.shape[1],:v.shape[0]] = v.T 318 for M in U, W, V: 319 for i in xrange(M.shape[0]): 320 for j in xrange(M.shape[1]): 321 if abs(M[i,j]) < 1e-10: 322 M[i,j] = 0.0 323 for i,x in enumerate(S): 324 if abs(x) < 1e-10: 325 S[i] = 0.0 326 Sinv = array([x and 1/x or 0 for x in S]) 327 Winv = diag(Sinv) 328 # A = U*W*V.T; A.I = V*Winv*U.T 329 if linalg.norm(identity(U.shape[0]) - U*U.T) > 1e-8: 330 raise Exception("Bad SVD: U's columns are not orthonormal") 331 if linalg.norm(identity(V.shape[0]) - V*V.T) > 1e-8: 332 raise Exception("Bad SVD: V's columns are not orthonormal") 333 return U, W, V, Winv, V*(Winv*U.T), N

334 335 # This function cleans up the constraint system prior to the transform, 336 # by deleting duplicate constraints. It uses two checks in case of numerical inaccuracy. 337 # It's no longer called automatically in linear_constraints().

338 -def reduce_constraints(A, B):

339 """Returns constraint system (A, B) with linearly dependent rows deleted.""" 340 Csrc = matrix(hstack((A, B)).T) # may not be square 341 N = max(Csrc.shape) 342 C = matrix(zeros(shape=(N,N))) # (14) 343 C[:Csrc.shape[0],:Csrc.shape[1]] = Csrc 344 Csub = C 345 cols = [C[:,i] for i in xrange(Csrc.shape[1])] 346 i = 0 347 rank = 0 348 while i < len(cols): 349 Csub = hstack(cols[:i+1]) # (15) 350 U, W, V, Winv, Cinv, Nsv = adjusted_svd(Csub) # (16) 351 if Nsv == rank: # check 1: added constraint must increase rank (21) 352 del cols[i] 353 else: 354 rank = Nsv 355 j = i+1 356 while j < len(cols): 357 # check 2: other constraints can't be lin. combinations of the prior 358 image = Cinv * cols[j] # (19) 359 if 1e-5 > linalg.norm(cols[j] - (Csub*image)): 360 del cols[j] 361 else: 362 j += 1 363 i += 1 364 return matrix(Csub[:-1,:].T), matrix(Csub[-1,:].T)

365 366 367 # Now that we have K0, K1, L, V and rate vector K, 368 # we do some magic with the SVD to change your linear constraints 369 # A * K = B 370 # into a transformation between K and the column vector P of free parameters 371 # K = A*P + B 372 # P = Ainv * K 373

374 -def linear_constraints(Ain, Bin, K0):

375 """Returns A, B, A.I, P0 [12].""" 376 A, B = Ain, Bin # reduce_constraints(Ain, Bin) 377 Nc = A.shape[0] 378 Nk = A.shape[1] 379 if (A.shape[0] == 0) or (A.shape[1] == 0): 380 return identity(len(K0)), matrix(zeros(shape=(len(K0),1))), matrix(array([x for x in K0])).T, K0 381 U, W, V, Winv, Ainv, rank = adjusted_svd(A) 382 if rank < Nc: 383 raise Exception('Bad SVD: Some constraints may be incompatible') 384 if rank != Nc: 385 raise Exception('Too many constraints. %i constraints may be linear combinations of the others.' % (Nc - rank)) 386 for i in xrange(Nc): 387 for j in xrange(Nc,Nk): 388 if abs(U[i,j] > 1e-5): 389 raise Exception('Bad SVD: 0 != U[c,non] == %.6g' % U.A[i,j]) 390 Acns = matrix(V[:,Nc:]) # (24) 391 Bcns = Ainv * B # (29) 392 Out0 = Acns.T*K0 # (13) 393 return Acns, Bcns, Acns.T, Out0

394 395 # And yet another linear transformation, this time to scale initial parameter values to 1.0. 396 # (In the case where p[i] == 0, we scale by 1 and translate by -1.) 397 # P = As*S + Bs 398 # S = As.I * (P - Bs) 399

400 -def start_at_ones(P0):

401 """Returns As, Bs, As.I, S0 [32].""" 402 scales = array([x or 1.0 for x in P0.A[:,0]]) 403 reverses = array([x and 1/x or 1.0 for i,x in enumerate(P0.A[:,0])]) 404 Ascal = matrix(diag(scales)) # (35) 405 Bscal = matrix(zeros(shape=P0.shape)) # (36) 406 for i,x in enumerate(P0.A[:,0]): 407 if not x: 408 Bscal[i] = -1.0 409 return Ascal, Bscal, matrix(diag(reverses)), matrix(ones(P0.shape))

410

411 -def KToPars(rates, A, Ainv, B): # (13)

412 k = matrix(zeros(shape=(prod(rates.shape),1))) 413 k[:,0] = matrix(rates.flatten()).T 414 return Ainv*k 415 416 # variant 1: pars is a python or nd- array, such as S provided by an optimizer

417 -def parsToK(pars, A, Ainv, B): # (12, 32)

418 pa = array(pars).flatten() 419 P = matrix(pa).T 420 K = ParsToK(P, A, Ainv, B) 421 return K 422 423 # variant 2: pars is a column vector, such as P = parsToK(S, ...)

424 -def ParsToK(P, A, Ainv, B): # (12, 32)

425 return A*P + B 426 427 428 # To set up, you provide K0,K1,L,V and constraints, and we derive 429 # P and S and their transforms. Then the optimizer adjusts S and we use it 430 # to build the Q matrices: 431 # S --parsToK---> P --ParsToK---> K --KtoK01---> K0,K1 --BuildQ---> Q 432 433 434

435 -def CalcStdOfK(Acns, Ascal, HessInv):

436 """Returns array([error limit of k for k in K]). 437 438 @param Acns: the Nk x Np matrix from L{linear_constraints}; K = Acns*P + B 439 @param Ascal: the Np x Np matrix from L{start_at_ones} 440 @param HessInv: inverse Hessian matrix, as returned by qubx.optimize 441 """ 442 hdiag = matrix(diag(HessInv).reshape((HessInv.shape[0],1))) 443 Asqr = matrix(array(Acns * Ascal)**2) 444 return array(sqrt(Asqr * hdiag)).flatten()

445 446 447 # mil confidence: 448 # dx[i] = sqrt( sum_j( A[i,j]**2 * H[j,j] ) ) { * exp(x[i]) } 449

Source Code for Module qubx.model_constraints