/[escript]/trunk/downunder/py_src/splitinversioncostfunctions.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/splitinversioncostfunctions.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5654 - (hide annotations)
Mon Jun 15 05:21:00 2015 UTC (4 years, 1 month ago) by jfenwick
File MIME type: text/x-python
File size: 42255 byte(s)
further
1 jfenwick 5622 from __future__ import division, print_function
2 jfenwick 5562 ##############################################################################
3     #
4 jfenwick 5593 # Copyright (c) 2003-2015 by The University of Queensland
5 jfenwick 5562 # http://www.uq.edu.au
6     #
7     # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10     #
11     # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development 2012-2013 by School of Earth Sciences
13     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14     #
15     ##############################################################################
16    
17     """Cost functions for inversions with one or more forward models"""
18    
19 jfenwick 5593 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland
20 jfenwick 5562 http://www.uq.edu.au
21     Primary Business: Queensland, Australia"""
22     __license__="""Licensed under the Open Software License version 3.0
23     http://www.opensource.org/licenses/osl-3.0.php"""
24     __url__="https://launchpad.net/escript-finley"
25    
26     __all__ = [ 'SplitInversionCostFunction']
27    
28     from .costfunctions import MeteredCostFunction
29     from .mappings import Mapping
30     from .forwardmodels import ForwardModel
31     from esys.escript.pdetools import ArithmeticTuple
32 jfenwick 5637 from esys.escript import Data, inner, addJobPerWorld, addVariable, makeLocalOnly, makeScalarReducer, makeDataReducer, FunctionJob, Job
33 jfenwick 5562 import numpy as np
34    
35    
36     class SplitInversionCostFunction(MeteredCostFunction):
37     """
38     Class to define cost function *J(m)* for inversion with one or more
39     forward models based on a multi-valued level set function *m*:
40    
41     *J(m) = J_reg(m) + sum_f mu_f * J_f(p)*
42    
43     where *J_reg(m)* is the regularization and cross gradient component of the
44     cost function applied to a level set function *m*, *J_f(p)* are the data
45     defect cost functions involving a physical forward model using the
46     physical parameter(s) *p* and *mu_f* is the trade-off factor for model f.
47    
48     A forward model depends on a set of physical parameters *p* which are
49     constructed from components of the level set function *m* via mappings.
50    
51     Example 1 (single forward model):
52     m=Mapping()
53     f=ForwardModel()
54     J=InversionCostFunction(Regularization(), m, f)
55    
56     Example 2 (two forward models on a single valued level set)
57     m0=Mapping()
58     m1=Mapping()
59     f0=ForwardModel()
60     f1=ForwardModel()
61    
62     J=InversionCostFunction(Regularization(), mappings=[m0, m1], forward_models=[(f0, 0), (f1,1)])
63    
64     Example 3 (two forward models on 2-valued level set)
65     m0=Mapping()
66     m1=Mapping()
67     f0=ForwardModel()
68     f1=ForwardModel()
69    
70     J=InversionCostFunction(Regularization(self.numLevelSets=2), mappings=[(m0,0), (m1,0)], forward_models=[(f0, 0), (f1,1)])
71    
72     :cvar provides_inverse_Hessian_approximation: if true the class provides an
73     approximative inverse of the Hessian operator.
74     """
75     provides_inverse_Hessian_approximation=True
76    
77     # Original params for InversionCostFunction
78     # regularization --- need to feed in later
79     # mappings --- could take a domain so need to be created later
80     # forward_models --- need to be split (do we have a function which takes a parameter to indicate which model to create?)
81     #
82     # New constructor
83     # num args, who many of each type
84     # splitw is the splitworld jobs are running on
85     # worldsinit_fn is run on each world at startup
86 jfenwick 5622 def __init__(self, numLevelSets=None, numModels=None, numMappings=None, splitworld=None, worldsinit_fn=None):
87 jfenwick 5562 """
88     fill this in.
89     """
90 jfenwick 5622 import math
91     if numLevelSets==None or numModels==None or numMappings==None or splitworld==None or worldsinit_fn==None:
92     raise ValueError("Please supply all required parameters")
93 jfenwick 5562 super(SplitInversionCostFunction, self).__init__()
94     if numModels<1 or numModels<1 or numMappings<1:
95 jfenwick 5563 raise ValueError("The inversion function requires at least one LevelSet, Mapping and Models.")
96 jfenwick 5562 self.numModels=numModels
97     self.numMappings=numMappings
98     self.numLevelSets=numLevelSets
99 jfenwick 5630 self.splitworld=splitworld
100 jfenwick 5622
101 jfenwick 5633 addVariable(splitworld,"regularization", makeLocalOnly)
102     addVariable(splitworld,"mappings", makeLocalOnly)
103     addVariable(splitworld,"fwdmodels", makeLocalOnly)
104     addVariable(splitworld,"initial_guess", makeLocalOnly) # Used to load the initial guess
105     addVariable(splitworld,"model_args", makeLocalOnly) # arguments for models stored on that world
106     addVariable(splitworld,"props", makeLocalOnly) # Properties for the current guess
107     addVariable(splitworld,"current_point", makeLocalOnly) # Current approximate solution. Starts out as initial_guess
108     addVariable(splitworld,"mu_model", makeLocalOnly)
109 jfenwick 5635
110 jfenwick 5648 addVariable(splitworld, "phi_a", makeScalarReducer, "SUM")
111 jfenwick 5654 addVariable(splitworld,"Jx_original", makeScalarReducer,"SET")
112 jfenwick 5635 addVariable(splitworld, "Jx", makeScalarReducer, "SUM")
113 jfenwick 5654 addVariable(splitworld, "Jx_old", makeScalarReducer,"SET")
114 jfenwick 5637 addVariable(splitworld, "g_Jx_0", makeDataReducer, "SUM")
115 jfenwick 5638 addVariable(splitworld, "g_Jx_1", makeLocalOnly) # This component is not merged with values from other worlds
116 jfenwick 5654
117     addVariable(splitworld, "old_g_Jx_0", makeDataReducer, "SUM")
118     addVariable(splitworld, "old_g_Jx_1", makeLocalOnly) # This component is not merged with values from other worlds
119    
120    
121 jfenwick 5643 addVariable(splitworld, "search_direction", makeDataReducer, "SET")
122 jfenwick 5645
123     addVariable(splitworld, "s_and_y", makeLocalOnly)
124 jfenwick 5648 addVariable(splitworld, "gphi0", makeLocalOnly)
125     addVariable(splitworld, "old_phi_a", makeLocalOnly)
126     addVariable(splitworld, "phi0", makeLocalOnly)
127     addVariable(splitworld, "base_point", makeLocalOnly)
128 jfenwick 5633
129 jfenwick 5654 addVariable(splitworld, "conv_flag", makeLocalOnly)
130    
131 jfenwick 5622 howmany=splitworld.getSubWorldCount()
132     rlen=int(math.ceil(numModels/howmany))
133     rstart=rlen*splitworld.getSubWorldID()
134     extraparams={'rangelen':rlen, 'rangestart':rstart, 'numLevelSets':numLevelSets}
135 jfenwick 5562 # sanity check
136 jfenwick 5622 addJobPerWorld(splitworld, FunctionJob, worldsinit_fn, **extraparams)
137     splitworld.runJobs()
138 jfenwick 5630 #reqd=["fwdmodels", "regularization", "mappings","mu_model"]
139     reqd=["fwdmodels", "regularization", "mappings", "initial_guess"] #For our script, mu_model appears not to be used
140 jfenwick 5622 knownvars=splitworld.getVarList()
141 jfenwick 5563 print(knownvars)
142 jfenwick 5562 for n in reqd:
143 jfenwick 5563 if [n,True] not in knownvars:
144     raise RuntimeError("Required variable "+n+" was not created by the world init function")
145 jfenwick 5633 if ['mu_model',True] not in knownvars:
146     self.setTradeOffFactorsModels()
147 jfenwick 5630 self.configured=True
148 jfenwick 5562
149     # Function to put the (possible list of) forward model(s) into the form expected by the rest of the system
150     @staticmethod
151     def formatModels(forward_models, numMappings):
152     if isinstance(forward_models, ForwardModel):
153     forward_models = [ forward_models ]
154     result=[]
155     for i in range(len(forward_models)):
156 jfenwick 5563 print("Doing iteration "+str(i))
157 jfenwick 5562 f=forward_models[i]
158     if isinstance(f, ForwardModel):
159     idx=[0]
160     fm=f
161     elif len(f) == 1:
162     idx=[0]
163     fm=f[0]
164     else:
165     if isinstance(f[1],int):
166     idx=[f[1]]
167     else:
168     idx=list(f[1])
169     for k in idx:
170     if k<0 or k> numMappings:
171     raise ValueError("mapping index %s in model %s is out of range."%(k,i))
172     fm=f[0]
173     result.append((fm,idx))
174     return result
175 jfenwick 5630
176     # Function to put the (possible list of) forward model(s) into a form expected by the rest of the system
177     @staticmethod
178     def formatMappings(mappings, numLevelSets):
179 jfenwick 5634 if isinstance(mappings, Mapping):
180     mappings = [ mappings ]
181 jfenwick 5630 newmappings = []
182     for i in range(len(mappings)):
183     mm=mappings[i]
184     if isinstance(mm, Mapping):
185     m=mm
186     if numLevelSets>1:
187     idx=[ p for p in range(numLevelSets)]
188     else:
189     idx=None
190     elif len(mm) == 1:
191     m=mm[0]
192     if numLevelSets>1:
193     idx=[ p for p in range(numLevelSets)]
194     else:
195     idx=None
196     else:
197     m=mm[0]
198     if isinstance(mm[1], int):
199     idx=[mm[1]]
200     else:
201     idx=list(mm[1])
202     if numLevelSets>1:
203     for k in idx:
204     if k < 0 or k > numLevelSets-1:
205     raise ValueError("level set index %s is out of range."%(k,))
206    
207     else:
208     if idx[0] != 0:
209     raise ValueError("Level set index %s is out of range."%(idx[0],))
210     else:
211     idx=None
212     newmappings.append((m,idx))
213     return newmappings
214    
215 jfenwick 5562
216     def getDomain(self):
217     """
218     returns the domain of the cost function
219    
220     :rtype: `Domain`
221     """
222     raise RuntimeError("Can't extract domains for SplitInversionCostFunctions")
223    
224     def getNumTradeOffFactors(self):
225     """
226     returns the number of trade-off factors being used including the
227     trade-off factors used in the regularization component.
228    
229     :rtype: ``int``
230     """
231     if self.configured:
232     return self.__num_tradeoff_factors
233     else:
234 jfenwick 5563 raise RuntimeError("This inversion function has not been configured yet")
235 jfenwick 5562
236     def getForwardModel(self, idx=None):
237     """
238     returns the *idx*-th forward model.
239    
240     :param idx: model index. If cost function contains one model only `idx`
241     can be omitted.
242     :type idx: ``int``
243     """
244     raise RuntimeError("Can't extract forward models for SplitInversionCostFunctions")
245    
246     def getRegularization(self):
247     """
248     returns the regularization
249    
250     :rtype: `Regularization`
251     """
252     if self.configured:
253     return self.regularization
254     else:
255 jfenwick 5563 raise RuntimeError("This inversion function has not been configured yet")
256 jfenwick 5562
257 jfenwick 5633 #Written to be executed inside a FunctionJob
258     @staticmethod
259     def subworld_setMu_model(self, **args):
260     if not isinstance(self, Job):
261     raise RuntimeError("This command should be run inside a Job")
262     extmu=args['mu']
263     chunksize=max(len(extmu)//self.swcount,1) #In case we have more worlds than models
264     minindex=self.swid*chunksize
265     maxindex=(self.swid+1)*chunksize # yes this could go off the end but I will slice
266     mymu=extmu[minindex:maxindex]
267     self.exportValue("mu_model", mymu)
268    
269 jfenwick 5562 def setTradeOffFactorsModels(self, mu=None):
270     """
271     sets the trade-off factors for the forward model components.
272    
273     :param mu: list of the trade-off factors. If not present ones are used.
274     :type mu: ``float`` in case of a single model or a ``list`` of
275     ``float`` with the length of the number of models.
276     """
277     if mu==None:
278     self.mu_model=np.ones((self.numModels, ))
279     else:
280     if self.numModels > 1:
281     mu=np.asarray(mu, dtype=float)
282     if min(mu) > 0:
283     self.mu_model= mu
284     else:
285     raise ValueError("All values for trade-off factor mu must be positive.")
286     else:
287     mu=float(mu)
288     if mu > 0:
289     self.mu_model= [mu, ]
290     else:
291     raise ValueError("Trade-off factor must be positive.")
292 jfenwick 5633 addJobPerWorld(self.splitworld, FunctionJob, self.subworld_setMu_model, mu=self.mu_model)
293     self.splitworld.runJobs()
294 jfenwick 5562
295 jfenwick 5633
296    
297 jfenwick 5562 def getTradeOffFactorsModels(self):
298     """
299     returns the trade-off factors for the forward models
300    
301     :rtype: ``float`` or ``list`` of ``float``
302     """
303     if not self.configured:
304 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
305 jfenwick 5562 if self.numModels>1:
306     return self.mu_model
307     else:
308     return self.mu_model[0]
309    
310     def setTradeOffFactorsRegularization(self, mu=None, mu_c=None):
311     """
312     sets the trade-off factors for the regularization component of the
313     cost function, see `Regularization` for details.
314    
315     :param mu: trade-off factors for the level-set variation part
316     :param mu_c: trade-off factors for the cross gradient variation part
317     """
318     self.regularization.setTradeOffFactorsForVariation(mu)
319     self.regularization.setTradeOffFactorsForCrossGradient(mu_c)
320    
321     def setTradeOffFactors(self, mu=None):
322     """
323     sets the trade-off factors for the forward model and regularization
324     terms.
325    
326     :param mu: list of trade-off factors.
327     :type mu: ``list`` of ``float``
328     """
329 jfenwick 5630 if not self.configured:
330     raise ValueError("This inversion function has not been configured yet")
331     raise ValueError("setTradeOffFactors not supported yet.")
332 jfenwick 5562 if mu is None:
333     mu=np.ones((self.__num_tradeoff_factors,))
334     self.setTradeOffFactorsModels(mu[:self.numModels])
335     self.regularization.setTradeOffFactors(mu[self.numModels:])
336    
337     def getTradeOffFactors(self, mu=None):
338     """
339     returns a list of the trade-off factors.
340    
341     :rtype: ``list`` of ``float``
342     """
343     if not self.configured:
344 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
345 jfenwick 5562 mu1=self.getTradeOffFactorsModels(mu[:self.numModels])
346     mu2=self.regularization.getTradeOffFactors()
347     return [ m for m in mu1] + [ m for m in mu2]
348    
349 jfenwick 5630 @staticmethod
350     def createLevelSetFunctionHelper(self, regularization, mappings, *props):
351     """
352     Returns an object (init-ed) with 0s.
353     Components can be overwritten by physical
354     properties `props`. If present entries must correspond to the
355     `mappings` arguments in the constructor. Use ``None`` for properties
356     for which no value is given.
357     """
358     if not isinstance(self, Job):
359     raise RuntimeError("This function is designed to be run inside a Job.")
360     m=regularization.getPDE().createSolution()
361     if len(props) > 0:
362     numMappings=len(mappings)
363     for i in range(numMappings):
364     if props[i]:
365     mp, idx=self.mappings[i]
366     m2=mp.getInverse(props[i])
367     if idx:
368     if len(idx) == 1:
369     m[idx[0]]=m2
370     else:
371     for k in range(idx): m[idx[k]]=m2[k]
372     else:
373     m=m2
374     return m
375    
376     @staticmethod
377     def calculatePropertiesHelper(self, m, mappings):
378     """
379     returns a list of the physical properties from a given level set
380     function *m* using the mappings of the cost function.
381    
382     :param m: level set function
383     :type m: `Data`
384     :rtype: ``list`` of `Data`
385     """
386     if not isinstance(self, Job):
387     raise RuntimeError("This function is designed to be run inside a Job.")
388     props=[]
389     for i in range(len(mappings)):
390     mp, idx=mappings[i]
391     if idx:
392     if len(idx)==1:
393     p=mp.getValue(m[idx[0]])
394     else:
395     m2=Data(0.,(len(idx),),m.getFunctionSpace())
396     for k in range(len(idx)): m2[k]=m[idx[k]]
397     p=mp.getValue(m2)
398     else:
399     p=mp.getValue(m)
400     props.append(p)
401     return props
402    
403 jfenwick 5562 def createLevelSetFunction(self, *props):
404     """
405     returns an instance of an object used to represent a level set function
406     initialized with zeros. Components can be overwritten by physical
407     properties `props`. If present entries must correspond to the
408     `mappings` arguments in the constructor. Use ``None`` for properties
409     for which no value is given.
410     """
411     if not self.configured:
412 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
413 jfenwick 5562 #Since this involves solving a PDE (and therefore a domain, it must be kept local
414     #to each subworld
415     raise RuntimeError("This needs to run inside the subworld --- create a function for it")
416     m=self.regularization.getPDE().createSolution()
417     if len(props) > 0:
418     for i in range(self.numMappings):
419     if props[i]:
420     mp, idx=self.mappings[i]
421     m2=mp.getInverse(props[i])
422     if idx:
423     if len(idx) == 1:
424     m[idx[0]]=m2
425     else:
426     for k in range(idx): m[idx[k]]=m2[k]
427     else:
428     m=m2
429     #return m
430    
431     def getProperties(self, m, return_list=False):
432     """
433     returns a list of the physical properties from a given level set
434     function *m* using the mappings of the cost function.
435    
436     :param m: level set function
437     :type m: `Data`
438     :param return_list: if ``True`` a list is returned.
439     :type return_list: ``bool``
440     :rtype: ``list`` of `Data`
441     """
442     if not self.configured:
443 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
444 jfenwick 5562 #Since this involves solving a PDE (and therefore a domain, it must be kept local
445     #to each subworld
446     raise RuntimeError("This needs to run inside the subworld --- create a function for it")
447     props=[]
448     for i in range(self.numMappings):
449     mp, idx=self.mappings[i]
450     if idx:
451     if len(idx)==1:
452     p=mp.getValue(m[idx[0]])
453     else:
454     m2=Data(0.,(len(idx),),m.getFunctionSpace())
455     for k in range(len(idx)): m2[k]=m[idx[k]]
456     p=mp.getValue(m2)
457     else:
458     p=mp.getValue(m)
459     props.append(p)
460     #if self.numMappings > 1 or return_list:
461     #return props
462     #else:
463     #return props[0]
464    
465     def _getDualProduct(self, x, r):
466     """
467     Returns the dual product, see `Regularization.getDualProduct`
468    
469     :type x: `Data`
470     :type r: `ArithmeticTuple`
471     :rtype: ``float``
472     """
473     if not self.configured:
474 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
475 jfenwick 5562 #This involves an 'x' which is a Data object.
476     #Does this only need to run on one subworld?
477     #Then shipped around using getDoubleValues?
478     raise RuntimeError("Still need to work this one out")
479     return self.regularization.getDualProduct(x, r)
480    
481 jfenwick 5642
482     @staticmethod
483     def update_point_helper(self, newpoint):
484     """
485     Call within a subworld to set 'current_point' to newpoint
486     and update all the cached args info
487     """
488     if not isinstance(self, Job):
489     raise RuntimeError("This function should only be called from within a Job")
490     mods=self.importValue("fwdmodels")
491     reg=self.importValue("regularization")
492     mappings=self.importValue("mappings")
493     props=[]
494     props=SplitInversionCostFunction.calculatePropertiesHelper(self, newpoint, mappings)
495     self.exportValue("props", props)
496     reg.setPoint(newpoint)
497     #Going to try this - each world stores the args for its
498     #models rather than going the setPoint route.
499     local_args=[]
500     for m,idx in mods:
501     pp=tuple( [props[k] for k in idx] ) # build up collection of properties used by this model
502     local_args.append(m.getArguments(*pp))
503     self.exportValue("current_point", newpoint)
504     self.exportValue("model_args", local_args)
505    
506 jfenwick 5562 def setPoint(self):
507 jfenwick 5630 self._setPoint()
508 jfenwick 5562
509     def _setPoint(self):
510     """
511     This should take in a value to set the point to, but that can wait
512 jfenwick 5637 ... It probably the shouldn't actually. We want all values to
513     be constructed inside the subworlds, so being able to (easily - we
514     can't stop closures) pass them
515     in from outside would defeat the purpose.
516 jfenwick 5630
517 jfenwick 5642 To modify the point, we probably want a separate move_point()
518 jfenwick 5637 function.
519    
520 jfenwick 5630 There is also the question of how this is expected to get its info.
521     Should it be passed in as a parameter or should it be read from
522     the environment?
523     We can expect the actuall initial guess to come from the world init
524     function, but what about later calls? (or are we hoping they won't
525     actually happen that often and that relative changes will be done instead?)
526    
527 jfenwick 5562 """
528     if not self.configured:
529 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
530 jfenwick 5562
531 jfenwick 5637 def load_guess_to_subworlds(self, **args):
532 jfenwick 5563 reg=self.importValue("regularization")
533 jfenwick 5630 mappings=self.importValue("mappings")
534 jfenwick 5637 # we are not passing in property values here because we don't have any yet
535     initguess=SplitInversionCostFunction.createLevelSetFunctionHelper(self, reg, mappings)
536 jfenwick 5642 SplitInversionCostFunction.update_point_helper(self, initguess)
537 jfenwick 5563
538 jfenwick 5637 addJobPerWorld(self.splitworld, FunctionJob, load_guess_to_subworlds, imports=["fwdmodels", "regularization", "mappings"])
539 jfenwick 5630 self.splitworld.runJobs()
540 jfenwick 5562
541     def _getArguments(self, m):
542     """
543     returns pre-computed values that are shared in the calculation of
544     *J(m)* and *grad J(m)*. In this implementation returns a tuple with the
545     mapped value of ``m``, the arguments from the forward model and the
546     arguments from the regularization.
547    
548     :param m: current approximation of the level set function
549     :type m: `Data`
550     :return: tuple of of values of the parameters, pre-computed values
551     for the forward model and pre-computed values for the
552     regularization
553     :rtype: ``tuple``
554     """
555     if not self.configured:
556 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
557 jfenwick 5562 raise RuntimeError("Call to getArguments -- temporary block to see where this is used")
558     args_reg=self.regularization.getArguments(m)
559     # cache for physical parameters:
560     props=self.getProperties(m, return_list=True)
561     args_f=[]
562     for i in range(self.numModels):
563     f, idx=self.forward_models[i]
564     pp=tuple( [ props[k] for k in idx] )
565     aa=f.getArguments(*pp)
566     args_f.append(aa)
567    
568     return props, args_f, args_reg
569    
570     def calculateValue(self, vnames):
571 jfenwick 5631 self._calculateValue(vnames)
572 jfenwick 5562
573     def _calculateValue(self, vnames):
574    
575     if not self.configured:
576 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
577 jfenwick 5562 #The props is already in each world as a variable
578 jfenwick 5631 #Each world has the arguments for the point for all of its models
579     # as a variable.
580     #Regularization already has its point set
581 jfenwick 5562
582 jfenwick 5631 def calculateValueWorker(self, **args):
583 jfenwick 5563 props=self.importValue("props")
584 jfenwick 5631 mods=self.importValue("fwdmodels")
585 jfenwick 5563 reg=self.importValue("regularization")
586     mu_model=self.importValue("mu_model")
587 jfenwick 5634 local_args=self.importValue("model_args")
588 jfenwick 5631 current_point=self.importValue("current_point")
589     try:
590 jfenwick 5634 vnames=args['vnames']
591     except KeyError as e:
592     raise RuntimeError("Function requires vnames as kwarg")
593 jfenwick 5563 J=None
594     for i in range(len(mods)): # note: iterating over local models not ones on other worlds
595     m,idx=mods[i]
596 jfenwick 5634 args=local_args[i]
597 jfenwick 5631 z=m.getDefect(current_point, *args)
598 jfenwick 5633 z*=mu_model[i];
599 jfenwick 5563 if J is None:
600     J=z
601     else:
602 jfenwick 5633 J+=z
603     if self.swid==0: # we only want to add the regularization term once
604 jfenwick 5563 J+=reg.getValueAtPoint() # We actually want to get a value here but
605     # I want to distiguish it from the other getValue call
606     if isinstance(vnames, str):
607 jfenwick 5648 self.exportValue(vnames, J)
608 jfenwick 5563 else:
609     for n in vnames:
610 jfenwick 5635 self.exportValue(n,J)
611 jfenwick 5648 # End calculateValueWorker
612    
613 jfenwick 5633 addJobPerWorld(self.splitworld,FunctionJob, calculateValueWorker, imports=["fwdmodels", "regularization", "props",
614     "model_args", "mu_model"], vnames=vnames)
615     self.splitworld.runJobs()
616     # The result will now be stored in the named variables
617     # The caller will need to execute splitworld.getDoubleVariable to extract them
618    
619 jfenwick 5562
620     def _getValue(self, m, *args):
621     """
622     Returns the value *J(m)* of the cost function at *m*.
623     If the pre-computed values are not supplied `getArguments()` is called.
624    
625     :param m: current approximation of the level set function
626     :type m: `Data`
627     :param args: tuple of values of the parameters, pre-computed values
628     for the forward model and pre-computed values for the
629     regularization
630     :rtype: ``float``
631     """
632     if not self.configured:
633 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
634 jfenwick 5562 raise RuntimeError("Call to getArguments -- temporary block to see where this is used")
635     if len(args)==0:
636     args=self.getArguments(m)
637    
638     props=args[0]
639     args_f=args[1]
640     args_reg=args[2]
641    
642     J = self.regularization.getValue(m, *args_reg)
643     self.logger.debug("J_R (incl. trade-offs) = %e"%J)
644    
645     for i in range(self.numModels):
646     f, idx=self.forward_models[i]
647     args=tuple( [ props[k] for k in idx] + list( args_f[i] ) )
648     J_f = f.getDefect(*args)
649     self.logger.debug("J_f[%d] = %e, mu_model[%d] = %e"%(i, J_f, i, self.mu_model[i]))
650     J += self.mu_model[i] * J_f
651    
652     return J
653    
654     def getComponentValues(self, m, *args):
655     if not self.configured:
656 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
657 jfenwick 5562 raise RuntimeError("Call to getComponentValues -- temporary block to see where this is used")
658     return self._getComponentValues(m, *args)
659    
660     def _getComponentValues(self, m, *args):
661     """
662     returns the values of the individual cost functions that make up *f(x)*
663     using the precalculated values for *x*.
664    
665     :param x: a solution approximation
666     :type x: x-type
667     :rtype: ``list<<float>>``
668     """
669     if not self.configured:
670 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
671 jfenwick 5562 raise RuntimeError("Call to getArguments -- temporary block to see where this is used")
672     if len(args)==0:
673     args=self.getArguments(m)
674    
675     props=args[0]
676     args_f=args[1]
677     args_reg=args[2]
678    
679     J_reg = self.regularization.getValue(m, *args_reg)
680     result = [J_reg]
681    
682     for i in range(self.numModels):
683     f, idx=self.forward_models[i]
684     args=tuple( [ props[k] for k in idx] + list( args_f[i] ) )
685     J_f = f.getValue(*args)
686     self.logger.debug("J_f[%d] = %e, mu_model[%d] = %e"%(i, J_f, i, self.mu_model[i]))
687    
688     result += [J_f] # self.mu_model[i] * ??
689    
690     return result
691    
692 jfenwick 5637 @staticmethod
693     def getModelArgs(self, fwdmodels):
694 jfenwick 5638 """
695     Attempts to import the arguments for forward models, if they are not available,
696     Computes and exports them
697     """
698     if not isinstance(self, Job):
699     raise RuntimeError("This function should only be called inside a Job")
700     args=self.importValue("model_args")
701     p=self.importValue("current_point")
702     if args is not None:
703     return args
704     args=[]
705     for mod in fwdmodels:
706     args.append(mod.getArguments(p))
707     self.exportValue("model_arguments",args)
708     return args
709    
710 jfenwick 5637 def calculateGradient(self, vnames1, vnames2):
711 jfenwick 5638 """
712     The gradient operation produces two components (designated (Y^,X) in the non-split version).
713     vnames1 gives the variable name(s) where the first component should be stored.
714     vnames2 gives the variable name(s) where the second component should be stored.
715     """
716 jfenwick 5637 return self._calculateGradient(vnames1, vnames2)
717 jfenwick 5635
718 jfenwick 5637 def _calculateGradient(self, vnames1, vnames2):
719 jfenwick 5562 if not self.configured:
720 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
721 jfenwick 5562
722 jfenwick 5563 numLevelSets=self.numLevelSets # pass in via closure
723 jfenwick 5635 def calculateGradientWorker(self, **args):
724 jfenwick 5563 """
725     vnames1 gives the names to store the first component of the gradient in
726     vnames2 gives the names to store the second component of the gradient in
727     """
728 jfenwick 5635 vnames1=args['vnames1']
729     vnames2=args['vnames2']
730 jfenwick 5563 props=self.importValue("props")
731 jfenwick 5637 mods=self.importValue("fwdmodels")
732 jfenwick 5563 reg=self.importValue("regularization")
733     mu_model=self.importValue("mu_model")
734     mappings=self.importValue("mappings")
735 jfenwick 5637 m=self.importValue("current_point")
736 jfenwick 5563
737 jfenwick 5637 model_args=SplitInversionCostFunction.getModelArgs(self, mods)
738    
739 jfenwick 5563 g_J = reg.getGradientAtPoint()
740     p_diffs=[]
741     # Find the derivative for each mapping
742     # If a mapping has a list of components (idx), then make a new Data object with only those
743     # components, pass it to the mapping and get the derivative.
744 jfenwick 5637 for i in range(len(mappings)):
745 jfenwick 5563 mm, idx=mappings[i]
746     if idx and numLevelSets > 1:
747     if len(idx)>1:
748     m2=Data(0,(len(idx),),m.getFunctionSpace())
749     for k in range(len(idx)): m2[k]=m[idx[k]]
750     dpdm = mm.getDerivative(m2)
751     else:
752     dpdm = mm.getDerivative(m[idx[0]])
753     else:
754     dpdm = mm.getDerivative(m)
755     p_diffs.append(dpdm)
756     #Since we are going to be merging Y with other worlds, we need to make sure the the regularization
757     #component is only added once. However most of the ops below are in terms of += so we need to
758     #create a zero object to use as a starting point
759 jfenwick 5637 if self.swid==0:
760 jfenwick 5563 Y=g_J[0] # Because g_J==(Y,X) Y_k=dKer/dm_k
761     else:
762     Y=Data(0, g_J[0].getShape(), g_J[0].getForwardModel())
763 jfenwick 5637 for i in range(len(mods)):
764     mu=mu_model[i]
765 jfenwick 5563 f, idx_f=mods[i]
766 jfenwick 5637 args=tuple( [ props[k] for k in idx_f] + list( model_args[i] ) )
767     Ys = f.getGradient(*args) # this d Jf/d props
768 jfenwick 5563 # in this case f depends on one parameter props only but this can
769     # still depend on several level set components
770     if Ys.getRank() == 0:
771     # run through all level sets k prop j is depending on:
772 jfenwick 5637 idx_m=mappings[idx_f[0]][1]
773 jfenwick 5563 # tmp[k] = dJ_f/d_prop * d prop/d m[idx_m[k]]
774     tmp=Ys * p_diffs[idx_f[0]] * mu
775     if idx_m:
776     if tmp.getRank()== 0:
777     for k in range(len(idx_m)):
778     Y[idx_m[k]]+=tmp # dJ_f /d m[idx_m[k]] = tmp
779     else:
780     for k in range(len(idx_m)):
781     Y[idx_m[k]]+=tmp[k] # dJ_f /d m[idx_m[k]] = tmp[k]
782     else:
783     Y+=tmp # dJ_f /d m[idx_m[k]] = tmp
784     else:
785     s=0
786     # run through all props j forward model f is depending on:
787     for j in range(len(idx_f)):
788     # run through all level sets k prop j is depending on:
789 jfenwick 5637 idx_m=mappings[j][1]
790 jfenwick 5563 if p_diffs[idx_f[j]].getRank() == 0 :
791     if idx_m: # this case is not needed (really?)
792     raise RuntimeError("something wrong A")
793     # tmp[k] = dJ_f/d_prop[j] * d prop[j]/d m[idx_m[k]]
794     tmp=Ys[s]*p_diffs[idx_f[j]] * mu
795     for k in range(len(idx_m)):
796     Y[idx_m[k]]+=tmp[k] # dJ_f /d m[idx_m[k]] = tmp[k]
797     else:
798     Y+=Ys[s]*p_diffs[idx_f[j]] * mu
799     s+=1
800     elif p_diffs[idx_f[j]].getRank() == 1 :
801     l=p_diffs[idx_f[j]].getShape()[0]
802     # tmp[k]=sum_j dJ_f/d_prop[j] * d prop[j]/d m[idx_m[k]]
803     tmp=inner(Ys[s:s+l], p_diffs[idx_f[j]]) * mu
804     if idx_m:
805     for k in range(len(idx_m)):
806     Y[idx_m[k]]+=tmp # dJ_f /d m[idx_m[k]] = tmp[k]
807     else:
808     Y+=tmp
809     s+=l
810     else: # rank 2 case
811     l=p_diffs[idx_f[j]].getShape()[0]
812     Yss=Ys[s:s+l]
813     if idx_m:
814     for k in range(len(idx_m)):
815     # dJ_f /d m[idx_m[k]] = tmp[k]
816     Y[idx_m[k]]+=inner(Yss, p_diffs[idx_f[j]][:,k])
817     else:
818     Y+=inner(Yss, p_diffs[idx_f[j]]) * mu
819     s+=l
820     if isinstance(vnames1, str):
821 jfenwick 5637 self.exportValue(vnames1, Y)
822 jfenwick 5563 else:
823     for n in vnames1:
824 jfenwick 5645 self.exportValue(n, Y)
825 jfenwick 5563 if isinstance(vnames2, str): #The second component should be strictly local
826 jfenwick 5637 self.exportValue(vnames2, g_J[1])
827 jfenwick 5563 else:
828     for n in vnames2:
829 jfenwick 5648 self.exportValue(n, g_J[1])
830     # End CalculateGradientWorker
831 jfenwick 5637
832     addJobPerWorld(self.splitworld, FunctionJob, calculateGradientWorker, vnames1=vnames1, vnames2=vnames2, imports=["models", "regularization", "props", "mu_models", "current_point"])
833 jfenwick 5635 self.splitworld.runJobs()
834 jfenwick 5562
835     def _getGradient(self, m, *args):
836     """
837     returns the gradient of the cost function at *m*.
838     If the pre-computed values are not supplied `getArguments()` is called.
839    
840     :param m: current approximation of the level set function
841     :type m: `Data`
842     :param args: tuple of values of the parameters, pre-computed values
843     for the forward model and pre-computed values for the
844     regularization
845    
846     :rtype: `ArithmeticTuple`
847    
848     :note: returns (Y^,X) where Y^ is the gradient from regularization plus
849     gradients of fwd models. X is the gradient of the regularization
850     w.r.t. gradient of m.
851     """
852     if not self.configured:
853 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
854 jfenwick 5562 raise RuntimeError("Call to getGradient -- temporary block to see where this is used")
855     if len(args)==0:
856     args = self.getArguments(m)
857    
858     props=args[0]
859     args_f=args[1]
860     args_reg=args[2]
861    
862     g_J = self.regularization.getGradient(m, *args_reg)
863     p_diffs=[]
864     # Find the derivative for each mapping
865     # If a mapping has a list of components (idx), then make a new Data object with only those
866     # components, pass it to the mapping and get the derivative.
867     for i in range(self.numMappings):
868     mm, idx=self.mappings[i]
869     if idx and self.numLevelSets > 1:
870     if len(idx)>1:
871     m2=Data(0,(len(idx),),m.getFunctionSpace())
872     for k in range(len(idx)): m2[k]=m[idx[k]]
873     dpdm = mm.getDerivative(m2)
874     else:
875     dpdm = mm.getDerivative(m[idx[0]])
876     else:
877     dpdm = mm.getDerivative(m)
878     p_diffs.append(dpdm)
879    
880     Y=g_J[0] # Because g_J==(Y,X) Y_k=dKer/dm_k
881     for i in range(self.numModels):
882     mu=self.mu_model[i]
883     f, idx_f=self.forward_models[i]
884     args=tuple( [ props[k] for k in idx_f] + list( args_f[i] ) )
885     Ys = f.getGradient(*args) # this d Jf/d props
886     # in this case f depends on one parameter props only but this can
887     # still depend on several level set components
888     if Ys.getRank() == 0:
889     # run through all level sets k prop j is depending on:
890     idx_m=self.mappings[idx_f[0]][1]
891     # tmp[k] = dJ_f/d_prop * d prop/d m[idx_m[k]]
892     tmp=Ys * p_diffs[idx_f[0]] * mu
893     if idx_m:
894     if tmp.getRank()== 0:
895     for k in range(len(idx_m)):
896     Y[idx_m[k]]+=tmp # dJ_f /d m[idx_m[k]] = tmp
897     else:
898     for k in range(len(idx_m)):
899     Y[idx_m[k]]+=tmp[k] # dJ_f /d m[idx_m[k]] = tmp[k]
900     else:
901     Y+=tmp # dJ_f /d m[idx_m[k]] = tmp
902     else:
903     s=0
904     # run through all props j forward model f is depending on:
905     for j in range(len(idx_f)):
906     # run through all level sets k prop j is depending on:
907     idx_m=self.mappings[j][1]
908     if p_diffs[idx_f[j]].getRank() == 0 :
909     if idx_m: # this case is not needed (really?)
910     self.logger.error("something wrong A")
911     # tmp[k] = dJ_f/d_prop[j] * d prop[j]/d m[idx_m[k]]
912     tmp=Ys[s]*p_diffs[idx_f[j]] * mu
913     for k in range(len(idx_m)):
914     Y[idx_m[k]]+=tmp[k] # dJ_f /d m[idx_m[k]] = tmp[k]
915     else:
916     Y+=Ys[s]*p_diffs[idx_f[j]] * mu
917     s+=1
918     elif p_diffs[idx_f[j]].getRank() == 1 :
919     l=p_diffs[idx_f[j]].getShape()[0]
920     # tmp[k]=sum_j dJ_f/d_prop[j] * d prop[j]/d m[idx_m[k]]
921     tmp=inner(Ys[s:s+l], p_diffs[idx_f[j]]) * mu
922     if idx_m:
923     for k in range(len(idx_m)):
924     Y[idx_m[k]]+=tmp # dJ_f /d m[idx_m[k]] = tmp[k]
925     else:
926     Y+=tmp
927     s+=l
928     else: # rank 2 case
929     l=p_diffs[idx_f[j]].getShape()[0]
930     Yss=Ys[s:s+l]
931     if idx_m:
932     for k in range(len(idx_m)):
933     # dJ_f /d m[idx_m[k]] = tmp[k]
934     Y[idx_m[k]]+=inner(Yss, p_diffs[idx_f[j]][:,k])
935     else:
936     Y+=inner(Yss, p_diffs[idx_f[j]]) * mu
937     s+=l
938     return g_J
939    
940     def _getInverseHessianApproximation(self, m, r, *args):
941     """
942     returns an approximative evaluation *p* of the inverse of the Hessian
943     operator of the cost function for a given gradient type *r* at a
944     given location *m*: *H(m) p = r*
945    
946     :param m: level set approximation where to calculate Hessian inverse
947     :type m: `Data`
948     :param r: a given gradient
949     :type r: `ArithmeticTuple`
950     :param args: tuple of values of the parameters, pre-computed values
951     for the forward model and pre-computed values for the
952     regularization
953     :rtype: `Data`
954     :note: in the current implementation only the regularization term is
955     considered in the inverse Hessian approximation.
956     """
957     if not self.configured:
958 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
959 jfenwick 5562 raise RuntimeError("Call to getInverseHessianApproximation -- temporary block to see where this is used")
960    
961     m=self.regularization.getInverseHessianApproximation(m, r, *args[2])
962     return m
963    
964     def updateHessian(self):
965     """
966     notifies the class that the Hessian operator needs to be updated.
967     """
968     if not self.configured:
969 jfenwick 5630 raise ValueError("This inversion function has not been configured yet")
970     addJobPerWorld(self.splitworld, FunctionJob, updateHessianWorker, imports=["regularization"])
971     self.splitworld.runJobs()
972 jfenwick 5562
973     def _getNorm(self, m):
974     """
975     returns the norm of `m`
976    
977     :param m: level set function
978     :type m: `Data`
979     :rtype: ``float``
980     """
981     if not self.configured:
982 jfenwick 5563 raise ValueError("This inversion function has not been configured yet")
983     raise RuntimeError("Need to have this in a subworld --- one or all?")
984 jfenwick 5562 return self.regularization.getNorm(m)
985    
986 jfenwick 5630 def updateHessianWorker(self, **kwargs):
987     reg=self.importValue("regularization")
988     reg.updateHessian()
989     #self.exportValue(reg, "regularization")

  ViewVC Help
Powered by ViewVC 1.1.26