text:
% Plot Simulation and On-board Model Output

sim_INS_Out = get(logsout, 'INS_Out');
ref_INS_Out = get(logsout, 'INS_Out_Ref');

sim_FMS_Out = get(logsout, 'FMS_Out');
ref_FMS_Out = get(logsout, 'FMS_Out_Ref');

sim_Control_Out = get(logsout, 'Control_Out');
ref_Control_Out = get(logsout, 'Control_Out_Ref');

%% Euler
figure;
ax1 = subplot(3,1,1);
plot(sim_INS_Out.Values.phi, 'Color', 'r');
hold on;
plot(ref_INS_Out.Values.phi, 'Color', 'b');
grid on;
legend('\phi', '\phi_r')

ax2 = subplot(3,1,2);
plot(sim_INS_Out.Values.theta, 'Color', 'r');
hold on;
plot(ref_INS_Out.Values.theta, 'Color', 'b');
grid on;
legend('\theta', '\theta_r')

ax3 = subplot(3,1,3);
plot(sim_INS_Out.Values.psi, 'Color', 'r');
hold on;
plot(ref_INS_Out.Values.psi, 'Color', 'b');
grid on;
legend('\psi', '\psi_r')

linkaxes([ax1,ax2,ax3], 'x');

%% Velocity
figure;
ax1 = subplot(3,1,1);
plot(sim_INS_Out.Values.vn, 'Color', 'r');
hold on;
plot(ref_INS_Out.Values.vn, 'Color', 'b');
grid on;
legend('vn', 'vn_r')

ax2 = subplot(3,1,2);
plot(sim_INS_Out.Values.ve, 'Color', 'r');
hold on;
plot(ref_INS_Out.Values.ve, 'Color', 'b');
grid on;
legend('ve', 've_r')

ax3 = subplot(3,1,3);
plot(sim_INS_Out.Values.vd, 'Color', 'r');
hold on;
plot(ref_INS_Out.Values.vd, 'Color', 'b');
grid on;
legend('vd', 'vd_r')

linkaxes([ax1,ax2,ax3], 'x');

%% Position
figure;
ax1 = subplot(3,1,1);
plot(sim_INS_Out.Values.x_R, 'Color', 'r');
hold on;
plot(ref_INS_Out.Values.x_R, 'Color', 'b');
grid on;
legend('x', 'x_r')

ax2 = subplot(3,1,2);
plot(sim_INS_Out.Values.y_R, 'Color', 'r');
hold on;
plot(ref_INS_Out.Values.y_R, 'Color', 'b');
grid on;
legend('y', 'y_r')

ax3 = subplot(3,1,3);
plot(sim_INS_Out.Values.h_R, 'Color', 'r');
hold on;
plot(ref_INS_Out.Values.h_R, 'Color', 'b');
grid on;
legend('h', 'h_r')

linkaxes([ax1,ax2,ax3], 'x');

%% FMS Output
figure;
ax1 = subplot(3,1,1);
plot(sim_FMS_Out.Values.u_cmd, 'Color', 'r');
hold on;
plot(ref_FMS_Out.Values.u_cmd, 'Color', 'b');
grid on;
legend('u\_cmd', 'u\_cmd_r')
ax2 = subplot(3,1,2);
plot(sim_FMS_Out.Values.v_cmd, 'Color', 'r');
hold on;
plot(ref_FMS_Out.Values.v_cmd, 'Color', 'b');
grid on;
legend('v\_cmd', 'v\_cmd_r')
ax3 = subplot(3,1,3);
plot(sim_FMS_Out.Values.w_cmd, 'Color', 'r');
hold on;
plot(ref_FMS_Out.Values.w_cmd, 'Color', 'b');
grid on;
legend('w\_cmd', 'w\_cmd_r')
linkaxes([ax1,ax2,ax3], 'x');

figure;
ax1 = subplot(3,1,1);
plot(sim_FMS_Out.Values.phi_cmd, 'Color', 'r');
hold on;
plot(ref_FMS_Out.Values.phi_cmd, 'Color', 'b');
grid on;
legend('\phi\_cmd', '\phi\_cmd_r')
ax2 = subplot(3,1,2);
plot(sim_FMS_Out.Values.theta_cmd, 'Color', 'r');
hold on;
plot(ref_FMS_Out.Values.theta_cmd, 'Color', 'b');
grid on;
legend('\theta\_cmd', '\theta\_cmd_r')
ax3 = subplot(3,1,3);
plot(sim_FMS_Out.Values.psi_rate_cmd, 'Color', 'r');
hold on;
plot(ref_FMS_Out.Values.psi_rate_cmd, 'Color', 'b');
grid on;
legend('\psi\_rate\_cmd', '\psi\_rate\_cmd_r')
linkaxes([ax1,ax2,ax3], 'x');

%% Control Output
figure;

ax1 = subplot(4,1,1);
plot(sim_Control_Out.Values.actuator_cmd.Time, sim_Control_Out.Values.actuator_cmd.Data(:,1), 'Color', 'r');
hold on;
plot(ref_Control_Out.Values.actuator_cmd.Time, ref_Control_Out.Values.actuator_cmd.Data(:,1), 'Color', 'b');
grid on;
legend('motor1', 'motor1_r')

ax2 = subplot(4,1,2);
plot(sim_Control_Out.Values.actuator_cmd.Time, sim_Control_Out.Values.actuator_cmd.Data(:,2), 'Color', 'r');
hold on;
plot(ref_Control_Out.Values.actuator_cmd.Time, ref_Control_Out.Values.actuator_cmd.Data(:,2), 'Color', 'b');
grid on;
legend('motor2', 'motor2_r')

ax3 = subplot(4,1,3);
plot(sim_Control_Out.Values.actuator_cmd.Time, sim_Control_Out.Values.actuator_cmd.Data(:,3), 'Color', 'r');
hold on;
plot(ref_Control_Out.Values.actuator_cmd.Time, ref_Control_Out.Values.actuator_cmd.Data(:,3), 'Color', 'b');
grid on;
legend('motor3', 'motor3_r')

ax4 = subplot(4,1,4);
plot(sim_Control_Out.Values.actuator_cmd.Time, sim_Control_Out.Values.actuator_cmd.Data(:,4), 'Color', 'r');
hold on;
plot(ref_Control_Out.Values.actuator_cmd.Time, ref_Control_Out.Values.actuator_cmd.Data(:,4), 'Color', 'b');
grid on;
legend('motor4', 'motor4_r')

linkaxes([ax1,ax2,ax3,ax4], 'x');

author: Firmament-Autopilot

repo: FMT-Model


text:
%SerialLink.jacob_dot Derivative of Jacobian
%
% JDQ = R.jacob_dot(Q, QD) is the product (6x1) of the derivative of the
% Jacobian (in the world frame) and the joint rates.
%
% Notes::
% - This term appears in the formulation for operational space control XDD = J(Q)QDD + JDOT(Q)QD
% - Written as per the reference and not very efficient.
%
% References::
% - Fundamentals of Robotics Mechanical Systems (2nd ed)
%   J. Angleles, Springer 2003.
% - A unified approach for motion and force control of robot manipulators: The operational space formulation
%  O Khatib, IEEE Journal on Robotics and Automation, 1987.
%
% See also SerialLink.jacob0, diff2tr, tr2diff.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function Jdot = jacob_dot(robot, q, qd)

	n = robot.n;
    links = robot.links;

    % Using the notation of Angeles:
    %   [Q,a] ~ [R,t] the per link transformation
    %   P ~ R   the cumulative rotation t2r(Tj) in world frame
    %   e       the last column of P, the local frame z axis in world coordinates
    %   w       angular velocity in base frame
    %   ed      deriv of e
    %   r       is distance from final frame
    %   rd      deriv of r
    %   ud      ??

    for i=1:n
        T = links(i).A(q(i));
        Q{i} = t2r(T);
        a{i} = transl(T)';
    end

    P{1} = Q{1};
    e{1} = [0 0 1]';
    for i=2:n
        P{i} = P{i-1}*Q{i};
        e{i} = P{i}(:,3);
    end

    % step 1
    w{1} = qd(1)*e{1};
    for i=1:(n-1)
        w{i+1} = qd(i+1)*[0 0 1]' + Q{i}'*w{i};
    end

    % step 2
    ed{1} = [0 0 0]';
    for i=2:n
        ed{i} = cross(w{i}, e{i});
    end

    % step 3
    rd{n} = cross( w{n}, a{n});
    for i=(n-1):-1:1
        rd{i} = cross(w{i}, a{i}) + Q{i}*rd{i+1};
    end

    r{n} = a{n};
    for i=(n-1):-1:1
        r{i} = a{i} + Q{i}*r{i+1};
    end

    ud{1} = cross(e{1}, rd{1});
    for i=2:n
        ud{i} = cross(ed{i}, r{i}) + cross(e{i}, rd{i});
    end

    % step 4
    %  swap ud and ed
    v{n} = qd(n)*[ud{n}; ed{n}];
    for i=(n-1):-1:1
        Ui = blkdiag(Q{i}, Q{i});
        v{i} = qd(i)*[ud{i}; ed{i}] + Ui*v{i+1};
    end

    Jdot = v{1};

author: petercorke

repo: robotics-toolbox-matlab


text:
function [stop, options] = msOutFcn(optimValues, ~, totIterExpect, hWaitBar, iniMessage, ticID)
    elapsedTime = toc(ticID);
    avgTimePerRun = elapsedTime/optimValues.localrunindex;
    runsRemaining = totIterExpect - optimValues.localrunindex;
    estTimeRemaining = runsRemaining*avgTimePerRun;

    [timeRemainingStr] = getDurationStr(estTimeRemaining);
%     timeRemainingStr =  num2str(estTimeRemaining, '%8.3f');
%     timeRemainingStr = paddStr(timeRemainingStr, 7);

    if(isnan(estTimeRemaining) || ~isempty(strfind(timeRemainingStr,'NaN')))
        estTimeStr = '';
    else
        estTimeStr = ['Est. Time Left: ', timeRemainingStr];
    end

    msgStr = {[iniMessage, ' [',num2str(optimValues.localrunindex),'/',num2str(totIterExpect),']'],...
              estTimeStr};

    if(isa(hWaitBar, 'matlab.ui.dialog.ProgressDialog'))
        hWaitBar.Message = msgStr;
        hWaitBar.Value = optimValues.localrunindex/totIterExpect;
    else
        waitbar(optimValues.localrunindex/totIterExpect, hWaitBar, msgStr);
    end

    stop = false;
    if(nargout == 2)
        options=[];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [dx,dy]=   centered_gradient(input,dx,dy, nx,ny)
dx = zeros(1,round(nx*ny*+nx));
dy = zeros(1,round(nx*ny*+nx));
nx = round(nx);
ny = round(ny);
size(input);
for i = 1:ny-1
        for  j = 1:  nx-1
           k = round(i * nx + j);
           if(nx+k < length(input))
            dx(k) = 0.5*(input(k+1) - input(k-1));
            dy(k) = 0.5*(input(k+nx) - input(k-nx));
           end
        end
end
z = 1;
for  j = 2: nx-1

      dx(z) = 0.5*(input(z+1) - input(j-1));
        dy(z) = 0.5*(input(z+nx) - input(z));

    k = (ny - 1) * nx + z;
if(k < length(input))
        dx(k) = 0.5*(input(k+1) - input(k-1));
        dy(k) = 0.5*(input(k) - input(k-nx));
end
        z = z +1;
end
for  i = 1: ny-1
        p = (i * nx)+1;
        if(p+nx < length(input))
        dx(p) = 0.5*(input(p+1) - input(p));
        dy(p) = 0.5*(input(p+nx) - input(p-nx));
        end
        k = (i+1) * nx - 1;
if(k+nx < length(input))
        dx(k) = 0.5*(input(k) - input(k-1));
        dy(k) = 0.5*(input(k+nx) - input(k-nx));
end

    dx(1) = 0.5*(input(2) - input(1));
    dy(1) = 0.5*(input(nx) - input(1));

    dx(nx-1) = 0.5*(input(nx-1) - input(nx-2));
    dy(nx-1) = 0.5*(input(2*nx-1) - input(nx-1));

    dx((ny-1)*nx) = 0.5*(input((ny-1)*nx + 1) - input((ny-1)*nx));
    dy((ny-1)*nx) = 0.5*(input((ny-1)*nx) - input((ny-2)*nx));
    if(ny*nx-2 <length(input))
    dx(ny*nx-2) = 0.5*(input(ny*nx-2) - input(ny*nx-1-1-1));
    dy(ny*nx-2) = 0.5*(input(ny*nx-2) - input((ny-2)*nx-2));
    disp('at least once');
    end
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function b = onLine(point, line)
%ONLINE test if a point belongs to a line.
%
%   B = onLine(POINT, LINE)
%   with POINT being [xp yp], and LINE being [x0 y0 dx dy].
%   Returns 1 if point lies on the line, 0 otherwise.
%
%   If POINT is an N*2 array of points, B is a N*1 array of booleans.
%
%   If LINE is a N*4 arrat of line, B is a 1*N array of booleans.
%
%   See also:
%   lines2d, points2d, onEdge, onRay, angle3Points

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2003-10-31
% Copyright 2003 INRA - TPV URPOI - BIA IMASTE

% deprecation warning
warning('geom2d:deprecated', ...
    '''onLine'' is deprecated, use ''isPointOnLine'' instead');

Nl = size(line, 1);
Np = size(point, 1);

x0 = repmat(line(:,1)', Np, 1);
y0 = repmat(line(:,2)', Np, 1);
dx = repmat(line(:,3)', Np, 1);
dy = repmat(line(:,4)', Np, 1);
xp = repmat(point(:,1), 1, Nl);
yp = repmat(point(:,2), 1, Nl);

% test if lines are colinear
b = abs((xp-x0).*dy-(yp-y0).*dx)./sqrt(dx.*dx+dy.*dy) < 1e-14;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function varargout = process_ssmooth( varargin )
% PROCESS_SSMOOTH: Spatial smoothing of the sources (DEPRECATED)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2010-2016

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Spatial smoothing (DEPRECATED)';
    sProcess.FileTag     = 'ssmooth';
    sProcess.Category    = 'Filter';
    sProcess.SubGroup    = 'Sources';
    sProcess.Index       = 0;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/VisualGroup';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'results', 'timefreq'};
    sProcess.OutputTypes = {'results', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Default values for some options
    sProcess.isSourceAbsolute = 1;
    sProcess.processDim  = 2;    % Process time by time
    % === GAUSSIAN PROP
    sProcess.options.label1.Comment = 'Gaussian filter properties:';
    sProcess.options.label1.Type    = 'label';
    % === FWHM (kernel size)
    sProcess.options.fwhm.Comment = '      FWHM (Full width at half maximum):  ';
    sProcess.options.fwhm.Type    = 'value';
    sProcess.options.fwhm.Value   = {10, 'mm', 0};
    % === METHOD
    sProcess.options.label2.Comment = '
Distance between vertices (v1,v2):';
    sProcess.options.label2.Type    = 'label';
    sProcess.options.method.Comment = {'Euclidean distance: norm(v1-v2)', ...
                                       'Path length: number of edges between v1 and v2', ...
                                       'Average: (euclidian distance + path length) / 2'};
    sProcess.options.method.Type    = 'radio';
    sProcess.options.method.Value   = 3;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    % Absolute values
    if isfield(sProcess.options, 'source_abs') && sProcess.options.source_abs.Value
        strAbs = ',abs';
    else
        strAbs = '';
    end
    % Method
    switch (sProcess.options.method.Value)
        case 1,    Method = 'euclidian';
        case 2,    Method = 'path';
        case 3,    Method = 'average';
        otherwise, error(['Unknown method: ' sProcess.options.method.Value]);
    end
    % Final comment
    Comment = sprintf('%s (%d%c%s)', sProcess.Comment, sProcess.options.fwhm.Value{1}, Method(1), strAbs);
end

%% ===== RUN =====
function sInput = Run(sProcess, sInput) %#ok<DEFNU>
    global GlobalData;
    % Get options
    FWHM = sProcess.options.fwhm.Value{1} / 1000;
    switch (sProcess.options.method.Value)
        case 1,    Method = 'euclidian';
        case 2,    Method = 'path';
        case 3,    Method = 'average';
        otherwise, error(['Unknown method: ' sProcess.options.method.Value]);
    end

    % ===== LOAD SURFACE =====
    % Load the surface filename from results file
    switch (file_gettype(sInput.FileName))
        case {'results', 'link'}
            FileMat = in_bst_results(sInput.FileName, 0, 'SurfaceFile', 'GridLoc', 'Atlas', 'nComponents', 'HeadModelType');
            nComponents = FileMat.nComponents;
        case 'timefreq'
            FileMat = in_bst_timefreq(sInput.FileName, 0, 'SurfaceFile', 'GridLoc', 'Atlas', 'HeadModelType');
            nComponents = 1;
        otherwise
            error('Unsupported file format.');
    end
    % Error: cannot smooth results on volume grids
    if ~strcmpi(FileMat.HeadModelType, 'surface') % ~isempty(FileMat.GridLoc)
        bst_report('Error', sProcess, sInput, 'Spatial smoothing is only supported for surface head models.');
        sInput = [];
        return;
    % Error: cannot smooth results that are already based on atlases
    elseif ~isempty(FileMat.Atlas)
        bst_report('Error', sProcess, sInput, 'Spatial smoothing is not supported for sources based on atlases.');
        sInput = [];
        return;
    % Error: only for constrained sources
    elseif (nComponents ~= 1)
        bst_report('Error', sProcess, sInput, ['This process is only available for source models with constrained orientations.' 10 'With unconstrained orientations, the source maps are usually already very smooth.']);
        sInput = [];
        return;
    end

    % ===== COMPUTE SMOOTHING OPERATOR =====
    % Get existing interpolation for this surface
    Signature = sprintf('ssmooth(%d,%s):%s', round(FWHM*1000), Method, FileMat.SurfaceFile);
    WInterp = [];
    if isfield(GlobalData, 'Interpolations') && ~isempty(GlobalData.Interpolations) && isfield(GlobalData.Interpolations, 'Signature')
        iInterp = find(cellfun(@(c)isequal(c,Signature), {GlobalData.Interpolations.Signature}), 1);
        if ~isempty(iInterp)
            WInterp = GlobalData.Interpolations(iInterp).WInterp;
        end
    end
    % Calculate new interpolation matrix
    if isempty(WInterp)
        % Load surface file
        SurfaceMat = in_tess_bst(FileMat.SurfaceFile);
        % Compute the smoothing operator
        WInterp = tess_smooth_sources(SurfaceMat.Vertices, SurfaceMat.Faces, SurfaceMat.VertConn, FWHM, Method);
        % Check for errors
        if isempty(WInterp)
            sInput = [];
            return;
        end
        % Save interpolation in memory for future calls
        sInterp = db_template('interpolation');
        sInterp.WInterp   = WInterp;
        sInterp.Signature = Signature;
        if isempty(GlobalData.Interpolations)
            GlobalData.Interpolations = sInterp;
        else
            GlobalData.Interpolations(end+1) = sInterp;
        end
    end

    % ===== APPLY TO THE DATA =====
    % Apply smoothing operator
    for iFreq = 1:size(sInput.A,3)
        sInput.A(:,:,iFreq) = WInterp * sInput.A(:,:,iFreq);
    end
    % Force the output comment
    sInput.CommentTag = [sProcess.FileTag, num2str(FWHM*1000), Method(1)];
    % Do not keep the Std field in the output
    if isfield(sInput, 'Std') && ~isempty(sInput.Std)
        sInput.Std = [];
    end
end

author: brainstorm-tools


text:
function K = exponential_correlation ( s, t )

%*****************************************************************************80
%
%% EXPONENTIAL_CORRELATION evaluates the exponential correlation function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real S(*), T(*), pairs of argument values.
%
%    Output, real K(*), the correlation function values
%
  K = exp ( - abs ( s - t ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p05_title ( )

%*****************************************************************************80
%
%% P05_TITLE prints a title for problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Problem 05\n' );
  fprintf ( 1, '  Name:       ST04\n' );
  fprintf ( 1, '              Stroud #4, page 26.\n' );
  fprintf ( 1, '  Region:     0 <= X(i) <= 1\n' );
  fprintf ( 1, '  Integrand:  F(X) = 1 / ( 1 + sum ( 2 * X(i) ) )\n' );

  return
end


text:
classdef P2Function < FeFunction

    properties (Access = public)
    end

    properties (Access = private)
        interpolation
        coord
        connec
    end

    methods (Access = public)

        function obj = P2Function(cParams)
            obj.init(cParams);
            obj.createInterpolation();
            obj.createDOFCoordConnec();
        end

        function fxV = evaluate(obj, xV)
            obj.interpolation.computeShapeDeriv(xV);
            shapes = obj.interpolation.shape;
            nNode  = size(shapes,1);
            nGaus  = size(shapes,2);
            nF     = size(obj.fValues,2);
            nElem  = size(obj.connec,1);
            fxV = zeros(nF,nGaus,nElem);
            for iGaus = 1:nGaus
                for iNode = 1:nNode
                    node = obj.connec(:,iNode);
                    Ni = shapes(iNode,iGaus);
                    fi = obj.fValues(node,:);
                    f(:,1,:) = Ni*fi';
                    fxV(:,iGaus,:) = fxV(:,iGaus,:) + f;
                end
            end

        end

        function N = computeShapeFunctions(obj, quad)
%             obj.mesh.computeInverseJacobian(quad,obj.interpolation);
            xV = quad.posgp;
            obj.interpolation.computeShapeDeriv(xV);
            N = obj.interpolation.shape;
        end

        function dNdx  = computeCartesianDerivatives(obj,quad)
            nElem = size(obj.connec,1);
            nNode = obj.interpolation.nnode;
            nDime = obj.interpolation.ndime;
            nGaus = quad.ngaus;
            invJ  = obj.mesh.computeInverseJacobian(quad,obj.interpolation);
            dShapeDx  = zeros(nDime,nNode,nElem,nGaus);
            for igaus = 1:nGaus
                dShapes = obj.interpolation.deriv(:,:,igaus);
                for jDime = 1:nDime
                    invJ_JI   = invJ(:,jDime,:,igaus);
                    dShape_KJ = dShapes(jDime,:);
                    dSDx_KI   = bsxfun(@times, invJ_JI,dShape_KJ);
                    dShapeDx(:,:,:,igaus) = dShapeDx(:,:,:,igaus) + dSDx_KI;
                end
            end
            dNdx = dShapeDx;
        end

        function gradFun = computeGradient(obj, quad)
            dNdx = obj.computeCartesianDerivatives(quad);
            nDimf = obj.ndimf;
            nDims = size(dNdx, 1); % derivX, derivY (mesh-related?)
            nNode = size(dNdx, 2);
            nElem = size(dNdx, 3);
            nGaus = size(dNdx, 4);

            grad = zeros(nDims,nDimf, nElem, nGaus);
            for iGaus = 1:nGaus
                dNdx_g = dNdx(:,:,:,iGaus);
                for iDims = 1:nDims
                    for iNode = 1:nNode
                        dNdx_i = squeeze(dNdx_g(iDims, iNode,:));
                        nodes = obj.connec(:,iNode);
                        f = obj.fValues(nodes,:);
                        p = (dNdx_i.*f)';
                        pp(1,:,:) = p;
                        grad(iDims,:,:,iGaus) = grad(iDims,:,:,iGaus) + pp;
                    end
                end
            end
            fVR = reshape(grad, [nDims*nDimf,nElem, nGaus]);
            s.fValues = permute(fVR, [1 3 2]);
%             s.ndimf      = nDimf;
            s.quadrature = quad;
            gradFun = FGaussDiscontinuousFunction(s);
        end

        function symGradFun = computeSymmetricGradient(obj,quad)
            grad = obj.computeGradient(quad);
            nDimf = obj.ndimf;
            nDims = size(grad.fValues, 1)/nDimf;
            nGaus = size(grad.fValues, 2);
            nElem = size(grad.fValues, 3);

            gradReshp = reshape(grad.fValues, [nDims,nDimf,nGaus,nElem]);
            gradT = permute(gradReshp, [2 1 3 4]);
            symGrad = 0.5*(gradReshp + gradT);

            s.fValues    = reshape(symGrad, [nDims*nDimf,nGaus,nElem]);
            s.quadrature = quad;
            symGradFun = FGaussDiscontinuousFunction(s);
        end

        function plot(obj, m) % 2D domains only
            s.mesh          = m;
            s.interpolation = obj.interpolation;
            c = ConnecCoordFromInterpAndMesh(s);
            c.compute();
            coord = c.coord;
            connec = obj.connec(:, [1 4 2 5 3 6]);
            x = coord(:,1);
            y = coord(:,2);
            figure()
            for idim = 1:obj.ndimf
                subplot(1,obj.ndimf,idim);
                z = obj.fValues(:,idim);
                a = trisurf(connec,x,y,z);
                view(0,90)
    %             colorbar
                shading interp
                a.EdgeColor = [0 0 0];
                title(['dim = ', num2str(idim)]);
            end
        end

        function dofConnec = computeDofConnectivity(obj)
            conne  = obj.connec;
            nDimf  = obj.ndimf;
            nNode  = size(conne, 2);
            nDofsE = nNode*nDimf;
            dofsElem  = zeros(nDofsE,size(conne,1));
            for iNode = 1:nNode
                for iUnkn = 1:nDimf
                    idofElem   = nDimf*(iNode - 1) + iUnkn;
                    globalNode = conne(:,iNode);
                    idofGlobal = nDimf*(globalNode - 1) + iUnkn;
                    dofsElem(idofElem,:) = idofGlobal;
                end
            end
            dofConnec = dofsElem;
        end

        function dof = getDofsFromCondition(obj, condition)
            nodes = condition(obj.coord);
            iNode = find(nodes==1);
            dofElem = repmat(1:obj.ndimf, [length(iNode) 1]);
            dofMat = obj.ndimf*(iNode - 1) + dofElem;
            dof = sort(dofMat(:));
        end

    end

    methods (Access = public, Static)

        function p1 = create(mesh, ndimf)
            s.fValues = zeros(mesh.nnodes, ndimf); % wrong
            s.mesh    = mesh;
            p1 = P2Function(s);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh = cParams.mesh;
            obj.fValues = cParams.fValues;
            obj.ndimf   = size(cParams.fValues,2);
        end

        function createInterpolation(obj)
            m.type = obj.mesh.type;
            obj.interpolation = Interpolation.create(m,'QUADRATIC');
        end

        function createDOFCoordConnec(obj)
            s.mesh          = obj.mesh;
            s.interpolation = obj.interpolation;
            c = ConnecCoordFromInterpAndMesh(s);
            c.compute();
            obj.coord  = c.coord;
            obj.connec = c.connec;
            nDimf = size(obj.fValues,2);
            if isequal(size(obj.mesh.coord,1), size(obj.fValues,1))
                obj.fValues = zeros(size(obj.coord,1),nDimf);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [model, B, elapse] = SH_learn(A, maxbits)
%   This is a wrapper function of Spectral Hashing learning.
%
%	Usage:
%	[model, B, elapse] = SH_learn(A, maxbits)
%
%	      A: Rows of vectors of data points. Each row is sample point
%   maxbits: Code length
%
%     model: Used for encoding a test sample point.
%	      B: The binary code of the input data A. Each row is sample point
%    elapse: The coding time (training time).
%
%
%
%   version 2.0 --Dec/2016
%   version 1.0 --Jan/2013
%
%   Written by  Yue Lin (linyue29@gmail.com)
%               Deng Cai (dengcai AT gmail DOT com)
%

tmp_T = tic;

model = SpectralHashing(A, maxbits);

B = SH_compress(A, model);

elapse = toc(tmp_T);

end

author: ZJULearning

repo: MatlabFunc


text:
function t_dmc_element(quiet)
%T_DMC_ELEMENT  Tests for MP.DMC_ELEMENT.

%   MATPOWER
%   Copyright (c) 2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

% define_constants;
if quiet
    verbose = 0;
else
    verbose = 1;
end

table_class = mp_table_class();

%%  1       2       3       4       5       6
%%  alpha   beta    uid     gamma   delta   epsilon
foo = [ ...
    11      5       1       0.1     9       0
    22      4       2       0.3     3       2
    33      3       3       0.5     7       0
    44      2       4       0.7     5       4
    55      1       5       0.9     1       6
];
baz = {
    'uno', 'one';
    'dos', 'two';
    'tres', 'three';
    'cuatro', 'four';
    'cinco', 'five';
};
d0 = struct('foo', foo, 'bar', struct('baz', {baz}));
ridx1 = 2;
ridx2 = [3;1];

r = [2;4;5];
var_names = {'uid', 'name', 'status', 'source_uid', ...
        'ids', 'twos', 'alpha', 'beta', 'gamma', 'delta', 'epsilon'};
ec = cell(5,1); [ec{:}] = deal('');
vals = {foo(:,3), baz(:,2), ones(5,1), ec, foo(:,3), 2*ones(5,1), ...
    foo(:,1), 10*foo(:,2), 2*foo(:,4), 5*foo(:,5), foo(:,6)};
tab0 = table_class(vals{:}, 'VariableNames', var_names);
tab1 = tab0(r,:); tab1.ids = [1:3]'; tab1.source_uid = r;

t_begin(24, quiet);

%% 1 - dme table rows correspond to original data rows
k = 'A';
t = sprintf('%s : dmce : ', k);
dmce = mp.dmce_test_widget();
t_ok(isa(dmce, 'mp.dmc_element'), [t 'class']);
t_is(dmce.use_r, 0, 12, [t 'dmce.use_r == 0']);

t = sprintf('%s : dme : ', k);
dme = mp.dme_test_widget();
t_ok(isa(dme, 'mp.dm_element'), [t 'class']);
t_ok(isempty(dme.tab), [t 'dme.tab empty']);

t = sprintf('%s : import(dme, d) : ', k);
dme = dmce.import(mp.dme_test_widget(), d0);
t_ok(~isempty(dme.tab), [t 'dme.tab not empty']);
t_ok(isequal(dme.tab, tab0), [t 'dme.tab']);

t = sprintf('%s : export(dme, d) : ', k);
d = struct();
d = dmce.export(dme, d);
ed = d0; ed.bar.baz(:,1) = ed.bar.baz(:,2); [ed.bar.baz{:,2}] = deal([]);
t_ok(isequal(d, ed), [t 'd']);

t = sprintf('%s : import(dme, d, var_names) : ', k);
d = d0; d.foo = d0.foo(:, end:-1:1);    %% reverse columns
tab = tab0; tmp = tab.alpha; tab.alpha = tab.epsilon; tab.epsilon = tmp;
dme = dmce.import(mp.dme_test_widget(), d0);
dme = dmce.import(dme, d, {'alpha', 'epsilon'});
t_ok(isequal(dme.tab, tab), [t 'dme.tab']);

t = sprintf('%s : export(dme, d, var_names) : ', k);
tmp = dme.tab.beta; dme.tab.beta = dme.tab.delta; dme.tab.delta = tmp;
ed = d0; ed.foo(:, [2 5]) = ed.foo(:, [5 2]);
ed.foo(:, 2) = ed.foo(:, 2)/2;
ed.foo(:, 5) = ed.foo(:, 5)*2;
d = dmce.export(dme, d0, {'beta', 'delta'});
t_ok(isequal(d, ed), [t 'd']);

t = sprintf('%s : import(dme, d, var_names, <scalar>) : ', k);
d = d0; d.foo = d0.foo(:, end:-1:1);    %% reverse columns
tab = tab0; tmp = tab.alpha(ridx1); tab.alpha(ridx1) = tab.epsilon(ridx1); tab.epsilon(ridx1) = tmp;
dme = dmce.import(mp.dme_test_widget(), d0);
dme = dmce.import(dme, d, {'alpha', 'epsilon'}, ridx1);
t_ok(isequal(dme.tab, tab), [t 'dme.tab']);

t = sprintf('%s : export(dme, d, var_names, <scalar>) : ', k);
tmp = dme.tab.beta; dme.tab.beta = dme.tab.delta; dme.tab.delta = tmp;
ed = d0; ed.foo(ridx1, [2 5]) = ed.foo(ridx1, [5 2]);
ed.foo(ridx1, 2) = ed.foo(ridx1, 2)/2;
ed.foo(ridx1, 5) = ed.foo(ridx1, 5)*2;
d = dmce.export(dme, d0, {'beta', 'delta'}, ridx1);
t_ok(isequal(d, ed), [t 'd']);

t = sprintf('%s : import(dme, d, var_names, <vector>) : ', k);
d = d0; d.foo = d0.foo(:, end:-1:1);    %% reverse columns
tab = tab0; tmp = tab.alpha(ridx2); tab.alpha(ridx2) = tab.epsilon(ridx2); tab.epsilon(ridx2) = tmp;
dme = dmce.import(mp.dme_test_widget(), d0);
dme = dmce.import(dme, d, {'alpha', 'epsilon'}, ridx2);
t_ok(isequal(dme.tab, tab), [t 'dme.tab']);

t = sprintf('%s : export(dme, d, var_names, <vector>) : ', k);
tmp = dme.tab.beta; dme.tab.beta = dme.tab.delta; dme.tab.delta = tmp;
ed = d0; ed.foo(ridx2, [2 5]) = ed.foo(ridx2, [5 2]);
ed.foo(ridx2, 2) = ed.foo(ridx2, 2)/2;
ed.foo(ridx2, 5) = ed.foo(ridx2, 5)*2;
d = dmce.export(dme, d0, {'beta', 'delta'}, ridx2);
t_ok(isequal(d, ed), [t 'd']);

%% 2 - dme table has r for indexing original data rows
k = 'B';
t = sprintf('%s : dmce : ', k);
dmce = mp.dmce_test_widget();
dmce.use_r = 1;
t_ok(isa(dmce, 'mp.dmc_element'), [t 'class']);
t_is(dmce.use_r, 1, 12, [t 'dmce.use_r == 1']);

t = sprintf('%s : import(dme, d) : ', k);
dme = dmce.import(mp.dme_test_widget(), d0);
t_ok(~isempty(dme.tab), [t 'dme.tab not empty']);
t_ok(isequal(dme.tab, tab1), [t 'dme.tab']);

t = sprintf('%s : export(dme, d) : ', k);
d = d0; d.foo = 0 * d.foo;
d = dmce.export(dme, d);
ed = d0;
ed.foo([1;3], :) = 0 * ed.foo([1;3], :);
ed.bar.baz(r,1) = ed.bar.baz(r,2);
t_ok(isequal(d, ed), [t 'd']);

t = sprintf('%s : import(dme, d, var_names) : ', k);
d = d0; d.foo(:, 1:end-2) = d0.foo(:, end-2:-1:1);  %% reverse columns
d.bar.baz(:, [1 2]) = d.bar.baz(:, [2 1]);          %% reverse columns
tab = tab1;
tmp = tab.alpha;
tab.alpha = tab.gamma / 2;
tab.gamma = tmp * 2;
tab.name = baz(r, 1);
dme = dmce.import(mp.dme_test_widget(), d0);
dme = dmce.import(dme, d, {'alpha', 'name', 'gamma'});
t_ok(isequal(dme.tab, tab), [t 'dme.tab']);

t = sprintf('%s : export(dme, d, var_names) : ', k);
tmp = dme.tab.beta; dme.tab.beta = dme.tab.delta; dme.tab.delta = tmp;
d = d0; d.foo = 0 * d.foo;
[d.bar.baz{:}] = deal('');
ed = d;
ed.foo(r, [2;5]) = d0.foo(r, [5;2]) * [0.5 0; 0 2];
ed.bar.baz(r,1) = d0.bar.baz(r,1);
d = dmce.export(dme, d, {'beta', 'delta', 'name'});
t_ok(isequal(d, ed), [t 'd']);

t = sprintf('%s : import(dme, d, var_names, <scalar>) : ', k);
d = d0; d.foo(:, 1:end-2) = d0.foo(:, end-2:-1:1);  %% reverse columns
d.bar.baz(:, [1 2]) = d.bar.baz(:, [2 1]);          %% reverse columns
tab = tab1;
tmp = tab.alpha(ridx1);
tab.alpha(ridx1) = tab.gamma(ridx1) / 2;
tab.gamma(ridx1) = tmp * 2;
tab.name(ridx1) = baz(r(ridx1), 1);
dme = dmce.import(mp.dme_test_widget(), d0);
dme = dmce.import(dme, d, {'alpha', 'name', 'gamma'}, ridx1);
t_ok(isequal(dme.tab, tab), [t 'dme.tab']);

t = sprintf('%s : export(dme, d, var_names, <scalar>) : ', k);
tmp = dme.tab.beta; dme.tab.beta = dme.tab.delta; dme.tab.delta = tmp;
d = d0; d.foo = 0 * d.foo;
[d.bar.baz{:}] = deal('');
ed = d;
ed.foo(r(ridx1), [2;5]) = d0.foo(r(ridx1), [5;2]) * [0.5 0; 0 2];
ed.bar.baz(r(ridx1),1) = d0.bar.baz(r(ridx1),1);
d = dmce.export(dme, d, {'beta', 'delta', 'name'}, ridx1);
t_ok(isequal(d, ed), [t 'd']);

t = sprintf('%s : import(dme, d, var_names, <vector>) : ', k);
d = d0; d.foo(:, 1:end-2) = d0.foo(:, end-2:-1:1);  %% reverse columns
d.bar.baz(:, [1 2]) = d.bar.baz(:, [2 1]);          %% reverse columns
tab = tab1;
tmp = tab.alpha(ridx2);
tab.alpha(ridx2) = tab.gamma(ridx2) / 2;
tab.gamma(ridx2) = tmp * 2;
tab.name(ridx2) = baz(r(ridx2), 1);
dme = dmce.import(mp.dme_test_widget(), d0);
dme = dmce.import(dme, d, {'alpha', 'name', 'gamma'}, ridx2);
t_ok(isequal(dme.tab, tab), [t 'dme.tab']);

t = sprintf('%s : export(dme, d, var_names, <vector>) : ', k);
tmp = dme.tab.beta; dme.tab.beta = dme.tab.delta; dme.tab.delta = tmp;
d = d0; d.foo = 0 * d.foo;
[d.bar.baz{:}] = deal('');
ed = d;
ed.foo(r(ridx2), [2;5]) = d0.foo(r(ridx2), [5;2]) * [0.5 0; 0 2];
ed.bar.baz(r(ridx2),1) = d0.bar.baz(r(ridx2),1);
d = dmce.export(dme, d, {'beta', 'delta', 'name'}, ridx2);
t_ok(isequal(d, ed), [t 'd']);

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%
% read_eccen_patch.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

% Should run 'left or right' scripts to set up fname_* variable,s
% and set nfig so that plotting starts in the right location.
fclose('all');
fid=fopen(fname_patch,'rt');
s = fgetl(fid);
s = fgetl(fid);
[poly,count] = sscanf(s,'%d');
numvert = poly(1);
numquad = poly(2);
vertex_coordinates = zeros(3,1);
face_index = zeros(4,1);
vertx_list = zeros(numvert,3);
face_list = zeros(numquad,5);
tic;
for vert = 1:1:(numvert),
    s = fgetl(fid);
    vertx = sscanf(s,'%d');
    s = fgetl(fid);
    vertx_coordinates = sscanf(s,'%f');
    vertx_list(vert,:) = [vertx_coordinates(1:2)' vertx];
end;
toc
tic;
for face = 1:1:(numquad),
    s = fgetl(fid);
    facenum = sscanf(s,'%d');
    s = fgetl(fid);
    face_vertx = sscanf(s,'%d');
    face_list(face,:) = [face_vertx' facenum];
end;
toc
fclose(fid);
full_vertx=zeros(max(max(vertx_list))+1,3);
full_vertx(vertx_list(:,3)+1,1)=vertx_list(:,1);
full_vertx(vertx_list(:,3)+1,2)=vertx_list(:,2);
mesh_real=File2Var(fname_real);
vertx_values=zeros(max(max(vertx_list))+1,1);
vertx_values(vertx_list(:,3)+1)=mesh_real(:,5);
mesh_imag=File2Var(fname_imag);
vertx_complex=zeros(max(max(vertx_list))+1,1);
vertx_complex(vertx_list(:,3)+1)=mesh_imag(:,5);
v_complex = vertx_values(:,1) + i*vertx_complex;
%v_phase=angle(v_complex');

strpwd=pwd;
%subplot(4,3,nfig+1);
figure;
title([strpwd((length(strpwd)-25):length(strpwd)) ' ' fname_real]);
p_eccen_handle=patch('Vertices',full_vertx,'Faces',face_list(:,1:4)+1,'FaceVertexCData',angle(v_complex),'FaceColor','interp','EdgeColor','none');
%colormap(rgb(256));
colorbar;
%tic; t_unwrap=unwrap(angle(v_complex)); toc

short_complex = mesh_real(:,5) + (i*mesh_imag(:,5));
short_unwrap=angle(short_complex);
[zmat_eccen,xvec,yvec]=ffgrid(vertx_list(:,1),vertx_list(:,2),short_unwrap,1.0,1.0);
figure;
title([strpwd((length(strpwd)-25):length(strpwd)) ' ' fname_real]);
%subplot(4,3,nfig+2);
imagesc(xvec,yvec,zmat_eccen);
axis xy;
%colormap(rgb(64));
colorbar;
%subplot(4,3,nfig+3);
figure;
title([strpwd((length(strpwd)-25):length(strpwd)) ' ' fname_real]);
contour(xvec,yvec,zmat_eccen,20);
%figure;

%ffgrid(vertx_list(:,1),vertx_list(:,2),short_unwrap,0.75,0.75);
[zmat_imag_eccen,xvec,yvec]=ffgrid(vertx_list(:,1),vertx_list(:,2),mesh_imag(:,5),1.0,1.0);
[zmat_real_eccen,xvec,yvec]=ffgrid(vertx_list(:,1),vertx_list(:,2),mesh_real(:,5),1.0,1.0);
smooth_filt = fspecial('average',3);
zmat_eccen_cplx = zmat_real_eccen + (i*zmat_imag_eccen);
zmat_imag_smooth = zmat_imag_eccen;
zmat_real_smooth = zmat_real_eccen;
maxiter=2;
for niter=1:maxiter,
  zmat_imag_smooth = filter2(smooth_filt,zmat_imag_smooth);
  zmat_real_smooth = filter2(smooth_filt,zmat_real_smooth);
end;
zmat_smooth_eccen = angle(zmat_real_smooth + (i*zmat_imag_smooth));
zmat_smooth_eccen_cplx = zmat_real_smooth + (i*zmat_imag_smooth);
%subplot(4,3,nfig+5);
figure;
title([strpwd((length(strpwd)-25):length(strpwd)) ' ' fname_real]);
imagesc(xvec,yvec,zmat_smooth_eccen); axis xy;
figure;
%subplot(4,3,nfig+6);
title([strpwd((length(strpwd)-25):length(strpwd)) ' ' fname_real]);
imcontour(xvec,yvec,zmat_smooth_eccen_cplx,20);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%Hand-eye calibration
%
%This add on allows to compute the hand to eye calibration based ont he
%calibration toolbox from Jean-Yves Bouguet
%see http://www.vision.caltech.edu/bouguetj/calib_doc/
%more information can be found here:
%http://www.vision.ee.ethz.ch/~cwengert/calibration_toolbox.php
%
%You need to have computed the intrinsic and extrinsic parameters of your
%camera (the extrinsic parameters with respect to the grid's local
%coordinate system). Then you also need to supply the pose / transformation
%of the robot arm / marker with repsect to the robot base / external
%tracking device in the following format:
%Pack into a 4x4xNumber_of_Views Matrix the following data
%Hmarker2world(:,:,i) = [Ri_3x3 ti_3x1;[ 0 0 0 1]]
%with
%i = number of the view,
%Ri_3x3 the rotation matrix
%ti_3x1 the translation vector.
%
%The following parameters can be set (if not set the default values are
%used):
%doShow [default=0]:
%Display the results from the calibration
%doSortHandEyeMovement [default=0]:
%Set this to 1 if you want your views sorted in a way that the interstation
%movement is ideal for hand-eye calibration, see [Tsai]
%HandEyeMethod	[default = 'Tsai']
%Here you specifiy which method you want to use for the
%hand-eye calibration, default is using Tsai's method. See
%ftp://ftp.vision.ee.ethz.ch/publications/proceedings/eth_biwi_00363.pdf
%Possible values are 'Tsai', 'Inria', 'Dual_quaternion', 'Navy'
%Tsai, Inria and Navy give the same results and are usually a bit better
%than the Dual quaternion approach
%
%Christian Wengert
%Computer Vision Laboratory
%ETH Zurich
%Sternwartstrasse 7
%CH-8092 Zurich
%www.vision.ee.ethz.ch/cwengert
%wengert@vision.ee.ethz.ch

%With the doShow flag you can show the results
if(~exist('doShow'))
    doShow = 0;
end

%Tsai states that the inter-station movement should be as large as possible
%for better accuracy. This flag will sort the movements for higher
%accuracy
if(~exist('doSortHandEyeMovement'))
    doSortHandEyeMovement = 0;
end

%If it does, check whether the Hrobot2hand exists
if(~exist('Hmarker2world'))
    disp(['handeye:: No Hmarker2world data available, Hand eye calibration aborted']);
    return
end

%If it does, check whether the Hrobot2hand exists
if(~exist('HandEyeMethod'))
    HandEyeMethod = 'Tsai';
end

%Now go on
disp(['handeye:: Hand-Eye Calibration']);
correctSets = 0;
for i=1:length(active_images)
    if(active_images(i)~=0 && Hmarker2world(1,1,i)~=0) %Make sure its a calibration image and it has tracker info
        %Extract extrinsic parameters
        correctSets = correctSets+1;
        eval(['Rc = Rc_' num2str(active_images(i)) ';']);%%%%%%%%%%%%%%%%%%%%%ACHTUNG
        eval(['Tc = Tc_' num2str(active_images(i)) ';']);
        Hgrid2cam(:,:,correctSets) = inv(([Rc Tc;[ 0 0 0 1]]));
        Hcam2grid(:,:,correctSets) = inv(Hgrid2cam(:,:,correctSets));
        Hm2w(:,:,correctSets) = Hmarker2world(:,:,active_images(i));
    end
end
if(correctSets)
    %Remove bad trackerdata
    badTrackerDataIndex = [];
    for i=1:correctSets
        if(abs(Hm2w(1,1,i))<1e-18) %Its bad
            badTrackerDataIndex = [badTrackerDataIndex;i];
            correctSets = correctSets-1;
        end
    end
    %Init
    Hm2w(:,:,badTrackerDataIndex) = [];
    Hcam2grid(:,:,badTrackerDataIndex) = [];
    Hgrid2cam(:,:,badTrackerDataIndex) = [];
    if(doSortHandEyeMovement)
        index = sortHandEyeMovement(Hm2w);
    else
        index = 1:size(Hm2w,3);
    end
    Hm2w2 = Hm2w(:,:,index);
    Hcam2grid2 = Hcam2grid(:,:,index);
    %Now calibrate
    switch(HandEyeMethod)
        case 'Tsai'
            [Hcam2marker_, err] = TSAIleastSquareCalibration(Hm2w2, Hcam2grid2)
        case 'Inria'
            [Hcam2marker_, err] = inria_calibration(Hm2w, Hcam2grid2);
        case 'Navy'
            [Hcam2marker_, err] = navy_calibration(Hm2w, Hcam2grid2);
        case 'Dual_quaternion'
            [Hcam2marker_, err] = hand_eye_dual_quaternion(Hm2w, Hcam2grid) ;
    end
    %Create the average Hworld2grid, givin an idea where the grid is
    %in the coordinate system of the tracker/robot
    for i=1:correctSets
        Hcam2world_(:,:,i) = Hm2w(:,:,i)*Hcam2marker_;  %Hc2m(:,:,k)
        Hworld2cam_(:,:,i) = inv(Hcam2world_(:,:,i));
        %The above is correct as it gives the same as in my
        %simulation
        Hgrid2world_(:,:,i) = Hcam2world_(:,:,i)*Hcam2grid(:,:,i);
        Hworld2grid_(:,:,i) = inv(Hgrid2world_(:,:,i));
    end
    %Average it, using the algorithm described in ...
    Hgrid2worldAvg = averageTransformation(Hgrid2world_);
    xd = [];

    %BACKPROJECT, this computes and displays the result
    errAvg = [];
    correctSets = 0;
    for i=1:length(active_images)
        if(active_images(i)~=0)
            correctSets  = correctSets + 1;
            if(correctSets<=length(Hm2w))
                err = [];
                eval(['tmp = X_' num2str(active_images(i)) ';'])
                eval(['ALLPTS.X_' num2str(correctSets) ' = X_' num2str(active_images(i)) ';'])
                tmp = [tmp, [[0;0;0],[3;0;0]]];
                eval(['xd_' num2str(correctSets) ' = x_' num2str(active_images(i)) ';']);
                eval(['allpts.xd_' num2str(correctSets) ' = x_' num2str(active_images(i)) ';']);
                Xworld = Hgrid2world_(:,:,correctSets)*[tmp;ones(1,length(tmp))];
                XworldAvg = Hgrid2worldAvg*[tmp;ones(1,length(tmp))];     %
                [P,x] = backprojectDistorted(KK,Hworld2cam_(1:3,1:3,correctSets),Hworld2cam_(1:3,4,correctSets),XworldAvg,kc);

                %Compute error
                eval(['xd = xd_' num2str(correctSets) ';']);
                m = length(xd);
                for j=1:m
                    err(j) = norm(x(:,j)-xd(:,j));
                end
                %per point
                errAvg(correctSets) = sum(err)/m;
                %Store it
                result.error_calib_avg(correctSets) = errAvg(correctSets);
                %GFX
                if(doShow)
                    hold off,figure(n_ima+1),hold on
                    if(isfloat(im))
                        eval(['imshow(I_' num2str(i) '/255);']);,hold on
                    else
                        eval(['imshow(I_' num2str(i) ');']);,hold on
                    end

                    hold on
                    draw2DPoints(x(:,1:m), '', 0, n_ima+1)
                    plot(x(1,end-1:end),x(2,end-1:end),'rx','MarkerSize',25,'LineWidth',3);
                    draw2DPoints(xd, 'blue(real), red guessed', 0, n_ima+1,'bd')
                    title(['Error = ' num2str(sum(err)/m) ' image ' num2str(i)]);
                    axis equal
                    pause
                end
            end
        end
    end
    backprojection_error = sum(errAvg)/correctSets
end

author: christianwengert


text:
function mpc = order_radial(mpc)
%ORDER_RADIAL  Performs oriented ordering to buses and branches.
%
%   mpc = order_radial(mpc)
%
%   orders the branches by the the principle of oriented ordering:
%   indicies of sending nodes are smaller than the indicies of the
%   receiving nodes. The branch index is equal to the index of their
%   receiving node. The ordering is taken from:
%   D. Rajicic, R. Ackovski and R. Taleski, "Voltage correction power flow,"
%   IEEE Transactions on Power Delivery, vol. 9, no. 2, pp. 1056-1062, Apr 1994.
%
%   See also RADIAL_PF.

%% define named indices into bus, gen, branch matrices
define_constants;
%% Initialize
slack = mpc.bus(mpc.bus(:,BUS_TYPE) == REF, 1);
[f, t] = deal(mpc.branch(:,F_BUS),mpc.branch(:,T_BUS));
nl = size(mpc.branch,1);
branch_number = (1:nl)';
mpc.branch_order = [];
mpc.loop = [];
mpc.bus_order = slack;
%% Bus and branch ordering
iter = 1;
while ~isempty(f) && iter <= nl
    % For each of the "from" buses that are in bus_order,
    % add the corresponding "to" buses to it. If both "from" and "to" are
    % in bus_order, the branch is forming loop. Add corresponding branch
    % numbers to branch_order or to loop.
    mf = ismember(f,mpc.bus_order);
    mt = ismember(t,mpc.bus_order);
    is_loop = mf & mt;
    % Add branch to loop and delete it from the list
    if any(is_loop)
        mpc.loop = [mpc.loop; branch_number(is_loop)];
        mf(is_loop) = [];
%         mt(is_loop) = [];
        f(is_loop) = [];
        t(is_loop) = [];
        branch_number(is_loop) = [];
    end
    % Add unique buses to bus_order
    if any(mf)
        u = unique(t(mf)); % unique buses
        [junk,i] = intersect(t.*mf,u); % first occurence of unique buses in t
        mpc.bus_order = [mpc.bus_order; t(i)];
        % Add branch to branch_order and delete it from the list
        mpc.branch_order = [mpc.branch_order; branch_number(i)];
        mf(i) = [];
%         mt(i) = [];
        f(i) = [];
        t(i) = [];
        branch_number(i) = [];
    end
    % Add any remaining branch to loop and delete it from the list
    if any(mf)
        mpc.loop = [mpc.loop; branch_number(mf)];
        f(mf) = [];
        t(mf) = [];
        branch_number(mf) = [];
    end
    % For each of the "to" buses that are in bus_order,
    % add the corresponding "from" buses to it. If both "from" and "to" are
    % in bus_order, the branch is forming loop. Add corresponding branch
    % numbers to branch_order or to loop.
    mf = ismember(f,mpc.bus_order);
    mt = ismember(t,mpc.bus_order);
    is_loop = mf & mt;
    % Add branch to loop and delete it from the list
    if any(is_loop)
        mpc.loop = [mpc.loop; branch_number(is_loop)];
        mt(is_loop) = [];
%         mf(is_loop) = [];
        f(is_loop) = [];
        t(is_loop) = [];
        branch_number(is_loop) = [];
    end
    % Add unique buses to bus_order
    if any(mt)
        u = unique(f(mt)); % unique buses
        [junk,i] = intersect(f.*mt,u); % first occurence of unique buses in f
        mpc.bus_order = [mpc.bus_order; f(i)];
        % Add branch to branch_order and delete it from the list
        mpc.branch_order = [mpc.branch_order; branch_number(i)];
%         mf(i) = [];
        mt(i) = [];
        f(i) = [];
        t(i) = [];
        branch_number(i) = [];
    end
    % Add any remaining branch to loop and delete it from the list
    if any(mt)
        mpc.loop = [mpc.loop; branch_number(mt)];
        f(mt) = [];
        t(mt) = [];
        branch_number(mt) = [];
    end
    iter = iter + 1;
end
if ~isempty(f)
    mpc.not_connected = branch_number;
else
    mpc.not_connected = [];
end
%% Reorder bus, branch and gen.
if isempty(mpc.loop)
    % Permute rows in branch
    mpc.branch = mpc.branch(mpc.branch_order,:);
    % Make an "inverse" vector out of bus_order
    mpc.bus_order_inv = sparse(mpc.bus_order,ones(nl+1,1),1:nl+1);
    % Swap indicies in "from" and "to" buses using bus_order_inv
    [f, t] = deal(mpc.branch(:,F_BUS),mpc.branch(:,T_BUS));
    f = mpc.bus_order_inv(f);
    t = mpc.bus_order_inv(t);
    % Reverse branch orientation of "from" is biger than "to"
    mpc.br_reverse = f > t;
    [f(mpc.br_reverse), t(mpc.br_reverse)] = deal(t(mpc.br_reverse), f(mpc.br_reverse));
    % Put new "from" and "to" indicies in branch
    mpc.branch(:,[F_BUS T_BUS]) = [f t];
    % Make an "inverse" vector out of branch_order
    mpc.branch_order_inv = sparse(mpc.branch_order,ones(nl,1),1:nl);
    % Permute rows in bus and replace bus indicies
    mpc.bus = mpc.bus(mpc.bus_order,:);
    mpc.bus(:,BUS_I) = mpc.bus_order_inv(mpc.bus(:,BUS_I));
    % Replace bus indicies in gen
    mpc.gen(:,GEN_BUS) = mpc.bus_order_inv(mpc.gen(:,GEN_BUS));
end


text:
%computes the maximum of the Harmonic Product Spectrum
%> called by ::ComputePitch
%>
%> @param X: spectrogram (dimension FFTLength X Observations)
%> @param f_s: sample rate of audio data
%>
%> @retval f_0 HPS maximum (in Hz)
% ======================================================================
function [f_0] = PitchSpectralHps (X, f_s)

    % initialize
    iOrder = 4;
    f_min = 300;
    k_min = round(f_min/f_s * 2 * (size(X, 1)-1)) + 1;

    afHps = X;

    % compute the HPS
    for j = 2:iOrder
        afHps = afHps .* [X(1:j:end, :); zeros(size(X, 1)-size(X(1:j:end,:), 1), size(X, 2))];
    end

    % find max index and convert to Hz
    [fDummy, f_0] = max(afHps(k_min:end, :), [], 1);
    f_0 = (f_0 + k_min - 2) / (size(X, 1)-1) * f_s / 2;
    f_0(sum(afHps, 1) == 0) = 0;
end

author: alexanderlerch

repo: ACA-Code


text:
function [alpha,xp] = StepSizeSW(f,x,d,alpha,params)
% Line search algorithm satisfying strong Wolfe conditions.
% Algorithms 3.5 on pages 60-61 in Nocedal and Wright.
% Requires x.p, x.f and x.g to be initialized.

alpha0 = params.stpmin;
c1 = params.ftol;
c2 = params.gtol;
% alpha is alpha_i
gxd = x.g'*d;
% alphap is alpha_{i-1}
alphap = alpha0;
fxp = x.f;
i=1;
while norm(alphap-alpha) > params.xtol
  xp.p = x.p + alpha*d;
  xp.f = feval(f,xp.p,1);
  if (xp.f > x.f + c1*alpha*gxd) | ((i > 1) & (xp.f >= fxp)),
    [alpha,xp] = zoom(f,x,d,alphap,alpha,fxp,c1,c2);
    return;
  end
  xp.g = feval(f,xp.p,2); gxpd = xp.g'*d;
  if abs(gxpd) <= -c2*gxd,
    return;
  end
  if gxpd >= 0,
    [alpha,xp] = zoom(f,x,d,alpha,alphap,xp.f,c1,c2);
    return;
  end
  alphap = alpha;
  fxp = xp.f;
  %  alpha = alpha + (params.stpmax-alpha)*rand(1);
  alpha = alpha + (params.stpmax-alpha)*0.5;
  i = i+1;
end
alpha,
error('No stepsize found');

function [alpha,xp] = zoom(f,x,d,alphal,alphah,fxl,c1,c2)
% function [alpha,xp] = zoom(f,x,d,alphal,alphah,fxl)
% Algorithm 3.6 on page 61 in Nocedal and Wright

gxd = x.g'*d;

while 1
   alpha = 1/2*(alphal+alphah);
   xp.p = x.p + alpha*d;
   xp.f = feval(f,xp.p,1);
   if ((xp.f > x.f + c1*alpha*gxd) | (xp.f >= fxl)),
      alphah = alpha;
   else
      xp.g = feval(f,xp.p,2); gxpd = xp.g'*d;
      if abs(gxpd) <= -c2*gxd,
        return;
      end
      if gxpd*(alphah-alphal) >= 0,
        alphah = alphal;
      end
      alphal = alpha;
      fxl = xp.f;
   end
end

author: clarkzinzow

repo: Nonlinear-Optimization-Algorithms


text:
function [NNF, nUpdate]= sc_update_NNF(trgPatch, img, NNF, modelPlane, modelReg, optS)

% SC_UPDATE_NNF
%
% Update the nearest neighbor field using the PatchMatch algorithm
%
% Input:
%   - trgPatch
%   - img:
%   - NNF
%   - modelPlane
%   - modelReg
%   - optS
% Output:
%   - NNF
%   - nUpdate

nUpdate = zeros(1,3);

for i = 1:optS.numPassPerIter
    % propagate along four directions
    for iDirect = 1:4
        [NNF, n] = sc_propagate(trgPatch, img, NNF, modelPlane, optS, iDirect);
        nUpdate(1) = nUpdate(1) + n;
    end

    % Random sampling
    %     [NNF, n] = sc_regular_search(trgPatch, img, NNF, modelPlane, modelReg, optS);
    %     nUpdate(3) = nUpdate(3) + n;

    % Regularity guided sampling
    [NNF, n] = sc_random_search(trgPatch, img, NNF, modelPlane, optS);
    nUpdate(2) = nUpdate(2) + n;
end

end

% PatchMatch - propagate step
function [NNF, nUpdateTotal] = sc_propagate(trgPatch, img, NNF, modelPlane, optS, iDirect)

% SC_PROPAGATE: update the nearest neighbor field using propagation
% Input:
%   - trgPatch:     target patch
%   - img:          image at current level
%   - NNF:          current nearest neighbor field
%   - modelPlane:   planar structure model
%   - optS
%   - indDirection: propagation direction
% Output:
%   - NNF, nUpdateTotal

nUpdateTotal = 0;

% The positions of neighboring pixels
uvPixN = NNF.uvPixN{iDirect};
uvPixActiveInd = uvPixN.validInd;

numUpdatePix = NNF.uvPix.numPix;
while(numUpdatePix ~= 0)
    % Prepare uvPix, uvPixNCur
    uvPix.sub     = NNF.uvPix.sub(uvPixActiveInd, :);  uvPix.ind     = NNF.uvPix.ind(uvPixActiveInd, :);
    uvPixNCur.sub = uvPixN.sub(uvPixActiveInd, :);     uvPixNCur.ind = uvPixN.ind(uvPixActiveInd, :);

    uvDtBdPixPosCur = NNF.uvDtBdPixPos(uvPixActiveInd);

    trgPatchCur   = trgPatch(:,:, uvPixActiveInd);
    srcPosCur     = NNF.uvTform.data(uvPixActiveInd, 7:8);
    uvCostCur     = NNF.uvCost.data(uvPixActiveInd);
    uvPlaneIDCur  = NNF.uvPlaneID.data(uvPixActiveInd);

    % Get srcPosMap for computing coherence cost
    srcPosMapCur  = NNF.uvTform.map(:,:,7:8);

    % Active pixel positions
    uvPixActivePos = find(uvPixActiveInd);

    % Get candidate uvTform candidates
    uvTformCand = sc_uvMat_from_uvMap(NNF.uvTform.map, uvPixNCur.ind);

    % Generate candidate transformation by propagation
    uvTformCand = sc_trans_tform(uvTformCand, optS.propDir(iDirect,:));

    % Check if the nearest neighbors are valid source patches
    uvValidSrcInd = sc_check_valid_uv(uvTformCand(:,7:8), NNF.validPix.mask);

    % Check if the nearest neighbors are already the same as the existing one
    diff = abs(uvTformCand(:,7:8) - srcPosCur);
    uvValidDistInd = ((diff(:,1) > 1 ) | (diff(:,2) > 1 ));

    % Valid pixel indices
    uvValidInd = uvValidSrcInd & uvValidDistInd;

    numUvValid = sum(uvValidInd);

    if(numUvValid ~= 0)
        uvPixValid.sub = uvPix.sub(uvValidInd,:);
        uvPixValid.ind = uvPix.ind(uvValidInd);

        uvDtBdPixPosCur = uvDtBdPixPosCur(uvValidInd);

        trgPatchCur    = trgPatchCur(:,:, uvValidInd);
        uvTformCand    = uvTformCand(uvValidInd, :);
        uvCostCur      = uvCostCur(uvValidInd);
        uvPlaneIDCand  = uvPlaneIDCur(uvValidInd);

        uvPixUpdatePos = uvPixActivePos(uvValidInd);

        % Grab source patches
        srcPatch = sc_prep_source_patch(img, uvTformCand, optS);

        % Compute patch matching cost
        [costPatchCandAll, uvBiasCand] = ...
            sc_patch_cost(trgPatchCur, srcPatch, modelPlane, uvPlaneIDCand, ...
            uvPixValid.sub, uvTformCand, srcPosMapCur, uvDtBdPixPosCur, optS);
        costPatchCand = sum(costPatchCandAll, 2);

        % Check which one to update
        updateInd = costPatchCand < uvCostCur;

        uvPixUpdatePos = uvPixUpdatePos(updateInd);
        numUpdatePix = size(uvPixUpdatePos, 1);
    else
        numUpdatePix = 0;
    end

    % Update NNF data
    if(numUpdatePix ~= 0)
        nUpdateTotal = nUpdateTotal + numUpdatePix;

        % === Update NNF data ===
        NNF.uvTform.data(uvPixUpdatePos, :) = uvTformCand(updateInd,:);
        NNF.uvCost.data(uvPixUpdatePos)     = costPatchCand(updateInd);
        NNF.uvPlaneID.data(uvPixUpdatePos)  = uvPlaneIDCand(updateInd);

        % Apply bias correction
        if(optS.useBiasCorrection)
            NNF.uvBias.data(:,:,uvPixUpdatePos)   = uvBiasCand(:,:,updateInd);
        end
        %         NNF.update.data(uvPixUpdatePos)   = 1; % updateInd;

        % === Update NNF map ===
        uvPixValidInd = uvPixValid.ind(updateInd);
        NNF.uvTform.map    = sc_update_uvMap(NNF.uvTform.map, uvTformCand(updateInd,:), uvPixValidInd);
        NNF.uvCost.map     = sc_update_uvMap(NNF.uvCost.map, costPatchCand(updateInd), uvPixValidInd);
        NNF.uvPlaneID.map  = sc_update_uvMap(NNF.uvPlaneID.map, uvPlaneIDCand(updateInd), uvPixValidInd);

        %         NNF.update.map  = sc_update_uvMap(NNF.update.map, 1, uvPixValidInd);

        % === Update uvPixActiveInd ===
        uvPixNextSub = uvPixValid.sub(updateInd,:);
        uvPixNextSub = bsxfun(@plus, uvPixNextSub, optS.propDir(iDirect,:));
        uvPixNextInd = sub2ind([NNF.imgH, NNF.imgW], uvPixNextSub(:,2), uvPixNextSub(:,1));

        updateMap = NNF.uvPix.mask;
        updateMap(uvPixNextInd) = 0;
        uvPixActiveInd = ~updateMap(NNF.uvPix.ind);
        uvPixActiveInd = uvPixActiveInd & uvPixN.validInd;
    end
end

end

% PatchMatch - random search step
function [NNF, nUpdateTotal] = sc_random_search(trgPatch, img, NNF, modelPlane, optS)

% SC_RANDOM_SEARCH: update the nearest neighbor using random sampling

[imgH, imgW, nCh] = size(img);

uvPix    = NNF.uvPix;
numUvPix = size(uvPix.sub, 1);

searchRad = max(imgH, imgW)/2;
nUpdateTotal = 0;

% uvPixActiveInd = true(numUvPix, 1);
iter = 1;
% while(iter <= optS.numRandSample)
while(searchRad > 1)
    iter = iter + 1;
    % Reduce search radius by half
    searchRad = searchRad/2;
    if(searchRad < 1)
        break;
    end

    % Get srcPosMap for computing coherence cost
    srcPosMapCur  = NNF.uvTform.map(:,:,7:8);

    uvTformCandCur = sc_uvMat_from_uvMap(NNF.uvTform.map, uvPix.ind);

    % Draw random samples
    srcPos = uvTformCandCur(:,7:8) + 2*searchRad*(rand(numUvPix, 2) - 0.5);

    % Draw plane ID candidate
    uvPlaneIDCand = sc_draw_plane_id(NNF.uvPlaneID.planeProbAcc);

    % Estimate the domain transformation
    uvTformCand = sc_src_domain_tform(uvPlaneIDCand, modelPlane, [], srcPos, NNF.uvPix.sub, 1);

    % === Reject invalid samples ===
    % Check if the scale of the source patch valid
    uvTformScale    = sc_scale_tform(uvTformCand);
    uvValidScaleInd = (uvTformScale > optS.minScale) & (uvTformScale < optS.maxScale);
    % Check if the souce patch is valid
    uvValidSrcInd   = sc_check_valid_uv(uvTformCand(:,7:8), NNF.validPix.mask);

    uvValidInd = uvValidSrcInd & uvValidScaleInd;

    % Active pixels
    uvPixActivePos = find(uvValidInd);
    numActPix = size(uvPixActivePos, 1);

    if(numActPix~=0)
        % Update
        trgPatchCur      = trgPatch(:,:,uvValidInd);
        uvCostDataCur    = NNF.uvCost.data(uvValidInd);
        uvTformCandCur   = uvTformCand(uvValidInd, :);
        uvPlaneIDCandCur = uvPlaneIDCand(uvValidInd);

        uvPixValid.sub   = uvPix.sub(uvValidInd,:);
        uvPixValid.ind   = uvPix.ind(uvValidInd);

        uvDtBdPixPosCur = NNF.uvDtBdPixPos(uvValidInd);

        % Grab source patches
        srcPatch = sc_prep_source_patch(img, uvTformCandCur, optS);

        [costPatchCandAll, uvBiasCand] = ...
            sc_patch_cost(trgPatchCur, srcPatch, modelPlane, uvPlaneIDCandCur, ...
            uvPixValid.sub, uvTformCandCur, srcPosMapCur, uvDtBdPixPosCur, optS);
        costPatchCand = sum(costPatchCandAll, 2);
        % Check which one to update
        updateInd = (costPatchCand < uvCostDataCur);
        nUpdate = sum(updateInd);

        if(nUpdate~=0)
            uvPixActivePos = uvPixActivePos(updateInd);

            nUpdateTotal = nUpdateTotal + nUpdate;

            % === Update NNF data ===
            NNF.uvTform.data(uvPixActivePos, :) = uvTformCandCur(updateInd,:);
            NNF.uvPlaneID.data(uvPixActivePos)  = uvPlaneIDCandCur(updateInd);
            NNF.uvCost.data(uvPixActivePos)     = costPatchCand(updateInd);
            if(optS.useBiasCorrection)
                NNF.uvBias.data(:,:,uvPixActivePos)   = uvBiasCand(:,:,updateInd);
            end

            % === Update NNF map ===
            uvPixValidInd     = uvPixValid.ind(updateInd);
            NNF.uvTform.map   = sc_update_uvMap(NNF.uvTform.map, uvTformCandCur(updateInd,:), uvPixValidInd);
            NNF.uvPlaneID.map = sc_update_uvMap(NNF.uvPlaneID.map, uvPlaneIDCandCur(updateInd), uvPixValidInd);
            NNF.uvCost.map    = sc_update_uvMap(NNF.uvCost.map, costPatchCand(updateInd), uvPixValidInd);
        end
    end
end
end

repo: StructCompletion


text:
function test_failed=test_involute

Lr=[9,10];

test_failed=0;

disp(' ===============  TEST_INVOLUTE ===========');

for ii=1:length(Lr)

  L=Lr(ii);
  f=tester_crand(L,1);

  r1=conj(dft(f));
  r2=dft(involute(f));

  res=norm(r1-r2);
  [test_failed,fail]=ltfatdiditfail(res,test_failed);
  s=sprintf('INVOLUTE  L:%3i %0.5g %s',L,res,fail);
  disp(s);

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function value = r8_psi ( xx )

%*****************************************************************************80
%
%% R8_PSI evaluates the function Psi(X).
%
%  Discussion:
%
%    This routine evaluates the logarithmic derivative of the
%    Gamma function,
%
%      PSI(X) = d/dX ( GAMMA(X) ) / GAMMA(X)
%             = d/dX LN ( GAMMA(X) )
%
%    for real X, where either
%
%      - XMAX1 < X < - XMIN, and X is not a negative integer,
%
%    or
%
%      XMIN < X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2008
%
%  Author:
%
%    Original FORTRAN77 version by William Cody.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    William Cody, Anthony Strecok, Henry Thacher,
%    Chebyshev Approximations for the Psi Function,
%    Mathematics of Computation,
%    Volume 27, Number 121, January 1973, pages 123-127.
%
%  Parameters:
%
%    Input, real XX, the argument of the function.
%
%    Output, real VALUE, the value of the function.
%
  four = 4.0;
  fourth = 0.25;
  half = 0.5;
  one = 1.0;
  p1(1:9) = [ ...
   4.5104681245762934160E-03, ...
   5.4932855833000385356, ...
   3.7646693175929276856E+02, ...
   7.9525490849151998065E+03, ...
   7.1451595818951933210E+04, ...
   3.0655976301987365674E+05, ...
   6.3606997788964458797E+05, ...
   5.8041312783537569993E+05, ...
   1.6585695029761022321E+05 ];
  p2(1:7) = [ ...
  -2.7103228277757834192, ...
  -1.5166271776896121383E+01, ...
  -1.9784554148719218667E+01, ...
  -8.8100958828312219821, ...
  -1.4479614616899842986, ...
  -7.3689600332394549911E-02, ...
  -6.5135387732718171306E-21 ];
  piov4 = 0.78539816339744830962;
  q1(1:8) = [ ...
   9.6141654774222358525E+01, ...
   2.6287715790581193330E+03, ...
   2.9862497022250277920E+04, ...
   1.6206566091533671639E+05, ...
   4.3487880712768329037E+05, ...
   5.4256384537269993733E+05, ...
   2.4242185002017985252E+05, ...
   6.4155223783576225996E-08 ];
  q2(1:6) = [ ...
   4.4992760373789365846E+01, ...
   2.0240955312679931159E+02, ...
   2.4736979003315290057E+02, ...
   1.0742543875702278326E+02, ...
   1.7463965060678569906E+01, ...
   8.8427520398873480342E-01 ];
  three = 3.0;;
  x01 = 187.0;
  x01d = 128.0;
  x02 = 6.9464496836234126266E-04;
  xinf = 1.70E+38;
  xlarge = 2.04E+15;
  xmax1 = 3.60E+16;
  xmin1 = 5.89E-39;
  xsmall = 2.05E-09;
  zero = 0.0;

  x = xx;
  w = abs ( x );
  aug = zero;
%
%  Check for valid arguments, then branch to appropriate algorithm.
%
  if ( xmax1 <= - x | w < xmin1 )

    if ( zero < x )
      value = - xinf;
    else
      value = xinf;
    end

    return
  end

  if ( x < half )
%
%  X < 0.5, use reflection formula: psi(1-x) = psi(x) + pi * cot(pi*x)
%  Use 1/X for PI*COTAN(PI*X)  when  XMIN1 < |X| <= XSMALL.
%
    if ( w <= xsmall )

      aug = - one / x;
%
%  Argument reduction for cotangent.
%
    else

      if ( x < zero )
        sgn = piov4;
      else
        sgn = - piov4;
      end

      w = w - floor ( w );
      nq = floor ( w * four );
      w = four * ( w - nq * fourth );
%
%  W is now related to the fractional part of 4.0 * X.
%  Adjust argument to correspond to values in the first
%  quadrant and determine the sign.
%
      n = floor ( nq / 2 );

      if ( n + n ~= nq )
        w = one - w;
      end

      z = piov4 * w;

      if ( mod ( n, 2 ) ~= 0 )
        sgn = - sgn;
      end
%
%  Determine the final value for  -pi * cotan(pi*x).
%
      n = floor ( ( nq + 1 ) / 2 );
      if ( mod ( n, 2 ) == 0 )
%
%  Check for singularity.
%
        if ( z == zero )

          if ( zero < x )
            value = -xinf;
          else
            value = xinf;
          end

          return
        end

        aug = sgn * ( four / tan ( z ) );

      else

        aug = sgn * ( four * tan ( z ) );

      end

    end

    x = one - x;

  end
%
%  0.5 <= X <= 3.0.
%
  if ( x <= three )

    den = x;
    upper = p1(1) * x;
    for i = 1 : 7
      den = ( den + q1(i) ) * x;
      upper = ( upper + p1(i+1) ) * x;
    end
    den = ( upper + p1(9) ) / ( den + q1(8) );
    x = ( x - x01 / x01d ) - x02;
    value = den * x + aug;
    return

  end
%
%  3.0 < X.
%
  if ( x < xlarge )
    w = one / ( x * x );
    den = w;
    upper = p2(1) * w;
    for i = 1 : 5
      den = ( den + q2(i) ) * w;
      upper = ( upper + p2(i+1) ) * w;
    end
    aug = ( upper + p2(7) ) / ( den + q2(6) ) - half / x + aug;
  end

  value = aug + log ( x );

  return
end


text:
% THE DIAGNOSTIC TOOLS (in the diag-folder):
%
% Covergence diagnostics
%   PSRF     - Potential Scale Reduction Factor
%   CPSRF    - Cumulative Potential Scale Reduction Factor
%   MPSRF    - Multivariate Potential Scale Reduction Factor
%   CMPSRF   - Cumulative Multivariate Potential Scale Reduction Factor
%   IPSRF    - Interval-based Potential Scale Reduction Factor
%   CIPSRF   - Cumulative Interval-based Potential Scale Reduction Factor
%   KSSTAT   - Kolmogorov-Smirnov goodness-of-fit hypothesis test
%   HAIR     - Brooks' hairiness convergence diagnostic
%   CUSUM    - Yu-Mykland convergence diagnostic for MCMC
%   SCORE    - Calculate score-function convergence diagnostic
%   GBINIT   - Initial iterations for Gibbs iteration diagnostic
%   GBITER   - Estimate number of additional Gibbs iterations
%
% Time series analysis
%   ACORR      - Estimate autocorrelation function of time series
%   ACORRTIME  - Estimate autocorrelation evolution of time series (simple)
%   GEYER_ICSE - Compute autocorrelation time tau using Geyer's
%                initial convex sequence estimator
%                (requires Optimization toolbox)
%   GEYER_IMSE - Compute autocorrelation time tau using Geyer's
%                initial monotone sequence estimator
%
% Kernel density estimation etc.:
%   KERNEL1  - 1D Kernel density estimation of data
%   KERNELS  - Kernel density estimation of independent components of data
%   KERNELP  - 1D Kernel density estimation, with automatic kernel width
%   NDHIST   - Normalized histogram of N-dimensional data
%   HPDI     - Estimates the Bayesian HPD intervals
%
% Manipulation of MCMC chains
%   THIN     - Delete burn-in and thin MCMC-chains
%   JOIN     - Join similar structures of arrays to one structure of arrays
%   BATCH    - Batch MCMC sample chain and evaluate mean/median of batches
%
% Misc:
%   CUSTATS   - Calculate cumulative statistics of data
%   BBPRCTILE - Bayesian bootstrap percentile
%   GRADCHEK  - Checks a user-defined gradient function using finite
%               differences.
%   DERIVATIVECHECK - Compare user-supplied derivatives to
%                     finite-differencing derivatives.
%

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [tSeries, trends] = detrendTSeries(tSeries,detrendOption,smoothFrames)
%
% [detrendedTSeries,  trends] = detrendTSeries(tSeries,detrendOption,[smoothFrames])
%
% detrendOption is one of the following:
%   0 no trend removal
%   1 highpass trend removal
%   2 quadratic removal
%   -1 linear trend removal
% default determined by calling 'detrendFlag' that uses blockedAnalysisParams.detrend
%
% smoothFrames only needed for detrendOption==1
%
% djh, 2/2001
% jw, 7/2012: For linear and quadratic detrending, set model to have a
%             maximum value of 1 to avoid rank deficient matrices when
%             calculating
%                   wgts = model\tSeries;

%disp('Detrending tSeries...');

% reshape into matrix if tSeries is 3D array
dims = size(tSeries);
nFrames = dims(1);
tSeries = reshape(tSeries, nFrames, []);

switch detrendOption
case 2
    % remove a quadratic function

    model = [(1:nFrames).*(1:nFrames);(1:nFrames);ones(1,nFrames)]';

    % Limit range of model to [0 1] to avoid the possibility of rank
    % deficiency in calculating model \ tSeries
    model = bsxfun(@rdivide, model, max(model));
    wgts = model\tSeries;

    trends = model*wgts;
    tSeries = tSeries - trends;

case -1
    % remove a linear function
    model = [(1:nFrames);ones(1,nFrames)]';
    model = bsxfun(@rdivide, model, max(model));

    wgts = model\tSeries;

    trends = model*wgts;
    tSeries = tSeries - trends;

case 1
    % Do high-pass baseline removal
    calcstep = 1e4;
    if size(tSeries,2) <= calcstep;
        [tSeries, trends] = removeBaseline2(tSeries, smoothFrames);
    else % solve the out of memory problem
        trends = NaN(size(tSeries));
        for ii = 1:calcstep:size(tSeries,2);
            curRange = ii:min(ii+calcstep-1,size(tSeries,2));
            [tSeries(:,curRange), trends(:,curRange)] = ...
                removeBaseline2(tSeries(:,curRange), smoothFrames);
        end
    end
otherwise
    % Do nothing
    trends = zeros(size(tSeries));

end

% reshape tSeries in case we changed it from 3D array to matrix
tSeries = reshape(tSeries, dims);
trends    = reshape(trends, dims);

return

%% DEBUG

% seed random stream
s = RandStream('mt19937ar','Seed',1);
RandStream.setGlobalStream(s);

% generate a time series
ts = single(smooth(randn(300,1), 5));

% detrend three ways
[~, fit1] = detrendTSeries(ts,-1); % linear
[~, fit2] = detrendTSeries(ts, 2); % quadratic
[~, fit3] = detrendTSeries(ts, 1, 20); % high pass

% plot
figure(101)
clf

plot(1:length(ts), ts, 'k-');
hold on
plot(1:length(ts), fit1, 'r-', 'LineWidth', 2);
plot(1:length(ts), fit2, 'g-', 'LineWidth', 2);
plot(1:length(ts), fit3, 'b-', 'LineWidth', 2);

legend({'linear', 'quadratic', 'highpass'}, 'Location', 'Best')

% Compare quadratic detrend using old vs new method

% old:
n = length(ts);
model = [(1:n).*(1:n);(1:n);ones(1,n)]';
wgts = model\ts;
trends = model*wgts;

% new (divide each column of model by its max)
n = length(ts);
model = [(1:n).*(1:n);(1:n);ones(1,n)]';
model = bsxfun(@rdivide, model, max(model));
wgts = model\ts;
fit2 = model*wgts;

figure(102);
plot(1:n, trends, 'r', 1:n, fit2, 'k')


text:
function [P_MPP,U_MPP,I_MPP,U_OC_New,I_SC_New,U_PV,I_PV,P_PV] = PV_Module_Simulator_Function(G,T_amb)
%------ Solar Module simulation according to DIN EN 50530:2010 ------------
% This function can be used to simulate an solar cell to a PV-field by
% entering the number of modules (cells) in parallel and the number of
% module in serise.
% The equations used in this function are taken from the DIN EN 50530:2010
% norm.
% This function need two inputs:
% 1- G: The new solar irradiation at which the new working points
% (U_MPP,P_MPP,...) will be calculated.
% 2- T_amb: The new ambiant temperature at which the new working points
% (U_MPP,P_MPP,...) will be calculated.
% This function outputs are:
% 1- P_MPP,U_MPP,I_MPP: are the power, voltage, current at the maximum
% power point at the new conditions respectivily
% 2- U_OC_New,I_SC_New: are the open circuit voltage and the short circuit
% current respectivily calculated for the new weather(input) conditions.
% 3- U_PV,I_PV,P_PV: are the new characteristic curves values for the
% studied input condition.
% This function will also ask the user for addtional informations which
% need to be given in order to be able to simulate the characteristics of
% the selected PV-mdouel.
% Those information cal be obtained from the datasheet of the selected PV
% modules. if not given use the default values already given.
%Author: Aubai Alkhatib, date: 25.09.2013

prompt1 = {'Please Enter The STC Temperatur in C','Please Enter The STC Open Voltage in V','Please Enter The STC Short circut Current in A','Please Enter The Voltage temperatur coefficient Beta in %/C','Please Enter The Curent temperatur coefficient Alpha in %/C ','Please Enter The STC Irradiantion in W/m2','Please Eneter The Voltage of the MPP in V','Please Eneter The Current of the MPP in A','Please Eneter The Technology correction factor CG in W/m2','Please Eneter The Technology correction factor CU in pu','Please Eneter The Technology correction factor CR in m2/W','Please Eneter The number of PV-Modules in Parallel P','Please Eneter The number of PV-Modules in Seris R','Please Enter The Open circut Voltage of the LS'};
name = 'Initial Input data';
numlines = 1;
%defaultanswer1 = {'25','36.7','8.18','-0.32','0.04','1000','29.9','7.53','2.514','8.593','1.088','22','384','1200','Oldenburg,Germany','20130701','X:\Loc_4631\2013\07'};
defaultanswer1 = {'25','36.7','8.18','-0.32','0.04','1000','29.9','7.53','2.514','8.593','1.088','22','450','1600'};
%defaultanswer1 = {'4','1700','105','20','400','50','1'};
answer1 = inputdlg(prompt1,name,numlines,defaultanswer1);
T_STC = str2num(answer1{1})+273;
Uoc_STC =  str2num(answer1{2});
Isc_STC =  str2num(answer1{3});
V_T_C = str2num(answer1{4})/100;
I_T_C = str2num(answer1{5})/100;
E0 = str2num(answer1{6});
Umpp_STC = str2num(answer1{7});
Impp_STC = str2num(answer1{8});
CG = str2num(answer1{9})/1000;
CU = str2num(answer1{10})/100;
CR = str2num(answer1{11})/10000;
R = str2num(answer1{12});
P = str2num(answer1{13});
Ts = str2num(answer1{14});%#ok<*ST2NM>
G_STC = E0;
Uoc_STC = Uoc_STC * R;
Isc_STC = Isc_STC * P;
Umpp_STC = Umpp_STC * R;
Impp_STC = Impp_STC * P;
if 0 == 1
    T_0 = -3;
    k = 0.03;%km^2/w
    Tau = 5;% in min
    T_PV = T_amb + T_0 + ((k*1000000)/(1+(Tau*60)))*G;
else
    T_NOCT = 45;
    T_PV = T_amb + ((T_NOCT-20) * (G/800));
end
Alpha = I_T_C;
Isc_New = (Isc_STC*(G/G_STC)) * (1 + (Alpha * (T_PV - T_STC)));
Beta = V_T_C;
Uoc_New = Uoc_STC*(1+Beta*(T_PV-T_STC))*(log((G/CG)+1)*CU-CR*G);
FF_U = Umpp_STC/Uoc_STC;
FF_I = Impp_STC/Isc_STC;
I_0 = (Isc_STC*(1-FF_I)^(1/(1-FF_U)))*(G/G_STC);
CAQ = (FF_U-1)/(log(1-FF_I));
U_PV = 0:1:Ts;
I_PV = Isc_New - I_0*(exp(U_PV/(Uoc_New*CAQ))-1);
P_PV = U_PV.*I_PV;
P_MPP = max(P_PV);
U_MPP = U_PV(P_PV == P_MPP);
I_MPP = I_PV(U_MPP);
indx = find(I_PV < 0);
U_OC_New = U_PV(indx(1));
I_SC_New = I_PV(1);
figure(1);
subplot(2,1,1);
plot(U_PV,I_PV);xlabel('DC-Voltage in V');ylabel('DC-Current in A');title('PV-Module Charactaristic curves');grid on;ylim([0,I_MPP+((20/100)*I_MPP)]);xlim([0,U_OC_New+10]);
subplot (2,1,2);
plot(U_PV,P_PV);xlabel('DC-Voltage in V');ylabel('DC-Power in W');grid on;ylim([0,P_MPP+((5/100)*P_MPP)]);xlim([0,U_OC_New+10]);
end

author: Sable

repo: mcbench-benchmarks


text:
function savefigs(imageResultRepository, imageFileType)
% This function allows you to quickly save all currently open figures with
% a custom filename for each in multiple formats.  To use the function
% simply call savefigs with no arguments, then follow the prompts
%
% Upon execution this function will one-by-one bring each currently open
% figure to the foreground.  Then it will supply a text prompt in the main
% console window asking you for a filename.  It will save that figure to
% that filename in the .fig, .emf, .png, and .eps formats.
%
% The formats that it saves in can be changed by commenting out or adding
% lines below.
%
% Copyright 2010 Matthew Guidry
% matt.guidry ATT gmail DOTT com  (Email reformatted for anti-spam)

mkdir(imageResultRepository); % ensure that the repository exist

hfigs = get(0, 'children');                          %Get list of figures

for m = 1:length(hfigs)
    figure(hfigs(m))                                %Bring Figure to foreground
    % filename = input('Filename? (0 to skip)\n', 's')%Prompt user
    filename = num2str(hfigs(m));
    filename = strcat(imageResultRepository, '\', filename);

    if strcmp(filename, '0')                        %Skip figure when user types 0
        continue
    else
        saveas(hfigs(m), [filename imageFileType]);
        %saveas(hfigs(m), [filename '.fig']) %Matlab .FIG file
        %saveas(hfigs(m), [filename '.emf']) %Windows Enhanced Meta-File (best for powerpoints)
        %saveas(hfigs(m), [filename '.png']) %Standard PNG graphics file (best for web)
        %eval(['print -depsc2 ' filename])   %Enhanced Postscript (Level 2 color) (Best for LaTeX documents)
    end
end


text:
%INCLINING21	Inclining Experiment data for Example 7.1.
% Data taken from Hansen (1985), `An analytical treatment of the accuracy of the
% results of the inclining experiment', Naval Engineers Journal, Vol. 97, No.4,
% May, pp. 97-115.
% Format is [ moment tangent ], initial units [ ft-tons - ]
% Companion file for Biran, A. (2003), Ship Hydrostatics and Stability,
% Oxford: Butterworth-Heinemann.

incldata = [
 3853.7	 0.01869
 3853.7	 0.01854
 3853.7	 0.01788
 2570.1	 0.01259
 2570.1  0.01255
 2570.1	 0.01214
 1286.8	 0.00648
 1286.8	 0.00641
 1286.8	 0.00615
    3.8	 0.00042
    3.8	 0.00047
    3.8	 0.00000
-3785.8	-0.01791
-3785.8 -0.01802
-3785.8 -0.01850
-2523.3 -0.01194
-2523.3 -0.01203
-2523.3 -0.01235
-1263.7 -0.00574
-1263.7 -0.00599
-1263.7 -0.00646
   -0.6 -0.00000
   -0.6 -0.00000
   -0.6 -0.00063 ];
format bank
% separate data and convert to SI units
moment  = (0.305/1.016)*incldata(:, 1);
tangent = incldata(:, 2);
plot(tangent, moment, 'k.'), grid
Hl = ylabel('Inclining moment, pd, tm');
set(Hl, 'FontSize', 14)
Hl = xlabel('Heel angle tangent, tan\theta');
set(Hl, 'FontSize', 14)

hold on
tmin = min(tangent);
tmax = max(tangent);
M    = sum(tangent.*moment)/sum(tangent.^2);
Mmin = M*tmin;
Mmax = M*tmax;
plot( [ tmin tmax ], [ Mmin Mmax ], 'k-')
Ht = text(-0.015, 1100, [ 'Average slope = ' num2str(M) ]);
set(Ht, 'FontSize', 14)

hold off


text:
classdef CEC2010_F6 < PROBLEM
% <single> <real> <constrained>
% CEC'2010 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% R. Mallipeddi and P. N. Suganthan, Problem definitions and evaluation
% criteria for the CEC 2010 competition on constrained real-parameter
% optimization, Nanyang Technological University, Singapore, 2010.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat;	% Rotation matrix
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2010.mat'),'Data');
            obj.O = Data{6}.O;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 30
                obj.D   = 10;
                obj.Mat = Data{6}.M_10;
            else
                obj.D   = 30;
                obj.Mat = Data{6}.M_30;
            end
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 600;
            obj.upper    = zeros(1,obj.D) + 600;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = max(Z,[],2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Y = (Z+483.6106)*obj.Mat - 483.6106;
            PopCon(:,1) = abs(mean((-Y.*sin(sqrt(abs(Y)))),2)) - 1e-4;
            PopCon(:,2) = abs(mean((-Y.*cos(0.5*sqrt(abs(Y)))),2)) - 1e-4;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function r = isinf(a)
%ISINF        Array of 1's for inf components
%
%   r = isinf(a)
%

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

  r = reshape(any(isinf(a.t)),a.size);

repo: OpenMAS


text:
function [A,B] = kmo(X)
%KMO Kaiser-Meyer-Olkin Measure of Sampling Adequacy.
% Factor analysis can be used as a guide to how coherently a set of variables
% relate to a hypothesized underlying dimension that they are all being used
% to measure. External validity analysis assesses whether the scale that has
% been constructed performs as theoretically expected in correlation with
% other variables to which it is expected to be related.
% There are some assumptions about the characteristics of factors that are
% extracted and defined that are unobserved common dimensions that may be
% listed to account for the correlations among observed variables. Sampling
% adequacy predicts if data are likely to factor well, based on correlation
% and partial correlation. Is used to assess which variables to drop from the
% model because they are too multicollinear.
% It has been suggested that inv(R) should be a near-diagonal matrix in order
% to successfully fit a factor analysis model.  To assess how close inv(R)
% is to a diagonal matrix, Kaiser (1970) proposed a measure of sampling
% adequacy, now called KMO (Kaiser-Meyer-Olkin) index. The common part, called
% the image of a variable, is defined as that part which is predictable by
% regressing each variable on all other variables.
% The anti-image is the specific part of the variable that cannot be predicted.
% Examining the anti-image of the correlation matrix. That is the negative of the
% partial correlations, partialling out all other variables.
% There is a KMO statistic for each individual variable and their sum is
% the overall statistic. If it is not > 0.6 drop the indicator variables with
% the lowest individual statistic value until the overall one rises above 0.6:
% factors which is meritorious. The diagonal elements on the Anti-image
% correlation matrix are the KMO individual statistics for each variable. A KMO
% index <= 0.5 indicates the correlation matrix is not suitable for factor
% analysis.
%
% Syntax: function kmo(X)
%
%     Input:
%          X - Input matrix can be a data matrix (size n-data x p-variables)
%     Output(s):
%            - Kaiser-Meyer-Olkin Index.
%            - Degree of Common Variance Report (shared by a set of variables
%              and thus assesses the degree to which they measure a common
%              underlying factor).
%        optional(s):
%            - Anti-image Covariance Matrix.
%            - Anti-image Correlation Matrix
%
%  Example: From the example given on the web page
%  http://www.ncl.ac.uk/iss/statistics/docs/factoranalysis.html
%  We are interested to calculate the Kaiser-Meyer-Olkin measure of sampling
%  adequacy in order to see if proceeds a satisfactory factor analysis to
%  investigate the reasons why customers buy a product such as a particular
%  brand of soft drink (e.g. coca cola). Several variables were identified
%  which influence customer to buy coca cola. Some of the variables identified
%  as being influential include availability of product (X1), cost of product
%  (X2), experience with product (X3), popularity of product (X4), prestige
%  attached to product (X5), quality of product (X6), quantity of product (X7),
%  and respectability of product (X8). From this, you designed a questionnaire
%  to solicit customers' view on a seven point scale, where 1 = not important
%  and 7 = very important. The results from your questionnaire are show on the
%  table below. Only the first twelve respondents (cases) are used in this
%  example.
%
%  Table 1: Customer survey
%  --------------------------------------------------
%     X1    X2    X3    X4    X5    X6    X7    X8
%  --------------------------------------------------
%      4     1     4     5     2     3     6     7
%      4     2     7     6     6     3     3     4
%      6     4     3     4     2     5     7     7
%      5     3     5     4     3     4     6     7
%      5     2     4     5     2     5     5     6
%      6     3     3     5     4     4     7     7
%      6     2     4     4     4     3     4     5
%      4     1     3     4     3     3     5     6
%      5     3     4     3     4     3     6     6
%      5     4     3     4     4     4     6     7
%      6     2     4     4     4     3     7     5
%      5     2     3     3     3     3     7     6
%  --------------------------------------------------
%
% Data matrix must be:
%  X=[4 1 4 5 2 3 6 7;4 2 7 6 6 3 3 4;6 4 3 4 2 5 7 7;5 3 5 4 3 4 6 7;
%  5 2 4 5 2 5 5 6;6 3 3 5 4 4 7 7;6 2 4 4 4 3 4 5;4 1 3 4 3 3 5 6;
%  5 3 4 3 4 3 6 6;5 4 3 4 4 4 6 7;6 2 4 4 4 3 7 5;5 2 3 3 3 3 7 6];
%
%  Calling on Matlab the function:
%            kmo(X)
%
%  Answer is:
%
%  Kaiser-Meyer-Olkin Measure of Sampling Adequacy: 0.4172
%  The KMO test yields a degree of common variance unacceptable (Don't Factor).
%
%
%  Created by A. Trujillo-Ortiz, R. Hernandez-Walls, A. Castro-Perez,
%             K. Barba-Rojo and A. Otero-Limon
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.mx
%
%  Copyright. October 10, 2006.
%
% To cite this file, this would be an appropriate format:
% Trujillo-Ortiz, A., R. Hernandez-Walls, A. Castro-Perez, K. Barba-Rojo
%   and A. Otero-Limon (2006). kmo:Kaiser-Meyer-Olkin Measure of Sampling
%   Adequacy. A MATLAB file. [WWW document]. URL http://www.mathworks.com/
%   matlabcentral/fileexchange/loadFile.do?objectId=12736
%
%  References:
%  Rencher, A. C. (2002), Methods of Multivariate Analysis. 2nd. ed.
%            New-Jersey:John Wiley & Sons. Chapter 13 (pp. 408-450).
%

error(nargchk(1,1,nargin));
msg = nargoutchk(1, 2, nargout);

X = corrcoef(X);

iX = inv(X);
S2 = diag(diag((iX.^-1)));
AIS = S2*iX*S2; %anti-image covariance matrix
IS = X+AIS-2*S2; %image covariance matrix
Dai = diag(diag(sqrt(AIS)));
IR = inv(Dai)*IS*inv(Dai); %image correlation matrix
AIR = inv(Dai)*AIS*inv(Dai); %anti-image correlation matrix
a = sum((AIR - diag(diag(AIR))).^2);
AA = sum(a);
b = sum((X - eye(size(X))).^2);
BB = sum(b);
MSA = b./(b+a); %measures of sampling adequacy
AIR = AIR-eye(size(AIR))+diag(MSA);
%Examine the anti-image of the correlation matrix. That is the negative of the partial correlations,
%partialling out all other variables.
N = BB;
D = AA+BB;
kmo = N/D;

disp(' ')
fprintf('Kaiser-Meyer-Olkin Measure of Sampling Adequacy: %3.4f\n', kmo);
if (kmo >= 0.00 && kmo < 0.50);
    disp('The KMO test yields a degree of common variance unacceptable (Don''t Factor).')
elseif (kmo >= 0.50 && kmo < 0.60);
    disp('The KMO test yields a degree of common variance miserable.')
elseif (kmo >= 0.60 && kmo < 0.70);
    disp('The KMO test yields a degree of common variance mediocre.')
elseif (kmo >= 0.70 && kmo < 0.80);
    disp('The KMO test yields a degree of common variance middling.')
elseif (kmo >= 0.80 && kmo < 0.90);
    disp('The KMO test yields a degree of common variance meritorious.')
else (kmo >= 0.90 && kmo <= 1.00);
    disp('The KMO test yields a degree of common variance marvelous.')
end

if nargout == 1;
    disp(' ')
    disp('A = Anti-image covariance matrix.');
    A = AIS;
elseif nargout > 1;
    disp(' ')
    disp('A = Anti-image covariance matrix.');
    A = AIS;
    disp('B = Anti-image correlation matrix.');
    B = AIR;
end

return,