util1_jan/read_array_tiles_sp.m

   1 function [varargout]=read_array_tiles_sp(root,dmtiles,mptiles,num1,num2);
   2 % [i,j,a]=read_array_tiles_sp(root,dmtiles,mptiles,num1,num2);
   3 % a=read_array_tiles_sp(root,dmtiles,mptiles,num1,num2);
   4 %
   5 % read array produced by matching calls write_array_m in module_fr_sfire_util.F
   6 % cycle tiles over dmtiles and mptiles
   7 % from files root[_num1[_num2]].nnnn.txt,  nnnn=1:dmntiles
   8 % if mptiles nonempry, add tile*10000 to num1
   9 % if dmtiles=[] mptiles=[] same as read_array_sp(root,num1,num2)
  10 % output in [i j a] format works even when some indices are <1
  11
  12 % Jan Mandel, October 2008
  13
  14 if ~exist('num1','var'),
  15     num1=-1;
  16 end
  17 if ~exist('num2','var'),
  18     num2=-1;
  19 end
  20 ff=file_name(root,num1,num2);
  21 if ~exist('dmtiles','var')
  22     dmtiles=[];
  23 end
  24 if ~exist('mptiles','var')
  25     mptiles=[];
  26 end
  27 if isempty(dmtiles) & isempty(mptiles)
  28    [i,j,a]=read_array_sp([ff,'.txt']);
  29 else
  30    i=[];j=[];a=[];tile=[];
  31    if isempty(dmtiles),
  32       dmtiles=-1;
  33    end
  34    if isempty(mptiles),
  35       mptiles=0;
  36    end
  37    for dmtile=dmtiles
  38       for mptile=mptiles
  39          if dmtile >= 0,
  40             f=sprintf('%s_%5.5i.%4.4i.txt',root,num1+10000*mptile,dmtile);
  41          else
  42             f=sprintf('%s_%5.5i.txt',root,num1+10000*mptile);
  43          end
  44          [ii,jj,aa]=read_array_sp(f);
  45          i=[i;ii]; j=[j;jj]; a=[a;aa];
  46          tile=[tile;ones(length(ii),1)*[dmtile,mptile]];
  47       end
  48    end
  49    % check for consistent duplicates
  50    isize=max(i)-min(i)+1;
  51    ij=i+isize*(j-1);    % coded pairs (i,j)
  52    [k,kix]=sort(ij);        % k=ij(kix)
  53    n=length(k);
  54    ikix=zeros(n,1);
  55    ikix(kix)=[1:n]';        % ij=k(ikix);
  56    same=find(k(2:end)==k(1:end-1));
  57    idiff=same(a(kix(same+1))~=a(kix(same)));
  58    id1=kix(idiff);
  59    id2=kix(idiff+1);
  60    if any(idiff)
  61       i1=i(id1);
  62       i2=i(id2);
  63       j1=j(id1);
  64       j2=j(id2);
  65       a1=a(id1);
  66       a2=a(id2);
  67       t1=tile(id1,:);
  68       t2=tile(id2,:);
  69       for m=1:length(idiff);
  70             fprintf('tile %i %i: a(%i,%i)=%g tile %i %i: a(%i,%i)=%g diff %g\n',...
  71                 t1(m,:),i1(m),j1(m),a1(m),t2(m,:),i2(m),j2(m),a2(m),a1(m)-a2(m))
  72       end
  73       warning('inconsistent values at overlap, taking the first seen')
  74    end
  75    % take out the duplicates
  76    i(kix(same))=[];
  77    j(kix(same))=[];
  78    a(kix(same))=[];
  79 end
  80 % create the output
  81 switch nargout
  82     case 1
  83         varargout{1}=sparse(i,j,a);
  84     case 3
  85         varargout{1}=i;
  86         varargout{2}=j;
  87         varargout{3}=a;
  88     otherwise
  89         error('bad number of output arguments')
  90 end
  91 end