Source code :: IDA

#!/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()