hooke/plugin/pcluster.py

   1 #!/usr/bin/env python
   2 # -*- coding: utf-8 -*-
   3 from ..libhooke import WX_GOOD, ClickedPoint, config_file_path
   4
   5 from mdp import pca
   6 import wxversion
   7 wxversion.select(WX_GOOD)
   8 from wx import PostEvent
   9 import numpy as np
  10 import scipy as sp
  11 import copy
  12 import os.path
  13 import time
  14 import pylab as pyl
  15
  16 from .. import curve as lhc
  17
  18
  19 import warnings
  20 warnings.simplefilter('ignore',np.RankWarning)
  21
  22
  23 class pclusterCommands(object):
  24
  25     def _plug_init(self):
  26         self.clustplot1=None
  27         self.clustplot2=None
  28
  29     def do_pcluster(self,args):
  30         '''
  31         pCLUSTER
  32         (pcluster.py)
  33         Automatically measures peaks and extracts informations for further clustering
  34         (c)Paolo Pancaldi, Massimo Sandal 2009
  35         '''
  36         blindw = str(self.convfilt_config['blindwindow'])
  37         pclus_dir = "pCluster_blind"+blindw+"_"+time.strftime("%Y%m%d_%H%M")
  38         self.my_work_dir = os.path.join(os.getcwd(), pclus_dir)
  39         self.my_curr_dir = os.path.basename(os.getcwd())
  40         os.mkdir(self.my_work_dir)
  41
  42         #--Custom persistent length
  43         pl_value=None
  44         for arg in args.split():
  45             #look for a persistent length argument.
  46             if 'pl=' in arg:
  47                 pl_expression=arg.split('=')
  48                 pl_value=float(pl_expression[1]) #actual value
  49             else:
  50                 pl_value=None
  51
  52         #configuration variables
  53         min_npks = self.convfilt_config['minpeaks']
  54         min_deviation = self.convfilt_config['mindeviation']
  55
  56         pclust_filename = "automeasure_"+self.my_curr_dir+"_blind"+blindw+".txt" #raw_input('Automeasure filename? ')
  57         realclust_filename = "coordinate_"+self.my_curr_dir+"_blind"+blindw+".txt" #raw_input('Coordinates filename? ')
  58         peackforce_filename = "peakforce_"+self.my_curr_dir+"_blind"+blindw+".txt"  #raw_input('Peacks and Forces filename? ')
  59
  60         f=open(self.my_work_dir+pclust_filename,'w+')
  61         f.write('Analysis started '+time.asctime()+'\n')
  62         f.write('----------------------------------------\n')
  63         f.write('; Contour length (nm)  ;  Persistence length (nm) ;  Max.Force (pN)  ;  Slope (N/m) ;  Sigma contour (nm) ; Sigma persistence (nm)\n')
  64         f.close()
  65
  66         f=open(self.my_work_dir+realclust_filename,'w+')
  67         f.write('Analysis started '+time.asctime()+'\n')
  68         f.write('----------------------------------------\n')
  69         f.write('; Peak number ; Mean delta (nm)  ;  Median delta (nm) ;  Mean force (pN)  ;  Median force (pN) ; First peak length (nm) ; Last peak length (nm) ; Max force (pN) ; Min force (pN) ; Max delta (nm) ; Min delta (nm) ; Peaks Diff\n')
  70         f.close()
  71
  72         f=open(self.my_work_dir+peackforce_filename,'w+')
  73         f.write('Analysis started '+time.asctime()+'\n')
  74         f.write('----------------------------------------\n')
  75         f.write('; Peak number  ;  1 peak Length (nm) ; 1 peak Force (pN) ;  2 peak Length (nm) ; 2 peak Force (pN) ;  3 peak Length (nm) ; 3 peak Force (pN) ;  4 peak Length (nm) ; 4 peak Force (pN) ;  5 peak Length (nm) ; 5 peak Force (pN) ;  6 peak Length (nm) ; 6 peak Force (pN) ;  7 peak Length (nm) ; 7 peak Force (pN) ;  8 peak Length (nm) ; 8 peak Force (pN)\n')
  76         f.close()
  77
  78
  79         def plot_informations(itplot,pl_value):
  80             '''
  81             OUR VARIABLES
  82             contact_point.absolute_coords               (2.4584142802103689e-007, -6.9647135616234017e-009)
  83             peak_point.absolute_coords                  (3.6047748250571423e-008, -7.7142802788854212e-009)
  84             other_fit_point.absolute_coords     (4.1666139243838867e-008, -7.3759393477579707e-009)
  85             peak_location                                                                               [510, 610, 703, 810, 915, 1103]
  86             peak_size                                                                                           [-1.2729111505202212e-009, -9.1632775347399312e-010, -8.1707438353929907e-010, -8.0335812578148904e-010, -8.7483955226387558e-010, -3.6269619757067322e-009]
  87             params                                                                                                      [2.2433999931959462e-007, 3.3230248825175678e-010]
  88             fit_errors                                                                                  [6.5817195369767644e-010, 2.4415923138871498e-011]
  89             '''
  90             fit_points=int(self.config['auto_fit_points']) # number of points to fit before the peak maximum <50>
  91
  92             T=self.config['temperature'] #temperature of the system in kelvins. By default it is 293 K. <301.0>
  93             cindex=self.find_contact_point(itplot[0]) #Automatically find contact point <158, libhooke.ClickedPoint>
  94             contact_point=self._clickize(itplot[0].vectors[1][0], itplot[0].vectors[1][1], cindex)
  95             self.basepoints=[]
  96             base_index_0=peak_location[-1]+self.fit_interval_nm(peak_location[-1], itplot[0], self.config['auto_right_baseline'],False)
  97             self.basepoints.append(self._clickize(itplot[0].vectors[1][0],itplot[0].vectors[1][1],base_index_0))
  98             base_index_1=self.basepoints[0].index+self.fit_interval_nm(self.basepoints[0].index, itplot[0], self.config['auto_left_baseline'],False)
  99             self.basepoints.append(self._clickize(itplot[0].vectors[1][0],itplot[0].vectors[1][1],base_index_1))
 100             self.basecurrent=self.current.path
 101             boundaries=[self.basepoints[0].index, self.basepoints[1].index]
 102             boundaries.sort()
 103             to_average=itplot[0].vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
 104             avg=np.mean(to_average)
 105             return fit_points, contact_point, pl_value, T, cindex, avg
 106
 107         def features_peaks(itplot, peak, fit_points, contact_point, pl_value, T, cindex, avg):
 108             '''
 109             calculate informations for each peak and add they in
 110             c_lengths, p_lengths, sigma_c_lengths, sigma_p_lengths, forces, slopes
 111             '''
 112             c_leng=None
 113             p_leng=None
 114             sigma_c_leng=None
 115             sigma_p_leng=None
 116             force=None
 117             slope=None
 118
 119             delta_force=10
 120             slope_span=int(self.config['auto_slope_span'])
 121
 122             peak_point=self._clickize(itplot[0].vectors[1][0],itplot[0].vectors[1][1],peak)
 123             other_fit_point=self._clickize(itplot[0].vectors[1][0],itplot[0].vectors[1][1],peak-fit_points)
 124
 125             points=[contact_point, peak_point, other_fit_point]
 126
 127             params, yfit, xfit, fit_errors = self.wlc_fit(points, itplot[0].vectors[1][0], itplot[0].vectors[1][1], pl_value, T, return_errors=True)
 128
 129             #Measure forces
 130             delta_to_measure=itplot[0].vectors[1][1][peak-delta_force:peak+delta_force]
 131             y=min(delta_to_measure)
 132             #Measure slopes
 133             slope=self.linefit_between(peak-slope_span,peak)[0]
 134             #check fitted data and, if right, add peak to the measurement
 135             if len(params)==1: #if we did choose 1-value fit
 136                 p_leng=pl_value
 137                 c_leng=params[0]*(1.0e+9)
 138                 sigma_p_leng=0
 139                 sigma_c_leng=fit_errors[0]*(1.0e+9)
 140                 force = abs(y-avg)*(1.0e+12)
 141             else: #2-value fit
 142                 p_leng=params[1]*(1.0e+9)
 143                 #check if persistent length makes sense. otherwise, discard peak.
 144                 if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
 145                     '''
 146                     p_lengths.append(p_leng)
 147                     c_lengths.append(params[0]*(1.0e+9))
 148                     sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
 149                     sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
 150                     forces.append(abs(y-avg)*(1.0e+12))
 151                     slopes.append(slope)
 152                     '''
 153                     c_leng=params[0]*(1.0e+9)
 154                     sigma_c_leng=fit_errors[0]*(1.0e+9)
 155                     sigma_p_leng=fit_errors[1]*(1.0e+9)
 156                     force=abs(y-avg)*(1.0e+12)
 157                 else:
 158                     p_leng=None
 159                     slope=None
 160             #return c_lengths, p_lengths, sigma_c_lengths, sigma_p_lengths, forces, slopes
 161             return  c_leng, p_leng, sigma_c_leng, sigma_p_leng, force, slope
 162
 163
 164         # ------ PROGRAM -------
 165         c=0
 166         for item in self.current_list:
 167             c+=1
 168             item.identify(self.drivers)
 169             itplot=item.curve.default_plots()
 170             flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
 171             itplot[0]=flatten(itplot[0], item, customvalue=1)
 172             try:
 173                 peak_location,peak_size=self.exec_has_peaks(item,min_deviation)
 174             except:
 175                 #We have troubles with exec_has_peaks (bad curve, whatever).
 176                 #Print info and go to next cycle.
 177                 print 'Cannot process ',item.path
 178                 continue
 179
 180             if len(peak_location)==0:
 181                 print 'No peaks!'
 182                 continue
 183
 184             fit_points, contact_point, pl_value, T, cindex, avg = plot_informations(itplot,pl_value)
 185
 186             print '\n\nCurve',item.path, 'is',c,'of',len(self.current_list),': found '+str(len(peak_location))+' peaks.'
 187
 188             #initialize output data vectors
 189             c_lengths=[]
 190             p_lengths=[]
 191             sigma_c_lengths=[]
 192             sigma_p_lengths=[]
 193             forces=[]
 194             slopes=[]
 195
 196             #loop each peak of my curve
 197             for peak in peak_location:
 198                 c_leng, p_leng, sigma_c_leng, sigma_p_leng, force, slope = features_peaks(itplot, peak, fit_points, contact_point, pl_value, T, cindex, avg)
 199                 for var, vector in zip([c_leng, p_leng, sigma_c_leng, sigma_p_leng, force, slope],[c_lengths, p_lengths, sigma_c_lengths, sigma_p_lengths, forces, slopes]):
 200                     if var is not None:
 201                         vector.append(var)
 202
 203             #FIXME: We need a dictionary here...
 204             allvects=[c_lengths, p_lengths, sigma_c_lengths, sigma_p_lengths, forces, slopes]
 205             for vect in allvects:
 206                 if len(vect)==0:
 207                     for i in range(len(c_lengths)):
 208                         vect.append(0)
 209
 210             print 'Measurements for all peaks detected:'
 211             print 'contour (nm)', c_lengths
 212             print 'sigma contour (nm)',sigma_c_lengths
 213             print 'p (nm)',p_lengths
 214             print 'sigma p (nm)',sigma_p_lengths
 215             print 'forces (pN)',forces
 216             print 'slopes (N/m)',slopes
 217
 218             '''
 219             write automeasure text file
 220             '''
 221             print 'Saving automatic measurement...'
 222             f=open(self.my_work_dir+pclust_filename,'a+')
 223             f.write(item.path+'\n')
 224             for i in range(len(c_lengths)):
 225                 f.write(' ; '+str(c_lengths[i])+' ; '+str(p_lengths[i])+' ; '+str(forces[i])+' ; '+str(slopes[i])+' ; '+str(sigma_c_lengths[i])+' ; '+str(sigma_p_lengths[i])+'\n')
 226             f.close()
 227
 228             peak_number=len(c_lengths)
 229
 230             '''
 231             write peackforce text file
 232             '''
 233             print 'Saving automatic measurement...'
 234             f=open(self.my_work_dir+peackforce_filename,'a+')
 235             f.write(item.path+'\n')
 236             peackforce_info = ''
 237             for i in range(len(c_lengths)):
 238                 peackforce_info = peackforce_info + ' ; ' + str(c_lengths[i]) + ' ; ' + str(forces[i])
 239             f.write(' ; '+str(peak_number)+peackforce_info+'\n')
 240             f.close()
 241
 242             '''
 243             calculate clustering coordinates
 244             '''
 245             if peak_number > 1:
 246                 deltas=[]
 247                 for i in range(len(c_lengths)-1):
 248                     deltas.append(c_lengths[i+1]-c_lengths[i])
 249
 250                 delta_mean=np.mean(deltas)
 251                 delta_median=np.median(deltas)
 252
 253                 force_mean=np.mean(forces)
 254                 force_median=np.median(forces)
 255
 256                 first_peak_cl=c_lengths[0]
 257                 last_peak_cl=c_lengths[-1]
 258
 259                 max_force=max(forces[:-1])
 260                 min_force=min(forces)
 261
 262                 max_delta=max(deltas)
 263                 min_delta=min(deltas)
 264
 265                 delta_stdev=np.std(deltas)
 266                 forces_stdev=np.std(forces[:-1])
 267
 268                 peaks_diff=(last_peak_cl-first_peak_cl)/peak_number
 269
 270                 print 'Coordinates'
 271                 print 'Peaks',peak_number
 272                 print 'Mean delta',delta_mean
 273                 print 'Median delta',delta_median
 274                 print 'Mean force',force_mean
 275                 print 'Median force',force_median
 276                 print 'First peak',first_peak_cl
 277                 print 'Last peak',last_peak_cl
 278                 print 'Max force',max_force
 279                 print 'Min force',min_force
 280                 print 'Max delta',max_delta
 281                 print 'Min delta',min_delta
 282                 print 'Delta stdev',delta_stdev
 283                 print 'Forces stdev',forces_stdev
 284                 print 'Peaks difference',peaks_diff
 285
 286                 '''
 287                 write clustering coordinates
 288                 '''
 289                 f=open(self.my_work_dir+realclust_filename,'a+')
 290                 f.write(item.path+'\n')
 291                 f.write(' ; '+str(peak_number)+     # non considerato
 292                         ' ; '+str(delta_mean)+      # 0
 293                         ' ; '+str(delta_median)+    # 1 -
 294                         ' ; '+str(force_mean)+      # 2
 295                         ' ; '+str(force_median)+    # 3 -
 296                         ' ; '+str(first_peak_cl)+   # 4 -
 297                         ' ; '+str(last_peak_cl)+    # 5 -
 298                         ' ; '+str(max_force)+       # 6
 299                         ' ; '+str(min_force)+       # 7
 300                         ' ; '+str(max_delta)+       # 8
 301                         ' ; '+str(min_delta)+       # 9
 302                         ' ; '+str(delta_stdev)+     # 10
 303                         ' ; '+str(forces_stdev)+    # 11
 304                         ' ; '+str(peaks_diff)+      # 12
 305                         '\n')
 306                 f.close()
 307
 308         # start PCA
 309         self.do_pca(pclus_dir+"/"+realclust_filename)
 310
 311
 312     def do_pca(self,args):
 313         '''
 314         PCA -> "pca gaeta_coor_blind50.txt 1,3,6"
 315         Automatically measures pca from coordinates filename and shows two interactives plots
 316         With the second argument (arbitrary) you can select the columns and the multiplier factor
 317         to use for the pca (for es "1,3*50,6,8x10,9"). Dont use spaces. "*" or "x" are the same thing.
 318         Without second argument it reads pca_config.txt file
 319         (c)Paolo Pancaldi, Massimo Sandal 2009
 320         '''
 321
 322         # reads the columns of pca
 323         conf=open(config_file_path("pca_config.txt"), 'r')
 324         config = conf.readlines()
 325         conf.close()
 326
 327         self.plot_myCoord = []          # tiene le coordinate prese direttamente dal file creato con pCluster
 328         self.plot_origCoord = []        # tiene le coordinate solo delle colonne scelte e moltiplicate per i valori scelti
 329         self.plot_pcaCoord = []         # tiene le due colonne della pca
 330         self.plot_pcaCoordTr = []       # tiene le due colonne della pca trasposta
 331         self.plot_FiltOrigCoord = []    # tiene le coordinate solo dei punti filtrati per densita
 332         self.plot_FiltPaths = []        # tiene i paths dei plot solo dei punti filtrati per densita
 333         self.plot_paths = []            # tiene i paths di tutti i plots
 334         self.plot_NewPcaCoord = []      # tiene le due colonne della pca filtrate per densita
 335         self.plot_NewPcaCoordTr=[]      # tiene le due colonne della pca trasposta filtrate per densita
 336         plot_path_temp = ""
 337
 338         # prende in inpunt un arg (nome del file)
 339         # e il secondo le colonne su cui lavorare (e' arbitrario, riceve x es "1,2,3")
 340         arg = args.split(" ")
 341         if arg[0]==args:
 342             file_name=args
 343         else:
 344             file_name=arg[0]
 345             config[0] = arg[1]
 346
 347         # creo l'array "plot_myCoord" con tutte le coordinate dei plots
 348         # e l'array plot_paths con tutti i percorsi dei plots
 349         nPlotTot = -3 #tolgo le prime 3 righe iniziali del file
 350         f=open(file_name)
 351         rows = f.readlines()
 352         for row in rows:
 353             if row[0]!=" " and row[0]!="":
 354                 nPlotTot = nPlotTot+1
 355                 plot_path_temp = row
 356             if row[0]==" " and row.find('nan')==-1 and row.find("-1.#IND")==-1:
 357                 row = row[row.index(";",2)+2:].split(" ; ")     # non considero la prima colonna col #picchi
 358                 row = [float(i) for i in row]
 359
 360                 #0:Mean delta, 1:Median delta, 2:Mean force, 3:Median force, 4:First peak length, 5:Last peak length
 361                 #6:Max delta 7:Min delta 8:Max force 9:Min force 10:Std delta 11:Std force
 362                 if (row[0]<500 and row[1]<500 and row[2]<500 and row[3]<500 and row[4]<500 and row[5]<500 and row[6]<500 and row[7]<500 and row[8]<500 and row[9]<500 and row[10]<500 and row[11]<500):
 363                     if (row[0]>0 and row[1]>0 and row[2]>0 and row[3]>0 and row[4]>0 and row[5]>0 and row[6]>0 and row[7]>0 and row[8]>0 and row[9]>0 and row[10]>0 and row[11]>0):
 364                         #row = row[0], row[2], row[3]*3, row[6], row[7]*56, row[8]
 365                         self.plot_myCoord.append(row)
 366                         self.plot_paths.append(plot_path_temp)
 367         f.close()
 368
 369         # creo l'array con alcune colonne e pure moltiplicate
 370         for row in self.plot_myCoord:
 371             res=[]
 372             for cols in config[0].split(","):
 373                 if cols.find("*")!=-1:
 374                     col = int(cols.split("*")[0])
 375                     molt = int(cols.split("*")[1])
 376                 elif cols.find("x")!=-1:
 377                     col = int(cols.split("x")[0])
 378                     molt = int(cols.split("x")[1])
 379                 else:
 380                     col = int(cols)
 381                     molt = 1
 382                 res.append(row[col]*molt)
 383             self.plot_origCoord.append(res)
 384
 385         # array convert, calculate PCA, transpose
 386         self.plot_origCoord = np.array(self.plot_origCoord,dtype='float')
 387         #print self.plot_origCoord.shape
 388         self.plot_pcaCoord = pca(self.plot_origCoord, output_dim=2)     #other way -> y = mdp.nodes.PCANode(output_dim=2)(array)
 389         self.plot_pcaCoordTr = np.transpose(self.plot_pcaCoord)
 390         pca_X=np.array(self.plot_pcaCoordTr[0],dtype='float')
 391         pca_Y=np.array(self.plot_pcaCoordTr[1],dtype='float')
 392
 393         '''
 394         # Start section of testing with good plots                                  # 4 TESTING!
 395         Xsyn_1=[]
 396         Ysyn_1=[]
 397         Xgb1_1=[]
 398         Ygb1_1=[]
 399         Xbad_1=[]
 400         Ybad_1=[]
 401         goodnamefile=open(file_name.replace("coordinate", "good"),'r')
 402         goodnames=goodnamefile.readlines()
 403         nPlotGood = len(goodnames)-2 #tolgo prima e ultima riga
 404         goodnames=[i.split()[0] for i in goodnames[1:]]
 405
 406         for index in range(len(self.plot_paths)):
 407             if self.plot_paths[index][:-1] in goodnames:
 408                 Xsyn_1.append(pca_X[index])
 409                 Ysyn_1.append(pca_Y[index])
 410             else:
 411                 Xbad_1.append(pca_X[index])
 412                 Ybad_1.append(pca_Y[index])
 413         # Stop section of testing with good plots                                   # 4 TESTING!
 414         '''
 415
 416         # print first plot
 417         clustplot1=lhc.PlotObject()
 418         clustplot1.add_set(pca_X,pca_Y)
 419         #clustplot1.add_set(Xbad_1,Ybad_1) # 4 TESTING!
 420         #clustplot1.add_set(Xsyn_1,Ysyn_1) # 4 TESTING!
 421         clustplot1.normalize_vectors()
 422         clustplot1.styles=['scatter', 'scatter','scatter']
 423         clustplot1.colors=[None,'red','green']
 424         clustplot1.destination=0
 425         self._send_plot([clustplot1])
 426         self.clustplot1=clustplot1
 427
 428         # density and filer estimation
 429         kernel = sp.stats.kde.gaussian_kde(sp.c_[pca_X,pca_Y].T)
 430         tallest = 0
 431         for i in range(len(pca_X)):
 432             kern_value = kernel.evaluate([pca_X[i],pca_Y[i]])
 433             if tallest < kern_value:
 434                     tallest = float(kern_value)
 435         if float(config[1]) == 0:
 436             my_filter = float(tallest / 3.242311147)
 437         else:
 438             my_filter = float(config[1])
 439         '''
 440         # section useful only for graphic printing
 441         xmin = pca_X.min()
 442         xmax = pca_X.max()
 443         ymin = pca_Y.min()
 444         ymax = pca_Y.max()
 445         mX, mY = sp.mgrid[xmin:xmax:100j, ymin:ymax:100j]
 446         Z = sp.rot90(sp.fliplr(sp.reshape(kernel(sp.c_[mX.ravel(), mY.ravel()].T).T, mX.T.shape)))
 447         axis_X = np.linspace(xmin,xmax,num=100)
 448         axis_Y = np.linspace(ymin,ymax,num=100)
 449         '''
 450
 451         # density filtering:
 452         # tramite "kernel.evaluate" trovo lo score (altezza) di ogni coordinata e decido se mantenerla o no
 453         filtered_pca_X = []
 454         filtered_pca_Y = []
 455         filtered_PcaCoordTr = []
 456         filtered_PcaCoord = []
 457         for i in range(len(pca_X)):
 458             kern_value = kernel.evaluate([pca_X[i],pca_Y[i]])
 459             if kern_value > my_filter:
 460                 filtered_pca_X.append(pca_X[i])
 461                 filtered_pca_Y.append(pca_Y[i])
 462         filtered_PcaCoordTr.append(filtered_pca_X)
 463         filtered_PcaCoordTr.append(filtered_pca_Y)
 464         filtered_PcaCoord = np.transpose(filtered_PcaCoordTr)
 465
 466         # creo i due array "plot_FiltOrigCoord" e "plot_FiltPaths" contenenti solo i dati filtrati con alta densita
 467         for index in range(len(self.plot_pcaCoord)):
 468             if self.plot_pcaCoord[index] in filtered_PcaCoord:
 469                 self.plot_FiltOrigCoord.append(self.plot_myCoord[index])
 470                 self.plot_FiltPaths.append(self.plot_paths[index])
 471
 472         '''
 473         # START PCA#2: USELESS!!!
 474
 475         # creo l array con alcune colonne e pure moltiplicate
 476         temp_coord = []
 477         for row in self.plot_FiltOrigCoord:
 478             res=[]
 479             for cols in config[2].split(","):
 480                 if cols.find("*")!=-1:
 481                     col = int(cols.split("*")[0])
 482                     molt = int(cols.split("*")[1])
 483                 elif cols.find("x")!=-1:
 484                     col = int(cols.split("x")[0])
 485                     molt = int(cols.split("x")[1])
 486                 else:
 487                     col = int(cols)
 488                     molt = 1
 489                 res.append(row[col]*molt)
 490             temp_coord.append(res)
 491         self.plot_FiltOrigCoord = temp_coord
 492
 493         # ricalcolo la PCA: array convert, calculate PCA, transpose
 494         self.plot_FiltOrigCoord = np.array(self.plot_FiltOrigCoord,dtype='float')
 495         #print self.plot_FiltOrigCoord.shape
 496         self.plot_NewPcaCoord = pca(self.plot_FiltOrigCoord, output_dim=2)      #other way -> y = mdp.nodes.PCANode(output_dim=2)(array)
 497         self.plot_NewPcaCoordTr = np.transpose(self.plot_NewPcaCoord)
 498         pca_X2=np.array(self.plot_NewPcaCoordTr[0],dtype='float')
 499         pca_Y2=np.array(self.plot_NewPcaCoordTr[1],dtype='float')
 500
 501         # Start section of testing with good plots                              # 4 TESTING!
 502         Xsyn_2=[]
 503         Ysyn_2=[]
 504         Xbad_2=[]
 505         Ybad_2=[]
 506         for index in range(len(self.plot_FiltPaths)):
 507             if self.plot_FiltPaths[index][:-1] in goodnames:
 508                 Xsyn_2.append(pca_X2[index])
 509                 Ysyn_2.append(pca_Y2[index])
 510             else:
 511                 Xbad_2.append(pca_X2[index])
 512                 Ybad_2.append(pca_Y2[index])
 513
 514         # print second plot
 515         clustplot2=lhc.PlotObject()
 516         #clustplot2.add_set(pca_X2,pca_Y2)
 517         clustplot2.add_set(Xbad_2,Ybad_2)                                       # 4 TESTING!
 518         clustplot2.add_set(Xsyn_2,Ysyn_2)                                       # 4 TESTING!
 519         clustplot2.normalize_vectors()
 520         clustplot2.styles=['scatter', 'scatter','scatter']
 521         clustplot2.colors=[None,'red','green']
 522         clustplot2.destination=1
 523         self._send_plot([clustplot2])
 524         self.clustplot2=clustplot2
 525         '''
 526
 527         # PRINT density plot
 528         clustplot2=lhc.PlotObject()
 529         clustplot2.add_set(filtered_pca_X,filtered_pca_Y)
 530         clustplot2.normalize_vectors()
 531         clustplot2.styles=['scatter', 'scatter','scatter']
 532         clustplot2.colors=[None,'red','green']
 533         clustplot2.destination=1
 534         self._send_plot([clustplot2])
 535         self.clustplot2=clustplot2
 536
 537         # printing results
 538         config_pca1 = config[0].replace("*", "x").rstrip("\n")
 539         config_pca2 = config[2].replace("*", "x").rstrip("\n")
 540         print ""
 541         print "- START: "+file_name
 542         print "Curve totali: ", nPlotTot
 543         #print "Curve totali good: ", nPlotGood                                  # 4 TESTING!
 544         print "- FILTRO 1: 0-500 e NaN"
 545         print "Curve totali rimaste: ", len(self.plot_origCoord)
 546         #print 'Curve good rimaste: ', len(Xsyn_1)                               # 4 TESTING!
 547         print "- FILTRO 2: PCA:"+config_pca1+" e DENSITA:"+str(my_filter)
 548         print "Curve totali rimaste: ", len(self.plot_FiltOrigCoord)
 549         #print 'Curve good rimaste: ', len(Xsyn_2)                               # 4 TESTING!
 550         print "Piu alta: ", tallest
 551         #print "- FILTRO 3: 2'PCA:"+config_pca2
 552         print ""
 553
 554         # -- exporting coordinates and plot of PCA in debug mode! --
 555         if config[3].find("true")!=-1:
 556             #1' PCA: save plot and build coordinate s file
 557             self.do_export(file_name.replace("coordinate_", "debug_pca1graph_").replace('.txt','_'+config_pca1) + " 0")
 558             f = open(file_name.replace("coordinate_", "debug_pca1coor_").replace('.txt','_'+config_pca1+'.txt'),'w')
 559             for i in range(len(pca_X)):
 560                 f.write (str(i) + "\t" + str(pca_X[i]) + "\t" + str(pca_Y[i]) + "\n")
 561             f.close()
 562             #2' PCA: save plot and build coordinate s file
 563             #self.do_export(file_name.replace("coordinate_", "debug_pca2graph_").replace('.txt','_'+config_pca2) + " 1")
 564             #f = open(file_name.replace("coordinate_", "debug_pca2coor_").replace('.txt','_'+config_pca2+'.txt'),'w')
 565             #for i in range(len(pca_X2)):
 566             #    f.write (str(i) + "\t" + str(pca_X2[i]) + "\t" + str(pca_Y2[i]) + "\n")
 567             #f.close()
 568             #DENSITY: save plot
 569             self.do_export(file_name.replace("coordinate_", "debug_densitygraph_").replace('.txt','_'+config_pca1+'_'+str(my_filter).replace(".",",")) + " 1")
 570             f = open(file_name.replace("coordinate_", "debug_densitycoor_").replace('.txt','_'+config_pca1+'_'+str(my_filter).replace(".",",")+'.txt'),'w')
 571             for i in range(len(filtered_pca_X)):
 572                 f.write (str(i) + "\t" + str(filtered_pca_X[i]) + "\t" + str(filtered_pca_Y[i]) + "\n")
 573             f.close()
 574             #ALL GOOD COORDINATES (without NaN and 0<x<500)
 575             f = open(file_name.replace("coordinate_", "debug_allgoodcoor_"),'w')
 576             for i in range(len(self.plot_myCoord)):
 577                 for cel in self.plot_myCoord[i]:
 578                     f.write (" ; " + str(cel))
 579                 f.write ("\n")
 580             f.close()
 581
 582         # pCLUSTER SAVING!!!
 583         import shutil
 584         pcl_name = file_name.replace("coordinate_", "goodplots_").replace('.txt','_'+config_pca1+'_'+str(my_filter).replace(".",","))
 585         if os.path.exists(pcl_name+slash): shutil.rmtree(pcl_name)
 586         os.mkdir(pcl_name+slash)
 587         f = open(pcl_name+'.txt','w')
 588         for i in range(len(self.plot_FiltPaths)):
 589             myfile = str(self.plot_FiltPaths[i]).rstrip("\n")
 590             f.write (myfile+"\n")
 591             shutil.copy2(myfile, pcl_name)
 592         f.close()
 593
 594
 595     def do_multipca(self,args):
 596         '''
 597         MULTIPCA -> "multipca gaeta_coor_blind50.txt 3"
 598         Automatically multiply the column suggest in second argument for value between 1-100 (step of 2),
 599         measures pca from coordinates filename and save the png plots.
 600         (c)Paolo Pancaldi, Massimo Sandal 2009
 601         '''
 602         # reads the columns of pca
 603         conf=open("pca_config.txt")
 604         config = conf.readlines() # config[0] = "1,2,3"
 605         conf.close()
 606         # cycling pca
 607         arg = args.split(" ")
 608         file_name=arg[0]
 609         column=str(arg[1])
 610         for i in range(1, 51, 1):
 611             self.do_pca(file_name + " " + config[0].replace(column,column+"*"+str(i),1))
 612
 613     def do_doublepca(self,args):
 614         '''
 615         DOUBLEPCA -> "doublepca gaeta_coor_blind50.txt"
 616         Automatically it launches the pca command for all combinations with two column
 617         (c)Paolo Pancaldi, Massimo Sandal 2009
 618         '''
 619         # cycling pca
 620         arg = args.split(" ")
 621         file_name=arg[0]
 622         for i in range(1, 13):
 623             for j in range(1, 13):
 624                 if i!=j:
 625                     self.do_pca(file_name + " " + str(i) + "," + str(j))
 626
 627     def do_triplepca(self,args):
 628         '''
 629         TRIPLEPCA -> "triplepca gaeta_coor_blind50.txt"
 630         Automatically it launches the pca command for all combinations with three column
 631         (c)Paolo Pancaldi, Massimo Sandal 2009
 632         '''
 633         # cycling pca
 634         arg = args.split(" ")
 635         file_name=arg[0]
 636         for i in range(1, 13):
 637             for j in range(1, 13):
 638                 for k in range(1, 13):
 639                     if i!=j and i!=k and j!=k:
 640                         self.do_pca(file_name + " " + str(i) + "," + str(j) + "," + str(k))
 641
 642     def do_pclick(self,args):
 643         '''
 644         It returns id, coordinates and file name of a clicked dot on a PCA graphic
 645         '''
 646
 647         self._send_plot([self.clustplot1]) #quick workaround for BAD problems in the GUI
 648         print 'Click point'
 649         point = self._measure_N_points(N=1, whatset=0)
 650         indice = point[0].index
 651         plot_file = self.plot_paths[indice]
 652         dot_coord = self.plot_pcaCoord[indice]
 653         print "file: " + str(plot_file).rstrip()
 654         print "id: " + str(indice)
 655         print "coord: " + str(dot_coord)
 656         self.do_genlist(str(plot_file))
 657         #self.do_jump(str(plot_file))
 658
 659     # indea iniziata e messa da parte...
 660     def do_peakforce(self, args):
 661         '''
 662         peackforce -> "peackforce peackforce_file.txt"
 663         Automatically measures peack and force plots
 664         (c)Paolo Pancaldi, Massimo Sandal 2009
 665         '''
 666
 667         # prende in inpunt un arg (nome del file)
 668         file_name=args
 669         f=open(file_name)
 670
 671         # scrivo un file temp
 672         g = open('_prove.txt','w')
 673
 674         plot_file = '';
 675         rows = f.readlines()
 676         for row in rows:
 677             if row[0]=="/":
 678                 plot_file = row
 679             if row[0]==" " and row.find('nan')==-1 and row.find("-1.#IND")==-1:
 680                 # FILTRO SUI 7 PICCHI
 681                 num_pic = int(row.split(" ; ")[1])
 682                 if num_pic==7:
 683                     width_force = row.split(" ; ")
 684                     w1 = float(width_force[2]); f1 = float(width_force[3]);
 685                     w2 = float(width_force[4]); f2 = float(width_force[5]);
 686                     w3 = float(width_force[6]); f3 = float(width_force[7]);
 687                     w4 = float(width_force[8]); f4 = float(width_force[9]);
 688                     w5 = float(width_force[10]); f5 = float(width_force[11]);
 689                     w6 = float(width_force[12]); f6 = float(width_force[13]);
 690                     w7 = float(width_force[14]); f7 = float(width_force[15]);
 691                     if w1>0 and w1<1000 and w2>0 and w2<1000 and w3>0 and w3<1000 and w4>0 and w4<1000 and w5>0 and w5<1000 and w6>0 and w6<1000 and w7>0 and w7<1000:
 692                         score_76 = abs(32 - (w7 - w6))
 693                         score_65 = abs(32 - (w6 - w5))
 694                         score_54 = abs(32 - (w5 - w4))
 695                         score_43 = abs(32 - (w4 - w3))
 696                         score_32 = abs(32 - (w3 - w2))
 697                         score_21 = abs(32 - (w2 - w1))
 698                         writeme = str(score_76) + " --- " + str(row)
 699                         g.write(writeme)
 700         g.close()
 701         f.close()