function crray = aggread(chr, reg, loc, slc, cmut, cstrn, strn, ustrn, zs )
% compute quantities (read, mutation counts, etc.) at each location from aligned reads
%   Detailed explanation goes here
[dat.slc, idx] = unique(slc);
dat.reg = reg(idx);
dat.loc = loc(idx);
dat.chr = repmat(chr,length(idx),1);
dat.strn = logical(strn(idx));

m = length(slc);
n = length(dat.slc);
dat.mut = zeros(n,3,'int16');

dat.read = int16(histc(slc,dat.slc));
czs = mat2cell(zs,m,1);
dat.zs = cellfun(@(x) histc(slc(x),dat.slc), czs,'UniformOutput',false);
dat.zs = int16(cat(2,dat.zs{:}));

for k = 1:2
    pmut = nonzeros(cmut(cstrn == ustrn(k),:));
    umut = unique(pmut);
    freq = int16(histc(pmut,umut));
    dat.mut(dat.strn == (k-1) & ismember(dat.loc,umut)) = freq;
    
    if ~issorted(dat.loc(dat.strn == (k-1)))
        error('Locations are not sorted.');
    end
end

crray = {dat.chr, dat.loc, dat.mut, dat.read, ...
    dat.reg, dat.slc, dat.strn, dat.zs};

end