function tracks = track(xyzs,maxdisp,param)

%;
% ; see http://glinda.lrsm.upenn.edu/~weeks/idl
% ;   for more information
% ;
% ;+
% ; NAME:
% ; track
% ; PURPOSE:
% ; Constructs n-dimensional trajectories from a scrambled list of
% ; particle coordinates determined at discrete times (e.g. in
% ; consecutive video frames).
% ; CATEGORY:
% ; Image Processing
% ; CALLING SEQUENCE:
% ; result = track( positionlist, maxdisp, param )
% ;  set all keywords in the space below
% ; INPUTS:
% ; positionlist: an array listing the scrambled coordinates and data 
% ;     of the different particles at different times, such that:
% ;  positionlist(0:d-1,*): contains the d coordinates and
% ;     data for all the particles, at the different times. must be positve
% ;  positionlist(d,*): contains the time t that the position 
% ;     was determined, must be integers (e.g. frame number.  These values must 
% ;               be monotonically increasing and uniformly gridded in time.
% ; maxdisp: an estimate of the maximum distance that a particle 
% ;     would move in a single time interval.(see Restrictions)
%  OPTIONAL INPUT:
%   param:  a structure containing a few tracking parameters that are
%       needed for many applications.  If param is not included in the
%       function call, then default values are used.  If you set one value
%       make sure you set them all:
% ;         param.mem: this is the number of time steps that a particle can be
% ;             'lost' and then recovered again.  If the particle reappears
% ;             after this number of frames has elapsed, it will be
% ;             tracked as a new particle. The default setting is zero.
% ;             this is useful if particles occasionally 'drop out' of
% ;             the data.
% ;         param.dim: if the user would like to unscramble non-coordinate data
% ;             for the particles (e.g. apparent radius of gyration for
% ;             the particle images), then positionlist should
% ;             contain the position data in positionlist(0:param.dim-1,*)
% ;             and the extra data in positionlist(param.dim:d-1,*). It is then
% ;             necessary to set dim equal to the dimensionality of the
% ;             coordinate data to so that the track knows to ignore the
% ;             non-coordinate data in the construction of the 
% ;             trajectories. The default value is two.
% ;         param.good: set this keyword to eliminate all trajectories with
% ;             fewer than param.good valid positions.  This is useful
% ;             for eliminating very short, mostly 'lost' trajectories
% ;             due to blinking 'noise' particles in the data stream.
%;          param.quiet: set this keyword to 1 if you don't want any text
% ; OUTPUTS:
% ; result:  a list containing the original data rows sorted 
% ;     into a series of trajectories.  To the original input 
% ;     data structure there is appended an additional column 
% ;     containing a unique 'id number' for each identified 
% ;     particle trajectory.  The result array is sorted so 
% ;     rows with corresponding id numbers are in contiguous 
% ;     blocks, with the time variable a monotonically
% ;     increasing function inside each block.  For example:
% ;     
% ;     For the input data structure (positionlist):
% ;         (x)      (y)      (t)
% ;     pos = 3.60000      5.00000      0.00000
% ;           15.1000      22.6000      0.00000
% ;           4.10000      5.50000      1.00000 
% ;           15.9000      20.7000      2.00000
% ;           6.20000      4.30000      2.00000
% ;
% ;     IDL> res = track(pos,5,mem=2)
% ;
% ;     track will return the result 'res'
% ;         (x)      (y)      (t)          (id)
% ;     res = 3.60000      5.00000      0.00000      0.00000
% ;           4.10000      5.50000      1.00000      0.00000
% ;           6.20000      4.30000      2.00000      0.00000
% ;           15.1000      22.6000      0.00000      1.00000
% ;           15.9000      20.7000      2.00000      1.00000
% ;
% ;     NB: for t=1 in the example above, one particle temporarily
% ;     vanished.  As a result, the trajectory id=1 has one time
% ;     missing, i.e. particle loss can cause time gaps to occur 
% ;     in the corresponding trajectory list. In contrast:
% ;
% ;     IDL> res = track(pos,5)
% ;
% ;     track will return the result 'res'
% ;         (x)      (y)      (t)          (id)
% ;     res = 15.1000      22.6000      0.00000      0.00000
% ;                   3.60000      5.00000      0.00000      1.00000
% ;               4.10000      5.50000      1.00000      1.00000
% ;               6.20000      4.30000      2.00000      1.00000
% ;               15.9000      20.7000      2.00000      2.00000
% ; 
% ;     where the reappeared 'particle' will be labelled as new
% ;     rather than as a continuation of an old particle since
% ;     mem=0.  It is up to the user to decide what setting of 
% ;     'mem' will yeild the highest fidelity .
% ; 
% ; SIDE EFFECTS:
% ; Produces informational messages.  Can be memory intensive for
% ; extremely large data sets.
% ; RESTRICTIONS:
% ; maxdisp should be set to a value somewhat less than the mean 
% ; spacing between the particles. As maxdisp approaches the mean
% ; spacing the runtime will increase significantly. The function 
% ; will produce an error message: "Excessive Combinatorics!" if
% ; the run time would be too long, and the user should respond 
% ; by re-executing the function with a smaller value of maxdisp.
% ; Obviously, if the particles being tracked are frequently moving
% ; as much as their mean separation in a single time step, this
% ; function will not return acceptable trajectories.
% ; PROCEDURE:
% ; Given the positions for n particles at time t(i), and m possible
% ; new positions at time t(i+1), this function considers all possible 
% ; identifications of the n old positions with the m new positions,
% ; and chooses that identification which results in the minimal total
% ; squared displacement. Those identifications which don't associate
% ; a new position within maxdisp of an old position ( particle loss )
% ; penalize the total squared displacement by maxdisp^2. For non-
% ; interacting Brownian particles with the same diffusivity, this
% ; algorithm will produce the most probable set of identifications 
% ; ( provided maxdisp >> RMS displacement between frames ).
% ; In practice it works reasonably well for systems with oscillatory,
% ; ballistic, correlated and random hopping motion, so long as single 
% ; time step displacements are reasonably small.  NB: multidimensional
% ; functionality is intended to facilitate tracking when additional
% ; information regarding target identity is available (e.g. size or 
% ; color).  At present, this information should be rescaled by the
% ; user to have a comparable or smaller (measurement) variance than 
% ; the spatial displacements.
% ;
% ; MODIFICATION HISTORY:
% ;  2/93 Written by John C. Crocker, University of Chicago (JFI).
% ;  7/93 JCC fixed bug causing particle loss and improved performance
% ;     for large numbers of (>100) particles.
% ; 11/93 JCC improved speed and memory performance for large
% ;     numbers of (>1000) particles (added subnetwork code).
% ;  3/94 JCC optimized run time for trivial bonds and d<7. (Added
% ;     d-dimensional raster metric code.)
% ;  8/94 JCC added functionality to unscramble non-position data
% ;     along with position data.
% ;  9/94 JCC rewrote subnetwork code and wrote new, more efficient 
% ;     permutation code.
% ;  5/95 JCC debugged subnetwork and excessive combinatorics code.
% ; 12/95 JCC added memory keyword, and enabled the tracking of
% ;     newly appeared particles.
% ;  3/96 JCC made inipos a keyword, and disabled the adding of 'new'
% ;     particles when inipos was set.
% ;  3/97 JCC added 'add' keyword, since Chicago users didn't like 
% ;     having particle addition be the default. 
% ;  9/97 JCC added 'goodenough' keyword to improve memory efficiency
% ;     when using the 'add' keyword and to filter out bad tracks.
% ;       10/97 JCC streamlined data structure to speed runtime for >200 
% ;               timesteps.  Changed 'quiet' keyword to 'verbose'. Made
% ;               time labelling more flexible (uniform and sorted is ok).
% ;  9/98 JCC switched trajectory data structure to a 'list' form,
% ;     resolving memory issue for large, noisy datasets.
% ;  2/99 JCC added Eric Weeks's 'uberize' code to post-facto 
% ;     rationalize the particle id numbers, removed 'add' keyword.
% ;  1/05 Transmuted to MATLAB by D. Blair
% ;  5/05  ERD Added the param structure to simplify calling.
%    6/05  ERD Added quiet to param structure
%    7/05  DLB Fixed slight bug in trivial bond code
%    3/07  DLB Fixed bug with max disp pointed out by Helene Delanoe-Ayari
%
% ; This code 'track.pro' is copyright 1999, by John C. Crocker. 
% ; It should be considered 'freeware'- and may be distributed freely 
% ; (outside of the military-industrial complex) in its original form 
% ; when properly attributed.
% ;
% ;-

dd = length(xyzs(1,:));

%use default parameters if none given
if nargin==2
    %default values
    memory_b=0; % if mem is not needed set to zero
    goodenough = 0;  % if goodenough is not wanted set to zero
    dim = dd - 1;
    quiet=0;
else
    memory_b    =   param.mem;
    goodenough  =   param.good;
    dim         =   param.dim;
    quiet       =   param.quiet;
end


% checking the input time vector
t = xyzs(:,dd);
st = circshift(t,1);
st = t(2:end) - st(2:end);
if  sum(st(find(st < 0))) ~= 0
    disp('The time vectors is not in order')
    return
end
info = 1;

w = find(st > 0);
z = length(w);
z = z +1;
if isempty(w)
    disp('All positions are at the same time... go back!')
    return
end

% partitioning the data with unique times

%res = unq(t);
% implanting unq directly
    indices = find(t ~= circshift(t,-1));
        count = length(indices);
        if count > 0
            res = indices;
        else  
            res = length(t) -1;
        end
 %%%%%%%%%%%%%%%%%%%%%%%       
        
res = [1,res',length(t)];
ngood = res(2) - res(1) + 1;
eyes = 1:ngood;
pos = xyzs(eyes,1:dim);
istart = 2;
n = ngood;

zspan = 50;
if n > 200 
    zspan = 20;
end
if n > 500 
    zspan = 10;
end
resx = zeros(zspan,n) - 1;

bigresx = zeros(z,n) - 1;
mem = zeros(n,1);
%  whos resx
%  whos bigresx
uniqid = 1:n;
maxid = n;
olist = [0.,0.];

if goodenough > 0 
    dumphash = zeros(n,1);
    nvalid = ones(n,1);
end

%  whos eyes;
resx(1,:) = eyes;
% setting up constants
maxdisq = maxdisp^2;

% John calls this the setup for "fancy code" ???
notnsqrd = (sqrt(n*ngood) > 200) & (dim < 7);
notnsqrd = notnsqrd(1);

if notnsqrd
    %;   construct the vertices of a 3x3x3... d-dimensional hypercube
    
    cube = zeros(3^dim,dim);
    
    
    for d=0:dim-1,
        numb = 0;
        for j=0:(3^d):(3^dim)-1,
            cube(j+1:j+(3^(d)),d+1) = numb;
            numb = mod(numb+1,3);
        end
    end    
    
    %   calculate a blocksize which may be greater than maxdisp, but which
    %   keeps nblocks reasonably small.  
    
    volume = 1;
    for d = 0:dim-1
        minn = min(xyzs(w,d+1));
        maxx = max(xyzs(w,d+1));
        volume = volume * (maxx-minn);
    end
    volume;
    blocksize = max( [maxdisp,((volume)/(20*ngood))^(1.0/dim)] );
end
%   Start the main loop over the frames.
for i=istart:z
    ispan = mod(i-1,zspan)+1;
    %disp(ispan)
    % get new particle positions
    m = res(i+1) - res(i);
    res(i);
    eyes = 1:m;
    eyes = eyes + res(i);
    
    if m > 0
        
        xyi = xyzs(eyes,1:dim);
        found = zeros(m,1);
        
        % THE TRIVIAL BOND CODE BEGINS   
        
        if notnsqrd
            %Use the raster metric code to do trivial bonds
            
            % construct "s", a one dimensional parameterization of the space 
            % which consists of the d-dimensional raster scan of the volume.)
            
            abi = fix(xyi./blocksize);
            abpos = fix(pos./blocksize);
            si = zeros(m,1);
            spos = zeros(n,1);
            dimm = zeros(dim,1);
            coff = 1.;
            
            for j=1:dim
                minn = min([abi(:,j);abpos(:,j)]);
                maxx = max([abi(:,j);abpos(:,j)]);
                abi(:,j) = abi(:,j) - minn;
                abpos(:,j) = abpos(:,j) - minn;
                dimm(j) = maxx-minn + 1;
                si = si + abi(:,j).*coff;
                spos = spos + abpos(:,j).*coff;
                coff = dimm(j).*coff;
            end
            nblocks = coff;
            % trim down (intersect) the hypercube if its too big to fit in the
            % particle volume. (i.e. if dimm(j) lt 3)
            
            cub = cube;
            deg = find( dimm < 3);
            if ~isempty(deg)
                for j = 0:length(deg)-1
                    cub = cub(find(cub(:,deg(j+1)) < dimm(deg(j+1))),:);
                end
            end 
            
            % calculate the "s" coordinates of hypercube (with a corner @ the origin)
            scube = zeros(length(cub(:,1)),1);
            coff = 1;
            for j=1:dim
                scube = scube + cub(:,j).*coff;
                coff = coff*dimm(j);      
            end
            
            % shift the hypercube "s" coordinates to be centered around the origin
            
            coff = 1;
            for j=1:dim
                if dimm(j) > 3
                    scube = scube - coff;
                end
                coff = dimm(j).* coff;
            end
            scube = mod((scube + nblocks),nblocks);
            % get the sorting for the particles by their "s" positions.
            [ed,isort] = sort(si);
            
            % make a hash table which will allow us to know which new particles
            % are at a given si.
            strt = zeros(nblocks,1) -1;
            fnsh = zeros(nblocks,1);
            h = find(si == 0);
            lh = length(h);
            if lh > 0
                
            si(h) = 1;  
            end
            
            for j=1:m
                if strt(si(isort(j))) == -1
                    strt(si(isort(j))) = j;
                    fnsh(si(isort(j))) = j;
                else
                    fnsh(si(isort(j))) = j;
                end
            end
            if lh > 0
            si(h) = 0;   
            end
            coltot = zeros(m,1);
            rowtot = zeros(n,1);
            which1 = zeros(n,1);
            for j=1:n
                
                
                map = fix(-1);
                
                scub_spos = scube + spos(j);
                s = mod(scub_spos,nblocks);
                whzero = find(s == 0 );
                if ~isempty(whzero)
                    nfk = find(s ~=0);
                    s = s(nfk);
                end
               
                w = find(strt(s) ~= -1);
               
                ngood = length(w);
                ltmax=0;
                if ngood ~= 0
 
                    s = s(w);
                    for k=1:ngood
                        map = [map;isort( strt(s(k)):fnsh(s(k)))];
                    end
                    map = map(2:end);
%                     if length(map) == 2
%                         if (map(1) - map(2)) == 0
%                             map = unique(map);
%                          end
%                     end
                    %   map = map(umap);
                    %end
                    % find those trival bonds
                    distq = zeros(length(map),1);
                    for d=1:dim     
                        distq = distq + (xyi(map,d) - pos(j,d)).^2;
                    end
                    ltmax = distq < maxdisq;
                    
                    rowtot(j) = sum(ltmax);
                    
                    if rowtot(j) >= 1 
                        w = find(ltmax == 1);
                        coltot( map(w) ) = coltot( map(w)) +1;
                        which1(j) = map( w(1) );
                    end
                end

            end
            
          
            ntrk = fix(n - sum(rowtot == 0));
        
            w = find( rowtot == 1);
            ngood = length(w);

           
            if ngood ~= 0 
                ww = find(coltot( which1(w) ) == 1);
                ngood = length(ww);
                if ngood ~= 0 
                     %disp(size(w(ww)))
                    resx(ispan,w(ww)) = eyes( which1(w(ww)));
                    found( which1( w(ww))) = 1;
                    rowtot( w(ww)) = 0;
                    coltot( which1(w(ww))) = 0;
                end
            end
            
            labely = find( rowtot > 0);
            ngood = length(labely);
            if ngood ~= 0 
                labelx = find( coltot > 0);
                
                nontrivial = 1;
            else
                nontrivial = 0;
            end

        else 
    
            %   or: Use simple N^2 time routine to calculate trivial bonds      
    
            % let's try a nice, loopless way!
            % don't bother tracking perm. lost guys.
            wh = find( pos(:,1) >= 0);
            ntrack = length(wh);
            if ntrack == 0 
                'There are no valid particles to track idiot!'
                break
            end
            xmat = zeros(ntrack,m);
            count = 0;
            for kk=1:ntrack
                for ll=1:m
                    xmat(kk,ll) = count;
                    count = count+1;
                end
            end
            count = 0;
            for kk=1:m
                for ll=1:ntrack
                    ymat(kk,ll) = count;
                    count = count+1;
                end
            end

            xmat = (mod(xmat,m) + 1);
            ymat = (mod(ymat,ntrack) +1)';
            [lenxn,lenxm] = size(xmat);
%            whos ymat
%            whos xmat
%            disp(m)

            for d=1:dim
                x = xyi(:,d);
                y = pos(wh,d);
                xm = x(xmat);
                ym = y(ymat(1:lenxn,1:lenxm));
                if size(xm) ~= size(ym)
                    xm = xm';
                end
                
                if d == 1
                    dq = (xm -ym).^2;
                    %dq = (x(xmat)-y(ymat(1:lenxn,1:lenxm))).^2;
                else
                    dq = dq + (xm-ym).^2;
                    %dq = dq + (x(xmat)-y(ymat(1:lenxn,1:lenxm)) ).^2;
                end
            end
            
            ltmax = dq < maxdisq;
            
            % figure out which trivial bonds go with which
            
            rowtot = zeros(n,1);
            rowtot(wh) = sum(ltmax,2);
            
            
            if ntrack > 1 
                coltot = sum(ltmax,1);
            else
                coltot = ltmax;
            end
            which1 = zeros(n,1);
            for j=1:ntrack 
                [mx, w] = max(ltmax(j,:));
                which1(wh(j)) = w;
            end
            
            ntrk = fix( n - sum(rowtot == 0));
            w= find( rowtot == 1) ;
            ngood = length(w);
            if ngood ~= 0
                ww = find(coltot(which1(w)) == 1);
                ngood = length(ww);
                if ngood ~= 0 
                    resx( ispan, w(ww) ) = eyes( which1( w(ww)));
                    found(which1( w(ww))) = 1;
                    rowtot(w(ww)) = 0;
                    coltot(which1(w(ww))) = 0;
                end
            end
            
            labely = find( rowtot > 0);
            ngood = length(labely);
            
            if ngood ~= 0
                labelx = find( coltot > 0);
                nontrivial = 1;
            else
                nontrivial = 0;
            end
        end
        
        %THE TRIVIAL BOND CODE ENDS
        
        if nontrivial
            
            xdim = length(labelx);
            ydim = length(labely);
            
            %  make a list of the non-trivial bonds            
            
            bonds = zeros(1,2);
            bondlen = 0;
            
            for j=1:ydim
                distq = zeros(xdim,1);
                
                for d=1:dim
                    %distq
                    distq = distq + (xyi(labelx,d) - pos(labely(j),d)).^2;
                    %distq    
                end
                
                w= find(distq <  maxdisq)' - 1;
                ngood = length(w);
                newb = [w;(zeros(1,ngood)+j)];
                
                
                bonds = [bonds;newb'];
               
                bondlen = [ bondlen;distq( w + 1) ];
                
            end
            bonds = bonds(2:end,:);
            
            bondlen = bondlen(2:end);
            numbonds = length(bonds(:,1));
            mbonds = bonds;
            max([xdim,ydim]);
                
                
            if max([xdim,ydim]) < 4
                nclust = 1;
                maxsz = 0;
                mxsz = xdim;
                mysz = ydim;
                bmap = zeros(length(bonds(:,1)+1),1) - 1;
               
            else
           

                %   THE SUBNETWORK CODE BEGINS            
                lista = zeros(numbonds,1);
                listb = zeros(numbonds,1);
                nclust = 0;
                maxsz = 0;
                thru = xdim;
                
                while thru ~= 0
                    %  the following code extracts connected 
                    %   sub-networks of the non-trivial 
                    %   bonds.  NB: lista/b can have redundant entries due to 
                    %   multiple-connected subnetworks      
                    
                    
                    w = find(bonds(:,2) >= 0);
   %                 size(w)
                    
                    lista(1) = bonds(w(1),2);
                    listb(1) = bonds(w(1),1);
                    bonds(w(1),:) = -(nclust+1);
                    bonds;
                    adda = 1;
                    addb = 1;
                    donea = 0;
                    doneb = 0;
                    if (donea ~= adda) | (doneb ~= addb)
                        true = 0;
                    else
                    true = 1;   
                    end
                                        
                    while ~true
                        
                        if (donea ~= adda)
                            w = find(bonds(:,2) == lista(donea+1));
                            ngood = length(w);
                            if ngood ~= 0 
                                listb(addb+1:addb+ngood,1) = bonds(w,1);
                                bonds(w,:) = -(nclust+1);
                                addb = addb+ngood;
                            end
                            donea = donea+1;
                        end
                        if (doneb ~= addb) 
                            w = find(bonds(:,1) == listb(doneb+1));
                            ngood = length(w);
                            if ngood ~= 0
                                lista(adda+1:adda+ngood,1) = bonds(w,2);
                                bonds(w,:) = -(nclust+1);
                                adda = adda+ngood;
                            end
                            doneb = doneb+1;
                        end
                      if (donea ~= adda) | (doneb ~= addb) 
                          true = 0;
                      else  
                          true = 1;
                      end
                    end
                    
                    [pp,pqx] = sort(listb(1:doneb));
                    %unx =  unq(listb(1:doneb),pqx);
                    %implanting unq directly
                        arr = listb(1:doneb);
                        q = arr(pqx);
                        indices = find(q ~= circshift(q,-1));
                        count = length(indices);
                        if count > 0
                            unx = pqx(indices);
                        else
                            unx = length(q) -1;
                        end
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    
                    xsz = length(unx);
                    [pp,pqy] = sort(lista(1:donea));
                    %uny =  unq(lista(1:donea),pqy);
                    %implanting unq directly
                        arr = lista(1:donea);
                        q = arr(pqy);
                        indices = find(q ~= circshift(q,-1));
                        count = length(indices);
                        if count > 0
                            uny = pqy(indices);
                        else
                            uny = length(q) -1;
                        end
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                   
                    
                    
                    
                    
                    ysz = length(uny);
                    if xsz*ysz > maxsz
                        maxsz = xsz*ysz;
                        mxsz = xsz;
                        mysz = ysz; 
                    end
                    
                    
                    thru = thru -xsz;
                    nclust = nclust + 1;
                end
                bmap = bonds(:,2);                    
            end
            % THE SUBNETWORK CODE ENDS
            % put verbose in for Jaci
            
            %   THE PERMUTATION CODE BEGINS
            
            for nc =1:nclust
                w = find( bmap == -1*(nc));
                
                nbonds = length(w);
                bonds = mbonds(w,:);
                lensq = bondlen(w);
                [pq,st] = sort( bonds(:,1));
                %un = unq(bonds(:,1),st);
                   %implanting unq directly     
                        arr = bonds(:,1);
                        q = arr(st);
                        indices = find(q ~= circshift(q,-1));
                        count = length(indices);
                        if count > 0
                            un = st(indices);
                        else
                            un = length(q) -1;
                        end
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
                
                uold = bonds(un,1);
                
                nold = length(uold);
                
                %un = unq(bonds(:,2));
                
                %implanting unq directly  
                indices = find(bonds(:,2) ~= circshift(bonds(:,2),-1));
                count = length(indices);
                    if count > 0
                        un = indices;
                    else  
                        un = length(bonds(:,2)) -1;
                    end
                 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
                    
                unew = bonds(un,2);
                nnew = length(unew);
                
                if nnew > 5
                    rnsteps = 1;
                    for ii =1:nnew
                        rnsteps = rnsteps * length( find(bonds(:,2) == ...
                            unew(ii)));
                        if rnsteps > 5.e+4
                            disp('Warning: difficult combinatorics encountered.')
                        end
                        if rnsteps > 2.e+5
                            disp(['Excessive Combinitorics you FOOL LOOK WHAT YOU HAVE' ...
                                    ' DONE TO ME!!!'])
                            return
                        end
                    end
                end
                st = zeros(nnew,1);
                fi = zeros(nnew,1);
                h = zeros(nbonds,1);
                ok = ones(nold,1);
                nlost = (nnew - nold) > 0;
                
                
                for ii=1:nold 
                    h(find(bonds(:,1) == uold(ii))) = ii;
                end
                st(1) = 1 ;
                fi(nnew) = nbonds; % check this later
                if nnew > 1 
                    sb = bonds(:,2);
                    sbr = circshift(sb,1);
                    sbl = circshift(sb,-1);
                    st(2:end) = find( sb(2:end) ~= sbr(2:end)) + 1;
                    fi(1:nnew-1) = find( sb(1:nbonds-1) ~= sbl(1:nbonds-1));
                end
%                if i-1 == 13
%                    hi
%                end
                checkflag = 0;
                while checkflag ~= 2
                    
                    pt = st -1;
                    lost = zeros(nnew,1);
                    who = 0;
                    losttot = 0;
                    mndisq = nnew*maxdisq;
                    
                    
                    while who ~= -1
                    
                        if pt(who+1) ~= fi(who+1)
                            
                            
                            w = find( ok( h( pt( who+1 )+1:fi( who+1 ) ) ) ); % check this -1
                            ngood = length(w);
                            if ngood > 0
                                if pt(who+1) ~= st(who+1)-1
                                    ok(h(pt(who+1))) = 1;
                                end
                                pt(who+1) = pt(who+1) + w(1);
                                ok(h(pt(who+1))) = 0;
                                if who == nnew -1
                                    ww = find( lost == 0);
                                    dsq = sum(lensq(pt(ww))) + losttot*maxdisq;
                                    
                                    if dsq < mndisq 
                                        minbonds = pt(ww);
                                        mndisq = dsq;
                                    end
                                else
                                    who = who+1;
                                end
                            else
                                if ~lost(who+1) & (losttot ~= nlost)
                                    lost(who+1) = 1;
                                    losttot = losttot + 1;
                                    if pt(who+1) ~= st(who+1) -1;
                                        ok(h(pt(who+1))) = 1;
                                    end
                                    if who == nnew-1
                                        ww = find( lost == 0);
                                        dsq = sum(lensq(pt(ww))) + losttot*maxdisq;
                                        if dsq < mndisq
                                            minbonds = pt(ww);
                                            mndisq = dsq;
                                        end
                                    else    
                                       who = who + 1;
                                    end
                                   
                                else
                                    if pt(who+1) ~= (st(who+1) -1) 
                                        ok(h(pt(who+1))) = 1;
                                    end
                                    pt(who+1) = st(who+1) -1;
                                    if lost(who+1) 
                                        lost(who+1) = 0;
                                        losttot = losttot -1;
                                    end
                                    who = who -1;
                                end
                            end
                        else  
                            if ~lost(who+1) & (losttot ~= nlost)
                                lost(who+1) = 1;
                                losttot = losttot + 1;
                                if pt(who+1) ~= st(who+1)-1
                                    ok(h(pt(who+1))) = 1;
                                end
                                if who == nnew -1
                                    ww = find( lost == 0);
                                    dsq = sum(lensq(pt(ww))) + losttot*maxdisq;
                                    
                                    if dsq < mndisq
                                        minbonds = pt(ww);
                                        mndisq = dsq;
                                    end
                                else
                                    who = who + 1;
                                end
                            else
                                if pt(who+1) ~= st(who+1) -1
                                    ok(h(pt(who+1))) = 1;
                                end
                                pt(who+1) = st(who+1) -1;
                                if lost(who+1) 
                                    lost(who+1) = 0;
                                    losttot = losttot -1;
                                end
                                who = who -1;
                            end
                        end
                    end
                    
                    checkflag = checkflag + 1;
                    if checkflag == 1
                        plost = min([fix(mndisq/maxdisq) , (nnew -1)]);
                        if plost > nlost 
                            nlost = plost; 
                        else
                            checkflag = 2;
                        end
                    end
                    
                end  
                %   update resx using the minimum bond configuration               
                
                resx(ispan,labely(bonds(minbonds,2))) = eyes(labelx(bonds(minbonds,1)+1));
                found(labelx(bonds(minbonds,1)+1)) = 1;

            end

            %   THE PERMUTATION CODE ENDS
        end
        
        w = find(resx(ispan,:) >= 0);
        nww = length(w);
        
        if nww > 0 
            pos(w,:) = xyzs( resx(ispan,w) , 1:dim);
            if goodenough > 0 
                nvalid(w) = nvalid(w) + 1;
            end
        end  %go back and add goodenough keyword thing   
        newguys = find(found == 0);
        
       nnew = length(newguys);
      
        if (nnew > 0) % & another keyword to workout inipos
            newarr = zeros(zspan,nnew) -1;
            resx = [resx,newarr];

            resx(ispan,n+1:end) = eyes(newguys);
            pos = [[pos];[xyzs(eyes(newguys),1:dim)]];
            nmem = zeros(nnew,1);
            mem = [mem;nmem];
            nun = 1:nnew;
            uniqid = [uniqid,((nun) + maxid)];
            maxid = maxid + nnew;
            if goodenough > 0 
                dumphash = [dumphash;zeros(1,nnew)'];
                nvalid = [nvalid;zeros(1,nnew)'+1];
            end
            % put in goodenough 
            n = n + nnew;
                
        end
 
    else
        ' Warning- No positions found for t='
    end
    w = find( resx(ispan,:) ~= -1);
    nok = length(w);
    if nok ~= 0
        mem(w) =0;
    end
    
    mem = mem + (resx(ispan,:)' == -1);
    wlost = find(mem == memory_b+1);
    nlost =length(wlost);

    if nlost > 0 
        pos(wlost,:) = -maxdisp;
        if goodenough > 0
            wdump = find(nvalid(wlost) < goodenough);
            ndump = length(wdump);
            if ndump > 0
                dumphash(wlost(wdump)) = 1;
            end
        end
        % put in goodenough keyword stuff if 
    end
    if (ispan == zspan) | (i == z)
        nold = length(bigresx(1,:));
        nnew = n-nold;
        if nnew > 0
            newarr = zeros(z,nnew) -1;
            bigresx = [bigresx,newarr];
        end
        if goodenough > 0  
            if (sum(dumphash)) > 0
                wkeep = find(dumphash == 0);
                nkeep = length(wkeep);
                resx = resx(:,wkeep);
                bigresx = bigresx(:,wkeep);
                pos = pos(wkeep,:);
                mem = mem(wkeep);
                uniqid = uniqid(wkeep);
                nvalid = nvalid(wkeep);
                n = nkeep;
                dumphash = zeros(nkeep,1);
            end
        end
        
        % again goodenough keyword
        if quiet~=1
            disp(strcat(num2str(i), ' of ' ,num2str(z), ' done.  Tracking  ',num2str(ntrk),' particles  ', num2str(n),' tracks total'));
        end
        bigresx(i-(ispan)+1:i,:) = resx(1:ispan,:);
        resx = zeros(zspan,n) - 1;

 
        wpull = find(pos(:,1) == -maxdisp);
        npull = length(wpull);
        
        if npull > 0
            lillist = zeros(1,2);
            for ipull=1:npull
                wpull2 = find(bigresx(:,wpull(ipull)) ~= -1);
                npull2 = length(wpull2);
                thing = [bigresx(wpull2,wpull(ipull)),zeros(npull2,1)+uniqid(wpull(ipull))];
                lillist = [lillist;thing];
                
            end
            olist = [[olist];[lillist(2:end,:)]];
             
        end

        
 
        wkeep = find(pos(:,1) >= 0);
        nkeep = length(wkeep);
        if nkeep == 0 
                'Were going to crash now, no particles....'
        end
        resx = resx(:,wkeep);
        bigresx = bigresx(:,wkeep);
        pos = pos(wkeep,:);
        mem = mem(wkeep);
        uniqid = uniqid(wkeep);
        n = nkeep;
        dumphash = zeros(nkeep,1);
        if goodenough > 0
            nvalid = nvalid(wkeep);
        end
    end
   
end

if goodenough > 0 
    nvalid = sum(bigresx >= 0 ,1);
    wkeep = find(nvalid >= goodenough);
    nkeep = length(wkeep);
    if nkeep == 0
        for i=1:10
        disp('You are not going any further, check your params and data')
        end
        disp('the code broke at line 1045')
        return
    end
    if nkeep < n
        bigresx = bigresx(:,wkeep);
        n = nkeep;
        uniqid = uniqid(wkeep);
        pos = pos(wkeep,:);
    end
end


wpull = find( pos(:,1) ~= -2*maxdisp);
npull = length(wpull);
if npull > 0
    lillist = zeros(1,2);
    for ipull=1:npull
        wpull2 = find(bigresx(:,wpull(ipull)) ~= -1);
        npull2 = length(wpull2);   
        thing = [bigresx(wpull2,wpull(ipull)),zeros(npull2,1)+uniqid(wpull(ipull))];
        lillist = [lillist;thing];
    end
    olist = [olist;lillist(2:end,:)];
end

olist = olist(2:end,:);
%bigresx = 0;
%resx = 0;

nolist = length(olist(:,1));
res = zeros(nolist,dd+1);
for j=1:dd
    res(:,j) = xyzs(olist(:,1),j);
end
res(:,dd+1) = olist(:,2);

% this is uberize included for simplicity of a single monolithic code

ndat=length(res(1,:));
newtracks=res;
    

%u=unq(newtracks(:,ndat));

% inserting unq
indices = find(newtracks(:,ndat) ~= circshift(newtracks(:,ndat),-1));
        count = length(indices);
        if count > 0
            u = indices;
        else  
            u = length(newtracks(:,ndat)) -1;
        end


ntracks=length(u);
u=[0;u];
for i=2:ntracks+1
    newtracks(u(i-1)+1:u(i),ndat) = i-1;
end

% end of uberize code

tracks = newtracks;



 
