Code covered by the BSD License  

Highlights from
CO2gui - lab control and automation

CO2gui - lab control and automation

by

 

06 Jan 2010 (Updated )

Software used for controlling and data logging lab equipment.

slmeval(x,slm,evalmode)
function y = slmeval(x,slm,evalmode)
% slmeval: evaluates a Hermite function, its derivatives, or its inverse
% usage 1: y = slmeval(x,slm)           % evaluates y = f(x)
% usage 2: y = slmeval(x,slm,evalmode)  % general form
%
% As opposed to extrapolation as I am, slmeval will not
% extrapolate. If you wanted to extrapolate, then you
% should have built the spline differently.
% 
% arguments: (input)
%  x       - array (any shape) of data to be evaluated through model
%
%            All points are assumed to lie inside the range of
%            the knots. Any which fall outside the first or last
%            knot will be assigned the value of the function
%            at the corresponding knot. NO extrapolation is
%            allowed. (If you really need to extrapolate, then
%            you should have done it when you built the model
%            in the first place. An alternative, if I ever write
%            it, is an extrapolator code, that builds a new,
%            extrapolated function.)
%
%  slm     - a shape language model structure, normally constructed
%            by slmfit or slmengine.
%
%            The fields in this struct are:
%            slm.type  = either 'cubic', 'linear', or 'constant'
%
%            slm.knots = a vector of knots (distinct & increasing)
%               There must be at least two knots.
%
%            slm.coef  = an array of coefficients for the hermite
%               function.
%
%            If the function is linear or constant Hermite, then
%            slm.coef will have only one column, composed of the
%            value of the function at each corresponding knot.
%            The 'constant' option uses a function value at x(i)
%            to apply for x(i) <= x < x(i+1).
%
%            If the function is cubic Hermite, then slm.coef
%            will have two columns. The first column will be the
%            value of the function at the corresponding knot,
%            the second column will be the correspodign first
%            derivative.
%
% evalmode - (OPTIONAL) numeric flag - specifies what evaluation
%            is to be done.
%
%            DEFAULT VALUE: 0
%
%            == 0 --> evaluate the function at each point in x.
%            == 1 --> the first derivative at each point in x.
%            == 2 --> the second derivative at each point in x.
%            == 3 --> the third derivative at each point in x.
%            == -1 --> evaluate the inverse of the function at
%             each point in x, thus y is returned such that x=f(y)
%
%            Note 1: Piecewise constant functions will return zero
%            for all order derivatives, since I ignore the delta
%            functions at each knot.
%            The inverse operation is also disabled for constant
%            functions.
%
%            Note 2: Linear hermite functions will return zero
%            for second and higher order derivatives. At a knot
%            point, while technically the derivative is undefined
%            at that location, the slope of the segment to the
%            right of that knot is returned.
%
%            Note 3: Inverse computations will return the
%            LEFTMOST zero (closest to -inf) in the event that
%            more than one solution exists.
%
%            Note 4: Inverse of points which fall above the
%            maximum or below the minimum value of the function
%            will be returned as a NaN.
%
% Arguments: (output)
%  y       - Evaluated result, the predicted value, i.e., f(x)
%            or f'(x), f''(x), f'''(x), or the functional inverse
%            such that x = f(y). y will have the same shape and
%            size as x.
%
% Example:
%




% See also: slmpar, ppval, slmengine
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 6/10/09

% get the shape of x, so we can restore it at the end
sizex = size(x);
x=x(:);
nx = length(x);

% check for the minimum required information
if ~isstruct(slm) || ~isfield(slm,'knots') || ...
    ~isfield(slm,'coef') || ~isfield(slm,'degree')
  error('slm must be a struct, with "knots" "coef" and "degree" fields')
end

knots = slm.knots(:);
coef = slm.coef;
nk = length(knots);
dx = diff(knots);

% default evalmode is 0, i.e., evaluate f(x)
if (nargin<3) || isempty(evalmode)
  evalmode = 0;
end

% evaluate the appropriate class of piecewise hermite function
switch slm.degree
  case 1
    % verify that coef is the right size
    sc=size(coef);
    if (sc(1) ~= nk)
      error('Improper size of coef field for these knots')
    elseif (sc(2)~=1)
      error('Improper size of coef field for a linear Hermite')
    end

    % which knot interval does each point fall in?
    % only treat evalmode == 0 or 1 right now
    if ismember(evalmode,[0 1])
      % clip the points first
      x = min(knots(end),max(x,knots(1)));

      % and use histc to bin the points.
      [junk,xbin] = histc(x,knots);

      % any point which falls at the top end, is said to
      % be in the last bin.
      xbin(xbin==nk)=nk-1;

    else
      % we need to treat the inverse carefully, and we
      % can ignore where a point falls for any higher
      % derivatives requested. do this part below.
      xbin = zeros(nx,1);

    end

    % Evaluation mode:
    switch evalmode
      case 0
        % f(x)
        t = (x-knots(xbin))./dx(xbin);
        y = coef(xbin) + (coef(xbin+1)-coef(xbin)).*t;

      case 1
        % first derivative
        y = (coef(xbin+1)-coef(xbin))./dx(xbin);

      case [2 3]
        % higher derivatives of a piecewise linear function
        % are all zero
        y = zeros(sizex);

      case -1
        % functional inverse
        % here the problem of which bin we fall into is not
        % trivial, since the function need not be monotone.
        % we might also have constant segments.
        % Identify the leftmost interval which contains the
        % point in question.

        % first, clip the data in terms of function value
        miny = min(coef);
        maxy = max(coef);
        y = repmat(NaN,nx,1);
        k = find((x<=maxy) & (x>=miny));

        % determine which knot interval to look in for
        % each point.
        ybin = nmbs(x(k),coef);

        % rule out the divide by zero cases for intervals
        % where the function was constant
        L = (coef(ybin+1) == coef(ybin));
        if any(L)
          y(k(L)) = knots(ybin(L));
          ybin(L)=[];
          k(L)=[];
        end

        % inverse interpolation
        y(k) = knots(ybin) + (x(k) - coef(ybin)).* ...
          (knots(ybin+1)-knots(ybin))./(coef(ybin+1)-coef(ybin));

      otherwise
        % anything else
        error('Evalmode must be one of [0 1 2 3 -1]')
    end

  case 3
    % verify shape of coef
    sc=size(coef);
    if (sc(1) ~= nk)
      error('Improper size of coef field for these knots')
    elseif (sc(2)~=2)
      error('Improper size of coef field for a cubic Hermite')
    end

    % which knot interval does each point fall in?
    % only treat evalmode == 0,1,2,3 right now
    if ismember(evalmode,[0 1 2 3])
      % clip the points first
      x = min(knots(end),max(x,knots(1)));

      % and use histc to bin the points.
      [junk,xbin] = histc(x,knots);

      % any point which falls at the top end is said to
      % be in the last bin.
      xbin(xbin==nk)=nk-1;

    else
      % we need to treat the inverse more carefully
      xbin = zeros(nx,1);

    end

    % Evaluation mode:
    switch evalmode
      case 0
        % f(x)
        t = (x-knots(xbin))./dx(xbin);
        t2 = t.^2;
        t3 = t.^3;
        s2 = (1-t).^2;
        s3 = (1-t).^3;
        y = (-coef(xbin,2).*(s3-s2) + ...
          coef(xbin+1,2).*(t3-t2)).*dx(xbin) + ...
          coef(xbin,1).*(3*s2-2*s3) + ...
          coef(xbin+1,1).*(3*t2-2*t3);
      case 1
        % first derivative for the cubic case
        t = (x-knots(xbin))./dx(xbin);
        t2 = t.^2;
        s = 1-t;
        s2 = (1-t).^2;
        y = -coef(xbin,2).*(-3*s2+2*s) + ...
          coef(xbin+1,2).*(3*t2-2*t) + ...
          (coef(xbin,1).*(-6*s+6*s2) + ...
          coef(xbin+1,1).*(6*t-6*t2))./dx(xbin);
      case 2
        % second derivative of a cubic
        t = (x-knots(xbin))./dx(xbin);
        s = 1-t;
        y = (-coef(xbin,2).*(6*s - 2) + ...
          coef(xbin+1,2).*(6*t - 2))./dx(xbin) + ...
          (coef(xbin,1).*(6 - 12*s) + ...
          coef(xbin+1,1).*(6 - 12*t))./(dx(xbin).^2);
      case 3
        % third derivative
        y = 6*(coef(xbin,2) + coef(xbin+1,2))./(dx(xbin).^2) + ...
          12*(coef(xbin,1) - coef(xbin+1,1))./(dx(xbin).^3);
      case -1
        % functional inverse
        % here the problem of which bin we fall into is not
        % trivial, since the function need not be monotone.
        % we might also have constant segments.
        
        % first, convert the spline into a pp form.
        pp = slm2pp(slm);
        coefs = pp.coefs;
        
        % scale the cubic polys so they live on [0,1].
        % this will stabilize things a bit, as well as make
        % the tests easier later on. I could do this in one
        % line with a variety of tricks, but feeling lazy...
        coefs(:,1) = coefs(:,1).*(dx.^3);
        coefs(:,2) = coefs(:,2).*(dx.^2);
        coefs(:,3) = coefs(:,3).*dx;
        
        % Identify the leftmost interval which contains the
        % point in question. The problem is the cubic segments
        % might not be monotone. nmbs will try though.
        binedges = slm.coef(:,1);
        ybin = nmbs(x,binedges);
        
        % a simple solution is to test every interval up to
        % and including the interval found in ybin, but no
        % further. We can stop there since we know a solution
        % must exist in that interval. If no interval was found,
        % then ybin will be a NaN. In that event, we must search
        % every interval, since the curve need not be monotone.
        ybin(isnan(ybin) | (ybin == nk)) = nk - 1;
        
        % just a loop over roots here. not terribly efficient,
        % but I don't terribly want to code a vectorized cubic
        % solver for this problem.
        y = NaN(size(x));
        tol = 1000*eps;
        for j = 1:nx
          flag = true;
          i = 1;
          while flag && (i <= ybin(j))
            Ci = coefs(i,:);
            
            % offset the constant term. a root of this
            % cubic is what we want.
            Ci(4) = Ci(4) - x(j);
            
            % scale the poly coefficients to improve things
            % yet some more.
            Ci = Ci./max(abs(Ci(1:3)));
            
            % get the roots. at least one must be real.
            Ri = roots(Ci);
            
            % of the real roots, is one of them strictly
            % inside [0,1]? If more than one is, take the
            % smallest root.
            k = (abs(imag(Ri)) == 0) & (real(Ri) >= 0) & (real(Ri) <= 1);
            if any(k)
              k = find(k,1,'first');
              y(j) = knots(i) + dx(i)*real(Ri(k));
              flag = false;
              continue
            end
            
            % if we did not succeed in the last test,
            % of the real roots, is one of them within 
            % inside [-tol,1+tol]? If more than one is,
            % take the smallest root.
            k = (abs(imag(Ri)) == 0) & (real(Ri) >= -tol) & (real(Ri) <= (1 + tol));
            if any(k)
              k = find(k,1,'first');
              y(j) = knots(i) + dx(i)*real(Ri(k));
              flag = false;
              continue
            end
            
            % increment i, check the next interval.
            i = i + 1;
          end
        end
        
      otherwise
        % anything else
        error('Evalmode must be one of [0 1 2 3 -1]')
    end
  case 0
    % piecewise constant function (discontinuous at the knots)
    % verify that coef is the right size
    sc=size(coef);
    if (sc(1) ~= (nk-1))
      error('Improper size of coef field for these knots')
    elseif (sc(2)~=1)
      error('Improper size of coef field for a piecewise constant function')
    end

    % Evaluation mode:
    switch evalmode
      case 0
        % f(x)

        % clip the points first
        x = min(knots(end),max(x,knots(1)));

        % which knot interval does each point fall in?
        % and use histc to bin the points.
        [junk,xbin] = histc(x,knots);

        % any point which falls at the top end, is said to
        % be in the last bin.
        xbin(xbin==nk)=nk-1;

        y = coef(xbin);

      case [1 2 3]
        % all derivatives of a piecewise constant function
        % will be zero if we ignore the delta functions
        % at each knot.
        y = zeros(sizex);

      case -1
        % functional inverse - no inverse exists for a
        % discontinuous function.
        y = reshape(NaN,sizex);

      otherwise
        % anything else
        error('Evalmode must be one of [0 1 2 3 -1]')
    end

  otherwise
    error('slmeval only handles the constant, linear, or cubic Hermite case')
end

% when all done, make sure the size of the output is the
% same as the input was.
y = reshape(y,sizex);


% ===================================================
% ================ begin subfunction ================
% ===================================================
function bind = nmbs(x,binedges)
% nmbs: non-monotone bin search
%
% Finds the leftmost bin that contains each x. Assumes
% that x has already been clipped so it must lie in
% some bin.
nb = length(binedges);
nx = length(x);

% if the bins are actually monotone, then just use histc
db = diff(binedges);
if all(db>0)
  % increasing bins
  [junk,bind] = histc(x,binedges);
  bind(bind==nb)=nb-1;

elseif all(db<0)
  % decreasing sequence of edges
  bind = histc(-x,-binedges);
  bind(bind==nb)=nb-1;

else
  % non-monotone sequence of edges. Do this one the
  % hard way. Find the first bin that fits.
  i = 1;
  j = 1:nx;
  bind = ones(nx,1);
  while ~isempty(j) && (i<nb)
    k = ((binedges(i)>=x(j)) & (binedges(i+1)<=x(j))) | ...
      ((binedges(i)<=x(j)) & (binedges(i+1)>=x(j)));

    if any(k)
      bind(j(k)) = i;
      j(k)=[];
    end

    % increment bin we will look in
    i=i+1;

  end

end






Contact us