diff src/mz_6_rescan_mz2.sci @ 0:245b2c2b3d75 draft

Uploaded
author jcb-mpl
date Tue, 27 Apr 2021 14:13:47 +0000
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/mz_6_rescan_mz2.sci	Tue Apr 27 14:13:47 2021 +0000
@@ -0,0 +1,102 @@
+function x_tr_mz2=mz_6_rescan_mz2(x_tr_mz,thresh)
+    
+    // 12juin19: cette fonction marche très bien!!
+    // mais on essaie de l'améliorer avec mz_6_rescan3
+    
+    // cette fonction réalise des boucles sur la matrice EIC, en regroupant les EIC qui se rapprochent
+    // stop quand sur 3 boucles successives on n'a pas regroupé plus de 5 EICs chaque fois 
+    // stop aussi quand on ne regroupe aucune EIC dans une boucle 
+
+    // thresh: des valeurs testées, pour chercher la valeur optimale 
+    // ex: thresh=0.001
+    
+    // normalement, cette fonction tourne rapidement, ~5 minutes 30jan20
+        
+    if argn(2)<3 then
+        coeff_thresh=1;
+    end;
+    
+    flag=0;
+
+    // verif qu'on a bien un div
+    if size(x_tr_mz.d,1)~=size(x_tr_mz.i,1) | size(x_tr_mz.d,2)~=size(x_tr_mz.v,1) then 
+        error('not a  div, mz_6_rescan_mz2')
+    end;
+
+
+    n_tr=size(x_tr_mz.i,1);
+    
+    // verif. que les m/z de x_tr_mz sont ordonnes 
+    x_tr_mz_v=strtod(x_tr_mz.v); 
+    diff_i=x_tr_mz_v(2:$)-x_tr_mz_v(1:$-1); // verif. que les m/z sont en ordre croissant
+    index=find(diff_i<0);
+    if max(size(index))>0 then              // pas ordonné -> il faut le faire 
+        [nul,tri0]=gsort(x_tr_mz_v,'g','i');
+        x_tr_mz=x_tr_mz(:,tri0);
+        x_tr_mz_v=x_tr_mz_v(tri0);
+    end;
+
+
+     while flag<2;    // enleve inferieur 31jan20
+        
+
+        clear x_tr_mz2;
+
+        flag=2*flag;
+        
+        x_tr_mz_v=strtod(x_tr_mz.v);                    
+        x_sum=sum(x_tr_mz.d,'r')';
+        n1=max(size(x_tr_mz.d,2));        
+
+//disp('index2')
+
+        for i=2:n1;
+            if x_tr_mz_v(i)-x_tr_mz_v(i-1)<thresh then                                                    // avant le 27mai 17h30: <=1.01
+                x_tr_mz_v(i-1)=x_tr_mz_v(i-1:i)'*x_sum(i-1:i) / (x_sum(i-1:i)'*ones(2,1));     // on repositionne la 1° valeur de m/z
+                x_tr_mz_v(i)=0;                                                                // et on annule la 2° valeur de m/z
+                x_tr_mz.d(:,i-1)=x_tr_mz.d(:,i-1)+x_tr_mz.d(:,i);                              // somme des signaux
+            end; 
+        end;
+        
+        x_tr_mz.v=string(x_tr_mz_v);           // on remet les nlles valeurs de m/z dans x_tr_mz 
+
+
+//disp('index3')
+
+        // tri
+        tri=find(x_tr_mz_v~=0);
+
+        // ligne suivante gardée: car elle ne prend que qq secondes avec x_tr_mz=3000x140000 30jan20
+        x_tr_mz2=x_tr_mz(:,tri);
+
+  
+//disp('index4')
+
+   
+        // verification des dimensions div 
+        if size(x_tr_mz2.d,1)~=size(x_tr_mz2.i,1) | size(x_tr_mz2.d,2)~=size(x_tr_mz2.v,1) then
+             disp('error in dimensions, mz_6_rescan_mz2')
+        end;
+
+//disp('index5')
+
+
+        n2=size(x_tr_mz2.d,2);
+        
+        x_tr_mz=x_tr_mz2;
+        
+        if flag <= 2 then                                                       // la boucle peut continuer
+            if abs(n1-n2)<=5 then
+                if n1==n2 then   
+                    flag=2;                                                    // stop immédiat -> sortie de la fonction
+                else
+                    flag=1;                                                    // pour bloquer à la boucle suivante
+                end
+            end
+        end
+       
+    end
+
+    
+    
+endfunction