Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/block_compose.py: 70%

1from __future__ import division 

2from __future__ import absolute_import 

3 

4from builtins import str 

5from builtins import map 

6from builtins import range 

7from builtins import object 

8"""Classes that define algebraic operations on matrices by deferring the actual 

9action until a vector is present. The classes are normally not used directly, 

10but instead returned implicitly by mathematical notation. For example, 

11 

12M = A*B+C 

13 

14returns the object 

15 

16M = block_add(block_mul(A,B), C), 

17 

18and the actual calculation is done later, once w=M*v is called for a vector v: 

19 

20w=block_add.__mul__(M, v) 

21--> a = block_mul.__mul__((A,B), v) 

22----> b = Matrix.__mul__(B, v) 

23----> a = Matrix.__mul__(A, b) 

24--> b = Matrix.__mul__(C, v) 

25--> w = a+b 

26""" 

27 

28from .block_base import block_base 

29 

30class block_mul(block_base): 

31 def __init__(self, A, B): 

32 """Args may be blockoperators or individual blocks (scalar or Matrix)""" 

33 A = A.chain if isinstance(A, block_mul) else [A] 

34 B = B.chain if isinstance(B, block_mul) else [B] 

35 self.chain = A+B 

36 

37 def __mul__(self, x): 

38 for op in reversed(self.chain): 

39 x = op * x 

40 if isinstance(x, type(NotImplemented)): 40 ↛ 41line 40 didn't jump to line 41, because the condition on line 40 was never true

41 return NotImplemented 

42 return x 

43 

44 def transpmult(self, x): 

45 from .block_util import isscalar 

46 for op in self.chain: 

47 if isscalar(op): 

48 if op != 1: 48 ↛ 46line 48 didn't jump to line 46, because the condition on line 48 was never false

49 x = op*x 

50 else: 

51 x = op.transpmult(x) 

52 return x 

53 

54 def create_vec(self, dim=1): 

55 """Create a vector of the correct size and layout for b (if dim==0) or 

56 x (if dim==1).""" 

57 

58 # Use try-except because even if op has a create_vec method, it may be 

59 # a composed operator and its create_vec method may fail if it cannot 

60 # find a component to create the vector. 

61 # 

62 # The operator may be scalar, in which case it doesn't change the 

63 # dimensions, but it may also be something that reduces to a 

64 # scalar. Hence, a simple is_scalar is not generally sufficient, and we 

65 # just try the next operator in the chain. This is not completely safe, 

66 # in particular it may produce a vector of the wrong size if non-square 

67 # numpy matrices are involved. 

68 

69 if dim==0: 69 ↛ 70line 69 didn't jump to line 70, because the condition on line 69 was never true

70 for op in self.chain: 

71 try: 

72 return op.create_vec(dim) 

73 except AttributeError: 

74 pass 

75 if dim==1: 75 ↛ 83line 75 didn't jump to line 83, because the condition on line 75 was never false

76 for op in reversed(self.chain): 76 ↛ 83line 76 didn't jump to line 83, because the loop on line 76 didn't complete

77 try: 

78 return op.create_vec(dim) 

79 except AttributeError: 

80 pass 

81 # Use AttributeError, because that's what the individual operators use, 

82 # and an outer create_vec will try the next candidate 

83 raise AttributeError('failed to create vec, no appropriate reference matrix') 

84 

85 def block_collapse(self): 

86 """Create a block_mat of block_muls from a block_mul of 

87 block_mats. See block_transform.block_collapse.""" 

88 from .block_mat import block_mat 

89 from .block_util import isscalar 

90 from .block_transform import block_collapse, block_simplify 

91 

92 # Reduce all composed objects 

93 ops = list(map(block_collapse, self.chain)) 

94 

95 # Do the matrix multiply, blockwise. Note that we use 

96 # block_mul(A,B) rather than A*B to avoid any implicit calculations 

97 # (e.g., scalar*matrix->matrix) -- the result will be transformed by 

98 # block_simplify() in the end to take care of any stray scalars. 

99 while len(ops) > 1: 

100 B = ops.pop() 

101 A = ops.pop() 

102 

103 if isinstance(A, block_mat) and isinstance(B, block_mat): 103 ↛ 111line 103 didn't jump to line 111, because the condition on line 103 was never false

104 m,n = A.blocks.shape 

105 p,q = B.blocks.shape 

106 C = block_mat(m,q) 

107 for row in range(m): 

108 for col in range(q): 

109 for i in range(n): 

110 C[row,col] += block_mul(A[row,i],B[i,col]) 

111 elif isinstance(A, block_mat) and isscalar(B): 

112 m,n = A.blocks.shape 

113 C = block_mat(m,n) 

114 for row in range(m): 

115 for col in range(n): 

116 C[row,col] = block_mul(A[row,col],B) 

117 elif isinstance(B, block_mat) and isscalar(A): 

118 m,n = B.blocks.shape 

119 C = block_mat(m,n) 

120 for row in range(m): 

121 for col in range(n): 

122 C[row,col] = block_mul(A,B[row,col]) 

123 else: 

124 C = block_mul(A,B) 

125 ops.append(C) 

126 return block_simplify(ops[0]) 

127 

128 def block_simplify(self): 

129 """Try to combine scalar terms and remove multiplicative identities, 

130 recursively. A fuller explanation is found in block_transform.block_simplify. 

131 """ 

132 from .block_util import isscalar 

133 from .block_transform import block_simplify 

134 operators = [] 

135 scalar = 1.0 

136 for op in self.chain: 

137 op = block_simplify(op) 

138 if isscalar(op): 

139 scalar *= op 

140 else: 

141 operators.append(op) 

142 if scalar == 0: 

143 return 0 

144 if scalar != 1 or len(operators) == 0: 144 ↛ 146line 144 didn't jump to line 146, because the condition on line 144 was never false

145 operators.insert(0, scalar) 

146 if len(operators) == 1: 

147 return operators[0] 

148 ret = block_mul(None, None) 

149 ret.chain = operators 

150 return ret 

151 

152 def __str__(self): 

153 return '{%s}'%(' * '.join(op.__str__() for op in self.chain)) 

154 

155 def __iter__(self): 

156 return iter(self.chain) 

157 def __len__(self): 

158 return len(self.chain) 

159 def __getitem__(self, i): 

160 return self.chain[i] 

161 

162class block_transpose(block_base): 

163 def __init__(self, A): 

164 self.A = A 

165 def matvec(self, x): 

166 return self.A.transpmult(x) 

167 def transpmult(self, x): 

168 return self.A.__mul__(x) 

169 

170 def create_vec(self, dim=1): 

171 return self.A.create_vec(1-dim) 

172 

173 def block_collapse(self): 

174 """See block_transform.block_collapse.""" 

175 from .block_transform import block_collapse, block_simplify 

176 from .block_mat import block_mat 

177 A = block_collapse(self.A) 

178 if not isinstance(A, block_mat): 178 ↛ 179line 178 didn't jump to line 179, because the condition on line 178 was never true

179 return block_transpose(A) 

180 m,n = A.blocks.shape 

181 ret = block_mat(n,m) 

182 for i in range(m): 

183 for j in range(n): 

184 ret[j,i] = block_transpose(A[i,j]) 

185 return block_simplify(ret) 

186 

187 def block_simplify(self): 

188 """Try to simplify the transpose, recursively. A fuller explanation is 

189 found in block_transform.block_simplify. 

190 """ 

191 from .block_util import isscalar 

192 from .block_transform import block_simplify 

193 A = block_simplify(self.A) 

194 if isscalar(A): 

195 return A 

196 if isinstance(A, block_transpose): 196 ↛ 197line 196 didn't jump to line 197, because the condition on line 196 was never true

197 return A.A 

198 return block_transpose(A) 

199 

200 def __str__(self): 

201 return '<block_transpose of %s>'%str(self.A) 

202 def __iter__(self): 

203 return iter([self.A]) 

204 def __len__(self): 

205 return 1 

206 def __getitem__(self, i): 

207 return [self.A][i] 

208 

209# It's probably best if block_sub and block_add do not allow coercion into 

210# block_mul, since that might mess up the operator precedence. Hence, they 

211# do not inherit from block_base. As it is now, self.A*x and self.B*x must be 

212# reduced to vectors, which means all composed multiplies are finished before 

213# __mul__ does anything. 

214 

215class block_sub(object): 

216 def __init__(self, A, B): 

217 self.A = A 

218 self.B = B 

219 

220 def _combine(self, y, z): 

221 y -= z 

222 

223 def _action(self, x, transposed): 

224 from .block_mat import block_vec 

225 from .block_util import mult 

226 from dolfin import GenericVector 

227 if not isinstance(x, (GenericVector, block_vec)): 227 ↛ 228line 227 didn't jump to line 228, because the condition on line 227 was never true

228 return NotImplemented 

229 y = mult(self.A, x, transposed) 

230 z = mult(self.B, x, transposed) 

231 if len(y) != len(z): 231 ↛ 232line 231 didn't jump to line 232, because the condition on line 231 was never true

232 raise RuntimeError( 

233 'incompatible dimensions in matrix subtraction -- %d != %d'%(len(y),len(z))) 

234 self._combine(y, z) 

235 return y 

236 

237 def __mul__(self, x): 

238 return self._action(x, transposed=False) 

239 

240 def transpmult(self, x): 

241 return self._action(x, transposed=True) 

242 

243 def __neg__(self): 

244 return block_sub(self.B, self.A) 

245 

246 def __add__(self, x): 

247 return block_add(self, x) 

248 

249 def __sub__(self, x): 

250 return block_sub(self, x) 

251 

252 def create_vec(self, dim=1): 

253 """Create a vector of the correct size and layout for b (if dim==0) or 

254 x (if dim==1).""" 

255 try: 

256 return self.A.create_vec(dim) 

257 except AttributeError: 

258 return self.B.create_vec(dim) 

259 

260 def block_simplify(self): 

261 """Try to combine scalar terms and remove additive identities, 

262 recursively. A fuller explanation is found in block_transform.block_simplify. 

263 """ 

264 from .block_util import isscalar 

265 from .block_transform import block_simplify 

266 A = block_simplify(self.A) 

267 B = block_simplify(self.B) 

268 if isscalar(A) and A==0: 

269 return -B 

270 if isscalar(B) and B==0: 

271 return A 

272 return A-B 

273 

274 def block_collapse(self): 

275 """Create a block_mat of block_adds from a block_add of block_mats. See block_transform.block_collapse.""" 

276 from .block_mat import block_mat 

277 from .block_util import isscalar 

278 from .block_transform import block_collapse, block_simplify 

279 

280 A = block_collapse(self.A) 

281 B = block_collapse(self.B) 

282 

283 # The self.__class__(A,B) used below works for both block_sub and 

284 # block_add, and any scalar terms are combined by the final call to 

285 # block_simplify(). 

286 if isinstance(A, block_mat) and isinstance(B, block_mat): 286 ↛ 292line 286 didn't jump to line 292, because the condition on line 286 was never false

287 m,n = A.blocks.shape 

288 C = block_mat(m,n) 

289 for row in range(m): 

290 for col in range(n): 

291 C[row,col] = self.__class__(A[row,col], B[row,col]) 

292 elif isinstance(A, block_mat) and isscalar(B): 

293 m,n = A.blocks.shape 

294 C = block_mat(m,n) 

295 for row in range(m): 

296 for col in range(n): 

297 C[row,col] = self.__class__(A[row,col], B) if row==col else A[row,col] 

298 elif isinstance(B, block_mat) and isscalar(A): 

299 m,n = B.blocks.shape 

300 C = block_mat(m,n) 

301 for row in range(m): 

302 for col in range(n): 

303 C[row,col] = self.__class__(A, B[row,col]) if row==col else B[row,col] 

304 else: 

305 C = self.__class__(A, B) 

306 return block_simplify(C) 

307 

308 def __str__(self): 

309 return '{%s - %s}'%(self.A.__str__(), self.B.__str__()) 

310 

311 def __iter__(self): 

312 return iter([self.A, self.B]) 

313 def __len__(self): 

314 return 2 

315 def __getitem__(self, i): 

316 return [self.A, self.B][i] 

317 

318class block_add(block_sub): 

319 def _combine(self, y, z): 

320 y += z 

321 

322 def block_simplify(self): 

323 """Try to combine scalar terms and remove additive identities, 

324 recursively. A fuller explanation is found in block_transform.block_simplify. 

325 """ 

326 from .block_util import isscalar 

327 from .block_transform import block_simplify 

328 A = block_simplify(self.A) 

329 B = block_simplify(self.B) 

330 if isscalar(A) and A==0: 

331 return B 

332 if isscalar(B) and B==0: 

333 return A 

334 return A+B 

335 

336 def __neg__(self): 

337 return block_mul(-1, self) 

338 

339 def __str__(self): 

340 return '{%s + %s}'%(self.A.__str__(), self.B.__str__())