Commit 79450fa2 authored by bepinat's avatar bepinat
Browse files

Inclusion of cap_mpfit from Michele Cappellari which works both with python2...

Inclusion of cap_mpfit from Michele Cappellari which works both with python2 and 3. This new CAMEL version also works both with python2 and 3.
parent ec6af203
#!/usr/bin/env python
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
......@@ -24,7 +24,8 @@ import os, sys
import copy
import numpy as np
import pyfits as pf
from mpfit import mpfit # python 2.7...
from cap_mpfit import mpfit # python 2.7 & 3!
#from mpfit import mpfit # python 2.7...
import ipdb
#from matplotlib import pyplot as plt
from scipy import ndimage
......@@ -134,10 +135,15 @@ def readconfig(filename, config):
if config != None:
keys = conf.keys()
for key in config:
if sys.version_info > (3,):
condkey = keys.isdisjoint([key])
else:
condkey = (not conf.has_key(key))
if config[key] != None:
conf[key] = config[key]
#elif keys.isdisjoint([key]): # c'est du python3
elif (not conf.has_key(key)): # c'est du python 2.7
#elif (not conf.has_key(key)): # c'est du python 2.7
elif condkey:
conf[key] = None
# XXX or give an error message but in that case, must test that the keywords have correct values
......@@ -145,10 +151,21 @@ def readconfig(filename, config):
needed_keys = ['FITSFILE', 'OUTPUT', 'SKYFILE', 'HALPHA', 'NII6548', 'NII6583', 'SII6716', 'SII6731', 'OIII4959', 'OIII5007', 'HBETA', 'OII3729', 'OII3726', 'EXTRAL', 'COMMW', 'REDSHIFT', 'REDMIN', 'REDMAX', 'INITW', 'WMIN', 'WMAX', 'DFIT', 'DGCTNUM', 'MFIT', 'SCLIP', 'XYCLIP', 'NCLIP', 'SPSF', 'WSMOOTH', 'SSMOOTH', 'THRES', 'MEDIAN', 'FITSPSF', 'XMIN', 'YMIN', 'ZMIN', 'XMAX', 'YMAX', 'ZMAX']
keys = conf.keys()
for key in needed_keys:
if sys.version_info > (3,):
condkey = keys.isdisjoint([key])
else:
condkey = (not conf.has_key(key))
#if keys.isdisjoint([key]): # c'est du python3
if (not conf.has_key(key)): # c'est du python 2.7
#if (not conf.has_key(key)): # c'est du python 2.7
if condkey:
conf[key] = None
if conf.has_key('OII'):
if sys.version_info > (3,):
condoii = (not keys.isdisjoint(['OII']))
else:
condoii = conf.has_key('OII')
#if conf.has_key('OII'):
if condoii:
conf['OII3729'] = conf['OII3726'] = conf.pop('OII')
if conf['SPSF'] == None: conf['SPSF'] = 0.
if conf['COMMW'] == None: conf['COMMW'] = False
......@@ -166,6 +183,7 @@ def readcubes(conf, multi_ext=False):
keyword to indicate if input fits has multiple extensions
"""
if multi_ext:
hdul = pf.open(conf['FITSFILE'])
cube = hdul[1].data
......@@ -186,7 +204,7 @@ def readcubes(conf, multi_ext=False):
cube = hdul[0].data
hdr = hdul[0].header
if conf['SKYFILE'] == conf['FITSFILE']:
var = np.zeros(cube.shape)
var = np.zeros(cube.shape, dtype='>f4')
if conf['XMIN'] == None:
conf['XMIN'] = 0
if conf['YMIN'] == None:
......@@ -265,15 +283,15 @@ def muse_whitelight_image(cube, hdr, filename):
data = np.sum(cube, axis=0)
hdrw = copy.deepcopy(hdr)
if 'CTYPE3' in hdr.ascard.keys(): del hdrw['CTYPE3']
if 'CRVAL3' in hdr.ascard.keys(): del hdrw['CRVAL3']
if 'CRPIX3' in hdr.ascard.keys(): del hdrw['CRPIX3']
if 'CUNIT3' in hdr.ascard.keys(): del hdrw['CUNIT3']
if 'CD3_3' in hdr.ascard.keys(): del hdrw['CD3_3']
if 'CD3_2' in hdr.ascard.keys(): del hdrw['CD3_2']
if 'CD3_1' in hdr.ascard.keys(): del hdrw['CD3_1']
if 'CD2_3' in hdr.ascard.keys(): del hdrw['CD2_3']
if 'CD1_3' in hdr.ascard.keys(): del hdrw['CD1_3']
if 'CTYPE3' in hdr.keys(): del hdrw['CTYPE3']
if 'CRVAL3' in hdr.keys(): del hdrw['CRVAL3']
if 'CRPIX3' in hdr.keys(): del hdrw['CRPIX3']
if 'CUNIT3' in hdr.keys(): del hdrw['CUNIT3']
if 'CD3_3' in hdr.keys(): del hdrw['CD3_3']
if 'CD3_2' in hdr.keys(): del hdrw['CD3_2']
if 'CD3_1' in hdr.keys(): del hdrw['CD3_1']
if 'CD2_3' in hdr.keys(): del hdrw['CD2_3']
if 'CD1_3' in hdr.keys(): del hdrw['CD1_3']
wlh = writedata(data, hdrw, filename)
return wlh
......@@ -361,7 +379,7 @@ def clipcube(cube, hdr, xy=3, sclip=3, npass=10):
"""
cube1 = np.zeros(cube.shape) # initialization of output cube
cube1 = np.zeros(cube.shape, dtype='>f4') # initialization of output cube
for z in range(cube.shape[0]): # loop on the z range
im = cube[z,:]
count = 0
......@@ -393,7 +411,7 @@ def spatialsmooth(cube, hdr, fwhm):
full width half maximum of the 2D Gaussian kernel
"""
cube1 = np.zeros(cube.shape) # initialization of output cube
cube1 = np.zeros(cube.shape, dtype='>f4') # initialization of output cube
sigma = fwhm / (2. * np.sqrt(2. * np.log(2.)))
for z in range(cube.shape[0]): # loop on the z range
cube1[z,:] = ndimage.gaussian_filter(cube[z,:], sigma)
......@@ -413,7 +431,7 @@ def spectralsmooth(cube, hdr, fwhm):
full width half maximum of the 1D Gaussian kernel
"""
cube1 = np.zeros(cube.shape) # initialization of output cube
cube1 = np.zeros(cube.shape, dtype='>f4') # initialization of output cube
sigma = fwhm / (2. * np.sqrt(2. * np.log(2.)))
for x in range(cube.shape[2]): # loop on the x range
for y in range(cube.shape[1]): # loop on the y range
......@@ -494,9 +512,9 @@ def waveindgen(hdr):
if cunit.strip().lower() in {'microns', 'micron', 'micrometers', 'mum'}: lconvfac=1e4
if cunit.strip().lower() in {'angstroms', 'angstrom'}: lconvfac=1
if cunit.strip().lower() in {'nanometers', 'nm'}: lconvfac=1e3
if 'CDELT3' in hdr.ascard.keys():
if 'CDELT3' in hdr.keys():
sclp = hdr['CDELT3']
elif 'CD3_3' in hdr.ascard.keys():
elif 'CD3_3' in hdr.keys():
sclp = hdr['CD3_3']
wave = lconvfac * (sclp * (np.arange(hdr['NAXIS3']) - hdr['CRPIX3'] + 1) + hdr['CRVAL3'])
return wave, lconvfac, sclp
......@@ -530,8 +548,13 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
els = lines()
for lline in els.lines:
#check the line name is set in the configuration file
if sys.version_info > (3,):
condline = (not conf.keys().isdisjoint([lline]))
else:
condline = conf.has_key(lline)
#if (not conf.keys().isdisjoint([line])): # c'est du python3
if conf.has_key(lline): # python 2.7
#if conf.has_key(lline): # python 2.7
if condline:
els.lines[lline].fit = conf[lline]
if conf['EXTRAL'] != None:
......@@ -554,10 +577,10 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
linewave = np.array([els.lines[i].wave for i in els.lines])[argsorted]
linename = np.array([els.lines[i].name for i in els.lines])[argsorted]
indok = np.zeros(0, dtype = int)
firstind = np.zeros(linefit.size)
lastind = np.zeros(linefit.size)
firstind0 = np.zeros(linefit.size)
lastind0 = np.zeros(linefit.size)
firstind = np.zeros(linefit.size, dtype='>f4')
lastind = np.zeros(linefit.size, dtype='>f4')
firstind0 = np.zeros(linefit.size, dtype='>f4')
lastind0 = np.zeros(linefit.size, dtype='>f4')
for i in range(linefit.size):
if (not linefit[i]): continue
linf = np.max([(zmin - zrange / 2. + 1) * linewave[i] , np.min(wave0)])
......@@ -604,7 +627,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#Weights
#if var == None: var = np.ones([indok.size, cube.shape[1], cube.shape[2]])
if var == None: var = np.ones(cube.shape)
if var == None: var = np.ones(cube.shape, dtype='>f4')
if var.ndim == 1: var = np.tile(var.reshape(var.size, 1, 1), (1, cube.shape[1], cube.shape[2]))
#if var.ndim == 3: weight = ((1. / var[indok,:,:]) / np.sum(1. / var[indok,:,:], axis=0)) * indok.size
# weight inverse variance
......@@ -622,7 +645,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
zzmin = zmean + np.ceil((zmin - zmean) / zstep) * zstep
#Parameters contraints and informations
params = np.zeros(1 + els.index * 3 + conf['DGCTNUM'] + 1)
params = np.zeros(1 + els.index * 3 + conf['DGCTNUM'] + 1, dtype='>f4')
parbase = {'value':0., 'fixed':1, 'limited':[0,0], 'limits':[0.,0.], 'tied':'', 'mpmaxstep':0, 'mpminstep':0, 'parname':''}
parinfo = []
......@@ -640,7 +663,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
if conf['COMMW']:
commwindex = np.min(lineindex[np.where(linefit)])
parinfo[1 + commwindex]['fixed'] = 0
linew = np.zeros(els.index)
linew = np.zeros(els.index, dtype='>f4')
for i in range(els.index):
name = linename[i]
......@@ -736,12 +759,12 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
parinfo[1 + els.index * 3 + i]['fixed'] = 0
#Output initialisation to be filled
paramscube = np.zeros(np.append(len(params), cube.shape[1:]))
perrorcube = np.zeros(np.append(len(params), cube.shape[1:]))
statusmap = np.zeros(cube.shape[1:])
dofmap = np.zeros(cube.shape[1:])
fnormmap = np.zeros(cube.shape[1:])
modcube = np.zeros(cube.shape)
paramscube = np.zeros(np.append(len(params), cube.shape[1:]), dtype='>f4')
perrorcube = np.zeros(np.append(len(params), cube.shape[1:]), dtype='>f4')
statusmap = np.zeros(cube.shape[1:], dtype='>f4')
dofmap = np.zeros(cube.shape[1:], dtype='>f4')
fnormmap = np.zeros(cube.shape[1:], dtype='>f4')
modcube = np.zeros(cube.shape, dtype='>f4')
#List of extra parameters
#xxx Ajuster wave et spectrum pour ne garder que ce qui nous intéresse?
......@@ -786,7 +809,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
p[0] = z
pstart = np.copy(p)
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)
fit0 = mpfit(myelspectrum, xall=pstart, functkw=fa, quiet=True, parinfo=parinfo, ftol=ftol, xtol=xtol, gtol=gtol, maxiter=maxiter, factor=factor)
#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
......@@ -814,8 +838,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#Continuum
#if we want the coefficient, we must multiply them by factor + give lambda[0]
contcube = np.zeros(cube.shape)
econtcube = np.zeros(cube.shape)
contcube = np.zeros(cube.shape, dtype='>f4')
econtcube = np.zeros(cube.shape, dtype='>f4')
for i in range(conf['DGCTNUM'] + 1):
contcube += paramscube[1 + els.index * 3 + i,:,:] * (wave0 - wave[0]).reshape(wave0.size, 1, 1) ** i
econtcube += perrorcube[1 + els.index * 3 + i,:,:] * (wave0 - wave[0]).reshape(wave0.size, 1, 1) ** i
......@@ -826,7 +850,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#Residual cube
#residual
residualcube = np.zeros(cube.shape)
residualcube = np.zeros(cube.shape, dtype='>f4')
residualcube[indok,:,:] = cube[indok,:,:] - modcube[indok,:,:]
#dispcont
dispcontmap = np.std(cube[indok,:,:] - modcube[indok,:,:], axis=0)
......@@ -841,17 +865,17 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
ervmap = compute_fwhm(ezmap, z=zmap, zcosmo=zcosmo)
#Lines
wavemaps = np.zeros(np.append(els.index, cube.shape[1:]))
ewavemaps = np.zeros(np.append(els.index, cube.shape[1:]))
wavedispmaps = np.zeros(np.append(els.index, cube.shape[1:]))
ewavedispmaps = np.zeros(np.append(els.index, cube.shape[1:]))
dispmaps = np.zeros(np.append(els.index, cube.shape[1:]))
edispmaps = np.zeros(np.append(els.index, cube.shape[1:]))
intmaps = np.zeros(np.append(els.index, cube.shape[1:]))
eintmaps = np.zeros(np.append(els.index, cube.shape[1:]))
fluxmaps = np.zeros(np.append(els.index, cube.shape[1:]))
efluxmaps = np.zeros(np.append(els.index, cube.shape[1:]))
snrmaps = np.zeros(np.append(els.index, cube.shape[1:]))
wavemaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
ewavemaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
wavedispmaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
ewavedispmaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
dispmaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
edispmaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
intmaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
eintmaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
fluxmaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
efluxmaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
snrmaps = np.zeros(np.append(els.index, cube.shape[1:]), dtype='>f4')
# dzmap and edzmap are intermediate products
if conf['COMMW']:
......@@ -913,16 +937,16 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
hh = writedata(econtcube, hdr, conf['OUTPUT'] + '_econtcube.fits')
#Images
if 'CTYPE3' in hdr.ascard.keys(): del hdr['CTYPE3']
if 'CRVAL3' in hdr.ascard.keys(): del hdr['CRVAL3']
if 'CRPIX3' in hdr.ascard.keys(): del hdr['CRPIX3']
if 'CUNIT3' in hdr.ascard.keys(): del hdr['CUNIT3']
if 'CD3_3' in hdr.ascard.keys(): del hdr['CD3_3']
if 'CD3_2' in hdr.ascard.keys(): del hdr['CD3_2']
if 'CD3_1' in hdr.ascard.keys(): del hdr['CD3_1']
if 'CD2_3' in hdr.ascard.keys(): del hdr['CD2_3']
if 'CD1_3' in hdr.ascard.keys(): del hdr['CD1_3']
if 'CDELT3' in hdr.ascard.keys(): del hdr['CDELT3']
if 'CTYPE3' in hdr.keys(): del hdr['CTYPE3']
if 'CRVAL3' in hdr.keys(): del hdr['CRVAL3']
if 'CRPIX3' in hdr.keys(): del hdr['CRPIX3']
if 'CUNIT3' in hdr.keys(): del hdr['CUNIT3']
if 'CD3_3' in hdr.keys(): del hdr['CD3_3']
if 'CD3_2' in hdr.keys(): del hdr['CD3_2']
if 'CD3_1' in hdr.keys(): del hdr['CD3_1']
if 'CD2_3' in hdr.keys(): del hdr['CD2_3']
if 'CD1_3' in hdr.keys(): del hdr['CD1_3']
if 'CDELT3' in hdr.keys(): del hdr['CDELT3']
suffixe = ''
if conf['COMMW']:
......@@ -1012,8 +1036,8 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False,
try:
testcubes(cube, var)
except:
logger.info("Variance shape: %i; Data cube shape: %i"%(np.shape(var), cube.shape))
logger.info("Variance size problem")
#logger.info("Variance shape: %i; Data cube shape: %i"%(np.shape(var), cube.shape))
return
......
This diff is collapsed.
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