hooke/plugin/pcluster.py

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