agd.Eikonal.run_detail

  1# Copyright 2020 Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
  2# Distributed WITHOUT ANY WARRANTY. Licensed under the Apache License, Version 2.0, see http://www.apache.org/licenses/LICENSE-2.0
  3
  4import numpy as np
  5from .LibraryCall import RunDispatch,GetBinaryDir
  6from .DictIn_detail import factoringPointChoice
  7from .. import Metrics
  8from .. import AutomaticDifferentiation as ad
  9
 10def GetGeodesics(output,suffix=''): 
 11	if suffix != '' and not suffix.startswith('_'): suffix='_'+suffix
 12	geodesics = np.vsplit(output['geodesicPoints'+suffix],
 13					 output['geodesicLengths'+suffix].cumsum()[:-1].astype(int))
 14	return [geo.T for geo in geodesics]
 15
 16class Cache(object):
 17	def __init__(self,needsflow=False):
 18		self.contents = dict()
 19		self.verbosity = None
 20		self.requested = None
 21		self.needsflow = needsflow
 22		self.dummy = False
 23	
 24	def empty(self):
 25		"""Wether the cache lacks data needed to bypass computation"""
 26		return not self.contents 
 27
 28#	def full(self): 
 29#		return self.contents and ('geodesicFlow' in self.contents or not needsflow)
 30	
 31	def PreProcess(self,hfmIn_raw):
 32		if self.dummy: return
 33		self.verbosity = hfmIn_raw.get('verbosity',1)
 34		self.requested = []
 35		if hfmIn_raw.get('exportValues',False): 		self.requested.append('values')
 36		if hfmIn_raw.get('exportActiveNeighs',False):	self.requested.append('activeNeighs')
 37		if hfmIn_raw.get('exportGeodesicFlow',False):	self.requested.append('geodesicFlow')
 38		if self.empty():
 39			if self.verbosity>=1: print("Requesting cacheable data")
 40			hfmIn_raw['exportValues']=True
 41			hfmIn_raw['exportActiveNeighs']=True
 42			if self.needsflow: hfmIn_raw['exportGeodesicFlow']=True
 43		else:
 44			if self.verbosity>=1: print("Providing cached data")
 45			for key in ('values','activeNeighs'):
 46				setkey_safe(hfmIn_raw,key,self.contents[key])
 47			# Note : all points set as accepted in HFM algorithm
 48			hfmIn_raw['exportValues']=False
 49			hfmIn_raw['exportActiveNeighs']=False
 50			hfmIn_raw['exportGeodesicFlow']=self.needsflow and ('geodesicFlow' not in self.contents)
 51
 52	def PostProcess(self,hfmOut_raw):
 53		if self.dummy: return
 54		if self.empty():
 55			if self.verbosity>=1 : print("Filling cache data")
 56			for key in ('values','activeNeighs'):
 57				self.contents[key] = hfmOut_raw[key]
 58			if self.needsflow:
 59				self.contents['geodesicFlow'] = hfmOut_raw['geodesicFlow']
 60		else:
 61			for key in self.requested:
 62				if key not in hfmOut_raw:
 63					setkey_safe(hfmOut_raw,key,self.contents[key])
 64				if self.needsflow and 'geodesicFlow' not in self.contents:
 65					self.contents['geodesicFlow'] = hfmOut_raw['geodesicFlow']
 66
 67	def geodesicFlow(self,hfmIn=None):
 68		if 'geodesicFlow' in self.contents:
 69			return self.contents['geodesicFlow']
 70		elif hfmIn is None:
 71			raise ValueError("geodesicFlow not cached and lacking hfm input data")
 72		else:
 73			self.dummy = False
 74			self.needsflow = True
 75
 76			hfmIn_ = {key:ad.remove_ad(value,iterables=(Metrics.Base,)) for key,value in hfmIn.items()
 77			if key not in [ # Keys useless for evaluating the geodesic flow
 78			'tips','tips_Unoriented',
 79			'costVariation','seedValueVariation','inspectSensitivity',
 80			'exportActiveNeighs','exportActiveOffsets']}
 81
 82			hfmOut_raw = RunSmart(hfmIn_,cache=self,returns='out_raw')
 83			if self.verbosity: 
 84				print("--- HFM call triggered above to compute geodesic flow ---")
 85			self.PostProcess(hfmOut_raw)
 86			return self.contents['geodesicFlow']
 87
 88
 89def RunRaw(hfmIn):
 90	"""Raw call to the HFM library"""
 91	rawIn = hfmIn.copy()
 92	mode = rawIn.pop('mode',None); assert mode is None or mode=='cpu_raw'
 93	float_t = rawIn.pop('float_t',None); assert float_t is None or float_t==np.float64 
 94	return RunDispatch(rawIn,GetBinaryDir("FileHFM","HFMpy"))
 95
 96
 97def RunSmart(hfmIn,co_output=None,cache=None,returns="out"):
 98	"""
 99	Calls the HFM library, a CPU eikonal solver written in C++,
100	with pre-processing and post-processing of data.
101
102	Main input : 
103	- hfmIn, a dictionary like structure, containing the eikonal solver data.
104
105	The C++ library embeds some help information, which can be accessed using the 
106	following key:value pairs in hfmIn.
107	* 'verbosity' : set to 1 or 2 to display information on run, including the defaulted keys.
108	 set to 0 to silence the run.
109	* 'keyHelp' : set as a string of keys separated by spaces, to print help on these.
110	 example 'keyHelp':"seeds tips"
111
112	Optional input:
113	- co_output (optional) : used for reverse automatic differentiation.
114	- cache (optional) : store some intermediate results to bypass computations at a later stage. 
115	- returns (optional) : early aborts the run (debug purposes). ('in_raw','out_raw','out')
116	"""
117	assert returns in ('in_raw','out_raw','out')
118
119	if factoringPointChoice(hfmIn) == 'Both':
120		hfmIn = hfmIn.copy()
121		hfmIn.SetFactor()
122
123	hfmIn_raw = {}
124
125	if cache is None:
126		cache = Cache()
127		cache.dummy = True
128
129	# Pre-process usual arguments
130	for key,value in hfmIn.items():
131		PreProcess(key,value,hfmIn,hfmIn_raw,cache)
132
133	# Reverse automatic differentiation
134	if co_output is not None:
135		assert ad.misc.reverse_mode(co_output)=="Reverse"
136		[co_hfmOut,co_value],_ = co_output
137		assert hfmIn.get('extractValues',False) and co_hfmOut is None
138		indices = np.nonzero(co_value)
139		positions = hfmIn.PointFromIndex(np.asarray(indices).T)
140		weights = co_value[indices]
141		setkey_safe(hfmIn_raw,'inspectSensitivity',positions)
142		setkey_safe(hfmIn_raw,'inspectSensitivityWeights',weights)
143		setkey_safe(hfmIn_raw,'inspectSensitivityLengths',[len(weights)])
144		if 'metric' in hfmIn or 'dualMetric' in hfmIn:
145			cache.needsflow = True
146			cache.dummy = False
147	
148	# Dealing with cached data
149	cache.PreProcess(hfmIn_raw)
150
151	# Call to the HFM library
152	if returns=='in_raw': return hfmIn_raw
153	hfmOut_raw = RunDispatch(hfmIn_raw,GetBinaryDir("FileHFM","HFMpy"))
154	if returns=='out_raw': return hfmOut_raw
155	
156	# Dealing with cached data
157	cache.PostProcess(hfmOut_raw)
158
159	# Post process
160	hfmOut = {}
161	for key,val in hfmOut_raw.items():
162		PostProcess(key,val,hfmOut_raw,hfmOut)
163
164	# Reverse automatic differentiation
165	if co_output is not None:
166		result=[]
167		for key,value in hfmIn.items(): 
168			if key in ('cost','speed'):
169				if isinstance(value,Metrics.Isotropic):
170					value = value.to_HFM()
171				result.append((value,
172					hfmOut['costSensitivity_0'] if key=='cost' else (hfmOut['costSensitivity_0']/value**2)))
173			if key in ('metric','dualMetric'):
174				shape_bound = value.shape
175				considered = co_output.second
176				value_ad = ad.Dense.register(value,iterables=(Metrics.Base,),shape_bound=shape_bound,considered=considered)
177				metric_ad = value_ad if key=='metric' else value_ad.dual()
178
179				costSensitivity = np.moveaxis(flow_variation(cache.geodesicFlow(hfmIn),metric_ad),-1,0)*hfmOut['costSensitivity_0']
180				shift = 0
181				size_bound = np.prod(shape_bound,dtype=int)
182				for x in value:
183					if any(x is val for val in considered):
184						xsize_free = x.size//size_bound
185						result.append((x,costSensitivity[shift:(shift+xsize_free)].reshape(x.shape)) )
186						shift+=xsize_free
187			elif key=='seedValues':
188				sens = hfmOut['seedSensitivity_0']
189				# Match the seeds with their approx given in sensitivity
190				corresp = np.argmin(ad.Optimization.norm(np.expand_dims(hfmIn_raw['seeds'],axis=2)
191					-np.expand_dims(sens[:,:2].T,axis=1),ord=2,axis=0),axis=0)
192				sens_ordered = np.zeros_like(value)
193				sens_ordered[corresp]=sens[:,2]
194				result.append((value,sens_ordered))
195			# TODO : speed
196		return result
197
198	return (hfmOut,hfmOut.pop('values')) if hfmIn.get('extractValues',False) else hfmOut	
199
200def setkey_safe(dico,key,value):
201	if key in dico:
202		if value is not dico[key]:
203			raise ValueError("Multiple values for key ",key)
204	else:
205		dico[key]=value
206
207def flow_variation(flow,metric):
208	flow = np.moveaxis(flow,-1,0)
209	zeros = np.all(flow==0.,axis=0)
210	flow.__setitem__((slice(None),zeros),np.nan)
211	norm = metric.norm(flow)
212	variation = norm.coef/np.expand_dims(norm.value,axis=-1)
213	variation.__setitem__((zeros,slice(None)),0.)
214	return variation
215
216# ----------------- Preprocessing ---------------
217def PreProcess(key,value,refined_in,raw_out,cache):
218	"""
219	copies key,val from refined to raw, with adequate treatment
220	"""
221
222	verbosity = refined_in.get('verbosity',1)
223	if   key=='mode': assert value in ('cpu','cpu_transfer'); return
224	elif key=='float_t': assert value==np.float64; return
225	elif key in ('cost','speed'):
226		if isinstance(value,Metrics.Isotropic):
227			value = value.to_HFM()
228		if isinstance(value,ad.Dense.denseAD):
229			setkey_safe(raw_out,'costVariation',
230				value.coef if key=='cost' else (1/value).coef)
231			value = ad.remove_ad(value)
232		setkey_safe(raw_out,key,value)
233	elif key in ('metric','dualMetric'):
234		if refined_in['model'].startswith('Isotropic') or refined_in['model'].startswith('Diagonal'):
235			assert isinstance(value,(Metrics.Isotropic,Metrics.Diagonal))
236			return PreProcess('cost' if key=='metric' else 'speed',value.to_HFM(),refined_in,raw_out,cache)
237		if isinstance(value,Metrics.Base): 
238			if ad.is_ad(value,iterables=(Metrics.Base,)):
239				metric_ad = value if key=='metric' else value.dual()
240				setkey_safe(raw_out,'costVariation',flow_variation(cache.geodesicFlow(refined_in),metric_ad))
241			value = ad.remove_ad(value.to_HFM())
242		setkey_safe(raw_out,key,value)
243
244	elif key=='seedValues':
245		if ad.is_ad(value):
246			setkey_safe(raw_out,'seedValueVariation',value.gradient())
247			value=ad.remove_ad(value)
248		setkey_safe(raw_out,key,value)
249
250	elif key=='extractValues':
251		setkey_safe(raw_out,'exportValues',value)
252	else:
253		setkey_safe(raw_out,key,value)
254
255#---------- Post processing ------------
256def PostProcess(key,value,raw_in,refined_out):
257	"""
258	copies key,val from raw to refined, with adequate treatment
259	"""
260	if key.startswith('geodesicPoints'):
261		suffix = key[len('geodesicPoints'):]
262		geodesics = GetGeodesics(raw_in,suffix=suffix)
263		setkey_safe(refined_out,"geodesics"+suffix,geodesics)
264	elif key.startswith('geodesicLengths'):
265		pass
266
267	elif key=='values':
268		if 'valueVariation' in raw_in:
269			value = ad.Dense.denseAD(value,raw_in['valueVariation'])
270		setkey_safe(refined_out,key,value)
271	elif key=='valueVariation':
272		pass
273	
274	elif key=='geodesicFlow':
275		setkey_safe(refined_out,'flow',np.moveaxis(value,-1,0))
276	elif key=='factoringGradients':
277		setkey_safe(refined_out,'factoringGradients',np.moveaxis(value,-1,0))
278	elif key=='activeOffsets':
279		setkey_safe(refined_out,'offsets',CastOffsets(value))
280	else:
281		setkey_safe(refined_out,key,value)
282
283def CastOffsets(raw_offsets):
284	"""Offsets are exported with their coefficients bundled together in a double.
285	This function unpacks them."""
286	raw_shape = raw_offsets.shape
287	d = len(raw_shape)-1
288	nOffsets = raw_shape[-1]
289	refined_shape = (d,nOffsets)+raw_shape[:-1]
290	refined_offsets = np.zeros(refined_shape,dtype=int)
291	def coef(x,i):
292		x = x//256**(d-i-1)
293		x = x%256
294		return np.where(x<128,x,x-256)
295
296	for j in range(nOffsets):
297		raw = raw_offsets[...,j].astype(np.int64)
298		for i in range(d):
299			refined_offsets[i,j] = coef(raw,i)
300	return refined_offsets