function [noisydata, indaffected] = simNoise(cleandata, nmissing, varargin)
% Simulates some statistical difficulties of a data set.
% Noise, elementwise outliers, and missing element values can be
% introduced in a formerly clean data. The alterations will be
% performed to randomly selected variables (As of yet, there's no way
% to select just some observations; they are treated equally here.)
%
% Input:
%
%   cleandata - the clean data to be made bad; columns are observations,
%               rows are variables.
%
%   nmissing  - number of missing values. This many randomly chosen values
%               will be replaced by NaN.
%
%   After the above-mentioned, you may give any number of sets of three
%   arguments, whose meanings are (in this order): "noise distribution",
%   "number affected", and "amplitude". The algorithm is as follows: 
%
%    "noise distribution" is a string; it must be either 'normal' to
%        induce normally distributed noise (randn(), centered on 0) or
%        'uniform' to induce uniformly distributed noise (rand(), range
%        [-1,1]). 
%
%    "number affected" must be an integer ranging from 0 to the number
%        of individual values in the clean data. Additionally, the sum of
%        all "numbers affected" plus the number of missing values
%        (nmissing; see above) must not exceed the total number of data
%        values.
%
%        You will get this many values randomly shifted by a random
%        value that is drawn from "noise distribution" and multiplied by
%        "amplitude".
%
%    "amplitude" can be either a scalar value (all variables shifted
%        with the given amplitude) or a row vector that gives an individual
%        amplitude for each of the variables. In this way you can have
%        different noise levels for different variables.
%
% Output:
%
%   noisydata   - a mangled data set that carries all of the requested
%                 complications. 
%
%   indaffected - a cell array of matrices (whose sizes are equal to that
%                 of the data) that remembers the affected data values as 
%                 binary matrices (0=not affected, 1=affected). Ordering of
%                 the cell array is the same as your final arguments to the
%                 function call. NOTE: Missing values are not included in
%                 indaffected; you can easily locate them by searching
%                 NaNs in the produced noisy data.
%
% Note:
%
%   The alterations are not cumulative: Any value is only affected by one
%   type of alteration. The "numbers affected" and the number of completely
%   missing values (nmissing) therefore must sum up to no more than the 
%   total number of values.
%
%   The indices of the altered values are randomly generated.
%
%   This affects only single variables within the data set. In order to
%   affect a number of observations (including all the variables), you must
%   use some other means.
%
%
%
% Example:
%
%   inputvecs = simnoise(inputvecs, 5, ...
%      'normal', 25, 0.3, ...
%      'uniform', 7, 0.6);
%
%   The above example would convert 5 of the elements of the matrix
%   *inputvecs* to NaN, perturb 25 elements with normally distributed
%   noise with amplification factor of 0.3, and 7 of the remaining
%   elements with uniformly distributed noise in the interval
%   [-0.6,+0.6].
%
%
%     TODO: Make it possible to mask missing values to be drawn from only 
%           some subset of the variables (or of observations?) - would be 
%           more realistic?


nvars = size(cleandata, 1);
ndata = size(cleandata, 2);
nvalues = nvars*ndata;

% Skim through varargin. Save info to indaffected{} and to distrib{} and to
% amp{}. Errors if arguments are badly formatted. Check that total number
% doesn't exceed that of the data.
va = varargin;

requestedtotal = nmissing;

j = 0;
for i = 1:3:length(va)
    s = va{i};
    num = va{i+1};
    a = va{i+2};
    if (~ischar(s))
        error('Check argument list. Expected a string (''normal'' or ''uniform'')');
    end
    if (~(strcmp(s, 'normal') || strcmp(s, 'uniform')))
        error('Check argument spelling. Expected either ''normal'' or ''uniform''');
    end
    if (~isnumeric(num) || ~isnumeric(a))
        error('Check argument list. Expected a numerical value');
    end
    if ((1 < size(a,1)) && (size(a,1) < nvars))
        error('A given amplitude is not scalar but its length doesn''t match the number of variables.');
    end

    % Enque a valid alteration (NOTE: those with zero amplitude are included!):
    j = j + 1;
    distrib{j} = s;
    numaff{j} = num;
    amp{j} = a;
    requestedtotal = requestedtotal + num;
end

if (requestedtotal > nvalues)
    error(sprintf('Total number of values is %d. Cannot alter as many as %d.', ...
        nvalues, requestedtotal));
end

% Prepare for selection of indices: make an array of alteration types.
% Warn if there is a request that actually makes no change to the data.
beg_i = 1;
for i = 1:j
    end_i = beg_i + numaff{i} - 1;
    typeofchange(beg_i:end_i) = i;
    beg_i = end_i + 1;
    if (numaff{i}==0)
        warning('Number of affected values for an alteration is set to zero.');
    elseif (amp{i}==0)
        warning('Noise amplitude for an alteration is set to zero.');
    end
end
if (nmissing > 0)
    end_i = beg_i + nmissing - 1;
    typeofchange(beg_i:end_i) = j+1; 
end
if (end_i < nvalues)
    typeofchange(end_i+1:nvalues) = 0;
end

% Shuffle requested alteration types and pack them into a matrix:
rp = randperm(nvalues);
typeofchange = typeofchange(rp);
typeofchange = reshape(typeofchange, nvars, ndata);

% Perform the alterations in a loop.
noisydata = cleandata;
for i = 1:j
    subjects = find(typeofchange == i);
    indaffected{i} = zeros(nvars, ndata);
    indaffected{i}(subjects) = 1;
    a(1:nvars) = amp{i};
    A = repmat(a', 1, ndata);

    if (strcmp(distrib{i},'normal'))
        noisydata(subjects) = noisydata(subjects) ...
            + A(subjects) .* randn(size(noisydata(subjects)));

    elseif (strcmp(distrib{i},'uniform'))
        rsize = size(noisydata(subjects));
        noisydata(subjects) = noisydata(subjects) ...
          + A(subjects) .* (2* (rand(rsize) - 0.5*ones(rsize)) );
	
    end
end

% And finally create missing values:
subjects = find(typeofchange == j+1);
noisydata(subjects) = NaN;
