staircut.m

用matalab开发的随机数生成程序包

function [ctimes, cval]=staircut(jmptimes, fval, cut_time)

% STAIRCUT Truncate piecewise constant functions at a given time. 
%   N piecewise constant functions "living" on the intervals [A, B_n],
%   n=1,...,N with the common start value A are given. The functions may
%   have different number of jumps. Truncate the functions at the time C,
%   so that they all live on the interval [A, C].
%
%   The functions are stored in two matrices containing the jumps points
%   and function values at these points respectively. If they have
%   different number of jumps, each vector is padded with the last value
%   to have the same length.
%
% [ctimes, cval] = staircut(jmptimes, fval, cut_time) 
%
% Inputs: 
%        jmptimes - a matrix with jump times stored columnwise
%        fval - a matrix with function values at the jump points stored 
%            columnwise 
%        cut_time - the cut-off time
%
% Outputs:
%        ctimes - a matrix with jump times of the truncated functions
%          stored columnwise 
%        cval - a matrix of function values at the jump points stored
%            columnwise 
%
% See also STAIRSUM, STAIRINTEGR, LINJPCUT.

% Authors: R.Gaigalas, I.Kaj
% v2.0 17-Oct-05

  nproc = size(jmptimes, 2);

  % cut off the jump times exceeding the limit  
  [ctimes, le_ij, cle_ij, exi, exj, cex_ij] = sttimescut(jmptimes, ...
                                                  cut_time);
  if (isempty(ctimes))  
    cval = [];
    return;
  end
 
  % copy the values satisfying the criterion
  cval = zeros(size(ctimes));
  cval(cle_ij) = fval(le_ij);
  
  % set the remaining values to the last value in each column
  % thanks to Peter J.Acklam
  % http://home.online.no/~pjacklam/matlab/doc/mtt/index.html  
  % indexes of the first exceeding jump time in each column 
  fexi = exi(logical(diff([0; exj]))); 
  % last not exceeding jump time
  lleij = sub2ind(size(jmptimes), max(fexi-1, 1), [1:nproc]');
  % if the k-th column contains j exceeding jump times then 
  % exj contains "k" replicated j times
  cval(cex_ij) = fval(lleij(exj));