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

Uploaded

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