text:
function [hc] = read_ctf_hc(filename);

% READ_CTF_HC reads the MEG headcoil marker positions from an ascii file
% and computes the coordinate transformation required to get from from
% dewar to head-coordinates
%
% the definition of head coordinates is according to CTF standard:
% - the origin is exactly between LPA and RPA
% - the positive x-axis goes throught NAS
% - the positive y-axis goes (approximately) through LPA
% - the positive z-axis goes up, orthogonal to the x- and y-axes
%
% hc = read_ctf_hc(filename)
%
% returns a structure with the following fields
%   hc.dewar.nas	marker positions relative to dewar
%   hc.dewar.lpa
%   hc.dewar.rpa
%   hc.head.nas		marker positions relative to head (measured)
%   hc.head.lpa
%   hc.head.rpa
%   hc.standard.nas	marker positions relative to head (expected)
%   hc.standard.lpa
%   hc.standard.rpa
% and
%   hc.affine		parameter for affine transformation (1x12)
%   hc.homogenous	homogenous transformation matrix (4x4, see warp3d)
%   hc.translation	translation vector (1x3)
%   hc.rotation		rotation matrix (3x3)
%
% Gradiometer positions can be transformed into head coordinates using the
% homogeneous transformation matrix, or using the affine parameters and
% the warp3d function from the WARPING toolbox
%
% WARNING: all returned positions are in mm, and the transformation matrices
% also should be applied on channel positions that are in mm.

% Copyright (C) 2002, Robert Oostenveld
%
global fb

hc.standard.nas = [0 0 0];
hc.standard.lpa = [0 0 0];
hc.standard.rpa = [0 0 0];
hc.dewar.nas    = [0 0 0];
hc.dewar.lpa    = [0 0 0];
hc.dewar.rpa    = [0 0 0];
hc.head.nas     = [0 0 0];
hc.head.lpa     = [0 0 0];
hc.head.rpa     = [0 0 0];

fid = fopen(filename, 'r');
if fid==-1
  error(sprintf('could not open file %s', filename));
end

fseek(fid, 0, 'bof');
line = [];

if fb
  fprintf('reading standard coil positions with respect to the dewar\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% standard coil positions with respect to the dewar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while ~strcmp(line, 'standard nasion coil position relative to dewar (cm):')
  line = fgetl(fid);
  if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.nas(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.nas(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.nas(3) = str2num(r(2:end));
% NOTE THAT THERE IS AN TYPING ERROR IN SOME CTF FILES WHICH I HAVE TO REPRODUCE HERE (staNdard)
while ~(strcmp(line, 'stadard left ear coil position relative to dewar (cm):') | ...
        strcmp(line, 'standard left ear coil position relative to dewar (cm):'))
  line = fgetl(fid);
  if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.lpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.lpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.lpa(3) = str2num(r(2:end));
while ~strcmp(line, 'standard right ear coil position relative to dewar (cm):')
  line = fgetl(fid);
  if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.rpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.rpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.standard.rpa(3) = str2num(r(2:end));

if fb
  fprintf('reading measured coil positions with respect to the dewar\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% measured coil positions with respect to the dewar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while ~strcmp(line, 'measured nasion coil position relative to dewar (cm):')
  line = fgetl(fid);
  if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.nas(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.nas(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.nas(3) = str2num(r(2:end));
while ~strcmp(line, 'measured left ear coil position relative to dewar (cm):')
  line = fgetl(fid);
  if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.lpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.lpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.lpa(3) = str2num(r(2:end));
while ~strcmp(line, 'measured right ear coil position relative to dewar (cm):')
  line = fgetl(fid);
  if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.rpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.rpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.dewar.rpa(3) = str2num(r(2:end));

if fb
  fprintf('reading measured coil positions with respect to the head\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% measured coil positions with respect to the head
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while ~strcmp(line, 'measured nasion coil position relative to head (cm):')
  line = fgetl(fid);
  if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.nas(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.nas(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.nas(3) = str2num(r(2:end));
while ~strcmp(line, 'measured left ear coil position relative to head (cm):')
  line = fgetl(fid);
  if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.lpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.lpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.lpa(3) = str2num(r(2:end));
while ~strcmp(line, 'measured right ear coil position relative to head (cm):')
  line = fgetl(fid);
  if ~ischar(line) & line==-1, error('premature end of file'), end
end
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.rpa(1) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.rpa(2) = str2num(r(2:end));
line = fgetl(fid); [t, r] = strtok(line, '='); hc.head.rpa(3) = str2num(r(2:end));

fclose(fid);

if fb
  fprintf('computing transformation matrix\n');
end

% change the marker positions from cm to mm
hc.standard.nas = 10*hc.standard.nas;
hc.standard.lpa = 10*hc.standard.lpa;
hc.standard.rpa = 10*hc.standard.rpa;
hc.dewar.nas    = 10*hc.dewar.nas;
hc.dewar.lpa    = 10*hc.dewar.lpa;
hc.dewar.rpa    = 10*hc.dewar.rpa;
hc.head.nas     = 10*hc.head.nas;
hc.head.lpa     = 10*hc.head.lpa;
hc.head.rpa     = 10*hc.head.rpa;

% compute the direction of the head coordinate axes in dewar coordinates
d_x = hc.dewar.nas - (hc.dewar.lpa + hc.dewar.rpa)/2;
d_z = cross(d_x, hc.dewar.lpa - hc.dewar.rpa);
d_y = cross(d_z, d_x);
d_x = d_x / norm(d_x);
d_y = d_y / norm(d_y);
d_z = d_z / norm(d_z);

% compute the translation and rotation which are neccessary to transform
% any given location from the dewar to the head coordinate system
hc.translation = -(hc.dewar.lpa + hc.dewar.rpa)/2;
hc.rotation = inv([d_x' d_y' d_z']);

% construct the homogenous coordinate transformation matrix
% from the subsequent translation and rotation
hc.homogenous = eye(4,4);
hc.homogenous(1:3,4)   = hc.rotation * hc.translation';
hc.homogenous(1:3,1:3) = hc.rotation;

% construct a vector with the 12 parameters for an affine warp
hc.affine = hc.homogenous';
hc.affine = hc.affine(1:12);

author: PatternRecognition

repo: OpenBMI


text:
% FUNCTION kdtree = kdtree_create(points)
%
% AUTHOR:     Steven Michael
%             (smichael@ll.mit.edu)
%
% DATE:       2/17/05
%
% DESCRIPTION:
%
%  This function creates a KD Tree from the given points
%  and outputs it in the abstract object "kdtree"
%  The "kdtree" object can then be used for range finding
%  and nearest neighbor searching.
%
% INPUTS:
%
%   points   :     A (npoints X ndim) array of points, where "npoints"
%                  is the number of points and "ndim" is the number
%                  of dimensions.  Note that the points, even if they
%                  are double precision, will be converted to single
%                  precision when the tree is populated.  This is for
%                  speed -- most kdtree search applications don't
%                  necessitate double precision data.
%
% OUTPUTS:
%
%   kdtree   :     The abstract KD Tree object.
%
%
% Example:
%
%    % Create a list of 1000 random points in 3d space
%    r = rand(1000,3);
%
%    % Create a tree from this list
%    tree = kdtree(r);
%
%    % Find the point closest to the origin
%    [pntidx,pntval] = kdtree_closestpoint(tree,[0 0 0]);
%
%    % Create a list "r2" of 100 random points in 3d space and
%    % find the points in "r" that are closest to each point in "r2"
%    [pntidx,pntval] = kdtree_closestpoint(tree,r2);
%
%    % Find all the points within the cube defined by "rng"
%    rng = [ [.45 .55]; [.45 .55]; ; [.45 .55] ];
%    pntidx = kdtree_range(tree,rng);

author: Sable

repo: mcbench-benchmarks


text:
function [dat, varargout] = proc_regressOutComponents(dat, s, varargin)
%PROC_REGRESSOUTCOMPONENTS - Artifact removal by regressing out components
%
% [DAT, W, A, S_HAT] = regress_out_component(DAT, S)
%
%
%Arguments:
% DAT - data structure of epoched data
% s   - time course of components to be removed, size(s) = [T, n_components]
%       n_components can be more than one
% OPT - struct or property/value list of optional properties:
%  .return  - determines what additional matrices are returned. Valid values are
%             'all' (filters, patterns, and estimated components are returned,
%             in that order'),
%             'none' (default, no other matrices are returned)
%
%Returns
% DAT               - Updated data structure. Note that this method reduces the rank of DAT.x
% W (optional)      - filter matrix of size(W) = [n_components, n_channels], which is used to
%                     estimate the components from the data.
% A (optional)      - pattern matrix of size(A) = [n_channels, n_components], which
%                     shows how each of the components in s project to the channels
%                     in X.
% s_hat (optional)  - the estimate of s, which is extracted from DAT.x using W
%
%
%Description:
% Uses regression to remove the time course(s) contained in s from the data
% contained in DAT. s could be the eye movements measured with EOG channels
% and DAT could contain the remaining EEG channels. In this case the EOG signal
% that is contained in DAT.x will be removed (as best as possible, assuming a
% stationary eye movement pattern).
%
%References:
% Parra, L. C., Spence, C. D., Gerson, A. D., & Sajda, P. (2005),
% "Recipes for the linear analysis of EEG. Neuroimage, 28(2), 326-341

% Author(s): Sven Daehne, Mihail Bogojeski

props= { 'return'    'none'      'CHAR'};

if nargin==0,
  dat = props; return
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

varargout={};
dat = misc_history(dat);
misc_checkType(dat, 'STRUCT(x clab)');
misc_checkType(s, 'DOUBLE[- -]');
X = dat.x;
misc_checkType(X, 'DOUBLE[- -]');

T = size(X,1);
if not(size(s,1) == T)
    error('X and s must have the same number of samples! size(X,1) = %d, size(s,1) = %d', size(X,1),size(s,1));
end

% regression filter and patterns for the components
Cx = X'*X;
Cs = s'*s;
W = Cx \ X' * s; % regression weights for s
A = Cx*W/Cs; % spatial patterns of s

% remove estimate of s from the data
s_hat = X * W;
X_new = X - s_hat*A';

dat.x = X_new;
if strcmp(opt.return, 'all')
  varargout{1} = W;
  varargout{2} = A;
  varargout{3} = s_hat;
end

author: bbci


text:
%  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 2 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, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function tfrrmsct
%TFRRMSCT Unit test for the function TFRRMSC.

%	O. Lemoine - April 1996.

N=64;

% Reality of the TFR
sig=noisecg(N);
[tfr,rtfr]=tfrrmsc(sig);
if sum(any(abs(imag(rtfr))>sqrt(eps)))~=0,
 error('tfrrmsc test 1 failed');
end

% Energy conservation
%sig=fmlin(N);
%[tfr,rtfr]=tfrrmsc(sig);
%Es=norm(sig)^2;
%Etfr=sum(mean(rtfr));
%if abs(Es-Etfr)>sqrt(eps),
% error('tfrrmsc test 2 failed');
%end

% Positivity
if any(any(rtfr<0)),
 error('tfrrsp test 3 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
[tfr,rtfr]=tfrrmsc(sig);
[ik,jk]=find(abs(rtfr)>sqrt(eps));
if any(jk~=t0)|any(ik'-(2:N)),
 error('tfrrmsc test 4 failed');
end;

% frequency localization
%f0=10;
%sig=fmconst(N+6,f0/N);
%[tfr rtfr]=tfrrmsc(sig,N/2+2,N);
%if any(find(rtfr>max(rtfr)/N)~=f0+1)|(abs(mean(rtfr)-1.0)>sqrt(eps)),
% error('tfrrmsc test 5 failed');
%end;

N=65;

% Reality of the TFR
sig=noisecg(N);
[tfr,rtfr]=tfrrmsc(sig);
if sum(any(abs(imag(rtfr))>sqrt(eps)))~=0,
 error('tfrrmsc test 6 failed');
end

% Energy conservation
%sig=fmlin(N);
%[tfr,rtfr]=tfrrmsc(sig);
%Es=norm(sig)^2;
%Etfr=sum(mean(rtfr));
%if abs(Es-Etfr)>sqrt(eps),
% error('tfrrmsc test 7 failed');
%end

% Positivity
if any(any(rtfr<0)),
 error('tfrrsp test 8 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
[tfr,rtfr]=tfrrmsc(sig);
[ik,jk]=find(abs(rtfr)>sqrt(eps));
if any(jk~=t0)|any(ik'-(2:N)),
 error('tfrrmsc test 9 failed');
end;

% frequency localization
%f0=10;
%sig=fmconst(N+6,f0/N);
%[tfr rtfr]=tfrrmsc(sig,round(N/2)+2,N);
%if any(find(rtfr>max(rtfr)/N)~=f0+1)|(abs(mean(rtfr)-1.0)>sqrt(eps)),
% error('tfrrmsc test 10 failed');
%end;

author: HeLiangHIT


text:
function [state,options,Aineq,bineq,Aeq,beq] = ...
    psocheckinitialpopulation(state,Aineq,bineq,Aeq,beq,LB,UB,nonlcon,...
    options)
% Checks initial population with respect to linear constraints. Requires
% optimization toolbox.

if exist('linprog','file') ~= 2
    msg = sprintf('Could not find a required function in Optimization ') ;
    msg = sprintf('%s Toolbox. Ignoring (non)linear constraints ',msg) ;
    msg = sprintf('%s for initial population distribution and',msg) ;
    warning('pso:linearconstraints:missingtoolbox',...
        '%s setting constraint behavior to ''soft''.',msg)
    options.ConstrBoundary = 'soft' ;
    return
end

% Some robustness for older versions of optimset
% vv = ver ;
% vernos = {vv.Version} ;
% if str2double(vernos{strcmpi({vv.Name},'MATLAB')}) < 7.7
%     state.LinprogOptions = optimset('Simplex','off',...
%         'LargeScale','off',...
%         'Display','off') ;
% else
    state.LinprogOptions = optimset('Simplex','off',...
        'LargeScale','off',...
        'Algorithm','active-set',...
        'Display','off') ;
% end

state.OutOfBounds = false(options.PopulationSize,1) ;

hw = waitbar(0,'Finding feasible initial positions...') ;
for i = 1:size(state.Population,1)
    if strcmpi(options.LinearConstr.type,'linearconstraints')
        if (~isempty([Aineq,bineq]) && any(Aineq*state.Population(i,:)' - bineq > options.TolCon)) ...
                || (~isempty([Aeq,beq]) && any(abs(Aeq*state.Population(i,:)' - beq) > ...
                options.TolCon))
            % Reposition the ith particle if it is outside of linear
            % constraint boundaries
            [newpoint,unused,exitflag] = ...
                linprog([],Aineq,bineq,Aeq,beq,LB,UB,...
                state.Population(i,:),...
                state.LinprogOptions) ;
            clear unused
            if exitflag == -2
                error('Problem is infeasible due to constraints')
            else
                state.Population(i,:) = reshape(newpoint,1,[]) ;
            end % if exitflag
        end
    elseif strcmpi(options.LinearConstr.type,'nonlinearconstraints')
        [c,ceq] = nonlcon(state.Population(i,:)) ;
        lineq = false ;
        if ~isempty(Aineq)
            lineq = Aineq*state.Population(i,:)' - bineq > ...
                    options.TolCon ;
        end

        % Check constraint boundary type
        if i==1 && ~isempty(ceq) && strcmpi(options.ConstrBoundary,'soft')
            msg = '''Soft'' boundaries don''t work with nonlinear' ;
            msg = sprintf('%s equality constraints.',msg) ;
            warning('pso:initpop:nonlcon',...
                '%s Changing options.ConstrBoundary to ''absorb''.',...
                msg)
            options.ConstrBoundary = 'absorb' ;
        end

        if isempty(ceq) && isempty(Aeq)
            % Keep trying random points within PopInitRange until we
            % find one that satisfies the nonlinear constraints.
            % This method is faster, but it depends on the user setting
            % appropriate values for PopInitRange. It doesn't work with
            % equality constraints, linear or nonlinear. If LB and UB
            % exist, PopInitRange should already have been set to match
            % them.
            while any(c > options.TolCon) || any(lineq)
                state.Population(i,:) = options.PopInitRange(1,:) + ...
                    rand(1,size(options.PopInitRange,2)) .* ...
                    (options.PopInitRange(2,:) - ...
                    options.PopInitRange(1,:)) ;
                c = nonlcon(state.Population(i,:)) ;
                if ~isempty(Aineq)
                    lineq = Aineq*state.Population(i,:)' - bineq > ...
                        options.TolCon ;
                end
            end % while any
        elseif any(c > options.TolCon) || any(abs(ceq) > options.TolCon)
%         if any(c > options.TolCon) || any(abs(ceq) > options.TolCon)
            % Random point rejection is much faster, and produces
            % a more uniform distribution of points within the bounded
            % region. However, using fmincon to find acceptable points will
            % allow simultaneous compliance with linear and boundary
            % constraints, as well as dealing well with nonlinear equality
            % constraints. Appropriate setting of PopInitRange will prevent
            % the initial points from accumulating near the boundaries.
            [newpoint,unused,exitflag] = ...
                fmincon(@void,state.Population(i,:),...
                Aineq,bineq,Aeq,beq,LB,UB,nonlcon,state.LinprogOptions) ;
            clear unused
            if exitflag == -2
                error('Problem is infeasible due to nonlinear constraints')
            else
                state.Population(i,:) = reshape(newpoint,1,[]) ;
            end % if exitflag
        end % if isempty
    end % if strcmpi
    waitbar(i/options.PopulationSize,hw)
end % for i
if ishandle(hw), close(hw), end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function obj = squeezenet_autonn_custom_fn(block, inputs, ~)
% SQUEEZENET_AUTONN_CUSTOM_FN autonn custom layer converter
%
% Copyright (C) 2017 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

  switch class(block)
    case 'dagnn.Permute'
      obj = Layer.create(@permute, {inputs{1}, block.order}) ;
    case 'dagnn.Flatten'
      obj = Layer.create(@vl_nnflatten, {inputs{1}, block.axis}) ;
    case 'dagnn.Reshape'
      obj = Layer.create(@vl_nnreshape, {inputs{1}, block.shape}) ;
  end

author: albanie

repo: convnet-burden


text:
function varargout = implies(varargin)
%IMPLIES Logical implication
%
% IMPLIES(X,Y) creates a mixed integer representation of
% the constraint X --> Y, i.e. Y is true if X is true.
%
% Syntax
%   F = implies(X,Y,tol)
%
% Input
%   X : binary SDPVAR variable or a constraint
%   Y : binary SDPVAR variable or a constraint
%  tol: Optional threshhold for defining zero (see NOTE)
%
% Output
%   F : SET object
%
% Examples
%
%  binvar X Y; F = implies(X,Y);
%  binvar X;sdpvar Y; F = [implies(X,Y>=5), -10 <= Y <= 100];
%  binvar X;Y=sdpvar(3,1); F = [implies(X,[sum(Y);Y(2)]>=[5;0]), -1<= Y <= 10];
%
% Note
%
%  All variables in the expressions have to be explicitly bounded somewhere
%  in the model (implicit constraints are not sufficients such as [Y <= Z,...,Z<= 10]
%
%  Using implies with X non-binary is highly sensitive numerically.
%  The problem comes from the definition of 0 in a floating-point
%  environment, and precision in the solver. To account for this,
%  the user can supply a third argument to define a dead-zone around
%  zero, i.e Implies(X<=0,Y) will be replaced with IMPLIES(X<=-tol,Y)
%  Note, you typically need to tweak this number for your
%  application/solver. By default, YALMIP uses tol = 0, which means you
%  easily can get garbage... A positive number means YALMIP is cautious in
%  terms of activating the condition, while a negative number means YALMIP
%  will be aggressive  in activating the condition.
%
%   See also @SDPVAR/AND, @SDPVAR/OR, IFF

% There are some cases to take care of...
%
% X --> Y     binary/binary :                     Implemented
% X --> Y     binary/(LP,equality,sdp)            Implemented
% X --> Y     (LP,equality,sdp)/binary            Not implemented
% X --> Y     (LP,equality,sdp)/(LP,equality,sdp) Not implemented

X = varargin{1};
Y = varargin{2};

if isempty(X)
    varargout{1} = [];
elseif isa(X,'double') && isequal(X,0)
    varargout{1} = [];
elseif isa(X,'double') && isequal(X,1)
    varargout{1} = Y;
end

switch class(X)

    case {'sdpvar','constraint','lmi'}

        if nargin > 2
            if ~isa(varargin{3},'double')
                error('Third argument in IMPLIES should be a numerical value (tolerance)');
            end
        end
        varargout{1} = setupMeta(lmi([]), mfilename,varargin{:});

    case 'char'
        varargout{1} = implies_internal(varargin{3:end});

    case 'logical'
        if length(X)==1
            if X
                varargout{1} = Y;
            else
                varargout{1} = [];
            end
        else
            if length(X) == length(Y)
                i = find(X);
                if isempty(i)
                    varargout{1} = [];
                else
                    varargout{1} = Y(i);
                end
            else
                error('Size mismatch in input arguments');
            end
        end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function verifyquadpts
%% VERIFYQUADPTS examples and verfication on quadrature rules.
%
%  error = verifyquadpts(n) computes the error of n-th order quadrature
%  rule in a triangle. This is an example on the usage of quadrature points
%  and verification of qudarture order for approximating integrals in a
%  triangle.
%
% See also quadpts
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

% a reference triangle
node = [0 0; 1 0; 0 1];
elem = [1 2 3];
area = 0.5;
err = zeros(9,2);
for n = 1:9
    % get quadrature points
    [lambda,weight] = quadpts(n);
    nQuad = size(lambda,1);
    t1 = 0;
    t2 = 0;
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        t1 = t1 + weight(p)*f1(pxy(1),pxy(2),n);
        t2 = t2 + weight(p)*f2(pxy(1),pxy(2));
    end
    t1 = t1*area;
    t2 = t2*area;
    err(n,1) = abs(t1 - 2/((n+1)*(n+2)));
    err(n,2) = abs(t2 - (sin(1) - cos(1)));
end
display('Table: Error of quadrature for two smooth functions')
colname = {'n','x^n + y^n','sin(x+y)'};
disptable(colname,(1:9)',[],err(:,1),'%0.5e',err(:,2),'%0.5e');
end

function z = f1(x,y,n)
z = x.^n + y.^n;
end

function z = f2(x,y)
z = sin(x+y);
end

%% Results
% Let T be the triangle formed by (0,0), (1,0), and (0,1).
%
% Error1 is for the integral
% $\int _{T} x^n + y^n \, dxdy$.
% It should be numerically exact.
%
% Error2 is for the integral
% $\int _{T} \sin(x+y) \, dxdy$.
% It decays as n increas.
%
% See the doc for qudrature rules in <matlab:ifemdoc('quadpts') quadpts>.

repo: ifem


text:
function [L, D] = ldlt_sytr(A)
%LDLT_SYTR  Block LDL^T factorization for a symmetric tridiagonal matrix.
%           [L, D] = LDLT_SYTR(A) factorizes A = L*D*L', where A is
%           Hermitian tridiagonal, L is unit lower triangular, and D is
%           block diagonal with 1x1 and 2x2 diagonal blocks.  It uses
%           Bunch's strategy for choosing the pivots.

%           References:
%           J. R. Bunch, Partial pivoting strategies for symmetric
%              matrices.  SIAM J. Numer. Anal., 11(3):521-528, 1974.
%           N. J. Higham, Accuracy and Stability of Numerical Algorithms,
%              Second edition, Society for Industrial and Applied Mathematics,
%              Philadelphia, PA, 2002; sec. 11.1.4.

n = length(A);
if norm( tril(A,-2), 1) | norm( triu(A,2), 1) | ~isequal(A,A')
   error('Matrix must be Hermitian tridiagonal.')
end

s = norm(A(:), inf);
a = (sqrt(5)-1)/2;
L = eye(n);
D = zeros(n);
k = 1;

while k < n

      if s*abs(A(k,k)) >= a*abs(A(k+1,k))^2

         % 1-by-1 pivot.
         D(k,k) = A(k,k);
         L(k+1,k) = A(k+1,k)/A(k,k);
         A(k+1,k+1) = A(k+1,k+1) - abs(A(k+1,k))^2/A(k,k);
         k = k+1;

      else

         % 2-by-2 pivot.
         E = A(k:k+1,k:k+1);
         D(k:k+1,k:k+1) = E;
         if k+2 <= n
            L(k+2:n,k:k+1) = A(k+2:n,k:k+1)/E;
            A(k+2,k+2) = A(k+2,k+2) - abs(A(k+2,k+1))^2*A(k,k)/det(E);
         end
         k = k+2;

      end

end
if k == n
   D(k,k) = A(k,k);
end

author: opencobra

repo: cobratoolbox


text:
% KDTREE_PRINT prints the structure of a kd-tree
%
% Added by Pablo Sala, November 25, 2008
%
% SYNTAX
% kdtree_print
%
% DESCRIPTION
% Displays the structure of a kd-tree
%

% Copyright (c) 2008 Andrea Tagliasacchi
% All Rights Reserved
% email: ata2@cs.sfu.ca
% $Revision: 1.0$  Created on: 2008/09/15

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function [tm]=diag(tt)
%Constructs diagonal TT-matrix from TT-tensor
%
%   [TM]=DIAG(TT) constructs diagonal TT-matrix from TT-tensor
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

d=tt.d;
n=tt.n;
r=tt.r;
cr=tt.core;
ps=tt.ps;
psm=cumsum([1;(n.^2).*r(1:d).*r(2:d+1)]);
crm=zeros(psm(d+1)-1,1);
if isa(cr, 'gpuArray')
  % In case the GPU data is single.
  crm=cast(crm, classUnderlying(cr));
  crm=gpuArray(crm);
end
for i=1:d
  cr1=cr(ps(i):ps(i+1)-1); cr1=reshape(cr1,[r(i),n(i),r(i+1)]);
  crm1=zeros(r(i),n(i),n(i),r(i+1));
  if isa(cr, 'gpuArray')
    % In case the GPU data is single.
    crm1=cast(crm1, classUnderlying(cr));
    crm1=gpuArray(crm1);
  end
  for s1=1:r(i)
    for s2=1:r(i+1)
        dd=cr1(s1,:,s2);
        dd=reshape(dd,[n(i),1]);
       crm1(s1,:,:,s2)=diag(dd);
    end
  end
  crm(psm(i):psm(i+1)-1)=crm1(:);
end
tt.ps=psm;
tt.core=crm;
tt.n=n.^2;
tm=tt_matrix;
tm.m=n;
tm.n=n;
tm.tt=tt;
return
end

author: oseledets

repo: TT-Toolbox


text:
function haar_test03 ( )

%*****************************************************************************80
%
%% HAAR_TEST03 tests HAAR_2D and HAAR_2D_INVERSE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HAAR_TEST03\n' );
  fprintf ( 1, '  HAAR_2D computes the Haar transform of an array.\n' );
  fprintf ( 1, '  HAAR_2D_INVERSE inverts the transform.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Apply this to a 128x128 matrix of 0 and 1 values\n' );
  fprintf ( 1, '  which is actually a bit map of the Sierpinski triangle.\n' );
%
%  Load data from a file.
%
  u = load ( 'sierpinski.txt' );

  [ m, n ] = size ( u );
%
%  Demonstrate successful inversion.
%
  r8mat_print_some ( m, n, u, 1, 1, 10, 10, '  Input array U:' );

  v = haar_2d ( u );

  r8mat_print_some ( m, n, v, 1, 1, 10, 10, '  Transformed array V:' );

  w = haar_2d_inverse ( v );

  r8mat_print_some ( m, n, w, 1, 1, 10, 10, '  Recovered array W:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function quadrule_test065 ( n )

%*****************************************************************************80
%
%% QUADRULE_TEST065 uses CHEBSHEV2_COMPUTE to integral over the semicircle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADRULE_TEST065\n' );
  fprintf ( 1, '  Approximate the integral of f(x,y) over the semicircle\n' );
  fprintf ( 1, '    -1 <= x <= 1, y = sqrt ( 1 - x^2 )\n' );
  fprintf ( 1, '  using N Chebyshev points.\n' );
  fprintf ( 1, '  If p(x,y) involves any term of odd degree in y,\n' );
  fprintf ( 1, '  the estimate will only be approximate.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Polynomial    N    Integral        Estimate       Error\n' );
  fprintf ( 1, '\n' );

  [ x, w ] = chebyshev2_compute ( n );
%
%  f(x,y) = 1
%
  exact = 1.5707963267948966192;
  f(1:n,1) = 1.0;
  q = w(1:n)' * f(1:n);
  error = abs ( q - exact );
  fprintf ( 1, '  1            %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x
%
  exact = 0.0;
  f(1:n,1) = x(1:n);
  q = w(1:n)' * f(1:n);
  error = abs ( q - exact );
  fprintf ( 1, '  x            %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = y = sqrt ( 1 - x^2 )
%
  exact = 0.66666666666666666667;
  f(1:n,1) = sqrt ( 1.0 - x(1:n).^2 );
  q = w(1:n)' * f(1:n) / 2.0;
  error = abs ( q - exact );
  fprintf ( 1, '     y         %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x^2
%
  exact = 0.39269908169872415481;
  f(1:n,1) = x(1:n).^2;
  q = w(1:n)' * f(1:n);
  error = abs ( q - exact );
  fprintf ( 1, '  x^2          %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = xy = x * sqrt ( 1 - x^2 )
%
  exact = 0.0;
  f(1:n,1) = x(1:n) .* sqrt ( 1.0 - x(1:n).^2 );
  q = w(1:n)' * f(1:n) / 2.0;
  error = abs ( q - exact );
  fprintf ( 1, '  x  y         %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = y^2 -> ( 1 - x^2 )
%
  exact = 0.39269908169872415481;
  f(1:n,1) = 1.0 - x(1:n).^2;
  q = w(1:n)' * f(1:n) / 3.0;
  error = abs ( q - exact );
  fprintf ( 1, '     y^2       %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x^3
%
  exact = 0.0;
  f(1:n,1) = x(1:n).^3;
  q = w(1:n)' * f(1:n);
  error = abs ( q - exact );
  fprintf ( 1, '  x^3          %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x^2 y = x^2 sqrt ( 1 - x^2 )
%
  exact = 0.13333333333333333333;
  f(1:n,1) = x(1:n).^2 .* sqrt ( 1.0 - x(1:n).^2 );
  q = w(1:n)' * f(1:n) / 2.0;
  error = abs ( q - exact );
  fprintf ( 1, '  x^2y         %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x y^2 = x * ( 1 - x^2 )
%
  exact = 0.0;
  f(1:n,1) = x(1:n) .* ( 1.0 - x(1:n).^2 );
  q = w(1:n)' * f(1:n) / 3.0;
  error = abs ( q - exact );
  fprintf ( 1, '  x  y^2       %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = y^3
%
  exact = 0.26666666666666666667;
  f(1:n,1) = ( 1.0 - x(1:n).^2 ).^(1.5);
  q = w(1:n)' * f(1:n) / 4.0;
  error = abs ( q - exact );
  fprintf ( 1, '     y^3       %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x^4
%
  exact = 0.19634954084936207740;
  f(1:n,1) = x(1:n).^4;
  q = w(1:n)' * f(1:n);
  error = abs ( q - exact );
  fprintf ( 1, '  x^4          %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x^2y^2 -> x^2( 1 - x^2 )
%
  exact = 0.065449846949787359135;
  f(1:n,1) = x(1:n).^2 .* ( 1.0 - x(1:n).^2 );
  q = w(1:n)' * f(1:n) / 3.0;
  error = abs ( q - exact );
  fprintf ( 1, '  x^2y^2       %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = y^4 -> ( 1 - x^2 )^2
%
  exact = 0.19634954084936207740;
  f(1:n,1) = ( 1.0 - x(1:n).^2 ).^2;
  q = w(1:n)' * f(1:n) / 5.0;
  error = abs ( q - exact );
  fprintf ( 1, '     y^4       %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x^4y = x^4 sqrt ( 1 - x^2 )
%
  exact = 0.057142857142857142857;
  f(1:n,1) = x(1:n).^4 .* sqrt ( 1.0 - x(1:n).^2 );
  q = w(1:n)' * f(1:n) / 2.0;
  error = abs ( q - exact );
  fprintf ( 1, '  x^4y         %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  x^2y^3 = x^2 ( 1 - x^2 )**(3/2)
%
  exact = 0.038095238095238095238;
  f(1:n,1) = x(1:n).^2 .* ( 1.0 - x(1:n).^2 ).^(1.5);
  q = w(1:n)' * f(1:n) / 4.0;
  error = abs ( q - exact );
  fprintf ( 1, '  x^2y^3       %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = y^5
%
  exact = 0.15238095238095238095;
  f(1:n,1) = ( 1.0 - x(1:n).^2 ).^(2.5);
  q = w(1:n)' * f(1:n) / 6.0;
  error = abs ( q - exact );
  fprintf ( 1, '     y^5       %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x^6
%
  exact = 0.12271846303085129838;
  f(1:n,1) = x(1:n).^6;
  q = w(1:n)' * f(1:n);
  error = abs ( q - exact );
  fprintf ( 1, '  x^6          %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x^4y^2 -> x^4( 1 - x^2 )
%
  exact = 0.024543692606170259675;
  f(1:n,1) = x(1:n).^4 .* ( 1.0 - x(1:n).^2 );
  q = w(1:n)' * f(1:n) / 3.0;
  error = abs ( q - exact );
  fprintf ( 1, '  x^4y^2       %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = x^2y^4 -> x^2( 1 - x^2 )^2
%
  exact = 0.024543692606170259675;
  f(1:n,1) = x(1:n).^2 .* ( 1.0 - x(1:n).^2 ).^2;
  q = w(1:n)' * f(1:n) / 5.0;
  error = abs ( q - exact );
  fprintf ( 1, '  x^2y^4       %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );
%
%  f(x,y) = y^6 -> ( 1 - x^2 )^3
%
  exact = 0.12271846303085129838;
  f(1:n,1) = ( 1.0 - x(1:n).^2 ).^3;
  q = w(1:n)' * f(1:n) / 7.0;
  error = abs ( q - exact );
  fprintf ( 1, '     y^6       %2d  %14.6g  %14.6g  %14.6g\n', n, exact, q, error );

  return
end


text:
function t_end
%T_END  Finish running tests and print statistics.
%   T_END checks the global counters that were updated by calls to
%   T_OK, T_IS and T_SKIP and prints out a summary of the test results.
%
%   Example:
%       quiet = 0;
%       t_begin(5, quiet);
%       t_ok(pi > 3, 'size of pi');
%       t_skip(3, 'not yet written');
%       t_is(2+2, 4, 12, '2+2 still equals 4');
%       t_end;
%
%   See also T_BEGIN, T_OK, T_IS, T_SKIP, T_RUN_TESTS.

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

global t_quiet;
global t_num_of_tests;
global t_counter;
global t_ok_cnt;
global t_not_ok_cnt;
global t_skip_cnt;
global t_clock;

t_counter = t_counter - 1;

if t_counter == t_num_of_tests && ...
        t_counter == t_ok_cnt + t_skip_cnt && ...
        t_not_ok_cnt == 0
    all_ok = 1;
else
    all_ok = 0;
end

if t_quiet
    if all_ok
        fprintf('ok');
        if t_skip_cnt
            fprintf(' (%d of %d skipped)', t_skip_cnt, t_num_of_tests);
        end
    else
        fprintf('not ok\n');
        fprintf('\t#####  Ran %d of %d tests: %d passed, %d failed', ...
            t_counter, t_num_of_tests, t_ok_cnt, t_not_ok_cnt);
        if t_skip_cnt
            fprintf(', %d skipped', t_skip_cnt);
        end
    end
    fprintf('\n');
else
    if all_ok
        if t_skip_cnt
            fprintf('All tests successful (%d passed, %d skipped of %d)', ...
                t_ok_cnt, t_skip_cnt, t_num_of_tests);
        else
            fprintf('All tests successful (%d of %d)', t_ok_cnt, t_num_of_tests);
        end
    else
        fprintf('Ran %d of %d tests: %d passed, %d failed', ...
            t_counter, t_num_of_tests, t_ok_cnt, t_not_ok_cnt);
        if t_skip_cnt
            fprintf(', %d skipped', t_skip_cnt);
        end
    end
    fprintf('\nElapsed time %.2f seconds.\n', toc(t_clock));
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.
function xyz = xy2xyz(xy, Y)

    x = xy(:,1);
    y = xy(:,2);
    iy = 1 ./ (y + 1e-6);

    X = x .* Y .* iy;
    Z = Y .* (1-x-y) .* iy;

    Y = ones(X) * Y;

    xyz = [X Y Z];

author: petercorke

repo: machinevision-toolbox-matlab


text:
function [loss, fchg, tchg, dloss_dV, dchg_dVm] = get_losses(baseMVA, bus, branch)
%GET_LOSSES   Returns series losses (and reactive injections) per branch.
%
%   LOSS = GET_LOSSES(RESULTS)
%   LOSS = GET_LOSSES(BASEMVA, BUS, BRANCH)
%
%   [LOSS, CHG] = GET_LOSSES(RESULTS)
%   [LOSS, FCHG, TCHG] = GET_LOSSES(RESULTS)
%   [LOSS, FCHG, TCHG, DLOSS_DV] = GET_LOSSES(RESULTS)
%   [LOSS, FCHG, TCHG, DLOSS_DV, DCHG_DVM] = GET_LOSSES(RESULTS)
%
%   Computes branch series losses, and optionally reactive injections from
%   line charging, as functions of bus voltages and branch parameters, using the
%   following formulae:
%
%       loss = abs( Vf / tau - Vt ) ^ 2 / (Rs - j Xs)
%       fchg = abs( Vf / tau ) ^ 2 * Bc / 2
%       tchg = abs( Vt ) ^ 2 * Bc / 2
%
%   Optionally, computes the partial derivatives of the line losses with
%   respect to voltage angles and magnitudes.
%
%   Input:
%       RESULTS - a MATPOWER case struct with bus voltages corresponding to
%                 a valid power flow solution.
%                 (Can optionally be specified as individual fields BASEMVA,
%                  BUS, and BRANCH.)
%
%   Output(s):
%       LOSS - complex NL x 1 vector of losses (in MW), where NL is the number
%              of branches in the system, representing only the losses in the
%              series impedance element of the PI model for each branch.
%       CHG -  NL x 1 vector of total reactive injection for each line
%              (in MVAr), representing the line charging injections of both
%              of the shunt elements of PI model for each branch.
%       FCHG - Same as CHG, but for the element at the "from" end of the
%              branch only.
%       TCHG - Same as CHG, but for the element at the "to" end of the branch.
%       DLOSS_DV - Struct with partial derivatives of LOSS with respect to bus
%              voltages, with fields:
%           .a  - Partial with respect to bus voltage angles.
%           .m  - Partial with respect to bus voltage magnitudes.
%       DCHG_DVM - Struct with partial derivatives of FCHG and TCHG with
%              respect to bus voltage magnitudes, with fields:
%           .f  - Partial of FCHG with respect to bus voltage magnitudes.
%           .t  - Partial of TCHG with respect to bus voltage magnitudes.
%
%   Example:
%       results = runpf(mycase);
%       [loss, chg] = get_losses(results);
%       total_system_real_losses = sum(real(loss));
%       total_system_reac_losses = sum(imag(loss)) - sum(chg);
%
%       [loss, fchg, tchg, dloss_dV] = get_losses(results);

%   MATPOWER
%   Copyright (c) 1996-2017, 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.

%% default arguments
if isstruct(baseMVA)
    mpc = baseMVA;
    [baseMVA, bus, branch] = deal(mpc.baseMVA, mpc.bus, mpc.branch);
end

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

%% create map of external bus numbers to bus indices
i2e = bus(:, BUS_I);
e2i = sparse(max(i2e), 1);
e2i(i2e) = (1:size(bus, 1))';
out = find(branch(:, BR_STATUS) == 0);          %% out-of-service branches

%% sizes of things
nb = size(bus, 1);      %% number of buses
nl = size(branch, 1);   %% number of branches

%% construct complex bus voltage vector
V = bus(:, VM) .* exp(1j * pi/180 * bus(:, VA));

%% parameters
Cf = sparse(1:nl, e2i(branch(:, F_BUS)), branch(:, BR_STATUS), nl, nb);
Ct = sparse(1:nl, e2i(branch(:, T_BUS)), branch(:, BR_STATUS), nl, nb);
tap = ones(nl, 1);                              %% default tap ratio = 1 for lines
xfmr = find(branch(:, TAP));                    %% indices of transformers
tap(xfmr) = branch(xfmr, TAP);                  %% include transformer tap ratios
tap = tap .* exp(1j*pi/180 * branch(:, SHIFT)); %% add phase shifters
A = spdiags(1 ./ tap, 0, nl, nl) * Cf - Ct;
Ysc = 1 ./ (branch(:, BR_R) - 1j * branch(:, BR_X));
Vdrop = A * V;      %% vector of voltage drop across series impedance element
loss = baseMVA * Ysc .* Vdrop .* conj(Vdrop);
% loss = baseMVA * abs(V(e2i(branch(:, F_BUS))) ./ tap - V(e2i(branch(:, T_BUS)))) .^ 2 ./ ...
%             (branch(:, BR_R) - 1j * branch(:, BR_X));
% loss(out) = 0;

if nargout > 1
    Vf = Cf * V;
    Vt = Ct * V;
    fchg = real(baseMVA / 2 * branch(:, BR_B) .* Vf .* conj(Vf) ./ (tap .* conj(tap)));
    tchg = real(baseMVA / 2 * branch(:, BR_B) .* Vt .* conj(Vt));
%     fchg = abs(V(e2i(branch(:, F_BUS))) ./ tap) .^ 2 .* branch(:, BR_B) * baseMVA / 2;
%     tchg = abs(V(e2i(branch(:, T_BUS)))       ) .^ 2 .* branch(:, BR_B) * baseMVA / 2;
    fchg(out) = 0;
    tchg(out) = 0;

    if nargout == 2
        fchg = fchg + tchg;
    end
end

if nargout > 3
    B = spdiags(A * V, 0, nl, nl) * conj(A) * spdiags(conj(V), 0, nb, nb);
    dYsc = spdiags(Ysc, 0, nl, nl);
    dloss_dV = struct(...
        'a', -1j * baseMVA * dYsc * (B - conj(B)), ...
        'm',       baseMVA * dYsc * (B + conj(B)) * spdiags(1 ./ abs(V), 0, nb, nb) ...
    );
    if nargout > 4
        Bc = spdiags(branch(:, BR_B), 0, nl, nl);
        tt = spdiags(1 ./ (tap .* conj(tap)), 0, nl, nl);
        dchg_dVm = struct(...
            'f', baseMVA * Bc * tt * spdiags(Cf * bus(:, VM), 0, nl, nl) * Cf, ...
            't', baseMVA * Bc      * spdiags(Ct * bus(:, VM), 0, nl, nl) * Ct ...
        );
    end
end


text:
function cvx_optval = det_inv( X, p )

%DET_INV   Internal cvx version.

narginchk(1,2);
n = size( X, 1 );
if ndims( X ) > 2, %#ok
    error( 'N-D arrays are not supported.' );
elseif size( X, 2 ) ~= n,
    error( 'Matrix must be square.' );
elseif nargin < 2,
    p = 1;
elseif ~isnumeric( p ) || ~isreal( p ) || numel( p ) ~=  1 || p <= 0,
    error( 'Second argument must be a positive scalar.' );
end

w = [ ones(n,1) ; p ];
if cvx_isconstant( X ),

    cvx_optval = cvx( det_inv( cvx_constant( X ), p ) );

elseif nnz( X ) <= n && nnz( diag( X ) ) == nnz( X ),

    y = [];
    cvx_begin
        epigraph variable y
        geo_mean( [ diag(X) ; y ], w ) >= 1; %#ok
    cvx_end

elseif isreal( X ),

	y = []; Z = [];
    cvx_begin
        epigraph variable y
        variable Z(n,n) lower_triangular
        D = diag( Z );
        [ diag( D ), Z' ; Z, X ] == semidefinite(2*n); %#ok
        geo_mean( [ D ; y ], [], w ) >= 1; %#ok
    cvx_end

else

	y = []; Z = [];
    cvx_begin
        epigraph variable y
        variable Z(n,n) lower_triangular complex
        D = diag( Z );
        [ diag( D ), Z' ; Z, X ] == hermitian_semidefinite(2*n); %#ok
        geo_mean( [ real( D ) ; y ], [], w ) >= 1; %#ok
    cvx_end

end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [newPsi, Stats] = sampleSplitMerge_NoQRev( Psi, data, algParams )
% Split Merge via Not-Quite-Valid Sequential Allocation
% Proposes new candidate configuration for feat matrix F and stateSeq z
%  and accepts or rejects via Metropolis-Hastings
% HOWEVER, this *does not* account for the probability of the reverse move
%  in the acceptance ratio... there is still an acceptance ratio, but
%  it is entirely based on the joint probability under the model,
%  not on the proposal probabilities.
% Proposal selects anchor objects and features to split/merge at random,
%   unless provided (see DEBUG NOTE below)
%INPUT:
%  Psi  : model config input
%  data : SeqData object
%  algParams : relevant params for this method
%    ** SM.featSelectDistr can be one of
%       'random' : choose features to merge at random
%       'splitBias+margLik' (default) : use smart criteria to
%                    make successful proposals
%    ** SM.doSeqUpdateThetaHatOnMerge
%       if 1, refine Theta after visiting every seq. during merge proposal
%          0 (default), just estimate Theta initially and stick with it
%OUTPUT:
%  newPsi : model configuration as result of MH proposal
%    equals a new configuration when proposal accepted
%    or the input config Psi when rejected
%  Stats  : struct indicating type of move and whether accepted/rejected
%ASSUMES:
%  Psi comes with correct sufficient statistics for observed data (Xstats)
%DEBUG NOTE:
%   Can provide input "Psi" with fields "anchorIDs" and "activeFeatIDs"
%     to indicatea specific proposal to try.
%   E.g. to merge features 1 and 5 using objects 101 and 202 as seeds, set
%    Psi.anchorIDs = [101 202];
%    Psi.activeFeatIDs = [1 5];  Note F(101,1) and F(202,5) must be ON.

if isfield( Psi, 'activeFeatIDs' )
    [anchorIDs, featIDs, qFWD] = sampleAnchorAndFeatIDsToSplitMerge( Psi, data, algParams );
else
    [anchorIDs, featIDs, qFWD] = sampleAnchorAndFeatIDsToSplitMerge( Psi, data, algParams );
    Psi.anchorIDs        = anchorIDs;
    Psi.activeFeatIDs = featIDs;
end

if featIDs(1) == featIDs(2)
    % =========================================== SPLIT
    moveDescrStr = 'ADD';
    [propPsi] = sampleSplitConfig( Psi, data, anchorIDs, featIDs(1), algParams );
    %[propPsi, logQ] = sampleSplitConfig( Psi, data, anchorIDs, featIDs(1), algParams );
    %[~, logQ_Rev]   = sampleMergeConfig( propPsi, data, anchorIDs, propPsi.activeFeatIDs, algParams, Psi );
else
    % =========================================== MERGE
    moveDescrStr = 'DEL';
    [propPsi] = sampleMergeConfig( Psi, data, anchorIDs, featIDs, algParams );

    %[propPsi, logQ] = sampleMergeConfig( Psi, data, anchorIDs, featIDs, algParams );
    %[~, logQ_Rev]   = sampleSplitConfig( propPsi, data, anchorIDs, propPsi.activeFeatIDs, algParams, Psi );
end

% Total up probabilities of FORWARD (Q) and REVERSE (Q_Rev) moves
[~, ~, qREV] = sampleAnchorAndFeatIDsToSplitMerge( propPsi, data, algParams );
logQ_Rev.all = log(qREV); % + logQ_Rev.F + logQ_Rev.z;
logQ.all     = log(qFWD); % + logQ.F + logQ.z;

% NB: not passing data as an arg here
%   means that we trust the stored X suff stats in Psi! Yay efficiency.
logPr_Cur  = calcJointLogPr_BPHMMState( Psi );
logPr_Prop = calcJointLogPr_BPHMMState( propPsi );

logPrAccept = logPr_Prop.all - logPr_Cur.all + logQ_Rev.all - logQ.all;
rho = exp( logPrAccept );
rho = min(1, rho);
doAccept = rand < rho;

if (  doAccept )
    newPsi = propPsi;
    % Remove empty columns of F, and rename state sequence appropriately
    newPsi = reallocateFeatIDs( newPsi );
else
    newPsi = Psi;
end

% Strip off info used internally to identify how to construct proposals
newPsi = rmfield( newPsi, 'anchorIDs');
newPsi = rmfield( newPsi, 'activeFeatIDs');

Stats.nAccept = doAccept;
Stats.rho = rho;
Stats.moveDescr = moveDescrStr;

end % main function

author: michaelchughes

repo: NPBayesHMM


text:
classdef StressNormSuperEllipseComputer < handle

    properties (Access = private)
        microProblem
        fName
        outputFolder
    end

    properties (Access = private)
        mx
        my
        q
        phi
        pNorm
        print
        fileName
        hasToCaptureImage
        testName
        iter
        meshBackground
        mesh
        hMesh
    end

    methods (Access = public)

        function obj = StressNormSuperEllipseComputer(cParams)
            obj.init(cParams);
        end

        function sPnorm = compute(obj)
            obj.createMicroProblem()
            sPnorm = obj.computePstressNorm();
        end

        function var = computeCellVariables(obj)
            obj.compute();
            mProblem = obj.microProblem;
            var.Ctensor = mProblem.variables.varFromCh.Chomog;
            var.tstress = mProblem.variables.varFromCh.tstress;
            var.tstrain = mProblem.variables.varFromCh.tstrain;
            var.displ   = mProblem.variables.varFromCh.tdisp;
            var.volume  = mProblem.mesh.computeVolume();
            var.mesh    = mProblem.mesh;

            quad = Quadrature.set(var.mesh.geometryType);
            quad.computeQuadrature('CONSTANT');
            volume = var.mesh.computeDvolume(quad);
            var.integrationVar.dV    = volume;
            var.integrationVar.nstre = size(var.tstress,1);
            var.integrationVar.ngaus = size(var.tstress,2);
            var.integrationVar.geoVol = obj.meshBackground.computeVolume();
        end

        function printImage(obj)
            outputName = [obj.fileName,'Print'];
            if obj.print
                    s.mesh       = obj.mesh;
                    s.outPutName = outputName;
                    printer = SuperEllipsePrinter(s);
                    printer.print();
                if obj.hasToCaptureImage
                    printer.captureImage()
                end
            end
        end

        function printStress(obj)
            if obj.print
                microP = obj.microProblem;
                outputName = [obj.fileName,'Print'];
                dI.mesh    =  microP.mesh;
                dI.outName = outputName;
                dI.pdim    = '2D';
                dI.ptype   = 'MICRO';
                ps = PostProcessDataBaseCreator(dI);
                dB = ps.getValue();
                postCase = 'ElasticityMicro';
                postProcess = Postprocess(postCase,dB);
                d.fields = microP.variables;
                d.quad   = microP.element.quadrature;
                postProcess.print(obj.iter,d);
            end
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mx       = cParams.mx;
            obj.my       = cParams.my;
            obj.q        = cParams.q;
            obj.phi      = cParams.phi;
            obj.pNorm    = cParams.pNorm;
            obj.print    = cParams.print;
            obj.fileName = cParams.fileName;
            obj.iter     = cParams.iter;
            obj.hMesh    = cParams.hMesh;
            obj.hasToCaptureImage = cParams.hasToCaptureImage;
            obj.computeFileNameAndOutputFolder();
        end

        function computeFileNameAndOutputFolder(obj)
            obj.fName = [obj.fileName];
            obj.outputFolder = fullfile(pwd,'Output',obj.fName);
        end

        function sPnorm2 = computePstressNorm(obj)
            stress = [sin(obj.phi) cos(obj.phi) 0];
            p = obj.pNorm;
            sPnorm2 = obj.microProblem.computeStressPnorm(stress,p);
        end

        function ls = createLevelSet(obj)
            sM.coord  = obj.meshBackground.coord;
            sM.connec = obj.meshBackground.connec;
            s.mesh = Mesh_Total(sM);

            s.widthH = obj.mx;
            s.widthV = obj.my;
            s.pnorm  = obj.q;
            s.type = 'smoothRectangle';

            s.levelSetCreatorSettings = s;
            s.type = 'LevelSet';

            s.scalarProductSettings.epsilon = 1;
            levelSet = LevelSet(s);
            ls = levelSet.value;
        end

        function createBackgroundMesh(obj)
            obj.testName = 'RVE_Square_Triangle_FineFine';
            %obj.testName = 'RVE_Square_Triangle_Fine';
            s.testName = obj.testName;
            obj.meshBackground = Mesh().createFromFile(s);
        end

        function createMesh(obj)
            s.fileName = obj.fileName;
            s.levelSet = obj.createLevelSet();
            s.meshBackground = obj.meshBackground;
            s.hMesh = obj.hMesh;
            mCreator = MeshCreatorFromLevelSetWithMMG(s);
            obj.mesh = mCreator.create();
        end

        function createMicroProblem(obj)
            obj.createBackgroundMesh();
            obj.createMesh();
            femSolver = Elastic_Problem_Micro.create(obj.testName);
            femSolver.setMesh(obj.mesh);
            props.kappa = .75;
            props.mu    = .375;
            femSolver.setMatProps(props);
            obj.microProblem = femSolver;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ quasi, seed ] = i4_sobol ( dim_num, seed )

%*****************************************************************************80
%
%% I4_SOBOL generates a new quasirandom Sobol vector with each call.
%
%  Discussion:
%
%    The routine adapts the ideas of Antonov and Saleev.
%
%    Thanks to Francis Dalaudier for pointing out that the range of allowed
%    values of DIM_NUM should start at 1, not 2!  17 February 2009.
%
%    This function was modified to use PERSISTENT variables rather than
%    GLOBAL variables, 13 December 2009.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2012
%
%  Author:
%
%    Original FORTRAN77 version by Bennett Fox.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Antonov, Saleev,
%    USSR Computational Mathematics and Mathematical Physics,
%    Volume 19, 1980, pages 252 - 256.
%
%    Paul Bratley, Bennett Fox,
%    Algorithm 659:
%    Implementing Sobol's Quasirandom Sequence Generator,
%    ACM Transactions on Mathematical Software,
%    Volume 14, Number 1, pages 88-100, 1988.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Ilya Sobol,
%    USSR Computational Mathematics and Mathematical Physics,
%    Volume 16, pages 236-242, 1977.
%
%    Ilya Sobol, Levitan,
%    The Production of Points Uniformly Distributed in a Multidimensional
%    Cube (in Russian),
%    Preprint IPM Akad. Nauk SSSR,
%    Number 40, Moscow 1976.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the number of spatial dimensions.
%    DIM_NUM must satisfy 1 <= DIM_NUM <= 40.
%
%    Input/output, integer SEED, the "seed" for the sequence.
%    This is essentially the index in the sequence of the quasirandom
%    value to be generated.  On output, SEED has been set to the
%    appropriate next value, usually simply SEED+1.
%    If SEED is less than 0 on input, it is treated as though it were 0.
%    An input value of 0 requests the first (0-th) element of the sequence.
%
%    Output, real QUASI(DIM_NUM), the next quasirandom vector.
%
  persistent atmost;
  persistent dim_max;
  persistent dim_num_save;
  persistent initialized;
  persistent lastq;
  persistent log_max;
  persistent maxcol;
  persistent poly;
  persistent recipd;
  persistent seed_save;
  persistent v;

  if ( isempty ( initialized ) )
    initialized = 0;
    dim_num_save = -1;
  end

  if ( ~initialized | dim_num ~= dim_num_save )

    initialized = 1;

    dim_max = 40;
    dim_num_save = -1;
    log_max = 30;
    seed_save = -1;
%
%  Initialize (part of) V.
%
    v(1:dim_max,1:log_max) = zeros(dim_max,log_max);

    v(1:40,1) = [ ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]';

    v(3:40,2) = [ ...
            1, 3, 1, 3, 1, 3, 3, 1, ...
      3, 1, 3, 1, 3, 1, 1, 3, 1, 3, ...
      1, 3, 1, 3, 3, 1, 3, 1, 3, 1, ...
      3, 1, 1, 3, 1, 3, 1, 3, 1, 3 ]';

    v(4:40,3) = [ ...
               7, 5, 1, 3, 3, 7, 5, ...
      5, 7, 7, 1, 3, 3, 7, 5, 1, 1, ...
      5, 3, 3, 1, 7, 5, 1, 3, 3, 7, ...
      5, 1, 1, 5, 7, 7, 5, 1, 3, 3 ]';

    v(6:40,4) = [ ...
                     1, 7, 9,13,11, ...
      1, 3, 7, 9, 5,13,13,11, 3,15, ...
      5, 3,15, 7, 9,13, 9, 1,11, 7, ...
      5,15, 1,15,11, 5, 3, 1, 7, 9 ]';

    v(8:40,5) = [ ...
                           9, 3,27, ...
     15,29,21,23,19,11,25, 7,13,17, ...
      1,25,29, 3,31,11, 5,23,27,19, ...
     21, 5, 1,17,13, 7,15, 9,31, 9 ]';

    v(14:40,6) = [ ...
              37,33, 7, 5,11,39,63, ...
     27,17,15,23,29, 3,21,13,31,25, ...
      9,49,33,19,29,11,19,27,15,25 ]';

    v(20:40,7) = [ ...
                                         13, ...
     33,115, 41, 79, 17, 29,119, 75, 73,105, ...
      7, 59, 65, 21,  3,113, 61, 89, 45,107 ]';

    v(38:40,8) = [ ...
                                7, 23, 39 ]';
%
%  Set POLY.
%
    poly(1:40)= [ ...
        1,   3,   7,  11,  13,  19,  25,  37,  59,  47, ...
       61,  55,  41,  67,  97,  91, 109, 103, 115, 131, ...
      193, 137, 145, 143, 241, 157, 185, 167, 229, 171, ...
      213, 191, 253, 203, 211, 239, 247, 285, 369, 299 ];

    atmost = 2^log_max - 1;
%
%  Find the number of bits in ATMOST.
%
    maxcol = i4_bit_hi1 ( atmost );
%
%  Initialize row 1 of V.
%
    v(1,1:maxcol) = 1;

  end
%
%  Things to do only if the dimension changed.
%
  if ( dim_num ~= dim_num_save )
%
%  Check parameters.
%
    if ( dim_num < 1 | dim_max < dim_num )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'I4_SOBOL - Fatal error!\n' );
      fprintf ( 1, '  The spatial dimension DIM_NUM should satisfy:\n' );
      fprintf ( 1, '    1 <= DIM_NUM <= %d\n', dim_max );
      fprintf ( 1, '  But this input value is DIM_NUM = %d\n', dim_num );
      return
    end

    dim_num_save = dim_num;
%
%  Initialize the remaining rows of V.
%
    for i = 2 : dim_num
%
%  The bits of the integer POLY(I) gives the form of polynomial I.
%
%  Find the degree of polynomial I from binary encoding.
%
      j = poly(i);
      m = 0;

      while ( 1 )

        j = floor ( j / 2 );

        if ( j <= 0 )
          break;
        end

        m = m + 1;

      end
%
%  Expand this bit pattern to separate components of the logical array INCLUD.
%
      j = poly(i);
      for k = m : -1 : 1
        j2 = floor ( j / 2 );
        includ(k) = ( j ~= 2 * j2 );
        j = j2;
      end
%
%  Calculate the remaining elements of row I as explained
%  in Bratley and Fox, section 2.
%
      for j = m + 1 : maxcol
        newv = v(i,j-m);
        l = 1;
        for k = 1 : m
          l = 2 * l;
          if ( includ(k) )
            newv = bitxor ( newv, l * v(i,j-k) );
          end
        end
        v(i,j) = newv;
      end
    end
%
%  Multiply columns of V by appropriate power of 2.
%
    l = 1;
    for j = maxcol-1 : -1 : 1
      l = 2 * l;
      v(1:dim_num,j) = v(1:dim_num,j) * l;
    end
%
%  RECIPD is 1/(common denominator of the elements in V).
%
    recipd = 1.0 / ( 2 * l );

    lastq(1:dim_num) = 0;

  end

  seed = floor ( seed );

  if ( seed < 0 )
    seed = 0;
  end

  if ( seed == 0 )

    l = 1;
    lastq(1:dim_num) = 0;

  elseif ( seed == seed_save + 1 )
%
%  Find the position of the right-hand zero in SEED.
%
    l = i4_bit_lo0 ( seed );

  elseif ( seed <= seed_save )

    seed_save = 0;
    l = 1;
    lastq(1:dim_num) = 0;

    for seed_temp = seed_save : seed - 1
      l = i4_bit_lo0 ( seed_temp );
      for i = 1 : dim_num
        lastq(i) = bitxor ( lastq(i), v(i,l) );
      end
    end

    l = i4_bit_lo0 ( seed );

  elseif ( seed_save + 1 < seed )

    for seed_temp = seed_save + 1 : seed - 1
      l = i4_bit_lo0 ( seed_temp );
      for i = 1 : dim_num
        lastq(i) = bitxor ( lastq(i), v(i,l) );
      end
    end

    l = i4_bit_lo0 ( seed );

  end
%
%  Check that the user is not calling too many times!
%
  if ( maxcol < l )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_SOBOL - Fatal error!\n' );
    fprintf ( 1, '  Too many calls!\n' );
    fprintf ( 1, '  MAXCOL = %d\n', maxcol );
    fprintf ( 1, '  L =      %d\n', l );
    return
  end
%
%  Calculate the new components of QUASI.
%
  for i = 1 : dim_num
    quasi(i) = lastq(i) * recipd;
    lastq(i) = bitxor ( lastq(i), v(i,l) );
  end

  seed_save = seed;
  seed = seed + 1;

  return
end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [ind,d] = find(lookup,ref,epsilon,varargin)
% return indece and distance of all nodes within a eps neighborhood
%
% Syntax
%   % find for each quaternion in ref the closest quaternion in table
%   [ind,d] = find(table, ref)
%
%   % find for each quaternion in ref all quaternion in lookup that have
%   % angle not larger then epsilon
%   [ind,d] = find(table, rot_ref, epsilon)
%
%   % find all quaternion which have distant to fibre f less then epsilon
%   [ind,d] = find(table, f_ref, epsilon)
%
% Input
%  table   - @quaternion
%  rot_ref - reference @quaternion
%  f_ref   - reference @fibre
%  epsilon - maximum distance
%
% Output
%  ind - index of the found @quaternions
%  d   - actual distances
%

if isa(ref,'fibre')

  ind = angle(ref,lookup) < epsilon;

else

  d = dot_outer(lookup,ref);

  if nargin == 2
    [d,ind] = max(d,[],1);
  else
    ind = d > cos(epsilon/2);
  end
end

author: mtex-toolbox

repo: mtex

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


text:
function timer_cputime_test02 ( )

%*****************************************************************************80
%
%% TEST02 times the vectorized EXP routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 December 2005
%
%  Author:
%
%    John Burkardt
%
  n_log_min = 12;
  n_log_max = 22;
  n_min = 2^n_log_min;
  n_max = 2^n_log_max;
  n_rep = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  Time vectorized operations:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    y(1:n) =        x(1:n)  \n' );
  fprintf ( 1, '    y(1:n) = PI *   x(1:n)  \n' );
  fprintf ( 1, '    y(1:n) = sqrt ( x(1:n) )\n' );
  fprintf ( 1, '    y(1:n) = exp  ( x(1:n) )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data vectors will be of minimum size %d\n', n_min );
  fprintf ( 1, '  Data vectors will be of maximum size %d\n', n_max );
  fprintf ( 1, '  Number of repetitions of the operation: %d\n', n_rep );

  for func = 1 : 4

    for i_rep = 1 : n_rep

      for n_log = n_log_min : n_log_max

        n = 2^n_log;

        x = rand(n,1);

        time1 = cputime;

        if ( func == 1 )
          y(1:n) = x(1:n);
        elseif ( func == 2 )
          y(1:n) = pi * x(1:n);
        elseif ( func == 3 )
          y(1:n) = sqrt ( x(1:n) );
        elseif ( func == 4 )
          y(1:n) = exp ( x(1:n) );
        end

        time2 = cputime;

        delta(n_log,i_rep) = time2 - time1;

      end

    end

    fprintf ( 1, '\n' );
    fprintf ( 1, 'Test02 Results:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Vector Size  Rep #1        Rep #2        ' );
    fprintf ( 1, 'Rep #3        Rep #4        Rep #5\n' );
    fprintf ( 1, '\n' );

    for n_log = n_log_min : n_log_max
      n = 2^n_log;
      fprintf ( 1, '%10d', n );
      for j = 1 : n_rep
        fprintf ( 1, '%14f', delta(n_log,1:n_rep) );
      end
      fprintf ( 1, '\n' );
    end

  end

  return
end


text:
function elem=delendelem(elem,mask)
%
% elem=delendelem(elem,mask)
%
% delete elements whose nodes are all edge nodes
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2007/11/24
%
% input/output:
%      elem: input/output, surface/volumetric element list
%      mask: of length of node number, =0 for internal nodes, =1 for edge nodes
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

badidx=sum(mask(elem)');
elem(find(badidx==size(elem,2)),:)=[];

repo: REST


text:
function segment_length = p11_boundary_segment_length ( segment_index, h )

%*****************************************************************************80
%
%% P11_BOUNDARY_SEGMENT_LENGTH returns boundary segment lengths in problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEGMENT_INDEX, the index of one of the boundary segments.
%
%    Input, real H, the suggested spacing between points.
%
%    Output, integer SEGMENT_LENGTH, the number of points in the segment.
%
  if ( h <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P11_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
    fprintf ( 1, '  Nonpositive H = %f\n', h );
    error ( 'P11_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
  end

  if ( segment_index == 1 )

    n = round ( 4.0 / h );
    n = max ( n, 5 );
    segment_length = n + mod ( 4 - mod ( n - 1, 4 ), 4 );

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P11_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
    fprintf ( 1, '  Illegal SEGMENT_INDEX = %d\n', segment_index );
    error ( 'P11_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );

  end

  return
end


text:
function fx = p46_fun ( n, x )

%*****************************************************************************80
%
%% P46_FUN evaluates the integrand for problem 46.
%
%  Discussion:
%
%    The problem has a parameter ALPHA that can be set by calling
%    P46_PARAM_SET.
%
%    The integrand is the radius of an ellipse as a function of angle.
%
%    The integral represents the arc length of the ellipse.
%
%    The suggested parameter range is 0.0 <= ALPHA < 1.0.  ALPHA is
%    the eccentricity of the ellipse.
%
%  Interval:
%
%    0 <= theta <= 2 pi
%
%  Integrand:
%
%    r(theta) = ( 1 - alpha^2 ) / ( 1 - alpha * cos ( theta ) )
%
%  Exact Integral:
%
%    When alpha = sin ( pi / 12 ), then
%
%      6.0690909595647754101
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Richard Crandall,
%    Projects in Scientific Computing,
%    Springer, 2000, page 47.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
%
  alpha = p46_param_get ( );

  fx(1:n) = ( 1.0 - alpha^2 ) ./ ( 1.0 - alpha * cos ( x(1:n) ) );

  return
end


text:
function [x,state] = struct_log(z,task)
%STRUCT_LOG Natural logarithm.
%   [x,state] = struct_log(z,[]) computes x as log(z). The structure state
%   stores information which is reused in computing the right and left
%   Jacobian-vector products.
%
%   struct_log(z,task) computes the right or left Jacobian-vector product
%   of this transformation, depending on the structure task. Use the
%   structure state and add the field 'r' of the same shape as z or the
%   field 'l' of the same shape as x to obtain the structure task for
%   computing the right and left Jacobian-vector products
%
%      (dF(:)/dz(:).')*task.r(:) and
%      (dF(:)/dz(:).')'*task.l(:) + conj((dF(:)/dconj(z(:)).')'*task.l(:)),
%
%   respectively. Here, F(z) represents this transormation, (:) signifies
%   vectorization and the derivative w.r.t. z (conj(z)) is a partial
%   derivative which treats conj(z) (z) as constant. The output has the
%   same shape as x or z for the right and left Jacobian-vector products,
%   respectively.
%
%   See also struct_power, struct_sqrt.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if nargin < 2, task = []; end
state = [];

if isempty(task) || (isempty(task.l) && isempty(task.r))
    x = log(z);
    state.deriv = 1./z;
elseif ~isempty(task.r)
    x = task.deriv.*task.r;
elseif ~isempty(task.l)
    x = conj(task.deriv).*task.l;
end

end

author: andrewssobral

repo: lrslibrary


text:
function res = indchannel(this, label)
% Method for getting channel indices based on channel labels
% FORMAT  res = indchannel(this, label)
% this       - MEEG object
% label      - string or cell array of labels
%
% res        - vector of channel indices matching labels
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: indchannel.m 3254 2009-07-07 15:18:54Z vladimir $

if ischar(label)
    label = {label};
end

[junk, res] = match_str(label, chanlabels(this));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clearvars; clearvars -global; clc; close all;

%% settings
addpath(genpath('Functions'));
addpath(genpath('./Libraries/'));

InitParam();

% fprintf('Preparing the training data\n');
% fprintf('***************************\n\n');
% PrepareTrainData();

fprintf('Preparing the test data\n');
fprintf('***********************\n\n');
PrepareTestData();

author: qingsenyangit

repo: AHDRNet


text:
function [idxStart, idxEnd] = find_ts_idx(ts, tStart, tEnd)
    % find indices for starting and ending time points
    if(tStart > ts.Time(end))
      warning(['Start time is greater than last point in time for timeseries: ' ts.Name]);
      idxStart = -1;
      idxEnd = -1;
      return;
    else
      idxStart = find(ts.Time >= tStart, 1);
    end
    if(tEnd < ts.Time(1))
      warning(['End time is smaller than first point in time for timeseries: ' ts.Name]);
      idxStart = -1;
      idxEnd = -1;
    else
      idxEnd = find(ts.Time <= tEnd, 1, 'last');
    end
    if(tStart >= tEnd)
      warning(['Start time is larger than end time']);
      idxStart = -1;
      idxEnd = -1;
    end
end

author: TUMFTM


text:
classdef FEM < handle

    properties (Access = protected)
        inputReader
    end

    methods (Access = public)

        function obj = FEM()
            obj.inputReader = FemInputReader_GiD();
        end

    end

    methods (Static, Access = public)

        function obj = create(s)
            switch s.type
                case 'ELASTIC'
                    switch s.scale
                        case 'MACRO'
                            obj = ElasticProblem(s);
                        case 'MICRO'
                            obj = ElasticProblemMicro(s);
                    end
                case 'THERMAL'
                    obj = ThermalProblem(s);
                case 'DIFF-REACT'
                    obj = DiffReactProblem(s);
                case 'HYPERELASTIC'
                    obj = Hyperelastic_Problem(s);
                case 'Stokes'
                    obj = StokesProblem(s);
            end
        end

    end

end


text:
function  gnss = gnss_m2r (lat, h, gnss)
% gnss_m2r: converts GPS standard deviation from meters to radians.
%
% INPUT
%   gnss, GNSS data structure with fields:
%       lat, 1x1 latitude (radians).
%       h,   1x1 altitude (meters).
%       stdm, 1x3 position error profile (m, m, m).
%
% OUTPUT
%   gnss.std, 1x3 position error profile (rad, rad, m).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%   R. Gonzalez, J. Giribet, and H. Patiño. NaveGo: a
% simulation framework for low-cost integrated navigation systems,
% Journal of Control Engineering and Applied Informatics, vol. 17,
% issue 2, pp. 110-120, 2015. Eq. 20.
%
%  	R. Gonzalez, J. Giribet, and H. Patiño. An approach to
% benchmarking of loosely coupled low-cost navigation systems,
% Mathematical and Computer Modelling of Dynamical Systems, vol. 21,
% issue 3, pp. 272-287, 2015. Eq. 7.
%
% Version: 003
% Date:    2021/03/09
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

gnss.std =  zeros(1,3);

[RM, RN] = radius(lat);

gnss.std(1) = gnss.stdm(1) / (RM + h);
gnss.std(2) = gnss.stdm(2) / (RN + h) / cos (lat);
gnss.std(3) = gnss.stdm(3);

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function opt = DatesCount(fo_year,lo_year,frequency,fo_period,lo_period)
% =======================================================================
% For given initial observation, final observation and frequency this
% function computes the number of observations between them and creates a
% cell aray with the dates (of the type 1979M2).
% =======================================================================
% opt = DatesCount(fo_year,lo_year,frequency,fo_period,lo_period)
% -----------------------------------------------------------------------
% INPUT
%	- fo_year: initial year of the timeline
%	- lo_year: final year of the timeline
%	- frequency: quarterly 'q' [def], monthly 'm', yearly 'y'
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%   - fo_period: initial quarter/month of the timeline (not for yearly frequency)
%   - lo_period: last quarter/month of the timeline (not for yearly frequency)
% -----------------------------------------------------------------------
% OUTPUT
%	- nobs: number of observations
%	- dates: cell array of dates
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

% Check inputs
if ~exist('frequency','var')
    frequency = 'q';
end

if strcmp(frequency,'m')
    if ~exist('fo_period','var'), error('You need to provide first period'), end
    if ~exist('lo_period','var'), error('You need to provide last period'), end
    fo = fo_year + fo_period*(1/12);
    lo = lo_year + lo_period*(1/12);
    aux = fo:(1/12):lo;
    nobs = length(aux);
elseif strcmp(frequency,'q')
    if ~exist('fo_period','var'), error('You need to provide first period'), end
    if ~exist('lo_period','var'), error('You need to provide last period'), end
    fo = fo_year + fo_period*(1/4);
    lo = lo_year + lo_period*(1/4);
    aux = fo:(1/4):lo;
    nobs = length(aux);
elseif strcmp(frequency,'y')
    fo = fo_year;
    lo = lo_year;
    aux = fo:lo;
    nobs = length(aux);
end

opt.dates = DatesCreate(fo_year,fo_period,nobs,frequency);
opt.nobs

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function visual_odo_main(visualDataFile,groundTruthFile)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     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/>.

    %% vo
    % for drawing images for estimating pitch and yaw
    global SUB_TRANS_IMG;
    global SUB_YAW_ROT_IMG;
    global SUB_HEIGHT_V_IMG;

    global ODO_IMG_TRANS_RESIZE_RANGE;
    global ODO_IMG_YAW_ROT_RESIZE_RANGE;
    global ODO_IMG_HEIGHT_V_RESIZE_RANGE;

    global PREV_TRANS_V_IMG_X_SUMS;
    global PREV_YAW_ROT_V_IMG_X_SUMS;
    global PREV_HEIGHT_V_IMG_Y_SUMS;

    global KEY_POINT_SET;
    global ODO_STEP;
    %% read imgs as VT input

    % get the visual data info
    [subFoldersPathSet, numSubFolders] = get_images_data_info(visualDataFile);

    % load ground truth data
    if nargin > 1
        [frameId, gt_x, gt_y, gt_z, gt_rx, gt_ry, gt_rz] = load_ground_truth_data(groundTruthFile);
        groundTruthData = [gt_x gt_y gt_z];
    end

    % define temp variables of previous profiles
    preProfilesTransImg = zeros(1,ODO_IMG_TRANS_RESIZE_RANGE(2));
    preProfilesYawRotImg = zeros(1,ODO_IMG_YAW_ROT_RESIZE_RANGE(2));
    preProfilesPitchRotImg = zeros(ODO_IMG_HEIGHT_V_RESIZE_RANGE(2), 1);

    % For drawing the current rot direction
    curRotDir = [0 0 0];
    curRotDir(1,3) = 0;

    % The total rotational angle
    sumRotAngle = 0;

    % for drawing the direction theta (rotational angle)
    theta = [0 0 0 0 0 0];
    rho = [0 0 0 0 0 0];
    startPoint =[0 0];
    endPoint = [0.8 0.8];

    % For drawing the velocity
    subRotVel = [];
    transVelVector = [];
    heightVelVector = [];
    sumHeight = [];
    transVelVector(1) = 0;
    heightVelVector(1) = 0;
    subRotVel(1) = 0;
    sumHeight = 0;

     % Transform the angle from degree to radian
    global DEGREE_TO_RADIAN;
    DEGREE_TO_RADIAN = pi / 180;

    % Transform the angle from radian to degree
    global RADIAN_TO_DEGREE;
    RADIAN_TO_DEGREE = 180 / pi;

    RENDER_RATE = 2;

    IMG_TYPE = '*.png';

    global OFFSET_YAW_ROT;
    offset_yaw_rot_vector = [];
    offset_yaw_rot_vector(1) = 0;

    global OFFSET_HEIGHT_V;
    offset_height_V_vector = [];
    offset_height_V_vector(1) = 0;

    odoMapTrajectory(1,1) = 0;
    odoMapTrajectory(1,2) = 0;
    odoMapTrajectory(1,3) = 0;
    odoMapTrajectory(1,4) = 0;

    ODO_MAP_X_SCALING = -1;
    ODO_MAP_Y_SCALING = 1;
    ODO_MAP_Z_SCALING = 1;
    GT_ODO_X_SCALING = 1;
    GT_ODO_Y_SCALING = 1;
    GT_ODO_Z_SCALING = 1;

    curFrame = 0;
    startFrame = 1;

    %% processing visual odometry
    for iSubFolder = 1 : numSubFolders

        [curFolderPath, imgFilesPathList, numImgs] = get_cur_img_files_path_list(subFoldersPathSet, IMG_TYPE, iSubFolder);

        if numImgs > 0
            for indFrame = startFrame :ODO_STEP: numImgs-1
                curFrame = curFrame + 1;
                [curImg] = read_current_image(curFolderPath, imgFilesPathList, indFrame);

                % visual templates and visual odo uses intensity so convert to grayscale
                curGrayImg = rgb2gray(curImg);
                curGrayImg = im2double(curGrayImg);

                % for drawing previous intensity
                preProfilesTransImg = PREV_TRANS_V_IMG_X_SUMS;
                preProfilesYawRotImg = PREV_YAW_ROT_V_IMG_X_SUMS;
                preProfilesHeightVImg = PREV_HEIGHT_V_IMG_Y_SUMS;

                %% step2: processing vo
                % get the odometry from the video

               % for common use
               %  [transV, yawRotV, heightV] = visual_odometry(curGrayImg);

                % for special use
                if length(KEY_POINT_SET) == 2
                    if indFrame < KEY_POINT_SET(1)
                        [transV, yawRotV, heightV] = visual_odometry(curGrayImg);
                    elseif indFrame < KEY_POINT_SET(2)
                        [transV, yawRotV, heightV] = visual_odometry_up(curGrayImg);
                    else
                        [transV, yawRotV, heightV] = visual_odometry(curGrayImg);
                    end
                    transV = 2;
                else
                    if indFrame < KEY_POINT_SET(1)
                        [transV, yawRotV, heightV] = visual_odometry(curGrayImg);
                    elseif indFrame < KEY_POINT_SET(2)
                        [transV, yawRotV, heightV] = visual_odometry_up(curGrayImg);
                    elseif indFrame < KEY_POINT_SET(3)
                        [transV, yawRotV, heightV] = visual_odometry(curGrayImg);
                    elseif indFrame < KEY_POINT_SET(4)
                        [transV, yawRotV, heightV] = visual_odometry_down(curGrayImg);
                    else
                        [transV, yawRotV, heightV] = visual_odometry(curGrayImg);
                    end
                end

                subRotVel(curFrame) = yawRotV;
                transVelVector(curFrame) = transV;
                heightVelVector(curFrame) = heightV;

                offset_yaw_rot_vector(curFrame) = OFFSET_YAW_ROT;
                offset_height_V_vector(curFrame) = OFFSET_HEIGHT_V;

                sumHeight(curFrame + 1) = sumHeight(curFrame) + heightV;

                % For drawing the total rotational angle
                if  sumRotAngle + yawRotV >= 360
                    sumRotAngle = mod(sumRotAngle + yawRotV, 360);
                elseif sumRotAngle + yawRotV <= -360
                    sumRotAngle = sumRotAngle + yawRotV + 360;
                else
                    sumRotAngle = sumRotAngle + yawRotV;
                end

                % Transform the vrot from degree to radian
                yawRotV = yawRotV * DEGREE_TO_RADIAN;

                % Direction by raw odometry
                curRotDir(curFrame + 1, 1) = cos(curRotDir(curFrame, 3) + yawRotV );
                curRotDir(curFrame + 1, 2) = sin(curRotDir(curFrame, 3) + yawRotV );
                curRotDir(curFrame + 1, 3) = curRotDir(curFrame, 3) + yawRotV;

                odoMapTrajectory(curFrame + 1,1) = odoMapTrajectory(curFrame,1) + transV * cos(sym(curRotDir(curFrame + 1, 3)));
                odoMapTrajectory(curFrame + 1,2) = odoMapTrajectory(curFrame,2) + transV * sin(sym(curRotDir(curFrame + 1, 3)));
                odoMapTrajectory(curFrame + 1,3) = odoMapTrajectory(curFrame,3) + heightV;
                odoMapTrajectory(curFrame + 1,4) = curRotDir(curFrame + 1, 3);

                % render results information
                if (mod(curFrame, RENDER_RATE) == 0)

                    % computing the normalized profile of horizontal translational image
                    profilesTransImg = sum(SUB_TRANS_IMG);
                    avgIntensity = sum(profilesTransImg) / size(profilesTransImg, 2);
                    profilesTransImg = profilesTransImg / avgIntensity;

                    % computing the normalized profile of rotational image
                    profilesYawRotImg = sum(SUB_YAW_ROT_IMG);
                    avgIntensity = sum(profilesYawRotImg) / size(profilesYawRotImg, 2);
                    profilesYawRotImg = profilesYawRotImg / avgIntensity;

%                     diffYawRotImgs = profilesYawRotImg - preProfilesYawRotImg;
                    % computing the normalized profile of vertical transV image
                    profilesHeightVImg = sum(SUB_HEIGHT_V_IMG,2);
                    avgIntensity = sum(profilesHeightVImg) / size(profilesHeightVImg, 1);
                    profilesHeightVImg = profilesHeightVImg / avgIntensity;

                    diffHeightVImgs = profilesHeightVImg - preProfilesHeightVImg;

                    figure(1)

                    % render the sub image for estimating horizontal translational velocity
                    subplot(3, 4, 1, 'replace');
                    imshow(SUB_TRANS_IMG, []);
                    title('Sub image (vHoriTrans)');
                    xlabel('Width');
                    ylabel('Height');
                    axis on

                    % draw the the normalized profile of horizontal translational image
                    subplot(3,4,5, 'replace');
                    hold on
                    plot(1: size(profilesTransImg,2),profilesTransImg, 'r');
                    plot(1: size(preProfilesTransImg,2),preProfilesTransImg, 'g');
                    hold off
%                     legend('curr', 'g: prev');
                    title('The intensity profiles');
                    xlabel('Width');
                    ylabel('Normalized profile');
                    axis on

                    % render the horizontal translational velocity
                    subplot(3, 4, 9, 'replace');
                    plot(transVelVector);
                    title('Translational velocity');
                    xlabel('Time');
                    ylabel('vTrans');
                    axis on

                    % render the sub image for estimating rotational velocity
                    subplot(3, 4, 2, 'replace');
                    imshow(SUB_YAW_ROT_IMG, []);
                    title('Sub image (vYawRot)');
                    xlabel('Width');
                    ylabel('Height');
                    axis on

                    % draw the the normalized profile of rotational image
                    subplot(3,4,6, 'replace');
                    hold on
                    plot(1: size(profilesYawRotImg,2),profilesYawRotImg, 'r');
                    plot(1: size(preProfilesYawRotImg,2),preProfilesYawRotImg, 'g');
                    hold off
%                     legend('r: curr', 'g: prev');
                    title('The intensity profiles');
                    xlabel('Width');
                    ylabel('Normalized profile');
                    axis on

                    % render the horizontal translational velocity
                    subplot(3, 4, 10, 'replace');
                    plot(subRotVel);
                    title('Yaw rotational velocity');
                    xlabel('Time');
                    ylabel('vYawRot');
                    axis on

                    % render the sub image for estimating pitch velocity
                    subplot(3, 4, 3, 'replace');
                    imshow(SUB_HEIGHT_V_IMG, []);
                    title('Sub image (vPitchRot)');
                    xlabel('Width');
                    ylabel('Height');
                    axis on

                    % draw the the normalized profile of height v image
                    subplot(3,4,7, 'replace');
                    hold on
                    plot(profilesHeightVImg, 1: size(profilesHeightVImg,1), 'r');
                    plot(preProfilesHeightVImg, 1: size(preProfilesHeightVImg,1), 'g');
                    hold off
%                     legend('r: curr', 'g: prev');
                    title('The intensity profiles');
                    xlabel('Width');
                    ylabel('Normalized profile');
                    axis on

                    % render the horizontal translational velocity
                    subplot(3, 4, 11, 'replace');
                    plot(heightVelVector);
                    title('Vertical translational velocity');
                    xlabel('Time');
                    ylabel('vVertTrans');
                    axis on

                    % render current rot direction
                    subplot(3, 4, 4, 'replace');
                    polar(theta,rho);
                    hold on
                    plot([startPoint(1) curRotDir(curFrame,1)],[startPoint(2) curRotDir(curFrame,2)],'linewidth',1,'color',[0.9 0 0]);
                    title('Current rotational direction');

                    % render the horizontal translational velocity
                    subplot(3, 4, 8, 'replace');
                    plot(sumHeight);
                    title('Height Change');
                    xlabel('Time');
                    ylabel('Height');
                    axis on

                    % render the horizontal translational velocity
                    subplot(3, 4, 12, 'replace');
                    hold on
                    plot(offset_yaw_rot_vector, 'r');
                    plot(offset_height_V_vector, 'g');
                    hold off
                    title('Offset yaw rotation');
                    xlabel('Time');
                    ylabel('Offset');
                    axis on

                    figure(2)

                    % draw odo map
                    subplot(1, 2, 1, 'replace');
                    hold on
                    plot3((odoMapTrajectory(:,2))* ODO_MAP_X_SCALING , ...
                        (odoMapTrajectory(:,1)) * ODO_MAP_Y_SCALING, ...
                        (odoMapTrajectory(:,3) ) * ODO_MAP_Z_SCALING, '.b');

%                     plot3((gt_x(1:curFrame) - gt_x(1)) * GT_ODO_X_SCALING, ...
%                         (gt_y(1:curFrame) - gt_y(1)) * GT_ODO_Y_SCALING, ...
%                         (gt_z(1:curFrame) - gt_z(1)) * GT_ODO_Z_SCALING, '.r');
                    hold off;
                    view(3);
                    grid on;
%                     view(0, 90);
                    xl = xlabel('odo-map-x');
                    yl = ylabel('odo-map-y');
                    zlabel('odo-map-z');
                    set(xl,'Rotation',15);
                    set(yl,'Rotation',-30);
                    title('Multilayered Odometry Map');
%                     legend('Result','Truth' ,'1');
%                     axis([-5 110 -5 100 -2 5]);
                    axis on
                    rotate3d on

                    % draw odo map
                    subplot(1, 2, 2, 'replace');
                    hold on
                    plot3((odoMapTrajectory(:,2))* ODO_MAP_X_SCALING , ...
                        (odoMapTrajectory(:,1)) * ODO_MAP_Y_SCALING, ...
                        (odoMapTrajectory(:,3) ) * ODO_MAP_Z_SCALING, '.b');

%                     plot3((gt_x(1:curFrame) - gt_x(1)) * GT_ODO_X_SCALING, ...
%                         (gt_y(1:curFrame) - gt_y(1)) * GT_ODO_Y_SCALING, ...
%                         (gt_z(1:curFrame) - gt_z(1)) * GT_ODO_Z_SCALING, '.r');
                    hold off;
                    grid on;
                    view(0, 90);
                    xlabel('odo-map-x');
                    ylabel('odo-map-y');
                    zlabel('odo-map-z');
                    title('Multilayered Odometry Map (top view)');
%                     legend('Result','Truth' ,'1');
%                      axis([-5 110 -5 100]);
                    axis on

%                     figure(3)
%
%                      subplot(1,1,1, 'replace');
%                     hold on
%                     plot(1: size(diffHeightVImgs,1),diffHeightVImgs, 'r');
%                     hold off
% %                     legend('r: curr', 'g: prev');
%                     title('The diff between two intensity profiles');
%                     xlabel('Width');
%                     ylabel('Diff');
%                     axis on

                    drawnow
                end
            end
        end
    end
end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text: %GETDATA Get data of datafile

author: marianux

repo: ecg-kit

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