#!/usr/bin/env python
# -*- coding: utf-8 -*-
# ====================== GEFIT 0.1 (Golem EFIT) ======================================
# plama position calculation is based on the least square fit of the signal from mirnov coils a plasma current.
#A position and new plasma current are calculated from the furmula for loop wire is calculeted mag, filed in the coils and compared
#with measured values. It should by robust and it was tested on the poor data from mirnov coils.
#magnetic filed from external coils can be calculated by setting extCoilsPos and extCoilsCurr, however it was not tested yet.
#
# Autor: Tomáš Odstrčil tomasodstrcil@gmail.com
#TODO problémy: při velké disrupci dojde k tak masivnímu skoku že selže
#integrace a nastane integrační drift na některé z cívek. Udělá to chybu v signálu i 20%
import time
#t = time.time()
from numpy import *
from numpy.linalg import norm,lstsq
#from scipy.interpolate import interpolate
from scipy.optimize import minimize
from multiprocessing import Process, Pool, cpu_count
#from Deconvolution import *
from pygolem_lite.modules import deconvolveExp
import sys
from MagFieldCalc import *
#print 'import time PP:', time.time()-t
#constants
a = 0.085 #[m]
R_0 = 0.4 #[m]
#http://www.phy.uct.ac.za/courses/python/examples/fitresonance.py
#http://www.scipy.org/doc/api_docs/SciPy.optimize.minpack.html
def directConductorApprox(B):
R = (B[3]-B[1])/(B[3]+B[1])*-0.093
Z = (B[4]-B[2])/(B[4]+B[2])*-0.093
r = hypot(Z,R)
if r > a:
Z *= a/r
R *= a/r
R+= R_0
return R,Z
def fun(x,param,B_coils,y,Bout):
#print x.shape, param
B = Bloop_analytic(x[:,0],x[:,1],param[0]+R_0,param[1],Bout)*param[2]
B += B_coils #NOTE musí být správná orientace proudů! meřených flukama
B *= x[:,2:]
y[1:] = sqrt(sum(square(B),1))
y[0] = abs(param[2])
return y
def MainCycle((tvec,data,x,sigma,Bext,extCoilsCurr)):
#print data.shape
#data = double(data)
#x = double(x)
#tvec = double(tvec)
#sigma = double(sigma)
#Bext = double(Bext)
#extCoilsCurr = double(extCoilsCurr)
#print data.shape
#exit()
N = size(data,0)
n_det = size(data,1)
position = empty((N, 2))
retrofit = empty((N,n_det))
reziduum = empty(N)
I0 = median(data[:,0])
steps = 0
y = empty(size(x,0)+1)
Bout = empty((size(x,0),2))
for t in xrange(N):
#print t
#sys.stdout.write('\r'+str(t))
#sys.stdout.flush()
if all(isfinite(sigma)):
R,Z = directConductorApprox(data[t,:])
R+=0.04
else: R,Z = R_0,0
scale = 1/array((1,1,0.5/I0)) #1/a/2
params0 = array([R-R_0,Z,data[t,0]])
if Bext == None:
B_coils = 0
else:
B_coils = dot(Bext,extCoilsCurr[t,:])
T = time.time()
#def fitfun(params):
#retrofit = fun(x,params*scale, B_coils,y,Bout)
#retrofit-= data
#retrofit/= sigma
#retrofit *=retrofit
#res = norm(retrofit)**2
#res+= norm((params*scale)[:2]/(a*1.0))**50
##res = norm((data[t,:] - fun(x,params*scale, B_coils,y,Bout))/sigma)**2+norm((params*scale)[:2]/(a*1.0))**50
#return res
fitfun = lambda params : norm((data[t,:] - fun(x,params*scale, B_coils,y,Bout))/sigma)**2+norm((params*scale)[:2]/(a*1.0))**50
res = minimize(fitfun, params0/scale, method='Nelder-Mead',options={'xtol':1e-4, 'ftol' : 1e-2})
#res = minimize(fitfun, params0*scale, method='BFGS',options={'xtol':1e-4, 'ftol' : 1e-2})
#print '\n', time.time()-T
#exit()
#sys.stdout.write('\r'+str(time.time()-T))
#sys.stdout.flush()
if not res.success:
res.x[:] = params0
res.fun = 1e6
print 'position calculation failure in ', str(tvec[t]), 's'
else:
params0 = res.x*scale
position[t,:] = (res.x*scale)[:2]
position[t,0]+= R_0
steps+= res.nfev
retrofit[t,:] = fun(x,params0,B_coils,y,Bout)
reziduum[t] = res.fun +sum(isinf(sigma))*1e4
return (position,retrofit,reziduum)
#http://www.phy.uct.ac.za/courses/python/examples/fitresonance.py
#http://www.scipy.org/doc/api_docs/SciPy.optimize.minpack.html CalcPlasmaPosition(tvec_ds,detectorPos, detectorSignal_ds, detectorDriftError, Ip_ds,IpDriftError, shot_num = shot_num)
def CalcPlasmaPosition(tvec,detectorPos, detectorSignal,detectorDriftError,
Ip,IpDriftError,extCoilsPos = None,extCoilsCurr= None ):
#očekávám stejné časové rozdělení vstupního signálu, všechny položky jsou už absolutně kalibrované, jednotky SI
print 'CalcPlasmaPosition'
t1 = time.time()
n = len(tvec)
coil_perpend=detectorPos-mean(detectorPos,axis=0)
coil_perpend[:,1]*= -1
coil_perpend=coil_perpend[:,::-1]/hypot(coil_perpend[:,1],coil_perpend[:,0])[:,None]
x = hstack((detectorPos,coil_perpend))
data = vstack((Ip, detectorSignal.T)).T
sigma = hstack((IpDriftError, detectorDriftError))
if extCoilsPos == None or extCoilsCurr == None:
Bext = None
else:
Bext = zeros(shape(detectorPos)+(size(extCoilsPos,0),))
for i in arange(size(extCoilsPos,0)):
Bext[...,i] = Bloop_analytic(x[:,0],x[:,1],extCoilsPos[i,0],extCoilsPos[i,0])
#main calculation , multiprocess
n_cpu = cpu_count()
p = Pool(n_cpu)
#split_ind = array_split(range(n), n_cpu*4)
#print data.shape
#exit()
data_spl = array_split(data, n_cpu*4,axis=0)
tvec_spl = array_split(tvec, n_cpu*4)
inputs = [(tvec_spl[i],data_spl[i],x,sigma,Bext,extCoilsCurr) for i in range(n_cpu*4)]
list_data = p.map(MainCycle,inputs )
p.close()
p.join()
list_position = list()
#list_position_dc = list()
list_retrofit = list()
list_reziduum = list()
for (position,retrofit,reziduum) in list_data:
list_position.append(position)
#list_position_dc.append(position_dc)
list_retrofit.append(retrofit)
list_reziduum.append(reziduum)
position = vstack(list_position)
#position_dc = vstack(list_position_dc)
retrofit = vstack(list_retrofit)
reziduum = hstack(list_reziduum)
radius = a-sqrt((position[:,0]-0.4)**2+position[:,1]**2)
chi2 = mean((data-retrofit)**2,axis=0)/sigma**2
print 'position calculated in %g s' % (time.time()-t1)
print 'chi2: ', chi2,' total ', norm(chi2)
return tvec, position, radius, reziduum, retrofit, data,chi2
def RemoveDriftsAuto(signal, Bt,Uloop,plasma_start,plasma_end,Bt_trigger,E_trigger, Bt_crosstalk=True, Uloop_crosstalk=False, Stabil_crosstalk=False, Trafo_effect=True ):
n_det = size(signal, 1)-1
index = 2
t_max = min(Bt[-1,0],Uloop[-1,0],signal[-1,0]) # čas po který to ty diagnostiky berou
t_min = max(Bt[0,0],Uloop[0,0],signal[0,0])
dt = mean(diff(signal[:,0]))
signal = signal[:int(t_max/dt+1),:]
tvec = signal[:,0]
n = size(tvec)
#resample to the same resolution
if Uloop_crosstalk:
i_Uloop = index
index+=1
Uloop = interp(tvec, Uloop[:,0],Uloop[:,1], left=0, right=None)
if Bt_crosstalk:
i_Bt = index
index+=1
Bt = interp(tvec, Bt[:,0],Bt[:,1], left=0, right=None)
if Trafo_effect:
t_exp = 21.2e-3 #[s]
i_Traf_res = index
Bt_conv,retrofit = deconvolveExp(Bt,-t_exp, dt,0,0)
index+=1
if Stabil_crosstalk:
i_St = index
index+=1
Bt = interp(tvec, fluke[:,0],fluke[:,1], left=0, right=None)
if isnan(plasma_start) or isnan(plasma_end):
plasma_start = Bt_trigger
plasma_end = Bt_trigger
#if removePlasmaDrift:
#drift_error = mean(cumsum((diff(signal[:,1:], axis = 0))**2, axis = 0), axis = 1)
#drift_error = hstack((0,drift_error))
#plot(drift_error)
#show()
#plot(diff(signal[:,1:], axis = 0))
#show()
#i_drift = index
#index+=1
Bt_trigger_n = argmin(abs(tvec-Bt_trigger))
E_trigger_n = argmin(abs(tvec- E_trigger))
plasma_start_n = max(argmin(abs(tvec-plasma_start)),E_trigger_n)
plasma_end_n = max(argmin(abs(tvec-plasma_end)),E_trigger_n)
interval0 = range(int(t_min/dt),Bt_trigger_n) #clean area before any trigger
interval1 = range(Bt_trigger_n+int(0.5e-3/dt), E_trigger_n) # short time between Bt and CD trigger
interval2 = range(E_trigger_n,plasma_start_n) #interval between CD trigger and breakdown
interval3 = range(plasma_end_n,n ) # interval beween plasma end and torroidal field end tyristor
interval_plasma = range(plasma_start_n,plasma_end_n) #plasma
interval = interval1+interval2+interval3 #interval used for removing of the drift
#print i
mainBasis = zeros((n, index))
if Bt_crosstalk:
mainBasis[:,i_Bt] = Bt
if Uloop_crosstalk:
mainBasis[:,i_Uloop] = Uloop
if Trafo_effect:
mainBasis[:,i_Traf_res] = Bt_conv
if Stabil_crosstalk:
mainBasis[:,i_St] = fluke
#if removePlasmaDrift:
#mainBasis[:,i_drift] = drift_error
#offset, integration drift
offset = zeros(n)
offset[interval0] = 1
mainBasis[:,0] = cumsum(offset, out = offset)
offset = ones(n)
offset[interval0] = 0
mainBasis[:,1] = cumsum(offset, out = offset)
#offset = zeros(n)
#offset[interval_plasma] = 1
#mainBasis[:,2] = cumsum(offset, out = offset)
x,res,rank,s = lstsq(mainBasis[interval,:],signal[interval,1:])
corr_signal = signal[:,1:]-dot(mainBasis, x)
res = sum((corr_signal[interval,:])**2,0)
res*= 1e15/len(interval)
res = sqrt(res)
#prepare graph
frames = list()
vlines = list()
hlines = list()
rectangles = list()
vlines.append((plasma_end*1000,'--' ))
vlines.append((plasma_start*1000,'--' ))
hlines.append(( 0, '-.'))
for i_det in arange(n_det):
curves = list()
curves.append((single(1e6*signal[:,i_det+1]),'Original signal', '--' ))
curves.append((single(1e6*corr_signal[:,i_det]), 'Corrected signal', 'k-'))
curves.append((single(1e6*dot(mainBasis[:,:2],x[:2,i_det])),'Integration offset','-' ))
if Bt_crosstalk:
curves.append((single(1e6*mainBasis[:,i_Bt]*x[i_Bt,i_det]), 'Toroidal mag. field', '-'))
if Trafo_effect:
curves.append((single(1e6*mainBasis[:,i_Traf_res]*x[i_Traf_res,i_det]),'Trafo mag. field','-'))
if Uloop_crosstalk:
curves.append((single(1e6*mainBasis[:,i_Uloop]*x[i_Uloop,i_det]),'Uloop','-'))
frames.append((1000*tvec,curves,'mc'+str(i_det*4+1)+' [ mV$\cdot$ms]','$\\chi^2$ = %2.1f'%res[i_det]))
intervals = (interval1,interval2,interval3 )
for inter in intervals:
if len(inter) == 0:
continue
(x_min,x_max) = (tvec[inter][0]*1e3, tvec[inter][-1]*1e3)
rectangles.append((x_min,x_max))
AutoRemoveGraph = (frames, vlines,hlines,rectangles)
#save projection constants
proj_file = 'mc1\tmc5\tmc9\tmc13\n'
proj_file+='drift1 [V] '+str(x[0,:])+'\n'
proj_file+='drift2 [V] '+str(x[1,:])+'\n'
if Bt_crosstalk:
proj_file+='Bt crosstalk [V/T]\t'+str(x[i_Bt,:])+'\n'
if Trafo_effect:
proj_file+='Trafo crosstalk (respons time %1.1f ms)[V/T]\t'%(1000*t_exp) +str(x[i_Traf_res,:])+'\n'
if Uloop_crosstalk:
proj_file+='Uloop (residual) crosstalk [V/V]\t'+str(x[i_Uloop,:])+'\n'
f = open('./constants/projections.txt', 'w')
f.write(proj_file)
f.close()
return vstack((tvec,corr_signal.T)).T,res,AutoRemoveGraph
