downunder/py_src/inversions.py


##############################################################################
#
# Copyright (c) 2003-2012 by University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development since 2012 by School of Earth Sciences
#
##############################################################################

"""Higher-level classes that allow running inversions with minimal set up"""

__copyright__="""Copyright (c) 2003-2012 by University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="https://launchpad.net/escript-finley"

__all__ = ['InversionBase','SingleParameterInversionBase','GravityInversion','MagneticInversion']

import logging
from esys.escript import *
import esys.escript.unitsSI as U
from esys.weipa import createDataset

from .inversioncostfunctions import SimpleInversionCostFunction
from .forwardmodels import GravityModel, MagneticModel
from .mappings import *
from .minimizers import *
from .regularizations import Regularization
from .datasources import DataSource

class InversionBase(object):
    """
    Base class for inversions.
    """
    def __init__(self):
        self.logger=logging.getLogger('inv.%s'%self.__class__.__name__)
        self._solver_callback = None
        self._solver_opts = {}
        self._solver_tol = 1e-9
        self._solver_maxiter = 200
        # set default mapping and solver
        self.setSolverClass()
        self.__domain=None
        self.__regularization=None
        self.__forwardmodel=None
        self.__mapping=None
        
    def setDomain(self, domain):
        """
        sets the domain of the inversion

        :param domain: domain of the inversion
        :type domain: `Domain`
        """
        self.__domain=domain

    def getDomain(self):
        """
        returns the domain of the inversion

        :rtype: `Domain`
        """
        return  self.__domain

    def setMapping(self, mapping):
        """
        Sets the mapping object to map between model parameters and the data.
        If no mapping is provided, an identity mapping is used
        (`ScalingMapping` with constant 1).

        :param mapping: parameter mapping
        :type mapping: `Mapping`
        """
        self.__mapping=mapping

    def getMapping(self):
        """
        return the mapping(s) used in the inversion

        :rtype: `Mapping`
        """
        return self.__mapping


    def setRegularization(self, regularization):
        """
        Sets the regularization for the inversion.

        :param regularization: regularization
        :type regularization: `Regularization`
        """
        self.__regularization=regularization

    def getRegularization(self):
        """
        returns the regularization method(s)

        :rtype: `Regularization`
        """
        return self.__regularization

    def setForwardModel(self, forwardmodel):
        """
        Sets the forward model(s) for the inversion.

        :param forwardmodel: forward model
        :type forwardmodel: `ForwardModel`
        """
        self.__forwardmodel=forwardmodel

    def getForwardModel(self):
        """
        returns the forward model

        :rtype: `ForwardModel`
        """
        return self.__forwardmodel

    def setSolverClass(self, solverclass=None):
        """
        The solver to be used in the inversion process. See the minimizers
        module for available solvers.
        By default, the L-BFGS minimizer is used.
        """
        if solverclass == None:
            self.solverclass=MinimizerLBFGS
        else:
            self.solverclass=solverclass

    def setSolverCallback(self, callback):
        """
        Sets the callback function which is called after every solver
        iteration.
        """
        self._solver_callback=callback

    def setSolverMaxIterations(self, maxiter):
        """
        Sets the maximum number of solver iterations to run.
        """
        if maxiter>0:
            self._solver_maxiter=maxiter
        else:
            raise ValueError("maxiter must be positive.")

    def setSolverOptions(self, **opts):
        """
        Sets additional solver options. The valid options depend on the
        solver being used.
        """
        self._solver_opts.update(**opts)

    def setSolverTolerance(self, tol):
        """
        Sets the error tolerance for the solver. An acceptable solution is
        considered to be found once the tolerance is reached.
        """
        if tol>0:
            self._solver_tol=tol
        else:
            raise ValueError("tolerance must be positive.")

    def isSetup(self):
        """
        returns True if the inversion is set up and ready to run.
        """
        raise NotImplementedError

    def run(self, *args):
        """
        This method starts the inversion process and must be overwritten.
        Preferably, users should be able to set a starting point so
        inversions can be continued after stopping.
        """
        raise NotImplementedError


class SingleParameterInversionBase(InversionBase):
    """
    Base class for inversions with a single parameter to be found.
    """
    def __init__(self):
        super(SingleParameterInversionBase,self).__init__()
        self.setWeights()

    def setWeights(self, mu_reg=None, mu_model=1.):
        """
        Sets the weighting factors for the cost function.
        """
        self.logger.debug("Setting weighting factors...")
        self.logger.debug("mu_reg = %s"%mu_reg)
        self.logger.debug("mu_model = %s"%mu_model)
        self._mu_reg=mu_reg
        self._mu_model=mu_model

    def siloWriterCallback(self, k, x, fx, gfx):
        """
        callback function that can be used to track the solution

        :param k: iteration count
        :param x: current m approximation
        :param fx: value of cost function
        :param gfx: gradient of f at x
        """
        fn='inv.%d'%k
        ds=createDataset(rho=self.getMapping().getValue(x))
        ds.setCycleAndTime(k,k)
        ds.saveSilo(fn)
        self.logger.debug("Jreg(m) = %e"%self.regularization.getValue(x))
        self.logger.debug("f(m) = %e"%fx)


    def isSetup(self):
        if self.getRegularization() and self.getMapping() \
                and self.getForwardModel() and self.getDomain():
            return True
        else:
            return False

    def run(self, initial_value=0.):
        if not self.isSetup():
            raise RuntimeError("Inversion is not setup properly.")

        f=SimpleInversionCostFunction(self.getRegularization(), self.getMapping(), self.getForwardModel())
        f.setWeights(mu_reg_1=self._mu_reg, mu_model=self._mu_model)

        solver=self.solverclass(f)
        solver.setCallback(self._solver_callback)
        solver.setMaxIterations(self._solver_maxiter)
        solver.setOptions(**self._solver_opts)
        solver.setTolerance(self._solver_tol)
        if not isinstance(initial_value, Data):
            initial_value=Scalar(initial_value, ContinuousFunction(self.getDomain()))
        m_init=self.getMapping().getInverse(initial_value)

        self.logger.info("Starting solver...")
        solver.run(m_init)
        m_star=solver.getResult()
        self.logger.info("m* = %s"%m_star)
        value_star=self.getMapping().getValue(m_star)
        self.logger.info("result* = %s"%value_star)
        solver.logSummary()
        return value_star


class GravityInversion(SingleParameterInversionBase):
    """
    Inversion of Gravity (Bouguer) anomaly data.
    """
    def setup(self, domainbuilder, 
                    rho0=None, drho=None, z0=None, beta=None, 
                    s0=None, s1=None):
        """
        Sets up the inversion parameters from a downunder.`DomainBuilder`.

        :param domainbuilder: Domain builder object with gravity source(s)
        :type domainbuilder: `DomainBuilder`
        
        
        :param rho0: reference density. If not specified, zero is used.
        :type rho0: ``float`` or `Scalar`
        :param w0: weighting factor for L2 term in the regularization
        :type w0: ``float`` or `Scalar`
        :param w: weighting factor for H1 term in the regularization
        :type w: ``list`` of float or `Vector`
        """
    self.logger.info('Retrieving domain...')
    self.setDomain(domainbuilder.getDomain())
    DIM=self.getDomain().getDim()
    #========================        
    self.logger.info('Set mapping ...')
    self.setMapping(DensityMapping(self.getDomain(), rho0=rho0, drho=drho, z0=z0, beta=beta))
    #========================
    self.logger.info("Setting up regularization...")
    if s1 is None: 
        s1=[1.]*DIM
    rho_mask = domainbuilder.getSetDensityMask()
    self.setRegularization(Regularization(self.getDomain(), numLevelSets=1, \
                   s0=s0, s1=s1, location_of_set_m=rho_mask))
    #====================================================================
    self.logger.info("Retrieving gravity surveys...")
    surveys=domainbuilder.getGravitySurveys()
    g=[]
    chi=[]
    for g_i,sigma_i in surveys:
        chi_i=safeDiv(1., sigma_i*sigma_i)
        if g_i.getRank()==0:
        g_i=g_i*kronecker(DIM)[DIM-1]
        if chi_i.getRank()==0:
        chi_i=chi_i*kronecker(DIM)[DIM-1]
        g.append(g_i)
        chi.append(chi_i)
        self.logger.debug("Added gravity survey:")
        self.logger.debug("g = %s"%g_i)
        self.logger.debug("sigma = %s"%sigma_i)
        self.logger.debug("chi = %s"%chi_i)
    #====================================================================
    
    self.logger.info("Setting up model...")
    self.setForwardModel(GravityModel(self.getDomain(), chi, g))

    # this is switched off for now:
    if self._mu_reg is None and False:
        print "SDSAD"
        x=self.getDomain().getX()
        l=0.
        for i in range(DIM-1):
        l=max(l, sup(x[i])-inf(x[i]))
        G=U.Gravitational_Constant
        mu_reg=0.5*(l*l*G)**2
        self.setWeights(mu_reg=mu_reg)


class MagneticInversion(SingleParameterInversionBase):
    """
    Inversion of magnetic data.
    """
    def setup(self, domainbuilder, 
                    k0=None, dk=None, z0=None, beta=None,
                    s0=None, s1=None):
        """
        Sets up the inversion parameters from a downunder.`DomainBuilder`.

        :param domainbuilder: Domain builder object with magnetic source(s)
        :type domainbuilder: `DomainBuilder`
        
        :param k_ref: reference susceptibility. If not specified, zero is used.
        :type k_ref: ``float`` or `Scalar`
        :param w0: weighting factor for the L2 term in the regularization
        :type w0: ``float`` or `Scalar`
        :param w: weighting factor for the H1 term in the regularization
        :type w: ``list`` of float or `Vector`
        """
        self.logger.info('Retrieving domain...')
        self.setDomain(domainbuilder.getDomain())
        DIM=self.getDomain().getDim()
        
        #========================
        self.logger.info('Set mapping ...')
        self.setMapping(SusceptibilityMapping(self.getDomain(), k0=k0, dk=dk, z0=z0, beta=beta))
        
        #========================
        self.logger.info("Setting up regularization...")
        if s1 is None:
            s1=[1.]*DIM
        k_mask = domainbuilder.getSetSusceptibilityMask()
        self.setRegularization(Regularization(self.getDomain(), numLevelSets=1, \
                               s0=s0, s1=s1, location_of_set_m=k_mask))

        #====================================================================
        self.logger.info("Retrieving magnetic field surveys...")
        surveys=domainbuilder.getMagneticSurveys()
        B=[]
        chi=[]
        for B_i,sigma_i in surveys:
            chi_i=safeDiv(1., sigma_i*sigma_i)
            if B_i.getRank()==0:
                B_i=B_i*kronecker(DIM)[DIM-1]
            if chi_i.getRank()==0:
                chi_i=chi_i*kronecker(DIM)[DIM-1]
            B.append(B_i)
            chi.append(chi_i)
            self.logger.debug("Added magnetic survey:")
            self.logger.debug("B = %s"%B_i)
            self.logger.debug("sigma = %s"%sigma_i)
            self.logger.debug("chi = %s"%chi_i)
        #====================================================================
        self.logger.info("Setting up model...")
        self.setForwardModel(MagneticModel(self.getDomain(), chi, B, domainbuilder.getBackgroundMagneticField()))

        # this is switched off for now:
        if self._mu_reg is None and False:
            x=self.getDomain().getX()
            l=0.
            for i in range(DIM-1):
                l=max(l, sup(x[i])-inf(x[i]))
            mu_reg=0.5*l**2
            self.setWeights(mu_reg=mu_reg)

1
2	##############################################################################
3	#
4	# Copyright (c) 2003-2012 by University of Queensland
5	# 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 since 2012 by School of Earth Sciences
13	#
14	##############################################################################
15
16	"""Higher-level classes that allow running inversions with minimal set up"""
17
18	__copyright__="""Copyright (c) 2003-2012 by University of Queensland
19	http://www.uq.edu.au
20	Primary Business: Queensland, Australia"""
21	__license__="""Licensed under the Open Software License version 3.0
22	http://www.opensource.org/licenses/osl-3.0.php"""
23	__url__="https://launchpad.net/escript-finley"
24
25	__all__ = ['InversionBase','SingleParameterInversionBase','GravityInversion','MagneticInversion']
26
27	import logging
28	from esys.escript import *
29	import esys.escript.unitsSI as U
30	from esys.weipa import createDataset
31
32	from .inversioncostfunctions import SimpleInversionCostFunction
33	from .forwardmodels import GravityModel, MagneticModel
34	from .mappings import *
35	from .minimizers import *
36	from .regularizations import Regularization
37	from .datasources import DataSource
38
39	class InversionBase(object):
40	"""
41	Base class for inversions.
42	"""
43	def __init__(self):
44	self.logger=logging.getLogger('inv.%s'%self.__class__.__name__)
45	self._solver_callback = None
46	self._solver_opts = {}
47	self._solver_tol = 1e-9
48	self._solver_maxiter = 200
49	# set default mapping and solver
50	self.setSolverClass()
51	self.__domain=None
52	self.__regularization=None
53	self.__forwardmodel=None
54	self.__mapping=None
55
56	def setDomain(self, domain):
57	"""
58	sets the domain of the inversion
59
60	:param domain: domain of the inversion
61	:type domain: `Domain`
62	"""
63	self.__domain=domain
64
65	def getDomain(self):
66	"""
67	returns the domain of the inversion
68
69	:rtype: `Domain`
70	"""
71	return self.__domain
72
73	def setMapping(self, mapping):
74	"""
75	Sets the mapping object to map between model parameters and the data.
76	If no mapping is provided, an identity mapping is used
77	(`ScalingMapping` with constant 1).
78
79	:param mapping: parameter mapping
80	:type mapping: `Mapping`
81	"""
82	self.__mapping=mapping
83
84	def getMapping(self):
85	"""
86	return the mapping(s) used in the inversion
87
88	:rtype: `Mapping`
89	"""
90	return self.__mapping
91
92
93	def setRegularization(self, regularization):
94	"""
95	Sets the regularization for the inversion.
96
97	:param regularization: regularization
98	:type regularization: `Regularization`
99	"""
100	self.__regularization=regularization
101
102	def getRegularization(self):
103	"""
104	returns the regularization method(s)
105
106	:rtype: `Regularization`
107	"""
108	return self.__regularization
109
110	def setForwardModel(self, forwardmodel):
111	"""
112	Sets the forward model(s) for the inversion.
113
114	:param forwardmodel: forward model
115	:type forwardmodel: `ForwardModel`
116	"""
117	self.__forwardmodel=forwardmodel
118
119	def getForwardModel(self):
120	"""
121	returns the forward model
122
123	:rtype: `ForwardModel`
124	"""
125	return self.__forwardmodel
126
127	def setSolverClass(self, solverclass=None):
128	"""
129	The solver to be used in the inversion process. See the minimizers
130	module for available solvers.
131	By default, the L-BFGS minimizer is used.
132	"""
133	if solverclass == None:
134	self.solverclass=MinimizerLBFGS
135	else:
136	self.solverclass=solverclass
137
138	def setSolverCallback(self, callback):
139	"""
140	Sets the callback function which is called after every solver
141	iteration.
142	"""
143	self._solver_callback=callback
144
145	def setSolverMaxIterations(self, maxiter):
146	"""
147	Sets the maximum number of solver iterations to run.
148	"""
149	if maxiter>0:
150	self._solver_maxiter=maxiter
151	else:
152	raise ValueError("maxiter must be positive.")
153
154	def setSolverOptions(self, **opts):
155	"""
156	Sets additional solver options. The valid options depend on the
157	solver being used.
158	"""
159	self._solver_opts.update(**opts)
160
161	def setSolverTolerance(self, tol):
162	"""
163	Sets the error tolerance for the solver. An acceptable solution is
164	considered to be found once the tolerance is reached.
165	"""
166	if tol>0:
167	self._solver_tol=tol
168	else:
169	raise ValueError("tolerance must be positive.")
170
171	def isSetup(self):
172	"""
173	returns True if the inversion is set up and ready to run.
174	"""
175	raise NotImplementedError
176
177	def run(self, *args):
178	"""
179	This method starts the inversion process and must be overwritten.
180	Preferably, users should be able to set a starting point so
181	inversions can be continued after stopping.
182	"""
183	raise NotImplementedError
184
185
186	class SingleParameterInversionBase(InversionBase):
187	"""
188	Base class for inversions with a single parameter to be found.
189	"""
190	def __init__(self):
191	super(SingleParameterInversionBase,self).__init__()
192	self.setWeights()
193
194	def setWeights(self, mu_reg=None, mu_model=1.):
195	"""
196	Sets the weighting factors for the cost function.
197	"""
198	self.logger.debug("Setting weighting factors...")
199	self.logger.debug("mu_reg = %s"%mu_reg)
200	self.logger.debug("mu_model = %s"%mu_model)
201	self._mu_reg=mu_reg
202	self._mu_model=mu_model
203
204	def siloWriterCallback(self, k, x, fx, gfx):
205	"""
206	callback function that can be used to track the solution
207
208	:param k: iteration count
209	:param x: current m approximation
210	:param fx: value of cost function
211	:param gfx: gradient of f at x
212	"""
213	fn='inv.%d'%k
214	ds=createDataset(rho=self.getMapping().getValue(x))
215	ds.setCycleAndTime(k,k)
216	ds.saveSilo(fn)
217	self.logger.debug("Jreg(m) = %e"%self.regularization.getValue(x))
218	self.logger.debug("f(m) = %e"%fx)
219
220
221	def isSetup(self):
222	if self.getRegularization() and self.getMapping() \
223	and self.getForwardModel() and self.getDomain():
224	return True
225	else:
226	return False
227
228	def run(self, initial_value=0.):
229	if not self.isSetup():
230	raise RuntimeError("Inversion is not setup properly.")
231
232	f=SimpleInversionCostFunction(self.getRegularization(), self.getMapping(), self.getForwardModel())
233	f.setWeights(mu_reg_1=self._mu_reg, mu_model=self._mu_model)
234
235	solver=self.solverclass(f)
236	solver.setCallback(self._solver_callback)
237	solver.setMaxIterations(self._solver_maxiter)
238	solver.setOptions(**self._solver_opts)
239	solver.setTolerance(self._solver_tol)
240	if not isinstance(initial_value, Data):
241	initial_value=Scalar(initial_value, ContinuousFunction(self.getDomain()))
242	m_init=self.getMapping().getInverse(initial_value)
243
244	self.logger.info("Starting solver...")
245	solver.run(m_init)
246	m_star=solver.getResult()
247	self.logger.info("m* = %s"%m_star)
248	value_star=self.getMapping().getValue(m_star)
249	self.logger.info("result* = %s"%value_star)
250	solver.logSummary()
251	return value_star
252
253
254	class GravityInversion(SingleParameterInversionBase):
255	"""
256	Inversion of Gravity (Bouguer) anomaly data.
257	"""
258	def setup(self, domainbuilder,
259	rho0=None, drho=None, z0=None, beta=None,
260	s0=None, s1=None):
261	"""
262	Sets up the inversion parameters from a downunder.`DomainBuilder`.
263
264	:param domainbuilder: Domain builder object with gravity source(s)
265	:type domainbuilder: `DomainBuilder`
266
267
268	:param rho0: reference density. If not specified, zero is used.
269	:type rho0: ``float`` or `Scalar`
270	:param w0: weighting factor for L2 term in the regularization
271	:type w0: ``float`` or `Scalar`
272	:param w: weighting factor for H1 term in the regularization
273	:type w: ``list`` of float or `Vector`
274	"""
275	self.logger.info('Retrieving domain...')
276	self.setDomain(domainbuilder.getDomain())
277	DIM=self.getDomain().getDim()
278	#========================
279	self.logger.info('Set mapping ...')
280	self.setMapping(DensityMapping(self.getDomain(), rho0=rho0, drho=drho, z0=z0, beta=beta))
281	#========================
282	self.logger.info("Setting up regularization...")
283	if s1 is None:
284	s1=[1.]*DIM
285	rho_mask = domainbuilder.getSetDensityMask()
286	self.setRegularization(Regularization(self.getDomain(), numLevelSets=1, \
287	s0=s0, s1=s1, location_of_set_m=rho_mask))
288	#====================================================================
289	self.logger.info("Retrieving gravity surveys...")
290	surveys=domainbuilder.getGravitySurveys()
291	g=[]
292	chi=[]
293	for g_i,sigma_i in surveys:
294	chi_i=safeDiv(1., sigma_i*sigma_i)
295	if g_i.getRank()==0:
296	g_i=g_i*kronecker(DIM)[DIM-1]
297	if chi_i.getRank()==0:
298	chi_i=chi_i*kronecker(DIM)[DIM-1]
299	g.append(g_i)
300	chi.append(chi_i)
301	self.logger.debug("Added gravity survey:")
302	self.logger.debug("g = %s"%g_i)
303	self.logger.debug("sigma = %s"%sigma_i)
304	self.logger.debug("chi = %s"%chi_i)
305	#====================================================================
306
307	self.logger.info("Setting up model...")
308	self.setForwardModel(GravityModel(self.getDomain(), chi, g))
309
310	# this is switched off for now:
311	if self._mu_reg is None and False:
312	print "SDSAD"
313	x=self.getDomain().getX()
314	l=0.
315	for i in range(DIM-1):
316	l=max(l, sup(x[i])-inf(x[i]))
317	G=U.Gravitational_Constant
318	mu_reg=0.5(llG)*2
319	self.setWeights(mu_reg=mu_reg)
320
321
322	class MagneticInversion(SingleParameterInversionBase):
323	"""
324	Inversion of magnetic data.
325	"""
326	def setup(self, domainbuilder,
327	k0=None, dk=None, z0=None, beta=None,
328	s0=None, s1=None):
329	"""
330	Sets up the inversion parameters from a downunder.`DomainBuilder`.
331
332	:param domainbuilder: Domain builder object with magnetic source(s)
333	:type domainbuilder: `DomainBuilder`
334
335	:param k_ref: reference susceptibility. If not specified, zero is used.
336	:type k_ref: ``float`` or `Scalar`
337	:param w0: weighting factor for the L2 term in the regularization
338	:type w0: ``float`` or `Scalar`
339	:param w: weighting factor for the H1 term in the regularization
340	:type w: ``list`` of float or `Vector`
341	"""
342	self.logger.info('Retrieving domain...')
343	self.setDomain(domainbuilder.getDomain())
344	DIM=self.getDomain().getDim()
345
346	#========================
347	self.logger.info('Set mapping ...')
348	self.setMapping(SusceptibilityMapping(self.getDomain(), k0=k0, dk=dk, z0=z0, beta=beta))
349
350	#========================
351	self.logger.info("Setting up regularization...")
352	if s1 is None:
353	s1=[1.]*DIM
354	k_mask = domainbuilder.getSetSusceptibilityMask()
355	self.setRegularization(Regularization(self.getDomain(), numLevelSets=1, \
356	s0=s0, s1=s1, location_of_set_m=k_mask))
357
358	#====================================================================
359	self.logger.info("Retrieving magnetic field surveys...")
360	surveys=domainbuilder.getMagneticSurveys()
361	B=[]
362	chi=[]
363	for B_i,sigma_i in surveys:
364	chi_i=safeDiv(1., sigma_i*sigma_i)
365	if B_i.getRank()==0:
366	B_i=B_i*kronecker(DIM)[DIM-1]
367	if chi_i.getRank()==0:
368	chi_i=chi_i*kronecker(DIM)[DIM-1]
369	B.append(B_i)
370	chi.append(chi_i)
371	self.logger.debug("Added magnetic survey:")
372	self.logger.debug("B = %s"%B_i)
373	self.logger.debug("sigma = %s"%sigma_i)
374	self.logger.debug("chi = %s"%chi_i)
375	#====================================================================
376	self.logger.info("Setting up model...")
377	self.setForwardModel(MagneticModel(self.getDomain(), chi, B, domainbuilder.getBackgroundMagneticField()))
378
379	# this is switched off for now:
380	if self._mu_reg is None and False:
381	x=self.getDomain().getX()
382	l=0.
383	for i in range(DIM-1):
384	l=max(l, sup(x[i])-inf(x[i]))
385	mu_reg=0.5l*2
386	self.setWeights(mu_reg=mu_reg)
387