text:
function [zgrid,xgrid,ygrid] = gridfit(x,y,z,xnodes,ynodes,varargin)
% gridfit: estimates a surface on a 2d grid, based on scattered data
%          Replicates are allowed. All methods extrapolate to the grid
%          boundaries. Gridfit uses a modified ridge estimator to
%          generate the surface, where the bias is toward smoothness.
%
%          Gridfit is not an interpolant. Its goal is a smooth surface
%          that approximates your data, but allows you to control the
%          amount of smoothing.
%
% usage #1: zgrid = gridfit(x,y,z,xnodes,ynodes);
% usage #2: [zgrid,xgrid,ygrid] = gridfit(x,y,z,xnodes,ynodes);
% usage #3: zgrid = gridfit(x,y,z,xnodes,ynodes,prop,val,prop,val,...);
%
% Arguments: (input)
%  x,y,z - vectors of equal lengths, containing arbitrary scattered data
%          The only constraint on x and y is they cannot ALL fall on a
%          single line in the x-y plane. Replicate points will be treated
%          in a least squares sense.
%
%          ANY points containing a NaN are ignored in the estimation
%
%  xnodes - vector defining the nodes in the grid in the independent
%          variable (x). xnodes need not be equally spaced. xnodes
%          must completely span the data. If they do not, then the
%          'extend' property is applied, adjusting the first and last
%          nodes to be extended as necessary. See below for a complete
%          description of the 'extend' property.
%
%          If xnodes is a scalar integer, then it specifies the number
%          of equally spaced nodes between the min and max of the data.
%
%  ynodes - vector defining the nodes in the grid in the independent
%          variable (y). ynodes need not be equally spaced.
%
%          If ynodes is a scalar integer, then it specifies the number
%          of equally spaced nodes between the min and max of the data.
%
%          Also see the extend property.
%
%  Additional arguments follow in the form of property/value pairs.
%  Valid properties are:
%    'smoothness', 'interp', 'regularizer', 'solver', 'maxiter'
%    'extend', 'tilesize', 'overlap'
%
%  Any UNAMBIGUOUS shortening (even down to a single letter) is
%  valid for property names. All properties have default values,
%  chosen (I hope) to give a reasonable result out of the box.
%
%   'smoothness' - scalar - determines the eventual smoothness of the
%          estimated surface. A larger value here means the surface
%          will be smoother. Smoothness must be a non-negative real
%          number.
%
%          Note: the problem is normalized in advance so that a
%          smoothness of 1 MAY generate reasonable results. If you
%          find the result is too smooth, then use a smaller value
%          for this parameter. Likewise, bumpy surfaces suggest use
%          of a larger value. (Sometimes, use of an iterative solver
%          with too small a limit on the maximum number of iterations
%          will result in non-convergence.)
%
%          DEFAULT: 1
%
%
%   'interp' - character, denotes the interpolation scheme used
%          to interpolate the data.
%
%          DEFAULT: 'triangle'
%
%          'bilinear' - use bilinear interpolation within the grid
%                     (also known as tensor product linear interpolation)
%
%          'triangle' - split each cell in the grid into a triangle,
%                     then linear interpolation inside each triangle
%
%          'nearest' - nearest neighbor interpolation. This will
%                     rarely be a good choice, but I included it
%                     as an option for completeness.
%
%
%   'regularizer' - character flag, denotes the regularization
%          paradignm to be used. There are currently three options.
%
%          DEFAULT: 'gradient'
%
%          'diffusion' or 'laplacian' - uses a finite difference
%              approximation to the Laplacian operator (i.e, del^2).
%
%              We can think of the surface as a plate, wherein the
%              bending rigidity of the plate is specified by the user
%              as a number relative to the importance of fidelity to
%              the data. A stiffer plate will result in a smoother
%              surface overall, but fit the data less well. I've
%              modeled a simple plate using the Laplacian, del^2. (A
%              projected enhancement is to do a better job with the
%              plate equations.)
%
%              We can also view the regularizer as a diffusion problem,
%              where the relative thermal conductivity is supplied.
%              Here interpolation is seen as a problem of finding the
%              steady temperature profile in an object, given a set of
%              points held at a fixed temperature. Extrapolation will
%              be linear. Both paradigms are appropriate for a Laplacian
%              regularizer.
%
%          'gradient' - attempts to ensure the gradient is as smooth
%              as possible everywhere. Its subtly different from the
%              'diffusion' option, in that here the directional
%              derivatives are biased to be smooth across cell
%              boundaries in the grid.
%
%              The gradient option uncouples the terms in the Laplacian.
%              Think of it as two coupled PDEs instead of one PDE. Why
%              are they different at all? The terms in the Laplacian
%              can balance each other.
%
%          'springs' - uses a spring model connecting nodes to each
%              other, as well as connecting data points to the nodes
%              in the grid. This choice will cause any extrapolation
%              to be as constant as possible.
%
%              Here the smoothing parameter is the relative stiffness
%              of the springs connecting the nodes to each other compared
%              to the stiffness of a spting connecting the lattice to
%              each data point. Since all springs have a rest length
%              (length at which the spring has zero potential energy)
%              of zero, any extrapolation will be minimized.
%
%          Note: I don't terribly like the 'springs' strategy.
%          It tends to drag the surface towards the mean of all
%          the data. Its been left in only because the paradigm
%          interests me.
%
%
%   'solver' - character flag - denotes the solver used for the
%          resulting linear system. Different solvers will have
%          different solution times depending upon the specific
%          problem to be solved. Up to a certain size grid, the
%          direct \ solver will often be speedy, until memory
%          swaps causes problems.
%
%          What solver should you use? Problems with a significant
%          amount of extrapolation should avoid lsqr. \ may be
%          best numerically for small smoothnesss parameters and
%          high extents of extrapolation.
%
%          Large numbers of points will slow down the direct
%          \, but when applied to the normal equations, \ can be
%          quite fast. Since the equations generated by these
%          methods will tend to be well conditioned, the normal
%          equations are not a bad choice of method to use. Beware
%          when a small smoothing parameter is used, since this will
%          make the equations less well conditioned.
%
%          DEFAULT: 'normal'
%
%          '\' - uses matlab's backslash operator to solve the sparse
%                     system. 'backslash' is an alternate name.
%
%          'symmlq' - uses matlab's iterative symmlq solver
%
%          'lsqr' - uses matlab's iterative lsqr solver
%
%          'normal' - uses \ to solve the normal equations.
%
%
%   'maxiter' - only applies to iterative solvers - defines the
%          maximum number of iterations for an iterative solver
%
%          DEFAULT: min(10000,length(xnodes)*length(ynodes))
%
%
%   'extend' - character flag - controls whether the first and last
%          nodes in each dimension are allowed to be adjusted to
%          bound the data, and whether the user will be warned if
%          this was deemed necessary to happen.
%
%          DEFAULT: 'warning'
%
%          'warning' - Adjust the first and/or last node in
%                     x or y if the nodes do not FULLY contain
%                     the data. Issue a warning message to this
%                     effect, telling the amount of adjustment
%                     applied.
%
%          'never'  - Issue an error message when the nodes do
%                     not absolutely contain the data.
%
%          'always' - automatically adjust the first and last
%                     nodes in each dimension if necessary.
%                     No warning is given when this option is set.
%
%
%   'tilesize' - grids which are simply too large to solve for
%          in one single estimation step can be built as a set
%          of tiles. For example, a 1000x1000 grid will require
%          the estimation of 1e6 unknowns. This is likely to
%          require more memory (and time) than you have available.
%          But if your data is dense enough, then you can model
%          it locally using smaller tiles of the grid.
%
%          My recommendation for a reasonable tilesize is
%          roughly 100 to 200. Tiles of this size take only
%          a few seconds to solve normally, so the entire grid
%          can be modeled in a finite amount of time. The minimum
%          tilesize can never be less than 3, although even this
%          size tile is so small as to be ridiculous.
%
%          If your data is so sparse than some tiles contain
%          insufficient data to model, then those tiles will
%          be left as NaNs.
%
%          DEFAULT: inf
%
%
%   'overlap' - Tiles in a grid have some overlap, so they
%          can minimize any problems along the edge of a tile.
%          In this overlapped region, the grid is built using a
%          bi-linear combination of the overlapping tiles.
%
%          The overlap is specified as a fraction of the tile
%          size, so an overlap of 0.20 means there will be a 20%
%          overlap of successive tiles. I do allow a zero overlap,
%          but it must be no more than 1/2.
%
%          0 <= overlap <= 0.5
%
%          Overlap is ignored if the tilesize is greater than the
%          number of nodes in both directions.
%
%          DEFAULT: 0.20
%
%
%   'autoscale' - Some data may have widely different scales on
%          the respective x and y axes. If this happens, then
%          the regularization may experience difficulties.
%
%          autoscale = 'on' will cause gridfit to scale the x
%          and y node intervals to a unit length. This should
%          improve the regularization procedure. The scaling is
%          purely internal.
%
%          autoscale = 'off' will disable automatic scaling
%
%          DEFAULT: 'on'
%
%
% Arguments: (output)
%  zgrid   - (nx,ny) array containing the fitted surface
%
%  xgrid, ygrid - as returned by meshgrid(xnodes,ynodes)
%
%
% Speed considerations:
%  Remember that gridfit must solve a LARGE system of linear
%  equations. There will be as many unknowns as the total
%  number of nodes in the final lattice. While these equations
%  may be sparse, solving a system of 10000 equations may take
%  a second or so. Very large problems may benefit from the
%  iterative solvers or from tiling.
%
%
% Example usage:
%
%  x = rand(100,1);
%  y = rand(100,1);
%  z = exp(x+2*y);
%  xnodes = 0:.1:1;
%  ynodes = 0:.1:1;
%
%  g = gridfit(x,y,z,xnodes,ynodes);
%
% Note: this is equivalent to the following call:
%
%  g = gridfit(x,y,z,xnodes,ynodes, ...
%              'smooth',1, ...
%              'interp','triangle', ...
%              'solver','normal', ...
%              'regularizer','gradient', ...
%              'extend','warning', ...
%              'tilesize',inf);
%
%
% Author: John D'Errico
% e-mail address: woodchips@rochester.rr.com
% Release: 2.0
% Release date: 5/23/06

% set defaults
params.smoothness = 1;
params.interp = 'triangle';
params.regularizer = 'gradient';
params.solver = 'normal';
params.maxiter = [];
params.extend = 'warning';
params.tilesize = inf;
params.overlap = 0.20;
params.mask = [];
params.autoscale = 'on';
params.xscale = 1;
params.yscale = 1;

% was the params struct supplied?
if ~isempty(varargin)
  if isstruct(varargin{1})
    % params is only supplied if its a call from tiled_gridfit
    params = varargin{1};
    if length(varargin)>1
      % check for any overrides
      params = parse_pv_pairs(params,varargin{2:end});
    end
  else
    % check for any overrides of the defaults
    params = parse_pv_pairs(params,varargin);

  end
end

% check the parameters for acceptability
params = check_params(params);

% ensure all of x,y,z,xnodes,ynodes are column vectors,
% also drop any NaN data
x=x(:);
y=y(:);
z=z(:);
k = isnan(x) | isnan(y) | isnan(z);
if any(k)
  x(k)=[];
  y(k)=[];
  z(k)=[];
end
xmin = min(x);
xmax = max(x);
ymin = min(y);
ymax = max(y);

% did they supply a scalar for the nodes?
if length(xnodes)==1
  xnodes = linspace(xmin,xmax,xnodes)';
  xnodes(end) = xmax; % make sure it hits the max
end
if length(ynodes)==1
  ynodes = linspace(ymin,ymax,ynodes)';
  ynodes(end) = ymax; % make sure it hits the max
end

xnodes=xnodes(:);
ynodes=ynodes(:);
dx = diff(xnodes);
dy = diff(ynodes);
nx = length(xnodes);
ny = length(ynodes);
ngrid = nx*ny;

% set the scaling if autoscale was on
if strcmpi(params.autoscale,'on')
  params.xscale = mean(dx);
  params.yscale = mean(dy);
  params.autoscale = 'off';
end

% check to see if any tiling is necessary
if (params.tilesize < max(nx,ny))
  % split it into smaller tiles. compute zgrid and ygrid
  % at the very end if requested
  zgrid = tiled_gridfit(x,y,z,xnodes,ynodes,params);
else
  % its a single tile.

  % mask must be either an empty array, or a boolean
  % aray of the same size as the final grid.
  nmask = size(params.mask);
  if ~isempty(params.mask) && ((nmask(2)~=nx) || (nmask(1)~=ny))
    if ((nmask(2)==ny) || (nmask(1)==nx))
      error 'Mask array is probably transposed from proper orientation.'
    else
      error 'Mask array must be the same size as the final grid.'
    end
  end
  if ~isempty(params.mask)
    params.maskflag = 1;
  else
    params.maskflag = 0;
  end

  % default for maxiter?
  if isempty(params.maxiter)
    params.maxiter = min(10000,nx*ny);
  end

  % check lengths of the data
  n = length(x);
  if (length(y)~=n) || (length(z)~=n)
    error 'Data vectors are incompatible in size.'
  end
  if n<3
    error 'Insufficient data for surface estimation.'
  end

  % verify the nodes are distinct
  if any(diff(xnodes)<=0) || any(diff(ynodes)<=0)
    error 'xnodes and ynodes must be monotone increasing'
  end

  % do we need to tweak the first or last node in x or y?
  if xmin<xnodes(1)
    switch params.extend
      case 'always'
        xnodes(1) = xmin;
      case 'warning'
        warning(['xnodes(1) was decreased by: ',num2str(xnodes(1)-xmin),', new node = ',num2str(xmin)])
        xnodes(1) = xmin;
      case 'never'
        error(['Some x (',num2str(xmin),') falls below xnodes(1) by: ',num2str(xnodes(1)-xmin)])
    end
  end
  if xmax>xnodes(end)
    switch params.extend
      case 'always'
        xnodes(end) = xmax;
      case 'warning'
        warning(['xnodes(end) was increased by: ',num2str(xmax-xnodes(end)),', new node = ',num2str(xmax)])
        xnodes(end) = xmax;
      case 'never'
        error(['Some x (',num2str(xmax),') falls above xnodes(end) by: ',num2str(xmax-xnodes(end))])
    end
  end
  if ymin<ynodes(1)
    switch params.extend
      case 'always'
        ynodes(1) = ymin;
      case 'warning'
        warning(['ynodes(1) was decreased by: ',num2str(ynodes(1)-ymin),', new node = ',num2str(ymin)])
        ynodes(1) = ymin;
      case 'never'
        error(['Some y (',num2str(ymin),') falls below ynodes(1) by: ',num2str(ynodes(1)-ymin)])
    end
  end
  if ymax>ynodes(end)
    switch params.extend
      case 'always'
        ynodes(end) = ymax;
      case 'warning'
        warning(['ynodes(end) was increased by: ',num2str(ymax-ynodes(end)),', new node = ',num2str(ymax)])
        ynodes(end) = ymax;
      case 'never'
        error(['Some y (',num2str(ymax),') falls above ynodes(end) by: ',num2str(ymax-ynodes(end))])
    end
  end

  % determine which cell in the array each point lies in
  [junk,indx] = histc(x,xnodes); %#ok
  [junk,indy] = histc(y,ynodes); %#ok
  % any point falling at the last node is taken to be
  % inside the last cell in x or y.
  k=(indx==nx);
  indx(k)=indx(k)-1;
  k=(indy==ny);
  indy(k)=indy(k)-1;
  ind = indy + ny*(indx-1);

  % Do we have a mask to apply?
  if params.maskflag
    % if we do, then we need to ensure that every
    % cell with at least one data point also has at
    % least all of its corners unmasked.
    params.mask(ind) = 1;
    params.mask(ind+1) = 1;
    params.mask(ind+ny) = 1;
    params.mask(ind+ny+1) = 1;
  end

  % interpolation equations for each point
  tx = min(1,max(0,(x - xnodes(indx))./dx(indx)));
  ty = min(1,max(0,(y - ynodes(indy))./dy(indy)));
  % Future enhancement: add cubic interpolant
  switch params.interp
    case 'triangle'
      % linear interpolation inside each triangle
      k = (tx > ty);
      L = ones(n,1);
      L(k) = ny;

      t1 = min(tx,ty);
      t2 = max(tx,ty);
      A = sparse(repmat((1:n)',1,3),[ind,ind+ny+1,ind+L], ...
        [1-t2,t1,t2-t1],n,ngrid);

    case 'nearest'
      % nearest neighbor interpolation in a cell
      k = round(1-ty) + round(1-tx)*ny;
      A = sparse((1:n)',ind+k,ones(n,1),n,ngrid);

    case 'bilinear'
      % bilinear interpolation in a cell
      A = sparse(repmat((1:n)',1,4),[ind,ind+1,ind+ny,ind+ny+1], ...
        [(1-tx).*(1-ty), (1-tx).*ty, tx.*(1-ty), tx.*ty], ...
        n,ngrid);

  end
  rhs = z;

  % Build regularizer. Add del^4 regularizer one day.
  switch params.regularizer
    case 'springs'
      % zero "rest length" springs
      [i,j] = meshgrid(1:nx,1:(ny-1));
      ind = j(:) + ny*(i(:)-1);
      m = nx*(ny-1);
      stiffness = 1./(dy/params.yscale);
      Areg = sparse(repmat((1:m)',1,2),[ind,ind+1], ...
        stiffness(j(:))*[-1 1],m,ngrid);

      [i,j] = meshgrid(1:(nx-1),1:ny);
      ind = j(:) + ny*(i(:)-1);
      m = (nx-1)*ny;
      stiffness = 1./(dx/params.xscale);
      Areg = [Areg;sparse(repmat((1:m)',1,2),[ind,ind+ny], ...
        stiffness(i(:))*[-1 1],m,ngrid)];

      [i,j] = meshgrid(1:(nx-1),1:(ny-1));
      ind = j(:) + ny*(i(:)-1);
      m = (nx-1)*(ny-1);
      stiffness = 1./sqrt((dx(i(:))/params.xscale).^2 + ...
        (dy(j(:))/params.yscale).^2);

      Areg = [Areg;sparse(repmat((1:m)',1,2),[ind,ind+ny+1], ...
        stiffness*[-1 1],m,ngrid)];

      Areg = [Areg;sparse(repmat((1:m)',1,2),[ind+1,ind+ny], ...
        stiffness*[-1 1],m,ngrid)];

    case {'diffusion' 'laplacian'}
      % thermal diffusion using Laplacian (del^2)
      [i,j] = meshgrid(1:nx,2:(ny-1));
      ind = j(:) + ny*(i(:)-1);
      dy1 = dy(j(:)-1)/params.yscale;
      dy2 = dy(j(:))/params.yscale;

      Areg = sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ...
        [-2./(dy1.*(dy1+dy2)), 2./(dy1.*dy2), ...
        -2./(dy2.*(dy1+dy2))],ngrid,ngrid);

      [i,j] = meshgrid(2:(nx-1),1:ny);
      ind = j(:) + ny*(i(:)-1);
      dx1 = dx(i(:)-1)/params.xscale;
      dx2 = dx(i(:))/params.xscale;

      Areg = Areg + sparse(repmat(ind,1,3),[ind-ny,ind,ind+ny], ...
        [-2./(dx1.*(dx1+dx2)), 2./(dx1.*dx2), ...
        -2./(dx2.*(dx1+dx2))],ngrid,ngrid);

    case 'gradient'
      % Subtly different from the Laplacian. A point for future
      % enhancement is to do it better for the triangle interpolation
      % case.
      [i,j] = meshgrid(1:nx,2:(ny-1));
      ind = j(:) + ny*(i(:)-1);
      dy1 = dy(j(:)-1)/params.yscale;
      dy2 = dy(j(:))/params.yscale;

      Areg = sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ...
        [-2./(dy1.*(dy1+dy2)), 2./(dy1.*dy2), ...
        -2./(dy2.*(dy1+dy2))],ngrid,ngrid);

      [i,j] = meshgrid(2:(nx-1),1:ny);
      ind = j(:) + ny*(i(:)-1);
      dx1 = dx(i(:)-1)/params.xscale;
      dx2 = dx(i(:))/params.xscale;

      Areg = [Areg;sparse(repmat(ind,1,3),[ind-ny,ind,ind+ny], ...
        [-2./(dx1.*(dx1+dx2)), 2./(dx1.*dx2), ...
        -2./(dx2.*(dx1+dx2))],ngrid,ngrid)];

  end
  nreg = size(Areg,1);

  % Append the regularizer to the interpolation equations,
  % scaling the problem first. Use the 1-norm for speed.
  NA = norm(A,1);
  NR = norm(Areg,1);
  A = [A;Areg*(params.smoothness*NA/NR)];
  rhs = [rhs;zeros(nreg,1)];
  % do we have a mask to apply?
  if params.maskflag
    unmasked = find(params.mask);
  end
  % solve the full system, with regularizer attached
  switch params.solver
    case {'\' 'backslash'}
      if params.maskflag
        % there is a mask to use
        % permute for minimum fill in for R (in the QR)
        p = colamd(A(:,unmasked));
        zgrid=nan(ny,nx);
        zgrid(unmasked(p)) = A(:,unmasked(p))\rhs;
      else
        % permute for minimum fill in for R (in the QR)
        p = colamd(A);
        zgrid=zeros(ny,nx);
        zgrid(p) = A(:,p)\rhs;
      end

    case 'normal'
      % The normal equations, solved with \. Can be fast
      % for huge numbers of data points.
      if params.maskflag
        % there is a mask to use
        % Permute for minimum fill-in for \ (in chol)
        APA = A(:,unmasked)'*A(:,unmasked);
        p = symamd(APA);
        zgrid=nan(ny,nx);
        zgrid(unmasked(p)) = APA(p,p)\(A(:,unmasked(p))'*rhs);
      else
        % Permute for minimum fill-in for \ (in chol)
        APA = A'*A;
        p = symamd(APA);
        zgrid=zeros(ny,nx);
        zgrid(p) = APA(p,p)\(A(:,p)'*rhs);
      end

    case 'symmlq'
      % iterative solver - symmlq - requires a symmetric matrix,
      % so use it to solve the normal equations. No preconditioner.
      tol = abs(max(z)-min(z))*1.e-13;
      if params.maskflag
        % there is a mask to use
        zgrid=nan(ny,nx);
        [zgrid(unmasked),flag] = symmlq(A(:,unmasked)'*A(:,unmasked), ...
          A(:,unmasked)'*rhs,tol,params.maxiter);
      else
        [zgrid,flag] = symmlq(A'*A,A'*rhs,tol,params.maxiter);
        zgrid = reshape(zgrid,ny,nx);
      end
      % display a warning if convergence problems
      switch flag
        case 0
          % no problems with convergence
        case 1
          % SYMMLQ iterated MAXIT times but did not converge.
          warning(['Symmlq performed ',num2str(params.maxiter), ...
            ' iterations but did not converge.'])
        case 3
          % SYMMLQ stagnated, successive iterates were the same
          warning 'Symmlq stagnated without apparent convergence.'
        otherwise
          warning(['One of the scalar quantities calculated in',...
            ' symmlq was too small or too large to continue computing.'])
      end

    case 'lsqr'
      % iterative solver - lsqr. No preconditioner here.
      tol = abs(max(z)-min(z))*1.e-13;
      if params.maskflag
        % there is a mask to use
        zgrid=nan(ny,nx);
        [zgrid(unmasked),flag] = lsqr(A(:,unmasked),rhs,tol,params.maxiter);
      else
        [zgrid,flag] = lsqr(A,rhs,tol,params.maxiter);
        zgrid = reshape(zgrid,ny,nx);
      end

      % display a warning if convergence problems
      switch flag
        case 0
          % no problems with convergence
        case 1
          % lsqr iterated MAXIT times but did not converge.
          warning(['Lsqr performed ',num2str(params.maxiter), ...
            ' iterations but did not converge.'])
        case 3
          % lsqr stagnated, successive iterates were the same
          warning 'Lsqr stagnated without apparent convergence.'
        case 4
          warning(['One of the scalar quantities calculated in',...
            ' LSQR was too small or too large to continue computing.'])
      end

  end

end  % if params.tilesize...

% only generate xgrid and ygrid if requested.
if nargout>1
  [xgrid,ygrid]=meshgrid(xnodes,ynodes);
end

% ============================================
% End of main function - gridfit
% ============================================

% ============================================
% subfunction - parse_pv_pairs
% ============================================
function params=parse_pv_pairs(params,pv_pairs)
% parse_pv_pairs: parses sets of property value pairs, allows defaults
% usage: params=parse_pv_pairs(default_params,pv_pairs)
%
% arguments: (input)
%  default_params - structure, with one field for every potential
%             property/value pair. Each field will contain the default
%             value for that property. If no default is supplied for a
%             given property, then that field must be empty.
%
%  pv_array - cell array of property/value pairs.
%             Case is ignored when comparing properties to the list
%             of field names. Also, any unambiguous shortening of a
%             field/property name is allowed.
%
% arguments: (output)
%  params   - parameter struct that reflects any updated property/value
%             pairs in the pv_array.
%
% Example usage:
% First, set default values for the parameters. Assume we
% have four parameters that we wish to use optionally in
% the function examplefun.
%
%  - 'viscosity', which will have a default value of 1
%  - 'volume', which will default to 1
%  - 'pie' - which will have default value 3.141592653589793
%  - 'description' - a text field, left empty by default
%
% The first argument to examplefun is one which will always be
% supplied.
%
%   function examplefun(dummyarg1,varargin)
%   params.Viscosity = 1;
%   params.Volume = 1;
%   params.Pie = 3.141592653589793
%
%   params.Description = '';
%   params=parse_pv_pairs(params,varargin);
%   params
%
% Use examplefun, overriding the defaults for 'pie', 'viscosity'
% and 'description'. The 'volume' parameter is left at its default.
%
%   examplefun(rand(10),'vis',10,'pie',3,'Description','Hello world')
%
% params =
%     Viscosity: 10
%        Volume: 1
%           Pie: 3
%   Description: 'Hello world'
%
% Note that capitalization was ignored, and the property 'viscosity'
% was truncated as supplied. Also note that the order the pairs were
% supplied was arbitrary.

npv = length(pv_pairs);
n = npv/2;

if n~=floor(n)
  error 'Property/value pairs must come in PAIRS.'
end
if n<=0
  % just return the defaults
  return
end

if ~isstruct(params)
  error 'No structure for defaults was supplied'
end

% there was at least one pv pair. process any supplied
propnames = fieldnames(params);
lpropnames = lower(propnames);
for i=1:n
  p_i = lower(pv_pairs{2*i-1});
  v_i = pv_pairs{2*i};

  ind = strmatch(p_i,lpropnames,'exact');
  if isempty(ind)
    ind = find(strncmp(p_i,lpropnames,length(p_i)));
    if isempty(ind)
      error(['No matching property found for: ',pv_pairs{2*i-1}])
    elseif length(ind)>1
      error(['Ambiguous property name: ',pv_pairs{2*i-1}])
    end
  end
  p_i = propnames{ind};

  % override the corresponding default in params
  params = setfield(params,p_i,v_i); %#ok

end

% ============================================
% subfunction - check_params
% ============================================
function params = check_params(params)

% check the parameters for acceptability
% smoothness == 1 by default
if isempty(params.smoothness)
  params.smoothness = 1;
else
  if (length(params.smoothness)>1) || (params.smoothness<=0)
    error 'Smoothness must be scalar, real, finite, and positive.'
  end
end

% regularizer  - must be one of 4 options - the second and
% third are actually synonyms.
valid = {'springs', 'diffusion', 'laplacian', 'gradient'};
if isempty(params.regularizer)
  params.regularizer = 'diffusion';
end
ind = find(strncmpi(params.regularizer,valid,length(params.regularizer)));
if (length(ind)==1)
  params.regularizer = valid{ind};
else
  error(['Invalid regularization method: ',params.regularizer])
end

% interp must be one of:
%    'bilinear', 'nearest', or 'triangle'
% but accept any shortening thereof.
valid = {'bilinear', 'nearest', 'triangle'};
if isempty(params.interp)
  params.interp = 'triangle';
end
ind = find(strncmpi(params.interp,valid,length(params.interp)));
if (length(ind)==1)
  params.interp = valid{ind};
else
  error(['Invalid interpolation method: ',params.interp])
end

% solver must be one of:
%    'backslash', '\', 'symmlq', 'lsqr', or 'normal'
% but accept any shortening thereof.
valid = {'backslash', '\', 'symmlq', 'lsqr', 'normal'};
if isempty(params.solver)
  params.solver = '\';
end
ind = find(strncmpi(params.solver,valid,length(params.solver)));
if (length(ind)==1)
  params.solver = valid{ind};
else
  error(['Invalid solver option: ',params.solver])
end

% extend must be one of:
%    'never', 'warning', 'always'
% but accept any shortening thereof.
valid = {'never', 'warning', 'always'};
if isempty(params.extend)
  params.extend = 'warning';
end
ind = find(strncmpi(params.extend,valid,length(params.extend)));
if (length(ind)==1)
  params.extend = valid{ind};
else
  error(['Invalid extend option: ',params.extend])
end

% tilesize == inf by default
if isempty(params.tilesize)
  params.tilesize = inf;
elseif (length(params.tilesize)>1) || (params.tilesize<3)
  error 'Tilesize must be scalar and > 0.'
end

% overlap == 0.20 by default
if isempty(params.overlap)
  params.overlap = 0.20;
elseif (length(params.overlap)>1) || (params.overlap<0) || (params.overlap>0.5)
  error 'Overlap must be scalar and 0 < overlap < 1.'
end

% ============================================
% subfunction - tiled_gridfit
% ============================================
function zgrid=tiled_gridfit(x,y,z,xnodes,ynodes,params)
% tiled_gridfit: a tiled version of gridfit, continuous across tile boundaries
% usage: [zgrid,xgrid,ygrid]=tiled_gridfit(x,y,z,xnodes,ynodes,params)
%
% Tiled_gridfit is used when the total grid is far too large
% to model using a single call to gridfit. While gridfit may take
% only a second or so to build a 100x100 grid, a 2000x2000 grid
% will probably not run at all due to memory problems.
%
% Tiles in the grid with insufficient data (<4 points) will be
% filled with NaNs. Avoid use of too small tiles, especially
% if your data has holes in it that may encompass an entire tile.
%
% A mask may also be applied, in which case tiled_gridfit will
% subdivide the mask into tiles. Note that any boolean mask
% provided is assumed to be the size of the complete grid.
%
% Tiled_gridfit may not be fast on huge grids, but it should run
% as long as you use a reasonable tilesize. 8-)

% Note that we have already verified all parameters in check_params

% Matrix elements in a square tile
tilesize = params.tilesize;
% Size of overlap in terms of matrix elements. Overlaps
% of purely zero cause problems, so force at least two
% elements to overlap.
overlap = max(2,floor(tilesize*params.overlap));

% reset the tilesize for each particular tile to be inf, so
% we will never see a recursive call to tiled_gridfit
Tparams = params;
Tparams.tilesize = inf;

nx = length(xnodes);
ny = length(ynodes);
zgrid = zeros(ny,nx);

% linear ramp for the bilinear interpolation
rampfun = inline('(t-t(1))/(t(end)-t(1))','t');

% loop over each tile in the grid
h = mrvWaitbar(0,'Relax and have a cup of JAVA. Its my treat.');
warncount = 0;
xtind = 1:min(nx,tilesize);
while ~isempty(xtind) && (xtind(1)<=nx)

  xinterp = ones(1,length(xtind));
  if (xtind(1) ~= 1)
    xinterp(1:overlap) = rampfun(xnodes(xtind(1:overlap)));
  end
  if (xtind(end) ~= nx)
    xinterp((end-overlap+1):end) = 1-rampfun(xnodes(xtind((end-overlap+1):end)));
  end

  ytind = 1:min(ny,tilesize);
  while ~isempty(ytind) && (ytind(1)<=ny)
    % update the mrvWaitbar
    mrvWaitbar((xtind(end)-tilesize)/nx + tilesize*ytind(end)/ny/nx)

    yinterp = ones(length(ytind),1);
    if (ytind(1) ~= 1)
      yinterp(1:overlap) = rampfun(ynodes(ytind(1:overlap)));
    end
    if (ytind(end) ~= ny)
      yinterp((end-overlap+1):end) = 1-rampfun(ynodes(ytind((end-overlap+1):end)));
    end

    % was a mask supplied?
    if ~isempty(params.mask)
      submask = params.mask(ytind,xtind);
      Tparams.mask = submask;
    end

    % extract data that lies in this grid tile
    k = (x>=xnodes(xtind(1))) & (x<=xnodes(xtind(end))) & ...
        (y>=ynodes(ytind(1))) & (y<=ynodes(ytind(end)));
    k = find(k);

    if length(k)<4
      if warncount == 0
        warning 'A tile was too underpopulated to model. Filled with NaNs.'
      end
      warncount = warncount + 1;

      % fill this part of the grid with NaNs
      zgrid(ytind,xtind) = NaN;

    else
      % build this tile
      zgtile = gridfit(x(k),y(k),z(k),xnodes(xtind),ynodes(ytind),Tparams);

      % bilinear interpolation (using an outer product)
      interp_coef = yinterp*xinterp;

      % accumulate the tile into the complete grid
      zgrid(ytind,xtind) = zgrid(ytind,xtind) + zgtile.*interp_coef;

    end

    % step to the next tile in y
    if ytind(end)<ny
      ytind = ytind + tilesize - overlap;
      % are we within overlap elements of the edge of the grid?
      if (ytind(end)+max(3,overlap))>=ny
        % extend this tile to the edge
        ytind = ytind(1):ny;
      end
    else
      ytind = ny+1;
    end

  end % while loop over y

  % step to the next tile in x
  if xtind(end)<nx
    xtind = xtind + tilesize - overlap;
    % are we within overlap elements of the edge of the grid?
    if (xtind(end)+max(3,overlap))>=nx
      % extend this tile to the edge
      xtind = xtind(1):nx;
    end
  else
    xtind = nx+1;
  end

end % while loop over x

% close down the mrvWaitbar
close(h)

if warncount>0
  warning([num2str(warncount),' tiles were underpopulated & filled with NaNs'])
end

author: vistalab

repo: vistasoft


text:
classdef NMPSO < ALGORITHM
% <multi/many> <real/integer>
% Novel multi-objective particle swarm optimization

%------------------------------- Reference --------------------------------
% Q. Lin, S. Liu, Q. Zhu, C. Tang, R. Song, J. Chen, C. A. Coello Coello,
% K. Wong, and J. Zhang, Particle swarm optimization with a balanceable
% fitness estimation for many-objective optimization problems, IEEE
% Transactions on Evolutionary Computation, 2018, 22(1): 32-46.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();
            Pbest      = Population;
            Archive    = UpdateArchive(Population(NDSort(Population.objs,1)==1),[],Problem.N);

            %% Optimization
            while Algorithm.NotTerminated(Archive)
                Population = Operator(Problem,Population,Pbest,Archive(randi(ceil(length(Archive)/10),1,Problem.N)));
                Pbest      = UpdatePbest(Pbest,Population);
                Archive    = UpdateArchive(Archive,Population,Problem.N);
                S          = OperatorGAhalf(Problem,Archive([1:length(Archive),randi(ceil(length(Archive)/2),1,length(Archive))]));
                Archive    = UpdateArchive(Archive,S,Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [gd,nlen]=ptpfundual(w,a,M,L,varargin)
%PTPFUNDUAL Sampled, periodized dual TP function of finite type
%   Usage: gd=ptpfundual(w,a,M,L)
%          gd=ptpfundual({w,width},a,M,L)
%          gd=ptpfundual(...,inc)
%          [gd,nlen]=ptpfundual(...)
%
%   Input parameters:
%         w      : Vector of reciprocals $w_j=1/\delta_j$ in Fourier representation of *g*
%         width  : Integer stretching factor for the essential support
%         a      : Length of time shift.
%         M      : Number of channels.
%         L      : Window length.
%         inc    : Extension parameter
%
%   Output parameters:
%         gd    : The periodized totally positive function dual window.
%         nlen  : Number of non-zero elements in gd.
%
%   `ptpfundual(w,a,M,L)` computes samples of a dual window of totally
%   positive function of finite type >=2 with weights *w*. Please see
%   |ptpfun| for definition of totally positive functions.
%   The lattice parameters $a,M$ must satisfy $M > a$ to ensure the
%   system is a frame.
%
%   `ptpfundual({w,width},a,M,L)` works as above but in addition the *width*
%   parameter determines the integer stretching factor of the original TP
%   function. For explanation see help of |ptpfun|.
%
%   `ptpfundual(...,inc)` or `ptpfundual(...,'inc',inc)` works as above,
%   but integer *inc* denotes number of additional columns to compute window
%   function *gd*. 'inc'-many are added at each side. It should be smaller
%   than 100 to have comfortable execution-time. The higher the number the
%   closer *gd* is to the canonical dual window.
%   The default value is 10.
%
%   `[gd,nlen]=ptpfundual(...)` as *gd* might have a compact support,
%   *nlen* contains a number of non-zero elements in *gd*. This is the case
%   when *gd* is symmetric. If *gd* is not symmetric, *nlen* is extended
%   to twice the length of the longer tail.
%
%   If $nlen = L$, *gd* has a 'full' support meaning it is a periodization
%   of a dual TP function.
%
%   If $nlen < L$, additional zeros can be removed by calling
%   `gd=middlepad(gd,nlen)`.
%
%   Examples:
%   ---------
%
%   The following example compares dual windows computed using 2 different
%   approaches.:::
%
%     w = [-3,-1,1,3];a = 25; M = 31; inc = 10;
%     L = 1e6; L = dgtlength(L,a,M);
%     width = M;
%
%     % Create the window
%     g = ptpfun(L,w,width);
%
%     % Compute a dual window using pebfundual
%     tic
%     [gd,nlen] = ptpfundual({w,width},a,M,L,inc);
%     ttpfundual=toc;
%
%     % We know that gd has only nlen nonzero samples, lets shrink it.
%     gd = middlepad(gd,nlen);
%
%     % Compute the canonical window using gabdual
%     tic
%     gdLTFAT = gabdual(g,a,M,L);
%     tgabdual=toc;
%
%     fprintf('PTPFUNDUAL elapsed time %f s\n',ttpfundual);
%     fprintf('GABDUAL elapsed time    %f s\n',tgabdual);
%
%     % Test on random signal
%     f = randn(L,1);
%
%     fr = idgt(dgt(f,g,a,M),gd,a,numel(f));
%     fprintf('Reconstruction error PTPFUNDUAL: %e\n',norm(f-fr)/norm(f));
%
%     fr = idgt(dgt(f,g,a,M),gdLTFAT,a,numel(f));
%     fprintf('Reconstruction error GABDUAL:    %e\n',norm(f-fr)/norm(f));
%
%   See also: dgt, idgt, ptpfun
%
%   References: grst13 kl12 bagrst14 klst14

% AUTHORS: Joachim Stoeckler, Tobias Kloos, 2012-2014

complainif_notenoughargs(nargin,4,upper(mfilename));
complainif_notposint(L,'L',upper(mfilename));
complainif_notposint(a,'a',upper(mfilename));
complainif_notposint(M,'M',upper(mfilename));

% Check lattice
if M<=a
    error('%s: Lattice parameters must satisfy M>a.',upper(mfilename));
end

% Check w
if iscell(w)
    if numel(w)~=2
        error('%s: w must be a 2 element cell array.',upper(mfilename));
    end
    width = w{2};
    w = w{1};
    complainif_notposint(width,'width',upper(mfilename));
else
    width = floor(sqrt(L));
end

if isempty(w) || ~isnumeric(w) || numel(w)<2
    error(['%s: w must be a nonempty numeric vector with at least',...
    ' 2 elements.'], upper(mfilename));
end

if any(w==0)
    error('%s: All weights w must be nonzero.', upper(mfilename));
    % TO DO: Also add a warning if w is very small or big?
end

% Define initial value for flags and key/value pairs.
%definput.import={'setnorm'};
definput.keyvals.inc = 10;
%definput.flags.scale = {'nomatchscale','matchscale'};
[flags,~,inc]=ltfatarghelper({'inc'},definput,varargin);
complainif_notnonnegint(inc,'inc',upper(mfilename));

% TP functions are scale invariant so we do scaling directly on w.
wloc = w/width;
% Converting a, M to alpha, beta
alpha = a;
beta = 1/M;

% check alpha beta
if (alpha<=0) || (beta<=0)
    error('lattice parameters alpha, beta must be positive')
end
if (width*beta > 10)
    warning('width/M should be smaller than 10: numerical instability may occur')
end

% compute m n and check that a has nonzero entries
if all(wloc<0)
    wloc = -wloc;
    case0 = 2;
else
    case0 = 0;
end
m = length(find(wloc>0));
n = length(find(wloc<0));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% preparations specially for computation of gd
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r = floor(1/(1-alpha*beta)+eps);

% check special cases according to n and m
if n == 0
    if m >= 2
        k2 = (m-1)*(r+1)-1;
        k1 = -k2;
        if case0 == 0
            case0 = 1;
        end
    end
elseif n == 1
    N = m*(r+1)+1;       % minimal column size
    k1 = -m*(r+1)+1;     % column index k1 from the paper
    k2 = k1+N-1;         % column index k2 from the paper
else
    N = (m+n-1)*(r+1);   % minimal column size
    k1 = -m*(r+1)+1;     % column index k1 from the paper
    k2 = k1+N-1;         % column index k2 from the paper
end

k1 = k1-inc;
k2 = k2+inc;

% minimal values for x and y
varl = floor((k1+m-1)/(alpha*beta))-1;
varr = ceil((k2-n+1)/(alpha*beta))+1;
x = varl*alpha:alpha:varr*alpha;
i0 = abs(varl)+1; % index of "central" row of P(x)
y = (k1-1)/beta:(1/beta):(k2+1)/beta;
k0 = abs(k1-1)+1; % index of "central" column of P(x)

[yy,xx] = meshgrid(y,x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% discretization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 0:(a-1); % for stepping through the interval [0,alpha)
tt = varl*alpha:varr*alpha; % choose same stepsize for t and tt
% left and right bounds large enough for the support of gamma
tt0 = abs(varl*a)+1; % index for tt == 0
gd = zeros(1,length(tt)); % dual window

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computation of gamma
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k1 = k1+k0;
k2 = k2+k0;

for k=1:length(t)
    % step through the interval [0,alpha)
    x0 = t(k); % compute dual window at points (x+j*alpha)

    % row indices for rectangular P0
    i1 = floor((k1-k0+m-1)/alpha/beta-x0/alpha)+1+i0;
    i2 = ceil((k2-k0-n+1)/alpha/beta-x0/alpha)-1+i0;

    % Computation of P0(x0)
    % z0 is the matrix of the abscissa x0+j*alpha-k/beta, j=i1:i2, k=k1:k2,
    % z1 puts all these abscissae into a row vector.
    % The computation of g(z1) is done as described above for the
    % vector tt.
    z0 = x0+xx(:,k1:k2)-yy(:,k1:k2);
    z1 = z0(:)';

    Y = zeros((n+m)-1,length(z1));
    for q = 1:(n+m)-1
        if wloc(q) == wloc(q+1)
            Y(q,:) = abs(wloc(q))^2*abs((z1.*((wloc(q)*z1)>=0))).*exp(-wloc(q)*(z1.*((wloc(q)*z1)>=0))).*((wloc(q)*z1)>=0);
        else
            if wloc(q)*wloc(q+1) < 0
                Y(q,:) = abs(wloc(q)*wloc(q+1))/(abs(wloc(q))+abs(wloc(q+1)))*(exp(-wloc(q)*(z1.*((wloc(q)*z1)>=0))).*((wloc(q)*z1)>=0) + ...
                exp(-wloc(q+1)*(z1.*((wloc(q+1)*z1)>=0))).*((wloc(q+1)*z1)>0));
            else
                Y(q,:) = wloc(q)*wloc(q+1)/(abs(wloc(q+1))-abs(wloc(q)))*(exp(-wloc(q)*(z1.*((wloc(q)*z1)>=0)))-exp(-wloc(q+1)*(z1.*((wloc(q)*z1)>=0)))).*((wloc(q)*z1)>=0);
            end
        end
    end

    for q = 2:(n+m)-1
        for j = 1:(n+m)-q
            if wloc(j) == wloc(j+q)
                Y(j,:) = Y(j,:).*abs(z1)/q*abs(wloc(j));
            else
                Y(j,:) = (wloc(j)*Y(j+1,:)-wloc(j+q)*Y(j,:))/(wloc(j)-wloc(j+q));
            end
        end
    end

    if (n+m) == 1
        A0 = abs(wloc)*exp(-wloc*(z1.*((wloc*z1)>=0))).*((wloc*z1)>=0);
    else
        A0 = Y(1,:);
    end

    A0 = reshape(A0,size(z0))*sqrt(width);%*L^(1/4);
    P0 = A0(i1:i2,:);

    % computation of pseudo-inverse matrix of P0
    P0inv = pinv(P0);
    gd(k-1+tt0-a*(i0-i1):a:k-1+tt0+a*(i2-i0)) = beta*P0inv(k0-k1+1,:); % row index k0-k1a+1
    % points to the "j=0" row of P0inv
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% periodization of gamma
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nlen = length(gd);
nr = ceil(nlen/L);
v = zeros(1,nr*L);
v(1:length(gd)) = gd;
v = [v(tt0:end),v(1:tt0-1)];

gd = sum(reshape(v,L,nr),2);
gd = gd(:);
if case0 == 2
    gd =  flipud(gd);
    gd = [gd(end);gd(1:end-1)];
end

% Determine nlen
if nlen<L
    negsupp = tt0-1;
    possupp = nlen-tt0; % excluding zero pos.
    nlen = 2*max([negsupp,possupp])+1;
end
nlen = min([L,nlen]);

% if flags.do_matchscale
%    g = ptpfun(L,w,width,flags.norm);
%    [scal,err] = gabdualnorm(g,gd,a,M,L);
%     assert(err<1e-10,sprintf(['%s: Assertion failed. This is not a valid ',...
%                               ' dual window.'],upper(mfilename)));
%    gd = gd/scal;
% else
%    gd = setnorm(gd,flags.norm);
% end

gd = gd(:);

end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%fsqfindewq(lambda,mu,c,m)
%   This function finds the average waiting time in the queue
%   (ewq) for a machine repair problem (finite source queue)

function out = fsqfindewq(lambda,mu,c,m)

ld = fsqfindld(lambda,mu,c,m);
lu = fsqfindlu(lambda,mu,c,m);

ewq = ld/(lambda*lu);

out = ewq;

author: Sable

repo: mcbench-benchmarks


text:
function results = vl_test_imdisttf(varargin)
% VL_TEST_DISTTF
vl_test_init ;

function test_basic()
for conv = {@single, @double}
  conv = conv{1} ;

  I = conv([0 0 0 ; 0 -2 0 ; 0 0 0]) ;
  D = vl_imdisttf(I);
  assert(isequal(D, conv(- [0 1 0 ; 1 2 1 ; 0 1 0]))) ;

  I(2,2) = -3 ;
  [D,map] = vl_imdisttf(I) ;
  assert(isequal(D, conv(-1 - [0 1 0 ; 1 2 1 ; 0 1 0]))) ;
  assert(isequal(map, 5 * ones(3))) ;
end

function test_1x1()
assert(isequal(1, vl_imdisttf(1))) ;

function test_rand()
I = rand(13,31) ;
for t=1:4
  param = [rand randn rand randn] ;
  [D0,map0] = imdisttf_equiv(I,param) ;
  [D,map] = vl_imdisttf(I,param) ;
  vl_assert_almost_equal(D,D0,1e-10)
  assert(isequal(map,map0)) ;
end

function test_param()
I = zeros(3,4) ;
I(1,1) = -1 ;

[D,map] = vl_imdisttf(I,[1 0 1 0]);
assert(isequal(-[1 0 0 0 ;
                 0 0 0 0 ;
                 0 0 0 0 ;], D)) ;

D0 = -[1 .9 .6 .1 ;
       0 0 0 0 ;
       0 0 0 0 ;] ;
[D,map] = vl_imdisttf(I,[.1 0 1 0]);
vl_assert_almost_equal(D,D0,1e-10);

D0 = -[1  .9 .6 .1 ;
       .9 .8 .5  0 ;
       .6 .5 .2  0 ;]  ;
[D,map] = vl_imdisttf(I,[.1 0 .1 0]);
vl_assert_almost_equal(D,D0,1e-10);

D0 = -[.9  1  .9  .6 ;
       .8 .9  .8  .5 ;
       .5 .6  .5  .2 ; ] ;
[D,map] = vl_imdisttf(I,[.1 1 .1 0]);
vl_assert_almost_equal(D,D0,1e-10);

function test_special()
I = rand(13,31) -.5 ;
D = vl_imdisttf(I, [0 0 1e5 0]) ;
vl_assert_almost_equal(D(:,1),min(I,[],2),1e-10);
D = vl_imdisttf(I, [1e5 0 0 0]) ;
vl_assert_almost_equal(D(1,:),min(I,[],1),1e-10);

function [D,map]=imdisttf_equiv(I,param)
D = inf + zeros(size(I)) ;
map = zeros(size(I)) ;
ur = 1:size(D,2) ;
vr = 1:size(D,1) ;
[u,v] = meshgrid(ur,vr) ;
for v_=vr
  for u_=ur
    E = I(v_,u_) + ...
        param(1) * (u - u_ - param(2)).^2 + ...
        param(3) * (v - v_ - param(4)).^2 ;
    map(E < D) = sub2ind(size(I),v_,u_) ;
    D = min(D,E) ;
  end
end

author: Cloud-CV

repo: object-proposals


text:
function [F,dF,G,H,varF,dH,varGss,varG,varH,I_sk,J_sjk] = negelcbo_vbmc(theta,beta,vp,gp,Ns,compute_grad,compute_var,altent_flag,thetabnd,entropy_alpha)
%NEGELCBO_VBMC Negative evidence lower confidence bound objective
%
% Note that THETA is a vector of *transformed* variational parameters:
% [MU_1,...,MU_K,log(SIGMA)] or
% [MU_1,...,MU_K,log(SIGMA),log(LAMBDA)] or
% [MU_1,...,MU_K,log(SIGMA),log(LAMBDA),log(W)]

if nargin < 5 || isempty(Ns); Ns = 0; end
if nargin < 6 || isempty(compute_grad); compute_grad = nargout > 1; end
if nargin < 7; compute_var = []; end
if nargin < 8 || isempty(altent_flag); altent_flag = false; end
if nargin < 9; thetabnd = []; end
if nargin < 10 || isempty(entropy_alpha); entropy_alpha = 0; end
if isempty(beta) || ~isfinite(beta); beta = 0; end
if isempty(compute_var); compute_var = beta ~=0 || nargout > 4; end
separate_K = nargout > 9; % Return expected log joint per component

% altent_flag and entropy_alpha are unused (kept here for retrocompatibility)

if compute_grad && beta ~= 0 && compute_var ~= 2
    error('negelcbo_vbmc:vargrad', ...
        'Computation of the gradient of ELBO with full variance not supported.');
end

D = vp.D;
K = vp.K;

avg_flag = 1;       % Average over multiple GP hyperparameters if provided
jacobian_flag = 1;  % Variational parameters are transformed

% Reformat variational parameters from THETA
if vp.optimize_mu
    vp.mu(:,:) = reshape(theta(1:D*K),[D,K]);
    idx_start = D*K;
else
    idx_start = 0;
end
if vp.optimize_sigma
    vp.sigma(1,:) = exp(theta(idx_start+(1:K)));
    idx_start = idx_start + K;
end
if vp.optimize_lambda; vp.lambda(:,1) = exp(theta(idx_start+(1:D))); end
if vp.optimize_weights
    vp.eta(1,:) = theta(end-K+1:end);
    vp.w(1,:) = exp(vp.eta);
    vp.w = vp.w/sum(vp.w);
end

% Which gradients should be computed, if any?
grad_flags = compute_grad*[vp.optimize_mu,vp.optimize_sigma,vp.optimize_lambda,vp.optimize_weights];

% Only weight optimization?
onlyweights_flag = vp.optimize_weights && ~vp.optimize_mu && ~vp.optimize_sigma && ~vp.optimize_lambda;

if separate_K
    if compute_grad
        error('Computing the gradient of variational parameters and requesting per-component results at the same time.');
    end

    if onlyweights_flag
        if compute_var
            [G,~,varG,~,~,I_sk,J_sjk] = gplogjoint_weights(vp,0,avg_flag,jacobian_flag,compute_var);
        else
            [G,dG,~,~,~,I_sk] = gplogjoint_weights(vp,compute_grad,avg_flag,jacobian_flag,0);
            J_sjk = [];
        end
        varGss = NaN;
    else
        if compute_var
            [G,~,varG,~,varGss,I_sk,J_sjk] = gplogjoint(vp,gp,grad_flags,avg_flag,jacobian_flag,compute_var);
        else
            [G,dG,~,~,~,I_sk] = gplogjoint(vp,gp,grad_flags,avg_flag,jacobian_flag,0);
            varGss = 0; varG = 0; J_sjk = [];
        end
    end
else
    if onlyweights_flag
        if compute_var
            if compute_grad
                [G,dG,varG,dvarG] = gplogjoint_weights(vp,1,avg_flag,jacobian_flag,compute_var);
            else
                [G,~,varG] = gplogjoint_weights(vp,0,avg_flag,jacobian_flag,compute_var);
            end
        else
            [G,dG] = gplogjoint_weights(vp,compute_grad,avg_flag,jacobian_flag,0);
        end
        varGss = NaN;
    else
        if compute_var
            if compute_grad
                [G,dG,varG,dvarG,varGss] = gplogjoint(vp,gp,grad_flags,avg_flag,jacobian_flag,compute_var);
            else
                [G,~,varG,~,varGss] = gplogjoint(vp,gp,grad_flags,avg_flag,jacobian_flag,compute_var);
            end
        else
            [G,dG] = gplogjoint(vp,gp,grad_flags,avg_flag,jacobian_flag,0);
            varGss = 0; varG = 0;
        end
    end
end

% Entropy term
if Ns > 0
    % Monte Carlo approximation
    [H,dH] = entmc_vbmc(vp,Ns,grad_flags,jacobian_flag);
else
    % Deterministic approximation via lower bound on the entropy
    [H,dH] = entlb_vbmc(vp,grad_flags,jacobian_flag);
end

%H_check = gmment_num(theta,lambda);
%[H - H_check, (H - H_check)/H_check ]

% Negative ELBO and its gradient
F = -G - H;
if compute_grad; dF = -dG - dH; else; dF = []; end

varH = 0;   % For the moment use zero variance for entropy
if compute_var
    varF = varG + varH;
else
    varF = 0;
end

% Negative ELCBO (add confidence bound)
if beta ~= 0; F = F + beta*sqrt(varF); end
if beta ~= 0 && compute_grad
    dF = dF + 0.5*beta*dvarG/sqrt(varF);
end

% Additional loss for variational parameter bound violation (soft bounds)
% and for weight size (if optimizing mixture weights)
% Only done when optimizing the variational parameters, but not when
% computing the EL(C)BO at each iteration
if ~isempty(thetabnd)
    if compute_grad
        [L,dL] = vpbndloss(theta,vp,thetabnd,thetabnd.TolCon);
        dF = dF + dL;
    else
        L = vpbndloss(theta,vp,thetabnd,thetabnd.TolCon);
    end
    F = F + L;

    % Penalty to reduce weight size
    if vp.optimize_weights
        Thresh = thetabnd.WeightThreshold;
        %L = sum(vp.w)*thetabnd.WeightPenalty;
        % L = sum(sqrt(vp.w))*thetabnd.WeightPenalty;
        L = sum(vp.w.*(vp.w<Thresh) + Thresh*(vp.w>=Thresh))*thetabnd.WeightPenalty;
        F = F + L;
        if compute_grad
            %w_grad = thetabnd.WeightPenalty*ones(K,1);
            % w_grad = 0.5./sqrt(vp.w(:))*thetabnd.WeightPenalty;
            w_grad = thetabnd.WeightPenalty.*(vp.w(:)<Thresh);
            eta_sum = sum(exp(vp.eta));
            J_w = bsxfun(@times,-exp(vp.eta)',exp(vp.eta)/eta_sum^2) + diag(exp(vp.eta)/eta_sum);
            w_grad = J_w*w_grad;
            dL = zeros(size(dF));
            dL(end-K+1:end) = w_grad;
            dF = dF + dL;
        end
    end
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function [w0, t] = adw_fan(ob, wi, mask, Phi, varargin)
% function [w0, phi, t] = adw_fan(ob, wi)
% Compute the angular-dependent weighting for 3rd gen fan-beam
% w0(\Phi) = w(x0, y0, \Phi)
%         for fully corrected penalty:
%          = w(s',\beta')*J(s')|\phi'=\Phi + ...
%            w(s',\beta')*J(s')|\phi'=\Phi-pi
% see fessler chapter 3.5.3
%
% in:
%    ob         a struct or object
%               including system parameters, such as pixel-size etc.
%    wi         cov(yi)
%    varargin   default: only compute phi' = Phi
%               'all': compute phi' = Phi + phi' = Phi- pi
% out:
%    w0     angular-dependent weighting
%    t      computation time
%
% Copyright 2005-4-07, Yingying Zhang, The University of Michigan

all = logical(0);
if ~isempty(varargin)
    arg = varargin{1};
    if isempty(arg)
			%	do nothing for empty arguments

	elseif ~ischar(arg)
        error 'unknown non-string argument?'
    elseif streq(arg, 'all')
        all = logical(1);
    else
        error 'unknown string argument?'
    end
end

% default argument
ob.Df = 0; % 3rd generation CT scanner

% image domain pixel locations
% only compute within mask

yflip = (-1); % since the top is positive in y direction %
% x = ([1:ob.nx] - (ob.nx+1)/2)* ob.dx;
% y = ([1:ob.ny] - (ob.ny+1)/2)* ob.dx * yflip;

x = ([0:ob.nx-1] - ob.nx/2)* ob.dx;
y = ([0:ob.ny-1] - ob.ny/2)* ob.dx * yflip;

[x y] = ndgrid(x, y);

rr = sqrt(x.^2 + y.^2);
theta = atan2(y,x); % use it for our purpose since atan will have pi flip

r0 = rr(mask(:));
phi0 = theta(mask(:));

clear rr theta

% dummy variables in freq. domain
phi = Phi;

% if ~isempty(varargin) & isnumeric(varargin{1})
%     phi = varargin{1};
% end

% sinogram domain locations
ob.beta = ([0:ob.na-1]'/ob.na * ob.orbit ...
    + ob.orbit_start) * pi / 180;
ob.s = ([-(ob.nb-1)/2:(ob.nb-1)/2]' ...
			- ob.offset_s) * ob.ds;
ob.gamma = ob.s / ob.dis_src_det;
[ss, beta] = meshgrid(ob.s,ob.beta); % meshgrid and ndgrid reverse

% loop over \phi_i
Dcf = ob.dis_src_det + ob.Df;
ratio_fcf = ob.Df/Dcf;
% w0 = zeros(ob.nx, ob.ny, ob.na);
w0 = zeros([size(r0,1) length(phi)]);

tic
for ia = 1:length(phi)
    % -------------------------
    % for \phi' = \phi
    % -------------------------

    % compute new variables in (3.5.2)
    r_p = r0 .* cos(phi(ia) - phi0); % use ob.beta, not outer_sum(ob.gamma, ob.beta)
%                                         since evaluate at \phi_p = \phi(view angle)
    s_p = Dcf.*(asin(r_p./ob.dis_src_iso) - asin(ratio_fcf.*r_p./ob.dis_src_iso));
    beta_p = phi(ia) - asin(r_p./ob.dis_src_iso);

    % compute Jacobian determinant (2.8.9)
    Jacob = ob.dis_src_det ./ sqrt(ob.dis_src_iso.^2 - r_p.^2);

%     % find corresponding wi: nearest neighbor method
%     s_cen = (ob.nb - 1)/2 + ob.offset_s; % for nb = even % <-----(yy) plus:(to me account for the right index) or (hugo) minus?
%     s_loc = round(s_p./ob.ds + s_cen);
%     s_loc = min(s_loc, ob.nb);
%     s_loc = max(s_loc, 1);
%     beta_loc = round((mod(beta_p,2*pi)-ob.orbit_start)./(2*pi/ob.na)) + 1;
%     beta_loc = min(beta_loc, ob.na);
%     beta_loc = max(beta_loc, 1);
%
%     % find the corresponding index in column-stack wi
%     loc = 1 + (s_loc(:)-1) + (beta_loc(:)-1)*ob.nb;

    % compute w0 at \phi' = \phi
%     w0(:,ia) =wi(loc(:)) .* Jacob;
    si = min(s_p, ob.s(end)); si = max(si, ob.s(1));
    betai = min(beta_p, ob.beta(end)); betai = max(betai, ob.beta(1));
    w0(:,ia) = interp2(ss, beta, wi.', si(:), betai(:)) .* Jacob;

    clear s_loc beta_loc loc
    % -------------------------
    % for \phi' = \phi - pi
    % -------------------------
    if all
        phi_pi = phi(ia) - pi;
        r_p = r0 .* cos(phi_pi - phi0);
        s_p = Dcf.*(asin(r_p./ob.dis_src_iso) - asin(ratio_fcf.*r_p./ob.dis_src_iso));
        beta_p = phi_pi - asin(r_p./ob.dis_src_iso);

% %         Jacob_p again: not needed since same
%         Jacob_p = abs(ob.dis_src_iso*cos(s_p./ob.dis_src_det) - ...
%             ob.source_offset*sin(s_p./ob.dis_src_det))./ob.dis_src_det;
%
%         % find corresponding wi: nearest neighbor method
%         s_loc = round(s_p./ob.ds + s_cen);
%         s_loc = min(s_loc, ob.nb);
%         s_loc = max(s_loc, 1);
%
%         beta_loc = round((mod(beta_p,2*pi)-ob.orbit_start)./(2*pi/ob.na)) + 1;
%         beta_loc = min(beta_loc, ob.na);
%         beta_loc = max(beta_loc, 1);
%
%         % find the corresponding index in column-stack wi
%         loc = 1 + (s_loc(:)-1) + (beta_loc(:)-1)*ob.nb;
%
%         % compute w0 at \phi' = \phi - pi and sum with \phi' = \phi
%         w0(:,ia) = w0(:,ia) + wi(loc(:)).*Jacob;

        si = min(s_p, ob.s(end)); si = max(si, ob.s(1));
        betai = min(beta_p, ob.beta(end)); betai = max(betai, ob.beta(1));
        w0(:,ia) = w0(:,ia) + interp2(ss, beta, wi.', si(:), betai(:)) .* Jacob;
        clear s_loc beta_loc loc
    end

%     % angle_dependent blur
%     b0(:,:,ia) = (ob.ray_spacing./Dcf) .* sqrt(ob.dis_src_iso^2 + r0^2 + 2 .* ob.dis_src_iso .* Dcf .* cos(phi(ia)-phi0));
end

t = toc;
sprintf('angular weighting computation time %.4f', t)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function b = back1(engine, bfuture, f, t)

if t ~= 1
  error('mixed up time stamps')
end
bnet = bnet_from_engine(engine);
ss = bnet.nnodes_per_slice;

int = engine.interface;
D = engine.in_clq; % from J2
C = engine.int_clq1; % from J1
phiD = marginalize_pot(bfuture.clpot{D}, int, engine.maximize);
phiC = marginalize_pot(f.clpot{C}, int, engine.maximize);
ratio = divide_by_pot(phiD, phiC);
f.clpot{C} = multiply_by_pot(f.clpot{C}, ratio);

[b.clpot, seppot] = distribute_evidence(engine.jtree_engine1, f.clpot, f.seppot);
for c=1:length(b.clpot)
  [b.clpot{c}, ll(c)] = normalize_pot(b.clpot{c});
end
b.t = t;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function FinalEx
%{
FinalEx.m

WASDQE for steering. V to switch to 1st person view. More below.

Combine stuff from all the previous examples to make:
1. Import a plane as a patch object.
2. Translate/rotate this through space as per user input. User changes
heading, constant velocity dynamics. Use Ws,AD, QE for pitch,roll,yaw
3. Put in horizon and ground with textures. If the plane COM tries to
penetrate, simply project it back up.
4. Have camera follow the plane in 2 ways, press v to switch between.

NOTE: I wouldn't start here if you want to understand each step. Start with
the other examples and see how this is built.

Matt Sheen, mws262@cornell.edu
%}

close all

fig = figure;

fig.UserData.firstPerson = false; %do we start in 1st person view, or not?
resetPositionOnBounds = true; %keeps the plane in the region with the matlab logo ground.

hold on
fig.Position = [600 600 1500 1200];
%Disable axis viewing, don't allow axes to clip the plane, turn on
%perspective mode.
fig.Children.Visible = 'off';
fig.Children.Clipping = 'off';
fig.Children.Projection = 'perspective';

%% Add the sky - Fancy environment folder ex
skymap = [linspace(1,0.4,100)',linspace(1,0.4,100)',linspace(1,0.99,100)'];

[skyX,skyY,skyZ] = sphere(50);
sky = surf(500000*skyX,500000*skyY,500000*skyZ,'LineStyle','none','FaceColor','interp');
colormap(skymap);

%% Add the plane - Fancy planes folder
fv = stlread('test.stl');
p1 = patch(fv,'FaceColor',       'red', ...
         'EdgeColor',       'none',        ...
         'FaceLighting',    'gouraud',     ...
         'AmbientStrength', 0.15);
vert = p1.Vertices;

material('metal')

%% Add the ground + textures - Fancy environments folder Ex
texture = imread('texture.jpg');

ground = 10000*membrane(1,40)-10000;
groundSurf = surf(linspace(-10000,10000,size(ground,1)),linspace(-10000,10000,size(ground,2)),ground,'LineStyle','none','AmbientStrength',0.3,'DiffuseStrength',0.8,'SpecularStrength',0,'SpecularExponent',10,'SpecularColorReflectance',1);
groundSurf.FaceColor = 'texturemap';
groundSurf.CData = texture;

%add some extra flat ground going off to (basically) infinity.
flatground = surf(linspace(-500000,500000,size(ground,1)),linspace(-500000,500000,size(ground,2)),-10001*ones(size(ground)));
flatground.FaceColor = 'texturemap';
flatground.CData = texture;
flatground.AlphaData = 0.8;

camlight('headlight');

camva(40); %view angle

%% Set keyboard callbacks and flags for movement.
set(fig,'WindowKeyPressFcn',@KeyPress,'WindowKeyReleaseFcn', @KeyRelease);
        fig.UserData.e = false;
         fig.UserData.q = false;
         fig.UserData.a = false;
         fig.UserData.d = false;
         fig.UserData.w = false;
         fig.UserData.s = false;

forwardVec = [1 0 0]'; %Vector of the plane's forward direction in plane frame
rot = eye(3,3); %Initial plane rotation
pos = [-8000,8000,-2000]; %Initial plane position
vel = 1000; %Velocity

hold off
axis([-10000 10000 -10000 10000 -10000 10000])

%% Animation loop:
tic
told = 0;
while(ishandle(fig))
  tnew = toc;

  %Check for user inputs:
  if fig.UserData.e
      rot = rot*angle2dcm(0.05,0,0);
  end
  if fig.UserData.q
      rot = rot*angle2dcm(-0.05,0,0);
  end
  if fig.UserData.s
      rot = rot*angle2dcm(0,-0.05,0);
  end
  if fig.UserData.w
      rot = rot*angle2dcm(0,0.05,0);
  end
  if fig.UserData.a
      rot = rot*angle2dcm(0,0,-0.05);
  end
  if fig.UserData.d
      rot = rot*angle2dcm(0,0,0.05);
  end

  %Update plane's center position.
  pos = vel*(rot*forwardVec*(tnew-told))' + pos;

  %If the plane wants to go under the ground, then bring it back up to the
  %ground surface.
  nearestGroundZ = interp2(groundSurf.XData,groundSurf.YData,groundSurf.ZData,pos(1),pos(2));
  if pos(3)<nearestGroundZ
      pos(3) = nearestGroundZ;
  end

  if resetPositionOnBounds
      % If we leave the ground area in the X direction, then snap the plane
      % back to the other side.
      if pos(1)>max(groundSurf.XData)
          pos(1) = min(groundSurf.XData);
      elseif pos(1)<min(groundSurf.XData)
          pos(1) = max(groundSurf.XData);
      end

      % If we leave the ground area in the y direction, then snap the plane
      % back to the other side.
      if pos(2)>max(groundSurf.YData)
          pos(2) = min(groundSurf.YData);
      elseif pos(2)<min(groundSurf.YData)
          pos(2) = max(groundSurf.YData);
      end
  end

  %Update the plane's vertices using the new position and rotation
  p1.Vertices = (rot*vert')' + repmat(pos,[size(vert,1),1]);

  %Camera updates:
  if fig.UserData.firstPerson %First person view -- follow the plane from slightly behind.
      camupvec = rot*[0 0 1]';
      camup(camupvec);
      campos(pos' - 1000*rot*[1 0 -0.25]');
      camtarget(pos' + 100*rot*[1 0 0]');
  else %Follow the plane from a fixed angle
    campos(pos + [-3000,3000,1000]);%3000*abs(pos-campos)/norm(pos-campos));
    camtarget(pos);

  end

     cam = campos;
    %Also keep the camera from going into the ground (could be done a
    %smarter way to also not look through the ground).
     nearestGroundZ = interp2(groundSurf.XData,groundSurf.YData,groundSurf.ZData,cam(1),cam(2));
     if cam(3)<nearestGroundZ
      	campos([cam(1),cam(2),nearestGroundZ]);
     end

    told = tnew;
    pause(0.01);

end
end

function KeyPress(varargin)
     fig = varargin{1};
     key = varargin{2}.Key;
     if strcmp(key,'e')
         fig.UserData.e = true;
     elseif strcmp(key,'q')
         fig.UserData.q = true;
     elseif strcmp(key,'a')
         fig.UserData.a = true;
     elseif strcmp(key,'d')
         fig.UserData.d = true;
     elseif strcmp(key,'w')
         fig.UserData.w = true;
     elseif strcmp(key,'s')
         fig.UserData.s = true;
     elseif strcmp(key,'v')
         fig.UserData.firstPerson = ~fig.UserData.firstPerson;
     end
end

function KeyRelease(varargin)
     fig = varargin{1};
     key = varargin{2}.Key;
     if strcmp(key,'e')
         fig.UserData.e = false;
     elseif strcmp(key,'q')
         fig.UserData.q = false;
     elseif strcmp(key,'a')
         fig.UserData.a = false;
     elseif strcmp(key,'d')
         fig.UserData.d = false;
     elseif strcmp(key,'w')
         fig.UserData.w = false;
     elseif strcmp(key,'s')
         fig.UserData.s = false;
     end
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
% DEMO_MGP ─ A practical example of the multivariate Gaussian process model
%            with different probabilistic models
%
% Description :
%
%   For many biological real-case scenarios two (or more)
%   animal populations can be strongly interwoven. The way in
%   which species dependent on each other are manifold. For instance,
%   predator-prey relationship, competition, symbiosis, etc.
%
%   It is fundamental to develop mechanistic understanding on the
%   underlying processes that governs the natural growth of the population.
%   However this goal is sometimes difficult and often the models are too
%   complex to derive from first principles. Nowadays, there are many
%   mechanistic models which are developed from first principles and that they
%   are able to describe the population dynamics in a very accurate manner
%   (system of ODEs or PDEs). Nevertheless, in many practical applications
%   it is easier to acquire data from the realization of such underlying biological
%   processes and infer their most probable form. Therefore Bayesian
%   nonparametric models become useful tools.
%
%   The example presented here uses a classical data-set from two animal
%   populations in Canada. The Lynx (Lynx canadensis) and the snowshoe hare
%   (Lepus americanus). The lynx population rise and fall according to the
%   variations in the populations of snowshoe hares over time. That is it, when
%   hares are abundant, lynx populations expand, and when the density of hares is
%   low, the population of lynx shrinks (there are many reasons why).
%   This is know as the predator-prey interaction.
%
%   In this example we want exemplify the joint modelling of the two species
%   population (abundance) aforementioned when data is missing in different
%   time intervals (for both species) and how one species can inform the
%   abundance of each another.
%
%   We assume that the number of lynx and hares follow the
%   negative-binomial distribution given unknown latent values (function
%   values of some unspecified function form for the regression), i.e.,
%
%       N1(t)|f1(t) ~ Neg-Binomial(exp(f1(t)), r1)  (hare)
%       N2(t)|f2(t) ~ Neg-Binomial(exp(f2(t)), r2)  (lynx)
%
%   and that N1(t)|f1(t) is independent of N2(t)|f2(t). The expression exp(f1(t))
%   is the expected value of N1 given f1 and exp(f2(t)) is the expected value
%   of N2 given f2. The correlation is now introduced through the multivariate
%   Gaussian process prior assuming the linear model of coregionalization, i.e.,
%
%       [f1(t), f2(t')] ~ MGP(0, K)
%
%   where K is a full covariance matrix formed by covariance function
%
%       k(fj(t), fj'(t')) = Sig_(c = 1, 2) u_c(j, j') k_c(t, t')
%
%   with u_c(j, j') the (j, j') entry of the matrix U_c = Lc Lc^T, where Lc
%   is cth column of the Cholesky decomposition of the coregionalization
%   matrix (covariance matrix) {Sig}j,j' = sig_j sig_j' rho_{j, j'}. We use the
%   analytical approximation expectation-propagation to carry out the inference
%   over the latent values f1, f2 and then calcule the unconditional expected values
%   of N1 and N2, i.e, E[N1(t)|Y1 = y1, Y2 = y2] and E[N2(t)|Y1 = y1, Y2 = y2].
%
% Additional references :
%
%   Murray, J, D (2002). Mathematical biology. Third Edition. Springer
%    Series.
%
%   Gelfand et al. (2004). Nonstationary multivariate process modelling
%     through spatially varying coregionalization.
%
%   Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian
%    Processes for Machine Learning. The MIT Press.
%
%   Bernardo, J and Smith, A (2008). Bayesian Theory. Wiley series in
%    probability and statistics
%
% -------------- Marcelo Hartmann

%%
% ------- data analysis with the multivariate Gaussian process model

% download the data
S = which('demo_multivariategp');
L = strrep(S, 'demo_multivariategp.m', 'demodata/predpreydata.txt');
data = importdata(L);

% colors;
col = [0 0 1; 1 0 0];

% all time points (yearly)
t = data.data(:, 1);

% full data
yy = [round(data.data(:, 2)); round(data.data(:, 3))];

% take sparse data (you can change the years);
yr1 = [1870, 1900];
% yr1 = [1846,  1935]; one observation
indt1 = ~logical((t >= yr1(1)) .* (t <= yr1(2)));
t1 = t(indt1);

yr2 = [1850, 1870];
% yr2 = [1845, 1934]; one observation
indt2 = ~logical((t >= yr2(1)) .* (t <= yr2(2)));
t2 = t(indt2);

% hare (thousands)
y1 = round(data.data(indt1, 2));

% lynx (thousands)
y2 = round(data.data(indt2, 3));

% predator-prey data;
y = [y1; y2];

% create species markers
c = [ones(size(t1)); 2*ones(size(t2))];

% data gp-format - species markers in the last column
x = [[t1; t2] c];
z = [ones(size(x, 1), 1) c];

% likelihood structure for each species (the likelihood can also be distinct)
lik1 = lik_negbin;
lik2 = lik_negbin;
likS = {lik1 lik2};
lik = lik_liks('likelihoods', likS, 'classVariables', 2);

% correlation function for each species (the correlation functions can also be distinct)
k1 = gpcf_sexp('magnSigma2', 1, 'magnSigma2_prior', prior_fixed, 'selectedVariables', 1);
k2 = gpcf_sexp('magnSigma2', 1, 'magnSigma2_prior', prior_fixed, 'selectedVariables', 1);

% the linear model of coregionalization (multivariate Gaussian process model)
k = gpcf_covar('numberClass', 2, 'classVariables', 2, ...
    'R_prior', prior_corrunif('nu', 3), 'corrFun', {k1 k2});

% set gp structure
% for the EP approximation
gp = gp_set('lik', lik, 'cf', k, 'latent_method', 'EP');

% For the Laplace approximation
%gp = gp_set('lik', lik, 'cf', k, 'latent_method', 'Laplace');

% optimzation options
opt = optimset('TolFun', 1e-3, 'TolX', 1e-5, 'Display', 'iter');

% random initialization (this is important ...)
% wini = 2.* randn(size(gp_pak(gp))); gp = gp_unpak(gp, wini);

% map ─ type-II maximum likelihood
gp = gp_optim(gp, x, y, 'z', z, 'opt', opt);

% check the gradient ...
% gp_g(gp_pak(gp), gp, x, y, 'z', z)';

% take the parameters, hyperparameters estimates
[th ss] = gp_pak(gp);

% correlation matrix estimate
%rho = k.fh.RealToRho(th(end - 2 - 2), 1, [])
Corr = gp.cf{1}.fh.sigma(gp.cf{1}, 'corr');
Corr{1}

% prediction
np = 200;
xp = repmat(linspace(min(x(:, 1)), max(x(:, 1)), np)', 2, 1);
xp = [xp repelem((1:2), np)'];
zp = ones(size(xp, 1), 1);

% leave-one-out
[~, ~, lpyt] = gp_loopred(gp, x, y,'z', z);

% prediction for new observations
[Ef, Varf, ~, Ey, Vary] = gp_pred(gp, x, y, xp, 'z', z, 'zt', [zp xp(:, 2)]);

% --- visualize predictions

figure; hold on;

subplt = subplot(2, 1, 1); hold on
subplt.XLim = [min(xp(:, 1)), max(xp(:, 1))];
psubplt = get(subplt, 'pos');
psubplt(4) = psubplt(4) + 0.03;
psubplt(3) = psubplt(3) + 0.04;
psubplt(1) = psubplt(1) - 0.04;
set(subplt, 'pos', psubplt);

for j = 1:2
    indj = xp(:, 2) == j;
    indX = xp(indj, 1);

    pl(j) = plot(xp(indj, 1), Ey(indj), 'color', col(j, :), 'LineWidth', 2);
 end

% visualize observed (training) predator-pray data
pl(3) = plot(t1(t1 <= yr1(1)), y1(t1 <= yr1(1)), '.', ...
    'MarkerSize', 23, 'lineWidth', 2, 'color', col(1, :));
plot(t1(t1 >= yr1(2)), y1(t1 >= yr1(1)), '.', ...
    'MarkerSize', 23, 'lineWidth', 2, 'color', col(1, :));
pl(4) = plot(t2(t2 <= yr2(1)), y2(t2 <= yr2(1)), '.', ...
    'MarkerSize', 23, 'lineWidth', 2, 'color', col(2, :));
plot(t2(t2 >= yr2(2)), y2(t2 >= yr2(2)), '.', ...
    'MarkerSize', 23, 'lineWidth', 2, 'color', col(2, :));

xlabel('years'); ylabel('Population size');
title('multivariate Gaussian process model');

legend(pl, 'E[N_1(t_*)|N_1 = n_1, N_2 = n_2]', 'E[N_2(t_*)|N_1 = n_1, N_2 = n_2]', ...
    'hare-data', 'lynx-data', 'Location', 'northeast');

%%
% -------- data analysis with independent Gaussian process model

EfI = {}; VarfI = {}; EyI = {}; VaryI = {};
gpI = {};

subplt = subplot(2, 1, 2); hold on
subplt.XLim = [min(xp(:, 1)), max(xp(:, 1))];
psubplt = get(subplt, 'pos');
psubplt(4) = psubplt(4) + 0.05;
psubplt(3) = psubplt(3) + 0.04;
psubplt(1) = psubplt(1) - 0.03;
set(subplt, 'pos', psubplt);

for j = 1:2
    inddj = x(:, 2) == j;
    cf = gpcf_sexp('selectedVariables', 1);

    gpI{j} = gp_set('lik', likS{j}, 'cf', cf);
    gpI{j} = gp_optim(gpI{j}, x(inddj, :), y(inddj), 'opt', opt);

    indj = xp(:, 2) == j;
    indX = xp(indj, 1);

    [EfI{j}, VarfI{j}, lpyt, EyI{j}, VaryI{j}] = ...
        gp_pred(gpI{j}, x(inddj, :), y(inddj), xp(indj, :));

    pl(j) = plot(xp(indj, 1), EyI{j}, 'color', col(j, :), 'LineWidth', 1.7);
end

% visualize observed (training) predator-pray data
pl(3) = plot(t1(t1 <= yr1(1)), y1(t1 <= yr1(1)), '.', ...
    'MarkerSize', 23, 'lineWidth', 2, 'color', col(1, :));
plot(t1(t1 >= yr1(2)), y1(t1 >= yr1(1)), '.', ...
    'MarkerSize', 23, 'lineWidth', 2, 'color', col(1, :));
pl(4) = plot(t2(t2 <= yr2(1)), y2(t2 <= yr2(1)), '.', ...
    'MarkerSize', 23, 'lineWidth', 2, 'color', col(2, :));
plot(t2(t2 >= yr2(2)), y2(t2 >= yr2(2)), '.', ...
    'MarkerSize', 23, 'lineWidth', 2, 'color', col(2, :));

xlabel('years'); ylabel('Population size');
title('independent Gaussian process models');

legend(pl, 'E[N1(t)|Y1 = y1, Y2 = y2]', 'E[N2(t)|Y1 = y1, Y2 = y2]', ...
    'hare-data', 'lynx-data', 'Location', 'northeast');

%%
% ------ all measurements (compare the difference ...)

figure; hold on;

for j = 1:2
    indj = x(:, 2) == j;
    plot(x(indj, 1), y(indj), '.', 'MarkerSize', 13, 'color', col(j, :))
    pd(j) =  plot(t, yy((j - 1) * 91 + [1:91]), '.-', ...
        'color', col(j, :), 'MarkerSize', 23, 'LineWidth', 2);
end

xlabel('years'); ylabel('Population size');
title('full data');

legend(pd, 'hare-data', 'lynx-data', 'Location', 'northeast');

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function linplus_test11 ( )

%*****************************************************************************80
%
%% TEST11 tests R83P_ML.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST11\n' );
  fprintf ( 1, '  R83P_ML computes A*x or A''*X\n' );
  fprintf ( 1, '    where A has been factored by R83P_FA.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );

  for job = 0 : 1
%
%  Set the matrix.
%
    [ a, seed ] = r83p_random ( n, seed );
%
%  Set the desired solution.
%
    x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
    if ( job == 0 )
      b = r83p_mxv ( n, a, x );
    else
      b = r83p_vxm ( n, a, x );
    end
%
%  Factor the matrix.
%
    [ a_lu, work2, work3, work4, info ] = r83p_fa ( n, a );

    if ( info ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TEST11 - Fatal error!\n' );
      fprintf ( 1, '  R83P_FA declares the matrix is singular!\n' );
      fprintf ( 1, '  The value of INFO is %d\n', info );
      return;
    end
%
%  Now multiply factored matrix times solution to get right hand side again.
%
    b2 = r83p_ml ( n, a_lu, x, job );

    if ( job == 0 )
      r8vec2_print_some ( n, b, b2, 10, '  A*x and PLU*x' );
    else
      r8vec2_print_some ( n, b, b2, 10, '  A''*x and (PLU)''*x' );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Fibres
%
%%
%
% The set of all rotations that rotate a certain vector u onto a certain
% vector v define a fibre in the rotation space. A discretisation of such a
% fibre is defined by

u = xvector;
v = yvector;
rot = rotation(fibre(u,v))

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [gradKeypointS,gradKeypointV,gradKeypointAlpha,gradOcclusionS,...
    gradOcclusionV,gradOcclusionAlpha,occlusionCount,gradSilS,gradSilV,...
    gradSilAlpha,gradShapePriorS,gradShapePriorV,gradShapePriorAlpha] ...
    = computeGrads(S,V,alpha,stateFiles,lambda,shapePriorInstance,corenum,numcores,params)

numpoints = size(S,1);
K = size(V,2);
numimages = length(stateFiles);

gradOcclusionS = zeros(size(S));
gradOcclusionV = zeros(size(V));
gradOcclusionAlpha = zeros(size(alpha));
occlusionCount = zeros(numpoints,1);

gradSilS = zeros(size(S));
gradSilV = zeros(size(V));
gradSilAlpha = zeros(size(alpha));

gradKeypointS = zeros(size(S));
gradKeypointV = zeros(size(V));
gradKeypointAlpha = zeros(size(alpha));

gradShapePriorS = zeros(size(S));
gradShapePriorV = zeros(size(V));
gradShapePriorAlpha = zeros(size(alpha));

%t_occ=0;t_sil=0;t_prior=0;
%t_kp = 0;

for i=1:numimages
    if(mod(i,numcores)==(corenum-1))

        load(stateFiles{i});
        if(params.opt.normalizeByViews)
            gradWeight = getGradWtFromClusters(state.viewClusters,state.globalID);
            gradWeight = gradWeight^(0.5);
        else
            gradWeight = 1;
        end
        S_i = S + reshape(V*alpha(:,i),numpoints,3);

        %% Keypoint Consistency
        %tic;
        kpGrad = keypointGrad(S_i,state,params);
        kpGrad = lambda(6)*kpGrad;
        gradKeypointS = gradKeypointS + kpGrad;
        gradKeypointV = gradKeypointV + kpGrad(:)*(alpha(:,i))';
        gradKeypointAlpha(:,i) = V'*kpGrad(:);
        %t_kp = t_kp + toc;

        %% Occlusion Reasoning
        %tic;
        [occGrad,occCount] = occlusionGradOptimal(S_i,state); %we also keep track of how many constraints each point violates
        occGrad = lambda(1)*occGrad;
        occGrad = occGrad*gradWeight;
        occlusionCount = occCount + occlusionCount;
        gradOcclusionS = gradOcclusionS + occGrad;
        gradOcclusionV = gradOcclusionV + occGrad(:)*(alpha(:,i))';
        gradOcclusionAlpha(:,i) = V'*occGrad(:);
        %t_occ = t_occ+toc;

        %% Silhouette Consistency
        %tic;
        silGrad = lambda(2)*silhouetteGrad(S_i,state,params); %we also keep track of how many constraints each point violates
        silGrad = silGrad*500; % because we divide by number of points on silhouette
        silGrad = silGrad*gradWeight;
        gradSilS = gradSilS + silGrad;
        gradSilV = gradSilV + silGrad(:)*(alpha(:,i))';
        gradSilAlpha(:,i) = V'*silGrad(:);
        %t_sil = t_sil + toc;

        %% Shape Prior
        %tic;
        shapePriorGrad = lambda(5)*shapePriorInstance(S_i);
        gradShapePriorS = gradShapePriorS + shapePriorGrad;
        gradShapePriorV = gradShapePriorV + shapePriorGrad(:)*(alpha(:,i))'; %%Verify
        gradShapePriorAlpha(:,i) = V'*shapePriorGrad(:); %% Verify
        %t_prior = t_prior + toc;
    end
end
end

repo: CategoryShapes


text:
function [vertex,face] = read_vtk(filename, verbose)

% read_vtk - read data from VTK file.
%
%   [vertex,face] = read_vtk(filename, verbose);
%
%   'vertex' is a 'nb.vert x 3' array specifying the position of the vertices.
%   'face' is a 'nb.face x 3' array specifying the connectivity of the mesh.
%
%   Copyright (c) Mario Richtsfeld

if nargin<2
    verbose = 1;
end

fid = fopen(filename,'r');
if( fid==-1 )
    error('Can''t open the file.');
    return;
end

str = fgets(fid);   % -1 if eof
if ~strcmp(str(3:5), 'vtk')
    error('The file is not a valid VTK one.');
end

%%% read header %%%
str = fgets(fid);
str = fgets(fid);
str = fgets(fid);
str = fgets(fid);
nvert = sscanf(str,'%*s %d %*s', 1);

% read vertices
[A,cnt] = fscanf(fid,'%f %f %f', 3*nvert);
if cnt~=3*nvert
    warning('Problem in reading vertices.');
end
A = reshape(A, 3, cnt/3);
vertex = A;

% read polygons
str = fgets(fid);
str = fgets(fid);

info = sscanf(str,'%c %*s %*s', 1);

if((info ~= 'P') && (info ~= 'V'))
    str = fgets(fid);
    info = sscanf(str,'%c %*s %*s', 1);
end

if(info == 'P')

    nface = sscanf(str,'%*s %d %*s', 1);

    [A,cnt] = fscanf(fid,'%d %d %d %d\n', 4*nface);
    if cnt~=4*nface
        warning('Problem in reading faces.');
    end

    A = reshape(A, 4, cnt/4);
    face = A(2:4,:)+1;
end

if(info ~= 'P')
    face = 0;
end

% read vertex ids
if(info == 'V')

    nv = sscanf(str,'%*s %d %*s', 1);

    [A,cnt] = fscanf(fid,'%d %d \n', 2*nv);
    if cnt~=2*nv
        warning('Problem in reading faces.');
    end

    A = reshape(A, 2, cnt/2);
    face = repmat(A(2,:)+1, 3, 1);
end

if((info ~= 'P') && (info ~= 'V'))
    face = 0;
end

fclose(fid);

return

author: gpeyre

repo: matlab-toolboxes


text:
%% vorMap2triSurf
% Below is a demonstration of the features of the |vorMap2triSurf| function

%%
clear; close all; clc;

%% Syntax
% |[Fd,Vd]=vorMap2triSurf(F,V,indSeed,seedIndex);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function w =  get_w(a,d)

w= (a.alpha' * a.child.dat.X)';

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [nb] = find_triangle_neighbours(pnt, tri)

% FIND_TRIANGLE_NEIGHBOURS determines the three neighbours for each triangle
% in a mesh. It returns NaN's if the triangle does not have a neighbour on
% that particular side.
%
% [nb] = find_triangle_neighbours(pnt, tri)

% Copyright (C) 2003, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

npnt = size(pnt,1);
ntri = size(tri,1);

% each triangle has maximally three neighbours, assuming that the
% surface mesh is not degenerate
nb = nan(size(tri));

% for i=1:ntri
%   for j=setdiff(1:ntri, i)
%     if length(intersect(tri(i,[1 2]), tri(j,:)))==2
%       nb(i,1) = j;
%       continue;
%     end
%     if length(intersect(tri(i,[2 3]), tri(j,:)))==2
%       nb(i,2) = j;
%       continue;
%     end
%     if length(intersect(tri(i,[3 1]), tri(j,:)))==2
%       nb(i,3) = j;
%       continue;
%     end
%   end
%   if all(~isnan(nb(i,:)))
%     continue;
%   end
% end

for i=1:ntri
  % find all neighbouring triangles
  tmp1 = (tri==tri(i,1));
  tmp2 = (tri==tri(i,2));
  tmp3 = (tri==tri(i,3));
  tmp  = (tmp1|tmp2|tmp3);
  sel = find(sum(tmp,2)==2);

  % ensure that each neighbour is assigned to the proper edge
  if length(sel)>3
    ft_error('more than three neighbours found for triangle %d', i);
  else
    for j=1:length(sel)
      if isempty(setdiff(intersect(tri(i,:), tri(sel(j),:)), tri(i,[1 2])))
        nb(i,1) = sel(j);
      elseif isempty(setdiff(intersect(tri(i,:), tri(sel(j),:)), tri(i,[2 3])))
        nb(i,2) = sel(j);
      elseif isempty(setdiff(intersect(tri(i,:), tri(sel(j),:)), tri(i,[3 1])))
        nb(i,3) = sel(j);
      end
    end
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%
% Copyright (c) 2019-present, Mahmoud Afifi
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
%
% Please, cite the following paper if you use this code:
%
% Mahmoud Afifi and Michael S. Brown. What else can fool deep learning?
% Addressing color constancy errors on deep neural network performance.
% ICCV, 2019
%
% Email: mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%%

function varargout = demo_GUI(varargin)
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @demo_GUI_OpeningFcn, ...
    'gui_OutputFcn',  @demo_GUI_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end

function demo_GUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
guidata(hObject, handles);
warning off
addpath('classes');

global WB_emulator
load('synthWBmodel.mat');
handles.save.Enable = 'off';

function varargout = demo_GUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

% --- Executes on button press in browse.
function browse_Callback(hObject, eventdata, handles)
% hObject    handle to browse (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global Path_Name
global File_Name
%global I_temp
global I
global outImages
global WB_emulator
global feature
Old_fileName = File_Name;
Old_pathName = Path_Name;
[File_Name, Path_Name] = uigetfile({'*.jpg';'*.png';'*.jpeg'},'Select input image',fullfile('..','..','images'));

if File_Name == 0
    File_Name = Old_fileName;
    Path_Name = Old_pathName;
    if sum(I(:)) ~= 0
        handles.save.Enable = 'on';
    else
        handles.save.Enable = 'off';
    end
else
    I = imread(fullfile(Path_Name,File_Name));
    axes(handles.image);
    handles.image.Visible = 'On';
    handles.input_text.Visible = 'On';
    imshow(I);
    pause(0.001);
    handles.status.String = 'Processing...';pause(0.001);
    k = handles.k.Value;
    sigma = handles.sigma.Value;
    WB_emulator.K = round(k);

%     if handles.wb_checkbox.Value==1
%         WBmodel = load('WBmodel.mat');
%         I_temp = WBmodel.WBmodel.correctImage(I);
%     else
%         I_temp = I ;
%     end
    hist = WB_emulator.RGB_UVhist(im2double(I));
    feature = WB_emulator.encode(hist);
    outImages = WB_emulator.generate_wb_srgb(I,10,feature,sigma);
    for i = 1 : size(outImages,4)
        tempstr= WB_emulator.wb_photo_finishing{i};
        eval(sprintf('axes(handles.%s);',tempstr(2:end)));
        eval(sprintf('handles.%s.Visible=''On'';',tempstr(2:end)));
        eval(sprintf('handles.%s_text.Visible=''On'';',tempstr(2:end)));
        imshow(outImages(:,:,:,i));
    end
    handles.save.Enable = 'on';
    handles.status.String = 'Done!';
    pause(0.01); handles.status.String = '';
end

% --- Executes on button press in save.
function save_Callback(hObject, eventdata, handles)
% hObject    handle to save (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

global Path_Name
global File_Name
global outImages
global WB_emulator
global I
[~,name,ext] = fileparts(File_Name);
outFile_Name = [name 'SynthWB' ext];
[file,path,~] = uiputfile({'*.jpg';'*.png'},'Save Image',fullfile(Path_Name,outFile_Name));
if file ~=0
    [~,name,ext] = fileparts(file);
    handles.status.String = 'Processing...';pause(0.001);
    imwrite(I,fullfile(path,sprintf('%s%s%s',name,'_original',ext)));
    for j =1 : size(outImages,4)
        imwrite(outImages(:,:,:,j),fullfile(path,...
            sprintf('%s%s%s',name,...
            WB_emulator.wb_photo_finishing{j},ext)));
    end
    handles.status.String = 'Done!';
    pause(0.01); handles.status.String = '';
end

% --- Executes on slider movement.
function k_Callback(hObject, eventdata, handles)
% hObject    handle to k (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
global WB_emulator
global I
global feature
global outImages
axes(handles.image);
handles.status.String = 'Processing...';pause(0.001);
k = handles.k.Value;
sigma = handles.sigma.Value;
WB_emulator.K = round(k);
outImages = WB_emulator.generate_wb_srgb(I,10,feature,sigma);
for i = 1 : size(outImages,4)
    tempstr= WB_emulator.wb_photo_finishing{i};
    eval(sprintf('axes(handles.%s);',tempstr(2:end)));
    eval(sprintf('handles.%s.Visible=''On'';',tempstr(2:end)));
    eval(sprintf('handles.%s_text.Visible=''On'';',tempstr(2:end)));
    imshow(outImages(:,:,:,i));
end
handles.status.String = 'Done!';
pause(0.01); handles.status.String = '';

% --- Executes during object creation, after setting all properties.
function k_CreateFcn(hObject, eventdata, handles)
% hObject    handle to k (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function sigma_Callback(hObject, eventdata, handles)
% hObject    handle to sigma (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

global WB_emulator
global I
global feature
global outImages
axes(handles.image);
handles.status.String = 'Processing...';pause(0.001);
k = handles.k.Value;
sigma = handles.sigma.Value;
WB_emulator.K = round(k);
outImages = WB_emulator.generate_wb_srgb(I,10,feature,sigma);
for i = 1 : size(outImages,4)
    tempstr= WB_emulator.wb_photo_finishing{i};
    eval(sprintf('axes(handles.%s);',tempstr(2:end)));
    eval(sprintf('handles.%s.Visible=''On'';',tempstr(2:end)));
    eval(sprintf('handles.%s_text.Visible=''On'';',tempstr(2:end)));
    imshow(outImages(:,:,:,i));
end
handles.status.String = 'Done!';
pause(0.01); handles.status.String = '';

% --- Executes during object creation, after setting all properties.
function sigma_CreateFcn(hObject, eventdata, handles)
% hObject    handle to sigma (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

%
% % --- Executes on button press in wb_checkbox.
% function wb_checkbox_Callback(hObject, eventdata, handles)
% % hObject    handle to wb_checkbox (see GCBO)
% % eventdata  reserved - to be defined in a future version of MATLAB
% % handles    structure with handles and user data (see GUIDATA)
%
% % Hint: get(hObject,'Value') returns toggle state of wb_checkbox
%
% %global I_temp
% global outImages
% global WB_emulator
% global feature
% global I
% pause(0.001);
% handles.status.String = 'Processing...';pause(0.001);
% k = handles.k.Value;
% sigma = handles.sigma.Value;
% WB_emulator.K = round(k);
% %
% % if handles.wb_checkbox.Value==1
% %     WBmodel = load('WBmodel.mat');
% %         I_temp = WBmodel.WBmodel.correctImage(I);
% % else
% %     I_temp = I ;
% % end
% hist = WB_emulator.RGB_UVhist(im2double(I));
% feature = WB_emulator.encode(hist);
% outImages = WB_emulator.generate_wb_srgb(I,10,feature,sigma);
% for i = 1 : size(outImages,4)
%     tempstr= WB_emulator.wb_photo_finishing{i};
%     eval(sprintf('axes(handles.%s);',tempstr(2:end)));
%     eval(sprintf('handles.%s.Visible=''On'';',tempstr(2:end)));
%     eval(sprintf('handles.%s_text.Visible=''On'';',tempstr(2:end)));
%     imshow(outImages(:,:,:,i));
% end
% handles.status.String = 'Done!';
% pause(0.01); handles.status.String = '';

author: mahmoudnafifi


text:
%% SAUSAGEPLOT
function sausageplot(obj, numBinsToUse, varargin)
    %sausageplot
    %    Under development, this function attempts to replicate the
    %    capabilities of sausageplot.xpy written by Glenn Thompson
    %    at AVO.
    p = inputParser;
    p.addRequired('numBinsToUse', @isnumeric);
    p.addOptional('numPoints', 1, @isnumeric);
    p.addOptional('radii', [], @isnumeric);
    p.parse(numBinsToUse, varargin{:});
    numBinsToUse = p.Results.numBinsToUse;
    numPoints = p.Results.numPoints;
    radii = p.Results.radii;
%             % Page size, dots-per-inch and scale settings
%             fh = figure()
%             if numBinsToUse > NUMPOINTS
%                 dpi = 6*(numBinsToUse+1);
%                 axes_width = 0.8 * (NUMPOINTS + 1) / (numBinsToUse + 1);
%                 axes_height = 0.8;
%             else
%                 dpi = 6*(NUMPOINTS+1);
%                 axes_width = 0.8;
%                 axes_height = 0.8 * (numBinsToUse + 1) / (NUMPOINTS + 1);
%             end
%             fh.set_dpi(dpi)
%             fh.set_size_inches((10.0,10.0),forward=True)
%             fig2ax1 = fig2.add_axes([0.1, 0.85-axes_height, axes_width, axes_height]);
%             print_pixels(fig2, fig2ax1, numBinsToUse, NUMPOINTS);
%             colormapname = 'hot_r';
    MAXMAGCOLORBAR = 5.0;
    MINMAGCOLORBAR = 1.0;
    % Marker size configuration
%             SCALEFACTOR = 84.0 / dpi;
%             MAXMARKERSIZE = 43.0 * SCALEFACTOR;
%             MINMARKERSIZE = 4.0 * SCALEFACTOR;
    PTHRESHOLD = 50;
    for c = 1:length(obj)
        % set up weekending list for yticklabels
        binendstr = {};
        for bin_index=numel(obj.time)-numBinsToUse+1: numel(obj.time)

            dstr = datestr(obj.time(bin_index))
            binendstr{bin_index} = dstr;
        end
        pointlabels = {};
        for i=1:length(numPoints)
            % filter here based on points and radii
            j = 1:length(obj(c).counts); % replace this
            counts = obj(c).counts(j);
            cummag = obj(c).cum_mag(j);
            prcntile = percentiles(obj(c).counts);
            %pointlabels{i} = sprintf('%s(%d)', point_label{i}, counts(-1));
            pointlabels{i}='';
            for bin_index=numel(obj.time)-numBinsToUse+1: numel(obj.time)%-numBinsToUse-1: 1: 0
                y = obj(c).counts(bin_index);
                magnitude = obj(c).cum_mag(bin_index);
                p = y2percentile(y,prcntile);
                if y>0
                    colorVal = scalarMap.to_rgba(magnitude)
                    msize = MINMARKERSIZE + (p-PTHRESHOLD) * (MAXMARKERSIZE - MINMARKERSIZE) / (100-PTHRESHOLD);
                    if msize<MINMARKERSIZE
                        msize=MINMARKERSIZE;
                    end
%                             fig2ax1.plot(i+0.5, w, 's', color=colorVal, markersize=msize, linewidth=0 );
                    scatter(i+0.5, bin_index, msize, y, 's', 'filled')
                    if msize > MAXMARKERSIZE * 0.3
                        text(i+0.5, bin_index, num2str(y))
%                                fig2ax1.text(i+0.5, w, '%d', y, horizontalalignment='center', verticalalignment='center', fontsize = 8 * SCALEFACTOR)
                    end
                end
            end
        end

        % Adding xticks, yticks, labels, grid
%                 fig2ax1.set_axisbelow(True) % I think this puts grid and tickmarks below actual data plotted

        % x-axis
%                 fig2ax1.set_xticks(np.arange(.5,NUMVOLCANOES+.5,1))
        set(gca, 'XTick', 0.5:1:numPoints+0.5);
%                 fig2ax1.set_xlim([-0.5, NUMVOLCANOES+0.5])
        set(gca, 'XLim', [-0.5 numPoints+0.5])
%                 fig2ax1.xaxis.grid(True, linestyle='-', color='gray')
        grid on;
        set(gca, 'XTickLabel', pointlabels)
%                 fig2ax1.xaxis.set_ticks_position('top')
%                 fig2ax1.xaxis.set_label_position('top')
        %plt.setp( fig2ax1.get_xticklabels(), rotation=45, horizontalalignment='left', fontsize=10*SCALEFACTOR )

        % y-axis
%                 fig2ax1.set_yticks(np.arange(-number_of_weeks_to_plot-0.5, 0, 1))
        set(gca, 'YTick', -numBinsToUse-0.5: 1: 0)
%                 fig2ax1.set_yticklabels(weekending)
        set(gca, 'YTickLabel', binendstr)
%                 fig2ax1.set_ylim([-number_of_weeks_to_plot - 0.5, -0.5])
        set(gca, 'YLim', -numBinsToUse -0.5 : 1 : -0.5)
%                 plt.setp( fig2ax1.get_yticklabels(), fontsize=10*SCALEFACTOR )
    end
end

function pvalue(obj)
    logt = log10(obj.time-obj.time(1));
    logc = log10(obj.counts);
    plot(logt, logc);
    p = logt \ logc'
%     p = polyfit(logt, logc, 1)
%     hold on
%     logt1 = min(logt);
%     logc1 = polyval(p, logt1);
%     logt2 = max(logt);
%     logc2 = polyval(p, logt2);
%     line([logt1 logt2], [logc1 logc2]);
end

%% PERCENTILES
function p=percentiles(vals)
    lenVals = length(vals);
    for i=1:100
        p(i) = vals(floor(i/100 * (lenVals-1))+1);
    end
end

%% PLOT_PERCENTILES
function plot_percentiles(p)
    figure(get(gcf,'Number')+1, 'Color', [1 1 1])
    plot(1:100, p);
end

author: geoscience-community-codes

repo: GISMO