#!/usr/bin/env python # -*- coding: utf-8 -*- """ calculate nonegative projections of the spectral basis to the measured spectra. Projection basis was found as nonnegative sparse decomposition of about 3500 spectra to 11 components. In the first part of the algoritm is remove stray light and made a corrections of the instrumental function It is very important because highly subpixel precision is necassary. First correction is made as a blind deconvolution of shape and then a zero and first order of the wavelength polynome are made. Than a projections even with realible estimate of the statistical errors are made. In the second part of the script is Integrated data analyzis of the radiation. Data from spectrometer with low time but very high wavelength resolution van be combinated with higt time, but low wavelength resoltion from othe diagnostics as photodiodes, cameras etc.. Advantages: extremly high sensitivity, can be identified even weak lines lost between strong surrounding lines of other ions. Immune agains line overburning. (Relatively) fast caculaton. It can never missidentified one of the strong lines. Disadvantage: In rare cases there can happend something interesting which can not be described by these limited number of components, one kind of event is here:7307,7960;9298,9156;10097 and also probably a new iont from unknow element can occure 9664,9665,9670,9839,9906,10131,10136. And finally a molecular spectrum 9504 On the other hand, without this algoritm it would by imposible to find that there is something interesting in this discharges. Autor: Tomas Odstrcil date:12.12.2012 """ from numpy import * import matplotlib matplotlib.rcParams['backend'] = 'Agg' from matplotlib.pyplot import * from SpectrometerControl import * from scipy.linalg import norm, inv,qr,pinv from scipy.special import erfinv, erf from scipy.stats.mstats import mquantiles from scikits.sparse.cholmod import cholesky, analyze, cholesky_AAt,CholmodError from scipy.signal import fftconvolve,convolve,gaussian,order_filter from numpy import linalg import time from scipy import sparse from scipy.sparse import * from scipy.optimize import minimize_scalar,nnls,fmin_powell,fmin_cg,fmin_bfgs,fmin_ncg from scipy.fftpack import fft,ifft from fftshift import fftshift from pygolem_lite.modules import save_adv,load_adv from scipy.sparse import block_diag from pygolem_lite import saveconst from scipy.interpolate import interp1d names = ('HI','OI','OII', 'OIII', 'HeI', 'CII', 'CIII', 'NI','NII', 'NIII','Mystery1', 'CIV+NIV+OIV' ) styles = ('b-','r-','r--', 'r-.' , 'g-' , 'k--', 'k-.' , 'y-','y--', 'y-.' , 'c-' , 'm:' ) def Shift(x,s,win,lam): #místo přímého posuvu udělám inverzi posuvu zpět, je to stabilnější s-= 0.5 # systematický posuv, aby se to shodovalo s lineární interpolací s*= -1 #invezní posuv win = (win/2)*2 diag_data = ones((2,win)) diag_data[1,:] *=s-floor(s) diag_data[0,:] *=floor(s)+1-s S = spdiags(diag_data, (-floor(s),-floor(s)-1), win, win,format='csr') f0 = zeros(win) f0[win/2] = 1 factor = cholesky_AAt(S.T,beta = lam) g = squeeze(factor(S.T*f0)) Sx = fftconvolve(x,g,mode = 'same') return Sx def GapsFilling(signal,win = 100, lam = 1e-1): #BUG loadovat to z pygolema """ =============================== Gaps Filling Filter 0.1 ===================== reconstruct the corrupted data (data with nans) by tikhonov-philips regularization with regulariting by laplace operator. And return smoothed retrofit Reconstruction is based on the invertation of the identical operator with zeros at the lines corresponding to the mising signal. due to memory and speed limitation the reconstruction is done on the overalaping intervals with width "win" signal - long data vector win - width of the recosntruction interval - it mas be much bigger than the gaps width lam - regularization parameter, dependes on the noise in data Autor: Tomas Odstrcil 2012 """ from scikits.sparse.cholmod import cholesky, analyze,cholesky_AAt from scipy.sparse import spdiags, eye n = len(signal) # signal extension -- avoid boundary effects n_ext = (n/win+3)*win ext_signal = ones(n_ext)*nan side = (n_ext-n)/2 ext_signal[side:-side-n%2] = signal ext_signal[:side+1] = median(signal[~isnan(signal)][:win/2]) ext_signal[-side-1:] = median(signal[~isnan(signal)][-win/2:]) intervals = arange(0,n_ext, win/2) ind_nan = isnan(ext_signal) ext_signal[ind_nan] = 0 recon = copy(ext_signal) diag_data = ones((2,win)) diag_data[1,:]*=-1 D = spdiags(diag_data, (0,1), win, win,format='csr') DD = D.T*D I = eye(win,win,format='csc') Factor = cholesky_AAt(DD, 1./lam) for i in range(len(intervals)-2): gaps = spdiags( int_(ind_nan[intervals[i]:intervals[i+2]]),0, win, win,format='csc') # use overlapping intervals !!! Factor.update_inplace(gaps/sqrt(lam), subtract=True) # speed up trick g = Factor(ext_signal[intervals[i]:intervals[i+2]]/lam) Factor.update_inplace(gaps/sqrt(lam), subtract=False) # speed up trick recon[(intervals[i]+intervals[i+1])/2:(intervals[i+1]+intervals[i+2])/2] = g[len(g)/4:-len(g)/4,0] #plot(recon,'r') #plot(ext_signal,'b') #savefig('gapsfill.png') #clf() chi2 = sum((ext_signal-recon)[~ind_nan]**2)/n recon = recon[(n_ext-n)/2:-(n_ext-n)/2] return recon, chi2 def Resample(tvec,vec,n,t_min = None,t_max = None,left = 0,right = 0): if t_min == None: t_min = amin(tvec) if t_max == None: t_max = amax(tvec) std = (t_max-t_min)/(tvec[1]-tvec[0])/n/2 std = max(std,1) gauss_win = gaussian(std*3,std) gauss_win/= sum(gauss_win) vec_smooth = fftconvolve(vec,gauss_win, mode = 'same') tvec_new = linspace(t_min,t_max,n) vec_new = interp(tvec_new,tvec, vec_smooth, left=left, right=right) return tvec_new,vec_new def downsample(data, sampleInteg = False, downsample = 1): downsampled = reshape(data, (-1,downsample), order = 'C') if sampleInteg: downsampled = sum(downsampled, axis = 1)/downsample else: downsampled = downsampled[:,0] return downsampled def upsample_smooth(data, win,lam, sampleInteg = False, upsample = 1): win = 2*(win/2) n = len(data) npix = upsample*win n_ext = (n/win+1)*win ext_signal = empty(n_ext) ext_signal[(n_ext-n)/2:-(n_ext-n)/2] = data ext_signal[:(n_ext-n)/2+1] = 0#median(data[:win/2]) ext_signal[-(n_ext-n)/2-1:] = 0#median(data[win/2:]) retrofit = zeros(n_ext*upsample) #solve the upsamling problem only for short interval win diag_data = ones((2,npix)) diag_data[1,:]*=-1 D = spdiags(diag_data, (0,1), npix, npix,format='csr') D = D.T*D DD = D.T*D if sampleInteg: ReductMatrix = kron(eye(win,win),ones((1,upsample))/upsample,format='csr') else: ReductMatrix = kron(eye(win,win),eye(1,upsample),format='csr') T = ReductMatrix f = ext_signal f0 =zeros(win) f0[win/2] = 1 factor = cholesky(T.T*T+lam*DD) g = squeeze(factor(T.T*f0)) #apply the single point solution g on the whole data vector T_Tf = zeros((n_ext,upsample)) T_Tf[:,0] = f T_Tf = squeeze(reshape(T_Tf, (-1,1), order = 'C')) if win*upsample > 100: #choose the faster way upsampled = fftconvolve(T_Tf,g,mode = 'same') else: upsampled = convolve(T_Tf,g,mode = 'same') retrofit[:-1] = upsampled[1:] #I don't know why but this must by done. retrofit = reshape(retrofit, (-1,upsample), order = 'C') if sampleInteg: retrofit = mean(retrofit, axis = 1) else: retrofit = retrofit[:,0] upsampled = upsampled[(n_ext-n)/2*upsample:-(n_ext-n)/2*upsample] retrofit = retrofit[(n_ext-n)/2:-(n_ext-n)/2] #print 'error ',sum(abs(data-retrofit)) return upsampled,retrofit def removeStrayLight(data): #experimentálně ověřeno, lepší algoritmus mě nenapadl m = len(data) win = 101 quantile = 0.1 shift = 1.4*(10*erfinv((quantile-0.5)*2)) domain = ones(win) rank = int(win*quantile) background = order_filter(data-shift,domain,rank) gauss = exp(-linspace(-1,1,win)**2) gauss-= min(gauss) gauss/= sum(gauss) background = fftconvolve(background,gauss , mode = 'same') data-= background*0.8 return data import pyfftw def deconv_opt(spectra,retrofit,overburned,lam,ker0=None): n = len(spectra) normal = fftconvolve(overburned, ones(5),mode='same') < 1./5 fr = fft(retrofit[normal], n)/sqrt(n) fs = fft( spectra[normal], n)/sqrt(n) ker_full = pyfftw.n_byte_align_empty(n, 16,dtype=complex) Fker_full = pyfftw.n_byte_align_empty(n, 16,dtype=complex) fft_forward = pyfftw.FFTW(ker_full,Fker_full, direction='FFTW_FORWARD', flags=['FFTW_ESTIMATE','FFTW_DESTROY_INPUT']) doF = sum(normal)*5. def f(ker): ker_full[:]=0 ker_full[n/2-(len(ker))/2+1:n/2+(len(ker))/2+1] = ker ker_full[:] = fftshift(ker_full) fft_forward() chi2 = norm(fs-fr*Fker_full)**2/doF return chi2+norm(ker)*lam if ker0 == None: L = 10 ker0 = zeros(L) ker0[L/2-1] = 1 else: L = len(ker0) #t = time.time() #y = fmin_powell(f, ker0, xtol=1e-5,ftol=1e-5,maxfun=1e6,disp=True) y= fmin_bfgs(f, ker0, gtol=1e-1,maxiter=1e6,disp=False) #print time.time()-t #exit() ker_full = zeros(n) ker_full[n/2-L/2:n/2+(L+1)/2] = real(y) return real(y),ker_full def optim_fun(noise,spectra,retrofit,overburned,lim): dec = deconv(noise,spectra,retrofit,overburned,lim) chi2 = norm(spectra-fftconvolve(retrofit,dec,mode='same'))**2/size(spectra)/5 return exp(abs(log(chi2))) def CalcProjections(spectrometr,shot, plasma_start,plasma_end): shot_number = 100000 #BUG!!!! tvec = shot.time_stamps/1000. wavelength, intensities = shot.getData() MaxIntensity = 32000 # hodnota po korekci nelinearity detektoru! if shot_number < 8500: MaxIntensity = 18200 overburned = intensities > MaxIntensity if shot_number < 7625: #koukali jsme se skrze sklo overburned[wavelength < 330,:] = True if shot_number > 10300: #koukali jsme se skrze sklo overburned[wavelength < 330,:] = True overburned[spectrometr.black_pixels,:] = False if spectrometr.hotPixels != []: overburned[spectrometr.hotPixels,:] = True data = load('SpectraComponents2.npz') #the most critical part, weeks of my work components = data['components'] intensities-= mean((intensities[spectrometr.black_pixels,:])) n_features = size(components,1) n_components = size(components,0) n_measurement = size(intensities,1) upsample = 10 projection = zeros((n_measurement, n_components)) projectionError = zeros((n_measurement, n_components)) dataSTD = shot.readoutNoiseRMS #detect plasma plasma = nanmax(intensities, axis = 0) > 10*dataSTD #překročení této hranice šumem je sttaisticky vellmi nepravděpodobné save_adv('./data/plasma', tvec,plasma) if not any(plasma): return True imin = amin(where(plasma)[0])-1 imax = amax(where(plasma)[0]) #print plasma #print 'tvec',imin,imax, isfinite(plasma_start), isfinite(plasma_end) #if isfinite(plasma_start) and isfinite(plasma_end): ##print 'tvec ', (arange(n_measurement)-imin)*(plasma_end-plasma_start)/(imax-imin)+plasma_start #tvec = (arange(n_measurement)-imin)*(plasma_end-plasma_start)/(imax-imin)+plasma_start+shot.integ_time/2000 plasma[1:] = plasma[1:] | plasma[:-1] plasma[:-1] = plasma[1:] | plasma[:-1] ##preprocess for i in where(plasma)[0]: intensities[:,i] = removeStrayLight(intensities[:,i]) intensities /= dataSTD #normalize by staistical error components_reshaped = reshape(components.T, (1,-1), order = 'F')[0,:] upsampled, retrofit = upsample_smooth(components_reshaped, 30,1, sampleInteg = True, upsample = upsample) upsampled_components = reshape(upsampled, (n_features*upsample,-1), order = 'F').T spect_lines = sum(components,axis = 0) != 0 #TODO 60% je nulových => jde to tím urychlit def estimateResid(components ): resid = zeros(n_measurement) for j in where(plasma)[0]: normalPx = ~overburned[:,j] _,resid[j] = nnls(components[:,normalPx].T,intensities[normalPx,j]) chi2 = norm(resid)**2/(sum(plasma)*size(intensities,0)) return chi2 def shiftSpectra(components,(s1,s2)): n = size(components,1) shift_components = interp1d(linspace(s2,n-s2-1,n),components, bounds_error=False,fill_value=0)(s1+arange(n)) return shift_components def fitfun(params,components_corr): shift_components = shiftSpectra(components_corr, params) chi2 = estimateResid(shift_components) #print chi2 return chi2 total_intens = sum( intensities[:,plasma], axis = 1) L = 10 ker = zeros(L) ker[L/2-1] = 1 ker_full = zeros(n_features) ker_full[n_features/2] = 1 normalPx = sum(overburned ,1) == 0 components_corr = zeros_like(components) print 'calc inst fun' for i in range(5): for j in range(n_components): components_corr[j,:] = fftconvolve(components[j,:],ker_full,mode='same') proj, resid = nnls(components_corr[:,normalPx].T,total_intens[normalPx]) retrofit = dot(proj.T,components).T ker,ker_full = deconv_opt(total_intens,retrofit,~normalPx,2,ker) t = time.time() par0 = (0,0) y = fmin_powell(fitfun,par0,args=[components_corr,], xtol=1e-4,ftol=1e-3,maxfun=1e6,disp=False) components_corr = shiftSpectra(components_corr, y) proj, resid = nnls(components_corr[:,normalPx].T,total_intens[normalPx]) #print resid**2/n_features/sum(plasma) print 'calc inst fun finished' components = components_corr components[components < 0] = 0 difference = zeros(n_features) chi2ion = zeros(n_components) resid = ones(n_measurement)*sqrt(n_features) q,r = qr(components.T,mode = 'economic') F = matrix(dot(inv(r),q.T)) for j in where(plasma)[0]: normalPx = ~overburned[:,j] projection[j, :], resid[j] = nnls(components[:,normalPx].T, intensities[normalPx,j]) #try to account error nonstatistical errors in fit difference[normalPx] = intensities[normalPx,j]-dot(components.T,projection[j, :])[normalPx] difference[~normalPx] = amax(difference[normalPx] ) #prostě tam bude velká chyba DataError = sqrt(difference**2+1)-(sqrt(2)-1) projectionError[j, :] = sqrt(diag(F*diag(DataError**2)*F.T)) for i in range(n_components): chi2ion[i] += norm(difference*(components[i,:]/(1e-3+sum(components,axis = 0))))**2/n_features #plot(wavelength, (dot(diag(projection[j, :]),components)).T,linewidth = 0.2) #plot(wavelength, intensities[:,j],'b-.') #plot(wavelength,dot(components.T,projection[j, :]), 'k') #plot(wavelength[:], intensities[:,j]-dot(components.T,projection[j, :]), 'r:' ) #ylim(-2, 100) #show() spectra = sum( intensities[:,plasma], axis = 1) retrofit = sum( dot(projection,components).T, axis = 1) chi2 = norm(spectra-retrofit)**2/n_features/sum(plasma) print 'standart chi2',chi2 normalPx = sum(overburned ,1) == 0 #calculate global charasteristics of the fit and spectra chi2 = sum(resid**2)/(n_measurement*n_features) print chi2, resid**2/n_features print chi2,chi2ion savetxt('./results/IonChi2'+'.txt',chi2ion) savetxt('./results/TimeChi2'+'.txt',resid**2/n_features) save('./results/TotalChi2',chi2) energy_constants = zeros(n_components) for i in range(n_components): norm_comp = spectrometr.convertCountToPhotons(shot.wavelength, copy(components[i,:]),shot.integ_time) energy_constants[i] = sum(norm_comp) abs_calibration = loadtxt('absolute_calibration.txt') energy_constants *= abs_calibration rel_proj = copy(projection) rel_proj_err = copy(projectionError) projection*= energy_constants projectionError *= energy_constants P_total = sum(projection,1) P_total_error = sqrt(sum(projectionError**2,1)) if isfinite(plasma_start) and isfinite(plasma_end): mean_tvec = sum(P_total*tvec)/sum(P_total) print mean_tvec print (plasma_start+plasma_end)/2 tvec = tvec-mean_tvec+(plasma_start+plasma_end)/2-shot.integ_time/1e3/2 save_adv('./data/TotalPower',tvec,P_total,data_err=P_total_error, tvec_err = 1e-3 ) n_frames = (plasma_end-plasma_start)/(shot.integ_time/1000) saveconst('./HistoricalAnalysis/meanPower',sum(P_total[plasma])/n_frames) save_adv('./data/projection', tvec,rel_proj, data_err=rel_proj_err, tvec_err=1e-3 ) save_adv('./data/intensities', tvec,intensities) save_adv('./data/components', wavelength,components) save_adv('./data/plasma', tvec,plasma) save('./data/energy_constants',energy_constants)