text:
function [dat] = volplot(x, y, z, dat, sel, cscale)

% VOLPLOT make 2D or 3D plot of volumetric data (e.g. MRI)
% that is defined on a regular orthogonal grid
%
% volplot(dat, sel) or
% volplot(x, y, z, dat, sel)
% volplot(x, y, z, dat, sel, caxis)
%
% where sel is one of
%   [x, y, z]     intersection through the three orthogonal directions
%   index         linear index of the voxel of interest
%   'min'         intersection at the minimum
%   'max'         intersection at the maximum
%   'center'      intersect at the center of each axis
%   'interactive' intersect at the center, then go into interactive mode
%   'maxproject'  project the maximum value along each orthogonal direction
%   'sumproject'  integrated value along each orthogonal direction (glassbrain)
%   'montage'     show all slices
% and caxis is the [min max] used for the color scaling
%
% See also TRIPLOT, LINEPLOT (in ~roberto/matlab/misc)
% See also NDGRID

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

if nargin<2
  dat = x;
  sel = 'interactive';
  x = 1:size(dat,1);
  y = 1:size(dat,2);
  z = 1:size(dat,3);
elseif nargin<3
  dat = x;
  sel = y;
  x = 1:size(dat,1);
  y = 1:size(dat,2);
  z = 1:size(dat,3);
elseif nargin>4
  if length(size(dat))==3
    if length(x)~=size(dat,1), ft_error('incorrect x-axis specification'), end
    if length(y)~=size(dat,2), ft_error('incorrect y-axis specification'), end
    if length(z)~=size(dat,3), ft_error('incorrect z-axis specification'), end
  end
else
  ft_error('incorrect number of input arguments');
end

if nargin==6
  cmin = cscale(1);
  cmax = cscale(2);
else
  % ensure same color scaling for all figures
  cmin = min(dat(:));
  cmax = max(dat(:));
end

% convert from 1-d to 3-d array
if any(size(dat)==1)
  dim(1) = length(x);
  dim(2) = length(y);
  dim(3) = length(z);
  dat = reshape(dat, dim);
else
  dim = size(dat);
end

% convert the selection to the indices of the x/y/z intersection
if ischar(sel) && strcmp(sel, 'min')
  [minval, minindx] = min(dat(:));
  [xi, yi, zi] = ind2sub(size(dat), minindx);
elseif ischar(sel) && strcmp(sel, 'max')
  [maxval, maxindx] = max(dat(:));
  [xi, yi, zi] = ind2sub(size(dat), maxindx);
elseif ischar(sel) && strcmp(sel, 'center')
  xi = round(length(x)/2);
  yi = round(length(y)/2);
  zi = round(length(z)/2);
elseif ischar(sel) && strcmp(sel, 'interactive')
  xi = round(length(x)/2);
  yi = round(length(y)/2);
  zi = round(length(z)/2);
elseif ~ischar(sel) && length(sel)==1
  [xi, yi, zi] = ind2sub(dim, sel);
else
  xi = nearest(x, sel(1));
  yi = nearest(y, sel(2));
  zi = nearest(z, sel(3));
end

% start the plotting
if strcmp(sel, 'interactive')
  xc = x(xi);
  yc = y(yi);
  zc = z(zi);
  nas = [];
  lpa = [];
  rpa = [];
  while(1)
    fprintf('============================================================\n');
    fprintf('click with mouse button to reslice the display to a new position\n');
    fprintf('press n/l/r on keyboard to record a fiducial position\n');
    fprintf('press q on keyboard to quit interactive mode\n');
    volplot(x, y, z, dat, [xc yc zc]);
    drawnow;
    try, [d1, d2, key] = ginput(1); catch, key='q'; end
    if key=='q'
      % rename the first output argument of this function
      h = nas;
      break;
    elseif key=='l'
      lpa = [xc yc zc];
    elseif key=='r'
      rpa = [xc yc zc];
    elseif key=='n'
      nas = [xc yc zc];
    elseif key=='x'
      selx= round((xc-0):(xc+0));
      sely= round((yc-0):(yc+0));
      selz= round((zc-0):(zc+0));
      dat(selx,sely,selz) = 0;
    elseif key=='X'
      selx= round((xc-1):(xc+1));
      sely= round((yc-1):(yc+1));
      selz= round((zc-1):(zc+1));
      dat(selx,sely,selz) = 0;
    else
      % update the view to a new position
      l1 = get(get(gca, 'xlabel'), 'string');
      l2 = get(get(gca, 'ylabel'), 'string');
      switch l1
        case 'x'
          xc = d1;
        case 'y'
          yc = d1;
        case 'z'
          zc = d1;
      end
      switch l2
        case 'x'
          xc = d2;
        case 'y'
          yc = d2;
        case 'z'
          zc = d2;
      end
    end
    if ~isempty(nas), fprintf('nas = [%f %f %f]\n', nas); else fprintf('nas = undefined\n'); end
    if ~isempty(lpa), fprintf('lpa = [%f %f %f]\n', lpa); else fprintf('lpa = undefined\n'); end
    if ~isempty(rpa), fprintf('rpa = [%f %f %f]\n', rpa); else fprintf('rpa = undefined\n'); end
  end

elseif strcmp(sel, 'montage')
  % make plot of x-y slices for all z values
  maxval = max(dat(:));
%  for z=1:size(dat,3)
%    % convert to 4D image for montage display
%    % transpose to correct for x-y axis change in MATLAB image function
%    img(:,:,1,z) = transpose(dat(:,:,z));
%  end
%  montage(img);
%  axis xy

  % this works if java fails
  siz    = size(dat);
  ndiv   = [ceil(sqrt(siz(3))) floor(sqrt(siz(3)))];
  map    = zeros(siz(2)*ndiv(2),siz(1)*ndiv(1));
  for k=1:siz(3)
    [nx,ny] = ind2sub(ndiv,k);
    map(siz(2)*(ny-1)+1:siz(2)*ny,siz(1)*(nx-1)+1:siz(1)*nx) = dat(:,:,k)';
  end
  imagesc(map);axis xy;axis equal;axis off;
  caxis([cmin cmax]);
  colormap jet

elseif strcmp(sel, 'sumproject')
  % make plot of integrated-value projection along the thee orthogonal directions
  delete(subplot(2,2,4));   % delete the old colorbar
  h1 = subplot(2,2,1);
  h2 = subplot(2,2,2);
  h3 = subplot(2,2,3);
  h4 = subplot(2,2,4);  % this will  be the new colorbar

  % change not-a-number values to zero
  dat(find(isnan(dat(:)))) = 0;

  % update cmin and cmax
  cmin = 0;
  cmax = squeeze(max(max(sum(dat,1))));
  cmax = max(cmax, squeeze(max(max(sum(dat,2)))));
  cmax = max(cmax, squeeze(max(max(sum(dat,3)))));

  subplot(h1);
  imagesc(x, z, squeeze(sum(dat, 2))'); set(gca, 'ydir', 'normal')
  axis equal; axis tight;
  xlabel('x'); ylabel('z');
  caxis([cmin cmax]);

  subplot(h2);
  imagesc(y, z, squeeze(sum(dat, 1))'); set(gca, 'ydir', 'normal')
  axis equal; axis tight;
  xlabel('y'); ylabel('z');
  caxis([cmin cmax]);

  subplot(h3);
  imagesc(x, y, squeeze(sum(dat, 3))'); set(gca, 'ydir', 'normal')
  axis equal; axis tight;
  xlabel('x'); ylabel('y');
  caxis([cmin cmax]);

  subplot(h4);
  imagesc(cmin:cmax);
  caxis([cmin cmax]);
  xlabel('colorscale')

elseif strcmp(sel, 'maxproject')
  % make plot of maximum-value projection along the thee orthogonal directions
  delete(subplot(2,2,4));   % delete the old colorbar
  h1 = subplot(2,2,1);
  h2 = subplot(2,2,2);
  h3 = subplot(2,2,3);
  h4 = subplot(2,2,4);  % this will  be the new colorbar

  subplot(h1);
  imagesc(x, z, squeeze(max(dat, [], 2))'); set(gca, 'ydir', 'normal')
  axis equal; axis tight;
  xlabel('x'); ylabel('z');
  caxis([cmin cmax]);

  subplot(h2);
  imagesc(y, z, squeeze(max(dat, [], 1))'); set(gca, 'ydir', 'normal')
  axis equal; axis tight;
  xlabel('y'); ylabel('z');
  caxis([cmin cmax]);

  subplot(h3);
  imagesc(x, y, squeeze(max(dat, [], 3))'); set(gca, 'ydir', 'normal')
  axis equal; axis tight;
  xlabel('x'); ylabel('y');
  caxis([cmin cmax]);

  subplot(h4);
  colorbar(h4, 'peer', h1);
  xlabel('colorscale')

else
  % make plot of three orthogonal slices intersecting at [xi yi zi]
  if ~exist('xi', 'var') || ~exist('yi', 'var') || ~exist('zi', 'var')
    ft_error('nothing to plot, no selection given')
  end

  fprintf('value of %f in voxel %d at [%.02f %.02f %.02f]\n', double(dat(xi, yi, zi)), sub2ind(dim, xi, yi, zi), x(xi), y(yi), z(zi));

  warning off
  delete(subplot(2,2,4));   % delete the old colorbar
  warning on
  h1 = subplot(2,2,1);
  h2 = subplot(2,2,2);
  h3 = subplot(2,2,3);
  h4 = subplot(2,2,4);  % this will  be the new colorbar

  subplot(h1);
  imagesc(x, z, squeeze(dat(:,yi,:))'); set(gca, 'ydir', 'normal')
  axis equal; axis tight;
  xlabel('x'); ylabel('z');
  caxis([cmin cmax]);
  ft_plot_crosshair([x(xi) z(zi)], 'color', 'yellow');

  subplot(h2);
  imagesc(y, z, squeeze(dat(xi,:,:))'); set(gca, 'ydir', 'normal')
  axis equal; axis tight;
  xlabel('y'); ylabel('z');
  caxis([cmin cmax]);
  ft_plot_crosshair([y(yi) z(zi)], 'color', 'yellow');

  subplot(h3);
  imagesc(x, y, squeeze(dat(:,:,zi))'); set(gca, 'ydir', 'normal')
  axis equal; axis tight;
  xlabel('x'); ylabel('y');
  caxis([cmin cmax]);
  ft_plot_crosshair([x(xi) y(yi)], 'color', 'yellow');

  subplot(h4);
  imagesc(cmin:((cmax-cmin)./64):cmax);
  caxis([cmin cmax]);
  xlabel('colorscale')
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Properties
%
%%

author: mtex-toolbox

repo: mtex

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


text:
function xyt = Xt(s)
x = s ;
y = 1-3*s+3*s^2 ;

xyt = [x ; y] ;

author: Sable

repo: mcbench-benchmarks


text:
function [ value, ifault ] = ppnd ( p )

%*****************************************************************************80
%
%% PPND produces the normal deviate value corresponding to lower tail area = P.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by J Beasley, S Springer.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    J Beasley, S Springer,
%    Algorithm AS 111:
%    The Percentage Points of the Normal Distribution,
%    Applied Statistics,
%    Volume 26, Number 1, 1977, pages 118-121.
%
%  Parameters:
%
%    Input, real P, the value of the cumulative probability
%    densitity function.  0 < P < 1.
%
%    Output, real VALUE, the normal deviate value with the property that
%    the probability of a standard normal deviate being less than or
%    equal to PPND is P.
%
%    Output, integer IFAULT, error flag.
%    0, no error.
%    1, P <= 0 or P >= 1.  PPND is returned as 0.
%
  a0 = 2.50662823884;
  a1 = -18.61500062529;
  a2 = 41.39119773534;
  a3 = -25.44106049637;
  b1 = -8.47351093090;
  b2 = 23.08336743743;
  b3 = -21.06224101826;
  b4 = 3.13082909833;
  c0 = -2.78718931138;
  c1 = -2.29796479134;
  c2 = 4.85014127135;
  c3 = 2.32121276858;
  d1 = 3.54388924762;
  d2 = 1.63706781897;
  split = 0.42;

  ifault = 0;
%
%  0.08 < P < 0.92
%
  if ( abs ( p - 0.5 ) <= split )

    r = ( p - 0.5 ) * ( p - 0.5 );

    value = ( p - 0.5 ) * ( ( ( ...
        a3   * r ...
      + a2 ) * r ...
      + a1 ) * r ...
      + a0 ) / ( ( ( ( ...
        b4   * r ...
      + b3 ) * r ...
      + b2 ) * r ...
      + b1 ) * r ...
      + 1.0 );
%
%  P < 0.08 or P > 0.92,
%  R = min ( P, 1-P )
%
  elseif ( 0.0 < p && p < 1.0 )

    if ( 0.5 < p )
      r = sqrt ( - log ( 1.0 - p ) );
    else
      r = sqrt ( - log ( p ) );
    end

    value = ( ( ( ...
        c3   * r ...
      + c2 ) * r ...
      + c1 ) * r ...
      + c0 ) / ( ( ...
        d2   * r ...
      + d1 ) * r ...
      + 1.0 );

    if ( p < 0.5 )
      value = - value;
    end
%
%  P <= 0.0 or 1.0 <= P
%
  else

    ifault = 1;
    value = 0.0;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function alpha_filter = getAlphaFilter(kgrid, medium, filter_cutoff, taper_ratio)
%GETALPHAFILTER Create filter for medium.alpha_filter.
%
% DESCRIPTION:
%       getAlphaFilter uses getWin to create a Tukey window via rotation to
%       pass to the medium.alpha_filter input field of the first order
%       simulation functions (kspaceFirstOrder1D, kspaceFirstOrder2D, and
%       kspaceFirstOrder3D). This parameter is used to regularise time
%       reversal image reconstruction when absorption compensation is
%       included.
%
% USAGE:
%       alpha_filter = getAlphaFilter(kgrid, medium, filter_cutoff)
%       alpha_filter = getAlphaFilter(kgrid, medium, filter_cutoff, taper_ratio)
%
% INPUTS:
%       kgrid           - k-Wave grid structure returned by makeGrid
%       medium          - k-Wave medium structure
%       filter_cutoff   - alpha_filter cutoff frequency [Hz], where
%                         filter_cutoff = [f_all] in 1D
%                         filter_cutoff = [f_all] or [f_x, f_y] in 2D
%                         filter_cutoff = [f_all] or [f_x, f_y, f_z] in 3D
%
%                         Any of the filter_cutoff inputs may be set to
%                         'max' to set the cutoff frequency to the maximum
%                         frequency supported by the grid
%
% OPTIONAL INPUTS
%       taper_ratio     - taper ratio for Tukey Window (default = 0.5)
%
% OUTPUTS:
%       alpha_filter    - filter
%
% ABOUT:
%       author          - Bradley Treeby
%       date            - 6th May 2010
%       last update     - 28th October 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also getWin

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% update the command line status
disp('Creating absorption filter...');

% check to see if the taper ratio is given
if nargin == 3
    taper_ratio = 0.5;
elseif nargin ~= 4
    error('incorrect number of inputs');
end

% update command line status
disp(['  taper ratio: ' num2str(taper_ratio)]);

% extract the dimensions
dim = numDim(kgrid.k);

% extract the alpha_filter cutoff frequencies
if numel(filter_cutoff) == 1
    filter_cutoff_x = filter_cutoff;
    if dim > 2
        filter_cutoff_z = filter_cutoff;
    end
    if dim > 1
        filter_cutoff_y = filter_cutoff;
    end
elseif numel(filter_cutoff) ~= dim
    error(['input filter_cutoff must have 1 or ' num2str(dim) ' elements for a ' num2str(dim) 'D grid']);
else
    if dim == 1
        filter_cutoff_x = filter_cutoff{1};
    end
    if dim == 2
        filter_cutoff_x = filter_cutoff{1};
        filter_cutoff_y = filter_cutoff{2};
    end
    if dim == 3
        filter_cutoff_x = filter_cutoff{1};
        filter_cutoff_y = filter_cutoff{2};
        filter_cutoff_z = filter_cutoff{3};
    end
end

% extract the maximium sound speed
c = max(medium.sound_speed(:));

% calculate the alpha_filter size in the z direction for 3D data
if dim > 2
    if strcmp(filter_cutoff_z, 'max')
        % set the alpha_filter size to be the same as the grid size
        filter_size_z = kgrid.Nz;
        filter_cutoff_z = kgrid.kz_max*c/(2*pi);
    else
        % convert the cutoff frequency to a wavenumber
        k_cutoff_z = 2*pi*filter_cutoff_z ./ c;

        % set the alpha_filter size
        filter_size_z = round(kgrid.Nz*k_cutoff_z/kgrid.kz(end));

        % check the alpha_filter size
        if filter_size_z > kgrid.Nz
            % set the alpha_filter size to be the same as the grid size
            filter_size_z = kgrid.Nz;
            filter_cutoff_z = kgrid.kz_max*c/(2*pi);
        end
    end
end

% calculate the alpha_filter size in the y direction for 2 and 3D data
if dim > 1
    if strcmp(filter_cutoff_y, 'max')
        % set the alpha_filter size to be the same as the grid size
        filter_size_y = kgrid.Ny;
        filter_cutoff_y = kgrid.ky_max*c/(2*pi);
    else
        % convert the cutoff frequency to a wavenumber
        k_cutoff_y = 2*pi*filter_cutoff_y ./ c;

        % set the alpha_filter size
        filter_size_y = round(kgrid.Ny*k_cutoff_y/kgrid.ky(end));

        % check the alpha_filter size
        if filter_size_y > kgrid.Ny
            % set the alpha_filter size to be the same as the grid size
            filter_size_y = kgrid.Ny;
            filter_cutoff_y = kgrid.ky_max*c/(2*pi);
        end
    end
end

% calculate the alpha_filter size in the x direction for 1, 2, and 3D data
if strcmp(filter_cutoff_x, 'max')
    % set the alpha_filter size to be the same as the grid size
    filter_size_x = kgrid.Nx;
    filter_cutoff_x = kgrid.kx_max*c/(2*pi);
else
    % convert the cutoff frequency to a wavenumber
    k_cutoff_x = 2*pi*filter_cutoff_x ./ c;

    % set the alpha_filter size
    filter_size_x = round(kgrid.Nx*k_cutoff_x/kgrid.kx(end));

    % check the alpha_filter size
    if filter_size_x > kgrid.Nx
        % set the alpha_filter size to be the same as the grid size
        filter_size_x = kgrid.Nx;
        filter_cutoff_x = kgrid.kx_max*c/(2*pi);
    end
end

% create the alpha_filter using getWin
switch dim
    case 1
        % create the alpha_filter
        filter_sec = getWin(filter_size_x, 'Tukey', 'Param', taper_ratio, 'Rotation', true);

        % enlarge the alpha_filter to the size of the grid
        alpha_filter  = zeros(kgrid.Nx, 1);
        x_index = round((kgrid.Nx - filter_size_x)/2) + 1;
        alpha_filter(x_index:x_index + filter_size_x - 1) = filter_sec;

        % update the command line status
        disp(['  filter cutoff: ' scaleSI(filter_cutoff_x) 'Hz']);
    case 2
        % create the alpha_filter
        filter_sec = getWin([filter_size_x, filter_size_y], 'Tukey', 'Param', taper_ratio, 'Rotation', true);

        % enlarge the alpha_filter to the size of the grid
        alpha_filter  = zeros(kgrid.Nx, kgrid.Ny);
        x_index = round((kgrid.Nx - filter_size_x)/2) + 1;
        y_index = round((kgrid.Ny - filter_size_y)/2) + 1;
        alpha_filter(x_index:x_index + filter_size_x - 1, y_index:y_index + filter_size_y - 1) = filter_sec;

        % update the command line status
        disp(['  filter cutoff: ' scaleSI(filter_cutoff_x) 'Hz by ' scaleSI(filter_cutoff_y) 'Hz']);
    case 3
        % create the alpha_filter
        filter_sec = getWin([filter_size_x, filter_size_y, filter_size_z], 'Tukey', 'Param', taper_ratio, 'Rotation', true);

        % enlarge the alpha_filter to the size of the grid
        alpha_filter  = zeros(kgrid.Nx, kgrid.Ny, kgrid.Nz);
        x_index = round((kgrid.Nx - filter_size_x)/2) + 1;
        y_index = round((kgrid.Ny - filter_size_y)/2) + 1;
        z_index = round((kgrid.Nz - filter_size_z)/2) + 1;
        alpha_filter(x_index:x_index + filter_size_x - 1, y_index:y_index + filter_size_y - 1, z_index:z_index + filter_size_z - 1) = filter_sec;

        % update the command line status
        disp(['  filter cutoff: ' scaleSI(filter_cutoff_x) 'Hz by ' scaleSI(filter_cutoff_y) 'Hz by ' scaleSI(filter_cutoff_z) 'Hz']);
end

repo: Aero-matlab


text:
%{
Authors:
Jonatas Lopes de Paiva
Claudio Fabiano Motta Toledo
Helio Pedrini

%}

function g = mutation(img)

r = randi(3);

switch r

    case 1
        tmp = double(img);
        rate = 0.7 + 0.6*rand;
        i = rate*tmp;
        i = floor(i + 0.5);
    case 2
        t = getFilterSize(2);
        sigma = rand * 5 + 0.05;
        h = fspecial('gaussian', t, sigma);
    case 3
        t = getFilterSize(2);
        h = fspecial('average', t);

end

if (r ==2 || r == 3)
    i = imfilter(img, h);
end

g = uint8(i);

end

function s = getFilterSize(maxSize)
n = randi(maxSize);

s = 2*n+1;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [maskedSlices, ROIRows, ROICols, ROIVoxels] = ROIMaskedSlices(inputSlices3M,maskM)
%  This function loads a 3D image array for masking with a single mask.
%  This could be an array of DCE images of one anatomic slice over time
%  or an array of varied flip angle images at one slice for T1 calculation
%  The calling program supplies the slice mask from the 3D mask array.
%  This function generates a masked version of the 3D array and
%  it also returns the coordinates of each voxel in the ROI
%  and the total number of voxels in the ROI
%
%  Kristen Zakian
%

% Get coordinates of in-ROI voxels
[ROIRows , ROICols] = find(maskM);
% Get no. voxels in ROI
ROIVoxels = nnz(maskM);

% Mask input array
%maskedSlices = bsxfun(@times, int16(maskM),int16(inputSlices));
maskedSlices = bsxfun(@times,maskM,double(inputSlices3M)); %CHANGED AI 2/7/17
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function xform = mrBaseXform(base, map);
%
% xform = mrBaseXform(base, map);
%
% Compute a 4x4 transform matrix to translate from
% data coords in a map volume to data coords in a
% base volume.
%
% base and map are both MR objects (see mrLoad). Will
% prompt if either is not specified.
%
%
% ras 07/05.
if ~exist('base', 'var') | isempty(base),  base = mrLoad; end
if ~exist('map', 'var') | isempty(map),  map = mrLoad; end

% first,  check if the two MR volumes cover the same space
% at the same resolution: then they're already in the same
% coordinate space and the xform is just an identity matrix:
if isequal(base.dims(1:3), map.dims(1:3)) & ...
        isequal(base.voxelSize(1:3), map.voxelSize(1:3))
    disp('Base and map volumes are coregistered')
    xform = eye(4);
    return
end

% next,  check that they don't specify the same extent
% in 3D space,  even if they may be at different resolutions:
if isequal(base.extent(1:3), map.extent(1:3)) | ...
        abs(base.extent(1:3)-map.extent(1:3)) < [1 1 1]
    disp('Base and map cover same extent -- scaling map to fit base')
    szRatio = map.voxelSize(1:3) ./ base.voxelSize(1:3);
    xform = affineBuild([0 0 0], [0 0 0], szRatio);
    return
end

% otherwise,  check if they share a common coordinate space.
% First,  remove the standard spaces (raw data in pix/mm,  L/R
% flipped -- see mrStandardSpaces):
spacesA = {base.spaces.name};
spacesB = {map.spaces.name};
std = {'Raw Data in Pixels' 'Raw Data in mm' 'L/R Flipped'};
spacesA = setdiff(spacesA, std);
spacesB = setdiff(spacesB, std);
common = intersect(spacesA, spacesB);

if isempty(common)
    error('base and map do not share any common coordinate spaces!');
end

% if we're here,  this means there is a common coord system.
% use the first common space (it shouldn''t matter which one
% is used) to derive the xform:
common = common{1};
fprintf('Coregistering via common space %s \n', common);
spacesA = {base.spaces.name};
spacesB = {map.spaces.name};
iA = cellfind(spacesA, common);
iB = cellfind(spacesB, common);
xform = inv(base.spaces(iA).xform) * map.spaces(iB).xform;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
filename='Cantilever_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'IPOPT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [aDeriv,aJacob,aHess,papt]=aCVRuv(x)
%ACVRUV The drift function for a continuous-time motion model where the
%       the target state is given in monostatic range and u-v direction
%       cosine coordinates and the motion is constant-velocity in 3D
%       Cartesian coordinates.
%
%INPUTS: x The 6XN state vector of N targets in the order of
%          [r;u;;v;rDot;uDot;vDot], where r is monostatic range and u and v
%          are direction cosines. The target is assumed to always be in
%          front of the radar and not pass behind the radar.
%
%OUTPUTS: aDeriv The 6XN set of time-derivatives of the N state vectors
%                under the Cartesian linear motion model in r-u-v
%                coordinates.
%         aJacob The 6X6XN set of Jacobians of aDeriv, which can be useful
%                in extended Kalman filters. This is the derivative of each
%                component of aDeriv (selected by row) with respect to the
%                elements of the state [r,u,v,rDot,uDot,vDot] selected by
%                column.
%
%A full derivation of the dynamic model is provided in [1]. Summarizing it
%here, let rVec be a position vector. We shall use the orthonormal basis
%vectors:
%uVec1=[u;v;sqrtw]
%uVec2=[sqrt(w/(1-v2));0;-u/sqrt(1-v2)]
%uVec3=[-u*v/sqrt(1-v2);sqrt(1-v2);-v*sqrt(w/(1-v2))];
%Note that uVec2 is a normalized version of the derivative of uVec1 with
%respect to u and uVec3 is cross(u1,u2). The time-derivatives of the basis
%vectors can themselves be expressed as basis vectors:
%uVec1Dot=c1*u2+c2*u3
%uVec2Dot=c3*u1+c4*c3
%uVec3Dot=c5*u1+c6*u2
%The coefficients c1 through c6 are omitted for brevity but can be found
%by taking the total derivative of each basis vector with respect to time
%and then taking the dot products of the result with the orthonormal basis
%vectors.
%The vector rVec is written
%rVec=r*u1
%The velocity vector is the derivative of rVec. Using the chain rule and
%substituting in the derivative parameters for the basis vectors,
%rVecDot=rDot*uVec1+r*c1*uVec2+r*c2*uVec3
%The decond derivative is more complicated, but simplifies to
%rVecDDot=a1*uVec1+a2*uVec2+a3*uVec3
%where two D's indicates two time-derivatives.
%a1 is a radial acceleration and a3 and a3 are two other Cartesian
%components of acceleration. The expressions for them are
% a1=rDDot+r*(c1*c3+c2*c5);
% a2=2*rDot*c1+r*(c1Dot+c2*c6);
% a3=2*rDot*c2+r*(c2Dot+c1*c4);
%The expression for the drift function comes from setting a1=a2=a3=0 and
%solving for the second derivative terms. The second derivatives in u and v
%arise through the derivatives of c1 and c2.
%
%EXAMPLE:
%Here, we verify that integrating forward with this model is equivalent to
%linear motion in Cartesian coordinates.
% xInitCart=[400;920;1000;-100;20;64];
% T=50;%The prediction time.
% F=FPolyKal(T,6,1);
% xInitRuv=stateCart2Ruv(xInitCart);
% xEndCart=F*xInitCart;
% RelTol=1e-10;
% AbsTol=1e-13;
% xStepsRuv=RKAdaptiveOverRange(xInitRuv,[0;T],@(x,t)aCVRuv(x),0.01,0,[],[],RelTol,AbsTol);
% xEndRuvRK=xStepsRuv(:,end);
% xEndRuvExact=stateCart2Ruv(xEndCart);
% max(abs(xEndRuvRK-xEndRuvExact)./xEndRuvExact)
%One will observe that the error is less than 1e-11, which is a good
%agreement.
%
%REFERENCES:
%[1] D. F. Crouse, "Basic Linear Dynamic Models in Local Coordinates,"
%    Naval Research Laboratory, Washington, D.C., Tech. Rep.
%    NRL/MR/5344--19-9882, 24 Aug. 2019.
%
%August 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

r=x(1,:);
u=x(2,:);
v=x(3,:);
rDot=x(4,:);
uDot=x(5,:);
vDot=x(6,:);

u2=u.^2;
v2=v.^2;
r2=r.^2;

w2=max(0,1-u2-v2);
w=sqrt(w2);

diffV2=1-v2;
diffV=sqrt(diffV2);

diffU2=1-u2;
diffU=sqrt(diffU2);

uDot2=uDot.*uDot;
vDot2=vDot.*vDot;
w3=w2.*w;

rDDot=-(r*(-uDot2*diffV2-2*u*uDot*v*vDot-diffU2*vDot2))/w2;
uDDot=-((2*rDot*uDot)/r)-u*(vDot2*diffU2+uDot2*diffV2+2*u*v*uDot*vDot)/w2;
vDDot=-((2*rDot*vDot)/r)-v*(uDot2*diffV2+vDot2*diffU2+2*u*v*uDot*vDot)/w2;

aDeriv=[rDot;uDot;vDot;rDDot;uDDot;vDDot];

if(nargout>1)
    N=size(x,2);
    if(N>1)
        error('Derivatives are only available for numPoints=1.')
    end

    w4=w2^2;

    ddrDDdr=-(-uDot2*diffV2-2*u*uDot*v*vDot-diffU2*vDot2)/w2;
    ddrDDdu=(2*r*(u*uDot+v*vDot)*(uDot-uDot*v2+u*v*vDot))/w4;
    ddrDDdv=-((2*r*(-u*uDot*v-diffU2*vDot)*(u*uDot+v*vDot))/w4);

    ddrDDduDot=((2*r*(uDot-uDot*v2+u*v*vDot))/w2);
    ddrDDdvDot=-(2*r*(-u*uDot*v-diffU2*vDot))/w2;

    dduDDdr=(2*rDot*uDot)/r2;
    dduDDdu=(-uDot2*(1+u2-v2)*diffV2-4*u*uDot*v*diffV2*vDot+(-1+v2-u2*(-2+u2+3*v2))*vDot2)/w4;
    dduDDdv=(2*u*(-u*uDot*v-diffU2*vDot)*(u*uDot+v*vDot))/w4;
    dduDDdrDot=-((2*uDot)/r);
    dduDDduDot=-((2*rDot)/r)-(2*u*(uDot-uDot*v2+u*v*vDot))/w2;
    dduDDdvDot=-(2*u*(u*uDot*v+vDot-u2*vDot))/w2;

    ddvDDdr=(2*rDot*vDot)/r2;
    ddvDDdu=-((2*v*(u*uDot+v*vDot)*(uDot-uDot*v2+u*v*vDot))/w4);
    ddvDDdv=(uDot2*(u2*(1-3*v2)-diffV2^2)-4*u*diffU2*uDot*v*vDot+diffU2*(-1+u2-v2)*vDot2)/w4;
    ddvDDdrDot=-((2*vDot)/r);
    ddvDDduDot=-(2*v*(uDot-uDot*v2+u*v*vDot))/w2;
    ddvDDdvDot=-((2*rDot)/r)-(2*v*(u*uDot*v+vDot-u2*vDot))/w2;

    aJacob=[0,               0,              0,              1,              0,              0;
            0,               0,              0,              0,              1,              0;
            0,               0,              0,              0,              0,              1;
            ddrDDdr,         ddrDDdu,        ddrDDdv,        0,              ddrDDduDot,     ddrDDdvDot;
            dduDDdr,         dduDDdu,        dduDDdv,        dduDDdrDot,     dduDDduDot,     dduDDdvDot;
            ddvDDdr,         ddvDDdu,        ddvDDdv,        ddvDDdrDot,     ddvDDduDot,     ddvDDdvDot];

    if(nargout>2)
        aHess=zeros(6,6,6);

        u4=u2*u2;
        v4=v2*v2;
        w6=w3*w3;
        diffU4=diffU2*diffU2;
        diffV4=diffV2*diffV2;

        ddrDDdrdr=0;
        ddrDDdudr=(2*(u*uDot+v*vDot)*(uDot-uDot*v2+u*v*vDot))/w4;
        ddrDDdvdr=-((2*(-u*uDot*v-diffU2*vDot)*(u*uDot+v*vDot))/w4);
        ddrDDduDotdr=((2*(uDot-uDot*v2+u*v*vDot))/w2);
        ddrDDdvDotdr=-(-2*u*uDot*v-2*diffU2*vDot)/w2;

        ddrDDdrdu=ddrDDdudr;
        ddrDDdudu=-(2*r*(-uDot2*(1+3*u2-v2)*diffV2-2*u*uDot*v*(3+u2-3*v2)*vDot+v2*(-1-3*u2+v2)*vDot2))/w6;
        ddrDDdvdu=(2*r*(uDot*(1+u2-v2)+2*u*v*vDot)*(2*u*uDot*v+vDot-u2*vDot+v2*vDot))/w6;
        ddrDDduDotdu=(-4*r*u*uDot*(-diffV2)+2*r*v*(1+u2-v2)*vDot)/w4;
        ddrDDdvDotdu=(2*r*v*(uDot*(1+u2-v2)+2*u*v*vDot))/w4;

        ddrDDdrdv=ddrDDdvdr;
        ddrDDdudv=ddrDDdvdu;
        ddrDDdvdv=-(2*r*(u2*uDot2*(-1+u2-3*v2)+2*u*uDot*v*(-3+3*u2-v2)*vDot+diffU2*(-1+u2-3*v2)*vDot2))/w6;
        ddrDDduDotdv=(2*r*u*(2*u*uDot*v+vDot-u2*vDot+v2*vDot))/w4;
        ddrDDdvDotdv=(2*r*u*uDot*(1-u2+v2)+4*r*diffU2*v*vDot)/w4;

        ddrDDdrdRDot=0;
        ddrDDdudRDot=0;
        ddrDDdvdRDot=0;
        ddrDDduDotdRDot=0;
        ddrDDdvDotdRDot=0;

        ddrDDdrduDot=ddrDDduDotdr;
        ddrDDduduDot=ddrDDduDotdu;
        ddrDDdvduDot=ddrDDduDotdv;
        ddrDDduDotduDot=(2*r*diffV2)/w2;
        ddrDDdvDotduDot=(2*r*u*v)/w2;

        ddrDDdrdvDot=ddrDDdvDotdr;
        ddrDDdudvDot=ddrDDdvDotdu;
        ddrDDdvdvDot=ddrDDdvDotdv;
        ddrDDduDotdvDot=ddrDDdvDotduDot;
        ddrDDdvDotdvDot=(2*r*diffU2)/w2;

        %%%
        dduDDdrdr=-((4*rDot*uDot)/r^3);
        dduDDdudr=0;
        dduDDdvdr=0;
        dduDDdrDotdr=(2*uDot)/r2;
        dduDDduDotdr=(2*rDot)/r2;
        dduDDdvDotdr=0;

        dduDDdrdu=dduDDdudr;
        dduDDdudu=-(2*u*uDot2*(3+u2-3*v2)*diffV2-4*uDot*v*diffV2*(-1-3*u2+v2)*vDot+2*u*v2*(3+u2-3*v2)*vDot2)/w6;
        dduDDdvdu=-(2*u2*uDot2*v*(3+u2-3*v2)+4*u*uDot*(1-v4+u2*(-1+3*v2))*vDot+2*v*(1-u4+(-1+3*u2)*v2)*vDot2)/w6;
        dduDDdrDotdu=0;
        dduDDduDotdu=-(2*diffV2*(uDot*(1+u2-v2)+2*u*v*vDot))/w4;
        dduDDdvDotdu=(-4*u*uDot*v*diffV2-2*(diffU4+(-1+3*u2)*v2)*vDot)/w4;

        dduDDdrdv=dduDDdvdr;
        dduDDdudv=dduDDdvdu;
        dduDDdvdv=-(2*u*(u2*uDot2*(1-u2+3*v2)+2*u*uDot*v*(3-3*u2+v2)*vDot-diffU2*(-1+u2-3*v2)*vDot2))/w6;
        dduDDdrDotdv=0;
        dduDDduDotdv=(2*u2*(-2*u*uDot*v+u2*vDot-(1+v2)*vDot))/w4;
        dduDDdvDotdv=(2*u*(u*uDot*(-1+u2-v2)-2*diffU2*v*vDot))/w4;

        dduDDdrdRDot=dduDDdrDotdr;
        dduDDdudRDot=dduDDdrDotdu;
        dduDDdvdRDot=dduDDdrDotdv;
        dduDDdrDotdRDot=0;
        dduDDduDotdRDot=-(2/r);
        dduDDdvDotdRDot=0;

        dduDDdrduDot=dduDDduDotdr;
        dduDDduduDot=dduDDduDotdu;
        dduDDdvduDot=dduDDduDotdv;
        dduDDdrDotduDot=dduDDduDotdRDot;
        dduDDduDotduDot=-((2*u*diffV2)/w2);
        dduDDdvDotduDot=-(2*u2*v)/w2;

        dduDDdrdvDot=dduDDdvDotdr;
        dduDDdudvDot=dduDDdvDotdu;
        dduDDdvdvDot=dduDDdvDotdv;
        dduDDdrDotdvDot=dduDDdvDotdRDot;
        dduDDduDotdvDot=dduDDdvDotduDot;
        dduDDdvDotdvDot=-(2*u*diffU2)/w2;

        %%%
        ddvDDdrdr=-((4*rDot*vDot)/r^3);
        ddvDDdudr=0;
        ddvDDdvdr=0;
        ddvDDdrDotdr=(2*vDot)/r2;
        ddvDDduDotdr=0;
        ddvDDdvDotdr=(2*rDot)/r2;

        ddvDDdrdu=ddvDDdudr;
        ddvDDdudu=-(2*v*(uDot2*(-1+v)*(1+v)*(-1-3*u2+v2)+2*u*uDot*v*(3+u2-3*v2)*vDot+v2*(1+3*u2-v2)*vDot2))/w6;
        ddvDDdvdu=-(2*u*uDot2*(1-v4+u2*(-1+3*v2))+4*uDot*v*(1-u4+(-1+3*u2)*v2)*vDot+2*u*v2*(3-3*u2+v2)*vDot2)/w6;
        ddvDDdrDotdu=0;
        ddvDDduDotdu=(2*v*(2*u*uDot*(-diffV2)+v*(-1-u2+v2)*vDot))/w4;
        ddvDDdvDotdu=-((2*v2*(uDot*(1+u2-v2)+2*u*v*vDot))/w4);

        ddvDDdrdv=ddvDDdvdr;
        ddvDDdudv=ddvDDdvdu;
        ddvDDdvdv=-(2*u2*uDot2*v*(3-3*u2+v2)-4*u*diffU2*uDot*(-1+u2-3*v2)*vDot+2*diffU2*v*(3-3*u2+v2)*vDot2)/w6;
        ddvDDdrDotdv=0;
        ddvDDduDotdv=(-2*uDot*(diffV4+u2*(-1+3*v2))-4*u*diffU2*v*vDot)/w4;
        ddvDDdvDotdv=(2*diffU2*(-2*u*uDot*v+u2*vDot-(1+v2)*vDot))/w4;

        ddvDDdrdRDot=ddvDDdrDotdr;
        ddvDDdudRDot=ddvDDdrDotdu;
        ddvDDdvdRDot=ddvDDdrDotdv;
        ddvDDdrDotdRDot=0;
        ddvDDduDotdRDot=0;
        ddvDDdvDotdRDot=-(2/r);

        ddvDDdrduDot=ddvDDduDotdr;
        ddvDDduduDot=ddvDDduDotdu;
        ddvDDdvduDot=ddvDDduDotdv;
        ddvDDdrDotduDot=ddvDDduDotdRDot;
        ddvDDduDotduDot=-((2*v*diffV2)/w2);
        ddvDDdvDotduDot=-(2*u*v2)/w2;

        ddvDDdrdvDot=ddvDDdvDotdr;
        ddvDDdudvDot=ddvDDdvDotdu;
        ddvDDdvdvDot=ddvDDdvDotdv;
        ddvDDdrDotdvDot=ddvDDdvDotdRDot;
        ddvDDduDotdvDot=ddvDDdvDotduDot;
        ddvDDdvDotdvDot=((2*(-1+u^2)*v)/w2);

        aHess(:,:,1)=[0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      ddrDDdrdr,       ddrDDdudr,      ddrDDdvdr,      0,              ddrDDduDotdr,   ddrDDdvDotdr;
                      dduDDdrdr,       dduDDdudr,      dduDDdvdr,      dduDDdrDotdr,   dduDDduDotdr,   dduDDdvDotdr;
                      ddvDDdrdr,       ddvDDdudr,      ddvDDdvdr,      ddvDDdrDotdr,   ddvDDduDotdr,   ddvDDdvDotdr];
        aHess(:,:,2)=[0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      ddrDDdrdu,       ddrDDdudu,      ddrDDdvdu,      0,              ddrDDduDotdu,   ddrDDdvDotdu;
                      dduDDdrdu,       dduDDdudu,      dduDDdvdu,      dduDDdrDotdu,   dduDDduDotdu,   dduDDdvDotdu;
                      ddvDDdrdu,        ddvDDdudu,     ddvDDdvdu,      ddvDDdrDotdu,   ddvDDduDotdu,   ddvDDdvDotdu];
        aHess(:,:,3)=[0,               0,              0,              0               0,              0;
                      0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      ddrDDdrdv,       ddrDDdudv,      ddrDDdvdv,      0,              ddrDDduDotdv,   ddrDDdvDotdv;
                      dduDDdrdv,       dduDDdudv,      dduDDdvdv,      dduDDdrDotdv,   dduDDduDotdv,   dduDDdvDotdv;
                      ddvDDdrdv,       ddvDDdudv,      ddvDDdvdv,      ddvDDdrDotdv,   ddvDDduDotdv,   ddvDDdvDotdv];
        aHess(:,:,4)=[0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      ddrDDdrdRDot,    ddrDDdudRDot,   ddrDDdvdRDot,   0,              ddrDDduDotdRDot,ddrDDdvDotdRDot;
                      dduDDdrdRDot,    dduDDdudRDot,   dduDDdvdRDot,   dduDDdrDotdRDot,dduDDduDotdRDot,dduDDdvDotdRDot;
                      ddvDDdrdRDot,    ddvDDdudRDot,   ddvDDdvdRDot,   ddvDDdrDotdRDot,ddvDDduDotdRDot,ddvDDdvDotdRDot];
        aHess(:,:,5)=[0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      ddrDDdrduDot,    ddrDDduduDot,   ddrDDdvduDot,   0,              ddrDDduDotduDot,ddrDDdvDotduDot;
                      dduDDdrduDot,    dduDDduduDot,   dduDDdvduDot,   dduDDdrDotduDot,dduDDduDotduDot,dduDDdvDotduDot;
                      ddvDDdrduDot,    ddvDDduduDot,   ddvDDdvduDot,   ddvDDdrDotduDot,ddvDDduDotduDot,ddvDDdvDotduDot];
        aHess(:,:,6)=[0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      0,               0,              0,              0,              0,              0;
                      ddrDDdrdvDot,    ddrDDdudvDot,   ddrDDdvdvDot,   0,              ddrDDduDotdvDot,ddrDDdvDotdvDot;
                      dduDDdrdvDot,    dduDDdudvDot,   dduDDdvdvDot,   dduDDdrDotdvDot,dduDDduDotdvDot,dduDDdvDotdvDot;
                      ddvDDdrdvDot,    ddvDDdudvDot,   ddvDDdvdvDot,   ddvDDdrDotdvDot,ddvDDduDotdvDot,ddvDDdvDotdvDot];

        if(nargout>3)
            papt=zeros(6,1);
        end
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function varargout = cgal_insurftri(varargin)
% CGAL_INSURFTRI  Find whether a point is inside or outside a closed
% surface.
%
%   This function evaluates whether one or more points belong inside a
%   closed surface. First, we check whether the point is on the surface
%   itself (in that case, it's considered inside). If not, a ray is
%   projected from the point and the intersections with the surface are
%   counted. An odd number means that the point is inside. This approach
%   fails if the point is not on the surface, but the ray lies on the
%   surface or crosses a vertex, because this spans many arbitrary
%   intersections. To solve this problem, a few rays are used for each
%   point, and the majority vote decides whether it's inside or outside.
%
% ISIN = cgal_insurftri(TRI, X, XI)
%
%   TRI is a 3-column matrix. Each row represents the indices of the three
%   vertices that form a triangle. TRI as a whole represents the closed
%   surface.
%
%   X is a 3-column matrix. Each row represents the Cartesian coordinates
%   of a vertex on the surface, indexed by TRI values.
%
%   XI is a 3-column matrix. Each row represents the Carterian coordinates
%   of a point for which we want to find whether it's inside or outside the
%   closed surface. Note that if you want to test all the voxels in an
%   image, it is very slow and memory intensive to generate coordinates for
%   each voxel. In that scenario, it is much better to use the cell array
%   CI syntax shown below, and provide only values for the coordinate axes.
%
%   ISIN is a boolean vector with one element per point in XI. True means
%   that the corresponding point is inside the closed surface (or on the
%   surface), and false means that it's outside.
%
% ISIN = cgal_insurftri(TRI, X, CI)
%
%   CI is a cell array CI={XI, YI, ZI}, where XI, YI and ZI are row vectors
%   that describe a rectangular grid. For example,
%
%     CI={linspace(-.25, .25, 5), ...
%         linspace(-.25, .25, 4), ...
%         linspace(-.25, .25, 3)};
%
%   describes a sampling grid of 4 rows x 5 columns x 3 slices (note that
%   rows correspond to YI and columns to XI), of a domain
%   [-0.25, 0.25] x [-0.25, 0.25] x [-0.25, 0.25].
%
% ISIN = cgal_insurftri(..., DIRECTIONS, TOL)
%
%   DIRECTIONS is a 3-column matrix. Each row represents a vector with a
%   ray direction. By default,
%
%          DIRECTIONS=[ 1.0,  0.0,  0.0; ...
%                      -1.0,  1.0,  1.0; ...
%                      -1.0, -1.0, -1.0]
%
%   This default can fail with regular voxels, as rays may cross vertices.
%   A good practical alternative is to use a few random directions, e.g.
%
%          DIRECTIONS=rand(5, 3);
%
%   Warning! For the voting system to make sense, select an odd number of
%   rays.
%
%   TOL is a scalar with the distance tolerance. Points at distance <= TOL
%   are considered to be on the surface, and thus "inside". By default,
%   TOL=1e-15.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2012-2013 University of Oxford
% Version: 0.2.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

error('MEX file not found')

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% sweep MasterThres param: percent cells included in Auto0.7, as a function
% of masterthres. for fig4 main.

author: xiuyechen

repo: FishExplorer


text:
function i4vec_permute_uniform_test ( )

%*****************************************************************************80
%
%% I4VEC_PERMUTE_UNIFORM_TEST tests I4VEC_PERMUTE_UNIFORM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2015
%
%  Author:
%
%    John Burkardt
%
  n = 12;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_PERMUTE_UNIFORM_TEST\n' );
  fprintf ( 1, '  I4VEC_PERMUTE_UNIFORM randomly reorders an I4VEC.\n' );

  a = 100 + ( 1 : n );

  i4vec_print ( n, a, '  A, before rearrangement:' );

  [ a, seed ] = i4vec_permute_uniform ( n, a, seed );

  i4vec_print ( n, a, '  A, after rearrangement:' );

  return
end


text:
% parse all fields of a solver proto
% transfer all FC layers in the existing model to CONV layers (equvelent)
function [ Solver ] = SolverParser( solver_def_file, resume_file )
if ~exist(solver_def_file,'file')||isempty(solver_def_file)
    error('Solver definition file %s is not found.',solver_def_file);
end
if ~exist('resume_file','var')
    fprintf('Creating a face lmk solver structure based on %s.', solver_def_file);
    Solver = [];
else
    load(resume_file);
end
fout = fopen(solver_def_file,'r');
tline = fgetl(fout);
while ischar(tline)
    disp(tline)
    pos = 1;
    while tline(pos)==' '
       pos = pos + 1;
    end
    if tline(pos) == '#'
            tline = fgetl(fout);
            continue;
    end
    ind = find(tline == '"',1);
    if  ~isempty(ind)
        field = tline(ind + 1 : end - 1);
        ind2 = find(tline == ':',1);
        name = tline(1:ind2-1);
    else
        ind2 = find(tline == ':',1);
        if isempty(ind2)
            error('incorrect format.')
        end
        ctr = tline(ind2+2:end);
        if isempty(str2num(ctr))
            ctr(isspace(ctr)) = [];
            field = ctr;
        else
            field = str2double(ctr);
        end
        name = tline(1:ind2-1);
    end
    Solver = setfield(Solver, name, field);
    tline = fgetl(fout);
end
fclose(fout);
if ~isfield(Solver, 'solver_mode')
    Solver.solver_mode = 'GPU';
end
if ~isfield(Solver, 'device_id') && strcmp(Solver.solver_mode, 'GPU')
        Solver.device_id = 0;
end
if isfield(Solver,'model')
    lnum = length(Solver.model);
    for ind = 1:lnum
        if strcmp(Solver.model(ind).layer_names,sprintf('fc%d',ind))
            Solver.model(ind).layer_names = sprintf('conv%d',ind);
            weights = Solver.model(ind).weights{1};
            [s1,s2] = size(weights);
            [~,~,~,ch] = size(Solver.model(ind-1).weights{1});
            filtersize = sqrt(s1/ch);
            weights = reshape(weights,[filtersize,filtersize,ch,s2]);
            Solver.model(ind).weights{1} = weights;
        end
    end
end
end

author: rwenqi

repo: GFN-dehazing

save_path: github-repos/MATLAB/rwenqi-GFN-dehazing


text:
%kkxv 'XV, image visualization and color edit, by J. Bradley '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kxv.pane file
%
% Parameters:
% InputFile: i 'Input ', optional: 'First Input data object'
% Toggle: kdfin 'KDF input compatibility', default: 0: 'accepts KDF input file'
% OutputFile: o 'Output', optional: 'Resulting output data object'
% Toggle: kdfout 'KDF output compatibility', default: 0: 'converts to KDF after calling xv'
%
% Example: o = kkxv(i, {'i','';'kdfin',0;'o','';'kdfout',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kxv - XV, image visualization and color edit, by J. Bradley
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc.  All rights reserved.
%

function varargout = kkxv(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,..] = kkxv(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'kdfin', 0;'o', '__output';'kdfout', 0};
maxval={1,0,1,0};
minval={1,0,1,0};
istoggle=[1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','Toggle','OutputFile','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; 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 'kxv"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function centroids = kMeansInitCentroids(X, K)
%KMEANSINITCENTROIDS This function initializes K centroids that are to be
%used in K-Means on the dataset X
%   centroids = KMEANSINITCENTROIDS(X, K) returns K initial centroids to be
%   used with the K-Means on the dataset X
%

% You should return this values correctly
centroids = zeros(K, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: You should set centroids to randomly chosen examples from
%               the dataset X
%
% Initialize the centroids to be random examples

% Randomly reorder the indices of examples
randidx = randperm(size(X, 1));
% Take the first K examples as centroids
centroids = X(randidx(1:K), :);

% =============================================================

end

author: loserChen

repo: Coursera-MachineLearning


text:
function vis_next_pred(joints, next_joints, graph)

    if nargin < 3
        next_joints = [];
    end

    colors = {'r','g','b','c','m','y'};
    markers = {'+','o','s', 'x','.'}; %,'-','s','d'};
    styles = {};
    j = 1;
    for m = 1:length(markers)
        for c = 1:length(colors)
            styles{j} = [colors{c} markers{m}];
            j = j + 1;
        end
    end

    visJoints = false;
    if visJoints
        for jj = 1:size(joints, 1)
            %fprintf('joint %d\n', jj);
            if isnan(joints(jj, 1))
                continue;
            end
            plot(joints(jj, 1), joints(jj, 2), styles{jj}, 'MarkerSize',5);
            text(joints(jj, 1), joints(jj, 2), num2str(jj), 'Color', 'g', 'FontSize', 12);
            %pause;
        end
    end

    starts = [2];
    ends = [1 3 5 6 14];
    %starts = [10];
    %ends = [3 9 14 11 1];

    if ~isempty(next_joints)
        %{
        for jj = 1:size(next_joints, 1)
            pts = zeros(2,2);
            pts(1,:) = joints(jj, :);
            pts(2,:) = joints(jj, :) + next_joints(jj,:);
            line(pts(:,1), pts(:,2), 'Color', 'r');
        end
        %}

        for jj = 1:size(next_joints, 1)
            pts = zeros(2,2);
            start = graph(jj, 1);
            end_  = graph(jj, 2);

            %if ~(start == 2 && (end_ == 1 || end_ == 3 || end_ == ))
            %    continue;
            %end
            if ~(any(start == starts) && any(end_ == ends))
                continue;
            end

            [~, idx] = find(ends == end_);
            color = colors{idx};

            %if ~((end_ == 7 ))
            %    continue;
            %end
            orig = joints(start, :);
            pts(1,:) = orig;
            pts(2,:) = orig + next_joints(jj,:);
            line(pts(:,1), pts(:,2), 'Color', color, 'LineWidth', 6);
            %{
            if jj <= 13
                line(pts(:,1), pts(:,2), 'Color', 'b');
            else
               line(pts(:,1), pts(:,2), 'Color', 'r');
            end
            %}
        end
    end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [zp,exitFlag]=nearestPointInEllipsoid(z,A,p,gammaVal,method,varargin)
%%NEARESTPOINTONELLIPSOID Given an ellipsoid such that a point zp on the
%        surface of the ellipsoid satisfies the equation
%        (z-zp)'*A*(z-zp)=gammaVal find the point on or inside the
%        ellipsoid that is closest to another point p. This function
%        differs from nearestPointOnEllipsoid in that it will return the
%        original point if the point is inside of the ellipsoid.
%
%INPUTS: z  The numDimX1 center of the ellipsoid.
%        A  A numDimXnumDim symmetric, positive definite matrix that
%           specifies the size and shape of the ellipse or ellipsoid, where
%           a point zp is on the ellipse/ellipsoid if
%           (zp-z)'*A*(zp-z)=gammaVal.
%        p  A numDimX1 point.
%  gammaVal The threshold for declaring a point to be in the ellipsoid. If
%           this parameter is omitted or an empty matrix is passed, the
%           default value of 1 is used.
%    method Two algorithms are available to solve the problem. Possible
%           values are:
%           0 (The default if omitted or an empty matrix is passed.
%             Determine whether the point is in the ellipse; if so, return
%             it. Otherwise, call the function nearestPointOnEllipsoid.
%           1 Solve the problem by refomrulating it as a quadratic
%             programming problem with quadratic constrains and use the
%             constrainedLSSpher function.
% varargin When using method 0, this just just be 'epsVal' followed by a
%          number to set the epsVal parameter of the
%          nearestPointOnEllipsoid function, if it is used. Otherwise, if
%          method =1, this is any parameters that one wishes to pass to the
%          fzero function. These are generally comma-separated things, such
%          as 'TolX',1e-12.
%
%OUTPUTS: zp The numDimX1 point on the ellipse that is closest to p. If the
%            algorithm fails, then an empty matrix is returned.
%  exitFlag The exit flag from the fzero function, if used. Otherwise, if
%           the point is inside of the ellipsoid, this is 0.
%
%If we say that B=A/gammaVal, then points on the surface of the ellipsoid
%satisfy the equation
%(zp-z)'*B*(zp-z)<=1
%We want to minimize the squared distance
%s^2=(zp-p)'*(zp-p)
%Subject to the point being on the ellipsoid. Let L be the lower-triangular
%Cholesky decomposition of B and say that
%y=L'*(zp-z)
%which is equivalent to
%zp=z+inv(L)'*y
%Then the constraint for a point being on the ellipsoid is  y'*y<=1 and,
%after dropping constant terms, the squared distance that we want to
%minimize becomes
%(inv(L)'*y+z-p)'*(inv(L)'*y+z-p)
%Thus, the problem is a quadratic programming problem with a quadratic
%constraint (that y is unit magnitude). If method=1, then this problem is
%solved using the constrainedLSSpher function.
%
%EXAMPLE 1:
%This is a 2D example that can be easily plotted. The point is outside the
%ellipse.
% A=[0.1,0;
%    0,10];
% z=[0;10];
% p=[6;10];
% M=[0.413074198133900,  0.910697373904216;
%    -0.910697373904216,   0.413074198133900];%A rotation matrix
% A=M*A*M';
% [zp,exitCode]=nearestPointInEllipsoid(z,A,p);
% figure(1)
% clf
% hold on
% axis([-6,6,-6+10,6+10])
% axis square
% drawEllipse(z,A,1,'g','linewidth',2)
% plot([p(1),zp(1)],[p(2),zp(2)],'--c')
% scatter(p(1),p(2),'ok','linewidth',2)
% scatter(zp(1),zp(2),'xr','linewidth',2)
%
%EXAMPLE 2:
%This is another 2D example. In this instance, the point is inside the
%ellipse, so the point itself is the closest point on or in the ellipse.
% A=[1,0;
%    0,5];
% z=[0;10];
% p=[0.5;10];
% [zp,exitCode]=nearestPointInEllipsoid(z,A,p);
% figure(1)
% clf
% hold on
% axis([-6,6,-6+10,6+10])
% axis square
% drawEllipse(z,A,1,'g','linewidth',2)
% plot([p(1),zp(1)],[p(2),zp(2)],'--c')
% scatter(p(1),p(2),'ok','linewidth',2)
% scatter(zp(1),zp(2),'xr','linewidth',2)
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(gammaVal))
    gammaVal=1;
end

if(nargin<5||isempty(method))
   method=0;
end

switch(method)
    case 0%Use the polynomial method.
        opts.epsVal=[];
        opts=addCelListToStruct(opts,varargin);
        epsVal=opts.epsVal;

        diff=(z-p);
        if(diff'*A*diff-gammaVal<=0)
            exitFlag=0;
            zp=z;
            return;
        else
            zp=nearestPointOnEllipsoid(z,A,p,gammaVal,epsVal);
            if(~isempty(zp))
                exitFlag=1;
            else
                exitFlag=-1;
            end
        end
    case 1%Use the constrainedLSSpher function rather than the polynomial
          %method.
        B=A/gammaVal;
        L=chol(B,'lower');

        A=inv(L)';
        b=p-z;

        [y,~,exitFlag]=constrainedLSSpher(A,b,1,varargin{:});

        if(exitFlag==0)
           exitFlag=1;
        end

        if(~isempty(y))
            zp=z+A*y;
        else
            zp=[];
        end
    otherwise
        error('Unknown method specified.')
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
function element_stress = elementStress(matrixD, element_strain, element_num)
% 该函数用于计算单元应力
% 输入为
%   D矩阵 matriD
%   单元节点应变 element_strain
%   单元编号element_num
% 输出为
%   单元应力 element_stress

element_stress = matrixD * element_strain(:, element_num);

author: Meelfy

repo: FEM

save_path: github-repos/MATLAB/Meelfy-FEM


text:
function n = normalDirection(L)
% normal direction

[~,~,n] = svd(matrix(L));

n = vector3d(n(:,1));


text:
function var = kernGetVariance(kern)

% KERNGETVARIANCE Get the signal associated with a the kernel.

% KERN

fhandleStr = ([kern.type 'KernGetVariance']);
if exist(fhandleStr)==2
  fhandle = str2func(fhandleStr);
  var = fhandle(kern);
  return
else
  if(strcmp(kern.type, 'white'))
    var = 0;
    return
  else
    if(isfield(kern, 'variance'))
      var = kern.variance;
    else
      error(['No field variance or file ' fhandleStr])
    end
    return
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ Vehicle ] = gen_trajectory( t )

     radius = 10;
    Vehicle.position = radius* [ 5*cos(0.3*t/4); 4*sin(0.2*t/4); 2*sin(0.2*t/4+1) ];
    %Vehicle.position = radius* [ 5*cos(0.03*t); 4*sin(0.02*t); 2*sin(0.02*t+1) ];
    Vehicle.euler = [0.5*t+2; -0.3*t; 0.4*t];

    % generate rotation matrix
    for i = 1:size(Vehicle.euler, 2)
        Vehicle.orientation((i-1)*3+1:i*3, 1:3) = euler2rotation_matrix(Vehicle.euler(:, i));
    end
end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function [nodeNo,nodeChildIdx] = depthIndex2NodeNo(d,k,wt)
%DEPTHINDEX2NODENO Get node from depth and index in the tree
%   Usage: [nodeNo,nodeChildIdx] = depthIndex2NodeNo(d,k,wt)
%
%   `[nodeNo,nodeChildIdx] = depthIndex2NodeNo(d,k,wt)` returns node
%   *nodeNo* and an array of its children nodes *nodeChildIdx* positioned
%   in depth *g* and index *k* in the tree *wt*.
%
if(d==0)
    nodeNo=0;
    nodeChildIdx=0;
    return;
end

% find ordered nodes at depth d-1
nodesNo = getNodesInDepth(d,wt);
if(isempty(nodesNo))
   error('%s: Depth of the tree is less than given d.',mfilename);
end

% k is index in children of ordered nodes at depth d

nodeNo = zeros(numel(k),1);
nodeChildIdx = zeros(numel(k),1);
chNo = cumsum(cellfun( @(nEl) length(nEl.g),wt.nodes(nodesNo)));
chNoZ = [0;chNo(:)];

for kIdx=1:numel(k)
    ktmp = k(kIdx);
    idx = find(chNo>ktmp,1);
    if isempty(idx)
       error('%s: Index k=%i out of bounds.',mfilename,ktmp);
    end
    nodeNo(kIdx) = nodesNo(idx);
    nodeChildIdx(kIdx) = ktmp-chNoZ(idx)+1;
end

function nodd = getNodesInDepth(d,wt)
% find all nodes with d steps to the root ordered
if d==1
    % return root
    nodd = find(wt.parents==0);
    return;
end

nbf = nodeBForder(0,wt);
nbfTmp = nbf;
tempd = 0;
while tempd<d
    nbf(nbfTmp==0) = [];
    nbfTmp(nbfTmp==0) = [];
    nbfTmp = wt.parents(nbfTmp);
    tempd = tempd+1;
end
nodd = nbf(nbfTmp==0);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function value = c8_log ( z )

%*****************************************************************************80
%
%% C8_LOG evaluates the logarithm of a C8.
%
%  Discussion:
%
%    A C8 is a complex value.
%
%    Here we use the relationship:
%
%      C8_LOG ( Z ) = LOG ( MAG ( Z ) ) + i * ARG ( Z )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, complex Z, the argument.
%
%    Output, complex VALUE, the function value.
%
  arg = c8_arg ( z );
  mag = c8_mag ( z );

  value = log ( mag ) + arg * i;

  return
end

repo: jburkardt-m


text:
function [trl] = read_fcdc_trl(fn);

% READ_FCDC_TRL reads trial definitions from a file
%
% Given a file which defines N trials, this function returns a Nx3
% matrix with the begin latency, end latency, and the latency offset
% of the first sample of each trial. The latencies are in seconds.
%
% [trl] = read_fcdc_trl(filename)
%
% An FCD trial definition file is formatted like
%   begin   end     offset
%   0.0000  1.0000  0.0000
%   3.0000  4.0000  0.0000
%   5.0000  5.5000  0.0000
%   ...
%
% The trial begin and end are given in seconds relative to the start
% of the recorded datafile. The offset is given in seconds and indicates
% the latency of the first sample, relative to the trial marker or
% trigger. E.g., given a trigger at 7000ms (relative to the recording
% begin), a trial of 1000ms with a pretrigger interval of 300ms would
% correspond to "6.700 7.700 -0.300".

% Copyright (C) 2003, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: read_fcdc_trl.m 945 2010-04-21 17:41:20Z roboos $

fid = fopen(fn, 'rt');
if fid<0
   error('could not open file');
end

trl = [];
while ~feof(fid)
  tmp = fscanf(fid, '%f %f %f', 3);
  if ~isempty(tmp)
    trl = [trl; tmp'];
  end
end

fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
function [intersects, edgeIndices] = intersectLinePolyline(line, poly, varargin)
%INTERSECTLINEPOLYLINE Intersection points between a line and a polyline
%
%   P = intersectLinePolyline(LINE, POLY)
%   Returns the intersection points of the lines LINE with polyline POLY.
%   LINE is a 1-by-4 row vector containing parametric representation of the
%   line (in the format [x0 y0 dx dy], see the function 'createLine' for
%   details).
%   POLY is a NV-by-2 array containing coordinates of the polyline vertices
%   P is a K-by-2 array containing the coordinates of the K intersection
%   points.
%
%   P = intersectLinePolyline(LINE, POLY, TOL)
%   Specifies the tolerance for geometric tests. Default is 1e-14.
%
%   [P INDS] = intersectLinePolyline(...)
%   Also returns the indices of edges involved in intersections. INDS is a
%   K-by-1 column vector, such that P(i,:) corresponds to intersection of
%   the line with the i-th edge of the polyline. If the intersection occurs
%   at a polyline vertex, the index of only one of the two neighbor edges
%   is returned.
%   Note that due to numerical approximations, the use of function
%   'isPointOnEdge' may give results not consistent with this function.
%
%
%   Examples
%   % compute intersections between a square and an horizontal line
%     poly = [0 0;10 0;10 10;0 10];
%     line = [5 5 1 0];
%     intersectLinePolyline(line, poly)
%     ans =
%           10     5
%     % also return indices of edges
%     [inters inds] = intersectLinePolyline(line, poly)
%     inters =
%           10     5
%     inds =
%           2
%
%   % compute intersections between a square and a diagonal line
%     poly = [0 0;10 0;10 10;0 10];
%     line = [5 5 1 1];
%     intersectLinePolyline(line, poly)
%     ans =
%            0     0
%           10    10
%
%   See Also
%   lines2d, polylines2d, intersectLines, intersectLinePolygon
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 31/10/2003.
%

%   HISTORY
%   2008-11-24 rename 'pi' as 'intersects', update doc
%   2009-07-23 removed forgotten occurence of 'pi' variable (thanks to Bala
%       Krishnamoorthy)
%   2010-01-26 rewrite using vectorisation
%   2011-05-20 returns unique results
%   2011-07-20 returns intersected edge indices
%   2012-11-33 add 'diag' option for linePosition

% get computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% create the array of edges
N = size(poly, 1);
edges = [poly(1:N-1, :) poly(2:N, :)];

% compute intersections with supporting lines of polyline edges
supportLines = edgeToLine(edges);
intersects = intersectLines(line, supportLines, tol);

% find edges that are not parallel to the input line
inds = find(isfinite(intersects(:, 1)));

% compute position of intersection points on corresponding lines
pos = linePosition(intersects(inds, :), supportLines(inds, :), 'diag');

% and keep only intersection points located on edges
b = pos > -tol & pos < 1+tol;
inds = inds(b);
intersects = intersects(inds, :);

% remove multiple vertices (can occur for intersections located at polyline
% vertices)
[intersects, I, J] = unique(intersects, 'rows'); %#ok<ASGLU>

if nargout > 1
    % return indices of edges involved in intersection
    % (in case of intersection located at a vertex, only one of the
    % neighbor edges is returned)
    edgeIndices = inds(I);
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function studyPath = studySave(studies);
% Save a set of studies in the local studies.mat file.
%
% studyPath = studySave(studies);
%
% Returns the path to the updated studies file. This should reside
% in the same directory as the study-related MATLAB functions. mrVista 2
% is intended to have a unique set of studies for each copy of the code,
% in order to be unambiguous and minimize the selection of files. (If it's
% decide to change this, these functions will needed to be significantly
% expanded.)
%
% ras, 07/06
if nargin<1, error('Not enough input args.'); end

parentDir = fileparts(which(mfilename)); % parent directory of this M-file.
studyPath = fullfile(parentDir, 'studies.mat');

save(studyPath, 'studies');

return


text:
close all;
clc;
clearvars;
rng('default');

try_spx = true;
try_ehsan = true;
try_cvx = false;

% dimension of ambient space
n = 100;
% number of subspaces = number of clusters
ns = 2;
% dimensions of individual subspaces 1 and 2
d1  = 2;
d2 = 2;
% number of signals in individual subspaces
s1 = 100;
s2 = 100;
% total number of signals
s = s1 + s2;
% A random basis for first subspace;
basis1 = randn(n,d1);
% An origin for first subspace
origin1 = randn(n, 1);
% coefficients for s1 vectors chosen randomly in subspace 1
coeffs1 = randn(d1,s1);
% Random signals from first subspace
Y1 = basis1 * coeffs1 + origin1;
% A random basis for second subspace
basis2 = randn(n,d2);
% An origin for second subspace
origin2 = randn(n, 1);
% coefficients for s2 vectors chosen randomly in subspace 2
coeffs2 = randn(d2,s2);
% Random signals from first subspace
Y2 = basis2 * coeffs2 + origin2;
% Prepare the overall set of signals
Y = [Y1 Y2];
% ground through clustering data
true_labels = [1*ones(s1,1) ; 2*ones(s2,1)];
% the largest dimension amongst all subspaces
K = max(d1, d2);
% All signals are expected to  have a K-sparse representation

% This is an affine problem
affine = true;

if try_spx
    fprintf('Attempting our implementation of SPR-ADMM-Affine\n');
    tstart = tic;
    options.affine = true;
    [C1, details] = spx.cluster.ssc.spr_admm(Y, options);
    elapsed_time = toc(tstart);
    fprintf('Maximum difference: %.4f\n', max(max(abs(Y - Y*C1))));
    fprintf('Number of iterations: %d\n', details.iterations);
    fprintf('Elapsed time: %.4f seconds\n', elapsed_time);
    disp(details);
    figure;
    imshow(C1);
end

if try_ehsan
    fprintf('Attempting admmLasso_mat_func\n');
    tstart = tic;
    [C2, details] = admmLasso_mat_func(Y, affine);
    elapsed_time = toc(tstart);
    fprintf('Maximum difference: %.4f\n', max(max(abs(Y - Y*C2))));
    fprintf('Number of iterations: %d\n', details.iterations);
    fprintf('Elapsed time: %.4f seconds\n', elapsed_time);
    disp(details);
    figure;
    imshow(C2);
end

if try_cvx
    fprintf('Attempting CVX\n');
    tstart = tic;
    options.verbose = 1;
    options.affine = affine;
    [C3, details] = spx.cluster.ssc.spr_cvx(Y, options);
    elapsed_time = toc(tstart);
    fprintf('Maximum difference: %.4f\n', max(max(abs(Y - Y*C3))));
    fprintf('Elapsed time: %.4f seconds\n', elapsed_time);
    disp(details);
    figure;
    imshow(C3);
end

author: indigits

repo: sparse-plex


text:
function a = transpose(a)
%TRANSPOSE    Implements  a.'  for Taylor
%
%   c = a.'
%

% written  05/21/09     S.M. Rump
%

  index = reshape(1:prod(a.size),a.size)';
  a.size = fliplr(a.size);
  a.t = a.t(:,index);

repo: OpenMAS


text:
  function sino = fbp2_sino_weight(sg, sino, varargin)
%|function sino = fbp2_sino_weight(sg, sino, varargin)
%|
%| Apply sinogram weighting for first step of 2D fan-beam FBP.
%| This matlab version is the backup alternative for users lacking mex routine.
%|
%| in
%|	sg			sino_geom()
%|	ig			image_geom()
%|	sino	[nb,na,(L)]	fan-beam sinogram(s) (line integrals)
%| out
%|	sino	[nb,na,(L)]	weighted sinogram
%|
%| Copyright 2006-4-19 by Jeff Fessler, University of Michigan

if nargin == 1 && streq(sg, 'test'), fbp2_sino_weight_test, return, end
if nargin < 2, ir_usage, end

arg.chat = 0;
arg = vararg_pair(arg, varargin);

idim = size(sino);
sino = reshape(sino, idim(1), idim(2), []); % [nb,na,*L]
sino = fbp2_sino_weight_do(sino, sg);
sino = reshape(sino, idim); % [nb,na,(L)]

%
% fbp2_sino_weight_do()
%
function sino = fbp2_sino_weight_do(sino, sg)

if isinf(sg.dfs)
	sino = fbp2_sino_weight_flat(sino, ...
		sg.s, sg.dsd, sg.dso, sg.source_offset);
elseif sg.dfs == 0
	sino = fbp2_sino_weight_arc(sino, ...
		sg.s, sg.dsd, sg.dso, sg.source_offset);
else
	error 'only flat and arc done'
end

%
% fbp2_sino_weight_arc()
%
function sino = fbp2_sino_weight_arc(sino, ss, dsd, dso, source_offset);
na = size(sino,2);
nz = size(sino,3);
gam = ss / dsd;
w1 = abs(dso * cos(gam) - source_offset * sin(gam)) / dsd; % 1D weighting
sino = sino .* repmat(w1, [1 na nz]);

%
% fbp2_sino_weight_flat()
%
function sino = fbp2_sino_weight_flat(sino, ss, dsd, dso, source_offset);
na = size(sino,2);
nz = size(sino,3);
gam = atan(ss / dsd);
w1 = abs(dso * cos(gam) - source_offset * sin(gam)) / dsd; % 1D weighting
sino = sino .* repmat(w1, [1 na nz]);

%
% fbp2_sino_weight_test()
%
function fbp2_sino_weight_test
sg = sino_geom('ge1', 'down', 4);
sino = sg.ones;
s1 = fbp2_sino_weight(sg, sino);

im pl 1 2
im(1, s1), cbar
%max_percent_diff(s1,s2)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function Model = buildTmodel(Model)
% Initiates the template model (Tmodel). It adds model
% designations before some of the cell arrays add creates indenty arrays
% that indicate which reactions and metabolties are contained in which
% models. This function should be used after `verifyModel`.
%
% USAGE:
%
%    Tmodel = buildTmodel(Model)
%
% INPUTS:
%    Model:     Cobra model.
%
% OUTPUTS:
%    Model:     Template model which to which another model can be
%               compared. Known as `Tmodel` in other scripts
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

modelName = Model.description ; % Declare variables.
nameFields = {'rxnID' 'metID' 'subSystems' 'rxnReferences' ...
              'rxnECNumbers' 'grRules'} ;

%% Put the model name in for appropriate fields.
for iF = 1:length(nameFields)
    % Determine which entries are non empty.
    haveData = find(~cellfun(@isempty, Model.(nameFields{iF}))) ;
    Model.(nameFields{iF})(haveData) = ...
        strcat(modelName, ':', Model.(nameFields{iF})(haveData)) ;
end

%% Create indentity structure.
Model.Models.(modelName).rxns = true(length(Model.rxns), 1) ;
Model.Models.(modelName).mets = true(length(Model.mets), 1) ;

%% Add genes.
Model.Models.(modelName).genes = Model.genes ;

%% Move bound info into structures.
lb = Model.lb ;
ub = Model.ub ;
c = Model.c ;

Model.lb = ([]);
Model.ub = ([]);
Model.c = ([]);

Model.lb.(modelName) = lb ;
Model.ub.(modelName) = ub ;
Model.c.(modelName) = c ;

%% Make S sparse.
Model.S = sparse(Model.S) ;

%% Remove extra fields that could get in the way.
unwantedFields = {'description' 'genes' 'rxnGeneMat' 'rules'} ;
for iField = 1:length(unwantedFields)
    try
        Model = rmfield(Model, unwantedFields{iField}) ;
        %fprintf(['Removed ' unwantedFields{iField} ' from in Tmodel\n'])
    catch
        %fprintf([unwantedFields{iField} ' not in Tmodel\n'])
    end
end

author: opencobra

repo: cobratoolbox


text:
function cmp = MS_complexity(x,n,binHow)
% Calculate the Lempel-Ziv complexity of the n-bit encoding of x.
%
% cmp is the normalised complexity, that is the number of distinct
% symbol sequences in x, divided by the expected number of distinct
% symbols for a noise sequence.
%
% Algorithm is implemented in complexitybs.c
%
% Michael Small
% michael.small@uwa.edu.au, http://school.maths.uwa.edu.au/~small/
% 7/10/04
% For further details, please see M. Small. Applied Nonlinear Time Series
% Analysis: Applications in Physics, Physiology and Finance. Nonlinear Science
% Series A, vol. 52. World Scientific, 2005. (ISBN 981-256-117-X) and the
% references therein.

if nargin < 2
    n = 2;
end
if nargin < 3
    binHow = 'equiprobable';
end
%-------------------------------------------------------------------------------

if length(n) > 1
    % Can run with multiple values of the number of symbols, n:
    for ni = 1:length(n)
        cmp(ni) = MS_complexity(x,n(ni),binHow);
    end
else

    % do the binning, with equiprobable bins
    switch binHow
    case 'equiprobable'
        x=x(:);
        nx=length(x);
        [xn,xi]=sort(x+eps*randn(size(x))); %introduce randomness for ties
        y=zeros(nx,1);
        y=1:nx;
        y=floor(y.*(n/(nx+1)));
        x(xi)=y;
    case 'equiwidth'
    % else,
    %     %do binning with equal width bins
        x=x(:);
        nx=length(x);
        minx=min(x);
        maxx=max(x);
        stepx=(maxx-minx)/n;
        y=zeros(nx,1);
        while minx<maxx
            minx=minx+stepx;
            y=y+double(x<minx);
        end
        x=floor(y);
    end

    %compute complexity with complexitybs
    cmp = MS_complexitybs(x);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function model = changelssvm(model,option, value)
% Change a field of the object oriented representation of the LS-SVM
%
%
% The different options of the fields are given in following table:
%
%     1. General options representing the kind of model:
%
%           type: 'classifier' ,'function estimation'
% implementation: 'CMEX' ,'CFILE' ,'MATLAB'
%         status: Status of this model ('trained'  or 'changed' )
%          alpha: Support values of the trained LS-SVM model
%              b: Bias term of the trained LS-SVM model
%       duration: Number of seconds the training lasts
%         latent: Returning latent variables ('no' ,'yes' )
%       x_delays: Number of delays of eXogeneous variables (by default 0 )
%       y_delays: Number of delays of responses (by default 0 )
%          steps: Number of steps to predict (by default 1 )
%            gam: Regularisation parameter
%    kernel_type: Kernel function
%    kernel_pars: Extra parameters of the kernel function
%
%
%     2. Fields used to specify the used training data:
%
%    x_dim: Dimension of input space
%    y_dim: Dimension of responses
%  nb_data: Number of training data
%   xtrain: (preprocessed) inputs of training data
%   ytrain: (preprocessed,coded) outputs of training data
% selector: Indexes of training data effectively used during training
%
%
%     3. Options used in the Conjugate Gradient (CG) algorithm:
%
%     cga_max_itr: Maximum number of iterations in CG
%         cga_eps: Stopcriterium of CG, largest allowed error
%    cga_fi_bound: Stopcriterium of CG, smallest allowed improvement
%        cga_show: Show the results of the CG algorithm (1 or 0)
% cga_startvalues: Starting values of the CG algorithm
%
%
%     4. Fields with the information for pre- and post-processing (only given if appropriate):
%
%  preprocess: 'preprocess'  or 'original'
%     schemed: Status of the preprocessing
%               ('coded' ,'original'  or 'schemed' )
% pre_xscheme: Scheme used for preprocessing the input data
% pre_yscheme: Scheme used for preprocessing the output data
%   pre_xmean: Mean of the input data
%    pre_xstd: Standard deviation of the input data
%   pre_ymean: Mean of the responses
%    pre_ystd: Standard deviation of the reponses
%
%
%     5. The specifications of the used encoding (only given if appropriate):
%
%          code: Status of the coding
%                 ('original' ,'changed'  or 'encoded')
%      codetype: Used function for constructing the encoding
%                  for multiclass classification (by default 'none')
% codetype_args: Arguments of the codetype function
%  codedist_fct: Function used to calculate to which class a
%                coded result belongs
% codedist_args: Arguments of the codedist function
%     codebook2: Codebook of the new coding
%     codebook1: Codebook of the original coding
%
% Full syntax
%
% >> model = changelssvm(model, field, value)
%
%       Outputs
%         model(*) : Obtained object oriented representation of the LS-SVM model
%       Inputs
%         model    : Original object oriented representation of the LS-SVM model
%         field    : Field of the model one wants to change (e.g. 'preprocess')
%         value    : New value of the field of the model one wants to change
%
% See also:
%   trainlssvm, initlssvm, simlssvm, plotlssvm.

% Copyright (c) 2010,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.ac.be/sista/lssvmlab

%
% alias sigma^2
%
if (strcmpi(option,'sig2')) option = 'kernel_pars'; end

%
% selector -> nb_data
% nb_data -> selector
%
if strcmp(option,'selector'),
  model.nb_data = length(value);
end
if strcmp(option,'nb_data'),
  model.selector = 1:value;
end

%
% xtrain
%
if strcmp(option,'xtrain'),
  [nb,model.x_dim] = size(value);
  model.nb_data = nb;%min(nb,model.nb_data);
  model.selector = 1:model.nb_data;
  if length(model.gam)>model.y_dim & length(model.gam)~=size(value,1),
    warning('Discarting  different gamma''s...');
    model.gam = max(model.gam);
  end
    eval('value=prelssvm(model,value);',...
       'warning(''new trainings inputdata not comform with used preprocessing'');');
end

%
% ytrain
%
if strcmp(option,'ytrain'),
  if size(value,2)~=size(model.ytrain,2),
    model.y_dim = size(value,2);
  end
  eval('value = codelssvm(model,[],value);',...
       'warning(''new trainings outputdata not comform with used encoding;'');');
  eval('[ff,value] = prelssvm(model,[],value);',...
       'warning(''new trainings outputdata not comform with used preprocessing;'');');
  [nb,model.y_dim] = size(value);
  model.nb_data = min(nb,model.nb_data);
  model.selector = 1:model.nb_data;
end

%
% switch between preprocessing - original data
% model.prestatus = {'changed','ok'}
%
if (strcmpi(option,'preprocess')) & model.preprocess(1)~=value(1),
  model.prestatus = 'changed';
end

%
% change coding
%
if strcmpi(option,'codetype') | strcmpi(option,'codebook2') | ...
      strcmpi(option, 'codeargs') | strcmpi(option, 'codedistfct'),
  model.code = 'changed';
elseif  strcmpi(option,'codebook1'),
  warning('change original format of the classifier; the toolbox will be unable to return results in the original format');
end

%
% final change
%
eval(['old_value = model.' lower(option) ';'],'old_value=[];');
eval(['model.' lower(option) '=value;']);

if (isempty(value) | isempty(old_value)),
  different = 1;
else
  eval('different = any(old_value~=value);','different=1;');
end

if different & ~strcmpi(option,'implementation'),
  model.status = 'changed';
end

author: peterhcharlton

repo: RRest


text:
%SIGM Sigmoid map
%
%   W = W*SIGM(SCALE)
%
%   B = SIGM(A,SCALE)
%   B = A*SIGM([],SCALE)
%   B = A*SIGM(SCALE)
%
% INPUT
%   A        Dataset (optional)
%   SCALE    Scaling parameter (optional, default: 1)
%
% OUTPUT
%   W        Sigmoid mapping, or
%   B        Dataset A mapped by sigmoid mapping
%
% DESCRIPTION
% Sigmoidal transformation, useful to transform a map to classifier,
% producing posterior probability estimates. The parameter SCALE scales the
% data first (A/SCALE), before the transformation. Default: SCALE = 1, i.e.
% no scaling.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, CLASSC

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function out = sigm (varargin)

  argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],1);

  if mapping_task(argin,'definition')
    out = define_mapping(argin,'fixed','Sigmoidal Mapping');

  else			% Evaluate

    [a,scale] = deal(argin{:});
		out = 1./(1+exp(-a/scale));

	end

	return

author: marianux

repo: ecg-kit

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


text:
function dat=readmptiff(fname)
%
% vol=readmptiff(fname)
%
% load a volume from a multi-page TIFF file
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%      fname: input file name
%
% output:
%      dat: output, data read from the TIFF file
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

info = imfinfo(fname);
nslice = length(info);
if(nslice<=0)
    error('no data found in the tiff');
end

for i = 1:nslice
    dat(:,:,i) = imread(fname, i);
end

author: fangq


text: %% unimodalComponent_index.m

author: mtex-toolbox

repo: mtex


text:
function  cplexProblem = buildCplexProblemFromCOBRAStruct(Problem)
% Build a cplex object from the given LP problem in COBRA Format
% USAGE:
%    cplexProblem = buildCplexProblemFromCOBRAStruct(LPproblem)
%
% INPUT:
%    LPproblem:         A COBRA style Problem with the following fields:
%                        * .A - The equality and in equality matrix
%                        * .b - The right hand side values of the constraints
%                        * .ub - the upper bounds of the variables
%                        * .lb - the lower bounds of the variables
%                        * .osense - The objective sense (-1 for max, 1 for min)
%                        * .c - the objective coefficient vector for the linear part
%                       OPTIONAL:
%                        * .F - The objective coefficient matrix for the quadratic part.
%                        * .varType - The variable types for mixed integer problems ('I', integer, 'C', continous,'B' binary)
%                        * .b_L - left hand sides of the constraints (will only be used if csense is empty)
%                        * .csense - The constraint senses, Default assumption is all 'E'
%                        * .x0 - Basis to use
%

try
    cplexProblem = Cplex();
    %cplexProblem = Cplex('COBRAProblem');
catch ME
    error('CPLEX not installed or licence server not up')
end
if (~isempty(Problem.csense))
    % build the rhs/lhs of the problem.
    boolE = Problem.csense == 'E';
    b_L(boolE) = Problem.b(boolE);
    b_U(boolE) = Problem.b(boolE);
    boolG = Problem.csense == 'G';
    b_L(boolG) = Problem.b(boolG);
    b_U(boolG) = inf;
    boolL = Problem.csense == 'L';
    b_L(boolL) = -inf;
    b_U(boolL) = Problem.b(boolL);
elseif isfield(Problem.csense, 'b_L') && isfield(Problem, 'b_U')
    % or extract them
    b_L = Problem.b_L;
    b_U = Problem.b_U;
else
    % or simply use the equality assumption.
    b_L = Problem.b;
    b_U = Problem.b;
end

cplexProblem.Model.A = Problem.A;
cplexProblem.Model.rhs = columnVector(b_U);
cplexProblem.Model.lhs = columnVector(b_L);
cplexProblem.Model.ub = Problem.ub;
cplexProblem.Model.lb = Problem.lb;

if isfield(Problem,'F')
    f = diag(Problem.F);
    bool0 = f==0;
    if any(bool0)
        Fbar = Problem.F - diag(f);
        if ~any(Fbar,'all')
            feasTol = getCobraSolverParams('LP', 'feasTol');
            f(bool0)=feasTol/10;
            f(~bool0)=0;
            Problem.F = Problem.F + spdiags(f,0,size(Problem.F,1),size(Problem.F,1));
            %fprintf('%s\n',['buildCplexProblemFromCOBRAStruct: Replacing zeros on the diagonal of QP problem.F with regularisation of ' num2str(feasTol/10)])
        end
    end
    cplexProblem.Model.Q = Problem.F;
end

if isfield(Problem,'c')
    cplexProblem.Model.obj = Problem.c;
end
if isfield(Problem,'osense')
    if Problem.osense == 1
        cplexProblem.Model.sense = 'minimize';
    else
        cplexProblem.Model.sense = 'maximize';
    end
else
    cplexProblem.Model.sense = 'minimize';
end

if isfield(Problem,'vartype')
    cplexProblem.Model.ctype = columnVector(Problem.vartype)';
end

if isfield(Problem,'x0')
    cplexProblem.Start.x = Problem.x0;
end


text:
function [results]=updateCD(text,parsed,note)
% Correct(Update) the species name according to a reference list of the names.
%
% USAGE:
%
%    [results] = updateCD(text, parsed, note)
%
% INPUTS:
%    text:            A matlab variable that saves the lines of the XML file.
%    parsed:          A parsed CD model structure generated by `parseCD` function.
%
% OPTIONAL INPUT:
%    note:            The default value is '0'; It can be set as 'true';
%                     Record the column number of the species that are to be
%                     replaced with a new name.
%
% OUTPUTS:
%    text:            The corrected lines of the XML file.
%    number:          The number of the line.
%    changdStrLines:  The changed line of text.
%    col:             Record the column number of the species that is to be
%                     replaced with a new species name.
%
% EXAMPLE:
%
%    [text, number, numberofchanges, col] = updateCD(annotedText, ref_corrected)

if nargin<3;
    ConductCol=0;
    col={};

elseif strcmp(note,'true')
    ConductCol=1;
    col={};
end

col_name=2; % the column 2 contains metabolite names.

ref=parsed.r_info;

pool={'<species (\w*)','<celldesigner:name>'}; % the format in CellDesiginer XML file
keywords={' name', ':nam'};
mark={'"','<'};

for i=1:length(pool);
     a=0;
    if i==2
        warning('starting the second round of the program');
    end

    str=pool{i}

    key=keywords{i}

    mark_key=mark{i}
    % str='<celldesigner:species(\w*) id=' ;

    for ln=1:length(text);

        line_species=regexpi(text(ln).str,str); % find the line that contains the string ('str').
        % dd=str2num({line_species})
        if ~isempty(line_species{1})
            % reText=cellstr(text(line).str)

            % if strcmpi(reText,ref.species{line,1})
            a=a+1;
            number(a,1)=ln;  % ln stores line number.
            % disp(line_species);%(number(a,1));
        end
    end
    if i==1;
        totalCount=length(number(:,1))
    else
        totalCount=totalCount;
    end

    changedStrLines={};
    nC=0;

    for d=1:totalCount  % d increase discretely from 1 to the end.

        nRxn=number(d,1); % d is the count, whereas nRxn contains the row number.
        if isfield(ref,'species');

            reText=cellstr(text(nRxn).str)
            %for nn=1:length(ref.species(:,1))
            sizeRef=size(ref.species);
            disp(sizeRef)
            disp(d);
            disp(col_name);
            disp(length(number(:,1)))

            if ~isempty(ref.species(d,col_name))
                %if strcmpi(text(line).str,ref.species{d,1});
                strd=strfind(reText,key) % ' name')
                if (~isempty(strd{1}))  % in some cases, the line does not contain ' name'.
                    %if (~isempty(str))

                    if iscell(reText) % convert cell type into string type
                        reText=reText{1};
                    else
                        error('the format of the reText is not right');
                    end
                    [p_st,p_ed]=position(reText,key,mark_key) % ' name')
                    fprintf(' the text: %s',reText);
                    ToBeRplaced=reText(p_st:p_ed);
                    % disp(ToBeRplaced);
                    %ToBeRplaced=ref.species{nn,3};
                    try
                        text(nRxn).str=strrep(text(nRxn).str,ToBeRplaced,ref.species{d,col_name})
                    catch
                        %                     disp(text(nRxn).str);
                        %                     disp(ToBeRplaced);
                        %                     dusp(ref.species{d,3});
                    end
                    nC=nC+1;
                    changedStrLines(nC,1)=text(nRxn).str;
                    fprintf('original text is %s and is replaced with %s',ToBeRplaced,ref.species{d,col_name});
                    if  ConductCol==1;
                    col{nC,1}=ToBeRplaced;
                    col{nC,2}=ref.species{d,col_name}
                    end
                end
                %end
            end
            %end
        else
            warning('cannot find the field name ''species''');
        end
    end
end

results.text=text;
results.number=number;
results.changedStrLines=changedStrLines;
results.col=col;

end

function [p_st,p_ed]=position(str_long,str_ID,mark)
%% name='metaid';

ind_pos=strfind(str_long,str_ID);
l=length(str_ID)+2;
try
    p_st=ind_pos(1)+l
catch
    error('cannot find the Keyword in the the line')
end
end_rem=strfind(str_long(p_st:end),mark) % '"');
p_ed=end_rem(1)+p_st-2;
% disp(ind_pos);
% disp(l);
% disp(p_ed)

end


text:
function [wh_predictors, betas, b, stat] = regress_best_subsets_ga(X, Y)
% :Usage:
% GA-based best subsets regression
%
% ::
%
%     [wh_predictors, betas, b_subset, stat_subset] = regress_best_subsets_ga(X, Y)
%
% :Inputs:
%
%   **Y:**
%        is outcome data
%
%   **X:**
%        is predictor matrix
%
%   **wh_predictors:**
%        is the primary outcome -- it is vector of which predictors to include in the model
%
% the objective criterion is AIC

seeds = get_seeds(X, Y);

% Example: get_AIC(X, Y, seeds(:, 1))

% The fitness function is higher when AIC is lower
% The 1/AIC transformation would introduce a nonlinear weighting function;
% so would -log(AIC)
% log(1/AIC) is approximately linear.

%fitness_fcn = @(wh_preds)  1 ./ get_AIC(X, Y, wh_preds)
fitness_fcn = @(wh_preds)  log(1 ./ get_AIC(X, Y, wh_preds));

% for i = 1:size(seeds, 2)
%     a(i) = get_AIC(X, Y, seeds(:, i));
%     f(i) = fitness_fcn(seeds(:, i));
% end

% Run GA

example_vec = {seeds(:, 1)}; % this is 'inputs', what will be passed to the objective function (ofun) in the GA

for i = 1:size(seeds, 2), gaseeds{i} = seeds(:, i); end

[best_params,fit,beff,in,isconverged] = tor_ga(300, 50, example_vec, fitness_fcn, 'discrete', 'seeds', gaseeds, 'genconverge', 5);

wh_predictors = logical(best_params{1});

 Xs = X(:, wh_predictors);
 [b, dev, stat] = glmfit(Xs, Y); % subset plus intercept first

 betas = zeros(size(X, 2) + 1, 1);
 betas(logical([1; wh_predictors])) = b;  % add the intercept always

end

function aic = get_AIC(X, Y, wh_preds)

    X = X(:, logical(wh_preds));
    [nobs, npreds] = size(X);

    [b, dev, stat] = glmfit(X, Y);
    RSS = stat.resid' * stat.resid;

    % AIC = 2k - 2ln(L),  k = npreds, L = likelihood
    % AIC = 2k + n(ln(2piRSS/n) + 1)
    %
    % simplifying to relative values (we need only relative AIC, not
    % absolute:

    npreds = sum(wh_preds) + 1;  % include intercept
    aic = 2 * npreds + nobs * log(RSS);

end

function seeds = get_seeds(X, Y)

[nobs, npreds] = size(X);
z = false(npreds, 1);

% seed the GA with stepwise regression results and univariate results (top
% 1, 2, ... k)

[b, dev, stat] = glmfit(X, Y);
b = abs(b(2:end));
t = abs(stat.t(2:end));

% Set of candidate subsets with descending order of abs magnitude of betas
[bs, wh] = sort(b, 1, 'descend');
bseeds = false(npreds);

for i = 1:npreds
    my_subset = z;
    my_subset(wh(1:i)) = true;
    bseeds(:, i) = my_subset;
end

% Set of candidate subsets with descending order of abs magnitude of
% t-values (last one is redundant with bseeds)
[ts, wh] = sort(t, 1, 'descend');
tseeds = false(npreds, npreds - 1);
for i = 1:npreds - 1
    my_subset = z;
    my_subset(wh(1:i)) = true;
    tseeds(:, i) = my_subset;
end

% single-predictor seeds with each variable
eseeds = logical(eye(npreds));

% stepwise seeds
pentervals = [.05:.05:.95];
for i = 1:length(pentervals)
    [b, se, p, inmodel] = stepwisefit(X, Y, 'penter', pentervals(i), 'display', 'off');
    fseeds(:, i) = inmodel';
end

seeds = [bseeds tseeds eseeds fseeds];

% random seeds
numseeds = max(100, 200 - size(seeds, 2)); % at least 100 seeds
rseeds = rand(npreds, numseeds);
rseeds = rseeds > .5;

seeds = [seeds rseeds];

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore