#!/usr/bin/python2
# -*- coding: utf-8 -*-
""" CREATED: 7/2012
AUTHOR: MICHAL ODSTRCIL
"""
from numpy import *
from config import *
import sys
#try:
#from matplotlib.pyplot import *
#except:
#pass
def import_matplotlib():
import matplotlib
matplotlib.rcParams['backend'] = 'Agg'
matplotlib.rc('font', size='10')
matplotlib.rc('text', usetex=True) # FIXME !! nicer but slower !!!
import matplotlib.pyplot as plt
return plt
try:
plt = import_matplotlib()
except:
pass
import os, re, datetime, time
from pygolem_lite import *
from copy import deepcopy
from pygolem_lite import isnone
###################### SIMPLE DATA LOADING ##############################
def cat(path, lines = [0], return_array = False):
""" read file from `path` and return string. If there is several lines, it can return array if `return_array == True`
"""
#try:
with open(path, 'r') as f:
content_tmp = f.readlines()
if return_array: return content_tmp
if len(lines) == 0:
lines = range(len(content_tmp))
content = ""
for i in lines:
content += content_tmp[i]
#except:
#content = 'N/A'
return content #return the raw string
####################### OTHER USEFUL ROUTINES #################################x
def nanmedian(data):
return median(data[~isnan(data)])
def list2array(L):
"""
convert list output from pygolem to array / scalar
"""
if type(L) is not list and type(L) is not tuple :
return L
if ndim(L[1]) == 1:
dim = 1
else:
dim = size(L[1],1)
out = empty([len(L[0]), dim+1])
out[:,0] = L[0]
out[:,1:] = reshape(L[1], (-1, dim))
return out
def random_number_downsample(tvec, data, N):
"""
downsample data without moire effects
"""
Ndata = len(tvec)
if Ndata <= N:
return tvec, data
ind = int_( linspace(0,Ndata-1, N))
ind_rand = ind + (random.randint(0,Ndata/N, size = (N)) if Ndata/N > 0 else 0)
ind_rand[-1] = Ndata-1
ind_rand = unique(ind_rand)
return tvec[ind_rand], data[ind_rand,:]
###################### PLOTTING ##############################
def multiplot(data_all, title_name, file_name, figsize = (9,None), dpi = 100, orientation = 'vertical', file_type = 'png', reduction = False, debug_info = True):
""" Plot several graphs to on multiplot. Use data `data_all` provided by function get_data !!
:var dict data_all: Input data or list of input data
:var str title_name: Title of the image
:var str file_name: Name of saved image (without filetype)
:var list figsize: Size of the output image
:var int dpi: DPI of output image
:var str orientation: Choose multiplot orientation: vertical/horizontal
:var str file_type: Image file type, recomended: png, svgz !! svgz is slow !!
:var bool reduction: Decrease number of plotted points => speed up plotting, use for all nice smooth lines !!
:var bool debug_info: Allow printing of some debug info
"""
from matplotlib.ticker import MultipleLocator, FormatStrFormatter,LogLocator, MaxNLocator
plt = import_matplotlib()
class MyFormatter(plt.ScalarFormatter): # improved format of axis
def __call__(self, x, pos=None):
self.set_scientific(True)
self.set_useOffset(True)
self.set_powerlimits((-2, 3))
#print x, pos, ScalarFormatter, ScalarFormatter()
if pos==0:
return ''
else:
return plt.ScalarFormatter.__call__(self, x, pos)
is_icon = dpi <= 50
if all(isnone(data_all)):
#print "data_all", data_all
raise IOError('Empty input data set for plotting')
else:
data_all = deepcopy(data_all) # avoid mysterious changes in plots
if type(data_all) is dict:
data_all = [data_all]
i = 0 # !! remove empty fields !!
while i < len(data_all):
if len(data_all[i]) == 0 or all(isnone( data_all[i] )) :
data_all.pop(i)
else:
i += 1
n_plots = len(data_all)
t0 = time.time()
if figsize[1] is None: # if the size is not determined, make it flexible
figsize = array(figsize)
figsize[1] = 3*n_plots
fig = plt.figure(num=None, figsize=figsize, edgecolor='k')
plt.subplots_adjust(wspace=0.3,hspace=0)
for i,data in enumerate(data_all):
if orientation is 'vertical':
ax = fig.add_subplot(n_plots,1,i+1)
ax.xaxis.set_major_formatter( plt.NullFormatter() )
ax.yaxis.set_major_formatter( MyFormatter() ) # improved format of axis
else:
ax = fig.add_subplot(1,n_plots,i+1)
if is_icon: # very small image
reduction = True
try:
_plot_adv(data, reduction, ax, dpi)
except Exception, e:
raise NotImplementedError("plotting failed", str(e), "\n" )# , data
try: # problem with log scale
ax.ticklabel_format(style='plain', axis='y')
ax.yaxis.set_major_locator( MaxNLocator(max(3,int(figsize[1]/double(n_plots)*5) )))
except:
pass
ax.tick_params(axis='both' if not is_icon else 'none',reset=False,which='both',length=2,width=0.5)
if i == 0 and title_name != "" and not is_icon:
ax.set_title(fix_str(title_name))
if orientation is 'horizontal':
plt.xlabel(data['xlabel'] if type(data) == dict else data[0]['xlabel'])
if orientation is 'vertical':
ax = fig.add_subplot(n_plots ,1,n_plots)
ax.xaxis.set_major_formatter( plt.ScalarFormatter())
#ax.yaxis.set_major_formatter( plt.ScalarFormatter())
#ax.xaxis.set_minor_locator( MultipleLocator(1) )
if not is_icon:
try:
plt.xlabel( fix_str(data['xlabel'] if type(data) == dict else data[0]['xlabel']) )
except Exception, e:
print "plotting xlabel failed", str(e), "\n", data
t0 = time.time()
file_name = re.sub('(.+)\.(.+)$', r'\1', re.sub('(.+)\.(.+)$', r'\1', file_name))
image_name = file_name+'.'+file_type
plt.savefig( image_name , dpi= dpi, bbox_inches='tight') # better but slower
plt.clf()
if file_type in ['png', 'jpg']:
os.system('convert '+ image_name +" -trim " + image_name )
if debug_info:
print "saving", time.time() - t0
print "plotting time %.2fs" % (time.time() - t0)
def paralel_multiplot(*kargs, **kwargs):
import multiprocessing
try:
multiprocessing.Process(target=multiplot, args=kargs, kwargs=kwargs).start()
except:
print "PLOTTING FAILED:" + kargs[0]
def _plot_adv(data, reduction, ax, dpi):
""" advanced plotter with dictionary input,
internal function, do not use directly
"""
if type(data) is dict:
data = [data]
is_icon = dpi <= 50
i = 0 # !! remove empty fields !!
while i < len(data):
if all(isnone(data[i])) or len(data[i]) == 0:
data.pop(i)
else:
i += 1
if len(data) == 0: # empty list
return
use_legend = False
N_data = len(data)
t0 = time.time()
for i in range(N_data):
if len(data[i]) == 0 or all(isnan(data[i]['data'])):
continue
d = data[i]
kwargs = d['kwargs']
ind_plot = ones(len(d['tvec']), dtype = bool)
if kwargs['label'] != "":
use_legend = True
if 'xlim' in d.keys(): # speed up plotting
ind_plot = where( (d['xlim'][0] <= d['tvec']) & ( d['xlim'][1] >= d['tvec']))
ind_plot = concatenate([ind_plot[0]-1, squeeze(ind_plot), ind_plot[-1]+1]) # extend the data by one point from the range
ind_plot = ind_plot[(ind_plot >= 0) & (ind_plot < len(d['tvec'])) ]
ind_plot = unique(ind_plot)
if reduction:
ind_plot = ind_plot[:: max(1, len(d['tvec'])/2000) ]
# downsampling of errobar data !!!
for name in ['xerr','yerr']:
if not kwargs[name] is None and not isscalar(kwargs[name]):
kwargs[name] = [kwargs[name][ind_plot,0], kwargs[name][ind_plot,1] ]
plt.errorbar(d['tvec'][ind_plot], d['data'][ind_plot], linewidth=0.5, capsize=0, **kwargs )
if data[0]['yscale'] != '':
ax.set_yscale(data[0]['yscale'])
if amax(data[0]['tvec']) - amin(data[0]['tvec']) < 100:
ax.xaxis.set_minor_locator(plt.MultipleLocator(1))
if 'vlines' in data[0].keys():
for line in data[0]['vlines']:
plt.axvline(x=line, color = 'b', linestyle='--')
for key in ['xstart', 'xend']:
if key in data[0].keys():
plt.axvline(x=data[0][key], color = 'k', linestyle='--')
if 'xlim' in data[0].keys():
Xrange =[min([data[i]['xlim'][0] for i in range(N_data)]),\
max([data[i]['xlim'][1] for i in range(N_data)])]
plt.xlim(Xrange)
if 'ylim' in data[0].keys():
ymin = inf ; ymax = -inf
for i in range(N_data):
ymin_tmp, ymax_tmp = _detect_range(data[i]['tvec'] , data[i]['data'], Xrange, data[i]['ylim'])
ymin = min(ymin_tmp,ymin)
ymax = max(ymax_tmp,ymax)
Yrange = [ymin, ymax]
plt.ylim(Yrange)
if not is_icon: # in case of small icon do no apply
if data[0]['ylabel'] != "" and data[0]['ylabel']:
plt.ylabel(data[0]['ylabel'])
if 'annotate' in data[0].keys():
for annot in data[0]['annotate']:
if annot['pos'][0] is None: annot['pos'][0] = Xrange[0] + (Xrange[1] - Xrange[0])*0.05
if annot['pos'][1] is None: annot['pos'][1] = Yrange[1] - (Yrange[1] - Yrange[0])*0.05
plt.annotate(annot['text'], annot['pos'], rotation=annot['angle'], fontsize=annot['fontsize'] )
if 'ygrid' in data[0].keys() and data[0]['ygrid']:
ax.yaxis.grid(color='gray', linestyle='dashed')
if 'xgrid' in data[0].keys() and data[0]['xgrid']:
ax.xaxis.grid(color='gray', linestyle='dashed')
if use_legend:
#leg = legend()
# FIXME !! nicer but slower !!!
leg = plt.legend(loc='best' , fancybox=True)
leg.get_frame().set_alpha(0.7)
return ax
def _detect_range(tvec, data, xlim, ylim):
""" smart way how to determine best ylim in the plots
internal fucntion do not use directly
"""
from scipy.stats.mstats import mquantiles
ind_ok = ~isnan(data) & ~isinf(data) & (tvec > xlim[0]) & (tvec < xlim[1])
yrange = [mquantiles(data[ind_ok], 0.005), mquantiles(data[ind_ok], 0.995)]
dy = abs(yrange[1] - yrange[0])
if dy/amax(abs(array(yrange))) < 1e-6: # constant line ...
yrange[0] = min(0, yrange[0]*1.1 )
yrange[1] = max(0, yrange[1]*1.1 )
else:
yrange[0] -= 0.2*dy
yrange[1] += 0.2*dy
for i in [0,1]:
if ylim[i] is None:
ylim[i] = yrange[i]
return ylim
def get_data(diagn, diagn_name = None, ylabel = None, xlabel = 'Time [ms]', shot_num = None, xlim = [], ylim = [], \
tvec_rescale = 1e3, data_rescale = 1, integrated = False, columns = [], yscale = '', \
reduction = False, plot_limits = True, smoothing = 0, line_format = "", \
vlines = [], annotate = [], xgrid= False, ygrid = False, **kwargs ):
"""
Load data/diagn and return object that is used for plotting.
:var str/list/Data class diagn: Input can be string with name of diagn from main database, name of file with locally saved data on disk, list with [tvec,data]
:var str diagn_name: Name of diagnostics that will be used in legend, if diagn has more signal than name will be `name 0`, `name 1`, ...., or diagn_name can be list with name for each column
:var str ylabel: Title of ylabel
:var str xlabel: Title of xlabel
:var list xlim: Range for X axis, use [1.132 ,None] to define only one limit
:var list ylim: Range for Y axis, use [1.132 ,None] to define only one limit
:var int shot_num: Number of shot, default is autodetect from current folder
:var double tvec_rescale: Rescale X axis to ms !! Be careful if the axis is not time in seconds !!
:var double data_rescale: Rescale Y axis
:var bool reduction: Try to reduce number of plotted points, automatically used for integrated signal
:var bool integrated: Integrate the input signal
:var list columns: Plot only the selected columns
:var str yscale: scale of y-axis, use `log`
:var bool plot_limits: If there is plasma, plot lines at start and end of the plasma
:var double smoothing: Smooth the signal
:var str line_format: Standard formating string
:var bool xgrid: Plot automatic xgrid in the plot
:var bool ygrid: Plot automatic ygrid in the plot
:var array xerr: error in X direction
:var array yerr: error in Y direction
:var **kwargs : keyword arguments to be passed to the standard plot() command, optional
this method accepts **kwargs as a standard :func:`matplotlib.pyplot.plot` function
Advanced:
:var list vlines : list of positions of vertical lines in plot (see historical analysis)
:var list annotate: list of dicts => plt.annotate(annot['text'], annot['pos'], rotation=annot['angle'] )
"""
try:
# process input setting
if line_format != "":
kwargs['fmt'] = line_format
cData = Shot(shot_num)
is_valid = False
#return "srs"
if type(diagn) is str: # if diagn
is_valid = cData._is_valid_name(diagn)
if is_valid == 0: # not in database
[tvec,data] = load_adv(diagn) # prefere locally saved data
elif is_valid in (1,2): # data
[tvec,data] = cData[diagn]
elif is_valid == 3: # DAS, channel_name
[tvec,data] = cData["any", diagn]
else:
raise NotImplementedError("ambiguous diagn name")
elif (type(diagn) is list or type(diagn) is tuple ) and len(diagn) == 2:
[tvec,data] = deepcopy(diagn)
elif type(diagn) is Data:
[tvec, data] = [diagn.tvec, diagn]
else:
raise NotImplementedError, "Unsupported input"
# Errobars ===================
if type(data) is Data or not 'yerr' in kwargs or 'xerr' in kwargs:
# 'xerr' not user defined in the function call
if 'xerr' not in kwargs : kwargs['xerr'] = data.tvec_err if type(data) is Data else None
if 'yerr' not in kwargs : kwargs['yerr'] = data.data_err if type(data) is Data else None
for name in ['xerr', 'yerr']:
err = kwargs[name]
if not any(isnone(err)) and ndim(err)==0: err = float(err)
if not any(isnone(err)) and not isscalar(err):
err = array(squeeze(err), ndmin=2)
if size(err,0) == 1: err = err.T
if size(err,1) == 1:
err = err.repeat(2, axis=1)
kwargs[name] = err
else:
kwargs['xerr'] = None; kwargs['yerr'] = None;
# =================================
if smoothing > 0:
_, data, _ = DiffFilter(data, mean(diff(tvec)), 2000, 1e6)
if len(columns) > 0:
data = data[:,columns]
data *= data_rescale
tvec *= tvec_rescale
if not any(isnone(kwargs['yerr'])): kwargs['yerr'] *= data_rescale
if not any(isnone(kwargs['xerr'])): kwargs['xerr'] *= tvec_rescale
Ncol = size(data,1) if ndim(data) == 2 else 1
data = reshape(data, (-1,Ncol)) # workaround for 1D arrays
if integrated:
dt = mean(diff(tvec))
data = cumsum(data,0)*dt
if reduction or integrated or smoothing > 10:
tvec, data = random_number_downsample(tvec, data, 4000) # downsample data to 4000 points
# default values
plasma = False
start = -inf
start_pl = nan
end_pl = nan
end = inf
try:
plasma = cData['plasma']
except:
pass
if plasma and plot_limits and tvec_rescale == 1e3: # tvec_rescale == 1e3 => x axis is probably time !!!
start_pl = cData['plasma_start']*1e3
end_pl = cData['plasma_end']*1e3
if cData['ucd'] > 10:
start = min(start_pl - 2, cData['tcd']*1e3) # 2ms before plasma start
else:
start = start_pl - 2
end = end_pl + 2
elif tvec_rescale == 1e3:
try:
start = cData['tb']*1e3-1
end = 35
except:
pass
t_range = amax(tvec) - amin(tvec)
if len(xlim) == 0:
xlim = array([max(start, amin(tvec)-t_range*0.05), min(end, amax(tvec)+t_range*0.05)])
if xlim[0] > xlim[1]: # some errror (no intersection)
xlim = array([amin(tvec)-t_range*0.05, amax(tvec)+t_range*0.05])
if type(data) is Data and ylabel is None:
if not all(isnone(data.ax_labels)):
ylabel = data.ax_labels[1]
if type(data) is Data and diagn_name is None:
diagn_name = data.info
output = list()
columns = range(Ncol)
if len(ylim) == 0:
ylim = [None, None]
### ============ auto-labeling =====================
for i in columns:
column_name = None
if len(diagn_name) == len(columns) and type(diagn_name) is list : # every column has its own label
column_name = diagn_name[i]
else:
if is_valid == 2: # use predefined name in config if it is DAS
column_name = cData.channel_name(diagn, i)
if column_name is None: # nothing matched ...
if len(columns) == 1: # only one variable
column_name = diagn_name
elif len(diagn_name) > 0: # any name
column_name = diagn_name+ ' '+str(i)
else:
column_name = ""
if "label" in kwargs: kwargs.pop( "label" ) ## BUG FIXME !!!!
kwargs_tmp = dict(label=fix_str(column_name), **kwargs)
o = {'tvec': tvec , 'data':data[:,i], 'ylabel':fix_str( ylabel ),
'xlabel': fix_str(xlabel), 'ylim':ylim, 'xlim':xlim, 'xstart':start_pl, 'xend':end_pl, 'yscale':yscale,
'vlines': vlines, 'annotate': annotate , 'xgrid':xgrid, 'ygrid':ygrid,
'kwargs': kwargs_tmp }
output.append(deepcopy(o))
if len(output) == 1:
return output[0]
else:
return output
except Exception, e:
print "Some mistake during get_data: " + str(e)
#return
#raise
#pass
def fix_str(string):
""" FIx the string to avoid any problems during plotting """
if isnone(string): return ""
if not re.match(".*\$.*", string): # hot fix of "latex font" if you can use $ solve it yourselfi
string = re.sub('([^\\\])_', r'\1\_', string) # basic latex friendly fix_string
return string
###################### WEB PAGES ##############################
def get_page_paths(shot, page, default_path = 'shots'):
"""
Return basic paths need for web pages generation
"""
base_path = default_path + '/' + str(shot)+'/'
[ page_path, page] = os.path.split(page)
page_path = base_path + page_path
return page_path,base_path, page
def emph(text):
return '<tt class="file docutils literal"><span class="pre">'+str(text)+'</span></tt>'
def modified(file):
try:
return datetime.datetime.fromtimestamp(os.stat(file)[-1]).strftime('%Y-%m-%d')
except:
return 'N/A'
def make_image(img_path, name = ""):
""" Show image from `img_path` with title `name`. Use loading animation of not avalible. Use svgz if possible and if the browser is capable to load.
"""
img_path = re.sub("(.+)(\.[\w]{3,4})$",r"\1", img_path) # remove file ending (.png)
out = ""
if name is not "":
out += "<h4>"+name+"</h4>\n"
#full_path = './'+img_path
full_path = img_path
rand_end = "?%s"%random.randint(1e6)
out += """<?php
$u_agent = $_SERVER['HTTP_USER_AGENT'];
if (file_exists("%s.png"))
{
// in firefox use svgz images !!
if (file_exists("%s.svgz") && (preg_match('/Firefox/i',$u_agent) || preg_match('/Chrome/i',$u_agent) ))
{
echo "<img src='%s.svgz%s' alt='%s' align='middle'/><br/>";
}
else
{
echo "<img src='%s.png%s' alt='%s' align='middle'/><br/>";
}
}
else
{
echo "<img src='/_static/loading.gif'/><br/>";
echo "<!-- Missing file %s.png --> " ;
}
?>""" % (full_path,full_path,full_path,rand_end,name,full_path,rand_end,name,full_path )
return out
def make_zoom_image(img_path, name = "", group = ""):
""" Load image from `img_path` with title 'name`. Make image zoomable, if possible use svgz version.
"""
rand_end = "?%s"%random.randint(1e6)
img_path = re.sub("(.+)(\.[\w]{3,4})$",r"\1", img_path) # remove file ending (.png)
out = ""
if name is not "":
out += "<h4>"+name+"</h4>\n"
out += """<?php
$u_agent = $_SERVER['HTTP_USER_AGENT'];
if (file_exists('""" + img_path + ".svgz" + """') && (preg_match('/Firefox/i',$u_agent) || preg_match('/Chrome/i',$u_agent) )){
$path='"""+ img_path +".svgz"+rand_end+"""'; }
else {
$path='"""+ img_path + ".png"+rand_end+"""'; }
?> """ + \
"<a href='<?php echo $path ?>' rel='lightbox"""+group+"'><img src='"+ img_path + ".png"+rand_end + "' alt='"+name + ("' title='"+name if name != "" else "") + "'/></a> "
return out
def make_config(path):
""" Try to nicely load and show config from text file in shape
value = 1
value2 = 3
"""
try:
cfg = cat(path, [], True)
out ="""
<table class="docutils field-list" frame="void" rules="none">
<col class="field-name" />
<col class="field-body" />
<tbody valign="top">
"""
for line in cfg:
# replace : or = by new cell
out += ' <tr class="field"><th class="field-name"> ' + re.sub('(.+)(\:|=)(.+)', r'\1 </th><td class="field-body"> \3', line) + ' </td> </tr>'
out += "</tbody></table>"
except:
out = "<h4> Missing config file " + path + "</h4>"
return out
def wiki(link):
"""
Short cat to geenrate link to wiki
"""
if re.match("^\/", link): # link to internal wiki, starts with /
return "<a href='http://golem.fjfi.cvut.cz:5001"+link+"'><font size=1><b><sup>WIKI</sup></b></font></a>"
else:
return "<a href='"+link+"'><font size=1><b><sup>WIKI</sup></b></font></a>"
def source_link(path , revision,default_path, link, name):
"""
Link to formatted link in bitbucket.org
"""
path = re.sub("^"+default_path+"\/[\d]+\/", "", path)
path = re.sub("\/\/", "\/", path + '/') # remove shots/ + shot_number from path
return '<a href="https://bitbucket.org/michalodstrcil/golem_velin/src/'+ revision +'/includes/'+link+ '">[' +name+ ']</a>'
################################# Historical analysis #############################################################
def get_history(diagn, shots, verbose=True, dtype = "float" ):
"""
Rapidly loads and returns values of scalar variables from selected shots
"""
# load default_path => either vshots or shots
N = len(shots)
Data = Shot(shots[-1])
shot = Data.shot_num
if os.path.exists(shotdir + '/' + str(shot) + '/' + diagn):
path = shotdir + '/' + str(shot) + '/' + diagn
else:
path = Data.get_data(diagn, return_path=True)
data = zeros(N, dtype=dtype)
try:
data *= NaN
except:
pass
for i in range(N):
path_tmp = re.sub(str(shot), str(shots[i]), path)
if os.path.exists(path_tmp) and os.path.getsize(path_tmp) > 0:
try:
f = open( path_tmp , 'r').read()
#if not re.match( ".*[A-Za-z\+].*", f): # no letters in the string
#data[i] = float( f )
#elif re.match('OK', f): # fix for plasma status
#data[i] = 1
try:
data[i] = float( f )
except:
data[i] = f
except Exception, e:
print diagn, "ERROR :", str(e), 'shot', shots[i]
if verbose:
sys.stdout.write("\r %2.1f %%" %(i/double(N)*100))
sys.stdout.flush()
sys.stdout.write("\r")
return data
############################# Advanced algorithms, do not touch :) #################################################
# Author: Tomas Odstrcil
# Date: 2012
def DiffFilter(signal,dt, win,lam):
from scipy.signal import medfilt, fftconvolve #loadování zabere většinu času
# signal - data; dt - sampling time of the signal, lam - regularization parameter, win - incresed length due to edge effets
signal = array(signal.T, copy = False,ndmin = 2).T
signal = medfilt(signal,(11,1))
#extend the signal due to edge effects
extended_signal = zeros((size(signal,0)+win,size(signal,1)))
extended_signal[win/2:-(win)/2,:] = signal
extended_signal[:win/2,:] = median(signal[ :win/4,:],axis=0)
extended_signal[-win/2-1:,:] = median(signal[-win/4:,:],axis=0)
N = size(extended_signal,0)
#calculate (II + lam*DD)^{-1}*I^T ale ve fourierovské doméně!
csignal = -cumsum(extended_signal, axis = 0,out = extended_signal)
fsig = fft.rfft(csignal,axis = 0).T
DDfft = (4*sin(linspace(0, pi/2, N/2+1))**2)
IIfft = 1/(DDfft+DDfft[1]*1e-4)
A = IIfft + lam*DDfft
fsig /= A
f = fft.irfft(fsig.T,axis = 0)
#remove edge effects
f = f[win/2-1:-(win+1)/2,:]
retrofit = cumsum(f, axis = 0)
chi2 = sum((retrofit-signal)**2, axis = 0)/size(signal,0)
return f/dt, retrofit, chi2
def DiffFilter_old(signal,dt, max_win= 2000, regularization = 1e6):
from scipy.signal import medfilt, fftconvolve
Nsig = size(signal,1) if ndim(signal) == 2 else 1
signal = reshape(signal, (-1, Nsig))
signal = medfilt(signal,(11, 1))
lam = regularization*(dt/5e-6)**-4
win = 15*(24*lam)**0.25
if win > max_win:
print 'Warning: too small window'
win = min(win, max_win) #to je otypovaná závislost,
win = int(win)+4
Integ = tril( ones((win, win)), 0)
II = dot(Integ.T,Integ)
D = diag(ones(win))-diag(ones(win-1),1)
DD = dot(D.T,D)
b = zeros(win)
b[win/2] = 1
B = dot(Integ.T,b)
A = II + lam*DD
c = array(linalg.solve( A,B.T))
extended_signal = zeros((size(signal,0)+win-1, Nsig))#ones(len(signal)+win-1)
extended_signal[win/2:-(win)/2+1,:] = signal
extended_signal[:win/2,:] += median(signal[:win/4],axis=0)
extended_signal[-win/2:,:] += median(signal[-win/4:],axis=0)
diffsig = empty(shape(signal))
for i in range(Nsig):
diffsig[:,i] = fftconvolve(extended_signal[:,i] , c, mode = 'valid')/dt
retrofit = cumsum(diffsig, axis = 0)*dt
retrofit += mean(signal[:int(4e-3/dt),:], axis = 0) #FIXME
#plot(diffsig)
##figure()
#plot(c)
#show()
#plot(signal)
#plot(retrofit)
#plot(signal-retrofit)
#show()
chi2 = sum((retrofit-signal)**2, axis = 0)/len(signal)
return diffsig, squeeze(retrofit), chi2
def GapsFilling(signal,win = 100, lam = 1e-1):
"""
=============================== 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+1)*win
ext_signal = zeros(n_ext)
ext_signal[(n_ext-n)/2:-(n_ext-n)/2] = signal
ext_signal[:(n_ext-n)/2+1] = median(signal[~isnan(signal)][:win/2])
ext_signal[-(n_ext-n)/2-1:] = median(signal[~isnan(signal)][win/2:])
intervals = arange(0,n+win/2, 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]
chi2 = sum((ext_signal-recon)[~ind_nan]**2)/n
recon = recon[(n_ext-n)/2:-(n_ext-n)/2]
#x = arange(n)
#plot(x,recon)
#plot(x,signal, '.')
#show()
return recon, chi2
def deconvolution( signal,responseFun, win, regularization):
from scipy.signal import fftconvolve
from scipy.sparse import spdiags
from numpy.matlib import repmat
"""
Hlavní algoritmu provádějící dekonvoluci
vrací to dekonvoluovaný signál a zpětnou rekonstrukci
je to zasložené na jednoduché metodě maximalizující hladkost
okno win musí být větší než "support" funkce g
"""
win = 2*(win/2)+1
responseFun/=sum(responseFun)
# vytvoří se konvoluční matice - matice odezvy systému na delta funkci
i_middle = argmax(responseFun)
diags = arange(i_middle-len(responseFun),+ i_middle) #zjistit kam to posouvá
diag_data = repmat(responseFun[::-1], win,1).T
ConvMatrix = spdiags(diag_data, diags, win, win).todense()
D = diag(ones(win))-diag(ones(win-1),1)
DD = dot(D.T,D)
#příprava dat do hlavného algoritmu na dekonvoluci, řeší se soustava Tg = f
T = matrix(ConvMatrix)
TT = dot(T.T,T)
f = zeros(win)
f[win/2] = 1
responseFun2 = zeros(2*len(responseFun))
responseFun2[len(responseFun):] = responseFun
responseFun = responseFun2
g = linalg.solve(TT +regularization*DD,dot(T.T,f).T )
g = squeeze(array(g))
extended_signal = zeros(len(signal)+win-1)
extended_signal[win/2:-win/2+1] = signal
extended_signal[:win/2] += median(signal[:win/4],axis = 0)
extended_signal[-win/2+1:] += median(signal[-win/4:],axis = 0)
deconv_sig = fftconvolve(extended_signal ,g, mode = 'valid')
retrofit = fftconvolve(deconv_sig ,responseFun, mode = 'same')
#plot(g)
#plot(responseFun)
#show()
#plot(deconv_sig)
#plot(retrofit, '--')
#plot(signal)
#show()
return deconv_sig, retrofit
#Ich = load('Ich.npy')
#tvec = Ich[:,0]
#signal = Ich[:,1]
#dt = tvec[1]-tvec[0]
#t_exp = 1e-4
#responseFun = exp(-tvec[:int(6*t_exp/dt)]/t_exp)
#win = 100
#deconvolution( signal,responseFun, win, 1e-2)
#deconvolve / convolve with exponential kernel
def deconvolveExp( signal,t_exp,dt,win,regularization):
from scipy.signal import fftconvolve
tvec = arange(len(signal))*dt
responseFun = exp(-tvec[:int(6*abs(t_exp)/dt)]/abs(t_exp))
responseFun/=sum(responseFun)
if abs(t_exp) <= 2*dt:
return signal, signal
elif t_exp < 0:
responseFun2 = zeros(2*len(responseFun))
responseFun2[len(responseFun):] = responseFun
responseFun = responseFun2
conv_sig = fftconvolve(signal,responseFun, mode = 'full')
nc = len(conv_sig)
n = len(signal)
conv_sig = conv_sig[(nc-n)/2 :-(nc-n)/2]
return conv_sig, signal
else:
return deconvolution( signal,responseFun, win, regularization)
