Commit 0214c6ae authored by bepinat's avatar bepinat

Use of a normalisation of data before fit to avoid failed fits. Outputs are...

Use of a normalisation of data before fit to avoid failed fits. Outputs are corrected so that they correspond to the initial input cubes
parent a74e92c9
......@@ -915,6 +915,9 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
number_of_pixels = np.float(cube.shape[2] * cube.shape[1])
tii = time.time()
tt = 0
normsp = 10. ** np.floor(np.log10(np.nanstd(cube[indok, :, :])))
for x in range(cube.shape[2]):
for y in range(cube.shape[1]):
if np.int((counter / number_of_pixels) * 100) != fac_prev:
......@@ -924,8 +927,9 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#if (counter / 10.) == np.int(counter / 10.): logger.debug('pixel %i / %i'%(counter, number_of_pixels))
counter += 1
fa['spectrum'] = cube[indok, y, x] # Spectrum
fa['err'] = np.sqrt(var[indok, y, x]) # Error
fa['spectrum'] = cube[indok, y, x] / normsp # Spectrum
fa['err'] = np.sqrt(var[indok, y, x]) / normsp # Error
#fa['weight'] = np.sqrt(weight[indok, y, x]) # Weight
if (np.min(fa['spectrum']) == np.max(fa['spectrum'])) & (np.min(fa['spectrum']) == 0):
continue # No need for fitting
......@@ -999,6 +1003,7 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#print(tt)
#Creating outputs
modcube *= normsp
#Error normalisation by reduced chi2
perrorcube = perrorcube * np.sqrt(fnormmap / dofmap)
......@@ -1013,6 +1018,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#; continuum is computed at the longest wavelength, not necessarily optimal
#cont[m,n]+=pfinal[2*n_elements(l)+1+i]*double(lambda[n_elements(lambda)-1]-lambda[0])^i/double(lambda[n_elements(lambda)-1]-lambda[0])
#econt[m,n]+=(paramerr[2*n_elements(l)+1+i]*double(lambda[n_elements(lambda)-1]-lambda[0])^i/double(lambda[n_elements(lambda)-1]-lambda[0]))^2
contcube *= normsp
econtcube *= normsp
#Residual cube
#residual
......@@ -1064,15 +1071,15 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
edispmaps[comp * els.index + i, :, :] = compute_fwhm(edzmap, z=zmap[comp * els.index + i, :, :], zcosmo=zcosmo)
#Flux
intmaps[comp * els.index + i, :, :] = paramscube[intind[comp * els.index + i], :, :] * paramscube[ratioind[comp * els.index + i], :, :]
intmaps[comp * els.index + i, :, :] = paramscube[intind[comp * els.index + i], :, :] * paramscube[ratioind[comp * els.index + i], :, :] * normsp
fluxmaps[comp * els.index + i, :, :] = np.sum(intmaps[comp * els.index + i, :, :] * np.exp(-0.5 * (wave.reshape(wave.size, 1, 1) - linewave[i] * (1 + zmap[comp * els.index + i, :, :])) ** 2 / ((dzmap * linewave[i]) ** 2 + conf['SPSF'] ** 2)) * sclp * lconvfac, axis=0)
snrmaps[comp * els.index + i, :, :] = np.sum(weight[firstind0[i]:lastind0[i], :, :] * (modcube[firstind0[i]:lastind0[i], :, :] - contcube[firstind0[i]:lastind0[i], :, :]), axis=0) * sclp * lconvfac / (np.sqrt(np.pi * 2) * np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2) * np.sqrt(np.sum(weight[firstind0[i]:lastind0[i], :, :] * (cube[firstind0[i]:lastind0[i], :, :] - modcube[firstind0[i]:lastind0[i], :, :]) ** 2, axis=0) / np.sum(weight[firstind0[i]:lastind0[i], :, :], axis=0))) * (lastind0[i] - firstind0[i]) / np.sum(weight[firstind0[i]:lastind0[i], :, :], axis=0)
if free_ratio:
eintmaps[comp * els.index + i, :, :] = perrorcube[intind[comp * els.index + i], :, :]
eintmaps[comp * els.index + i, :, :] = perrorcube[intind[comp * els.index + i], :, :] * normsp
elif i != j:
eintmaps[comp * els.index + i, :, :] = np.sqrt((paramscube[intind[comp * els.index + i], :, :] * perrorcube[ratioind[comp * els.index + i], :, :]) ** 2 + (perrorcube[intind[comp * els.index + j], :, :] * paramscube[ratioind[comp * els.index + i], :, :]) ** 2)
eintmaps[comp * els.index + i, :, :] = np.sqrt((paramscube[intind[comp * els.index + i], :, :] * perrorcube[ratioind[comp * els.index + i], :, :]) ** 2 + (perrorcube[intind[comp * els.index + j], :, :] * paramscube[ratioind[comp * els.index + i], :, :]) ** 2) * normsp
else:
eintmaps[comp * els.index + i, :, :] = perrorcube[intind[comp * els.index + i], :, :]
eintmaps[comp * els.index + i, :, :] = perrorcube[intind[comp * els.index + i], :, :] * normsp
efluxmaps[comp * els.index + i, :, :] = np.sqrt(np.pi * 2) * np.sqrt((eintmaps[comp * els.index + i, :, :] * np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2)) ** 2 + (linewave[i] * dzmap / np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2) * intmaps[comp * els.index + i, :, :] * edzmap) ** 2)
......@@ -1303,7 +1310,7 @@ def camel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul
if var.ndim == 3:
var, varhdr = spectralsmooth(var, varhdr, conf['WSMOOTH'])
# variance is divided by the integral of the kernel to account for the variance lowering compared to the signal when smoothed
var /= ((conf['SSMOOTH'] / (2. * np.sqrt(2. * np.log(2.)))) * np.sqrt(2 * np.pi))
var /= ((conf['WSMOOTH'] / (2. * np.sqrt(2. * np.log(2.)))) * np.sqrt(2 * np.pi))
varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var.fits')
if (conf['SSMOOTH'] != 0) & (conf['SSMOOTH'] is not None):
......
......@@ -718,8 +718,10 @@ def buildmaps(conf, data, hdr, var=None, plot=False, debug=False, free_ratio=Fal
tii = time.time()
tt = 0
fa['spectrum'] = data[indok] # Spectrum
fa['err'] = np.sqrt(var[indok]) # Error
normsp = 10. ** np.floor(np.log10(np.nanstd(cube[indok, :, :])))
fa['spectrum'] = data[indok] / normsp # Spectrum
fa['err'] = np.sqrt(var[indok]) / normsp # Error
if (np.min(fa['spectrum']) == np.max(fa['spectrum'])) & (np.min(fa['spectrum']) == 0):
return # No need for fitting
......@@ -787,6 +789,7 @@ def buildmaps(conf, data, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#print(tt)
#Creating outputs
moddata *= normsp
#Error normalisation by reduced chi2
perrordata = perrordata * np.sqrt(fnorm / dof)
......@@ -797,6 +800,8 @@ def buildmaps(conf, data, hdr, var=None, plot=False, debug=False, free_ratio=Fal
for i in range(conf['DGCTNUM'] + 1):
contdata += paramsdata[contind[i]] * (wave0 - wave[0]) ** i
econtdata += perrordata[contind[i]] * (wave0 - wave[0]) ** i
contdata *= normsp
econtdata *= normsp
#Residual data
#residual
......@@ -848,15 +853,15 @@ def buildmaps(conf, data, hdr, var=None, plot=False, debug=False, free_ratio=Fal
edispmaps[comp * els.index + i] = compute_fwhm(edzmap, z=zmap[comp * els.index + i], zcosmo=zcosmo)
#Flux
intmaps[comp * els.index + i] = paramsdata[intind[comp * els.index + i]] * paramsdata[ratioind[comp * els.index + i]]
intmaps[comp * els.index + i] = paramsdata[intind[comp * els.index + i]] * paramsdata[ratioind[comp * els.index + i]] * normsp
fluxmaps[comp * els.index + i] = np.sum(intmaps[comp * els.index + i] * np.exp(-0.5 * (wave.reshape(wave.size, 1, 1) - linewave[i] * (1 + zmap[comp * els.index + i])) ** 2 / ((dzmap * linewave[i]) ** 2 + conf['SPSF'] ** 2)) * sclp * lconvfac, axis=0)
snrmaps[comp * els.index + i] = np.sum(weight[firstind0[i]:lastind0[i]] * (moddata[firstind0[i]:lastind0[i]] - contdata[firstind0[i]:lastind0[i]]), axis=0) * sclp * lconvfac / (np.sqrt(np.pi * 2) * np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2) * np.sqrt(np.sum(weight[firstind0[i]:lastind0[i]] * (data[firstind0[i]:lastind0[i]] - moddata[firstind0[i]:lastind0[i]]) ** 2, axis=0) / np.sum(weight[firstind0[i]:lastind0[i]], axis=0))) * (lastind0[i] - firstind0[i]) / np.sum(weight[firstind0[i]:lastind0[i]], axis=0)
if free_ratio:
eintmaps[comp * els.index + i] = perrordata[intind[comp * els.index + i]]
eintmaps[comp * els.index + i] = perrordata[intind[comp * els.index + i]]normsp
elif i != j:
eintmaps[comp * els.index + i] = np.sqrt((paramsdata[intind[comp * els.index + i]] * perrordata[ratioind[comp * els.index + i]]) ** 2 + (perrordata[intind[comp * els.index + j]] * paramsdata[ratioind[comp * els.index + i]]) ** 2)
eintmaps[comp * els.index + i] = np.sqrt((paramsdata[intind[comp * els.index + i]] * perrordata[ratioind[comp * els.index + i]]) ** 2 + (perrordata[intind[comp * els.index + j]] * paramsdata[ratioind[comp * els.index + i]]) ** 2)normsp
else:
eintmaps[comp * els.index + i] = perrordata[intind[comp * els.index + i]]
eintmaps[comp * els.index + i] = perrordata[intind[comp * els.index + i]]normsp
efluxmaps[comp * els.index + i] = np.sqrt(np.pi * 2) * np.sqrt((eintmaps[comp * els.index + i] * np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2)) ** 2 + (linewave[i] * dzmap / np.sqrt((linewave[i] * dzmap) ** 2 + conf['SPSF'] ** 2) * intmaps[comp * els.index + i] * edzmap) ** 2)
......@@ -1047,7 +1052,7 @@ def samel(filename, plot=False, debug=False, free_ratio=False, multi_comp=1, mul
if var is not None:
var, varhdr = spectralsmooth(var, varhdr, conf['WSMOOTH'])
# variance is divided by the integral of the kernel to account for the variance lowering compared to the signal when smoothed
var /= ((conf['SSMOOTH'] / (2. * np.sqrt(2. * np.log(2.)))) * np.sqrt(2 * np.pi))
var /= ((conf['WSMOOTH'] / (2. * np.sqrt(2. * np.log(2.)))) * np.sqrt(2 * np.pi))
varhdr = writedata(var, varhdr, conf['OUTPUT'] + '_var.fits')
if smooth:
......
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