linjpcut.m

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

function [ctimes, cval]=linjpcut(jmptimes, fval, cut_time) 
    
% LINJPCUT Truncate piecewise linear functions at a given time. 
%   N piecewise linear 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 break points. 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 break or jump
%   points and function values at these points respectively. If they have
%   different number of breaks (jumps), each vector is padded with the
%   last value to have the same length.
%
% [ctimes, cval] = linjpcut(jmptimes, fval, cut_time) 
%
% Inputs: 
%        jmptimes - a matrix with break or jump times stored columnwise
%        fval - a matrix with function values at the break (jump) points
%          stored columnwise 
%        cut_time - the cut-off time
%
% Outputs:
%        ctimes - a matrix with break or jump times of the truncated
%          functions stored columnwise
%        cval - a matrix of function values at the break (jump) points
%          stored columnwise 
%

% Authors: R.Gaigalas, I.Kaj
% v1.0 16-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);

  % compute the values at the cut-off time by linear interpolation 
  % indices of the first exceeding jump time in each column 
  % thanks to Peter J.Acklam
  % http://home.online.no/~pjacklam/matlab/doc/mtt/index.html 
  fexi = exi(logical(diff([0; exj]))); 
  fexij = sub2ind(size(jmptimes), fexi, [1:nproc]');
  lleij = sub2ind(size(jmptimes), max(fexi-1, 1), [1:nproc]');
  cutval = lineatpoints([jmptimes(lleij) jmptimes(fexij)], ...
                        [fval(lleij) fval(fexij)], cut_time);
  
  % fill each column with the last value
  cval(cex_ij) = cutval(exj);

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [yi] = lineatpoints(X, Y, xi)
%   Draw a line between each pair (x_1j, y_1j) and (x_2j, y_2j) and
%   compute the values at the points xi. A simple linear interpolation.
%

yi = (repmat(xi, size(X(:, 1)))-X(:, 1)).*(Y(:, 2)-Y(:, 1))./...
      (X(:, 2)-X(:, 1))+Y(:, 1);