text:
function [h, T2] = ft_plot_slice(dat, varargin)

% FT_PLOT_SLICE plots a 2-D cut through a 3-D volume and interpolates if needed
%
% Use as
%   ft_plot_slice(dat, ...)
% or
%   ft_plot_slice(dat, mask, ...)
% where dat and mask are equal-sized 3-D arrays.
%
% Additional options should be specified in key-value pairs and can be
%   'transform'    = 4x4 homogeneous transformation matrix specifying the mapping from
%                    voxel coordinates to the coordinate system in which the data are plotted.
%   'location'     = 1x3 vector specifying a point on the plane which will be plotted
%                    the coordinates are expressed in the coordinate system in which the
%                    data will be plotted. location defines the origin of the plane
%   'orientation'  = 1x3 vector specifying the direction orthogonal through the plane
%                    which will be plotted (default = [0 0 1])
%   'unit'         = string, can be 'm', 'cm' or 'mm (default is automatic)
%   'resolution'   = number (default = 1 mm)
%   'datmask'      = 3D-matrix with the same size as the data matrix, serving as opacitymap
%                    If the second input argument to the function contains a matrix, this
%                    will be used as the mask
%   'maskstyle'    = string, 'opacity' or 'colormix', defines the rendering
%   'background'   = needed when maskstyle is 'colormix', 3D-matrix with
%                    the same size as the data matrix, serving as
%                    grayscale image that provides the background
%   'opacitylim'   = 1x2 vector specifying the limits for opacity masking
%   'interpmethod' = string specifying the method for the interpolation, see INTERPN (default = 'nearest')
%   'style'        = string, 'flat' or '3D'
%   'colormap'     = string, see COLORMAP
%   'clim'         = 1x2 vector specifying the min and max for the colorscale
%
% You can plot the slices from the volume together with an intersection of the slices
% with a triangulated surface mesh (e.g. a cortical sheet) using
%   'intersectmesh'       = triangulated mesh, see FT_PREPARE_MESH
%   'intersectcolor'      = string, color specification
%   'intersectlinestyle'  = string, line specification
%   'intersectlinewidth'  = number
%
% See also FT_PLOT_ORTHO, FT_PLOT_MONTAGE, FT_SOURCEPLOT

% Undocumented options
%   'plotmarker'     = Nx3 matrix with points to be plotted as markers, e.g. dipole positions
%   'markersize'
%   'markercolor'

% Copyrights (C) 2010-2014, Jan-Mathijs Schoffelen
% Copyrights (C) 2014-2016, Robert Oostenveld and Jan-Mathijs Schoffelen
%
% 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$

persistent dim X Y Z

if isequal(dim, size(dat(:,:,:,1,1)))
  % reuse the persistent variables to speed up subsequent calls with the same input
else
  dim       = size(dat);
  [X, Y, Z] = ndgrid(1:dim(1), 1:dim(2), 1:dim(3));
end

% parse first input argument(s). it is either
% (dat, varargin)
% (dat, msk, varargin)
% (dat, [], varargin)
if numel(varargin)>0 && (isempty(varargin{1}) || isnumeric(varargin{1}) || islogical(varargin{1}))
  msk      = varargin{1};
  varargin = varargin(2:end);
end

% get the optional input arguments
transform           = ft_getopt(varargin, 'transform', eye(4));
loc                 = ft_getopt(varargin, 'location');
ori                 = ft_getopt(varargin, 'orientation', [0 0 1]);
unit                = ft_getopt(varargin, 'unit');       % the default will be determined further down
resolution          = ft_getopt(varargin, 'resolution'); % the default depends on the units and will be determined further down
datmask             = ft_getopt(varargin, 'datmask');
maskstyle           = ft_getopt(varargin, 'maskstyle', 'opacity');
background          = ft_getopt(varargin, 'background');
opacitylim          = ft_getopt(varargin, 'opacitylim');
interpmethod        = ft_getopt(varargin, 'interpmethod', 'nearest');
cmap                = ft_getopt(varargin, 'colormap');
clim                = ft_getopt(varargin, 'clim');
doscale             = ft_getopt(varargin, 'doscale', true); % only scale when necessary (time consuming), i.e. when plotting as grayscale image & when the values are not between 0 and 1
h                   = ft_getopt(varargin, 'surfhandle', []);

mesh                = ft_getopt(varargin, 'intersectmesh');
intersectcolor      = ft_getopt(varargin, 'intersectcolor', 'yrgbmyrgbm');
intersectlinewidth  = ft_getopt(varargin, 'intersectlinewidth', 2);
intersectlinestyle  = ft_getopt(varargin, 'intersectlinestyle');

plotmarker          = ft_getopt(varargin, 'plotmarker');
markersize          = ft_getopt(varargin, 'markersize', 'auto');
markercolor         = ft_getopt(varargin, 'markercolor', 'w');

% convert from yes/no/true/false/0/1 into a proper boolean
doscale = istrue(doscale);

if ~isa(dat, 'double')
  dat = cast(dat, 'double');
end

if exist('msk', 'var') && isempty(datmask)
  ft_warning('using the second input argument as mask rather than the one from the varargin list');
  datmask = msk; clear msk;
end

% normalise the orientation vector to one
ori = ori./sqrt(sum(ori.^2));

% set the default location
if isempty(loc) && (isempty(transform) || isequal(transform, eye(4)))
  loc = (dim+1)./2;
elseif isempty(loc)
  loc = [0 0 0];
end

% shift the location to be along the orientation vector
loc = ori*dot(loc,ori);

% it should be a cell-array
if isstruct(mesh)
  tmp = mesh;
  mesh = cell(size(tmp));
  for i=1:numel(tmp)
    mesh{i} = tmp(i);
  end
elseif iscell(mesh)
  % do nothing
else
  mesh = {};
end

% replace pnt by pos
for k = 1:numel(mesh)
  mesh{k} = fixpos(mesh{k});
end

dointersect = ~isempty(mesh);
if dointersect
  for k = 1:numel(mesh)
    if ~isfield(mesh{k}, 'pos') || ~isfield(mesh{k}, 'tri')
      % ft_error('the mesh should be a structure with pos and tri');
      mesh{k}.pos = [];
      mesh{k}.tri = [];
    end
  end
end

% check whether the mask is ok
domask = ~isempty(datmask);
if domask
  if ~isequal(size(dat), size(datmask)) && ~isequal(cmap, 'rgb')
    % the exception is when the functional data is to be interpreted as rgb
    ft_error('the mask data should have the same dimensions as the functional data');
  end
end

dobackground = ~isempty(background);
if dobackground
  if ~isequal(size(dat), size(background))
    error('the background data should have the same dimensions as the functional data');
  end
end

% determine the voxel center
% voxel_center_vc = [X(:) Y(:) Z(:)];
% voxel_center_hc = ft_warp_apply(transform, voxel_center_vc);

% determine the edges, i.e. the corner points of each voxel
% [Xe, Ye, Ze] = ndgrid(0:dim(1), 0:dim(2), 0:dim(3));
% Xe = Xe+0.5;
% Ye = Ye+0.5;
% Ze = Ze+0.5;
% voxel_edge_vc = [Xe(:) Ye(:) Ze(:)];
% voxel_edge_hc = ft_warp_apply(transform, voxel_edge_vc);

% determine the corner points of the box encompassing the whole data block
% extend the box with half a voxel  in all directions to get the outer edge
corner_vc = [
  0.5        0.5        0.5
  0.5+dim(1) 0.5        0.5
  0.5+dim(1) 0.5+dim(2) 0.5
  0.5        0.5+dim(2) 0.5
  0.5        0.5        0.5+dim(3)
  0.5+dim(1) 0.5        0.5+dim(3)
  0.5+dim(1) 0.5+dim(2) 0.5+dim(3)
  0.5        0.5+dim(2) 0.5+dim(3)
  ];
corner_hc = ft_warp_apply(transform, corner_vc);

if isempty(unit)
  if ~isequal(transform, eye(4))
    % estimate the geometrical units we are dealing with
    unit = ft_estimate_units(norm(range(corner_hc)));
  else
    % units are in voxels, these are assumed to be close to mm
    unit = 'mm';
  end
end
if isempty(resolution)
  % the default resolution is 1 mm
  resolution = ft_scalingfactor('mm', unit);
end

% determine whether interpolation is needed
dointerp = false;
dointerp = dointerp || sum(sum(transform-eye(4)))~=0;
dointerp = dointerp || ~all(round(loc)==loc);
dointerp = dointerp || sum(ori)~=1;
dointerp = dointerp || ~(resolution==round(resolution));
% determine the caller function and toggle dointerp to true, if ft_plot_slice has been called from ft_plot_montage
% this is necessary for the correct allocation of the persistent variables
st = dbstack;
if ~dointerp && numel(st)>1 && strcmp(st(2).name, 'ft_plot_montage'), dointerp = true; end

% define 'x' and 'y' axis in projection plane, the definition of x and y is more or less arbitrary
[x, y] = projplane(ori);
% z = ori;

% project the corner points onto the projection plane
corner_pc = zeros(size(corner_hc));
for i=1:8
  corner   = corner_hc(i, :) - loc(:)';
  corner_pc(i,1) = dot(corner, x);
  corner_pc(i,2) = dot(corner, y);
  corner_pc(i,3) = 0;
end

% get the transformation matrix from the projection plane to head coordinates
T2 = [x(:) y(:) ori(:) loc(:); 0 0 0 1];

% get the transformation matrix from projection plane to voxel coordinates
T3 = transform\T2;

min_corner_pc = min(corner_pc, [], 1);
max_corner_pc = max(corner_pc, [], 1);
% round the bounding box limits to the nearest mm
switch unit
  case 'm'
    min_corner_pc = ceil(min_corner_pc*100)/100;
    max_corner_pc = floor(max_corner_pc*100)/100;
  case 'cm'
    min_corner_pc = ceil(min_corner_pc*10)/10;
    max_corner_pc = floor(max_corner_pc*10)/10;
  case 'mm'
    min_corner_pc = ceil(min_corner_pc);
    max_corner_pc = floor(max_corner_pc);
end

% determine a grid of points in the projection plane
xplane = min_corner_pc(1):resolution:max_corner_pc(1);
yplane = min_corner_pc(2):resolution:max_corner_pc(2);
zplane = 0;
[Xi, Yi, Zi]      = ndgrid(xplane, yplane, zplane);
siz               = [size(squeeze(Xi)) size(dat,4)];
interp_center_pc  = [Xi(:) Yi(:) Zi(:)];
% interp_center_hc = ft_warp_apply(T2, interp_center_pc);

% get the positions of the points in the projection plane in voxel coordinates
interp_center_vc = ft_warp_apply(T3, interp_center_pc);

Xi = reshape(interp_center_vc(:, 1), siz(1:2));
Yi = reshape(interp_center_vc(:, 2), siz(1:2));
Zi = reshape(interp_center_vc(:, 3), siz(1:2));

% check whether the values in the axes are close enough to integer
tol = nanmean([diff(unique(Xi(:)));diff(unique(Yi(:)))])./100;
isintegerXi = issufficientlyinteger(Xi(:),tol);
isintegerYi = issufficientlyinteger(Yi(:),tol);
isintegerZi = issufficientlyinteger(Zi(:),tol);

% check whether it's possible to select an orthogonal plane
[islineXi, lineXi] = isline(Xi);
[islineYi, lineYi] = isline(Yi);
[islineZi, lineZi] = isline(Zi);

use_interpn = ~isequal(transform, eye(4)) || ~isequal(interpmethod, 'nearest') || ~all([isintegerXi isintegerYi isintegerZi]);
get_slice   = ~use_interpn && all([islineXi islineYi islineZi]) && all([isintegerXi isintegerYi isintegerZi]);
if use_interpn
  V  = interpn(X, Y, Z, dat, Xi, Yi, Zi, interpmethod);
  if domask,       Vmask = interpn(X, Y, Z, datmask,       Xi, Yi, Zi, interpmethod); end
  if dobackground, Vback = interpn(X, Y, Z, background, Xi, Yi, Zi, interpmethod); end
elseif get_slice
  %something more efficient than an interpolation can be done
  % just select the appropriate plane, and permute to get the orientation
  % right in the plots, something to do with ndgrid vs meshgrid I think
  permutevec = [2 1];
  if ndims(dat)>3
    permutevec = [permutevec 3:ndims(dat)];
  end
  if numel(unique(lineXi(:)))==1
    lineXi = lineXi(1);
  elseif numel(unique(lineYi(:)))==1
    lineYi = lineYi(1);
  elseif numel(unique(lineZi(:)))==1
    lineZi = lineZi(1);
  end
  V = permute(reshape(dat(lineXi,lineYi,lineZi,:), siz(permutevec(1:ndims(dat)-1))), permutevec);
  if domask,       Vmask = permute(reshape(datmask(lineXi,lineYi,lineZi,:),       siz(permutevec(1:2))), [2 1]); end
  if dobackground, Vback = permute(reshape(background(lineXi,lineYi,lineZi,:), siz(permutevec(1:2))), [2 1]); end
else
  % use sub2ind in the unlikely case that it's an oblique plane, parallel
  % to one of the axes with only integer indices
  % this fails for rgb data
  V = dat(sub2ind(dim, Xi(:), Yi(:), Zi(:)));
  V = reshape(V, siz);
end

if all(isnan(V(:)))
  % the projection plane lies completely outside the box spanned by the data
else
  % trim the edges of the projection plane
  [sel1, sel2] = tight(V(:,:,1));
  V  = V (sel1,sel2,:);
  Xi = Xi(sel1,sel2);
  Yi = Yi(sel1,sel2);
  Zi = Zi(sel1,sel2);
  if domask
    Vmask = Vmask(sel1,sel2);
  end
  if dobackground
    Vback = Vback(sel1,sel2);
  end
end

if dobackground
  % convert the background plane to a grayscale image
  bmin  = nanmin(background(:));
  bmax  = nanmax(background(:));
  Vback = (Vback-bmin)./(bmax-bmin);
  Vback(~isfinite(Vback)) = 0;
  Vback = cat(3, Vback, Vback, Vback);
end

interp_center_vc = [Xi(:) Yi(:) Zi(:)]; clear Xi Yi Zi
interp_center_pc = ft_warp_apply(inv(T3), interp_center_vc);

% determine a grid of points in the projection plane
% this reconstruction is needed since the edges may have been trimmed off
xplane = min(interp_center_pc(:, 1)):resolution:max(interp_center_pc(:, 1));
yplane = min(interp_center_pc(:, 2)):resolution:max(interp_center_pc(:, 2));
zplane = 0;

[Xi, Yi, Zi] = ndgrid(xplane, yplane, zplane); % 2D cartesian grid of projection plane in plane voxels
siz          = size(squeeze(Xi));

% extend with one voxel along dim 1
Xi = cat(1, Xi, Xi(end,:)+mean(diff(Xi,[],1),1));
Yi = cat(1, Yi, Yi(end,:)+mean(diff(Yi,[],1),1));
Zi = cat(1, Zi, Zi(end,:)+mean(diff(Zi,[],1),1));
% extend with one voxel along dim 2
Xi = cat(2, Xi, Xi(:,end)+mean(diff(Xi,[],2),2));
Yi = cat(2, Yi, Yi(:,end)+mean(diff(Yi,[],2),2));
Zi = cat(2, Zi, Zi(:,end)+mean(diff(Zi,[],2),2));
% shift with half a voxel along dim 1 and 2
Xi = Xi-0.5*resolution;
Yi = Yi-0.5*resolution;
% Zi = Zi; % do not shift along this direction

interp_edge_pc = [Xi(:) Yi(:) Zi(:)]; clear Xi Yi Zi
interp_edge_hc = ft_warp_apply(T2, interp_edge_pc);

if false
  % plot all objects in head coordinates
  ft_plot_mesh(voxel_center_hc,   'vertexmarker', 'o')
  ft_plot_mesh(voxel_edge_hc,     'vertexmarker', '+')
  ft_plot_mesh(corner_hc,         'vertexmarker', '*')
  ft_plot_mesh(interp_center_hc,  'vertexmarker', 'o', 'vertexcolor', 'r')
  ft_plot_mesh(interp_edge_hc,    'vertexmarker', '+', 'vertexcolor', 'r')
  axis on
  grid on
  xlabel('x')
  ylabel('y')
  zlabel('z')
end

if isempty(cmap)
  % treat as gray value: scale and convert to rgb
  if doscale
    dmin = min(dat(:));
    dmax = max(dat(:));
    V    = (V-dmin)./(dmax-dmin);
    clear dmin dmax
  end
  V(~isfinite(V)) = 0;

  % deal with clim for RGB data here, where the purpose is to increase the
  % contrast range, rather than shift the average grey value
  if ~isempty(clim)
    V = (V-clim(1))./clim(2);
    V(V>1)=1;
  end

  % convert into RGB values, e.g. for the plotting of anatomy
  V = cat(3, V, V, V);
end

% get positions of the voxels in the interpolation plane in head coordinates
Xh = reshape(interp_edge_hc(:,1), siz+1);
Yh = reshape(interp_edge_hc(:,2), siz+1);
Zh = reshape(interp_edge_hc(:,3), siz+1);

% do the actual plotting of the slice
if ~domask
  % no masked slice to be plotted
  if isempty(h)
    % create surface object
    h = surface(Xh, Yh, Zh, V);
    set(h, 'linestyle', 'none');
  else
    % update the colordata in the surface object
    set(h, 'Cdata', V);
    set(h, 'Xdata', Xh);
    set(h, 'Ydata', Yh);
    set(h, 'Zdata', Zh);
  end
elseif domask
  % what should be done depends on the maskstyle
  switch maskstyle
    case 'opacity'
      if dobackground
        ft_warning('specifying maskstyle = ''opacity'' causes the supplied background image not to be used');
      end
      if isempty(h)
        % create surface object
        h = surface(Xh, Yh, Zh, V);
        set(h, 'linestyle', 'none');
      else
        % update the colordata in the surface object
        set(h, 'Cdata', V);
        set(h, 'Xdata', Xh);
        set(h, 'Ydata', Yh);
        set(h, 'Zdata', Zh);
      end
      if islogical(Vmask), Vmask = double(Vmask); end
      set(h, 'FaceColor', 'texture');
      set(h, 'FaceAlpha', 'texturemap'); %flat
      set(h, 'AlphaDataMapping', 'scaled');
      set(h, 'AlphaData', Vmask);
      if ~isempty(opacitylim)
        alim(opacitylim)
      end

    case 'colormix'
      if isempty(cmap), error('using ''colormix'' as maskstyle requires an explicitly defined colormap'); end
      V = bg_rgba2rgb(Vback,V,cmap,clim,Vmask,'rampup',opacitylim);
      if isempty(h)
        % create surface object
        h = surface(Xh, Yh, Zh, V);
        set(h, 'linestyle', 'none');
      else
        % update the colordata in the surface object
        set(h, 'Cdata', V);
        set(h, 'Xdata', Xh);
        set(h, 'Ydata', Yh);
        set(h, 'Zdata', Zh);
      end
  otherwise
    error('unsupported maskstyle');
  end
end

% plot the intersection with a mesh
if dointersect
  % determine three points on the plane
  inplane = eye(3) - (eye(3) * ori') * ori;
  v1 = loc + inplane(1,:);
  v2 = loc + inplane(2,:);
  v3 = loc + inplane(3,:);

  for k = 1:numel(mesh)
    [xmesh, ymesh, zmesh] = intersect_plane(mesh{k}.pos, mesh{k}.tri, v1, v2, v3);

    % draw each individual line segment of the intersection
    if ~isempty(xmesh)
      p = patch(xmesh', ymesh', zmesh', nan(1, size(xmesh,1)));
      if ~isempty(intersectcolor),     set(p, 'EdgeColor', intersectcolor(k)); end
      if ~isempty(intersectlinewidth), set(p, 'LineWidth', intersectlinewidth); end
      if ~isempty(intersectlinestyle), set(p, 'LineStyle', intersectlinestyle); end
    end
  end
end

if ~isempty(cmap) && ~isequal(cmap, 'rgb')
  colormap(cmap);
  if ~isempty(clim)
    caxis(clim);
  end
end

if ~isempty(plotmarker)
  % determine three points on the plane
  inplane = eye(3) - (eye(3) * ori') * ori;
  v1 = loc + inplane(1,:);
  v2 = loc + inplane(2,:);
  v3 = loc + inplane(3,:);
  pr = nan(size(plotmarker,1), 3);
  d  = nan(size(plotmarker,1), 1);
  for k = 1:size(plotmarker,1)
    [pr(k,:), d(k,:)] = ptriprojn(v1, v2, v3, plotmarker(k,:));
  end
  sel = d<eps*1e8;
  if sum(sel)>0
    ft_plot_dipole(pr(sel,:), repmat([0;0;1], 1, size(pr,1)), 'length', 0, 'color', markercolor, 'diameter', markersize);
  end
end

% update the axes to ensure that the whole volume fits
ax = [min(corner_hc) max(corner_hc)];
axis(ax([1 4 2 5 3 6])); % reorder into [xmin xmax ymin ymaz zmin zmax]

st = dbstack;
if numel(st)>1
  % ft_plot_slice has been called from another function
  % assume the remainder of the axis settings to be handled there
else
  set(gca,'xlim',[min(Xh(:))-0.5*resolution max(Xh(:))+0.5*resolution]);
  set(gca,'ylim',[min(Yh(:))-0.5*resolution max(Yh(:))+0.5*resolution]);
  set(gca,'zlim',[min(Zh(:))-0.5*resolution max(Zh(:))+0.5*resolution]);

  set(gca,'dataaspectratio',[1 1 1]);
  % axis equal; % this for some reason does not work robustly when drawing intersections, replaced by the above
  axis vis3d
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, y] = projplane(z)
[u, s, v] = svd([eye(3) z(:)]);
x = u(:, 2)';
y = u(:, 3)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [sel1, sel2] = tight(V)
% make a selection to cut off the nans at the edges
sel1 = sum(~isfinite(V), 2)<size(V, 2);
sel2 = sum(~isfinite(V), 1)<size(V, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bool = issufficientlyinteger(X, tolerance)
%isinteger only checks for integer class, so will always return false with
%double integers
bool = all(abs(X-round(X))<tolerance);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bool, lineX] = isline(X)
%isline returns an array of the values in X are columnwise or rowwise the
%same, otherwise returns false
if isequal(X(ones(1,size(X,1)),:),X)
  lineX = X(1,:);
  bool = true;
elseif isequal(X(:,ones(1,size(X,2))),X)
  lineX = X(:,1)';
  bool = true;
else
  lineX = [];
  bool = false;
end

author: spm


text:
function uTest_FiltFiltM(doSpeed)
% Unit test: FiltFiltM
% This is a routine for automatic testing. It is not needed for processing and
% can be deleted or moved to a folder, where it does not bother.
%
% uTest_FiltFiltM(doSpeed)
% INPUT:
%   doSpeed: If ANY(doSpeed) is TRUE, the speed is measured. Optional.
% OUTPUT:
%   On failure the test stops with an error.
%
% Tested: Matlab 6.5, 7.7, 7.8, WinXP
% Author: Jan Simon, Heidelberg, (C) 2009-2011 matlab.THISYEAR(a)nMINUSsimon.de

% $JRev: R-k V:010 Sum:VKPp090JCPL7 Date:19-Jul-2011 19:40:07 $
% $License: BSD (use/copy/change/redistribute on own risk, mention the author) $
% $File: Tools\UnitTests_\uTest_FiltFiltM.m $

% Init: ========================================================================
ErrHead = ['JSimon:', mfilename, ':TestFailed'];

if nargin == 0 || any(doSpeed)
   TestTime = 2.0;  % sec
else
   TestTime = 0.5;  % sec
end

% ==============================================================================
disp(['==== Test FiltFiltM  ', datestr(now, 0)]);
disp(['  Version: ', which('FiltFiltM')]);
pause(0.01);  % Flush events

fprintf('\nKnown answer tests:\n');
B = [];
A = [];
Y = FiltFiltM(B, A, []);
if ~isempty(Y)
   error(ErrHead, ['*** ', mfilename, ': Bad reply for empty signal']);
end
disp('  ok: all input empty');

B = [];
A = [];
X = rand(1000, 1);
Y = FiltFiltM(B, A, X);
if ~isempty(Y)
   error(ErrHead, ['*** ', mfilename, ': Bad reply for empty parameters']);
end
disp('  ok: empty parameters');

% [B, A] = ButterParam(2, 0.2, 'low')
B = [0.06745527, 0.1349105, 0.06745527];
A = [1, -1.142981, 0.4128016];
Y = FiltFiltM(B, A, []);
if ~isempty(Y)
   error(ErrHead, ...
      ['*** ', mfilename, ': Bad reply for empty signal with parameters']);
end
disp('  ok: empty signal');

Y = FiltFiltM(B, A, X);
Z = filtfilt(B, A, X);
if ~isequal(size(Y), size(Z)) || any(abs(Y - Z) > 10 * eps)
   error(ErrHead, ...
      ['*** ', mfilename, ': Bad result for [1000 x 1]']);
end
disp('  ok: [N x 1] vector');

Xt = reshape(X, 1, []);
Y  = FiltFiltM(B, A, Xt);
if ~isequal(size(Y), size(Xt)) || any(abs(Y(:) - Z(:)) > 10 * eps)
   error(ErrHead, ...
      ['*** ', mfilename, ': Bad result for [1 x 1000]']);
end
disp('  ok: [1 x N] vector');

X = [X, X];
Y = FiltFiltM(B, A, X);
Z = filtfilt(B, A, X);
if ~isequal(size(Y), size(Z)) || any(abs(Y(:) - Z(:)) > 10 * eps)
   error(ErrHead, ...
      ['*** ', mfilename, ': Bad result for [1000 x 2]']);
end
disp('  ok: 2D array, 1st dim');

X = transpose(X);
Y = FiltFiltM(B, A, X, 2);
Z = transpose(Z);
if ~isequal(size(Y), size(Z)) || any(abs(Y(:) - Z(:)) > 10 * eps)
   error(ErrHead, ...
      ['*** ', mfilename, ': Bad result for [2 x 1000], 2nd dim']);
end
disp('  ok: 2D array, 2nd dim');

X = rand(20, 30, 40);
Y = FiltFiltM(B, A, X);
Z = filtfilt(B, A, X);
if ~isequal(size(Y), size(Z)) || any(abs(Y(:) - Z(:)) > 10 * eps)
   error(ErrHead, ...
      ['*** ', mfilename, ': Bad result for [20 x 30 x 40]']);
end

Y = FiltFiltM(B, A, X, 2);
Z = ipermute(filtfilt(B, A, permute(X, [2,1,3])), [2,1,3]);
if ~isequal(size(Y), size(Z)) || any(abs(Y(:) - Z(:)) > 10 * eps)
   error(ErrHead, ...
      ['*** ', mfilename, ': Bad result for [20 x 30 x 40], 2nd dim']);
end

Y = FiltFiltM(B, A, X, 3);
Z = ipermute(filtfilt(B, A, permute(X, [3,1,2])), [3,1,2]);
if ~isequal(size(Y), size(Z)) || any(abs(Y(:) - Z(:)) > 10 * eps)
   error(ErrHead, ...
      ['*** ', mfilename, ': Bad result for [20 x 30 x 40], 3rd dim']);
end
disp('  ok: 3D array, 1st, 2nd, 3rd dim');

X = rand(1000, 1);
Y = FiltFiltM(B, A, X);
Z = FiltFiltM(B, A, reshape(X, 1, 1, []), []);
if ~isequal(size(Z), [1, 1, numel(X)]) || any(abs(Y(:) - Z(:)) > 10 * eps)
   error(ErrHead, ['*** ', mfilename, ...
      ': Bad result for [1 x 1 x 1000], 1st non-singelton']);
end
disp('  ok: [1 x 1 x N] array, 1st non-singelton');

% FILTER does not support SINGLES in Matlab 6.5:
try
   y = filter(B, A, single(ones(1, 100)));   %#ok<NASGU>
   FILTERforSINGLE = true;
catch
   FILTERforSINGLE = false;
   lasterr('');
end

if FILTERforSINGLE
   X  = rand(1000, 2);
   Y  = FiltFiltM(B, A, X);
   Xs = single(X);
   Ys = FiltFiltM(B, A, Xs);
   Zs = filtfilt(B, A, Xs);
   if ~isequal(size(Ys), size(Zs))
      error(ErrHead, ['*** ', mfilename, ': Bad result for SINGLE']);
   end
   epsSingle    = 1.192092895507813e-007;
   errFiltFiltM = max(abs(double(Ys(:)) - Y(:))) / epsSingle;
   errFILTFILT  = max(abs(double(Zs(:)) - Y(:))) / epsSingle;
   if any(errFiltFiltM > 10)
      error(ErrHead, ['*** ', mfilename, ': Bad result for SINGLE']);
   end
   disp('  ok: SINGLE array');
   fprintf(['      Difference to DOUBLE divided by EPS(SINGLE(1)):\n', ...
      '      FilterM: %g, FILTER: %g\n'], errFiltFiltM, errFILTFILT);

else
   v = sscanf(version, '%d.', 2);
   fprintf('  ok: FILTER does not accept SINGLE in Matlab %d.%d\n', v);
end

% ------------------------------------------------------------------------------
tooLazy = false;
fprintf('\nProvoke problems:\n');
try
   Y = FiltFiltM(B, A, rand(4));  %#ok<NASGU>
   tooLazy = true;
catch
   disp('  ok: Too short signal in dim 1 caught.');
end
if tooLazy
   error(ErrHead, ['*** ', mfilename, ': Too short signal accepted.']);
end

try
   Y = FiltFiltM(B, A, rand(4), 2);  %#ok<NASGU>
   tooLazy = true;
catch
   disp('  ok: Too short signal in dim 2 caught.');
end
if tooLazy
   error(ErrHead, ['*** ', mfilename, ': Too short signal accepted.']);
end

try
   Y = FiltFiltM(B, A, rand(1, 100), 3);  %#ok<NASGU>
   tooLazy = true;
catch
   disp('  ok: Dim out of range caught.');
end
if tooLazy
   error(ErrHead, ['*** ', mfilename, ': Dim out of range accepted.']);
end

try
   Y = FiltFiltM(B, A, rand(1, 100), -3.1);  %#ok<NASGU>
   tooLazy = true;
catch
   disp('  ok: Not integer Dim caught.');
end
if tooLazy
   error(ErrHead, ['*** ', mfilename, ': Not integer Dim accepted.']);
end

% Time test: -------------------------------------------------------------------
% Results on a Pentium-M 1500 MHz, 512MB RAM, WinXP, Matlab 2009a:
% Number of loops to get about 2 seconds run time
% 1000000 x 1:  6 loops,  FILTFILT: 1.930 sec,  FiltFiltM: 0.385 sec  ==> 20.0%
% 100000 x 1:  70 loops,  FILTFILT: 2.077 sec,  FiltFiltM: 0.488 sec  ==> 23.5%
% 10000 x 10:  82 loops,  FILTFILT: 1.996 sec,  FiltFiltM: 0.600 sec  ==> 30.0%
% 1000 x 100:  50 loops,  FILTFILT: 2.761 sec,  FiltFiltM: 0.518 sec  ==> 18.7%
% 100 x 100:  100 loops,  FILTFILT: 1.935 sec,  FiltFiltM: 0.079 sec  ==>  4.1%
% 100 x 1000:  12 loops,  FILTFILT: 2.303 sec,  FiltFiltM: 0.130 sec  ==>  5.6%

% A 4th order Butterworth filter: [B, A] = Butter(4, 0.25, 'low')
B = [0.0102094807912032, 0.0408379231648126, 0.0612568847472189, ...
   0.0408379231648126, 0.0102094807912032];
A = [1.0, -1.96842778693852, 1.73586070920889, -0.724470829507363, ...
   0.120389599896245];
fprintf('\n== Speed test with Butterworth filter of order %d:\n', ...
   length(B) - 1);

fprintf('Number of loops determined for about %g seconds run time\n', TestTime);

TestSize = {[1e6, 1], [1e5, 1], [1e4, 10], [1e3, 100], [100, 100], [100, 1e3]};
for iTest = 1:length(TestSize)
   % Create test data:
   aSize = TestSize{iTest};
   x     = rand(aSize(1), aSize(2));
   msg   = sprintf('  %d x %d:', aSize);
   fprintf('%-14s', msg);
   drawnow;

   % Run test of original function:
   iTime = cputime;
   nLoop = 0;
   while cputime - iTime < 1
      nLoop = nLoop + 1;
      origy = filtfilt(B, A, x);  %#ok<NASGU>
      clear('origy');
   end
   nLoop = max(2, ceil(nLoop * TestTime));
   fprintf('%5d loops\n', nLoop);

   % Run test for vectorized function:
   tic;
   for i = 1:nLoop
      origy = filtfilt(B, A, x);  %#ok<NASGU>
      clear('origy');
   end
   elapsed0 = toc;
   fprintf('    FILTFILT:  %8.3f sec\n', elapsed0);

   % Run test for vectorized function:
   tic;
   for i = 1:nLoop
      newy = FiltFiltM(B, A, x);  %#ok<NASGU>
      clear('newy');
   end
   elapsed1 = toc;
   fprintf('    FiltFiltM: %8.3f sec', elapsed1);
   fprintf('  ==> %.1f%%\n', 100 * elapsed1 / elapsed0);

   % Compare the results - usually the results are identical, but small round
   % off errors would be ok:
   origy = filtfilt(B, A, x);
   newy  = FiltFiltM(B, A, x);
   if isequal(origy, newy)
      fprintf('    Results are identical\n\n');
   elseif all(abs(origy(:) - newy(:)) < 10 * eps)
      fprintf('    Results are equal within 10*eps\n\n');
   else
      error(ErrHead, 'Different results - DO NOT USE FiltFiltM!');
   end
end

% Bye:
disp('== FiltFiltM passed the tests.');

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% MAPTEXT   Text annotations on maps.
%
%    MAPTEXT(C, 'string') adds the text in the quotes to location
%    specified by the coordinates C. C(1,:) specifies the longitudes and
%    C(2,:) the latitudes.
%
%    MAPTEXT(C, S) adds the texts in the cell array S to locations specified
%    by the coordinates C. C(1,:) specifies the longitudes and C(2,:) the
%    latitudes.
%
%    MAPTEXT(..., property/value pairs) can be used to change, e.g.,
%    fontsize, weight, and color using the standard TEXT properties.
%
%    See also TEXT and M_TEXT.

% Copyright (c) 2010 Jaakko Luttinen

function maptext(C, str, varargin)

if iscell(str)
  for i=1:numel(str)
    m_text(C(1,i), C(2,i), str{i}, varargin{:});
  end
else
  m_text(C(1,:), C(2,:), str, varargin{:});
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [data, wAxis, dAxis, xAxis, yAxis, misc] = lscload(filename)
% Reads in IR image data from PerkinElmer block structured files.
% This version is compatible with '1994' standard LSC files.
%
% [data, wAxis, dAxis, xAxis, yAxis, misc] = lscload(filename):
%   data:  2D array, [W, D]     wavelength x distance
%   wAxis: vector of e.g. wavenumbers
%   dAxis: vector of distance along line
%   xAxis: vector for stage x positions (e.g. micrometers)
%   yAxis: vector for stage y positions (e.g. micrometers)
%   misc: miscellanous information in name,value pairs

% Copyright (C)2007 PerkinElmer Life and Analytical Sciences
% Stephen Westlake, Seer Green
%
% History
% 2007-05-01 SW     Initial version

% Block IDs
DS4C3IData            = 5100;   % 4D DS: header info
DS4C3IPts             = 5101;   % 4D DS: point coordinates (CvCoOrdArray)
DS4C3IPtsCont         = 5102;   % 4D DS: point coordinates (continuator for large CvCoOrdArray)
DS4C3ITemp2DData      = 5103;   % 4D DS: temp block for storing part of component 2d dataset
DS4C3I2DData          = 5104;   % 4D DS: release sw block for storing ALL of component 2d data
DS4C3IFloatPts        = 5105;   % 4D DS: point coordinates (FloatArray)
DS4C3IFloatPtsCont    = 5106;   % 4D DS: point coordinates (continuator for large FloatArray)
DS4C3IHistory         = 5107;   % 4D DS: history record

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

% Fixed file header of signature and description
signature = setstr(fread(fid, 4, 'uchar')');
if ~strcmp(signature, 'PEPE')
    error('This is not a PerkinElmer block structured file.');
    return
end
description = setstr(fread(fid, 40, 'uchar')');

% Initialize a variable so we can tell if we have read it.
dLen = int32(0);
qLen = int32(0);

% The rest of the file is a list of blocks
while ~feof(fid)
    blockID = fread(fid,1,'int16');
    blockSize = fread(fid,1,'int32');

    % feof does not go true until after the read has failed.
    if feof(fid)
        break
    end

    switch blockID
        case DS4C3IData
            len = fread(fid,1,'int16');
            alias = setstr(fread(fid, len, 'uchar')');

            dDelta = fread(fid, 1, 'double');       % distance delta
            angle = fread(fid, 1, 'double');        % angle (radians)
            wDelta = fread(fid, 1, 'double');       % wavenumber delta
            fread(fid, 1, 'double');      % zstart
            fread(fid, 1, 'double');      % zend
            fread(fid, 1, 'double');      % min data value
            fread(fid, 1, 'double');      % max data value
            x0 = fread(fid, 1, 'double');           % stage x0
            y0 = fread(fid, 1, 'double');           % stage y0
            w0 = fread(fid, 1, 'double');           % wavenumber start
            dLen = fread(fid, 1, 'int32');          % distance points
            qLen = fread(fid, 1, 'int32');          % should be 1
            wLen = fread(fid, 1, 'int32');          % wavelength points

            len = fread(fid, 1, 'int16');
            dLabel = setstr(fread(fid, len, 'uchar')');
            len = fread(fid, 1, 'int16');
            setstr(fread(fid, len, 'uchar')');
            len = fread(fid, 1, 'int16');
            wLabel = setstr(fread(fid, len, 'uchar')');
            len = fread(fid, 1, 'int16');
            zLabel = setstr(fread(fid, len, 'uchar')');

            % imshow row 1 is the top, i.e. fsm row N.
            d = dLen;

        case DS4C3IFloatPts         % the next spectrum
            data(d, :) = fread(fid, wLen, 'float');
            if d == 0
                error('There are too many data points.');
                return
            end
            % set up the index for the next point
            d = d - 1;

        otherwise               % unknown block, just seek past it
            fseek(fid, blockSize, 'cof');
    end
end
fclose(fid);

if dLen == 0
    error('The file does not contain spectral image data.');
    return
end
if qLen ~= 1
    error('The file does not contain line scan data.');
    return
end

% d0 is implicitly 0
dEnd = (dLen - 1) * dDelta;
wEnd = w0 + (wLen - 1) * wDelta;

% Calculate end stage positions from angle and hypoteneuse
xEnd = x0 + dEnd * cos(angle);
yEnd = y0 + dEnd * sin(angle);
xDelta = dDelta * cos(angle);
yDelta = dDelta * sin(angle);

% Expand the axes specifications into vectors
% D axis is reversed to match the image data
xAxis = x0 : xDelta : xEnd;
yAxis = y0 : yDelta : yEnd;
dAxis = dEnd : -dDelta : 0;
wAxis = w0 : wDelta : wEnd;

% Return the other details as name,value pairs
misc(1,:) = {'dLabel', dLabel};
misc(2,:) = {'wLabel', wLabel};
misc(3,:) = {'zLabel', zLabel};
misc(4,:) = {'alias', alias};
misc(5,:) = {'angle', angle};

author: Sable

repo: mcbench-benchmarks


text:
function [Population,FrontNo,DWeight] = EnvironmentalSelection(Population,N)
% The environmental selection of DWU

%------------------------------- 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".
%--------------------------------------------------------------------------

% This function is written by Gladston Moreira

    %% Non-dominated sorting
    [FrontNo] = NDSort(Population.objs,Population.cons,N);

    %% Calculate the dominance information each solution
    DWeight = InfoDominance(Population.objs);

    %% Environment Select
    Next = ReplacementUniformity(Population,N,FrontNo,DWeight);

    %% Population for next generation
    Population = Population(Next);
    FrontNo    = FrontNo(Next);
    DWeight    = DWeight(Next);
end

function InfoD = InfoDominance(PopObj)
% Calculate the information dominance each solution

    N = size(PopObj,1);

    %% Dominance count each solution
    D = false(N);
    for i = 1 : N-1
        for j = i+1 : N
            k = any(PopObj(i,:)<PopObj(j,:)) - any(PopObj(i,:)>PopObj(j,:));
            if k == 1
                D(i,j) = true;
            elseif k == -1
                D(j,i) = true;
            end
        end
    end
    CountDominance = sum(D,2);

    %% Calculate information dominance each solution
    InfoD = D'*CountDominance;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
 function ob = Gtomo2_wtmex(arg1, varargin)
%function ob = Gtomo2_wtmex(sg, ig, options)
%function ob = Gtomo2_wtmex(file, options)
%function ob = Gtomo2_wtmex(arg_pairs, options)
%
% Generate a tomographic system model based on aspire's wtfmex() routine.
%
% in
%	sg	strum		sino_geom()
%	ig	strum		image_geom()
%
% options
%	'grouped'		row or col (default: row)
%	'nthread'		pthreads for multiple processors (default: 1)
%	'mask'	[nx,ny]		logical support mask, has precedence
%	'pairs' {}		cell array of name/value pairs for aspire_pair()
%
% out
%	ob	[nb*na,np]	Fatrix system "matrix", np = sum(ig.mask(:))
%
% Copyright 2006-3-3, Jeff Fessler, The University of Michigan

if nargin == 1 && streq(arg1, 'test'), Gtomo2_wtmex_test, return, end
if nargin < 1, ir_usage, end

% options
arg.grouped = 'row';
arg.nthread = 1;
arg.mask = [];
arg.chat = 0;
arg.pairs = {};

% given wtf or argument pairs
if ischar(arg1)

	arg = vararg_pair(arg, varargin);

	if exist(arg1, 'file') % file.wtf
		if ~isempty(arg.mask), error 'no mask option given .wtf', end
		arg.file = arg1;
		[nx ny nb na] = wtfmex('read', arg.file);

	else % arg_pair
		arg.arg = arg1;
		if isempty(arg.mask)
			mask_arg = {};
		else
			if ~islogical(arg.mask), error 'need logical mask', end
			mask_arg = uint8(arg.mask);
		end

		wtfmex('gensys', arg.arg', arg.grouped, mask_arg{:});

		nx = str2num(arg_get(arg.arg, 'nx'));
		ny = str2num(arg_get(arg.arg, 'ny'));
		nb = str2num(arg_get(arg.arg, 'nb'));
		na = str2num(arg_get(arg.arg, 'na'));
	end

	arg.ig = image_geom('nx', nx, 'ny', ny, 'dx', -1, ...
		'mask', wtfmex('mask') > 0);
	arg.sg = sino_geom('par', 'nb', nb, 'na', na);

% given sino_geom() and image_geom()
else

	arg.sg = arg1;
	arg.ig = varargin{1};
	varargin = {varargin{2:end}};
	arg = vararg_pair(arg, varargin);
	if ~isempty(arg.mask)
		arg.ig.mask = arg.mask;
	end

	arg.aspire_arg = aspire_pair(arg.sg, arg.ig, 'support', 'array', ...
		arg.pairs{:}); % trick: pairs can override default support

	wtfmex('gensys', arg.aspire_arg', arg.grouped, uint8(arg.ig.mask));

end

arg.power = 1;
arg.nd = arg.sg.nb * arg.sg.na;
arg.np = sum(arg.ig.mask(:));
dim = [arg.nd arg.np]; % trick: make it masked by default!

%
% build Fatrix object
%
ob = Fatrix(dim, arg, 'caller', 'Gtomo2_wtmex', ...
	'forw', @Gtomo2_wtmex_forw, 'back', @Gtomo2_wtmex_back, ...
	'free', @Gtomo2_wtmex_free, ...
	'mtimes_block', @Gtomo2_wtmex_mtimes_block);

if arg.chat
	wtfmex('print');
end

%
% Gtomo2_wtmex_forw(): y = G * x
%
function y = Gtomo2_wtmex_forw(arg, x)
y = Gtomo2_wtmex_mtimes_block(arg, 0, x, 1, 1);

%
% Gtomo2_wtmex_back(): x = G' * y
% full backprojection
%
function x = Gtomo2_wtmex_back(arg, y)
if arg.power == 2
	[flag_column y] = Gtomo2_wtmex_y_shape(arg, y, istart, nblock);
	x = wtfmex('chat', arg.chat, 'back2', single(y));
	if flag_column, x = arg.ig.maskit(x); end
else
	x = Gtomo2_wtmex_mtimes_block(arg, 1, y, 1, 1);
end

%
% Gtomo2_wtmex_mtimes_block()
%
function y = Gtomo2_wtmex_mtimes_block(arg, is_transpose, x, istart, nblock)

if is_transpose
	y = Gtomo2_wtmex_block_back(arg, x, istart, nblock);
else
	y = Gtomo2_wtmex_block_forw(arg, x, istart, nblock);
end

%
% Gtomo2_wtmex_block_forw()
%
function y = Gtomo2_wtmex_block_forw(arg, x, istart, nblock)

if arg.power ~= 1, error('power=%d not done', arg.power), end

% if needed, convert concise column to 3d array
flag_column = 0;
if size(x,1) == arg.np
	flag_column = 1;
	x = arg.ig.embed(x);
end

if nblock == 1
	y = wtfmex('chat', arg.chat, 'mult', single(x));
else
	y = wtfmex('chat', arg.chat, 'proj,block', single(x), ...
			int32(istart-1), int32(nblock));

	% fix: extract the relevant columns - should do in wtfmex?
	ia = istart:nblock:arg.sg.na;
	y = y(:,ia,:);
end

y = double6(y);

if flag_column % column in yields column out.
	y = reshape(y, arg.sg.nb*size(y,2), []);
end

%
% Gtomo2_wtmex_block_back()
%
function x = Gtomo2_wtmex_block_back(arg, y, istart, nblock)

if arg.power ~= 1, error('power=%d not done', arg.power), end

[flag_column y] = Gtomo2_wtmex_y_shape(arg, y, istart, nblock);
if nblock == 1
	x = wtfmex('chat', arg.chat, 'back', single(y));
else
	x = wtfmex('chat', arg.chat, 'back,block', single(y), ...
			int32(istart-1), int32(nblock));
end
x = double6(x);

if flag_column, x = arg.ig.maskit(x); end

%
% Gtomo2_wtmex_y_shape(arg, y, istart, nblock)
%
function [flag_column, y] = Gtomo2_wtmex_y_shape(arg, y, istart, nblock)

flag_column = 0;
nb = arg.sg.nb;
na = arg.sg.na;

if nblock == 1
	if size(y,1) == nb * na % [nb*na,?]
		flag_column = 1;
		y = arg.sg.shape(y); % [nb,na,?]
	end
else
	ia = istart:nblock:na;
	nv = length(ia);
	if size(y,1) == nb * nv % [nb*nv,?]
		flag_column = 1;
		tmp = zeros(nb, na, size(y,2));
		tmp(:,ia,:) = reshape(y, nb, nv, []);
		y = reshape(tmp, nb*na, []); % [nb*na,?]
	elseif size(y,1) == nb
		error todo
	else
		error bug
	end
end

%
% Gtomo2_wtmex_free()
%
function Gtomo2_wtmex_free(arg)
printm('freeing wtfmex static memory')
wtfmex('free');

%
% Gtomo2_wtmex_power()
%
function ob = Gtomo2_wtmex_power(ob, sup)
ob.arg.power = ob.arg.power * sup;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%Documentation
% Calls m2html to generate documentation of m code
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 6/8/2011
classdef Documentation
    methods (Static)
        function generate()
			import edu.stanford.covert.cell.sim.util.Documentation;

            [~, folder] = fileparts(pwd);
			try
				Documentation.cleanup(['doc' filesep 'm2html']);
			end
			if ~exist(['doc' filesep 'm2html'], 'dir')
				mkdir(['doc' filesep 'm2html']);
			end
            directories = [{folder}; cellfun(@(x) [folder filesep x], [
                edu.stanford.covert.util.IOUtil.getDirectoryNamesRecursively('src')], ...
                'UniformOutput', false)];
            cd('..');
            try
                m2html(...
                    'mFiles', directories, ...
                    'htmlDir', [folder filesep 'doc' filesep 'm2html'], ...
                    'recursive', 'off',...
                    'graph', 'on', ...
                    'todo', 'on', ...
                    'verbose', 'off');
				cd(folder);
            catch exception
				cd(folder);
				exception.rethrow();
			end
        end

		function cleanup(folder)
			import edu.stanford.covert.cell.sim.util.Documentation;

			files = dir(folder);
			for i = 1:numel(files)
				if isequal(files(i).name, '.') || isequal(files(i).name, '..') || isequal(files(i).name, '.svn')
					continue;
				end
				if isdir([folder filesep files(i).name])
					Documentation.cleanup([folder filesep files(i).name]);
				else
					delete([folder filesep files(i).name]);
				end
			end
		end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function CPD = discrete_CPD(clamped, dom_sizes)
% DISCRETE_CPD Virtual constructor for generic discrete CPD
% CPD = discrete_CPD(clamped, dom_sizes)

CPD.dom_sizes = dom_sizes;
CPD = class(CPD, 'discrete_CPD', generic_CPD(clamped));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%kvskelthinP 'Skeleton by Thinning Mult. P'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vskelthinP.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image'
% InputFile: str1 'Struct Element 1', required: 'Structuring Element'
% InputFile: str2 'Struct Element 2', required: 'Structuring Element'
% InputFile: str3 'Struct Element 3', required: 'Structuring Element'
% InputFile: str4 'Struct Element 4', required: 'Structuring Element'
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvskelthinP({i, str1, str2, str3, str4}, {'i','';'str1','';'str2','';'str3','';'str4','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vskelthinP - Skeleton by Thinning Mult. P
%
%  DESCRIPTION
% .I vskelthinP
% will do the skeleton by thinning with mult. str. sets of an input image by a set of structuring elements. This program is based on an algorithm proposed by Roland T. Chin and Ben-Kwei Jang.
%
% The output image data type will be the same as that of the highest data type
% of the input image.
%
% Useful set of structuring elements:
%
%                         0 0 0    1 1 0    0 0 0    0 1 0
%   Homotopic Skeleton    0 1 1    1 0 0    1 1 1    0 0 0
%                         0 1 0    0 0 0    0 1 0    0 0 0
%
%
%
%  EXAMPLES
% vskelthinP -i ball.xv -str1 b272 -str2 b273 -str3 b274 -str4 b275 -o outimage.xv
%
% Will do the skeleton by thinning with mult. str. sets of image "ball.xv" by the structuring sets b272, b273, b274 and b275, with the resulting image written to "outimage.xv".
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% .I vskelthinP
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
% The structuring elements are subsets of the 3x3 matrix and the origin is always at the center of this matrix.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvskelthinP(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvskelthinP(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'str1', '__input';'str2', '__input';'str3', '__input';'str4', '__input';'o', '__output'};
maxval={0,0,0,0,0,0};
minval={0,0,0,0,0,0};
istoggle=[0,0,0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','InputFile','InputFile','InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vskelthinP"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef GenericLinearTangentSteeringModel < AbstractSteeringModel
    %GenericLinearTangentSteeringModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        gammaAngleModel(1,1) PolynominalModel = PolynominalModel(0,0,0,0);
        betaAngleModel(1,1) LinearTangentModel = LinearTangentModel(0,0,0,0,0);
        alphaAngleModel(1,1) PolynominalModel = PolynominalModel(0,0,0,0);

        gammaContinuity(1,1) logical = false;
        betaContinuity(1,1) logical = false;
        alphaContinuity(1,1) logical = false;

        refFrame AbstractReferenceFrame
        controlFrame AbstractControlFrame
    end

    methods
        function dcm = getBody2InertialDcmAtTime(obj, ut, rVect, vVect, bodyInfo)
            gammaAng = obj.gammaAngleModel.getValueAtTime(ut);
            betaAng = obj.betaAngleModel.getValueAtTime(ut);
            alphaAng = obj.alphaAngleModel.getValueAtTime(ut);

%             elemSet = CartesianElementSet(ut, rVect(:), vVect(:), bodyInfo.getBodyCenteredInertialFrame());
%             if(not(isempty(obj.refFrame.getOriginBody())))
%                 elemSet = elemSet.convertToFrame(obj.refFrame);
%             end

            dcm = real(obj.controlFrame.computeDcmToInertialFrame(ut, rVect, vVect, bodyInfo, gammaAng, betaAng, alphaAng, obj.refFrame));
        end

        function [angleModel, continuity] = getAngleNModel(obj, n)
            angleModel = PolynominalModel.empty(1,0);

            switch n
                case 1
                    angleModel = obj.gammaAngleModel;
                    continuity = obj.gammaContinuity;
                case 2
                    angleModel = obj.betaAngleModel;
                    continuity = obj.betaContinuity;
                case 3
                    angleModel = obj.alphaAngleModel;
                    continuity = obj.alphaContinuity;
            end
        end

        function t0 = getT0(obj)
            t0 = obj.alphaAngleModel.t0;
        end

        function setT0(obj, newT0)
            obj.gammaAngleModel.t0 = newT0;
            obj.betaAngleModel.t0 = newT0;
            obj.alphaAngleModel.t0 = newT0;
        end

        function setConstTerms(obj, gammaConst, alphaConst)
            obj.gammaAngleModel.constTerm = gammaConst;
            obj.alphaAngleModel.constTerm = alphaConst;
        end

        function setLinearTerms(obj, gammaRate, alphaRate)
            obj.gammaAngleModel.linearTerm = gammaRate;
            obj.alphaAngleModel.linearTerm = alphaRate;
        end

        function setAccelTerms(obj, gammaRateRate, alphaRateRate)
            obj.gammaAngleModel.accelTerm = gammaRateRate;
            obj.alphaAngleModel.accelTerm = alphaRateRate;
        end

        function setLinearTangentTerms(obj, a, a_dot, b, b_dot)
            obj.betaAngleModel.a = a;
            obj.betaAngleModel.a_dot = a_dot;
            obj.betaAngleModel.b = b;
            obj.betaAngleModel.b_dot = b_dot;
        end

        function setTimeOffsets(obj, timeOffset)
            obj.alphaAngleModel.tOffset = timeOffset;
            obj.betaAngleModel.tOffset = timeOffset;
            obj.gammaAngleModel.tOffset = timeOffset;
        end

        function [angle1Cont, angle2Cont, angle3Cont] = getContinuityTerms(obj)
            angle1Cont = obj.gammaContinuity;
            angle2Cont = obj.betaContinuity;
            angle3Cont = obj.alphaContinuity;
        end

        function setContinuityTerms(obj, angle1Cont, angle2Cont, angle3Cont)
            obj.gammaContinuity = angle1Cont;
            obj.betaContinuity = angle2Cont;
            obj.alphaContinuity = angle3Cont;
        end

        function setConstsFromDcmAndContinuitySettings(obj, dcm, ut, rVect, vVect, bodyInfo)
            if(obj.gammaContinuity || obj.betaContinuity || obj.alphaContinuity)
                elemSet = CartesianElementSet(ut, rVect(:), vVect(:), bodyInfo.getBodyCenteredInertialFrame());

%                 if(not(isempty(obj.refFrame.getOriginBody())))
%                     elemSet = elemSet.convertToFrame(obj.refFrame);
%                 end

                [gammaAngle, betaAngle, alphaAngle] = obj.controlFrame.getAnglesFromInertialBodyAxes(LaunchVehicleAttitudeState(dcm), elemSet.time, elemSet.rVect(:), elemSet.vVect(:), elemSet.frame.getOriginBody(), obj.refFrame);

                if(obj.gammaContinuity)
                    obj.gammaAngleModel.constTerm = gammaAngle;
                end

                if(obj.betaContinuity)
                    obj.betaAngleModel.setBForContinuity(betaAngle);
                end

                if(obj.alphaContinuity)
                    obj.alphaAngleModel.constTerm = alphaAngle;
                end
            end
        end

        function setInitialAttitudeFromState(obj, stateLogEntry, tOffsetDelta)
            t0 = stateLogEntry.time;
            obj.setT0(t0);

            obj.gammaAngleModel.tOffset = obj.gammaAngleModel.tOffset + tOffsetDelta;
            obj.betaAngleModel.tOffset = obj.betaAngleModel.tOffset + tOffsetDelta;
            obj.alphaAngleModel.tOffset = obj.alphaAngleModel.tOffset + tOffsetDelta;
        end

        function [angle1Name, angle2Name, angle3Name] = getAngleNames(obj)
            angleNames = obj.controlFrame.getControlFrameEnum().angleNames;
            angle1Name = angleNames{1};
            angle2Name = angleNames{2};
            angle3Name = angleNames{3};
        end

        function newSteeringModel = deepCopy(obj)
            newSteeringModel = GenericPolySteeringModel(obj.gammaAngleModel.deepCopy(), obj.betaAngleModel.deepCopy(), obj.alphaAngleModel.deepCopy());
            newSteeringModel.gammaContinuity = obj.gammaContinuity;
            newSteeringModel.betaContinuity = obj.betaContinuity;
            newSteeringModel.alphaContinuity = obj.alphaContinuity;
            newSteeringModel.refFrame = obj.refFrame;
            newSteeringModel.controlFrame = obj.controlFrame;
        end

        function optVar = getNewOptVar(obj)
            optVar = SetGenericLinearTangentSteeringModelActionOptimVar(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function tf = usesRefFrame(~)
            tf = true;
        end

        function refFrame = getRefFrame(obj)
            refFrame = obj.refFrame;
        end

        function setRefFrame(obj, refFrame)
            obj.refFrame = refFrame;
        end

        function tf = usesControlFrame(~)
            tf = true;
        end

        function cFrame = getControlFrame(obj)
            cFrame = obj.controlFrame;
        end

        function setControlFrame(obj, cFrame)
            obj.controlFrame = cFrame;
        end

        function tf = requiresReInitAfterSoIChange(obj)
            tf = false;
        end

        function enum = getSteeringModelTypeEnum(~)
            enum = SteerModelTypeEnum.LinearTangentAngles;
        end
    end

    methods(Access=private)
        function obj = GenericLinearTangentSteeringModel(gammaAngleModel, betaAngleModel, alphaAngleModel)
            obj.gammaAngleModel = gammaAngleModel;
            obj.betaAngleModel = betaAngleModel;
            obj.alphaAngleModel = alphaAngleModel;
        end
    end

    methods(Static)
        function model = getDefaultSteeringModel()
            gammaAngleModel = PolynominalModel(0,0,0,0);
            betaAngleModel = LinearTangentModel(0,0,0,0,0);
            alphaAngleModel = PolynominalModel(0,0,0,0);

            model = GenericLinearTangentSteeringModel(gammaAngleModel, betaAngleModel, alphaAngleModel);
        end

        function typeStr = getTypeNameStr()
            typeStr = SteeringModelEnum.GenericLinTan.nameStr;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function s=csnr(A,B,row,col)

[n,m,ch]=size(A);

if ch==1
   e=A-B;
   e=e(row+1:n-row,col+1:m-col);
   me=mean(mean(e.^2));
   s=10*log10(255^2/me);
else
   e=A-B;
   e=e(row+1:n-row,col+1:m-col,:);
   e1=e(:,:,1);e2=e(:,:,2);e3=e(:,:,3);
   me1=mean(mean(e1.^2));
   me2=mean(mean(e2.^2));
   me3=mean(mean(e3.^2));
   mse=(me1+me2+me3)/3;
   s  = 10*log10(255^2/mse);
%    s(1)=10*log10(255^2/me1);
%    s(2)=10*log10(255^2/me2);
%    s(3)=10*log10(255^2/me3);
end

return;

author: csjunxu


text:
function mprob = convBaron(prob,opts)
%CONVBARON Convert OPTI problem to BARON Nonlinear problem
%
%   mprob = convBaron(prob,opts)

%   Copyright (C) 2013 Jonathan Currie (I2C2)

%Ensure all args passed
if(nargin < 2)
    error('You must supply both the problem + options');
end
if(~isstruct(prob) || ~isstruct(opts))
    error('Both prob + opts must be structures!');
end

%Get Warning Level
if(strcmpi(opts.warnings,'all'))
    warn = 2;
elseif(strcmpi(opts.warnings,'critical'))
    warn = 1;
else
    warn = 0;
end

%Assign Common Args
mprob.A = prob.A;
mprob.rl = prob.rl;
mprob.ru = prob.ru;
mprob.lb = prob.lb;
mprob.ub = prob.ub;
mprob.xtype = prob.int.str;

%Check if the user wants to use BARON LP/QP/QCQP solver
switch(lower(prob.type))
    case {'lp','milp','bilp','qp','qcqp','miqp','miqcqp'}
        if(warn > 1)
            optiwarn('OPTI:QCQP2NLP','Quadratic Program converted to (MI)NLP to suit BARON interface');
        end
        [prob,opts] = QCQP2NLP(prob,opts);
    case {'uno','nlp','minlp'}
        %nothing
    otherwise
        error('Problem type ''%s'' is not supported by BARON',prob.type);
end

%Assign remaining args
mprob.fun = prob.fun;
mprob.nlcon = prob.nlcon;
mprob.cl = prob.cl;
mprob.cu = prob.cu;
mprob.x0 = prob.x0;

% Set BARON options (optiset options override common ones)
%Build baronset struct, adding any passed options (ignored by baronset if not recognised)
options = baronset(opts.solverOpts);
%Add OPTISET options
options.epsa = opts.tolafun;
options.epsr = opts.tolrfun;
options.inttol = opts.tolint;
options.maxiter = opts.maxnodes;
options.maxtime = opts.maxtime;
options.prlevel = dispLevel(opts.display);
%Save Options Structure
mprob.options = options;

function  print_level = dispLevel(lev)
%Return IPOPT compatible display level
switch(lower(lev))
    case'off'
        print_level = 0;
    case 'iter'
        print_level = 1;
    case 'final'
        print_level = 1;
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function jac = p03_jac ( neqn, t, y )

%*****************************************************************************80
%
%% P03_JAC evaluates the jacobian for problem p03.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the jacobian.
%
%    Output, real JAC(NEQN,NEQN), the jacobian matrix.
%
  jac(1,1) = cos ( t );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ pwIdxsAllrel ] = build_joint_pairs( p )

pidxs = p.pidxs;

if p.person_part
    pidxs = [pidxs pidxs(end)+1];
end

pwIdxsAllrel = cell(1, 1);
n = 0;
for i = 1:length(pidxs)-1
    for j = i+1:length(pidxs)
        n = n + 1;
        pwIdxsAllrel{n} = [i j];
    end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function asa314_test ( )

%*****************************************************************************80
%
%% MAIN is the main program for ASA314_TEST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Roger Payne,
%    Inversion of matrices with contents subject to modulo arithmetic,
%    Applied Statistics,
%    Volume 46, Number 2, 1997, pages 295-298.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA314_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the ASA314 library.\n' );

  asa314_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA314_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% Fugacity coeficient calculation using Soave equation of state
% Author: Arnulfo Rosales-Quintero
% Email: arnol122@gmail.com
% This function calculates either liquid or vapor fugacity coeficients
%-------------------------------------------------------------------------%
function [Phi,Vmix]=Phi(P,T,xy,a,b,phase01)
%-------------------------------------------------------------------------%
    global delta1 delta2 R
    global Kij
%-----------------------------------------------------------------------%
nc = length(xy);
%Mixture rule
amix=0.;	bmix=0.;
for i=1:nc
    dda=0.;
    for j=1:nc
        aij(i,j)= sqrt(a(i)*a(j))*(1.0-Kij(i,j));
        dda     = 2.0 * xy(j) * aij(i,j) + dda;
        amix    = xy(i)*xy(j) * aij(i,j) + amix;	%"a" mixture constant
    end
    dna(i) = dda;						%First derivative   d(an)/dni
    dnb(i) = b(i);						%First derivative   d(bn)/dni
    bmix   = xy(i)*b(i) + bmix;         %"b" mixture constant
end
%-------------------------------------------------------------------------%
%Mixture Volume
    Bb = (P*bmix*(delta1+delta2-1.0) - R*T)/P;
    Cc = (amix - (P*bmix^2+R*T*bmix)*(delta1+delta2) + P*bmix^2*delta1*delta2)/P;
    Dd = (-amix*bmix - (P*bmix^3 + R*T*bmix^2)*delta1*delta2)/P;

switch phase01
    case ('vapor')
        Vmix = max(roots([1 Bb Cc Dd]));
    case ('liquid')
        Vmix = min(roots([1 Bb Cc Dd]));
end
%-------------------------------------------------------------------------%
%Fugacity calculation
for i=1:nc
	Phi(i) = dnb(i)/bmix*(P*Vmix/R/T-1.) - log(P*(Vmix-bmix)/R/T);
    Phi(i) = Phi(i) - amix/bmix*(dna(i)/amix-dnb(i)/bmix)/R/T*log((Vmix+delta1*bmix)/(Vmix+delta2*bmix))/(delta1-delta2);
end
	Phi = exp(Phi);
%-------------------------------------------------------------------------%


text:
function [xyTV,xyLV,xyBV,xyRV] = dividePolygon(xV,yV,xc,yc)
% function [xyTV,xyLV,xyBV,xyRV] = dividePolygon(xV,yV,xc,yc)
%
% This function divides the passed polygon into Top, Left, Bottom and Right
% parts.
%
% APA, 08/23/2012

%xV = planC{indexS.structures}(2).contour(71).segments.points(:,1);
%yV = planC{indexS.structures}(2).contour(71).segments.points(:,2);

if nargin < 3
    xc = mean(xV);
    yc = mean(yV);
end
theta1 = 45*pi/180;
x1 = xc+100*cos(theta1);
y1 = yc+100*sin(theta1);
theta2 = 135*pi/180;
x2 = xc+100*cos(theta2);
y2 = yc+100*sin(theta2);

% Get 45 degree cut
xyTLV = cutpolygon([xV(:) yV(:)], [[xc; x1],[yc;y1]], 'B');
xyBRV = cutpolygon([xV(:) yV(:)], [[xc; x1],[yc;y1]], 'T');

% Top part
xyTV = cutpolygon([xyTLV(:,1) xyTLV(:,2)], [[xc; x2],[yc;y2]], 'B');

% Left part
xyLV = cutpolygon([xyTLV(:,1) xyTLV(:,2)], [[xc; x2],[yc;y2]], 'T');

% Bottom part
xyBV = cutpolygon([xyBRV(:,1) xyBRV(:,2)], [[xc; x2],[yc;y2]], 'T');

% Right part
xyRV = cutpolygon([xyBRV(:,1) xyBRV(:,2)], [[xc; x2],[yc;y2]], 'B');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ temp ] = SAH(temp, vitalslab_hold)

% Matthieu Komorowski - Imperial College London 2017
% will copy a value in the rows below if the missing values are within the
% hold period for this variable (e.g. 48h for weight, 2h for HR...)
% vitalslab_hold = 2x55 cell (with row1 = strings of names ; row 2 = hold time)

h = waitbar(0,'Initializing waitbar...');

hold=table2array(cell2table(vitalslab_hold(2,:)));
nrow=size(temp,1);
ncol=size(temp,2);

lastcharttime=zeros(1,ncol);
lastvalue=zeros(1,ncol);
oldstayid=temp(1,2);

for i=4:ncol
     waitbar(i/ncol,h,i/ncol*100)

    for j=1:nrow

        if oldstayid~=temp(j,2)
            lastcharttime=zeros(1,ncol);
            lastvalue=zeros(1,ncol);
            oldstayid=temp(j,2);
        end

        if isnan(temp(j,i))==0
            lastcharttime(i)=temp(j,3);
            lastvalue(i)=temp(j,i);
        end

        if j>1
        if isnan(temp(j,i)) &&  temp(j,2)==oldstayid && (temp(j,3)-lastcharttime(i))<=hold(i-3)*3600 %note : hold has 53 cols, temp has 55
            temp(j,i)=lastvalue(i);
        end

        end
    end
end

close(h);

end

author: matthieukomorowski


text:
%%
 %  Copyright (c) 2014, Facebook, Inc.
 %  All rights reserved.
 %
 %  This source code is licensed under the BSD-style license found in the
 %  LICENSE file in the root directory of this source tree. An additional grant
 %  of patent rights can be found in the PATENTS file in the same directory.
 %
 %%
function features = prepare_level1_global_features(patches, config)
dims = config.CNN_GLOBAL_MODEL_DIMS;
crop_dims = dims - 2 * config.CNN_GLOBAL_FRINGE;
N = size(patches,4);

% Change the input_dim to be 100 in the deploy.prototxt
matcaffe_init(config.USE_GPU, config.GLOBAL_DEF, config.GLOBAL_MODEL);

d = load(config.GLOBAL_MEAN_FILE);
IMAGE_MEAN = d.image_mean;

% permute from RGB to BGR (IMAGE_MEAN is already BGR)
patches = patches(:, :, [3 2 1], :);
features = zeros(N, 4096 * 2);
input = cell(1, 1);
for i = 1 : ceil(N/ 100)
    fprintf('Extract global features for batch %d\n',i);
    idx = [100 * (i - 1) + 1 : min(N, 100 * i)];
    input{1} = single(zeros(crop_dims(1), crop_dims(2), 3, 100));
    % the top half
    patches_1 = single(patches(1:dims(1),:,:,idx)) - repmat(IMAGE_MEAN, [1 1 1 length(idx)]);
    % flip width and height to make width the fastest dimension
    patches_1 = permute(patches_1, [2 1 3 4]);
    input{1}(:, :, :, 1:length(idx)) = patches_1(config.CNN_GLOBAL_FRINGE + (1:crop_dims(1)), ...
             config.CNN_GLOBAL_FRINGE + (1:crop_dims(2)), :, :);
    results = caffe('forward', input);
    fea_top = reshape(results{1}, [4096 100]);
    clear results;
    fea_top = fea_top(:, 1:length(idx));

    % the bottom half
    patches_2 = single(patches((dims(1)+1): (dims(1)*2),:,:,idx)) - repmat(IMAGE_MEAN, [1 1 1 length(idx)]);
    patches_2 = permute(patches_2, [2 1 3 4]);
    input{1}(:, :, :, 1:length(idx)) = patches_2(config.CNN_GLOBAL_FRINGE + (1:crop_dims(1)),...
        config.CNN_GLOBAL_FRINGE + (1:crop_dims(2)), :, :);
    results = caffe('forward', input);
    fea_bottom = reshape(results{1}, [4096 100]);
    fea_bottom = fea_bottom(:, 1:length(idx));
    features(idx, :)  = [fea_top' fea_bottom'];
end

author: facebookarchive

repo: pose-aligned-deep-networks


text:
function [beta, bo] = svm_multi_predK(X,Y,C,K)
% SVM_MULTI_PREDK
%
% Support Vector Multi Classification
%
% USAGE: [beta, bo] = svm_multi_pred(X,Y,C,K)
%
% PARAMETERS:  X      - (m,d) matrix of m Training inputs in R^d
%
%              Y      - m vector of m Training targets in {1,..,Q}
%
%              C      - Trade-off between regularization and empirical error
%
%              K      - (m,m) Gram matrix
%
%              beta   - (m,Q) Coefficients matrix of the expansion of
%                       the vectors w_i over the x_p's
%                       w_i = \sum_p beta(p,i)x_p
%
%              bo     - (Q,1) bias terms
%
% SUBROUTINES: svm_multi_init.m -> initialization of the optimization problem
%             squadsolve.m          -> optimization
%             compute_kernel.m -> computation of the gram matrix
%
% DESCRIPTION:
%
%              This procedure is the implementation of a multiclass SVM corresponding to
%              the following problem
%
%              \min_{w_i,w_j} \sum_{i\neq j} \|w_i - w_j\|^2 + C\sum_{p,j} \xi_{pj}
%
%              subject to :
%               for all p=1,..,m   ( w_{c(p)} - w_j ).x_p + b_{c(p)} - b_j \geq 1 - \xi_{pj}
%                                  \xi_{pj} \geq 0
%
%      where:  m is the number of training points
%              Q is the number of classes
%              i and j are indices between 1 and Q
%              p is an index between 1 and m
%              w.x is the dot product between w and x
%
%              To solve this problem, a direct approach is used. All the gram matrix K (K_{i,j} =x_i.x_j)
%              is computed before beginning the optimization. The memory cost is thus proportionnal to Q^2*m^2.
%              The solution is not approximated but is provided by a quadratic programming method implemented in
%              squadsolve.m. Instead of minimizing the primal objective function, we maximize the dual. Its
%              formulation is not given here.
%
%
% ERRORS AND BUGS:
%              It is recommended not to use this code for large problems: unless your computer has a huge memory,
%              the procedure will stop with the message: 'memory exhausted'. Use the code msvm_fw.m or chunk_msvm.m
%              instead. Sometimes, the bias may not be computable because the term 'bias_eps' (see below) is too large.
%              In that case, decrease 'bias_eps' and it should work. In all cases, the program will output the bias and
%              try to give the best solution.
%
% EXAMPLE OF USE:
%
%      > [beta,bo] = svm_multi_pred(inputs,targets,10,K);
%
%       -> solve the multiclass optimization problem for a learning set S = \cup_i {(inputs(i,:),targets(i))}
%          with C=10, and with the gram matrix K.
%
% Description of variables:
%
%      kernel = 0,1,2 (integer)    -> type of kernel
%      m              (integer)    -> number of input points
%      Q              (integer)    -> number of classes
%      bias_eps       (real)       -> precision for the computation of the bias
%      reg_eps        (real)       -> to avoid ill conditionning, a diagonal matrix
%                                     is added to the hessian of the minimization problem,
%                                     reg_eps scales this diagonal matrix. (usual < 10^(-6))
%      A              ((Q,Q*m) matrix) -> equality constraint matrix
%      H            ((Q*m,Q*m) matrix) -> hessian of the objective function
%      c               (Q*m vector)    -> linear part of the objectif function
%      alpha           (Q*m vector)    -> variables of the optimization problem
%      constraints     (Q*m vector)    -> store 1-(w_{c(p)} - w_j).x_p
%      beta            ((m,Q) matrix)  -> vector w_i is computed as w_i = sum_p beta_{pi} x_p
%      b               ((Q,Q) matrix)  -> store the differences between bias: b(i,j) = b_i - b_j
% Avoid warning messages
warning off;
disp('Optimizing...');
% Test if nargin is correct
if (nargin < 2) | (nargin > 6),
help svm_multi_pred;
else

% Init
 empty_list_elements_ok = 1;
 m = size(X,1);
 Q = max(Y);
 if (nargin<3),
     C=Inf;
 end;
 if (C==Inf)
    tol = 1e-5;
  else
    tol = C*1e-6;
 end;
 bias_eps=10^(-8)/C;
 reg_eps = 10^(-6);
 sv_eps=10^(-6);
% Message
 % disp(sprintf('\nMulti Support Vector Classification\n'));
 % disp(sprintf('-----------------------------------\n'));

% Set up the parameters for the Optimisation problem
%disp(sprintf('Initialization...\n'));
[H,A,c]=svm_multi_init(Y,K);
% Add small amount of zero order regularization to
% avoid problems when Hessian is badly conditioned.
H = 1/(2*Q^2)*H+reg_eps*eye(size(H));% the problem is to minimize (1/4Q^2)x'Hx - c.x
c = -c;
% Some variables are irrelevant. Consider only relevant variables
mul = Q*(0:m-1);
indx = mul' + Y;
s1 = (1:Q*m);
s2 = (indx);
indx = setdiff(s1,s2)';
clear s1;clear s2;clear mul;
% Solve the optimization problem
% the solution is stored in alpha
%disp(sprintf('Begin the optimization ...\n'));
bo = zeros(size(A,1),1);
[x1,y] = quadsolve(H(indx,indx),c(indx),A(:,indx),bo,C);
x1=x1(1:length(indx));
obj = 0.5*x1'*(H(indx,indx)-reg_eps*eye(length(indx)))*x1 + c(indx)'*x1;
clear H;clear c;
alpha = zeros(Q*m,1);
alpha(indx)=x1;
clear x1;
% Output the objective value
%disp(sprintf('Final objective function: %f \n',-obj));
% if the alpha's are all greater than C then, C can be considered as infinite
% (for the computation of the bias)
if (max(abs(alpha)) < C*0.95)
    C==Inf;
end;
% Compute the number of support vector
% Can be removed if not desired
M = max(alpha);
nsv=0;
for i=1:m,
  if ~(isempty(find(alpha(Q*(i-1)+1:Q*i)>100*eps)))
    nsv = nsv+1;
  end;
end;
%disp(sprintf('Support vectors : %d (%3.1f)\n',nsv,100*nsv/(m)));
% Compute the coefficients beta_ip of the vector w_i = sum_p beta_ip x_p
% clean alpha
beta=zeros(m,Q);
notvoid = find(alpha>=100*eps);
alpha_tmp=zeros(Q*m,1);
alpha_tmp(notvoid) = alpha(notvoid);
alpha=alpha_tmp;
clear notvoid;clear alpha_tmp;
% compute beta
for i =1:Q,
 for p=1:m,
  tmp=0;
  if (i==Y(p))
   tmp = tmp + sum(alpha((p-1)*Q + 1: (p-1)*Q + Q))/Q;
  end;
  beta(p,i) = tmp - alpha((p-1)*Q + i)/Q;
 end;
end;
% Computation of the bias
%disp(sprintf('\n Computation of the bias........\n'));
% Computation of the outputs 1-(w_{c(i)} - w_j).x_i
errorcache = zeros(Q*m,1);
for p=1:m,
  for i=1:Q,
    temp = ((beta(:,Y(p))-beta(:,i))'*K(:,p));
    errorcache((p-1)*Q+i)=temp;
  end;
end;

  % Computation of the output 1-(w_Y(p)-w_i).x_p
  constraints=1-errorcache;
   if C==Inf,
     b=-Inf*ones(Q,Q);
     for i=1:m,
         for k=1:Q,
            b(Y(i),k) = max([b(Y(i),k),(constraints(Q*(i-1)+k) - 0.05)]);
         end;
     end;
     for k=1:Q,
         b(k,k)=0;
         b(k,find(b(k,:)==-Inf))=0;
     end;
   end; %if C>max(alpha),
%%% For finite C
if C<Inf,
% computation of alpha_bias
  % here: computing of the bias that minimize the \sum_{ip} \xi_{pi}
  % subject to linear constraints:
  % (w_{c(p)} - w_i).x_p + b_{c(p)} - b_i >= 1 - \xi_{pi}
  % \xi_{pi}
  % This is a linear program whose dual is easy to compute. It can be
  % solved with the slinearsolve method.
  alpha_bias=zeros(Q*m,1);
  c_un = ones(Q*m,1);
  [alpha_bias(indx)] = slinearsolve(errorcache(indx)-c_un(indx),A(1:(Q-1),indx),zeros(Q-1,1),1);

  % clean alpha_bias
  tmp = 0:m-1;
  tmp = tmp*Q;
  tmp = tmp + Y';
  alpha_bias(tmp)=0;

  % if problem, consider the alpha instead
  if (flag~=1)
    alpha_bias=alpha;
  end;

  b=zeros(Q,Q); % matrix of the difference b_i-b_j
  % Computation of the output 1-(w_Y(p)-w_i).x_p
  constraints=1-errorcache;

  % consider the alpha_ip s.t. 0 < alpha_bias_ip < 1
  % here, i consider alpha_new, because i compute the bias
  % that minimize the l_1 error of the current solution
  % and not of the optimal.
  if C == Inf,
    Cm =1;
  else
    Cm=C;
  end;
  if (flag~=1)
   indice_fix = find(alpha_bias < C*(1-bias_eps) & alpha_bias > Cm*bias_eps);
  else
   indice_fix = find(alpha_bias < (1-bias_eps) & alpha_bias > bias_eps);
  end;

  % if there are no such alpha -> problem
  if isempty(indice_fix)
    disp('WARNING : Unable to compute the bias, use another method\n');
  end;
  % the bias are computed by averaging
  % the matrix 'compteur' counts when  b_i-b_j is computed
  compteur=zeros(Q,Q);
  for ind=indice_fix',
    % find the index Q and p corresponding to ind
    indQ = rem(ind,Q);
    if (indQ==0), indQ=Q;end;
    indp = (ind-indQ)/Q;
    compteur(Y(indp+1),indQ)=compteur(Y(indp+1),indQ)+1;
    b(Y(indp+1),indQ) = b(Y(indp+1),indQ) + constraints(ind);
  end;
  % average by dividing by compteur
  for i=1:Q,
    b(i,i) = 0;
    for j=1:Q,
      if (compteur(i,j))
       b(i,j)=b(i,j)/compteur(i,j);
      end;
    end;
  end;
  end; %if C...,
  clear indice_fix;clear constraints;

  % transform b in order to be symetric
  for i=1:Q,
    for j=i+1:Q,
      if (b(i,j)==0)
        b(i,j)=-b(j,i);
      end;
      if (b(j,i)==0)
        b(j,i)=-b(i,j);
      end;
    end;
  end;

  % average between both symetric elements
  if C==Inf,
   for i=1:Q,
    for j=i+1:Q,
        b(i,j)=max(b(i,j),-b(j,i));
        b(j,i)=-b(i,j);
    end;
  end;
  else,
  for i=1:Q,
    for j=i+1:Q,
        b(i,j)=(b(i,j)-b(j,i))/2;
        b(j,i)=-b(i,j);
    end;
  end;
  end;

  % compute the values of b_i s.t. sum_i b_i = 0
  ok=0;% ok > 0, if computation is possible otherwise ok=0
  b_guess = zeros(Q,1); % if no computation is possible still give a guess
  bo = zeros(Q,1); % the bias vector
  for i=1:Q,
    z = find(b(i,:)~=0);
    if (length(z)==Q-1)
      ok=ok+1;
      temp = sum(b(i,:))/Q;
      for j=1:Q,
        bo(j) =bo(j)+temp-b(i,j);
      end;
    else
      %disp(sprintf('Warning : bias %d -> problem... %d\n',i,length(z)));
      if length(z),
        b_guess(i) = sum(b(i,:))/length(z);
      end;
      for j=1:length(z),
        b_guess(z(j)) =b_guess(i) -b(i,z(j));
      end;
    end;
  end;

  if ~ok
    % then one should compute the bias with another method
    disp('WARNING : Problem with the bias, use another method');
    bo=b_guess;
  else
    bo =bo/ok;
  end;

end; %if nargin...

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [words, k] = feature_color(I, I_type, opts)
% Clusters RGB values of pixels into words of a dictionary using kmeans.

I_lin = reshape(I, size(I,1)*size(I,2), 3);

if opts.load_color_dict
    % sets color_dict
    switch I_type
        case 'rgb'
            load('dicts/rgb_dict_k150');
        case 'nrgb'
            load('dicts/nrgb_dict_k150');
        case 'opp'
            load('dicts/opp_nocut_dict_k150');
            %load('dicts/opp_dict_k150');
        case 'hsv'
            load('dicts/hsv_dict_k150');
    end
    words = vl_ikmeanspush(I_lin', color_dict);
    k = size(color_dict, 2); % number of clusters
else % generate color dictionary for the image
    error('It is recommended to load a color dictionary.');
    k = 50; % number of clusters
    [color_dict, words] = vl_ikmeans(I_lin', k); % integer k-means clustering
    %[color_dict, color_words] = vl_ikmeans(I_lin', k, 'Method', 'elkan'); % integer k-means clustering
    %save('car_color_dict', 'color_dict'); % create example color dictionary
end

words = words(:); % transpose. Using (:) instead of ' to emphasize that the words variable for each feature has the same format

% Visualize color clustering
% color_dict = uint8(color_dict);
% Ic = color_dict(:,words')';
% Ic = reshape(Ic, size(I,1), size(I,2), 3);
% imshow(Ic);

% Shows that there is no indexing error
% words(1:20)
% words = reshape(words, size(I,1), size(I,2));
% image(words) % pixels in correct places, original image distinguishable
% words2 = words(:); % this reverses the above reshape...
% words2(1:20) % ... because this is equal to above words(1:20)

author: Cloud-CV

repo: object-proposals


text:
%DATGAUSS Apply Gaussian filter on images in a dataset
%
% This command is deprecated, use IM_GAUSS instead.

function a = datgauss(varargin)

  a = im_gauss(varargin{:});

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
classdef IMOP6 < PROBLEM
% <multi> <real> <expensive/none>
% Benchmark MOP with irregular Pareto front
% a1 --- 0.05 --- Parameter a1
% a2 ---   10 --- Parameter a2
% K  ---    5 --- Parameter K

%------------------------------- Reference --------------------------------
% Y. Tian, R. Cheng, X. Zhang, M. Li, and Y. Jin, Diversity assessment of
% multi-objective evolutionary algorithms: Performance metric and benchmark
% problems, IEEE Computational Intelligence Magazine, 2019, 14(3): 61-74.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties(Access = private)
        a1 = 0.05;  % Parameter a1
        a2 = 10;    % Parameter a2
        K  = 5;     % Parameter K
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            [obj.a1,obj.a2,obj.K] = obj.ParameterSet(0.05,10,5);
            obj.M = 3;
            if isempty(obj.D); obj.D = 10; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            y1 = mean(PopDec(:,1:2:obj.K),2).^obj.a1;
            y2 = mean(PopDec(:,2:2:obj.K),2).^obj.a2;
            g  = sum((PopDec(:,obj.K+1:end)-0.5).^2,2);
            r  = max(0,min(sin(3*pi*y1).^2,sin(3*pi*y2).^2)-0.05);
            PopObj(:,1) = (1+g).*y1 + ceil(r);
            PopObj(:,2) = (1+g).*y2 + ceil(r);
            PopObj(:,3) = (0.5+g).*(2-y1-y2) + ceil(r);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            [x,y] = meshgrid(linspace(0,1,ceil(sqrt(N))));
            R = [x(:),y(:)];
            r = max(0,min(sin(3*pi*R(:,1)).^2,sin(3*pi*R(:,2)).^2)-0.05);
            R(:,3) = 1 - sum(R,2)/2;
            R = R + repmat(ceil(r),1,3);
            R = R(NDSort(R,1)==1,:);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            [x,y] = meshgrid(linspace(0,1,50));
            z     = 1 - x/2 - y/2;
            R     = [x(:),y(:),z(:)];
            r     = max(0,min(sin(3*pi*R(:,1)).^2,sin(3*pi*R(:,2)).^2)-0.05);
            R     = R + repmat(ceil(r),1,3);
            fes   = NDSort(R,1) == 1;
            z(reshape(~fes,size(z))) = nan;
            R     = {x,y,z};
        end
    end
end


text:
% pop_apcluster() - Clusters STUDY ICs using Affinity Product method.
%                   In this method, IC measures, except equiv. dipoles, (ERP, ERSP...)
%                   are compared for each IC pair and their dissimilarity is multiplied
%                   together to form a combined pairwise dissimilarity matrix. This matrix
%                   is then normalized, weighted and added to the normalized and weighted
%                   IC equiv. dipole distance matrix. The final dissimilarity matrix is
%                   then clustered using affinity clustering  method.
%                   You can control the effect of equiv. dipole distances in
%                   the clustering by setting the 'Relative dipole weight'
%                   parameter in the pop-uo GUI. For example, by setting
%                   this value to 0.8, the final dissimilarity matrix will consist of 80%
%                   distance dissimilarity and 20% of other measures combined together.
%                   Please note that the number of returned clusters may slighty
%                   differ from the number requested in the GUI.
%
%
% Usage:
%     >> STUDY = pop_apcluster(STUDY, ALLEEG) % popup window
%
% See also:  std_apcluster(), std_apreclust(), , pop_apreclust(), apclusterK()
%
% Author: Nima Bigdely-Shamlo, SCCN/INC/UCSD, 2009

function [STUDY ALLEEG command] = pop_apcluster(STUDY, ALLEEG)
% command is used for keeping a history.
% disable measure checkboxes which are not present (calculated) in
% pre-clustering.

% ERP
if isfield(STUDY.preclust,'similarity') && isfield(STUDY.preclust.similarity,'erpCorr')
    erpEnable = 'on';
else
    erpEnable = 'off';
end;
erpChecked = strcmp(erpEnable, 'on'); % only check items if they are enabled.

% ERSP
if isfield(STUDY.preclust,'similarity') && isfield(STUDY.preclust.similarity,'erspCorr')
    erspEnable = 'on';
else
    erspEnable = 'off';
end;
erspChecked = strcmp(erspEnable, 'on'); % only check items if they are enabled.

% ITC
if isfield(STUDY.preclust,'similarity') && isfield(STUDY.preclust.similarity,'itcCorr')
    itcEnable = 'on';
else
    itcEnable = 'off';
end;
itcChecked = strcmp(itcEnable, 'on'); % only check items if they are enabled.

% dipole
if isfield(STUDY.preclust,'similarity') && isfield(STUDY.preclust.similarity,'compDistance')
    dipoleEnable = 'on';
else
    dipoleEnable = 'off';
end;
dipoleChecked = strcmp(dipoleEnable, 'on'); % only check items if they are enabled.

% spectra
if isfield(STUDY.preclust,'similarity') && isfield(STUDY.preclust.similarity,'specCorr')
    specEnable = 'on';
else
    specEnable = 'off';
end;
specChecked = strcmp(specEnable, 'on'); % only check items if they are enabled.

% scalp map
if isfield(STUDY.preclust,'similarity') && isfield(STUDY.preclust.similarity,'mapCorr')
    scalpEnable = 'on';
else
    scalpEnable = 'off';
end;

% popup GUI

returnedFromGui = inputgui( 'geometry', { [1 0.5]  [1 0.5] 1 1 [1 1 1] [1 1 1] [1 1 1] [1 1 1] [1 1 1] [1 1 1] 1 [3 1.5]}, ...
    'geomvert', [], 'uilist', { ...
    { 'style', 'text', 'string', 'Number of clusters to compute:'}, ...
    { 'style', 'edit', 'string', '10' 'tag' 'numberOfClusters' } , ...
       { 'style', 'text', 'string', 'Relative dipole weight (between 0 and 1):'}, ...
    { 'style', 'edit', 'string', '0.8' 'tag' 'numberOfClusters' } , ...
        { 'style', 'text', 'string', [ 'Select measuretures to be used in the clustering:' ] }, {}, ...
    {},{ 'Style', 'checkbox', 'string' 'Dipole' 'tag' 'scale' 'value' 1} , {}, ...
    {},{ 'Style', 'checkbox', 'string' 'ERP' 'tag' 'scale' 'value' erpChecked 'enable' erpEnable}, {},...
    {},{ 'Style', 'checkbox', 'string' 'ERSP' 'tag' 'scale' 'value' erspChecked 'enable' erspEnable}, {},...
    {},{ 'Style', 'checkbox', 'string' 'ITC' 'tag' 'scale' 'value' itcChecked 'enable' itcEnable}, {},...
    {},{ 'Style', 'checkbox', 'string' 'Spectra' 'tag' 'scale' 'value' specChecked 'enable' specEnable}, {}, ...
    {}, { 'Style', 'checkbox', 'string' 'Scalp map' 'tag' 'scale' 'value' 0 'enable' scalpEnable}, {},...
    {}, { 'Style', 'checkbox', 'string' 'Separate outliers (enter std.)' 'tag' 'scale' 'value' 1}, { 'style', 'edit', 'string', '3' 'tag' 'outlierSTD' }, ...

    }, 'helpcom','pophelp(''pop_mpcluster'');', 'title', 'Affinity Product clustering -- pop_apcluster()');

if isempty(returnedFromGui) % an empty returnedFromGui means the Cancel button has been pressed so nothing should be done.
    command = '';
    return; % Cancel button is pressed, so do nothing.
else

    % analysze answers returned from the GUI
    numberOfClusters = str2num(returnedFromGui{1});
    methodParameter = str2num(returnedFromGui{2});

    answers = cell2mat(returnedFromGui(3:end-1));
    measureNamesInGUIorder = {'dipole', 'erp', 'ersp', 'itc', 'spec', 'map'};
    measuresToUseInClustering = measureNamesInGUIorder(find(answers(1:end-1))); %#ok<FNDSB>

    if answers(end) % checkbox for outlier
        outlierSTD = str2num(returnedFromGui{end});
    else
        outlierSTD = Inf;
    end;

    STUDY = std_apcluster(STUDY, ALLEEG, numberOfClusters, outlierSTD, measuresToUseInClustering, methodParameter);

    % prepare 'command' variable for placing both in eeglab histry (accessible with eegh() ) and also
    % adding to  STUDY.history

    measuresInOneString = [];
    for i=1:length(measuresToUseInClustering)
        if i>1
            measuresInOneString = [measuresInOneString ' , ' '''' measuresToUseInClustering{i} ''''];
        else
            measuresInOneString = ['''' measuresToUseInClustering{1} ''''];
        end;
    end;

    % pop up the cluster edit and visualization .
    [STUDY commandFromPop_clustedit] = pop_clustedit(STUDY, ALLEEG);

    command = ['STUDY = std_apcluster(STUDY, ALLEEG, ' num2str(numberOfClusters) ', ' num2str(outlierSTD) ', {' measuresInOneString '} , ' num2str(methodParameter) ');'];
    command = [command '\n' commandFromPop_clustedit];     % add the command from pop_clustedit() to the history too.
    STUDY.history =  sprintf('%s\n%s',  STUDY.history, command);

end;

author: PatternRecognition

repo: OpenBMI


text:
function dist = circle_imp_point_dist_signed_2d ( r, center, p )

%*****************************************************************************80
%
%% CIRCLE_IMP_POINT_DIST_SIGNED_2D: signed distance ( implicit circle, point ) in 2D.
%
%  Discussion:
%
%    The signed distance is zero if the point is on the circle.
%    The signed distance is negative if the point is inside the circle.
%
%    An implicit circle in 2D satisfies the equation:
%
%      ( X - CENTER(1) )**2 + ( Y - CENTER(2) )**2 = R**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real CENTER(2,1), the center of the circle.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, real DIST, the signed distance of the point
%    to the circle.  If the point is inside the circle, the signed distance
%    is negative.
%
  dim_num = 2;

  r2 = sqrt ( sum ( ( p(1:dim_num,1) - center(1:dim_num,1) ).^2 ) );

  dist = r2 - r;

  return
end


text:
function hessfd = getHessianFD(problem, x, d, storedb, key)
% Computes an approx. of the Hessian w/ finite differences of the gradient.
%
% function hessfd = getHessianFD(problem, x, d)
% function hessfd = getHessianFD(problem, x, d, storedb)
% function hessfd = getHessianFD(problem, x, d, storedb, key)
%
% Returns a finite difference approximation of the Hessian at x along d of
% the cost function described in the problem structure. The finite
% difference is based on computations of the gradient.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% If the gradient cannot be computed, an exception is thrown.
%
% See also: approxhessianFD

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   Feb. 19, 2015 (NB):
%       It is sufficient to ensure positive radial linearity to guarantee
%       (together with other assumptions) that this approximation of the
%       Hessian will confer global convergence to the trust-regions method.
%       Formerly, in-code comments referred to the necessity of having
%       complete radial linearity, and that this was harder to achieve.
%       This appears not to be necessary after all, which simplifies the
%       code.
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    if ~canGetGradient(problem)
        up = MException('manopt:getHessianFD:nogradient', ...
            'getHessianFD requires the gradient to be computable.');
        throw(up);
    end

	% Step size
    norm_d = problem.M.norm(x, d);

    % First, check whether the step d is not too small
    if norm_d < eps
        hessfd = problem.M.zerovec(x);
        return;
    end

    % Parameter: how far do we look?
    % (Use approxhessianFD explicitly to gain access to this parameter.)
    epsilon = 1e-4;

    c = epsilon/norm_d;

    % Compute the gradient at the current point.
    grad = getGradient(problem, x, storedb, key);

    % Compute a point a little further along d and the gradient there.
    % Since this is a new point, we need a new key for it, for the storedb.
    x1 = problem.M.retr(x, d, c);
    key1 = storedb.getNewKey();
    grad1 = getGradient(problem, x1, storedb, key1);

    % Transport grad1 back from x1 to x.
    grad1 = problem.M.transp(x1, x, grad1);

    % Return the finite difference of them.
    hessfd = problem.M.lincomb(x, 1/c, grad1, -1/c, grad);

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function mu = truncgaussmf(z,a,b,s)
%TRUNCGAUSSMF Truncated Gaussian membership function.
%   MU = TRUNCGAUSSMF(Z,A,B,S) computes a truncated Gaussian fuzzy
%   membership function. Z is the input variable and can be a vector of
%   any length. A, B, and S are scalar shape parameters. A and B have to
%   be ordered such that A <= B.
%
%       MU = exp(-(Z - B).^2/(2*s^2))   abs(Z - B) <= (B - A)
%       MU = 0                          otherwise
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

mu = zeros(size(z));

c = a + 2*(b - a);
range = (a <= z) & (z <= c);
mu(range) = exp(-(z(range) - b).^2/(2*s^2));

author: dipum

repo: dipum-toolbox

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


text:
function [bcx,bcy] = specific_flow(xbd,ybd)
%leakycavity_flow   Reference problem 5.3 inflow condition
%   [bcx,bcy] = specific_flow(xbd,ybd);
%   input
%          xbd          x coordinate vector
%          ybd          y coordinate vector
%
%   specifies leaky cavity flow boundary condition
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bcx=0*xbd; bcy=0*xbd;
k=find(ybd==1); bcx(k)=1;
return


text:
%KF_LOOP  Performs the prediction and update steps of the Kalman filter
%         for a set of measurements.
%
% Syntax:
%   [MM,PP] = KF_LOOP(X,P,H,R,Y,A,Q)
%
% In:
%   X - Nx1 initial estimate for the state mean
%   P - NxN initial estimate for the state covariance
%   H - DxN measurement matrix
%   R - DxD measurement noise covariance
%   Y - DxM matrix containing all the measurements.
%   A - Transition matrix of the discrete model (optional, default identity)
%   Q - Process noise of the discrete model     (optional, default zero)
%
% Out:
%   MM - Filtered state mean sequence
%   PP - Filtered state covariance sequence
%
%  Description:
%    Calculates state estimates for a set measurements using the
%    Kalman filter. This function is for convience, as it basically consists
%    only of a space reservation for the estimates and of a for-loop which
%    calls the predict and update steps of the KF for each time step in
%    the measurements.
%
%  See also:
%    KF_PREDICT, KF_UPDATE

%  History:
%
%    12.2.2007 JH Initial version.
%
% Copyright (C) 2007 Jouni Hartikainen
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function [MM,PP] = kf_loop(X,P,H,R,Y,A,Q)

  % Check the input parameters.
  if nargin < 5
     error('Too few arguments');
  end
  if nargin < 6
     A = [];
  end
  if nargin < 7
     Q = [];
  end

  % Apply the defaults
  if isempty(A)
    A = eye(size(X,1));
  end
  if isempty(Q)
    Q = zeros(size(X,1));
  end

  % Space for the estimates.
  MM = zeros(size(X,1), size(Y,2));
  PP = zeros(size(X,1), size(X,1), size(Y,2));

  % Filtering steps.
  for i = 1:size(Y,2)
     [X,P] = kf_predict(X,P,A,Q);
     [X,P] = kf_update(X,P,Y(:,i),H,R);
     MM(:,i) = X;
     PP(:,:,i) = P;
  end

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcCos(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

MergeM=get(handles.Matrix_name_edit,'String');
if isfloat(handles.TMatrix(1))
    set(handles.Matrix_name_edit,'String',['cos([' MergeM '])']);
else
    set(handles.Matrix_name_edit,'String',['cos(double([' MergeM ']))']);
end
MU_calc_matrix(handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab