#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Algorithm for integrated data analyze of the visible light signal from GOLEM tokamak.
Many different diagnotics with different wavelength resolution, time resolution, time shift error,
and signal error is availible on the tokamak.This algirithm puts them together and finds simplest
solution which corespondes to the measured signal. Usually most of the wavelength information is provided
by spectrometer and the higth freqency part is provided by fast cameras and photodiodes.
The problem is solved as least squares fit of the data by nonparametrical model sim maximal simplisity.
Autor: Tomas Odstrcil
date:12.12.2012
"""
#TODO jak jsou ošetřená přepálená data z foťáku?? dát tam inf???
from numpy import *
from matplotlib.pyplot import *
from SpectrometerControl import *
from scikits.sparse.cholmod import cholesky, analyze, cholesky_AAt,CholmodError
from scipy.signal import fftconvolve,gaussian
from numpy import linalg
import time
from scipy import sparse
from numpy.matlib import repmat
from scipy.sparse import *
from scipy.linalg import norm, inv,qr,pinv
from pygolem_lite.modules import save_adv,load_adv
from CalcIonProjections import Resample,upsample_smooth,GapsFilling
import numexpr as ne
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 WeightedMean(values,error,correlMatrix = None): #implicit expectation of the gaussin error distribution
n = len(values)
if correlMatrix == None:
correlMatrix = np.identity(n)
covarMatrix = (diag(error)*matrix(correlMatrix**2)*diag(error))
icovarMatrix = linalg.pinv(covarMatrix)
s = 1/sum(icovarMatrix)
m = s*sum(values*icovarMatrix)
chi2 = double((values-m)*((values-m)*icovarMatrix).T)
s*= chi2/(n-1)
return m,sqrt(s)
def WeightedMeanVec(vectors,vec_errors,correlMatrixes): #implicit expectation of the gaussin error distribution
n = len(vectors[0])
#http://en.wikipedia.org/wiki/Weighted_mean#Vector-valued_estimates
s = zeros((n,n))
x = zeros((n,1))
for (cor,err,vec) in zip(correlMatrixes,vec_errors,vectors):
cov = diag(err)*matrix(cor**2)*diag(err)
icov = linalg.pinv(cov)
s += icov
x += dot(icov,vec).T
s = linalg.pinv(s)
x = dot(s,x).T
chi2 = 0
for (cor,err,vec) in zip(correlMatrixes,vec_errors,vectors):
cov = diag(err)*matrix(cor**2)*diag(err)
icov = linalg.pinv(cov)
chi2 += double((vec-x)*((vec-x)*icov).T)
s*= chi2/len(vectors)
return squeeze(x),sqrt(diag(s))
def IntegLightAnalyz(wavelength,plasma_start,plasma_end,SpectrometerData,PhotodiodesData, upsample = 10):
print 'IDA start'
wavelength,components = load_adv('./data/components') #the most critical part, weeks of my work
(spec_tvec,spec_tvecError,SpectrSensitiv,projection,projectionError) = SpectrometerData
(photoData,photoError,correlList,photoTvec,photoTvecPrecis,photoFilter,photoSensitivity,photoLabel) = PhotodiodesData
#photoFilter - musí být všechny stejně normované aby se nepokazily vzájemé citlivosti
SpectrSensitiv = interp(wavelength, SpectrSensitiv[:,0], SpectrSensitiv[:,1])
for i,(pfilter,rsPhotodiod) in enumerate(zip(photoFilter,photoSensitivity )):
photoFilter[i] = interp(wavelength,pfilter[:,0],pfilter[:,1],left=0,right = 0)
photoSensitivity[i]=interp(wavelength,rsPhotodiod[:,0],rsPhotodiod[:,1],left=0, right=0)
n_features = size(components,1)
n_components = size(components,0)
n_measurement = size(projection,0)
n_photodiodes = len(photoData)
dt = mean(diff(spec_tvec))/upsample
#this matrix make a projection of the spectrometer data to the photodiodes data
FilterProjMatrix = empty((n_components, n_photodiodes))
for i,(pfilter,rsPhotodiod) in enumerate(zip(photoFilter,photoSensitivity)):
FilterProjMatrix[:,i] = dot(components,1/SpectrSensitiv*pfilter*rsPhotodiod)
photoData_down = list()
photoError_down = list()
photoTvec_down = list()
#downsample data from photodiods
for i,(pdata,perror,ptvec) in enumerate(zip(photoData,photoError,photoTvec)):
#nanind = array(~isfinite(pdata))
if not all(isfinite(pdata)):
perror,_ = GapsFilling(perror)
pdata,_ = GapsFilling(pdata)
t_min = amin(spec_tvec)
t_max = amax(spec_tvec+dt*(upsample-1))
minerr= amin(perror)
#plot(perror)
#savefig('err%d.png'%i)
#clf()
tvec,data = Resample(ptvec, pdata,n_measurement*upsample,t_min=t_min,t_max=t_max)
tvec,err = Resample(ptvec,perror,n_measurement*upsample,
t_min=t_min,t_max=t_max,left=minerr,right=minerr)
#plot(ptvec,nanind)
#plot(ptvec,pdata)
#plot(tvec,data)
#savefig('data'+str(i)+'.png')
#clf()
photoTvec_down.append(tvec)
photoData_down.append(data)
photoError_down.append(err)
##_,nanind = Resample(ptvec, nanind,n_measurement*upsample,t_min=t_min,t_max=t_max)
#plot(photoTvec[i],nanind)
#nanind = nanind > 0.3
#plot(photoTvec[i],nanind)
#xlim(plasma_start,plasma_end)
#savefig('nanind.pdf')
#(photoError[i])[nanind] = 1e6 #big number insted of infinity
#upsample spectrometer data
upsampled_proj = empty((n_components, upsample*n_measurement))
for i, proj in enumerate(projection.T ):
upsampled_proj[i,:],_ = upsample_smooth(proj,15,1,sampleInteg = True,upsample = upsample)
#calculate expected photodiods signal from spectrometer
PhotodFromSpec = dot(FilterProjMatrix.T,upsampled_proj )
PhotodFromSpec2 = dot(FilterProjMatrix.T,projection.T )
#calculate the mutual renormalization of the intensity and time shift of the time axis if the diagnostics
shift = zeros(n_photodiodes)
amplitude = zeros(n_photodiodes) #TODO měly by být taky všechny stejné, až to bude dobře zkalibrované
for i, (photod,spect) in enumerate(zip(photoData_down,PhotodFromSpec)):
convolution = fftconvolve(photod,spect[::-1], mode = 'same')
shift[i] = len(convolution)/2-argmax(convolution)
amplitude[i] = sum(spect)/sum(photod)
tshift = shift*dt
#estimate correction of the spectrometer timeaxis
correlList.insert(0,1) # spectrometer
M = sparse.block_diag(correlList)
correlMatrix = M.todense()
spec_tshift,spec_tshift_err = WeightedMean(np.hstack((0,tshift)),
np.hstack((spec_tvecError,photoTvecPrecis)),
correlMatrix=correlMatrix)
spec_tvec += -spec_tshift+dt*upsample/2 #korekce subixelové chyby
tshift -= spec_tshift
#correstion of the timeaaxis of the photodiods
tshift,tshift_err = WeightedMeanVec((zeros(n_photodiodes),tshift),
(photoTvecPrecis,spec_tshift_err*ones(n_photodiodes)),
(correlMatrix[1:,1:],correlMatrix[1:,1:]))
tshift = array(tshift,ndmin=1)
tshift_err = array(tshift_err,ndmin=1)
#make a coresponding shift of the photodiod signal
for i, (shift,shift_err,photod,photodE,photoT,spect) \
in enumerate(zip(tshift,tshift_err,photoData_down,photoError_down,photoTvec_down,PhotodFromSpec)):
err_shift = hypot(shift_err,spec_tshift)
photoTvec[i] += shift
photoData[i] *= amplitude[i]
shift += spec_tshift
photod[:] = interp(photoT,photoT+shift,photod )*amplitude[i]
photodE[:] = interp(photoT,photoT+shift,photodE)*amplitude[i]
photodE[1:]+= err_shift*abs(diff(photodE))/dt #error in x variable
photoT -= spec_tshift
#plot(photoData[i])
#plot(spect)
#savefig('_'+str(i)+'.png')
#clf()
#errorbar(photoT,photod,photodE)
#title(photoLabel[i])
#plot(spec_tvec,PhotodFromSpec2[i,:])
#axvline(x = plasma_start,c = 'r')
#axvline(x = plasma_end,c = 'r')
#savefig(str(i)+'.png')
#clf()
# x is solution
n_b = n_components*n_measurement
n_x = n_b*upsample
x = zeros(n_x)
# b are measured data
b1 = reshape(projection, (-1,1),order='F') #data from spectrometer
b2 = array(photoData_down,copy=False) #data from photodiode
b2 = reshape(b2, (-1,1),order='C')
b = squeeze(np.vstack((b1,b2)))
#e is expected error in data_file
e1 = reshape( projectionError, (-1,1), order='F') #error from spectrometer
e2 = array(photoError_down,copy=False) #error from photodiode
e2 = reshape( e2, (-1,1),order='C')
e = squeeze(np.vstack((e1,e2)))
e+= 1e-4 #nonzeroes errors
#The vogel matrix
ReductMatrix = kron(eye(n_b,n_b),ones((1,upsample))/upsample,format='csr')
FiltMatrix = kron(FilterProjMatrix.T,
identity(n_measurement*upsample,format='dia'),format='csr')
T = vstack((ReductMatrix, FiltMatrix))
iE = diags(1/e, 0)
T = iE*T
b = iE*b
tvec = photoTvec_down[0]
#Smoothing matrix
interv_noplasma = (tvec< plasma_start)+(tvec> plasma_end)
interv_noplasma = squeeze(repmat(interv_noplasma,1,n_components))
diag_data = ones((2,n_x))
diag_data[1,:] *= -1
diag_data[1,interv_noplasma] = 0
std = sum(~interv_noplasma)/50
gauss_win = gaussian(std*3,std)
gauss_win/= sum(gauss_win)
diag_data[1,:] = fftconvolve(diag_data[1,:],gauss_win,mode='same')
diag_data[1,::(n_measurement*upsample)] = 0
D = spdiags(diag_data, (0,1), n_x, n_x,format='csr')
D = D.T*D
DD = D.T*D
TT = T.T*T
W = ones(n_x) #weight matrix
factor = analyze(TT+DD)
for j in range(5):
W = D.T*sparse.spdiags(W,0, n_x,n_x)*D
lam = 1e4 #initial guess, just estimated "by eye"
for i in range(5):
factor.cholesky_inplace(TT +lam*W)
g = squeeze(factor(T.T*b))
chi2 = norm(T*g-b)**2/len(b)
print lam, chi2
if chi2 < 1.3: #very simple root finder
lam*= 2
else:
break
#elif chi2 < 1:
#lam*= 3
#else:
#break
g_tmp = copy(g)
g_tmp[g_tmp < 0.001] = 0.001 # remove negative points (result must be positive)
W = mean(g_tmp)/(g_tmp) # new weight matrix
f = T*g
resid = f-b
chi2 = norm(resid)**2/len(b)
print 'chi2',chi2
E = diags(e, 0)
g = reshape(g, (-1,n_components),order='F')
save_adv('./data/IDA_ions_projections',tvec,g )
#exit()
#plot results
fig = figure(num=None, figsize=(8, 6), dpi=80, facecolor='w', edgecolor='k')
ax = fig.add_subplot(111)
for i in range(n_components):
plot(tvec*1e3,g[:,i], styles[i], label=names[i])
ax.text(0.9,0.9,'$\\chi^2/doF = %2.1f$'%chi2,horizontalalignment='right',
verticalalignment='bottom',transform=ax.transAxes,backgroundcolor='w')
leg = legend(loc='upper left', fancybox=True)
leg.get_frame().set_alpha(0.7)
axvline(x = plasma_start*1e3,ls = '-.',c = 'g',)
axvline(x = plasma_end*1e3 ,ls = '-.',c = 'g',)
ax.set_xlabel('t [ms]')
ax.set_ylabel('relative intensity [-]')
ax.axis('tight')
ax.set_xlim(plasma_start*1e3-1,plasma_end*1e3+1)
fig.savefig('./graphs/IDA_ions_projections.png',bbox_inches='tight')
fig.clf()
#plot retrofit
retroSpec = (E*f)[:n_measurement*n_components]
retroSpec = reshape(retroSpec, (n_components,-1),order='C')
chi2 = norm(resid[:n_b])**2/n_b
space = 0.1
ax = fig.add_subplot(111)
for i in range(n_components):
plot(spec_tvec*1e3,retroSpec[i,:]+i*space, styles[i], label=names[i])
ax.text(0.9,0.9,'$\\chi^2/doF = %2.1f$'%chi2,horizontalalignment='right',
verticalalignment='bottom',transform=ax.transAxes,backgroundcolor='w')
ax.set_xlabel('t [ms]')
ax.set_ylabel('shift+intensity [a.u.]')
leg = legend(loc='upper left', fancybox=True)
leg.get_frame().set_alpha(0.7)
for i in range(n_components):
errorbar(spec_tvec*1e3,projection[:,i]+i*space,
yerr=projectionError[:,i],fmt=styles[i][0]+'.')#intensity with errorbars
ax.axis('tight')
ax.set_xlim(plasma_start*1e3-1,plasma_end*1e3+1)
fig.savefig('./graphs/spectrometer_retrofit.png',bbox_inches='tight')
fig.clf()
retroDiods = (E*f)
cm = get_cmap('gist_rainbow')
color = (cm(i/(1.+n_photodiodes)) for i in range(n_photodiodes+1))
ax = fig.add_subplot(111)
chi2 = norm(resid[n_b:])**2/(len(b)-n_b)
ymax = 0
index = n_measurement*n_components
for (ptvec,pdata,perr,ptvec_full,pdata_full,lab,pc) in zip(photoTvec_down,photoData_down,photoError_down,photoTvec,photoData,photoLabel,color):
plot(ptvec_full*1e3,pdata_full, color=pc, linestyle='-')
fill(np.concatenate([ptvec*1e3, ptvec[::-1]*1e3]), \
np.concatenate([pdata - perr,(pdata + perr)[::-1]]), \
alpha=.4, fc=pc, ec='None')
ymax = amax(pdata) if amax(pdata) > ymax else ymax
plot(ptvec*1e3,retroDiods[index:index+len(pdata)]+index/100.,'--', color=pc,label=lab)
index+= len(pdata)
#fig.savefig('./graphs/photodiodes_retrofit'+str(index)+'.png',bbox_inches='tight')
ax.text(0.9,0.9,'$\\chi^2/doF = %2.1f$'%chi2,horizontalalignment='right',
verticalalignment='bottom',transform=ax.transAxes,backgroundcolor='w')
leg = legend(loc='upper left', fancybox=True)
leg.get_frame().set_alpha(0.7)
ax.set_xlim(plasma_start*1e3-1,plasma_end*1e3+1)
ax.set_ylim(-ymax/20,ymax*1.1)
ax.set_xlabel('t [ms]')
ax.set_ylabel('intensity [a.u.]')
fig.savefig('./graphs/photodiodes_retrofit.png',bbox_inches='tight')
fig.clf()
#exit()
