From ed2f6e21511120078e98437ea12bcb420f484c81 Mon Sep 17 00:00:00 2001 From: Yutong Zhao Date: Thu, 8 Oct 2020 10:48:40 -0400 Subject: [PATCH 1/2] Make python3 friendly. --- A_opt.py | 84 +++++++++++++++++++++++++++--------------------------- diffnet.py | 2 +- 2 files changed, 43 insertions(+), 43 deletions(-) diff --git a/A_opt.py b/A_opt.py index f9f191c2..afb849b8 100644 --- a/A_opt.py +++ b/A_opt.py @@ -15,9 +15,9 @@ def solution_to_nij( sol, K, measure_indices=None): CVXOPT solution. ''' if sol['status'] != 'optimal': - raise ValueError, sol['status'] + raise ValueError(sol['status']) x = sol['x'] - # print x + # print(x) n = matrix( 0., (K, K)) for i in xrange(K): if measure_indices is not None: @@ -492,35 +492,35 @@ def kkt_solver( x, y, z): dy = np.max(np.abs(y - yp)) tol = 1e-5 if dx > tol: - print 'dx=' - print dx - print x - print xp + print('dx=') + print(dx) + print(x) + print(xp) if dy > tol: - print 'dy=' - print dy - print y - print yp + print('dy=') + print(dy) + print(y) + print(yp) if dz > tol: - print 'dz=' - print dz - print z - print zp + print('dz=') + print(dz) + print(z) + print(zp) if dx > tol or dy > tol or dz > tol: for i, (r, rti) in enumerate( zip(ris, rtis)): - print 'r[%d]=' % i - print r - print 'rti[%d]=' % i - print rti - print 'rti.T*r=' - print rti.T*r + print('r[%d]=' % i) + print(r) + print('rti[%d]=' % i) + print(rti) + print('rti.T*r=') + print(rti.T*r) for i, d in enumerate( ds): - print 'd[%d]=%g' % (i, d) - print 'x0, y0, z0=' - print x0 - print y0 - print z0 - print Bm0 + print('d[%d]=%g' % (i, d)) + print('x0, y0, z0=') + print(x0) + print(y0) + print(z0) + print(Bm0) ### # END of kkt_solver. @@ -1064,10 +1064,10 @@ def my_solver( W): dx, dy, dz = np.max(np.abs(dx)), np.max(np.abs(dy)), np.max(np.abs(dz)) if tol < np.max( [dx, dy, dz]): - print 'KKT solver FAILS: max(dx=%g, dy=%g, dz=%g) > tol = %g' % \ - (dx, dy, dz, tol) + print('KKT solver FAILS: max(dx=%g, dy=%g, dz=%g) > tol = %g' % \ + (dx, dy, dz, tol)) success = False - print 'KKT solver succeeds: dx=%g, dy=%g, dz=%g' % (dx, dy, dz) + print('KKT solver succeeds: dx=%g, dy=%g, dz=%g' % (dx, dy, dz)) return success def test_Gfunc( ntrials=10, tol=1e-10): @@ -1095,9 +1095,9 @@ def test_Gfunc( ntrials=10, tol=1e-10): dy = np.max(np.abs(y - yp)) if (dy > tol): success = False - print 'G function FAILS for trans=N: dy=%g' % dy + print('G function FAILS for trans=N: dy=%g' % dy) else: - print 'G function succeeds for trans=N: dy=%g' % dy + print('G function succeeds for trans=N: dy=%g' % dy) trans = 'T' nx = K*(K+1)/2 + K*(K+1)*(K+1) @@ -1118,9 +1118,9 @@ def test_Gfunc( ntrials=10, tol=1e-10): dy = np.max(np.abs(y - yp)) if (dy > tol): success = False - print 'G function FAILS for trans=T: dy=%g' % dy + print('G function FAILS for trans=T: dy=%g' % dy) else: - print 'G function succeeds for trans=T: dy=%g' % dy + print('G function succeeds for trans=T: dy=%g' % dy) return success @@ -1149,11 +1149,11 @@ def test_sumdR2( ntrials=10, tol=1e-9): delta = np.max(np.abs(ddR2 - ddR2p)) if (delta > tol): success = False - print 'sum dR test FAILED: delta=%g > tol=%g' % (delta, tol) + print('sum dR test FAILED: delta=%g > tol=%g' % (delta, tol)) else: - print 'sum dR test succeeds: delta=%g' % delta - print 'Timing for naive sum dR: %f seconds per call.' % (tnaive/ntrials) - print 'Timing for aligned sum dR: %f seconds per call.' % (tfast/ntrials) + print('sum dR test succeeds: delta=%g' % delta) + print('Timing for naive sum dR: %f seconds per call.' % (tnaive/ntrials)) + print('Timing for aligned sum dR: %f seconds per call.' % (tfast/ntrials)) return success def test_relative_only(): @@ -1163,7 +1163,7 @@ def test_relative_only(): sij = 0.5*(sij + sij.T) for i in range(K): sij[i,i] = np.inf nij = A_optimize_fast( sij) - print nij + print(nij) def unit_test(): test_Gfunc( ntrials=100) @@ -1195,8 +1195,8 @@ def unit_test(): nij = A_optimize_fast( sij, N, nsofar, only_include_measurements) tend = time.time() tlapse = tend - tstart - print 'Fast A-optimize took %g seconds.' % tlapse - # print nij + print('Fast A-optimize took %g seconds.' % tlapse) + # print(nij) if (K>=80): import sys @@ -1206,6 +1206,6 @@ def unit_test(): nij0 = update_A_optimal_sdp( sij, N, nsofar, only_include_measurements) tend = time.time() tlapse = tend - tstart - print 'SDP A-optimize took %g seconds.' % tlapse + print('SDP A-optimize took %g seconds.' % tlapse) - print 'dn=', np.max(np.abs( nij0 - nij)) + print('dn=', np.max(np.abs( nij0 - nij))) diff --git a/diffnet.py b/diffnet.py index 2c2137b5..93d2925e 100644 --- a/diffnet.py +++ b/diffnet.py @@ -173,7 +173,7 @@ def A_optimize( sij, nadd=1., nsofar=None, delta=None, only_include_measurements) else: if delta is not None: - raise ValueError, 'Currently delta values are only supported in A-optimal by the conelp method.' + raise ValueError('Currently delta values are only supported in A-optimal by the conelp method.') if nsofar is None: nij = A_optimize_sdp( sij) nij *= nadd From cf1e76397eced2141a3575c5a54e79405a4be5b6 Mon Sep 17 00:00:00 2001 From: Yutong Zhao Date: Thu, 8 Oct 2020 12:05:32 -0400 Subject: [PATCH 2/2] Replace xrange with range --- A_opt.py | 144 ++++++++++++++++++++++++------------------------ README.md | 4 +- analysis.ipynb | 44 +++++++-------- diffnet.py | 70 +++++++++++------------ examples.py | 22 ++++---- graph.py | 4 +- netbfe.ipynb | 2 +- netbfe.py | 8 +-- test_diffnet.py | 18 +++--- 9 files changed, 158 insertions(+), 158 deletions(-) diff --git a/A_opt.py b/A_opt.py index afb849b8..b741f79c 100644 --- a/A_opt.py +++ b/A_opt.py @@ -19,14 +19,14 @@ def solution_to_nij( sol, K, measure_indices=None): x = sol['x'] # print(x) n = matrix( 0., (K, K)) - for i in xrange(K): + for i in range(K): if measure_indices is not None: m = measure_indices.get( (i,i), None) if m is not None: n[i,i] = x[m] else: n[i,i] = x[i] - for j in xrange(i+1, K): + for j in range(i+1, K): if measure_indices is not None: m = measure_indices.get( (i,j), None) if m is None: continue @@ -55,8 +55,8 @@ def conelp_solution_to_nij( x, K): ''' n = matrix( 0.0, (K, K)) p = 0 - for i in xrange(K): - for j in xrange(i, K): + for i in range(K): + for j in range(i, K): n[i,j] = n[j,i] = x[p] p += 1 return n @@ -94,11 +94,11 @@ def pairwise_diff( x, y, n): ''' k = 0 r = x.size[0] - for i in xrange(n): + for i in range(n): #y[:,k] = x[:,i] blas.copy( x, y, n=r, offsetx=i*r, offsety=k*r) k+=1 - for j in xrange(i+1, n): + for j in range(i+1, n): #y[:,k] = x[:,i] - x[:,j] blas.copy( x, y, n=r, offsetx=i*r, offsety=k*r) blas.axpy( x, y, alpha=-1, n=r, offsetx=j*r, offsety=k*r) @@ -119,11 +119,11 @@ def congruence_matrix( r, W, offset=0): ''' K = r.size[0] ij = 0 - for i in xrange( K): - for j in xrange( K): + for i in range( K): + for j in range( K): ab = 0 - for a in xrange( K): - for b in xrange( K): + for a in range( K): + for b in range( K): W[offset+ij, offset+ab] = r[a,i]*r[b,j] ab += 1 ij += 1 @@ -134,8 +134,8 @@ def symmetrize_matrix( x, n, xstart=0): Given a sequence x representing the lower triangle of a symmetric square nxn matrix, fill in the upper triangle with symmetric values. ''' - for a in xrange(n): - for b in xrange(a+1, n): + for a in range(n): + for b in range(a+1, n): x[xstart+b*n+a] = x[xstart+a*n+b] def tri2symm( x, n): @@ -146,7 +146,7 @@ def tri2symm( x, n): ''' y = matrix( 0., (n, n)) p = 0 - for i in xrange( n): + for i in range( n): y[i:,i] = x[p:p+n-i] y[i, i:] = y[i:, i].T p += (n-i) @@ -202,7 +202,7 @@ def sumdR2_aligned( Ris, K): # RK[:,1] R2[:,2] ... RK[:,K] KR = matrix( 0., (K*K, K)) start = 0 - for i in xrange(K): + for i in range(K): KR[start:start+K, :] = Ris[i][:K,:K] start += K @@ -245,13 +245,13 @@ def sumdR2_aligned( Ris, K): ddR2 = matrix( 0., (K*(K+1)/2, K*(K+1)/2)) start = 0 bsize = K*(K+1)/2 - for i in xrange(K): + for i in range(K): ddR2 += ddRs[i::K, :] return ddR2 def sumdR2( Ris, K): ddR2 = matrix( 0., (K*(K+1)/2, K*(K+1)/2)) - for i in xrange(K): + for i in range(K): Ri = Ris[i] # dR[:(a,b)] = R[:,a] - R[:,b] dR = matrix( 0., (K, K*(K+1)/2)) @@ -288,8 +288,8 @@ def Aopt_KKT_solver( si2, W): di2s = dis**2 # R_i = r_i^{-t}r_i^{-1} - Ris = [ matrix(0.0, (K+1, K+1)) for i in xrange(K) ] - for i in xrange(K): + Ris = [ matrix(0.0, (K+1, K+1)) for i in range(K) ] + for i in range(K): blas.gemm( rtis[i], rtis[i], Ris[i], transB = 'T') ddR2 = sumdR2( Ris, K) @@ -297,7 +297,7 @@ def Aopt_KKT_solver( si2, W): # upper triangular representation si2ab[(a,b)] := si2[a,b] si2ab = matrix( 0., (K*(K+1)/2, 1)) p = 0 - for i in xrange(K): + for i in range(K): si2ab[p:p+(K-i)] = si2[i:,i] p += (K-i) @@ -321,8 +321,8 @@ def Aopt_KKT_solver( si2, W): Bm[:K*(K+1)/2,:K*(K+1)/2] = cvxopt.mul( si2q, ddR2) row = 0 - for a in xrange(K): - for b in xrange(a, K): + for a in range(K): + for b in range(a, K): Bm[row, row] += di2s[row] # d_{ab}^{-2} n_{ab} row += 1 assert(K*(K+1)/2 == row) @@ -333,19 +333,19 @@ def Aopt_KKT_solver( si2, W): # g_i^t F g_i + R_{i,K+1,K+1}^2 u_i = pi - L_{i,K+1,K+1} dg = matrix( 0., (K, K*(K+1)/2)) g = matrix( 0., (K, K)) - for i in xrange(K): + for i in range(K): g[i,:] = Ris[i][K,:K] # dg[:,(a,b)] = g[a] - g[b] if a!=b else g[a] pairwise_diff( g, dg, K) dg2 = dg**2 # dg2 := s[(a,b)]^{-2} dg[(a,b)]^2 - for i in xrange( K): + for i in range( K): dg2[i,:] = cvxopt.mul( si2ab.T, dg2[i,:]) Bm[K*(K+1)/2:K*(K+1)/2+K,:-K] = dg2 # Diagonal coefficients for u_i. uoffset = K*(K+1)/2 - for i in xrange(K): + for i in range(K): RiKK = Ris[i][K,K] Bm[uoffset+i,uoffset+i] = RiKK**2 @@ -397,7 +397,7 @@ def kkt_solver( x, y, z): zp[:] = z[:] default_solver( xp, yp, zp) offset = K*(K+1)/2 - for i in xrange(K): + for i in range(K): symmetrize_matrix( zp, K+1, offset) offset += (K+1)*(K+1) @@ -432,10 +432,10 @@ def kkt_solver( x, y, z): dL = matrix( 0., (K*(K+1)/2, 1)) ab = 0 - for a in xrange(K): + for a in range(K): dL[ab] = Ls[a,a] ab += 1 - for b in xrange(a+1,K): + for b in range(a+1,K): dL[ab] = Ls[a,a] + Ls[b,b] - 2*Ls[b,a] ab += 1 @@ -470,7 +470,7 @@ def kkt_solver( x, y, z): F = Fisher_matrix( si2, nab) offset = K*(K+1)/2 - for i in xrange( K): + for i in range( K): start, end = i*(K+1)*(K+1), (i+1)*(K+1)*(K+1) Fu = matrix( 0.0, (K+1, K+1)) Fu[:K,:K] = F @@ -484,7 +484,7 @@ def kkt_solver( x, y, z): if (TEST_KKT): offset = K*(K+1)/2 - for i in xrange(K): + for i in range(K): symmetrize_matrix( z, K+1, offset) offset += (K+1)*(K+1) dz = np.max(np.abs(z - zp)) @@ -567,7 +567,7 @@ def Aopt_Gfunc( si2, x, y, alpha=1.0, beta=0.0, trans='N'): Fu = matrix( 0., (K+1, K+1)) Fu[:K,:K] = F start = hkkp1 - for i in xrange(K): + for i in range(K): Fu[K,K] = u[i] y[start:start+kp1sqr] = -Fu[:] + beta*y[start:start+kp1sqr] start += kp1sqr @@ -576,16 +576,16 @@ def Aopt_Gfunc( si2, x, y, alpha=1.0, beta=0.0, trans='N'): xKK = alpha*x[hkkp1+kp1sqr-1::kp1sqr] start = hkkp1 xsum = matrix( 0., (kp1sqr, 1)) - for i in xrange(K): + for i in range(K): xsum += x[start:start+kp1sqr] start += kp1sqr xsum *= alpha xsum = matrix( xsum, (K+1, K+1)) ab = 0 - for a in xrange(K): + for a in range(K): xab[ab] += si2[a,a]*xsum[a,a] ab += 1 - for b in xrange(a+1, K): + for b in range(a+1, K): xab[ab] += si2[b,a]*(xsum[a,a] + xsum[b,b] - 2*xsum[b,a]) ab += 1 y[:hkkp1] = -xab + beta*y[:hkkp1] @@ -625,35 +625,35 @@ def G( x, y, alpha=1., beta=0., trans='N'): Gs = [] # Gs = np.zeros( nrows*ncols) # -I_{K*(K+1)/2} identity matrix - Gs.extend( [ ( i, i, -1. ) for i in xrange(K*(K+1)/2) ]) + Gs.extend( [ ( i, i, -1. ) for i in range(K*(K+1)/2) ]) # Gs[:K*(K+1)/2*nrows:nrows+1] = -1. skip = K*(K+1)/2 # vec( [ V_{ij}, 0; 0, 0 ]) col = 0 - for i in xrange(K): + for i in range(K): # Skip the first K*(K+1)/2 rows Gs.extend( [ (skip+i*(K+2)+t*(K+1)*(K+1), col, -si2[i,i]) - for t in xrange( K) ]) + for t in range( K) ]) # offset = col*nrows + K*(K+1)/2 # Gs[offset+i*(K+2) : (col+1)*nrows : (K+1)*(K+1)] = si2[i,i] col += 1 - for j in xrange(i+1, K): + for j in range(i+1, K): Gs.extend( [ (skip+i*(K+2) + t*(K+1)*(K+1), col, -si2[i,j]) - for t in xrange(K) ]) + for t in range(K) ]) Gs.extend( [ (skip+j*(K+2) + t*(K+1)*(K+1), col, -si2[i,j]) - for t in xrange(K) ]) + for t in range(K) ]) Gs.extend( [ (skip+i*(K+1) + j + t*(K+1)*(K+1), col, si2[i,j]) - for t in xrange(K) ]) + for t in range(K) ]) Gs.extend( [ (skip+j*(K+1) + i + t*(K+1)*(K+1), col, si2[i,j]) - for t in xrange(K) ]) + for t in range(K) ]) col += 1 # vec( [ 0, 0; 0, 1 ]) Gs.extend( [ (skip+(i+1)*(K+1)*(K+1)-1, K*(K+1)/2+i, -1.) - for i in xrange(K) ]) - I, J, X = [ [ ijx[p] for ijx in Gs ] for p in xrange(3) ] + for i in range(K) ]) + I, J, X = [ [ ijx[p] for ijx in Gs ] for p in range(3) ] G = spmatrix( X, I, J, (nrows, ncols)) # h vector. @@ -675,8 +675,8 @@ def G( x, y, alpha=1., beta=0., trans='N'): F += omega row = K*(K+1)/2 - for i in xrange(K): - for j in xrange(K): + for i in range(K): + for j in range(K): h[row + j*(K+1) : row + (j+1)*(K+1)-1] = F[:,j] h[row + (i+1)*(K+1)-1] = 1. # e_i^t h[row + K*(K+1) + i] = 1. # e_i @@ -737,13 +737,13 @@ def A_optimize_fast( sij, N=1., nsofar=None, delta=None, if delta is not None: di2 = np.array( [ 1./delta[i]**2 if delta[i] is not None else 0. - for i in xrange(K)]) + for i in range(K)]) else: di2 = None if only_include_measurements is not None: - for i in xrange(K): - for j in xrange(i, K): + for i in range(K): + for j in range(i, K): if not (i,j) in only_include_measurements: # Set the s[i,j] to infinity, thus excluding the pair. si2[i,j] = si2[j,i] = 0. @@ -817,16 +817,16 @@ def A_optimize_sdp( sij): # G matrix, of dimension ((K+1)*(K+1), (M+K)). Each column is a # column-major vector representing the KxK matrix of U_m augmented # by a length K vector, hence the dimension (K+1)x(K+1). - # Gs = [ matrix( 0., ((K+1)*(K+1), (M+K))) for k in xrange( K) ] + # Gs = [ matrix( 0., ((K+1)*(K+1), (M+K))) for k in range( K) ] G0 = [] - hs = [ matrix( 0., (K+1, K+1)) for k in xrange( K) ] + hs = [ matrix( 0., (K+1, K+1)) for k in range( K) ] - for i in xrange( K): + for i in range( K): # The index of matrix element (i,i) in column-major representation # of a (K+1)x(K+1) matrix is i*(K+1 + 1) # Gs[0][i*(K+2), i] = 1./(sij[i,i]*sij[i,i]) G0.append( (i*(K+2), i, -1./(sij[i,i]*sij[i,i]))) - for j in xrange( i+1, K): + for j in range( i+1, K): m = measurement_index( i, j, K) # The index of matrix element (i,j) in column-major representation # of a (K+1)x(K+1) matrix is j*(K+1) + i @@ -842,7 +842,7 @@ def A_optimize_sdp( sij): # Gs[0] *= -1. Gs = [] - for k in xrange( K): + for k in range( K): # if (k>0): Gs[k][:,:M] = Gs[0][:,:M] # for the term u_k [ [0, 0], [0, 1] ] # Gs[k][-1, M+k] = -1. @@ -929,11 +929,11 @@ def update_A_optimal_sdp( sij, nadd, nsofar, only_include_measurements=None): # G matrix, of dimension ((K+1)*(K+1), (M+K)). Each column is a # column-major vector representing the KxK matrix of U_m augmented # by a length K vector, hence the dimension (K+1)x(K+1). - # Gs = [ matrix( 0., ((K+1)*(K+1), (M+K))) for k in xrange( K) ] + # Gs = [ matrix( 0., ((K+1)*(K+1), (M+K))) for k in range( K) ] G0 = [] - hs = [ matrix( 0., (K+1, K+1)) for k in xrange( K) ] + hs = [ matrix( 0., (K+1, K+1)) for k in range( K) ] - for i in xrange( K): + for i in range( K): # The index of matrix element (i,i) in column-major representation # of a (K+1)x(K+1) matrix is i*(K+1 + 1) v2 = 1./(sij[i,i]*sij[i,i]) @@ -946,7 +946,7 @@ def update_A_optimal_sdp( sij, nadd, nsofar, only_include_measurements=None): # Gs[0][i*(K+2), i] = v2 G0.append( (i*(K+2), i, -v2)) hs[0][i,i] += nsofar[i,i]*v2 - for j in xrange( i+1, K): + for j in range( i+1, K): # The index of matrix element (i,j) in column-major representation # of a (K+1)x(K+1) matrix is j*(K+1) + i v2 = 1./(sij[i,j]*sij[i,j]) @@ -970,7 +970,7 @@ def update_A_optimal_sdp( sij, nadd, nsofar, only_include_measurements=None): # Gs[0] *= -1. Gs = [] - for k in xrange( K): + for k in range( K): if (k>0): # Gs[k][:,:M] = Gs[0][:,:M] hs[k][:K,:K] = hs[0][:K,:K] @@ -1017,17 +1017,17 @@ def my_solver( W): return Aopt_KKT_solver( si2, W) success = True - for t in xrange( ntrials): + for t in range( ntrials): x = matrix( 1*(np.random.rand( K*(K+1)/2+K) - 0.5), (K*(K+1)/2+K, 1)) y = matrix( np.random.rand( 1), (1,1)) z = matrix( 0.0, (K*(K+1)/2 + K*(K+1)*(K+1), 1)) z[:K*(K+1)/2] = 5.*(np.random.rand( K*(K+1)/2) - 0.5) offset = K*(K+1)/2 - for i in xrange(K): + for i in range(K): r = 10*(np.random.rand( (K+1)*(K+2)/2) - 0.3) p = 0 - for a in xrange(K+1): - for b in xrange(a, K+1): + for a in range(K+1): + for b in range(a, K+1): z[offset + a*(K+1) + b] = r[p] z[offset + b*(K+1) + a] = r[p] p+=1 @@ -1035,7 +1035,7 @@ def my_solver( W): ds = matrix( 10*np.random.rand( K*(K+1)/2), (K*(K+1)/2, 1)) rs = [ matrix(np.random.rand( (K+1)*(K+1)) - 0.3, (K+1, K+1)) - for i in xrange(K) ] + for i in range(K) ] W = dict( d=ds, di=cvxopt.div(1., ds), r=rs, @@ -1055,7 +1055,7 @@ def my_solver( W): dx = xp - x dy = yp - y offset = K*(K+1)/2 - for i in xrange(K): + for i in range(K): symmetrize_matrix( zp, K+1, offset) symmetrize_matrix( z, K+1, offset) offset += (K+1)*(K+1) @@ -1081,7 +1081,7 @@ def test_Gfunc( ntrials=10, tol=1e-10): success = True - for i in xrange( ntrials): + for i in range( ntrials): trans = 'N' nx = K*(K+1)/2+K ny = K*(K+1)/2+K*(K+1)*(K+1) @@ -1105,10 +1105,10 @@ def test_Gfunc( ntrials=10, tol=1e-10): x = matrix( np.random.rand( nx), (nx, 1)) y = matrix( 1.e6*np.random.rand( ny), (ny, 1)) - for i in xrange(K): + for i in range(K): start = K*(K+1)/2 + i*(K+1)*(K+1) - for a in xrange(K+1): - for b in xrange(a+1, K+1): + for a in range(K+1): + for b in range(a+1, K+1): x[start+a*(K+1)+b] = x[start+b*(K+1)+a] yp = y[:] @@ -1131,9 +1131,9 @@ def test_sumdR2( ntrials=10, tol=1e-9): tnaive = tfast = 0. success = True - for t in xrange(ntrials): - Ris = [ matrix(np.random.rand(K*K), (K,K)) for i in xrange(K) ] - for i in xrange(K): + for t in range(ntrials): + Ris = [ matrix(np.random.rand(K*K), (K,K)) for i in range(K) ] + for i in range(K): Ris[i] = 0.5*(Ris[i].T + Ris[i]) tstart = time.time() @@ -1181,7 +1181,7 @@ def unit_test(): if (False): connectivity = 5 only_include_measurements = set() - for i in xrange( K): + for i in range( K): js = i + np.floor((K-i)*np.random.rand(connectivity)).astype('int') for j in js: only_include_measurements.add( (i,j)) diff --git a/README.md b/README.md index 0fb91dce..b3e418fe 100644 --- a/README.md +++ b/README.md @@ -72,10 +72,10 @@ while not converged: # nij = A_optimize( sij) # call this if not using iterative optimization. # nij[nij < ncut] = 0 # Omit pairs with impractically small allocations. nij = round_to_integers( nij) - for i in xrange(nmols): + for i in range(nmols): if nij[i,i] > 0: __individual_free_energy__( mols[i], nsamples=nij[i,i], ...) - for j in xrange(i+1, nmols): + for j in range(i+1, nmols): if nij[i,j] == 0: continue # Compute the relative free energy between i and j, # using nij[i,j] number of samples. diff --git a/analysis.ipynb b/analysis.ipynb index 67924cf3..c6c96187 100644 --- a/analysis.ipynb +++ b/analysis.ipynb @@ -186,8 +186,8 @@ "# K0 = 5\n", "# x0 = np.arange( 1., K0+1, 1)\n", "sij0 = np.diag( x0)\n", - "for i in xrange(K0):\n", - " for j in xrange(i+1, K0):\n", + "for i in range(K0):\n", + " for j in range(i+1, K0):\n", " sij0[i,j] = sij0[j,i] = x0[j] - x0[i]\n", "sij0 = matrix( sij0) " ] @@ -267,7 +267,7 @@ " u = np.zeros( K)\n", " u[0] = x0[0]/np.sqrt(K)\n", " s = np.sqrt(K)*x0[0]\n", - " for i in xrange( 1, K):\n", + " for i in range( 1, K):\n", " u[i] = u[i-1] + (x0[i] - x0[i-1])/np.sqrt(K-i)\n", " s += (x0[i] - x0[i-1])*np.sqrt(K-i)\n", " return u*s" @@ -284,7 +284,7 @@ "def distnet_minTrC( xs):\n", " K = len(xs)\n", " trC = np.sqrt(K)*xs[0]\n", - " for i in xrange( 1, K):\n", + " for i in range( 1, K):\n", " trC += np.sqrt(K-i)*(xs[i] - xs[i-1])\n", " return trC**2" ] @@ -765,7 +765,7 @@ " else:\n", " sCOX2rel[i,i] = sCOX2[a, ROF]\n", " origins[i] = 'R'\n", - " for j in xrange(i+1, len(allmols)):\n", + " for j in range(i+1, len(allmols)):\n", " b = allmols[j]\n", " sCOX2rel[i,j] = sCOX2rel[j,i] = sCOX2[a,b]\n", " return matrix(sCOX2rel), origins\n", @@ -1285,7 +1285,7 @@ " stats2 = results[val][o2]\n", " ratio = stats2/stats1\n", " bl = len(ratio)/blocks\n", - " bavg = [ np.mean( ratio[b*bl:(b+1)*bl]) for b in xrange(blocks)]\n", + " bavg = [ np.mean( ratio[b*bl:(b+1)*bl]) for b in range(blocks)]\n", " return np.mean(ratio), np.std( bavg)/np.sqrt(blocks)" ] }, @@ -1829,7 +1829,7 @@ } ], "source": [ - "for j in xrange(4): \n", + "for j in range(4): \n", " plt.plot( x0[j::4], xML[j::4], 'o')\n", "plt.gca().set_aspect( 1)\n", "plt.gca().set_xlabel( r'$x_0$')\n", @@ -1945,9 +1945,9 @@ " '''\n", " K = len(xtrue)\n", " xijtrue = np.zeros( (K, K), dtype=float)\n", - " for i in xrange( K):\n", + " for i in range( K):\n", " xijtrue[i,i] = xtrue[i]\n", - " for j in xrange( i+1, K):\n", + " for j in range( i+1, K):\n", " xijtrue[i,j] = xtrue[i] - xtrue[j]\n", " xijtrue[j,i] = -xijtrue[i,j]\n", " xij = np.zeros( (K, K), dtype=float)\n", @@ -1958,8 +1958,8 @@ " for t, nsofar in enumerate( ntotals):\n", " nactual = np.asarray(fij*(nsofar - nprev), dtype=int)\n", " invsij2 = np.zeros( (K, K), dtype=float)\n", - " for i in xrange(K):\n", - " for j in xrange(i, K):\n", + " for i in range(K):\n", + " for j in range(i, K):\n", " if nactual[i,j] > 0:\n", " gij = np.random.normal( xijtrue[i,j], strue[i,j], nactual[i,j])\n", " xijnew = np.mean( gij)\n", @@ -1994,9 +1994,9 @@ " np.random.seed( rseed)\n", " K = len(xtrue)\n", " xijtrue = np.zeros( (K, K), dtype=float)\n", - " for i in xrange( K):\n", + " for i in range( K):\n", " xijtrue[i,i] = xtrue[i]\n", - " for j in xrange( i+1, K):\n", + " for j in range( i+1, K):\n", " xijtrue[i,j] = xtrue[i] - xtrue[j]\n", " xijtrue[j,i] = -xijtrue[i,j]\n", "\n", @@ -2020,8 +2020,8 @@ " \n", " smin, smax = np.inf, 0\n", " invsij2 = np.zeros( (K, K), dtype=float)\n", - " for i in xrange( K):\n", - " for j in xrange( i, K):\n", + " for i in range( K):\n", + " for j in range( i, K):\n", " if nactual[i,j] > 0:\n", " gij = np.random.normal( xijtrue[i,j], strue[i,j], nactual[i,j])\n", " xijnew = np.mean( gij)\n", @@ -2041,8 +2041,8 @@ " if sij[i,j] < smin: smin = sij[i,j]\n", " \n", " # For pairs where a difference has not been estimated, we set sij to be a random number between smin and smax\n", - " for i in xrange( K):\n", - " for j in xrange( i, K):\n", + " for i in range( K):\n", + " for j in range( i, K):\n", " if nij[i,j] <= 1: # Need at 2 samples to have meaningful sij estimate\n", " z = np.random.rand( 1)\n", " sij[i,j] = smin*(1 - z) + smax*z\n", @@ -3285,7 +3285,7 @@ "nijCOX2trj = []\n", "dGCOX2MSTtrj = []\n", "dGCOX2Aopttrj = []\n", - "for t in xrange( TRIALS):\n", + "for t in range( TRIALS):\n", " _dgtrj, _sijtrj, _trCtrj, _nijtrj = iterative_diffnet( dGCOX2, sCOX2, N=10000, rseed=2019 + t)\n", " dGCOX2trj.append( _dgtrj)\n", " sCOX2trj.append( _sijtrj)\n", @@ -3314,7 +3314,7 @@ } ], "source": [ - "[ dn.sum_upper_triangle( matrix(nijCOX2trj[i][-1])) for i in xrange(TRIALS) ]" + "[ dn.sum_upper_triangle( matrix(nijCOX2trj[i][-1])) for i in range(TRIALS) ]" ] }, { @@ -3351,8 +3351,8 @@ " for t, nij in enumerate( nijs):\n", " nij = nij/dn.sum_upper_triangle( matrix(nij))\n", " d = 0.\n", - " for i in xrange( K):\n", - " for j in xrange( i, K):\n", + " for i in range( K):\n", + " for j in range( i, K):\n", " if nij[i,j]>0: \n", " d += nij[i,j]*np.log( nij[i,j]/(nopt[i,j] + epsilon))\n", " D[t] = d\n", @@ -3503,7 +3503,7 @@ " sstd = np.zeros( len(ns))\n", " for k, n in enumerate(ns):\n", " n = int(n)\n", - " sigma = np.array([ np.std( x[b*n:(b+1)*n]) for b in xrange( blocks) ])\n", + " sigma = np.array([ np.std( x[b*n:(b+1)*n]) for b in range( blocks) ])\n", " smean[k] = np.mean(sigma)\n", " sstd[k] = np.std( sigma)\n", " return ns, smean, sstd\n", diff --git a/diffnet.py b/diffnet.py index 93d2925e..ce0a3905 100644 --- a/diffnet.py +++ b/diffnet.py @@ -50,8 +50,8 @@ def sum_upper_triangle( x): ''' if not isinstance(x, matrix): x = matrix( x) s = 0. - for i in xrange( x.size[0]): - for j in xrange( i, x.size[1]): + for i in range( x.size[0]): + for j in range( i, x.size[1]): s += x[i,j] return s @@ -78,10 +78,10 @@ def lndetC( sij, x, hessian=False): K = sij.size[0] M = K*(K+1)/2 F = matrix( 0., (K, K)) - for i in xrange( K): + for i in range( K): # n_{ii}*v_{ii}.v_{ii}^t F[i,i] += x[i]/(sij[i,i]*sij[i,i]) - for j in xrange( i+1, K): + for j in range( i+1, K): m = measurement_index( i, j, K) v2 = x[m]/(sij[i,j]*sij[i,j]) F[i,i] += v2 @@ -90,32 +90,32 @@ def lndetC( sij, x, hessian=False): C = linalg.inv( F) fval = -np.log(linalg.det( F)) df = matrix( 0., (1, M)) - for i in xrange( K): + for i in range( K): df[i] = -C[i,i]/(sij[i,i]*sij[i,i]) - for j in xrange( i+1, K): + for j in range( i+1, K): m = measurement_index( i, j, K) df[m] = (2*C[i,j] - C[i,i] - C[j,j])/(sij[i,j]*sij[i,j]) if not hessian: return (fval, df) # Compute the Hessian d2f = matrix( 0., (M, M)) - for i in xrange( K): - for j in xrange( i, K): + for i in range( K): + for j in range( i, K): # d^2/dx_i dx_j = C_{ij}^2/(s_{ii}^2 s_{jj}^2) d2f[i, j] = C[i,j]*C[i,j]/(sij[i,i]*sij[i,i]*sij[j,j]*sij[j,j]) d2f[j, i] = d2f[i, j] - for i2 in xrange( K): - for j2 in xrange( i2+1, K): + for i2 in range( K): + for j2 in range( i2+1, K): m2 = measurement_index( i2, j2, K) # d^2/dx_id_x(i',j') = (C_{ii'}-C_{ji'})^2/(s_{i'i'}^2 s_{ij}^2) dC = C[i2,i] - C[j2,i] d2f[i, m2] = dC*dC/(sij[i,i]*sij[i,i]*sij[i2,j2]*sij[i2,j2]) d2f[m2, i] = d2f[i, m2] - for j in xrange( i+1, K): + for j in range( i+1, K): m = measurement_index( i, j, K) invs2 = 1/(sij[i,j]*sij[i,j]) - for i2 in xrange( i, K): - for j2 in xrange( i2+1, K): + for i2 in range( i, K): + for j2 in range( i2+1, K): m2 = measurement_index( i2, j2, K) # d^2/dx_{ij}dx_{i'j'} = # (C_{ii'}+C_{jj'}-C_{ji'}-C_{ij'})^2/(s_{i'j'}^2 s_{ij}^2) @@ -212,9 +212,9 @@ def D_optimize( sij): def F( x=None, z=None): if x is None: - x0 = matrix( [ sij[i,i] for i in xrange( K) ] + - [ sij[i,j] for i in xrange( K) - for j in xrange( i+1, K) ], (M, 1)) + x0 = matrix( [ sij[i,i] for i in range( K) ] + + [ sij[i,j] for i in range( K) + for j in range( i+1, K) ], (M, 1)) return (0, x0) if z is None: return lndetC( sij, x) @@ -246,8 +246,8 @@ def constant_relative_error( si): K = len(si) si = np.sort( si) sij = np.diag( si) - for i in xrange( K): - for j in xrange( i+1, K): + for i in range( K): + for j in range( i+1, K): sij[i,j] = sij[j,i] = si[j] - si[i] return matrix( sij) @@ -260,7 +260,7 @@ def A_optimize_const_relative_error( si): nij = np.zeros( (K, K), dtype=float) N = nij[0,0] = np.sqrt( K)*si[0] - for i in xrange(K-1): + for i in range(K-1): nij[i+1, i] = nij[i, i+1] = np.sqrt(K - (i+1))*(si[i+1] - si[i]) N += nij[i, i+1] @@ -278,7 +278,7 @@ def D_optimize_const_relative_error( si): iK = 1./K nij = np.zeros( (K, K), dtype=float) nij[0,0] = iK - for i in xrange(K-1): + for i in range(K-1): nij[i,i+1] = nij[i+1,i] = iK return matrix( nij) @@ -323,10 +323,10 @@ def E_optimize( sij): # G[i*K + j] = (v_m.v_m^t)[i,j] G = matrix( 0., (K*K, M+1)) h = matrix( 0., (K, K)) - for i in xrange( K): + for i in range( K): G[i*(K+1), i] = 1./(sij[i,i]*sij[i,i]) G[i*(K+1), M] = 1. # The column-major identity matrix for t. - for j in xrange( i+1, K): + for j in range( i+1, K): m = measurement_index( i, j, K) v2 = 1./(sij[i,j]*sij[i,j]) G[j*K + i, m] = G[i*K + j, m] = -v2 @@ -366,7 +366,7 @@ def Dijkstra_shortest_path( sij): while q: d, u = heapq.heappop( q) - for v in xrange( K): + for v in range( K): suv = sij[u,v] if u != K else sij[v,v] dp = d + suv if dp < dist[v]: @@ -407,7 +407,7 @@ def E_optimal_tree( sij): # For each node, compute \sum_{j \elem T_i} a_j, where j runs over # the set of nodes in the subtree T_i rooted at i, including i itself. suma = np.zeros( K, dtype=float) - for v in xrange( K): + for v in range( K): suma[v] += a[v] u = prev[v] # Follow up the tree until the root at the origin @@ -417,7 +417,7 @@ def E_optimal_tree( sij): # n_{i\mu_i} = s_{i\mu_i} \lambda \sum_{j\elem T_i} a_j nij = matrix( 0., (K, K)) - for i in xrange( K): + for i in range( K): j = prev[i] if j!=K: # not the origin nij[i,j] = sij[i,j]*suma[i] @@ -513,7 +513,7 @@ def MLestimate( xij, invsij2, x0=None, di2=None): # z_i = \sigma_i^{-2} x_i + \sum_{j\neq i} \sigma_{ij}^{-2} x_{ij} z = np.sum( invsij2*xij, axis=1) if di2 is not None and x0 is not None: - for i in xrange( len(z)): + for i in range( len(z)): if x0[i] is not None and x0[i] is not np.nan and di2[i] != 0: z[i] += di2[i]*x0[i] @@ -574,7 +574,7 @@ def covariance( invv, delta=None): if np.all( np.diag(f) == 0) and delta is None: return covariance_singular_Fisher( f) F = np.diag( np.sum(f, axis=1) ) - for k in xrange(K): f[k,k] = 0 + for k in range(K): f[k,k] = 0 F -= f if delta is not None: di2 = np.square( 1/delta) @@ -632,15 +632,15 @@ def round_to_integers( n): nceilsum = 0 edges = [] nint = np.zeros( (K, K), dtype=int) - for i in xrange(K): - for j in xrange(i, K): + for i in range(K): + for j in range(i, K): nsum += n[i,j] nint[i,j] = nint[j,i] = int(np.ceil( n[i,j])) nceilsum += nint[i,j] heapq.heappush( edges, (-n[i,j], (i,j))) nsum = int(np.floor(nsum + 0.5)) surplus = nceilsum - nsum - for k in xrange( surplus): + for k in range( surplus): d, (i, j) = heapq.heappop( edges) nint[i,j] -= 1 nint[j,i] = nint[i,j] @@ -747,9 +747,9 @@ def weight( i,j): G.add_node( 'O') edges = [] - for i in xrange(K): + for i in range(K): edges.append( ('O', i, weight( i,i))) - for j in xrange(i+1, K): + for j in range(i+1, K): edges.append( (i, j, weight(i,j))) edges = list(nx.k_edge_augmentation( G, k=connectivity, partial=True)) @@ -771,13 +771,13 @@ def weight( i,j): # network. if (len(only_include_measurements) < n_measure): indices = [] - for i in xrange(K): - for j in xrange(i, K): + for i in range(K): + for j in range(i, K): if (i,j) in only_include_measurements: continue heapq.heappush( indices, (weight(i,j), (i,j))) addition = [] - for m in xrange(n_measure - len(only_include_measurements)): + for m in range(n_measure - len(only_include_measurements)): _w, (i,j) = heapq.heappop( indices) addition.append( (i,j)) only_include_measurements.update( addition) diff --git a/examples.py b/examples.py index 2ddf76ae..418bd2ac 100644 --- a/examples.py +++ b/examples.py @@ -28,14 +28,14 @@ def distance_net( K, dim=2, rmin=0): x = 2.*(np.random.rand( K, dim) - 0.5) dij = matrix( 0., (K, K)) - for i in xrange( K): + for i in range( K): dij[i,i] = np.sqrt(x[i].dot(x[i])) if rmin>0 and dij[i,i] 0: # results for absolute binding free energy. dG[i,i] = dG0[i] dG[i,i] += 2*np.sqrt(1/isigma2[i,i])*(np.random.rand() - 0.5) - for j in xrange(i+1, K): + for j in range(i+1, K): if n[i,j] > 0: ddG = dG0[i] - dG0[j] ddG += 2*np.sqrt(1./isigma2[i,j])*(np.random.rand() - 0.5) @@ -153,7 +153,7 @@ def test_A_optimality_with_reference( s, n, delta, dn=1E-1, ntimes=10): cov = covariance( cvxopt.div( n, s**2), delta) f = np.trace( cov) df = np.zeros( ntimes) - for t in xrange( ntimes): + for t in range( ntimes): zeta = matrix( 1. + 2*dn*(np.random.rand( K, K) - 0.5)) n1 = cvxopt.mul( n, zeta) n1 = 0.5*(n1 + n1.trans()) # Symmetrize @@ -185,7 +185,7 @@ def unit_test(): # The experimental values of the molecules not in the references # will be unavailable. if references is not None: - for i in xrange( K): + for i in range( K): if i not in references: dG0p[i] = None delta[i] = np.infty diff --git a/test_diffnet.py b/test_diffnet.py index afdd82db..a44ffc4c 100644 --- a/test_diffnet.py +++ b/test_diffnet.py @@ -22,7 +22,7 @@ def check_optimality( sij, nij, optimality='A', delta=1E-1, ntimes=10): fC['Etree'] = np.max( linalg.eig( C)[0]).real df = np.zeros( ntimes) - for t in xrange( ntimes): + for t in range( ntimes): zeta = matrix( 1. + 2*delta*(np.random.rand( K, K) - 0.5)) nijp = cvxopt.mul( nij, zeta) nijp = 0.5*(nijp + nijp.trans()) # Symmetrize @@ -80,7 +80,7 @@ def check_update_A_optimal( sij, delta=5e-1, ntimes=10, tol=1e-5): trC = np.trace( C) dtr = np.zeros( ntimes) - for t in xrange( ntimes): + for t in range( ntimes): zeta = matrix( 1. + 2*delta*(np.random.rand( K, K) - 0.5)) nnextp = cvxopt.mul( nnext, zeta) nnextp = 0.5*(nnextp + nnextp.trans()) @@ -134,7 +134,7 @@ def check_sparse_A_optimal( sij, ntimes=10, delta=1e-1, tol=1e-5): trC = np.trace( covariance( cvxopt.div( nij, sij**2))) dtr = np.zeros( ntimes) - for t in xrange( ntimes): + for t in range( ntimes): zeta = matrix( 1. + 2*delta*(np.random.rand( K, K) - 0.5)) nijp = cvxopt.mul( nij, zeta) nijp = 0.5*(nijp + nijp.trans()) # Symmetrize @@ -180,7 +180,7 @@ def check_hessian( dF, d2F, x0): N = len(x0) esqr = 0. - for i in xrange( N): + for i in range( N): def func( x): return dF(x)[i] def dfunc( x): @@ -196,21 +196,21 @@ def fabricate_measurements( K=10, sigma=0.1, noerror=True, disconnect=False): invsij2 = 0.5*(invsij2 + np.transpose( invsij2)) sij = np.sqrt( 1./invsij2) if noerror: sij *= 0. - for i in xrange(K): + for i in range(K): xij[i][i] = x0[i] + sij[i,i]*np.random.randn() - for j in xrange(i+1, K): + for j in range(i+1, K): xij[i][j] = x0[i] - x0[j] + sij[i][j]*np.random.randn() xij[j][i] = -xij[i][j] if (disconnect >= 1): # disconnect the origin and thus eliminate the individual measurements - for i in xrange(K): invsij2[i][i] = 0 + for i in range(K): invsij2[i][i] = 0 if (disconnect >= 2): # disconnect the network into the given number of disconnected # components. - for i in xrange( K): + for i in range( K): c1 = i % disconnect - for j in xrange( i+1, K): + for j in range( i+1, K): c2 = j % disconnect if (c1 != c2): invsij2[i][j] = invsij2[j][i] = 0