@@ -42,6 +42,7 @@ def __init__(self, support):
4242 self .A = [] # sparse matrix storage coo format
4343 self .col = []
4444 self .row = [] # sparse matrix storage
45+ self .w = []
4546 self .solver = None
4647 self .eq_const_C = []
4748 self .eq_const_row = []
@@ -135,7 +136,7 @@ def reset(self):
135136 self .B = []
136137 self .n_constraints = 0
137138
138- def add_constraints_to_least_squares (self , A , B , idc , name = 'undefined' ):
139+ def add_constraints_to_least_squares (self , A , B , idc , w = 1. , name = 'undefined' ):
139140 """
140141 Adds constraints to the least squares system. Automatically works
141142 out the row
@@ -166,49 +167,51 @@ def add_constraints_to_least_squares(self, A, B, idc, name='undefined'):
166167
167168 if len (A .shape ) > 2 :
168169 nr = A .shape [0 ] * A .shape [1 ]
170+ w = np .tile (w ,(A .shape [1 ]))
169171 A = A .reshape ((A .shape [0 ]* A .shape [1 ],A .shape [2 ]))
170172 idc = idc .reshape ((idc .shape [0 ]* idc .shape [1 ],idc .shape [2 ]))
173+ B = B .reshape ((A .shape [0 ]))
174+ # w = w.reshape((A.shape[0]))
171175 # normalise by rows of A
172- length = norm (A ,axis = 1 )#.getcol(0).norm()
176+ length = np . linalg . norm (A ,axis = 1 )#.getcol(0).norm()
173177 B [length > 0 ]/= length [length > 0 ]
174- A = normalize (A ,axis = 1 )
175178 # going to assume if any are nan they are all nan
176179 mask = np .any (np .isnan (A ),axis = 1 )
177180 A [mask ,:] = 0
181+ A [length > 0 ,:] /= length [length > 0 ,None ]
182+ if isinstance (w ,float ):
183+ w = np .ones (A .shape [0 ])* w
184+
185+ if isinstance (w ,np .ndarray ) == False :
186+ raise BaseException ('w must be a numpy array' )
187+
188+ if w .shape [0 ] != A .shape [0 ]:
189+ # # make w the same size as A
190+ # w = np.tile(w,(A.shape[1],1)).T
191+ # else:
192+ raise BaseException ('Weight array does not match number of constraints' )
178193 if np .any (np .isnan (idc )) or np .any (np .isnan (A )) or np .any (np .isnan (B )):
179194 logger .warning ("Constraints contain nan not adding constraints: {}" .format (name ))
180195 # return
181-
182196 rows = np .arange (0 , nr ).astype (int )
183197 rows += self .c_
184198 constraint_ids = rows .copy ()
185-
186- if name in self .constraints :
199+ base_name = name
200+ while name in self .constraints :
187201 count = 0
188202 if '_' in name :
189- count = int (name .split ('_' )[0 ])+ 1
190- name = name + '_{}' .format (count )
203+ count = int (name .split ('_' )[1 ])+ 1
204+ name = base_name + '_{}' .format (count )
191205
192206 # self.constraints[name]['A'] = A#np.vstack([self.constraints[name]['A'],A])
193207 # self.constraints[name]['B'] = B#np.hstack([self.constraints[name]['B'], B])
194208 # self.constraints[name]['idc'] = idc#np.vstack([self.constraints[name]['idc'],
195209 # idc])
196-
197- self .constraints [name ] = {'node_indexes' :constraint_ids ,'A' :A ,'B' :B .flatten (),'idc' :idc }
198210 rows = np .tile (rows , (A .shape [- 1 ], 1 )).T
211+ self .constraints [name ] = {'node_indexes' :constraint_ids ,'A' :A ,'B' :B .flatten (),'col' :idc ,'w' :w ,'row' :rows }
199212
200213 self .c_ += nr
201- if self .shape == 'rectangular' :
202- # don't add operator where it is = 0 to the sparse matrix!
203- A = A .flatten ()
204- rows = rows .flatten ()
205- idc = idc .flatten ()
206- B = B .flatten ()
207- mask = A == 0
208- self .A .extend (A [~ mask ].tolist ())
209- self .row .extend (rows [~ mask ].tolist ())
210- self .col .extend (idc [~ mask ].tolist ())
211- self .B .extend (B .tolist ())
214+
212215
213216 def calculate_residual_for_constraints (self ):
214217 residuals = {}
@@ -318,11 +321,41 @@ def build_matrix(self, square=True, damp=True):
318321
319322 logger .info ("Interpolation matrix is %i x %i" % (self .c_ ,self .nx ))
320323 cols = np .array (self .col )
321- A = coo_matrix ((np .array (self .A ), (np .array (self .row ), \
324+ # To keep the solvers consistent for different model scales the range of the constraints should be similar.
325+ # We normalise the row vectors for the interpolation matrix
326+ # Each constraint can then be weighted separately for the least squares problem
327+ # The weights are normalised so that the max weight is 1.0
328+ # This means that the tolerance and other parameters for the solver
329+ # are kept the same between iterations.
330+ # #TODO currently the element size is not incorporated into the weighting.
331+ # For cartesian grids this is probably ok but for tetrahedron could be more problematic if
332+ # the tetras have different volumes. Would expect for the size of the element to influence
333+ # how much it contributes to the system.
334+ # It could be implemented by multiplying the weight array by the element size.
335+ # I am not sure how to integrate regularisation into this framework as my gut feeling is the regularisation
336+ # should be weighted by the area of the element face and not element volume, but this means the weight decreases with model scale
337+ # which is not ideal.
338+ max_weight = 0
339+ for c in self .constraints .values ():
340+ if c ['w' ].max () > max_weight :
341+ max_weight = c ['w' ].max ()
342+ a = []
343+ b = []
344+ rows = []
345+ cols = []
346+ for c in self .constraints .values ():
347+ aa = (c ['A' ]* c ['w' ][:,None ]/ max_weight ).flatten ()
348+ b .extend ((c ['B' ]* c ['w' ]/ max_weight ).tolist ())
349+ mask = aa == 0
350+ a .extend (aa [~ mask ].tolist ())
351+ rows .extend (c ['row' ].flatten ()[~ mask ].tolist ())
352+ cols .extend (c ['col' ].flatten ()[~ mask ].tolist ())
353+
354+ A = coo_matrix ((np .array (a ), (np .array (rows ), \
322355 cols )), shape = (self .c_ , self .nx ),
323356 dtype = float ) # .tocsr()
324- B = np . array ( self . B )
325-
357+
358+ B = np . array ( b )
326359 if not square :
327360 logger .info ("Using rectangular matrix, equality constraints are not used" )
328361 return A , B
0 commit comments