Commit 4ac1cb92 authored by bepinat's avatar bepinat

additionnal convergence options

parent 8670d89b
......@@ -460,12 +460,15 @@ def elspectrum(p, wave=None, lines=None, psf=None, degc=None):
return (continuum + slines)
def myelspectrum(p, fjac=None, wave=None, spectrum= None, err=None, lines=None, psf=None, degc=None):
model = elspectrum(p, wave=wave, lines=lines, psf=psf, degc=degc)
def myelspectrum(p, fjac=None, wave=None, spectrum=None, err=None, lines=None, psf=None, degc=None):
#t1=time.time()
model = elspectrum(p, wave=wave, lines=lines, psf=psf, degc=degc)
status = 0
#t2=time.time()
#print(t2-t1)
return [status, (spectrum - model) / err]
def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=False):
def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=False, factor=10, ftol=1e-7, gtol=1e-7, xtol=1e-7, maxiter=200):
"""
"""
......@@ -554,11 +557,11 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#Redshift step
zstep = (conf['DFIT'] * 1e3 / ct.c + zcosmo) / (1 - conf['DFIT'] * 1e3 / ct.c) - zcosmo
if zstep >= (zmax - zmin) / 2: zstep = (zmax - zmin) / 2.
#if zstep >= (zmax - zmin) / 2: zstep = (zmax - zmin) / 2.
zmean = (zmax + zmin) / 2.
zzmax = zmean + (np.round(zmean / zstep) + 1) * zstep
zzmin = zmean - np.round(zmean / zstep) * zstep
zzmax = zmean + np.floor(zmax / zstep + 1) * zstep
zzmin = zmean + np.ceil(zmin / zstep) * zstep
#Parameters contraints and informations
params = np.zeros(1 + els.index * 3 + conf['DGCTNUM'] + 1)
parbase = {'value':0., 'fixed':1, 'limited':[0,0], 'limits':[0.,0.], 'tied':'', 'mpmaxstep':0, 'mpminstep':0, 'parname':''}
......@@ -687,6 +690,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
counter = 0
fac_prev = -1
number_of_pixels = np.float(cube.shape[2] * cube.shape[1])
tii = time.time()
tt = 0
for x in range(cube.shape[2]):
for y in range(cube.shape[1]):
if np.int((counter / number_of_pixels) * 100) != fac_prev:
......@@ -718,13 +723,16 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
for z in np.arange(zzmin, zzmax, zstep):
p[0] = z
pstart = np.copy(p)
fit0 = mpfit.mpfit(myelspectrum, xall=pstart, functkw=fa, quiet=True, parinfo=parinfo)
ti = time.time()
fit0 = mpfit.mpfit(myelspectrum, xall=pstart, functkw=fa, quiet=True, parinfo=parinfo, ftol=ftol, xtol=xtol, gtol=gtol, maxiter=maxiter, factor=factor)
tf = time.time()
tt += tf - ti
if fit0.status == 0: continue
if fit0.status == -16: continue
if fit0.fnorm <= minfnorm:
minfnorm = fit0.fnorm
fit = copy.copy(fit0)
#print(z,tf - ti)
if minfnorm == np.infty: continue
#print(fit.params[2 * els.index + 1 + 12], parinfo[2 * els.index + 1 + 12]['limits'], fit.params[2 * els.index + 1 + 13])
paramscube[:, y, x] = fit.params[:]
......@@ -734,7 +742,10 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
statusmap[y, x] = fit.status
modcube[indok, y, x] = elspectrum(fit.params, wave=wave, lines=linewave, psf=conf['SPSF'], degc=conf['DGCTNUM'])
tff = time.time()
#print(tff - tii)
#print(tt)
#Creating outputs
#Error normalisation by reduced chi2
perrorcube = perrorcube * np.sqrt(fnormmap / dofmap)
......@@ -872,7 +883,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
return rvmap
def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False, config=None):
def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False, config=None, ftol=1e-7, xtol=1e-7, gtol=1e-7, factor=10, maxiter=200):
"""CAMEL (Cube Analysis: Moment maps of Emission Lines) computes moment maps
Parameters
......@@ -962,7 +973,7 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False,
varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var.fits')
#Creating the maps
rv = buildmaps(conf, cube, hdr, var=var, plot=plot, debug=debug, free_ratio=free_ratio)
rv = buildmaps(conf, cube, hdr, var=var, plot=plot, debug=debug, free_ratio=free_ratio, ftol=ftol, xtol=xtol, gtol=gtol, factor=factor, maxiter=maxiter)
t1 = time.time()
#print(t1-t0)
......@@ -1018,6 +1029,11 @@ def main(argv):
parser.add_argument('--debug', '-d', action="store_true", dest="debug", default=False, help="keyword to show debug information (not implemented yet)")
parser.add_argument('--free_line_ratio', '-r', action="store_true", dest="free_ratio", default=False, help="keyword to allow unconstrained line ratio")
parser.add_argument('--multi_ext', '-m', action="store_true", dest="multi_ext", default=False, help="keyword to indicate if input fits has multiple extensions")
parser.add_argument('--ftol', action="store", dest="ftol", default=1e-7, type=float, help="ftol argument in mpfit")
parser.add_argument('--xtol', action="store", dest="xtol", default=1e-7, type=float, help="xtol argument in mpfit")
parser.add_argument('--gtol', action="store", dest="gtol", default=1e-7, type=float, help="gtol argument in mpfit")
parser.add_argument('--maxiter', action="store", dest="maxiter", default=200, type=int, help="maxiter argument in mpfit")
parser.add_argument('--factor', action="store", dest="factor", default=10, type=float, help="factor argument in mpfit (between 0.1 and 100)")
parser.add_argument('--FITSFILE', action="store", dest="FITSFILE", default=None, help="name of the input cube")
parser.add_argument('--OUTPUT', action="store", dest="OUTPUT", default=None, help="generic output name")
parser.add_argument('--SKYFILE', action="store", dest="SKYFILE", default=None, help="name of the variance cube (can be a sky cube)")
......@@ -1077,7 +1093,7 @@ def main(argv):
opts = vars(options) # this is a dictionnary
opts['OII3726'] = opts['OII3729']
camel(options.filename, plot=options.plot, debug=options.debug, free_ratio=options.free_ratio, multi_ext=options.multi_ext, config=opts)
camel(options.filename, plot=options.plot, debug=options.debug, free_ratio=options.free_ratio, multi_ext=options.multi_ext, ftol=options.ftol, xtol=options.xtol, gtol=options.gtol, maxiter=options.maxiter, factor=options.factor, config=opts)
return
if __name__ == "__main__":
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment