Commit 6b4c6e82 authored by bepinat's avatar bepinat

bugs corrections EXTRAL

parent 1964b3e4
...@@ -83,8 +83,9 @@ def readconfig(filename, config): ...@@ -83,8 +83,9 @@ def readconfig(filename, config):
conf = {} conf = {}
if filename != None: if filename != None:
data = open(filename) data = open(filename)
for line in data: count_extral = 0
keyvalcom = line.split('= ') for raw in data:
keyvalcom = raw.split('= ')
key = keyvalcom[0].rstrip() # on vire les espaces de fin de keyword key = keyvalcom[0].rstrip() # on vire les espaces de fin de keyword
value = keyvalcom[1].split('/ ')[0].rstrip('\t').rstrip() # on vire les tabultation et les espaces de fin de keyword value = keyvalcom[1].split('/ ')[0].rstrip('\t').rstrip() # on vire les tabultation et les espaces de fin de keyword
value = value.replace("'", '') value = value.replace("'", '')
...@@ -97,11 +98,28 @@ def readconfig(filename, config): ...@@ -97,11 +98,28 @@ def readconfig(filename, config):
value = False value = False
elif value.upper() == 'TRUE': elif value.upper() == 'TRUE':
value = True value = True
elif (value.count('.') == 1) & (value.replace('.', '').isdigit()): #elif (value.count('.') == 1) & (value.replace('.', '').isdigit()):
value = float(value) #value = float(value)
elif value.isdigit(): elif (value.count('.') == 1):
value = int(value) try:
conf[key] = value value = float(value)
except:
value = value
#elif value.isdigit():
#value = int(value)
elif (value.count('.') == 0):
try:
value = int(value)
except:
value = value
if (key == 'EXTRAL') & (value != None):
if count_extral == 0:
conf[key] = np.array([value])
else:
conf[key] = np.append(conf[key], value)
count_extral += 1
else:
conf[key] = value
# add input options to configuration dictionnary + check needed keywords # add input options to configuration dictionnary + check needed keywords
if config != None: if config != None:
...@@ -147,7 +165,7 @@ def readcubes(conf, multi_ext=False): ...@@ -147,7 +165,7 @@ def readcubes(conf, multi_ext=False):
var = hdul[2].data var = hdul[2].data
varhdr = hdul[2].header varhdr = hdul[2].header
elif (not conf['SKYFILE']): elif (not conf['SKYFILE']):
var = 0 var = None
varhdr = '' varhdr = ''
else: else:
hdul = pf.open(conf['SKYFILE']) hdul = pf.open(conf['SKYFILE'])
...@@ -158,7 +176,7 @@ def readcubes(conf, multi_ext=False): ...@@ -158,7 +176,7 @@ def readcubes(conf, multi_ext=False):
cube = hdul[0].data cube = hdul[0].data
hdr = hdul[0].header hdr = hdul[0].header
if (not conf['SKYFILE']): if (not conf['SKYFILE']):
var = 0 var = None
varhdr = '' varhdr = ''
else: else:
hdul = pf.open(conf['SKYFILE']) hdul = pf.open(conf['SKYFILE'])
...@@ -177,6 +195,8 @@ def testcubes(cube, var): ...@@ -177,6 +195,8 @@ def testcubes(cube, var):
input variance data input variance data
""" """
if var == None:
return
if var.ndim == 3: if var.ndim == 3:
if var.shape != cube.shape: if var.shape != cube.shape:
sys.exit(2) sys.exit(2)
...@@ -451,11 +471,11 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal ...@@ -451,11 +471,11 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#Initialisation of lines #Initialisation of lines
els = lines() els = lines()
for line in els.lines: for lline in els.lines:
#check the line name is set in the configuration file #check the line name is set in the configuration file
#if (not conf.keys().isdisjoint([line])): # c'est du python3 #if (not conf.keys().isdisjoint([line])): # c'est du python3
if conf.has_key(line): # python 2.7 if conf.has_key(lline): # python 2.7
els.lines[line].fit = conf[line] els.lines[lline].fit = conf[lline]
if conf['EXTRAL'] != None: if conf['EXTRAL'] != None:
for i in np.arange(np.shape(conf['EXTRAL'])[0]): for i in np.arange(np.shape(conf['EXTRAL'])[0]):
...@@ -469,7 +489,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal ...@@ -469,7 +489,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
zcosmo = conf['REDSHIFT'] zcosmo = conf['REDSHIFT']
zmax = conf['REDMAX'] zmax = conf['REDMAX']
zmin = conf['REDMIN'] zmin = conf['REDMIN']
zrange = zmax-zmin zrange = zmax - zmin
argsorted = np.argsort(np.array([els.lines[i].index for i in els.lines])) argsorted = np.argsort(np.array([els.lines[i].index for i in els.lines]))
lineindex = np.array([els.lines[i].index for i in els.lines])[argsorted] lineindex = np.array([els.lines[i].index for i in els.lines])[argsorted]
linefit = np.array([els.lines[i].fit for i in els.lines])[argsorted] linefit = np.array([els.lines[i].fit for i in els.lines])[argsorted]
...@@ -525,14 +545,19 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal ...@@ -525,14 +545,19 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#lastind[i] = condo[0][np.size(condo) - 1] #lastind[i] = condo[0][np.size(condo) - 1]
#Weights #Weights
if var == None: var = np.ones([indok.size, cube.shape[1], cube.shape[2]]) #if var == None: var = np.ones([indok.size, cube.shape[1], cube.shape[2]])
if var == None: var = np.ones(cube.shape)
if var.ndim == 1: var = np.tile(var.reshape(var.size, 1, 1), (1, cube.shape[1], cube.shape[2])) 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 #if var.ndim == 3: weight = ((1. / var[indok,:,:]) / np.sum(1. / var[indok,:,:], axis=0)) * indok.size
if var.ndim == 3: weight = ((1. / var) / np.sum(1. / var, axis=0)) * var.shape[0] if var.ndim == 3: weight = ((1. / var) / np.sum(1. / var, axis=0)) * var.shape[0]
#Redshift step #Redshift step
zstep = (conf['DFIT'] * 1e3 / ct.c + zcosmo) / (1 - conf['DFIT'] * 1e3 / ct.c) - zcosmo zstep = (conf['DFIT'] * 1e3 / ct.c + zcosmo) / (1 - conf['DFIT'] * 1e3 / ct.c) - zcosmo
if zstep >= (zmax - zmin): zstep = (zmax - zmin) * 0.99
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
#Parameters contraints and informations #Parameters contraints and informations
params = np.zeros(1 + els.index * 3 + conf['DGCTNUM'] + 1) params = np.zeros(1 + els.index * 3 + conf['DGCTNUM'] + 1)
...@@ -553,6 +578,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal ...@@ -553,6 +578,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
commwindex = np.min(lineindex[np.where(linefit)]) commwindex = np.min(lineindex[np.where(linefit)])
parinfo[1 + commwindex]['fixed'] = 0 parinfo[1 + commwindex]['fixed'] = 0
linew = np.zeros(els.index) linew = np.zeros(els.index)
for i in range(els.index): for i in range(els.index):
name = linename[i] name = linename[i]
j = els.lines[els.lines[name].ref].index j = els.lines[els.lines[name].ref].index
...@@ -689,8 +715,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal ...@@ -689,8 +715,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#Inititalisation of continuum (degree 0) #Inititalisation of continuum (degree 0)
p[1 + els.index * 3] = np.median(fa['spectrum']) p[1 + els.index * 3] = np.median(fa['spectrum'])
for z in np.arange(zmin, zmax, zstep): for z in np.arange(zzmin, zzmax, zstep):
#print(z)
p[0] = z p[0] = z
pstart = np.copy(p) pstart = np.copy(p)
fit0 = mpfit.mpfit(myelspectrum, xall=pstart, functkw=fa, quiet=True, parinfo=parinfo) fit0 = mpfit.mpfit(myelspectrum, xall=pstart, functkw=fa, quiet=True, parinfo=parinfo)
...@@ -878,7 +903,10 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False, ...@@ -878,7 +903,10 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False,
return return
#Reading data #Reading data
print("Reading data cube " + conf['FITSFILE'] + " and variance " + conf['SKYFILE']) try:
print("Reading data cube " + conf['FITSFILE'] + " and variance " + conf['SKYFILE'])
except:
print("Reading data cube " + conf['FITSFILE'] + ". No variance data.")
try: try:
cube, hdr, var, varhdr = readcubes(conf, multi_ext=multi_ext) cube, hdr, var, varhdr = readcubes(conf, multi_ext=multi_ext)
except: except:
...@@ -888,7 +916,7 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False, ...@@ -888,7 +916,7 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_ext=False,
try: try:
testcubes(cube, var) testcubes(cube, var)
except: except:
print("Variance shape: %i; Data cube shape: %i"%(var.shape, cube.shape)) print("Variance shape: %i; Data cube shape: %i"%(np.shape(var), cube.shape))
print("Variance size problem") print("Variance size problem")
return return
......
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