text:
function lspec = lyapspec(ts, dim, delay, NNR, eps, Nref)

ts = ts(:);
N = length(ts) - dim*delay;
ref = 1:delay:min(N,(1+(Nref-1)*delay));

[nn, dists, points] = emb_nn_search(ts(1:(end-delay)), [dim delay], ref, NNR+1, -1);

Q = diag(ones(dim,1));
lspec = [];zeros(dim, 1);

for i=1:length(ref)
	X = points(nn(i,:), :);
	X0 = mean(X);		% center of gravity
	Y = X - repmat(X0, NNR+1, 1);
	[U,W,V] = svd(Y);
	w = diag(W) / sum(diag(W));
	k = max(find(w >= eps));
	V = V(:,1:k);
	Z = Y * V;

	images = ts(dim*delay+nn(i,:));
	%all(points(nn(i,:)+delay,1) == images);

	a = [Z ones(NNR+1,1)] \ images;

	b = a(end);
	a = a(1:end-1);
	c = V * a;

	%mean(abs(images-(X * c + ( b - X0 * c))) ./ abs(images))

	JAC = sparse([ones(1,dim) 2:dim], [1:dim 1:dim-1], [c' ones(1, dim-1)]);
	%full(JAC);

	A =  JAC * Q;
	[Q,R,E] = qr(A);

	lspec = [lspec; diag(R)'];
end

lspec = mean(log(abs(lspec(min(length(ref)/10,500):end,:))));

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function C = getAutoCovMat(hmm,k,verbose)
% Get the autocovariance from a TDE-HMM model, for state k,
%
% Diego Vidaurre, OHBA, University of Oxford (2019)

if nargin < 3, verbose = 1; end
C = getFuncConn(hmm,k,verbose);

end

author: OHBA-analysis

repo: HMM-MAR


text:
function area = polygon_area_3d_2 ( n, v, area )

%*****************************************************************************80
%
%% POLYGON_AREA_3D_2 computes the area of a polygon in 3D.
%
%  Discussion:
%
%    The computation is not valid unless the vertices of the polygon
%    lie in a plane, so that the polygon that is defined is "flat".
%
%    The polygon does not have to be "regular", that is, neither its
%    sides nor its angles need to be equal.
%
%    The area is computed as the sum of the areas of the triangles
%    formed by the last node with consecutive pairs of nodes (1,2),
%    (2,3), ..., and (N-2,N-1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%
%    Input, real V(3,N), the coordinates of the vertices.
%
%    Output, real AREA, the area of the polygon.
%
  area_vector(1:3) = 0.0;

  for i = 1 : n - 2

    t(1:3,1:3) = [ v(1:3,i)'; v(1:3,i+1)'; v(1:3,n)' ]';

    area_vector_triangle = triangle_area_vector_3d ( t );

    area_vector(1:3) = area_vector(1:3) + area_vector_triangle(1:3);

  end

  area = 0.5 * sqrt ( sum ( area_vector(1:3).^2 ) );

  return
end

author: johannesgerer

repo: jburkardt-m


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

%%
clear; close all; clc;

%% Syntax
% |[Vs]=triSurfSmoothFourier(F,V,numFreq);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
siz = 512;
rad = 36;

% Find the maximum

imm = im;
imt = imm;
fil = zeros( siz, siz );

[val ind] = max( imt(:) );

while (val > 0)

x = floor( (ind-1) / siz ) + 1;
y = mod  ( (ind-1),  siz ) + 1;

fil = fil + mkCircMask( rad, [siz siz], [x y] ) .* val;

imt = imm - fil;
imt = imt.*(imt > 0);

[val ind] = max( imt(:) );

end

fmax = fil;

% Find the minimum

imm = max(max(im)) - im;
imt = imm;
fil = zeros( siz, siz );

[val ind] = max( imt(:) );

while (val > 0)

x = floor( (ind-1) / siz ) + 1;
y = mod  ( (ind-1),  siz ) + 1;

fil = fil + mkCircMask( rad, [siz siz], [x y] ) .* val;

imt = imm - fil;
imt = imt.*(imt > 0);

[val ind] = max( imt(:) );

end

fmin = max(max(im)) - fil;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = sumk(varargin)
% SUMK  Returns sum of k largest (eigen-)values.
%
% s = SUMK(X,k,w)
%
% For a vector X, SUMK returns the sum of the k largest elements.
%
% For a symmetric matrix X, SUMK returns the sum of the k largest eigen-values.
%
% A third argument can be used to generalize to weighted sum.
% The weights have to be non-negative and non-increasing.
%
% See also SUMABSK

if nargin == 1
    error('sumk needs at least two arguments');
end

switch class(varargin{1})

    case 'double' % What is the numerical value of this argument (needed for displays etc)

        X = varargin{1};
        [n,m] = size(X);
        if (min(n,m) > 1 & ~issymmetric(X))
            error('sumk can only be applied on vectors and symmetric matrices');
        else
            k = min(length(X),varargin{2});
            w = ones(k,1);
            if nargin == 3
                w = varargin{3};
                if length(w)==1
                    w = ones(k,1)*w;
                end
            end
            if min(n,m)==1
                sorted = sort(X,'descend');
            else
                sorted = sort(eig(X),'descend');
            end
            sorted = sorted(:);
            w = w(:);
            varargout{1} = sum(sorted(1:k).*w(1:k));
        end

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.
        X = varargin{1};
        [n,m] = size(X);
        if nargin < 3
            w = 1;
        else
            w = varargin{3};
            if length(w)>1
                if any(diff(w)>0) || any(w<0)
                    error('The weights have to be non-negative non-increasing')
                end
            end
        end
        if (min(n,m) > 1 & ~issymmetric(X))
            error('sumk can only be applied on vectors and symmetric matrices');
        else
            varargout{1} = yalmip('define',mfilename,varargin{:});
        end

    case 'char' % YALMIP send 'model' when it wants the epigraph or hypograph
        if isequal(varargin{1},'graph')
            t = varargin{2}; % Second arg is the extended operator variable
            X = varargin{3}; % Third arg and above are the args user used when defining t.
            k = min(varargin{4},length(X));
            if nargin == 5
                w = varargin{5};
            else
                w = 1;
            end
            [Model,Properties] = sumk_generator(X,k,t,w);
            varargout{1} = Model;
            varargout{2} = Properties;
            varargout{3} = X;
        else
        end
    otherwise
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function value = triangle_orientation_2d ( t )

%*****************************************************************************80
%
%% TRIANGLE_ORIENTATION_2D determines the orientation of a triangle in 2D.
%
%  Discussion:
%
%    Three distinct non-colinear points in the plane define a circle.
%    If the points are visited in the order P1, P2, and then
%    P3, this motion defines a clockwise or counterclockwise
%    rotation along the circle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 July 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Output, integer VALUE, reports if the three points lie
%    clockwise on the circle that passes through them.  The possible
%    return values are:
%    0, the points are distinct, noncolinear, and lie counterclockwise
%    on their circle.
%    1, the points are distinct, noncolinear, and lie clockwise
%    on their circle.
%    2, the points are distinct and colinear.
%    3, at least two of the points are identical.
%
  dim_num = 2;

  if ( all ( t(1:dim_num,1) == t(1:dim_num,2) ) | ...
       all ( t(1:dim_num,2) == t(1:dim_num,3) ) | ...
       all ( t(1:dim_num,3) == t(1:dim_num,1) ) )
    value = 3;
    return
  end

  det = ( t(1,1) - t(1,3) ) * ( t(2,2) - t(2,3) ) ...
      - ( t(1,2) - t(1,3) ) * ( t(2,1) - t(2,3) );

  if ( det == 0.0 )
    value = 2;
  elseif ( det < 0.0 )
    value = 1;
  elseif ( 0.0 < det )
    value = 0;
  end

  return
end


text:
function hsi = rgb2hsi(rgb)
%RGB2HSI Converts an RGB image to HSI.
%   HSI = RGB2HSI(RGB) converts an RGB image to HSI. The input image is
%   assumed to be of size M-by-N-by-3, where the third dimension
%   accounts for three image planes: red, green, and blue, in that
%   order. If all RGB component images are equal, the HSI conversion is
%   undefined. The input image can be of class double (with values in
%   the range [0, 1]), uint8, or uint16.
%
%   The output image, HSI, is of class double, where:
%     HSI(:,:,1) = hue image normalized to the range [0,1] by dividing
%                  all angle values by 2*pi.
%     HSI(:,:,2) = saturation image, in the range [0,1].
%     HSI(:,:,3) = intensity image, in the range [0,1].
%
%   All the conversion equations are explained in Chapter 7 of DIPUM3E.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Extract the individual component images.
rgb = im2double(rgb);
r = rgb(:,:,1);
g = rgb(:,:,2);
b = rgb(:,:,3);

% Implement the conversion equations.
num = 0.5*((r - g) + (r - b));
den = sqrt((r - g).^2 + (r - b).*(g - b));
theta = acos(num./(den + eps));

H = theta;
H(b > g) = 2*pi - H(b > g);
H = H/(2*pi);

num = min(min(r,g),b);
den = r + g + b;
den(den == 0) = eps;
S = 1 - 3.* num./den;

H(S == 0) = 0;

I = (r + g + b)/3;

% Combine all three results into an hsi image.
hsi = cat(3,H,S,I);

author: dipum

repo: dipum-toolbox

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


function results = OTB_HC_settings(seq, res_path, bSaveImage, parameters)

close all

s_frames = seq.s_frames;

% Feature specific parameters
hog_params.cell_size = 6;
hog_params.compressed_dim = 10;

cn_params.tablename = 'CNnorm';
cn_params.useForGray = false;
cn_params.cell_size = 4;
cn_params.compressed_dim = 3;

ic_params.tablename = 'intensityChannelNorm6';
ic_params.useForColor = false;
ic_params.cell_size = 4;
ic_params.compressed_dim = 3;

% Which features to include
params.t_features = {
    struct('getFeature',@get_fhog,'fparams',hog_params),...
    struct('getFeature',@get_table_feature, 'fparams',cn_params),...
    struct('getFeature',@get_table_feature, 'fparams',ic_params),...
};

% Global feature parameters1s
params.t_global.normalize_power = 2;    % Lp normalization with this p
params.t_global.normalize_size = true;  % Also normalize with respect to the spatial size of the feature
params.t_global.normalize_dim = true;   % Also normalize with respect to the dimensionality of the feature

% Image sample parameters
params.search_area_shape = 'square';    % The shape of the samples
params.search_area_scale = 4.0;         % The scaling of the target size to get the search area
params.min_image_sample_size = 150^2;   % Minimum area of image samples
params.max_image_sample_size = 200^2;   % Maximum area of image samples

% Detection parameters
params.refinement_iterations = 1;       % Number of iterations used to refine the resulting position in a frame
params.newton_iterations = 5;           % The number of Newton iterations used for optimizing the detection score
params.clamp_position = false;          % Clamp the target position to be inside the image

% Learning parameters
params.output_sigma_factor = 1/16;		% Label function sigma
params.learning_rate = 0.007;	 	 	% Learning rate
params.nSamples = 30;                   % Maximum number of stored training samples
params.sample_replace_strategy = 'lowest_prior';    % Which sample to replace when the memory is full
params.lt_size = 0;                     % The size of the long-term memory (where all samples have equal weight)
params.train_gap = 5;                   % The number of intermediate frames with no training (0 corresponds to training every frame)
params.skip_after_frame = 10;           % After which frame number the sparse update scheme should start (1 is directly)
params.use_detection_sample = true;     % Use the sample that was extracted at the detection stage also for learning

% Factorized convolution parameters
params.use_projection_matrix = true;    % Use projection matrix, i.e. use the factorized convolution formulation
params.update_projection_matrix = true; % Whether the projection matrix should be optimized or not
params.proj_init_method = 'pca';        % Method for initializing the projection matrix
params.projection_reg = 1e-6;	 	 	% Regularization paremeter of the projection matrix

% Generative sample space model parameters
params.use_sample_merge = true;                 % Use the generative sample space model to merge samples
params.sample_update_criteria = 'Merge';        % Strategy for updating the samples
params.weight_update_criteria = 'WeightedAdd';  % Strategy for updating the distance matrix
params.neglect_higher_frequency = false;        % Neglect hiigher frequency components in the distance comparison for speed

% Conjugate Gradient parameters
params.CG_iter = 5;                     % The number of Conjugate Gradient iterations in each update after the first frame
params.init_CG_iter = 10*20;            % The total number of Conjugate Gradient iterations used in the first frame
params.init_GN_iter = 10;               % The number of Gauss-Newton iterations used in the first frame (only if the projection matrix is updated)
params.CG_use_FR = false;               % Use the Fletcher-Reeves (true) or Polak-Ribiere (false) formula in the Conjugate Gradient
params.CG_standard_alpha = true;        % Use the standard formula for computing the step length in Conjugate Gradient
params.CG_forgetting_rate = 50;	 	 	% Forgetting rate of the last conjugate direction
params.precond_data_param = 0.75;	 	% Weight of the data term in the preconditioner
params.precond_reg_param = 0.25;	 	% Weight of the regularization term in the preconditioner
params.precond_proj_param = 60;	 	 	% Weight of the projection matrix part in the preconditioner

% Regularization window parameters
params.use_reg_window = true;           % Use spatial regularization or not
params.reg_window_min = 1e-4;			% The minimum value of the regularization window
params.reg_window_edge = 4e-3;          % The impact of the spatial regularization
params.reg_window_power = 2;            % The degree of the polynomial to use (e.g. 2 is a quadratic window)
params.reg_sparsity_threshold = 0.05;   % A relative threshold of which DFT coefficients that should be set to zero

% Interpolation parameters
params.interpolation_method = 'bicubic';    % The kind of interpolation kernel
params.interpolation_bicubic_a = -0.75;     % The parameter for the bicubic interpolation kernel
params.interpolation_centering = true;      % Center the kernel at the feature sample
params.interpolation_windowing = false;     % Do additional windowing on the Fourier coefficients of the kernel

% Scale parameters for the translation model
% Only used if: params.use_scale_filter = false
% params.number_of_scales = 7;            % Number of scales to run the detector
% params.scale_step = 1.01;               % The scale factor

% Scale filter parameters
% Only used if: params.use_scale_filter = true
params.use_scale_filter = true;        % Use the fDSST scale filter or not (for speed)
params.scale_sigma_factor = 1/16;       % Scale label function sigma
params.scale_learning_rate = 0.025;		% Scale filter learning rate
params.number_of_scales_filter = 17;    % Number of scales
params.number_of_interp_scales = 33;    % Number of interpolated scales
params.scale_model_factor = 1.0;        % Scaling of the scale model
params.scale_step_filter = 1.02;        % The scale factor for the scale filter
params.scale_model_max_area = 32*16;    % Maximume area for the scale sample patch
params.scale_feature = 'HOG4';          % Features for the scale filter (only HOG4 supported)
params.s_num_compressed_dim = 'MAX';    % Number of compressed feature dimensions in the scale filter
params.lambda = 1e-2;					% Scale filter regularization
params.do_poly_interp = true;           % Do 2nd order polynomial interpolation to obtain more accurate scale

% Other parameters
params.visualization = 0;               % Visualiza tracking and detection scores
params.debug = 0;                       % Do full debug visualization

% Initialize
params.init_sz = [seq.init_rect(1,4), seq.init_rect(1,3)];
params.init_pos = [seq.init_rect(1,2), seq.init_rect(1,1)] + (params.init_sz - 1)/2;
params.s_frames = s_frames;

% Run tracker
results = tracker(params);


text:
function [x_sigma, x_gam] = gaussian_para_esti(x)
% x = x(:);
gam = 0.2:0.001:10;
r_gam = (gamma(1./gam).*gamma(3./gam))./((gamma(2./gam)).^2);
x_mu = mean(x);
x_sigma_sq = mean((x - x_mu).^2);
x_sigma = sqrt(x_sigma_sq);
E_x = mean(abs(x - x_mu));
rho_x = x_sigma_sq/E_x^2;
[x_diff, x_ind] = min(abs(rho_x - r_gam));
x_gam = gam(x_ind);

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
function K = kmgraph(kern,d1,d2,ind1,ind2,kerparam),
% Marginalized Graph Kernel by Koji Tsuda
% The actual implementation uses the \delta version of the kernel (see publication)
% and a tuned C implementation by Alex Zien.
% Despite the paper we are using a normed K(X1,X2)/sqrt(K(X1,X1)*K(X2,X2))
% version.
% To use this kernel ensure that the function graphkerneldelta.cpp
% is compiled in the 'functions' folder.
% Furthermore note that d1.X must contain !cells! of the form :
% Let x,y,z be a graph, then
% struct(x)
%     adjacency: [NxN double]
%     vertices: [Nx1 double]
% and d=data({x,y,z}') is a valid data object
%
% gb 03-Mar-2004

X1=get_x(d1,ind1);
X2=get_x(d2,ind2);

K=zeros(length(X2),length(X1));
eps=1e-5;

lambda = kerparam;

if( size(ind1)==size(ind2))
    if( ind1==ind2)
        for i=ind1
            for j=ind2
                Ks(j,i)=graphkerneldelta(X1{i}.adjacency,X2{j}.adjacency,X1{i}.vertices,X2{j}.vertices,lambda,eps);

            end
        end

        for i=ind1
            for j=ind2
                K(j,i)=Ks(j,i)/sqrt(Ks(j,j)*Ks(i,i));
            end
        end
    end
else

    for i=ind1
        for j=ind2
            k12=graphkerneldelta(X1{i}.adjacency,X2{j}.adjacency,X1{i}.vertices,X2{j}.vertices,lambda,eps);
            k11=graphkerneldelta(X1{i}.adjacency,X1{i}.adjacency,X1{i}.vertices,X1{i}.vertices,lambda,eps);
            k22=graphkerneldelta(X2{j}.adjacency,X2{j}.adjacency,X2{j}.vertices,X2{j}.vertices,lambda,eps);
            K(j,i)=k12/sqrt(1e-5+k11*k22);
        end
    end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef MimMarkerMagenta < MimGuiPlugin
    % MimMarkerMagenta. Gui Plugin for switching transparency of zero
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimMarkerMagenta is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Magenta'
        SelectedText = 'Magenta'
        ToolTip = 'Place magenta markers'
        Category = 'Marker colour'
        Visibility = 'Dataset'
        Mode = 'Markers'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'single_marker.png'
        IconColour = GemMarkerPoint.DefaultColours{5}
        Location = 5
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.NewMarkerColour = 5;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && isequal(gui_app.GetCurrentModeName, MimModes.MarkerMode);
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.NewMarkerColour == 5;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function tests = test_spaces
  tests = functiontests(localfunctions);
end

function test_chol_update1(testCase)
    n = 5;
    A  = gallery('moler', n);
    L = sqrt(A(1,1));
    for c = 2:n
        % pick c elements from c-th column
        atom = A(1:c, c);
        L = spx.la.chol.chol_update(L, atom);
        estimate = L * L';
        verifyEqual(testCase, A(1:c, 1:c), estimate, 'AbsTol', 1e-6);
    end
end

function test_chol_update2(testCase)
    n = 10;
    D  = randn(n, n);
    A = D'*D;
    L = [];
    for c = 1:n
        % pick c elements from c-th column
        atom = A(1:c, c);
        L = spx.la.chol.chol_update(L, atom);
        estimate = L * L';
        verifyEqual(testCase, A(1:c, 1:c), estimate, 'AbsTol', 1e-6);
    end
end

author: indigits

repo: sparse-plex


text:
function element_neighbor = triangulation_neighbor_elements ( ...
  element_order, element_num, element_node )

%*****************************************************************************80
%
%% TRIANGULATION_NEIGHBOR_ELEMENTS determines element neighbors.
%
%  Discussion:
%
%    A triangulation of a set of nodes can be completely described by
%    the coordinates of the nodes, and the list of nodes that make up
%    each element.  However, in some cases, it is necessary to know
%    element adjacency information, that is, which element, if any,
%    is adjacent to a given element on a particular side.
%
%    This routine creates a data structure recording this information.
%
%    The primary amount of work occurs in sorting a list of 3 * ELEMENT_NUM
%    data items.
%
%    This routine was modified to use columns instead of rows.
%
%  Example:
%
%    The input information from ELEMENT_NODE:
%
%    Element    Nodes
%    --------   ---------------
%     1         3      4      1
%     2         3      1      2
%     3         3      2      8
%     4         2      1      5
%     5         8      2     13
%     6         8     13      9
%     7         3      8      9
%     8        13      2      5
%     9         9     13      7
%    10         7     13      5
%    11         6      7      5
%    12         9      7      6
%    13        10      9      6
%    14         6      5     12
%    15        11      6     12
%    16        10      6     11
%
%    The output information in ELEMENT_NEIGHBOR:
%
%    Element   Neighboring Elements
%    --------  ---------------------
%
%     1        -1     -1      2
%     2         1      4      3
%     3         2      5      7
%     4         2     -1      8
%     5         3      8      6
%     6         5      9      7
%     7         3      6     -1
%     8         5      4     10
%     9         6     10     12
%    10         9      8     11
%    11        12     10     14
%    12         9     11     13
%    13        -1     12     16
%    14        11     -1     15
%    15        16     14     -1
%    16        13     15     -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ELEMENT_ORDER, the order of the element.
%
%    Input, integer ELEMENT_NUM, the number of element.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), the nodes that
%    make up each element.
%
%    Output, integer ELEMENT_NEIGHBOR(3,ELEMENT_NUM), the three elements that are direct
%    neighbors of a given element.  ELEMENT_NEIGHBOR(1,I) is the index of the element
%    which touches side 1, defined by nodes 2 and 3, and so on.  ELEMENT_NEIGHBOR(1,I)
%    is negative if there is no neighbor on that side.  In this case, that
%    side of the element lies on the boundary of the triangulation.
%

%
%  Step 1.
%  From the list of nodes for element E, of the form: (I,J,K)
%  construct the three neighbor relations:
%
%    (I,J,3,E) or (J,I,3,E),
%    (J,K,1,E) or (K,J,1,E),
%    (K,I,2,E) or (I,K,2,E)
%
%  where we choose (I,J,3,E) if I < J, or else (J,I,3,E)
%
  col = zeros ( 4, 3 * element_num );

  for element = 1 : element_num

    i = element_node(1,element);
    j = element_node(2,element);
    k = element_node(3,element);

    if ( i < j )
      col(1:4,1+3*(element-1)) = [ i, j, 3, element ]';
    else
      col(1:4,1+3*(element-1)) = [ j, i, 3, element ]';
    end

    if ( j < k )
      col(1:4,2+3*(element-1)) = [ j, k, 1, element ]';
    else
      col(1:4,2+3*(element-1)) = [ k, j, 1, element ]';
    end

    if ( k < i )
      col(1:4,3+3*(element-1)) = [ k, i, 2, element ]';
    else
      col(1:4,3+3*(element-1)) = [ i, k, 2, element ]';
    end

  end
%
%  Step 2. Perform an ascending dictionary sort on the neighbor relations.
%  We only intend to sort on rows 1 and 2; the routine we call here
%  sorts on rows 1 through 4 but that won't hurt us.
%
%  What we need is to find cases where two elements share an edge.
%  Say they share an edge defined by the nodes I and J.  Then there are
%  two columns of COL that start out ( I, J, ?, ? ).  By sorting COL,
%  we make sure that these two columns occur consecutively.  That will
%  make it easy to notice that the elements are neighbors.
%
  col = i4col_sort_a ( 4, 3*element_num, col );
%
%  Step 3. Neighboring elements show up as consecutive columns with
%  identical first two entries.  Whenever you spot this happening,
%  make the appropriate entries in ELEMENT_NEIGHBOR.
%
  element_neighbor(1:3,1:element_num) = -1;

  icol = 1;

  while ( 1 )

    if ( 3 * element_num <= icol )
      break
    end

    if ( col(1,icol) ~= col(1,icol+1) || col(2,icol) ~= col(2,icol+1) )
      icol = icol + 1;
      continue
    end

    side1 = col(3,icol);
    element1 =  col(4,icol);
    side2 = col(3,icol+1);
    element2 =  col(4,icol+1);

    element_neighbor(side1,element1) = element2;
    element_neighbor(side2,element2) = element1;

    icol = icol + 2;

  end

  return
end


text:
function test_ft_rejectcomponent

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_rejectcomponent

fs = 500;
nchan = 32;
start_time = -1; % seconds
end_time = 2.5; % seconds
nsamples = (end_time - start_time) * fs + 1;

data = [];
data.time{1} = linspace(start_time, end_time, nsamples);
data.trial{1} = randn(nchan,nsamples);
data.label = cellstr(num2str((1:nchan).'));

cfg = [];
comp = ft_componentanalysis(cfg, data);

cfg = [];
cfg.component = [1 4 14];
dataout = ft_rejectcomponent(cfg,comp,data);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [source] = fixinside(source, target)

% FIXINSIDE ensures that the region of interest (which is indicated by the
% field "inside") is consistently defined for source structures and volume
% structures. Furthermore, it solves backward compatibility problems.
%
% Use as
%   [source] = fixinside(source, 'logical');
% or
%   [source] = fixinside(source, 'index');

% Copyright (C) 2006, 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
  target = 'logical';
end

if ~isfield(source, 'inside')
  if isfield(source, 'leadfield')
    % determine the positions outside the brain on basis of the empty cells
    source.inside = ~cellfun(@isempty, source.leadfield);
  elseif isfield(source, 'filter')
    % determine the positions outside the brain on basis of the empty cells
    source.inside = ~cellfun(@isempty, source.filter);
  elseif isfield(source, 'pos')
    % assume that all positions are inside the region of interest
    source.inside  = [1:size(source.pos,1)]';
    source.outside = [];
  elseif isfield(source, 'dim') && isfield(source, 'transform')
    % assume that all positions are inside the region of interest
    source.inside  = [1:prod(source.dim(1:3))]';
    source.outside = [];
  end
end

if ~isfield(source, 'inside')
  % nothing to do
  return;
end

% determine the format
if isa(source.inside, 'logical')
  current = 'logical';
elseif all(source.inside(:)==0 | source.inside(:)==1)
  source.inside = logical(source.inside);
  current = 'logical';
else
  current = 'indexed';
end

if strcmp(current, 'indexed') && strcmp(target, 'indexed')
  % nothing to do
elseif strcmp(current, 'logical') && strcmp(target, 'logical')
  % nothing to do
elseif strcmp(current, 'indexed') && strcmp(target, 'logical')
  % remove outside
  if isfield(source, 'outside')
    source = rmfield(source, 'outside');
  end
  % convert inside to a logical array
  if isfield(source, 'pos')
    tmp = false(size(source.pos,1), 1);
  elseif isfield(source, 'dim') && isfield(source, 'transform')
    tmp = false(source.dim(1:3));
  end
  tmp(source.inside) = true;
  source.inside = tmp;
elseif strcmp(current, 'logical') && strcmp(target, 'index')
  % convert to a vectors with indices
  tmp = source.inside;
  source.inside  = find( tmp(:));
  source.outside = find(~tmp(:));
else
  ft_error('incorrect specification of the insidestyle')
end

if strcmp(target, 'logical')
  if isfield(source, 'pos')
    % reshape it so that it matches the positions
    source.inside = reshape(source.inside, size(source.pos,1), 1);
  elseif isfield(source, 'dim') && isfield(source, 'transform')
    % reshape it so that it matches the volume dimensions
    source.inside = reshape(source.inside, source.dim(1:3));
  end
end


text:
function [coeff, score, latent, tsquared, explained, mu] = pca(x,varargin)
%PCA Principal Component Analysis (PCA) on raw data.
%   COEFF = PCA(X) returns the principal component coefficients for the N
%   by P data matrix X. Rows of X correspond to observations and columns to
%   variables. Each column of COEFF contains coefficients for one principal
%   component. The columns are in descending order in terms of component
%   variance (LATENT). PCA, by default, centers the data and uses the
%   singular value decomposition algorithm. For the non-default options,
%   use the name/value pair arguments.
%
%   [COEFF, SCORE] = PCA(X) returns the principal component score, which is
%   the representation of X in the principal component space. Rows of SCORE
%   correspond to observations, columns to components. The centered data
%   can be reconstructed by SCORE*COEFF'.
%
%   [COEFF, SCORE, LATENT] = PCA(X) returns the principal component
%   variances, i.e., the eigenvalues of the covariance matrix of X, in
%   LATENT.
%
%   [COEFF, SCORE, LATENT, TSQUARED] = PCA(X) returns Hotelling's T-squared
%   statistic for each observation in X. PCA uses all principal components
%   to compute the TSQUARED (computes in the full space) even when fewer
%   components are requested (see the 'NumComponents' option below). For
%   TSQUARED in the reduced space, use MAHAL(SCORE,SCORE).
%
%   [COEFF, SCORE, LATENT, TSQUARED, EXPLAINED] = PCA(X) returns a vector
%   containing the percentage of the total variance explained by each
%   principal component.
%
%   [COEFF, SCORE, LATENT, TSQUARED, EXPLAINED, MU] = PCA(X) returns the
%   estimated mean, MU, when 'Centered' is set to true; and all zeros when
%   set to false.
%
%   [...] = PCA(..., 'PARAM1',val1, 'PARAM2',val2, ...) specifies optional
%   parameter name/value pairs to control the computation and handling of
%   special data types. Parameters are:
%
%    'Algorithm' - Algorithm that PCA uses to perform the principal
%                  component analysis. Choices are:
%        'svd'   - Singular Value Decomposition of X (the default).
%        'eig'   - Eigenvalue Decomposition of the covariance matrix. It
%                  is faster than SVD when N is greater than P, but less
%                  accurate because the condition number of the covariance
%                  is the square of the condition number of X.
%        'als'   - Alternating Least Squares (ALS) algorithm which finds
%                  the best rank-K approximation by factoring a X into a
%                  N-by-K left factor matrix and a P-by-K right factor
%                  matrix, where K is the number of principal components.
%                  The factorization uses an iterative method starting with
%                  random initial values. ALS algorithm is designed to
%                  better handle missing values. It deals with missing
%                  values without listwise deletion (see {'Rows',
%                  'complete'}).
%
%     'Centered' - Indicator for centering the columns of X. Choices are:
%         true   - The default. PCA centers X by subtracting off column
%                  means before computing SVD or EIG. If X contains NaN
%                  missing values, NANMEAN is used to find the mean with
%                  any data available.
%         false  - PCA does not center the data. In this case, the original
%                  data X can be reconstructed by X = SCORE*COEFF'.
%
%     'Economy'  - Indicator for economy size output, when D the degrees of
%                  freedom is smaller than P. D, is equal to M-1, if data
%                  is centered and M otherwise. M is the number of rows
%                  without any NaNs if you use 'Rows', 'complete'; or the
%                  number of rows without any NaNs in the column pair that
%                  has the maximum number of rows without NaNs if you use
%                  'Rows', 'pairwise'. When D < P, SCORE(:,D+1:P) and
%                  LATENT(D+1:P) are necessarily zero, and the columns of
%                  COEFF(:,D+1:P) define directions that are orthogonal to
%                  X. Choices are:
%         true   - This is the default. PCA returns only the first D
%                  elements of LATENT and the corresponding columns of
%                  COEFF and SCORE. This can be significantly faster when P
%                  is much larger than D. NOTE: PCA always returns economy
%                  size outputs if 'als' algorithm is specifed.
%         false  - PCA returns all elements of LATENT. Columns of COEFF and
%                  SCORE corresponding to zero elements in LATENT are
%                  zeros.
%
%     'NumComponents' - The number of components desired, specified as a
%                  scalar integer K satisfying 0 < K <= P. When specified,
%                  PCA returns the first K columns of COEFF and SCORE.
%
%     'Rows'     - Action to take when the data matrix X contains NaN
%                  values. If 'Algorithm' option is set to 'als, this
%                  option is ignored as ALS algorithm deals with missing
%                  values without removing them. Choices are:
%         'complete' - The default action. Observations with NaN values
%                      are removed before calculation. Rows of NaNs are
%                      inserted back into SCORE at the corresponding
%                      location.
%         'pairwise' - If specified, PCA switches 'Algorithm' to 'eig'.
%                      This option only applies when 'eig' method is used.
%                      The (I,J) element of the covariance matrix is
%                      computed using rows with no NaN values in columns I
%                      or J of X. Please note that the resulting covariance
%                      matrix may not be positive definite. In that case,
%                      PCA terminates with an error message.
%         'all'      - X is expected to have no missing values. All data
%                      are used, and execution will be terminated if NaN is
%                      found.
%
%     'Weights'  - Observation weights, a vector of length N containing all
%                  positive elements.
%
%     'VariableWeights' - Variable weights. Choices are:
%          - a vector of length P containing all positive elements.
%          - the string 'variance'. The variable weights are the inverse of
%            sample variance. If 'Centered' is set true at the same time,
%            the data matrix X is centered and standardized. In this case,
%            PCA returns the principal components based on the correlation
%            matrix.
%
%   The following parameter name/value pairs specify additional options
%   when alternating least squares ('als') algorithm is used.
%
%      'Coeff0'  - Initial value for COEFF, a P-by-K matrix. The default is
%                  a random matrix.
%
%      'Score0'  - Initial value for SCORE, a N-by-K matrix. The default is
%                  a matrix of random values.
%
%      'Options' - An options structure as created by the STATSET function.
%                  PCA uses the following fields:
%          'Display' - Level of display output.  Choices are 'off' (the
%                      default), 'final', and 'iter'.
%          'MaxIter' - Maximum number of steps allowed. The default is
%                      1000. Unlike in optimization settings, reaching
%                      MaxIter is regarded as convergence.
%           'TolFun' - Positive number giving the termination tolerance for
%                      the cost function.  The default is 1e-6.
%             'TolX' - Positive number giving the convergence threshold
%                      for relative change in the elements of L and R. The
%                      default is 1e-6.
%
%
%   Example:
%       load hald;
%       [coeff, score, latent, tsquared, explained] = pca(ingredients);
%
%   See also PPCA, PCACOV, PCARES, BIPLOT, BARTTEST, CANONCORR, FACTORAN,
%   ROTATEFACTORS.

% References:
%   [1] Jolliffe, I.T. Principal Component Analysis, 2nd ed.,Springer,2002.
%   [2] Krzanowski, W.J., Principles of Multivariate Analysis, Oxford
%       University Press, 1988.
%   [3] Seber, G.A.F., Multivariate Observations, Wiley, 1984.
%   [4] Jackson, J.E., A User's Guide to Principal Components, Wiley, 1988.
%   [5] Sam Roweis, EM algorithms for PCA and SPCA, In Proceedings of the
%       1997 conference on Advances in neural information processing
%       systems 10 (NIPS '97), MIT Press, Cambridge, MA, USA, 626-632,1998.
%   [6] Alexander Ilin and Tapani Raiko. Practical Approaches to Principal
%       Component Analysis in the Presence of Missing Values. J. Mach.
%       Learn. Res. 11 (August 2010), 1957-2000, 2010.

%   Copyright 2012 The MathWorks, Inc.

[n, p] = size(x);
internal.stats.checkSupportedNumeric('X',x,false,false,true); % complex is accepted here

% Parse arguments and check if parameter/value pairs are valid
paramNames = {'Algorithm','Centered','Economy','NumComponents','Rows',...
    'Weights','VariableWeights','Coeff0','Score0','Options'};
defaults   = {'svd',       true,      true,    p,           'complete',...
    ones(1,n,'like',x) ,ones(1,p,'like',x),        [],      [], statset('pca')};

[vAlgorithm, vCentered, vEconomy, vNumComponents, vRows,vWeights,...
    vVariableWeights, c0, s0, opts, setFlag]...
    = internal.stats.parseArgs(paramNames, defaults, varargin{:});
% Validate String value for  Algorithm value
AlgorithmNames = {'svd','eig','als'};
vAlgorithm = internal.stats.getParamVal(vAlgorithm,AlgorithmNames,...
    '''Algorithm''');
% Validate boolean value for 'Centered' option
vCentered = internal.stats.parseOnOff(vCentered,'''Centered''');
% Validate boolean value for 'Economy' option
vEconomy = internal.stats.parseOnOff(vEconomy,'''Economy''');
% Validate the number of components option 'NumComponents'
if ~isempty(x) && ~internal.stats.isScalarInt(vNumComponents,1,p)
    error(message('stats:pca:WrongNumComponents',p));
end
% Validate value for 'Rows' option
RowsNames = {'complete','pairwise','all'};
vRows = internal.stats.getParamVal(vRows,RowsNames,'''Rows''');

switch vAlgorithm
    case 'svd'
        % Switch 'Algorithm' to 'eig' if 'Rows' set to 'pairwise'
        if strcmp(vRows,'pairwise')
            if setFlag.Algorithm
                warning(message('stats:pca:NoPairwiseSVD'));
            end
            vAlgorithm = 'eig';
        end
        % Switch algorithm to 'als' if user specify 'Coeff0' and 'Score0'.
        if setFlag.Coeff0 || setFlag.Score0
            vAlgorithm = 'als';
        end
    case 'als'
        % If 'als' is chosen, force PCA to use ALS and to ignore the
        % Rows' option
        if setFlag.Rows
            warning(message('stats:pca:NoALSRows'));
        end
end

% Validate Weights Vectors
if isvector(vWeights) && isequal(numel(vWeights),n)
    vWeights = reshape(vWeights,1,n); % make sure it is a row vector
else
    error(message('stats:pca:WrongObsWeights', n));
end

if internal.stats.isString(vVariableWeights)
    WeightsNames = {'variance'};
    internal.stats.getParamVal(vVariableWeights,WeightsNames,...
        '''VariableWeights''');
    vVariableWeights = 1./classreg.learning.internal.wnanvar(x,vWeights,1);
elseif isnumeric(vVariableWeights) && isvector(vVariableWeights)...
        && (isequal(numel(vVariableWeights), p))
    vVariableWeights = reshape(vVariableWeights,1,p);
else
    error(message('stats:pca:WrongVarWeights', p));
end

if any(vWeights <= 0) || any(vVariableWeights <= 0)
    error(message('stats:pca:NonPositiveWeights'));
end
% end of checking input arguments

% Handling special empty case
if isempty(x)
    pOrZero = ~vEconomy * p;
    coeff = zeros(p, pOrZero);
    coeff(1:p+1:end) = 1;
    score = zeros(n, pOrZero);
    latent = zeros(pOrZero, 1);
    tsquared = zeros(n, 1);
    explained = [];
    mu = [];
    return;
end

nanIdx = isnan(x);
numNaN = sum(nanIdx, 2); % number of NaNs in each row
wasNaN = any(numNaN,2); % Rows that contain NaN

% Handling special cases where X is all NaNs:
if all(nanIdx(:))
    coeff = NaN;
    score = NaN;
    latent = NaN;
    tsquared = NaN;
    explained = NaN;
    mu = NaN;
    return;
end
% Handling special scalar case;
if isscalar(x)
    coeff = 1;
    score = (~vCentered)*x;
    latent = (~vCentered)*x^2;
    tsquared = ~vCentered;
    explained = 100;
    mu = vCentered*x;
    return;
end

if strcmp(vRows,'all') && (~strcmp(vAlgorithm,'als'))
    if any(wasNaN)
        error(message('stats:pca:RowsAll'));
    else
        vRows = 'complete';
    end
end

if strcmp(vRows,'complete')
    % Degrees of freedom (DOF) is n-1 if centered and n if not centered,
    % where n is the numer of rows without any NaN element.
    DOF = max(0,n-vCentered-sum(wasNaN));
elseif strcmp(vRows,'pairwise')
    % DOF is the maximum number of element pairs without NaNs
    notNaN = double(~nanIdx);
    nanC = notNaN'*notNaN;
    nanC = nanC.*(~eye(p));
    DOF = max(nanC(:));
    DOF = DOF-vCentered;
else
    DOF = max(0,n-vCentered);
end

if vCentered
    % Weighted sample mean:
    mu = classreg.learning.internal.wnanmean(x, vWeights);
else
    mu = zeros(1,p,'like',x);
end

% Compute by EIG if no weights are given
switch vAlgorithm
case 'eig'
    % Center the data if 'Centered' is true.
    if vCentered
        x = bsxfun(@minus,x,mu);
    end

    % Use EIG to compute.
    [coeff, eigValues] = localEIG(x, vCentered, vRows, vWeights,...
        vVariableWeights);

    % When 'Economy' value is true, nothing corresponding to zero
    % eigenvalues should be returned.
    if (DOF<p)
        if vEconomy
            coeff(:, DOF+1:p) = [];
            eigValues(DOF+1:p, :) = [];
        else % make sure they are zeros.
            eigValues(DOF+1:p, :) = 0;
        end
    end

    % Check if eigvalues are all postive
    if any(eigValues<0)
        error(message('stats:pca:CovNotPositiveSemiDefinite'));
    end

    if nargout > 1
        score = x/coeff';
        latent = eigValues; % Output Eigenvalues
        if nargout > 3
            tsquared = localTSquared(score, latent, n, p);
        end
    end

case 'svd' % Use SVD to compute
    % Center the data if 'Centered' is true.
    if vCentered
        x = bsxfun(@minus,x,mu);
    end

    [U,sigma, coeff, wasNaN] = localSVD(x, n,...
        vEconomy, vWeights, vVariableWeights);
    if nargout > 1
        score =  bsxfun(@times,U,sigma');
        latent = sigma.^2./DOF;
        if nargout > 3
            tsquared = localTSquared(score,latent,DOF,p);
        end
        %Insert NaNs back
        if any(wasNaN)
            score = internal.stats.insertnan(wasNaN, score);
            if nargout >3
                tsquared = internal.stats.insertnan(wasNaN,tsquared);
            end
        end
    end

    if DOF < p
        % When 'Economy' value is true, nothing corresponding to zero
        % eigenvalues should be returned.
        if vEconomy
            coeff(:, DOF+1:end) = [];
            if nargout > 1
                score(:, DOF+1:end)=[];
                latent(DOF+1:end, :)=[];
            end
        elseif nargout > 1
        % otherwise, eigenvalues and corresponding outputs need to pad
        % zeros because svd(x,0) does not return columns of U corresponding
        % to components of (DOF+1):p.
            score(:, DOF+1:p) = 0;
            latent(DOF+1:p, 1) = 0;
        end
    end

case 'als' % Alternating Least Square Algorithm

    vNumComponents = min([vNumComponents,n-vCentered,p]);  % ALS always return economy sized outputs

    if isempty(s0);
        s0 = randn(n,vNumComponents,'like',x);
    elseif ~isequal(size(s0),[n,vNumComponents])|| any(isnan(s0(:)))
        error(message('stats:pca:BadInitialValues','Score0',n,vNumComponents));
    end
    if isempty(c0);
        c0 = randn(p,vNumComponents,'like',x);
    elseif ~isequal(size(c0),[p,vNumComponents]) || any(isnan(c0(:)))
        error(message('stats:pca:BadInitialValues','Coeff0',p,vNumComponents));
    end

    [score,coeff,mu,latent]=alsmf(x,vNumComponents,'L0',s0,'R0',c0,...
        'Weights',vWeights,'VariableWeights',vVariableWeights,...
        'Orthonormal',true,'Centered',vCentered,'Options',opts);

    if nargout > 3
        % T-squared values are in reduced space.
        tsquared = localTSquared(score, latent,n-vCentered, vNumComponents);
    end
end % end of switch vAlgorithm

% Calcuate the percentage of the total variance explained by each principal
% component.
if nargout > 4
    explained = 100*latent/sum(latent);
end

% Output only the first k principal components
if (vNumComponents<DOF)
    coeff(:, vNumComponents+1:end) = [];
    if nargout > 1
    score(:, vNumComponents+1:end) = [];
    end
end

% Enforce a sign convention on the coefficients -- the largest element in
% each column will have a positive sign.
[~,maxind] = max(abs(coeff), [], 1);
[d1, d2] = size(coeff);
colsign = sign(coeff(maxind + (0:d1:(d2-1)*d1)));
coeff = bsxfun(@times, coeff, colsign);
if nargout > 1
    score = bsxfun(@times, score, colsign); % scores = score
end

end % End of main function

%----------------Subfucntions--------------------------------------------

function [coeff, eigValues]=localEIG(x,vCentered, vRows,vWeights,...
        vVariableWeights)
% Compute by EIG. vRows are the options of handing NaN when compute
% covariance matrix

% Apply observation and variable weights
OmegaSqrt = sqrt(vWeights);
PhiSqrt = sqrt(vVariableWeights);
x = bsxfun(@times, x, OmegaSqrt');
x = bsxfun(@times, x, PhiSqrt);

xCov = ncnancov(x, vRows, vCentered);

[coeff, eigValueDiag] = eig(xCov);
[eigValues, idx] = sort(diag(eigValueDiag), 'descend');
coeff = coeff(:, idx);

coeff = bsxfun(@times, coeff,1./PhiSqrt');
end

function [U,sigma, coeff, wasNaN] = localSVD(x, n,...,
    vEconomy, vWeights, vVariableWeights)
% Compute by SVD. Weights are supplied by vWeights and vVariableWeights.

% Remove NaNs missing data and record location
[~,wasNaN,x] = internal.stats.removenan(x);
if n==1  % special case because internal.stats.removenan treats all vectors as columns
    wasNaN = wasNaN';
    x = x';
end

% Apply observation and variable weights
vWeights(wasNaN) = [];
OmegaSqrt = sqrt(vWeights);
PhiSqrt = sqrt(vVariableWeights);
x = bsxfun(@times, x, OmegaSqrt');
x = bsxfun(@times, x, PhiSqrt);

if vEconomy
    [U,sigma,coeff] = svd(x,'econ');
else
    [U,sigma, coeff] = svd(x, 0);
end

U = bsxfun(@times, U, 1./OmegaSqrt');
coeff = bsxfun(@times, coeff,1./PhiSqrt');

if n == 1     % sigma might have only 1 row
    sigma = sigma(1);
else
    sigma = diag(sigma);
end
end

function tsquared = localTSquared(score, latent,DOF, p)
% Subfunction to calulate the Hotelling's T-squared statistic. It is the
% sum of squares of the standardized scores, i.e., Mahalanobis distances.
% When X appears to have column rank < r, ignore components that are
% orthogonal to the data.

if isempty(score)
    tsquared = score;
    return;
end

r = min(DOF,p); % Max possible rank of x;
if DOF > 1
    q = sum(latent > max(DOF,p)*eps(latent(1)));
    if q < r
        warning(message('stats:pca:ColRankDefX', q));
    end
else
    q = 0;
end
standScores = bsxfun(@times, score(:,1:q), 1./sqrt(latent(1:q,:))');
tsquared = sum(standScores.^2, 2);
end

function c = ncnancov(x,Rows,centered)
%   C = NCNANCOV(X) returns X'*X/N, where N is the number of observations
%   after removing rows missing values.
%
%   C = NCNANCOV(...,'pairwise') computes C(I,J) using rows with no NaN
%   values in columns I or J.  The result may not be a positive definite
%   matrix. C = NCNANCOV(...,'complete') is the default, and it omits rows
%   with any NaN values, even if they are not in column I or J.
%
%   C = NCNANCOV(...,true), C is normalized by N-1 if data X is already
%   centered. The default is false.

if nargin <2
    Rows = 'complete';
end

d = 0;
if nargin>2
    d =  d + centered;
end

idxnan = isnan(x);

[n, p] = size(x);

if ~any(any(idxnan))
    c = x'*x/(n-d);
elseif strcmp(Rows,'complete')
    nanrows = any(idxnan,2);
    xNotNaN = x((~nanrows),:);
    denom = max(0, (size(xNotNaN,1)-d) );
    c = xNotNaN'*xNotNaN/denom;
elseif strcmp(Rows,'pairwise')
    c = zeros(p,class(x));
    for i = 1:p
        for j = 1:i
            NonNaNRows = ~any(idxnan(:,[i, j]),2);
            denom = max(0,(sum(NonNaNRows)-d));
            c(i,j) = x(NonNaNRows,i)'*x(NonNaNRows,j)/denom;
        end
    end
    c = c + tril(c,-1)';
end
end

author: Daikenan

repo: ASRCF

save_path: github-repos/MATLAB/Daikenan-ASRCF


text:
function u = tapas_datagen_categorical
% This function generates categorical input data for the hgf_categorical model
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2015 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% First set of outcomes
u = mnrnd(1, [0.8, 0.1, 0.1], 64);

% Second set of outcomes
u = [u; mnrnd(1, [1/3, 1/3, 1/3], 64)];

% Third set of outcomes
u = [u; mnrnd(1, [0.1, 0.1, 0.8], 64)];

% Add next set of outcomes (...or don't)

% Turn u into a single column of natural numbers indicating outcome category
u = sum(u*diag([1 2 3]),2);

author: translationalneuromodeling

repo: tapas


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
%COFICOSTFUNC Collaborative filtering cost function
%   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
%   num_features, lambda) returns the cost and gradient for the
%   collaborative filtering problem.
%

% Unfold the U and W matrices from params
X = reshape(params(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(params(num_movies*num_features+1:end), ...
                num_users, num_features);

% You need to return the following values correctly
J = 0;
X_grad = zeros(size(X));
Theta_grad = zeros(size(Theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost function and gradient for collaborative
%               filtering. Concretely, you should first implement the cost
%               function (without regularization) and make sure it is
%               matches our costs. After that, you should implement the
%               gradient and use the checkCostFunction routine to check
%               that the gradient is correct. Finally, you should implement
%               regularization.
%
% Notes: X - num_movies  x num_features matrix of movie features
%        Theta - num_users  x num_features matrix of user features
%        Y - num_movies x num_users matrix of user ratings of movies
%        R - num_movies x num_users matrix, where R(i, j) = 1 if the
%            i-th movie was rated by the j-th user
%
% You should set the following variables correctly:
%
%        X_grad - num_movies x num_features matrix, containing the
%                 partial derivatives w.r.t. to each element of X
%        Theta_grad - num_users x num_features matrix, containing the
%                     partial derivatives w.r.t. to each element of Theta
%

% X.shape=[nm, n]; Theta.shape=[nu, n]
% Y.shape=[nm, nu];
J = 1.0 / 2 * sum(sum(R.*(X*Theta'-Y).^2));
J = J + lambda/2*sum(sum(Theta.^2)) + lambda/2*sum(sum(X.^2));

% X_grad.shape=[nm, n]; Theta.grad.shape=[nu, n]
X_grad = R.*(X*Theta'-Y)*Theta;
X_grad = X_grad + lambda*X;
Theta_grad = (R.*(X*Theta'-Y))'*X;
Theta_grad = Theta_grad + lambda*Theta;

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

grad = [X_grad(:); Theta_grad(:)];

end

author: imLogM


text:
% Uniform on a Moebius'string
function [data] = mobius_rotate(N)
 r = 3; R = 10;
%r = 5; R = 15;
%np = 100; nr = 20;
%N = 2000;

xi = rand(N,1);
eta = rand(N,1);
rho = (1-2*xi)*r;
theta = pi*eta;
%for k=1:10
theta = pi*eta+rho./R.*cos(theta);
%end
x = (R+rho.*sin(theta)).*cos(2*theta);
y = (R+rho.*sin(theta)).*sin(2*theta);
z = rho.*cos(theta);
x_r=x+8;
y_r=y+8;
z_r=z;

% figure;
% plot3(x,y,z,'b.')
% hold on
% plot3(x_r,y_r,z_r,'r.')
% axis('equal');
X1=[x';y';z'];
X2=[x_r';y_r';z_r'];
X_P=[X1,X2];
data=X_P;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [fo, ip] = func_excludeoutlier_cor( fi, g, g_s );
%======================================================================
%
% Version 1.0
%
% This program excludes the points which below g_s of g
%
% Input
%   fi : input x data
%
% Output
%   fo : excluded x data
%   ip : excluded array element number in xi and yi
%
% Example:
%   [fo, ip] = func_excludeoutlier_cor( fi, g, g_s );
%
%
%======================================================================
% Terms:
%
%       Distributed under the terms of the terms of the BSD License
%
% Copyright:
%
%       Nobuhito Mori, Kyoto University
%
%========================================================================
%
% Update:
%       1.00    2005/01/12 Nobuhito Mori
%
%========================================================================

n = max(size(fi));
fo = [];
ip = find(g<g_s);
%fo = fi(find(g>=g_s));
fo =fi;
fo(ip) = NaN;

author: Sable

repo: mcbench-benchmarks


text:
function drawSurfPatch(u, v, z, varargin)
%DRAWSURFPATCH Draw a 3D surface patch, with 2 parametrized surfaces
%
%   usage:
%   drawSurfPatch(u, v, zuv)
%   where u, v, and zuv are three matrices the same size, u and
%   corresponding to each parameter, and zuv being equal to a function of u
%   and v.
%
%   drawSurfPatch(u, v, zuv, p0)
%   If p0 is specified, two lines with u(p0(1)) and v(p0(2)) are drawn on
%   the surface, and corresponding tangent are also shown.
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 24/05/2005.
%

%   HISTORY
%   2005-06-08 add doc.
%   2007-01-04 remove unused variables and change function name
%   2010-03-08 code cleanup, use drawPolyline3d

% prepare figure
hold on;

% draw the surface interior
surf(u, v, z, 'FaceColor', 'g', 'EdgeColor', 'none');

% draw the surface boundaries
drawPolyline3d(u(1,:), v(1,:), z(1,:))
drawPolyline3d(u(end,:), v(end,:), z(end,:))
drawPolyline3d(u(:,end), v(:,end), z(:,end))
drawPolyline3d(u(:,1), v(:,1), z(:,1))

% eventually draw two perpendicular lines on the surface
if ~isempty(varargin)
    pos = varargin{1};
    drawPolyline3d(u(pos(1),:), v(pos(1),:), z(pos(1),:));
    drawPolyline3d(u(:,pos(2)), v(:,pos(2)), z(:,pos(2)));
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%AKSHAY GORE
% email id:akshaygore@live.com
%MOB NO 9464894829
%Chandigarh university
%EDGE DETECTION  BASED ON
% 'average'   averaging filter
%     'disk'      circular averaging filter
%     'gaussian'  Gaussian lowpass filter
%     'laplacian' filter approximating the 2-D Laplacian operator
%     'log'       Laplacian of Gaussian filter
%     'motion'    motion filter
%     'prewitt'   Prewitt horizontal edge-emphasizing filter
%     'sobel'     Sobel horizontal edge-emphasizing filter
%     'unsharp'   unsharp contrast enhancement filter
function EDGE_DETECTION_NEW_2013(userImage)
% Clean up.
clc; % Clear the command window.
clear all;
close all; % Close all figures (except those of imtool.)
workspace;
fontSize = 18;

if nargin == 0
    % No image passed in on the command line.
    % Read in one of the standard MATLAB demo images
    % as our original gray scale image and display it.
   promptMessage = sprintf('Which image do you want to use');
    button = questdlg(promptMessage, 'Select Image', 'Coins', 'Cameraman','Coins');
     if strcmp(button, 'Coins')
        grayImage = double(imread('coins.png')); % Cast to double.
    else
        grayImage = double(imread('cameraman.tif')); % Cast to double.
    end
else
    % Use the image array passed in on the command line.
    grayImage = double(userImage); % Cast to double.
end
% Start timing.
%startTime = tic;

figure(1),subplot(2, 3, 1);
imshow(grayImage, []);
title('Original Image', 'FontSize', fontSize);
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.

% Blur the image with a 5 by 5 averaging (box filter) window.
blurredImage = conv2(grayImage, ones(5,5)/25);
figure(1),subplot(2, 3, 2);
imshow(blurredImage, []);
title('Blurred Image', 'FontSize', fontSize);
% Perform a variance filter.
% Output image is the variance of the input image in a 3 by 3 sliding window.
VarianceFilterFunction = @(x) var(x(:));
varianceImage = nlfilter(grayImage, [3 3], VarianceFilterFunction);
% An alternate way of doing the variance filter is on the next line:
% varianceImage = reshape(std(im2col(originalImage,[3 3],'sliding')),
%size(Original Image-2);
figure(1),subplot(2, 3, 3);
imshow(varianceImage, [])
title('Variance Image', 'FontSize', fontSize);

% Compute the square root of the variance image to get the standard deviation.
standardDeviationImage = sqrt(varianceImage);
figure(1),subplot(2, 3, 4);
imshow(standardDeviationImage, [])
title('Standard Deviation Image', 'FontSize', fontSize);

% Compute the standard deviation filter with MATLAB's built-in stdfilt() filter.
standardDeviationImage2 = stdfilt(grayImage);
figure(1),subplot(2, 3, 5);
imshow(standardDeviationImage2, [])
title('Built-in stdfilt() filter', 'FontSize', fontSize);

% Perform Sobel filter on image in both direction
% h = fspecial('sobel') returns a 3-by-3 filter h (shown below) thatemphasizes horizontal edges
% using the smoothing effect by approximating a vertical gradient.
% If you need to emphasize vertical edges, transpose the filter h'.
% [ 1 2 1
% 0 0 0
% -1 -2 -1 ]
verticalSobel = fspecial('sobel')';
sobelImage = imfilter(blurredImage, verticalSobel);
figure(2),subplot(3, 3, 1);
imshow(sobelImage, [])
title('Sobel edge filter', 'FontSize', fontSize);
%Perform UNSHARP filter
verticalunsharp = fspecial('unsharp');
unsharpImage = imfilter(blurredImage, verticalunsharp);
figure(2),subplot(3,3,2);
imshow(unsharpImage, [])
title('Unsharp image', 'FontSize', fontSize);
set(gcf, 'Position', get(0,'Screensize'));
%Perform LAPLACIAN filter
verticallaplacian = fspecial('laplacian');
laplacianImage = imfilter(blurredImage, verticallaplacian);
figure(2),subplot(3,3,3);
imshow(laplacianImage, [])
title('Laplacian image', 'FontSize', fontSize);
%Perform LOG filter
verticallog = fspecial('log');
logImage = imfilter(blurredImage, verticallog);
figure(2),subplot(3,3,4);
imshow(logImage, [])
title('Log image', 'FontSize', fontSize);
%Perform MOTION FILTER ON IMAGE
verticalmotion = fspecial('motion');
motionImage = imfilter(blurredImage, verticalmotion);
figure(2),subplot(3,3,5);
imshow(motionImage, [])
title('Motion image', 'FontSize', fontSize);
%PERFORM PREWITT FILTER ON IMAGE
verticalprewitt = fspecial('prewitt');
prewittImage = imfilter(blurredImage, verticalprewitt);
figure(2),subplot(3,3,6);
imshow(prewittImage, [])
title('Prewitt image', 'FontSize', fontSize);
% PERFORM AVERAGE FILTER ON  IMAGE
verticalaverage = fspecial('average');
averageImage = imfilter(blurredImage, verticalaverage);
figure(2),subplot(3,3,7);
imshow(averageImage, [])
title('Average image', 'FontSize', fontSize);
%PERFORM DISK  FILTER ON IMAGE
verticaldisk = fspecial('disk');
diskImage = imfilter(blurredImage, verticaldisk);
figure(2),subplot(3,3,8);
imshow(diskImage, [])
title('Disk image', 'FontSize', fontSize);

%CHANDIGARH UNIVERSITY
end


text:
function pass = test_interp1(pref)

%% Linear interpolation:
% Test a scalar function:
seedRNG(6178);
x = (0:10).';
y = sin(x);
f = chebfun.interp1(x, y, 'linear');
tol = 10*eps;
pass(1) = norm(feval(f, x) - y) < tol;
pass(2) = numel(f.funs) == 10;
pass(3) = length(f) == 20;

% Test an array-valued function:
x = (0:10)';
y = [sin(x), cos(x)];
f = chebfun.interp1(x, y, 'linear');
tol = 10*eps;
pass(4) = norm(feval(f, x) - y) < tol;
pass(5) = numel(f.funs) == 10;
pass(6) = length(f) == 20;

% Test a different domain:
dom = [.01, 10.01];
x = (0:10)';
y = [sin(x), cos(x)];
f = chebfun.interp1(x, y, 'linear', dom);
err = norm(feval(f, x(2:end-1)) - y(2:end-1,:));
tol = 10*eps;
pass(7) = all(f.domain == [dom(1), 1:10, dom(2)]);
pass(8) = err < tol;
pass(9) = numel(f.funs) == 11;
pass(10) = length(f) == 22;

%% Polynomial interpolation:
% Test an scalar-valued function:
x = (0:10)';
y = sin(x);
f = chebfun.interp1(x, y);
tol = 100*eps;
pass(11) = norm(feval(f, x) - y) < tol;
pass(12) = numel(f.funs) == 1;
pass(13) = length(f) == length(x);

% Test an array-valued function:
x = (0:10)';
y = [sin(x), cos(x)];
f = chebfun.interp1(x, y);
tol = 100*eps;
pass(14) = norm(feval(f, x) - y) < tol;
pass(15) = numel(f.funs) == 1;
pass(16) = length(f) == length(x);

% Test a different domain:
dom = [.01, 10.01];
f = chebfun.interp1(x, y, dom);
tol = 100*eps;
pass(17) = norm(feval(f, x) - y) < tol;
pass(18) = numel(f.funs) == 1;
pass(19) = length(f) == length(x);
pass(20) = all(f.domain == dom);

%% Check that things work for row vector input.
x = (0:10).';
y = sin(x);
f1 = chebfun.interp1(x.', y, 'linear');
f2 = chebfun.interp1(x, y.', 'linear');
f3 = chebfun.interp1(x.', y.', 'linear');
pass(21) = (norm(feval(f1, x) - y) < 10*vscale(f1)*eps) && ...
           (norm(feval(f2, x) - y) < 10*vscale(f2)*eps) && ...
           (norm(feval(f3, x) - y) < 10*vscale(f3)*eps);

% Test random points:
x = rand(11,1);
y = sin(x);
f = chebfun.interp1(x, y, 'linear');
tol = 10*eps;
pass(22) = norm(feval(f, x) - y) < tol;
pass(23) = numel(f.funs) == 10;
pass(24) = length(f) == 20;

% Test a chebfun
x = chebfun('x', [0, 1]);
y = exp(x);
x = rand(11, 1);
f = chebfun.interp1(x, y, 'poly', [0, 1]);
tol = 100*eps;
pass(25) = norm(f(x)-y(x), inf) < tol;

%% Test trigonometric interpolation:
x = (0:10)';
y = cos(pi*x);
f = chebfun.interp1(x, y, 'trig');
tol = 100*eps;
pass(26) = norm(feval(f, x(2:end-1)) - y(2:end-1)) < tol && ...
    norm(feval(f, x(1))-(y(1)+y(end))/2) < tol ;
pass(27) = numel(f.funs) == 1;
pass(28) = length(f) == length(x)-1;

%%
% Test an array-valued function:
x = (0:10)';
y = [exp(sin(x)), cos(pi*x)];
f = chebfun.interp1(x, y, 'trig', [0, 11]);
tol = 100*eps;
pass(29) = norm(feval(f, x) - y) < tol;
pass(30) = numel(f.funs) == 1;
pass(31) = length(f) == length(x);
%%
% Test a different domain:
dom = [-.01, 10.01];
f = chebfun.interp1(x, y, dom);
tol = 1e3*eps;
pass(32) = norm(feval(f, x) - y) < tol;
pass(33) = numel(f.funs) == 1;
pass(34) = length(f) == length(x);
pass(35) = all(f.domain == dom);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [tp, fp, tn, fn, rawcounts] = segmentationMetrics(prob, gt, threshPR, threshIU)
	% prob in HxWxnumClass

	numClass = size(prob, 3);

	% Calculate category wise precision and recall
	if(length(threshPR) == 0)
		tp = zeros(0, numClass);
		fp = tp; tn = tp; fn = tp;
	else
		for i = 1:numClass,
			gtI = double(gt == i);
			gtI(gt == 0) = NaN;
			probI = prob(:,:,i);
			[tp(:,i), fp(:,i), tn(:,i), fn(:,i)] = getPR(gtI(:), probI(:), threshPR);
		end
	end

	rawcounts = getRawCounts(gt, prob, threshIU);
end

function rawcounts = getRawCounts(gt, prob, thresh)
	% Calculate Intersection over Union for different thresholds of projection...
	num = size(prob, 3) + 1;
	rawcounts = zeros([num, num, length(thresh)]);
	[conf, ind] = max(prob, [], 3);
	locs = gt >= 0;
	for i = 1:length(thresh),
		resim = ind .* (conf >= thresh(i));
		sumim = 1+gt+resim*num;
		hs = histc(sumim(locs), 1:num*num);
    	rawcounts(:,:,i) = rawcounts(:,:,i) + reshape(hs(:), size(rawcounts(:,:,i)));
	end
end

function a = reverse(a)
	a = a(end:-1:1);
end

function [tp fp tn fn] = getPR(gt, out, thresh)
	assert(issorted(thresh));
	out = out+1;
	thresh = thresh+1;
	tp = histc(out(:).*gt(:), thresh);
	tp = reverse(cumsum(reverse(tp)));
	fp = histc(out(:).*(1-gt(:)), thresh);
	fp = reverse(cumsum(reverse(fp)));
	fn = nnz(gt(:) == 1) - tp;
	tn = nnz(gt(:) == 0) - fp;
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
% Extended Kalman Filter.
%
% Files
%   innovation      - Innovation of an observation.
%   correctBlockEkf - Correct in block-defined EKF.
%   predictBlockEkf - Covariance predict in block-defined EKF.

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function GA
% Genetic Algorithm(real coding)
% By: Javad Ivakpour
% E-mail: javad7@gmail.com
% May 2006
% Goal: find maximum of function that introduced in fun00.m file in current
% directory and can be plotted in plot00
% this file is also include the random serach for comparision
tic
clc
figure(1)
clf
clear all
format long

%------------------------        parameters        ------------------------
% befor using this function you must specified your function in fun00.m
% file in current directory and then set the parameters
var=2;            % Number of variables (this item must be equal to the
                  %   number of variables that is used in the function in
                  %   fun00.m file)
n=100;            % Number of population

m0=20;            % Number of generations that max value remains constant
                  %   (use for termination criteria)
nmutationG=20;                  %number of mutation children(Gaussian)
nmutationR=20;                  %number of mutation children(random)
nelit=2;                        %number of elitism children
valuemin=ones(1,var)*-5*pi;     % min possible value of variables
valuemax=ones(1,var)*5*pi;      % max possible value of variables

%-------------------------------------------------------------------------
nmutation=nmutationG+nmutationR;
sigma=(valuemax-valuemin)/10;    %Parameter that related to Gaussian
                                 %   function and used in mutation step
max1=zeros(nelit,var);
parent=zeros(n,var);
cu=[valuemin(1) valuemax(1) valuemin(2) valuemax(2)];
for l=1:var
    p(:,l)=valuemin(l)+rand(n,1).*(valuemax(l)-valuemin(l));
end
initial=p;
m=m0;
maxvalue=ones(m,1)*-1e10;
maxvalue(m)=-1e5;
g=0;
meanvalue(m)=0;
%-------------   ****    termination criteria   ****-------------
while abs(maxvalue(m)-maxvalue(m-(m0-1)))>0.001*maxvalue(m) &...
        (abs(maxvalue(m))>1e-10 & abs(maxvalue(m-(m0-1)))>1e-10)...
        & m<10000 & abs(maxvalue(m)-meanvalue(m))>1e-5 | m<20
    sigma=sigma./(1.05);% reducing the sigma value
    %  ------  **** % reducing the number of mutation()random   **** ----
    g=g+1;
    if g>10 & nmutationR>0
        g=0;
        nmutationR=nmutationR-1;
        nmutation=nmutationG+nmutationR;
    end

    %-------------   ****    function evaluation   ****-------------
    for i=1:n
        y(i)=fun00(p(i,:));
    end
    s=sort(y);
    maxvalue1(1:nelit)=s(n:-1:n-nelit+1);
    if nelit==0
        maxvalue1(1)=s(n);
        for i=1:n
            if y(i)==maxvalue1(1)
                max1(1,:)=p(i,:);
            end
        end
    end
    for k=1:nelit
        for i=1:n
            if y(i)==maxvalue1(k)
                max1(k,:)=p(i,:);
            end
        end
    end
    if var==2
        figure(1)
        subplot(2,2,1)
        hold off
        plot00(cu)
        hold on
        plot3(p(:,1),p(:,2),y,'ro')
        plot3(max1(1,1),max1(1,2),maxvalue1(1),'bh')
        title({' Genetic Algorithm '...
            ,'Performance of GA ( o : each individual)'},'color','b')

    end
    y=y-min(y)*1.02;
    sumd=y./sum(y);
    meanvalue=y./(sum(y)/n);

    %-------------   ****   Selection: Rolette wheel   ****-------------
    for l=1:n
        sel=rand;
        sumds=0;
        j=1;
        while sumds<sel
            sumds=sumds+sumd(j);
            j=j+1;
        end
        parent(l,:)=p(j-1,:);
    end
    p=zeros(n,var);

    %-------------   ****    regeneration   ****-------------
    for l=1:var

        %-------------   ****    cross-over   ****-------------
        for j=1:ceil((n-nmutation-nelit)/2)
            t=rand*1.5-0.25;
            p(2*j-1,l)=t*parent(2*j-1,l)+(1-t)*parent(2*j,l);
            p(2*j,l)=t*parent(2*j,l)+(1-t)*parent(2*j-1,l);
        end

        %-------------   ****    elitism   ****-------------
        for k=1:nelit
            p((n-nmutation-k+1),l)=max1(k,l);
        end

        %-------------   ****    mutation   ****-------------
        for i=n-nmutation+1:n-nmutationR
            phi=1-2*rand;
            z=erfinv(phi)*(2^0.5);
            p(i,l)=z*sigma(l)+parent(i,l);

        end
        for i=n-nmutationR+1:n
            p(i,1:var)=valuemin(1:var)+rand(1,var).*(valuemax(1:var)-...
                valuemin(1:var));
        end
        for i=1:n
            for l=1:var
                if p(i,l)<valuemin(l)
                    p(i,l)=valuemin(l);
                elseif p(i,l)>valuemax(l)
                    p(i,l)=valuemax(l);
                end
            end
        end
    end
    p;
    m=m+1;
    max1;
    maxvalue(m)=maxvalue1(1);
    maxvalue00(m-m0)=maxvalue1(1);
    mean00(m-m0)=sum(s)/n;
    meanvalue(m)=mean00(m-m0);
    figure(1)
    if var==2
        subplot(2,2,2)
    end
    hold off
    plot(maxvalue00,'b')
    hold on
    plot(mean00,'r')
    hold on
    title({'Performance of GA',...
        'best value GA:blue, best value RS:black, mean value GA:red',''}...
        ,'color','b')
    xlabel('number of generations')
    ylabel('value')

    %-------------   ****    Random search   ****-------------
    %-------------   ****   for comparision  ****-------------
    p00=zeros(n,var);
    for l=1:var
        p00(:,l)=valuemin(l)+rand(n,1).*(valuemax(l)-valuemin(l));
    end
    for i=1:n
        y(i)=fun00(p00(i,:));
    end
    s=sort(y);
    maxvalueRAND(m-m0)=s(n);
    if m>(m0+1)
        if maxvalueRAND(m-m0)<maxvalueRAND(m-(m0+1))
            maxvalueRAND(m-m0)=maxvalueRAND(m-(m0+1));
        else
            for i=1:n
                if y(i)==maxvalueRAND(m-m0)
                    maxRand=p00(i,:);
                end
            end
        end
    else
        for i=1:n
            if y(i)==maxvalueRAND(m-m0)
                maxRand=p00(i,:);
            end
        end
    end
    plot(maxvalueRAND,'k')
    if var==2
        figure(1)
        subplot(2,2,3)
        plot00(cu)
        hold on
        plot3(maxRand(1,1),maxRand(1,2),maxvalueRAND(m-m0),'k*')
        plot3(max1(1,1),max1(1,2),maxvalue00(m-m0),'bo')
        title({'Best solution found by GA(: o) & RS(:*)'...
            'in each generation ',''},'color','b')
    end
    pause(0.001)
end

clc
disp('     Genetic Algorithm(real coding)   ')
disp('          By: Javad Ivakpour          ')
disp('       E-mail: javad7@gmail.com       ')
disp('**************************************')
num_of_fun_evaluation=n*m
max_point_GA=max1(1,:)
maxvalue_GA=maxvalue00(m-m0)
max_point_RS=maxRand
maxvalue_RS=maxvalueRAND(m-m0)
if var==2
    figure(1)
    subplot(2,2,4)
    hold off
    plot3(max1(1,1),max1(1,2),maxvalue1,'o')
    hold on
    plot00(cu)
    hold on
    plot3(maxRand(1,1),maxRand(1,2),maxvalueRAND(m-m0),'*')
    title({'Best solution ';'GA: o & RS: *'},'color','b')
end
figure(2)
title('Performance of GA(best value)','color','b')
xlabel('number of generations')
ylabel('max value of best solution')
hold on
plot(maxvalue00)
hold on
toc


text:
function [ SV ] = moveFV(V,F,S)
% MOVEFV  Move a scalar field defined on faces to vertices by averaging
%
% [ SV ] = moveFV(V,F,S)
%
% Input:
%   V,F  mesh
%   S  scalar field defined on faces, Fx1
%
% Output:
%   SV  scalar field defined on vertices

SV = zeros(size(V,1),size(S,2));
COUNT = zeros(size(V,1),1);

for i=1:size(F,1)
    SV(F(i,:)',:) = SV(F(i,:)',:) + repmat(S(i,:),size(F,2),1);
    COUNT(F(i,:)') = COUNT(F(i,:)') + 1;
end

SV = SV ./ repmat(COUNT, 1, size(S,2));

end

author: alecjacobson

repo: gptoolbox


text:
function [obXM,obYM,lp] = mrSpreadInplanes(numofanats,obXM,obYM,lp)
%NAME:   [obXM,obYM,lp] = mrSpreadInplanes(numofanats,obXM,obYM,lp)
%AUTHOR:  Poirson
%DATE:	  08.04.96
%PURPOSE: One of a set of routines that allows user
%         to set and select a set of oblique planes in saggital slice.
%	  The routines are mrTransInplanes.m, mrRotInplanes.m,
%	  mrClipInplanes.m, mrSelInPlane.m, mrSetupInplanes.m,
%         mrSpreadInplanes.m
%HISTORY: Started with mrGetOblPlane from G. Boynton 4/6/96
%NOTES:
%BUGS:
%  Loses the cropping information when changing inplane spread.

INCDELTA = 0.25;

% One more pair of points for the perpendicular line
nPtPairs = numofanats + 1;

% See if the user has already been working on some inPlanes
if size(obXM,1) == 0
	disp('You must first create a candidate set of Inplanes');
	disp('Choose Set_Up_Inplanes first');
	return
end

% Use X and Y travel of perpendicular line
% to find incX and incY between lines.
incX = (obXM(nPtPairs,2)- obXM(nPtPairs,1))/(numofanats-1);
incY = (obYM(nPtPairs,2)- obYM(nPtPairs,1))/(numofanats-1);

% percent of change in X and Y direction;
fractionX = incX/(incX+incY);
fractionY = incY/(incX+incY);

disp('');
disp('--Setting Distance Between Inplanes--');
disp('  Left Button=Decrease, Middle Button=Increase, Right Button=Quit');
button = 0;
while(button~=3)
 [tempx,tempy,button]=mrGinput(1,'cross');
 if (button~=3)
  if (button == 1)
	delta = -1*INCDELTA;
  end
  if (button == 2)
	delta = INCDELTA;
  end

  deltaX = fractionX * delta;
  deltaY = fractionY * delta;

  incX = incX + deltaX;
  incY = incY + deltaY;

  for i=1:numofanats
	obXM(i,:) = obXM(1,:) + (i-1)*incX;
	obYM(i,:) = obYM(1,:) + (i-1)*incY;
  end

  % the perpendicular line, (starting point is the same)
  obXM(nPtPairs,2) = obXM(nPtPairs,1)+(numofanats-1)*incX;
  obYM(nPtPairs,2) = obYM(nPtPairs,1)+(numofanats-1)*incY;

  % draw the new lines on the screen
  for i=1:nPtPairs
     delete(lp(i));
     lp(i)=line(obXM(i,:),obYM(i,:),'Color','w');
  end
 end
end

disp('--Done Setting Distance Between Inplanes--');
disp('');

return;

author: vistalab

repo: vistasoft


text:
function fem_basis_q4_display ( prefix )

%*****************************************************************************80
%
%% FEM_BASIS_Q4_DISPLAY displays a finite element Q4 basis function.
%
%  Discussion:
%
%    This program reads a data file defining a set of nodes, and a
%    data file defining the triangulation of those nodes using 3 node triangles
%    (or 6 node triangles, as long as the vertices are listed first).
%
%    The program then asks the user interactively to select one of the
%    nodes.  It computes the basis function associated with that node
%    and displays it over the entire mesh.  Of course, the basis function
%    will only be nonzero over a small number of the elements, but it
%    is instructive to see the entire mesh.
%
%    The display is initially "flat", but by using the manipulator
%    on the graphics menu, the user can easily get some dramatic images
%    of the basis function.
%
%  Usage:
%
%    fem_basis_q4_display ( 'prefix' )
%
%    where 'prefix' is the common prefix for the FEM files:
%
%    * 'prefix'_nodes.txt,    the node coordinates.
%    * 'prefix'_elements.txt, the nodes that make up each element;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string PREFIX, the common file prefix.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM_BASIS_Q4_DISPLAY:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Display basis functions associated with \n' );
  fprintf ( 1, '  a finite element grid of linear quadrilaterals ("Q4").\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The user specifies a basis function by node index, and\n' );
  fprintf ( 1, '  the program displays a surface plot of that basis function.\n' );
  fprintf ( 1, '  (Use the 3D ROTATE option to see the full picture!\n' );
%
%  Get the prefix if missing.
%
  if ( nargin < 1 )
    prefix = input ( '  Enter the common file prefix:  ' );
  end
%
%  Construct the file names.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  element_filename = strcat ( prefix, '_elements.txt' );
%
%  Read the nodes.
%
  [ dim_num, node_num ] = r8mat_header_read ( node_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension = %d\n', dim_num );
  fprintf ( 1, '  Number of nodes = %d\n', node_num );

  if ( dim_num ~= 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FEM_BASIS_Q4_DISPLAY - Fatal error!\n' );
    fprintf ( 1, '  The spatial dimension is not 2.\n' );
    error ( 'FEM_BASIS_Q4_DISPLAY - Fatal error!' );
  end

  node_xy = r8mat_data_read ( node_filename, dim_num, node_num );
%
%  Read the elements.
%
  [ element_order, element_num ] = i4mat_header_read ( element_filename );

  if ( element_order ~= 4 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FEM_BASIS_Q4_DISPLAY\n' );
    fprintf ( 1, '  This program requires that the elements have order 4.\n' );
    fprintf ( 1, '  Your elements have order %d\n', element_order );
    error ( 'FEM_BASIS_Q4_DISPLAY - Fatal error!' );
  end

  element_node = i4mat_data_read ( element_filename, element_order, ...
    element_num );
%
%  Set up the graph.
%
  x_min = min ( node_xy(1,:) );
  x_max = max ( node_xy(1,:) );
  y_min = min ( node_xy(2,:) );
  y_max = max ( node_xy(2,:) );
  z_min = 0.0;
  z_max = 1.0;

  node_min = min ( min ( element_node ) );
  node_max = max ( max ( element_node ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Every basis function is associated with a node.\n' );
  fprintf ( 1, '  To chooose a basis function, you specify a node.\n' );
  fprintf ( 1, '  Nodes range in value from %d to %d.\n', node_min, node_max );

  while ( 1 )

    fprintf ( 1, '\n' );
    prompt = sprintf ( 'Enter a node between %d and %d: ', node_min, node_max );
    node_index = input ( prompt );

    if ( node_index < node_min | node_max < node_index )
      break
    end
%
%  Clear the graphics page.
%
    clf

    fprintf ( 1, '\n' );

    for element = 1 : element_num

      local = 0;
      for j = 1 : element_order
        if ( element_node(j,element) == node_index )
          fprintf ( 1, '  Node %d occurs as local node %d in element %d.\n', ...
            node_index, j, element );
          local = j;
        end
      end

      z(1:2,1:2) = 0.0;

      x(1,1) = node_xy(1,element_node(1,element));
      y(1,1) = node_xy(2,element_node(1,element));
      x(2,1) = node_xy(1,element_node(2,element));
      y(2,1) = node_xy(2,element_node(2,element));
      x(1,2) = node_xy(1,element_node(4,element));
      y(1,2) = node_xy(2,element_node(4,element));
      x(2,2) = node_xy(1,element_node(3,element));
      y(2,2) = node_xy(2,element_node(3,element));

      if ( local ~= 0 )
        if ( local == 1 )
          z(1,1) = 1.0;
        elseif ( local == 2 )
          z(2,1) = 1.0;
        elseif ( local == 3 )
          z(2,2) = 1.0;
        else
          z(1,2) = 1.0;
        end
      end

      caxis ( [ -0.4, 1.2 ] )
      surface ( x, y, z, 'FaceColor', 'interp' )

    end

    axis ( [ x_min, x_max, y_min, y_max, z_min, z_max ] );
    axis equal

    xlabel ( '--X axis--' );
    ylabel ( '--Y axis--' );
    zlabel ( '--Z axis--' );

    title_string = sprintf ( 'Q4 basis function for node %d', node_index );
    title ( title_string );

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

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
    column_num = -1;
    return;
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
    row_num = -1;
    return;
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function table = i4mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% I4MAT_DATA_READ reads data from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns in the data.
%
%    Output, integer TABLE(M,N), the point coordinates.
%
  table = [];
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %d' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'I4MAT_DATA_READ - Error!' );
    return;
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
      fprintf ( 1, '  End of input while reading data.\n' );
      error ( 'I4MAT_DATA_READ - Error!' );
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = i4mat_header_read ( input_filename )

%*****************************************************************************80
%
%% I4MAT_HEADER_READ reads the header from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = [];
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
    return;
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LENGTH, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dirSynch( root1, root2, showOnly, flag, ignDate )
% Synchronize two directory trees (or show differences between them).
%
% If a file or directory 'name' is found in both tree1 and tree2:
%  1) if 'name' is a file in both the pair is considered the same if they
%     have identical size and identical datestamp (or if ignDate=1).
%  2) if 'name' is a directory in both the dirs are searched recursively.
%  3) if 'name' is a dir in root1 and a file in root2 (or vice-versa)
%     synchronization cannot proceed (an error is thrown).
% If 'name' is found only in root1 or root2 it's a difference between them.
%
% The parameter flag controls how synchronization occurs:
%  flag==0: neither tree1 nor tree2 has preference (newer file is kept)
%  flag==1: tree2 is altered to reflect tree1 (tree1 is unchanged)
%  flag==2: tree1 is altered to reflect tree2 (tree2 is unchanged)
% Run with showOnly=1 and different values of flag to see its effect.
%
% By default showOnly==1. If showOnly, displays a list of actions that need
% to be performed in order to synchronize the two directory trees, but does
% not actually perform the actions. It is highly recommended to run
% dirSynch first with showOnly=1 before running it with showOnly=0.
%
% USAGE
%  dirSynch( root1, root2, [showOnly], [flag], [ignDate] )
%
% INPUTS
%  root1    - root directory of tree1
%  root2    - root directory of tree2
%  showOnly - [1] show but do NOT perform actions
%  flag     - [0] 0: synchronize; 1: set root2=root1; 2: set root1==root2
%  ignDate  - [0] if true considers two files same even if have diff dates
%
% OUTPUTS
%  dirSynch( 'c:\toolbox', 'c:\toolbox-old', 1 )
%
% EXAMPLE
%
% See also
%
% Piotr's Image&Video Toolbox      Version 2.10
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if(nargin<3 || isempty(showOnly)), showOnly=1; end;
if(nargin<4 || isempty(flag)), flag=0; end;
if(nargin<5 || isempty(ignDate)), ignDate=0; end;
% get differences between root1/root2 and loop over them
D = dirDiff( root1, root2, ignDate );
roots={root1,root2}; ticId = ticStatus;
for i=1:length(D)
  % get action
  if( flag==1 )
    if( D(i).in1 ), act=1; src1=1; else act=0; src1=2; end
  elseif( flag==2 )
    if( D(i).in2 ), act=1; src1=2; else act=0; src1=1; end
  else
    act=1;
    if(D(i).in1 && D(i).in2)
      if( D(i).new1 ), src1=1; else src1=2; end
    else
      if( D(i).in1 ), src1=1; else src1=2; end
    end
  end
  src2=mod(src1,2)+1;
  % perform action
  if( act==1 )
    if( showOnly )
      disp(['COPY ' int2str(src1) '->' int2str(src2) ': ' D(i).name]);
    else
      copyfile( [roots{src1} D(i).name],  [roots{src2} D(i).name], 'f' );
    end;
  else
    if( showOnly )
      disp(['DEL in  ' int2str(src1) ': ' D(i).name]);
    else
      fName = [roots{src1} D(i).name];
      if(D(i).isdir), rmdir(fName,'s'); else delete(fName); end
    end
  end
  if(~showOnly), tocStatus( ticId, i/length(D) ); end;
end
end

function D = dirDiff( root1, root2, ignDate )
% get differences from root1 to root2
D1 = dirDiff1( root1, root2, ignDate, '/' );
% get differences from root2 to root1
D2 = dirDiff1( root2, root1, ignDate, '/' );
% remove duplicates (arbitrarily from D2)
D2=D2(~([D2.in1] & [D2.in2]));
% swap 1 and 2 in D2
for i=1:length(D2),
  D2(i).in1=0; D2(i).in2=1; D2(i).new1=~D2(i).new1;
end
% merge
D = [D1 D2];
end

function D = dirDiff1( root1, root2, ignDate, subdir )
if(root1(end)~='/'), root1(end+1)='/'; end
if(root2(end)~='/'), root2(end+1)='/'; end
if(subdir(end)~='/'), subdir(end+1)='/'; end
fs1=dir([root1 subdir]); fs2=dir([root2 subdir]);
D=struct('name',0,'isdir',0,'in1',0,'in2',0,'new1',0);
D=repmat(D,[1 length(fs1)]); n=0; names2={fs2.name}; Dsub=[];
for i1=1:length( fs1 )
  name=fs1(i1).name; isdir=fs1(i1).isdir;
  if( any(strcmp(name,{'.','..'})) ), continue; end;
  i2 = find(strcmp(name,names2));
  if(~isempty(i2) && isdir)
    % cannot handle this condition
    if(~fs2(i2).isdir), disp([root1 subdir name]); assert(false); end;
    % recurse and record possible differences
    Dsub=[Dsub dirDiff1(root1,root2,ignDate,[subdir name])]; %#ok<AGROW>
  elseif( ~isempty(i2) && fs1(i1).bytes==fs2(i2).bytes && ...
      (ignDate || fs1(i1).datenum==fs2(i2).datenum))
    % nothing to do - files are same
    continue;
  else
    % record differences
    n=n+1;
    D(n).name=[subdir name]; D(n).isdir=isdir;
    D(n).in1=1; D(n).in2=~isempty(i2);
    D(n).new1 = ~D(n).in2 || (fs1(i1).datenum>fs2(i2).datenum);
  end
end
D = [D(1:n) Dsub];
end

author: zouchuhang

repo: LayoutNet