text:
function [S] = ukfPropagation(dt,Rot,RotAnt,vAnt,omega_b,a_b,S,omega,acc,Qc,g)
q = length(S);
S_aug = blkdiag(S,Qc);
N_aug = length(S_aug);

% scaled unsented transform
W0 = 1-N_aug/3;
Wj = (1-W0)/(2*N_aug);
gamma = sqrt(N_aug/(1-W0));

ichi = Rot';
omega = omega-omega_b; % unbiased inputs
acc = acc-a_b;
X = gamma*[zeros(N_aug,1) S_aug' -S_aug'];% sigma-points
X(10:15,:) = X(10:15,:)+X(q+7:q+12,:)*dt;
for j = 2:2*N_aug+1
    xi_j = X(1:3,j);
    xi_v = X(4:6,j);
    xi_x = X(7:9,j);
    w_j = X(q+1:N_aug,j);
    omega_bj = X(10:12,j);
    a_bj = X(13:15,j);
    Rot_j = RotAnt*expSO3(xi_j)*expSO3((omega+w_j(1:3)-omega_bj)*dt);
    v_j = vAnt + xi_v + (Rot_j*(acc+w_j(4:6)-a_bj)+g)*dt;
    x_j = xi_x + v_j*dt;
    Xi_j = ichi*Rot_j;
    X(1:3,j) = logSO3(Xi_j); % propagated sigma points
    X(4:6,j) = v_j - vAnt;
    X(7:9,j) = x_j + (v_j-vAnt)*dt;
end
X = sqrt(Wj)*X;
[~,Rs] = qr(X(1:q,2:2*N_aug+1)');
S = Rs(1:q,1:q);
end

author: mbrossar


text:
function  [ gx,dgdx,dgdP ] = g_softmax(x,P,u,in )
% softmax decision rule for Q-learning (2-armed bandit task)
% function  [ gx,dgdx,dgdP ] = g_softmax(x,P,u,in )
% IN:
%   - x : Q-values
%   - P : inverse (log-) temperature and bias
%   - u : [useless]
%   - in : [useless]
% OUT:
%   - gx : P(a=1|x)

% inverse temperature
% -------------------------------------------------------------------------
beta = exp(P(1)); % exp: [-Inf,Inf] -> [0 Inf]

dQ = (x(1)-x(2));
if length(P)>1
    gx = VBA_sigmoid( beta*dQ + P(2));
else
    gx = VBA_sigmoid( beta*dQ );
end
dgdx = zeros(size(x,1),1);
dgdx(1) = beta*gx*(1-gx);
dgdx(2) = -beta*gx*(1-gx);
if length(P)>1
    dgdP = [beta*dQ*gx*(1-gx),gx*(1-gx)];
else
    dgdP = [beta*dQ*gx*(1-gx)];
end

author: MBB-team

repo: VBA-toolbox


text:
function enlarged_tt = add_non_essential_dims(small_tt, new_n, old_vars_pos)
    % enlarged_tt = add_non_essential_dims(small_tt, new_n, old_vars_pos)
    % creates enlarged version of tt_tensor, expanding it by ones.
    % Mode sizes of the enlarged_tt are equal to new_n.
    % First dimension of the small_tt is old_vars_pos(1) dimension of the enlarged_tt.
    % This procedure doesn't increase the maximal TT-rank.
    %
    % Input:
    % 	small_tt     -- P-dimensional tensor in the TT-format.
    % 	new_n        -- Q by 1 vector with enlarged_tt mode sizes (Q > P).
    % 	old_vars_pos -- P by 1 vector in ascending order. Consists of positions
    % 							of the small_tt dimensions in the enlarged_tt tensor.
    %
    % Output:
    % 	enlarged_tt  -- Q-dimensional tensor in the TT format.
    % 			enlarged_tt(i1, ..., iQ)
    % 			i1 take values from {1, ..., new_n(1)}
    % 			...
    % 			iQ take values from {1, ..., new_n(Q)}
    %
    % Example:
    %       m = rand(3, 4);
    %       small_tt = tt_tensor(m); % small_tt(i2, i4)
    %       enlarged_tt = add_non_essential_dims(small_tt, [2, 3, 6, 4, 4], [2, 4]);
    %       % enlarged_tt(i1, i2, i3, i4, i5, i6) == small_tt(i2, i4) for any i1, i3, i5, i6
    %

    if (~issorted(old_vars_pos))
        error('Dimensions (third argument) must be in ascending order.');
    end

    if ~isvector(new_n) | ~isvector(old_vars_pos) | small_tt.d ~= length(old_vars_pos)
    	error('Wrong usage, see help tt_tensor.add_non_essential_dims.');
    end

    d = length(new_n);
    enlarged_tt = tt_tensor;
    enlarged_tt.d = d;
    enlarged_tt.r = zeros(d + 1, 1);
    enlarged_tt.r(1) = 1;
    enlarged_tt.n = new_n(:);
    enlarged_tt.core = [];
    enlarged_tt.ps = zeros(d + 1, 1);
    enlarged_tt.over = 0;
    curr_var_i = 0;
    vars_length = length(old_vars_pos);
    before_first_var = true;
    after_last_var = false;
    for dimension_i = 1:d
        if curr_var_i == vars_length
            after_last_var = true;
        end

        if curr_var_i < vars_length && old_vars_pos(curr_var_i + 1) == dimension_i
            curr_var_i = curr_var_i + 1;
            before_first_var = false;
            if new_n(dimension_i) ~= small_tt.n(curr_var_i)
                error('Dimensions must agree!');
            end
        end

        if before_first_var || after_last_var
            % Dimension_i is before old_vars_pos(1) or after old_vars_pos(end).
            curr_core = ones(new_n(dimension_i), 1);
            curr_rank = 1;
        elseif old_vars_pos(curr_var_i) == dimension_i
            % It's real dimension.
            curr_core = core(small_tt, curr_var_i);
            curr_rank = small_tt.r(curr_var_i);
        else
            % Non-essential dimension between two real ones.
            curr_rank = small_tt.r(curr_var_i + 1);
            curr_core = core_eye(curr_rank, new_n(dimension_i));
        end
        enlarged_tt.r(dimension_i) = curr_rank;
        enlarged_tt.ps(dimension_i) = length(enlarged_tt.core) + 1;
        enlarged_tt.core = [enlarged_tt.core; curr_core(:)];
    end
    enlarged_tt.r(end) = 1;
    enlarged_tt.ps(end) = length(enlarged_tt.core) + 1;
end

function core = core_eye(r, n)
    % Return 3-dimensional array filled with eye matrices,
    % which corresponds to non-essential dimension core.
    core = zeros(r, n, r);
    for i = 1:n
        core(:, i, :) = eye(r);
    end
end

author: oseledets

repo: TT-Toolbox


text:
function K=fitTwoLines(x,y, showimage)
% K=fitTwoLines(x,y, showimage)
% This fits two adjecent lines into data [x,y] which are partitioed as
% x=[x1,x2] and y=[y1,y2]. It determines K which gives best fit for
% partition x1=x(1:K), x2=x(K:end) ...
% Used for estimation of the 'kink' in the graph.
if ~exist('showimage','var')
    showimage = 0;
end
n=length(y);

for ind_includeMidPoint=0:1 % when 0 then middle point is included (shared by both lines)
    index =0;
    for ii=1:n-ind_includeMidPoint
        index = index + 1;
        x1=x(1:ii);
        y1=y(1:ii);
        x2=x(ii+ind_includeMidPoint:n);
        y2=y(ii+ind_includeMidPoint:n);

        [pbest,perror,nchi]=nonlinft('linearfunction' ,x1,y1,ones(size(x1)),[-1 1e-5],[1 1]);
        pbest1(index,:,ind_includeMidPoint+1)=pbest;
        perror1(index,:,ind_includeMidPoint+1)=perror;
        nchi1(index,ind_includeMidPoint+1)=nchi*(numel(y1)-1); % to compensate for degrees of freedom
        [pbest,perror,nchi]=nonlinft('linearfunction' ,x2,y2,ones(size(x2)),[-1 1e-5],[1 1]);
        pbest2(index,:,ind_includeMidPoint+1)=pbest;
        perror2(index,:,ind_includeMidPoint+1)=perror;
        nchi2(index,ind_includeMidPoint+1)=nchi*(numel(y2)-1); % to compensate for degrees of freedom
    end
end
nn=squeeze((sqrt(perror2(:,1,:).^2+perror2(:,2,:).^2)));
objectiveFunction = (nchi1.^2+nchi2.^2);
mo=min(objectiveFunction,[],2); % minimum the two either including (1) or not (2) the middle point
[del, ixmo]=min(objectiveFunction./(1./nn),[],2); %whether to choose (1) or (2). This is weighted by error of the estimation...
% [m,mx]=min(objectiveFunction);
% [m,mx]=min(mo)+ixmo-1;
% K=x(mx+1);
[m,mx]=min(mo(1:end-1));
K=x(mx+ixmo(mx)-1);

if showimage
    figure;
    plot(x(1:end), mo)
    vline2(K,'r','K');
    xlabel ('x')
    ylabel ('objective function')

    figure;
    plot(x,y,'-o')
    hold on
    plot(x(1:mx), pbest1(mx,1,ixmo(mx))*x(1:mx)+pbest1(mx,2,ixmo(mx)),'g')
    plot(x(mx+ixmo(mx)-1:n), pbest2(mx,1, ixmo(mx))*x(mx+ixmo(mx)-1:n)+pbest2(mx,2,ixmo(mx)),'r')
    xlabel('x')
    ylabel('y')
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function ap_multi = compute_AP_multiCam(good_image, junk_image, index, queryCam, testCam)
good_cam = testCam(good_image);
good_cam_uni = unique(good_cam);
ap_multi = zeros(1, 6);

% on the same camera
good_cam_now = queryCam;
ngood = length(junk_image);
junk_image_now = good_image;
good_image_now = junk_image;
old_recall = 0;
old_precision = 1.0;
ap = 0;
intersect_size = 0;
j = 0;
good_now = 0;
for n = 1:length(index)
    flag = 0;
    if ~isempty(find(good_image_now == index(n), 1))
        flag = 1; % good image
        good_now = good_now+1;
    end
    if ~isempty(find(junk_image_now == index(n), 1))
        continue; % junk image
    end

    if flag == 1%good
        intersect_size = intersect_size + 1;
    end
    if ngood == 0
        ap_multi(good_cam_now) = 0;
        break;
    end
    recall = intersect_size/ngood;
    precision = intersect_size/(j + 1);
    ap = ap + (recall - old_recall)*((old_precision+precision)/2);
    old_recall = recall;
    old_precision = precision;
    j = j+1;

    if good_now == ngood
        ap_multi(good_cam_now) = ap;
        break;
    end
end

for k = 1:length(good_cam_uni)
    good_cam_now = good_cam_uni(k);
    ngood = length(find(good_cam == good_cam_now));
    pos_junk = find(good_cam ~= good_cam_now);
    junk_image_now = [junk_image good_image(pos_junk)];
    pos_good = find(good_cam == good_cam_now);
    good_image_now = good_image(pos_good);
    old_recall = 0;
    old_precision = 1.0;
    ap = 0;
    intersect_size = 0;
    j = 0;
    good_now = 0;
    for n = 1:length(index)
        flag = 0;
        if ~isempty(find(good_image_now == index(n), 1))
            flag = 1; % good image
            good_now = good_now+1;
        end
        if ~isempty(find(junk_image_now == index(n), 1))
            continue; % junk image
        end

        if flag == 1%good
            intersect_size = intersect_size + 1;
        end
        recall = intersect_size/ngood;
        precision = intersect_size/(j + 1);
        ap = ap + (recall - old_recall)*((old_precision+precision)/2);
        old_recall = recall;
        old_precision = precision;
        j = j+1;

        if good_now == ngood
            ap_multi(good_cam_now) = ap;
            break;
        end
    end
end

end

author: naiq


text:
function a = rutis1_inverse ( )

%*****************************************************************************80
%
%% RUTIS1_INVERSE returns the inverse of the RUTIS1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joan Westlake,
%    A Handbook of Numerical Matrix Inversion and Solution of
%    Linear Equations,
%    John Wiley, 1968,
%    ISBN13: 978-0471936756,
%    LC: QA263.W47.
%
%  Parameters:
%
%    Output, real A(4,4), the matrix.
%

%
%  Note that the matrix entries are listed by row.
%
  a(1:4,1:4) = [ ...
    -2.0,  4.0,  4.0, -5.0; ...
     4.0, -2.0, -5.0,  4.0; ...
     4.0, -5.0, -2.0,  4.0; ...
    -5.0,  4.0,  4.0, -2.0 ];

  a(1:4,1:4) = a(1:4,1:4) / 15.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function r8po_print_some ( n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8PO_PRINT_SOME prints some of a R8PO matrix.
%
%  Discussion:
%
%    The R8PO storage format is appropriate for a symmetric positive definite
%    matrix and its inverse.  (The Cholesky factor of a R8PO matrix is an
%    upper triangular matrix, so it will be in R8GE storage format.)
%
%    Only the diagonal and upper triangle of the square array are used.
%    This same storage scheme is used when the matrix is factored by
%    R8PO_FA, or inverted by R8PO_INVERSE.  For clarity, the lower triangle
%    is set to zero.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N,N), the R8PO matrix.
%
%    Input, integer ILO, JLO, IHI, JHI, the first row and
%    column, and the last row and column to be printed.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  incx = 5;
%
%  Print the columns of the matrix, in strips of 5.
%
  for j2lo = jlo : incx : jhi

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    inc = j2hi + 1 - j2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d       ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '  ---\n' );
%
%  Determine the range of the rows in this strip.
%
    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, n );

    for i = i2lo : i2hi

      fprintf ( 1, '%4d  ', i );
%
%  Print out (up to) 5 entries in row I, that lie in the current strip.
%
      for j2 = 1 : inc

        j = j2lo - 1 + j2;

        if ( i <= j )
          aij = a(i,j);
        else
          aij = a(j,i);
        end

        fprintf ( 1, '%12g  ', aij );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function net = resnet50_intialization_down4(warpsize,cdim,ndim,RoIRoD)

netStruct = load('data/pretrained_models/imagenet-resnet-50-dag.mat') ;
net = dagnn.DagNN.loadobj(netStruct) ;
net.layers = net.layers(1:140);
net.renameVar('data','input');
net.rebuild();

% Dimension reduction
net.addLayer('dimred', convBlocknoPad(1,1,1024,cdim), {'res4fx'}, {'xRed'}, {'dimred_filter','dimred_bias'});
net = initLayerParams(net);

% RoIAlign
if strcmp(RoIRoD,'RoIRoD') || strcmp(RoIRoD,'RoIOnly')
    roiSize = [warpsize,warpsize];
    net.addLayer('roialign', dagnn.RoiAlign('subdivisions',roiSize,'transformation',1/16), {'xRed','rois'}, 'xRoi');
    net.addLayer('roipool', dagnn.Pooling('method', 'avg', 'poolSize', [2 2], 'stride', 1), {'xRoi'}, {'pool_xRoi'}, {});
end

% RoDAlign
if strcmp(RoIRoD,'RoIRoD') || strcmp(RoIRoD,'RoDOnly')
    rodSize = [warpsize+1,warpsize+1];
    net.addLayer('rodalign', dagnn.RodAlign('subdivisions',rodSize,'transformation',1/16), {'xRed','rois'}, 'xRod');
    net.addLayer('rodpool', dagnn.Pooling('method', 'avg', 'poolSize', [2 2], 'stride', 1), {'xRod'}, {'pool_xRod'}, {});
end

switch RoIRoD
    case 'RoIRoD'
        net.addLayer('concat', dagnn.Concat('dim', 3), {'pool_xRoi', 'pool_xRod'}, {'xConc'}, {});
        net.addLayer('fc1', convBlocknoPad(warpsize,warpsize,cdim*2,ndim), {'xConc'}, {'fc1'}, {'fc1_filter', 'fc1_bias'});
    case 'RoIOnly'
        net.addLayer('fc1', convBlocknoPad(warpsize,warpsize,cdim,ndim), {'pool_xRoi'}, {'fc1'}, {'fc1_filter', 'fc1_bias'});
    case 'RoDOnly'
        net.addLayer('fc1', convBlocknoPad(warpsize,warpsize,cdim,ndim), {'pool_xRod'}, {'fc1'}, {'fc1_filter', 'fc1_bias'});
    otherwise
        error('Unknown RoIRoD choice.');
end
net = initLayerParams(net);
net.addLayer('relu_fc1', dagnn.ReLU(), {'fc1'}, {'relu_fc1'}, {});
net.addLayer('fc2', convBlocknoPad(1,1,ndim,ndim), {'relu_fc1'}, {'fc2'}, {'fc2_filter', 'fc2_bias'});
net = initLayerParams(net);
net.addLayer('relu_fc2', dagnn.ReLU(), {'fc2'}, {'relu_fc2'}, {});

net.addLayer('dropout',dagnn.DropOut('rate',0.5),'relu_fc2','dropout',{});

net.addLayer('predcls',convBlocknoPad(1,1,ndim,1), 'dropout','predcls',{'predcls_filter','predcls_bias'});
net = initLayerParams(net);

% net.print({'input', [256 256 3]}, 'all', true);
net.addLayer('losscls',dagnn.RegressionLoss('lossType','Huber'), {'predcls','label'}, 'losscls',{});

% Meta parameters
net.meta.trainOpts.learningRate = logspace(-4,-4,20);%logspace(-4,-5,10)
net.meta.trainOpts.weightDecay = 0.0001;
net.meta.trainOpts.batchSize = 4 ;
net.meta.trainOpts.numSubBatches = 4 ;
net.meta.trainOpts.solver = @solver.adam ;
net.meta.trainOpts.numEpochs = numel(net.meta.trainOpts.learningRate) ;
net.meta.normalization.averageImage = reshape([122.7717 102.9801 115.9465],[1 1 3]);

end

function convObj = convBlocknoPad(fh,fw,fc,k)
    convObj = dagnn.Conv('size', [fh fw fc k], 'hasBias', true);
end

function net = initLayerParams(net)
    p = net.getParamIndex(net.layers(end).params) ;
    params = net.layers(end).block.initParams();
    [net.params(p).value] = deal(params{:}) ;
end

author: HuiZeng

repo: Grid-Anchor-based-Image-Cropping


text:
function angles = polyangles(x,y)
%POLYANGLES Computes internal polygon angles.
%	ANGLES = POLYANGLES(X,Y) computes the interior angles (in degrees) of
%	an arbitrary polygon whose vertices have x- and y-coordinates given
%	in column vectors X and Y. The vertices must be arranged in a
%	clockwise manner.  The program eliminates duplicate adjacent rows in
%	[X,Y], except that the first row may equal the last, so that the
%	polygon is closed.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Preliminaries.
xy = [x(:),y(:)];
if isempty(xy)
    % No vertices!
    angles = zeros(0,1);
    return
end

% Close the polygon if necessary.
if (size(xy,1) == 1) || ~isequal(xy(1,:),xy(end,:))
    xy = [xy;xy(1,:)];
end

% Eliminate duplicate vertices.
if size(xy,1) > 2
    xy(all(diff(xy,1,1) == 0,2),:) = [];
end

% Form a set of vectors by taking the differences between adjacent
% vertices in the polygon.  These vectors form a chain directed in the
% clockwise direction because of the way the data is input.
v2 = diff(xy,1,1);

% Each angle of the polygon is formed by the head of a vector joining
% the tail of the vector following it in the sequence. Each vector in
% the following array precedes its corresponding vector in v2:
v1 = circshift(v2,[1 0]);

% Get the x- and y-components of each vector.
v1x = v1(:,1);
v1y = v1(:,2);
v2x = v2(:,1);
v2y = v2(:,2);

% The angles between the vectors in v2 relative to the vectors in v1 is
% given by the following expression. See, for example,
% http://www.euclideanspace.com/maths/algebra/vectors/ and (select
% "angle between" for an explanation of this expression).
angles = (180/pi)*(atan2(v2y,v2x) - atan2(v1y,v1x));

% Because the head of each vector in v1 points to the tail of the
% subsequent vector in v2 (see above), the preceding expression gives
% "outer" angles. We need the interior angles of the polygon. To do
% this, we have to reverse the direction of one of the vectors or,
% equivalently, add (or subtract) 180 degrees from the previous result.
% We also want to make sure that all angles are in the range 0 to 360
% degrees. The following expression implements both requirements:
angles = mod(angles + 180,360);

author: dipum

repo: dipum-toolbox

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


text:
function [v,x,t,m,ze]=v_quadpeak(z)
%V_PEAK2DQUAD find quadratically-interpolated peak in a N-D array
%
%  Inputs:  Z(m,n,...)   is the input array (ignoring trailing singleton dimensions)
%                        Note: a row vector will have 2 dimensions
%
% Outputs:  V        is the peak value
%           X(:,1)  is the position of the peak (in fractional subscript values)
%           T        is -1, 0, +1 for maximum, saddle point or minimum
%           M        defines the fitted quadratic: z = [x y ... 1]*M*[x y ... 1]'
%           ZE       the estimated version of Z

%	   Copyright (C) Mike Brookes 2008
%      Version: $Id: v_quadpeak.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

persistent wz a
% first calculate the fixed matrix, a (can be stored if sz is constant)
sz=size(z);         % size of input array
psz=prod(sz);       % number of elements in input array
dz=numel(sz);       % number of input dimensions
mz=find(sz>1);      % non-singleton dimension indices
nm=numel(mz);       % number of non-singleton dimensions
vz=sz(mz);          % size of squeezed input array
dx=max(mz);         % number of output dimensions
if ~nm              % if the input array is a scalar
    error('Cannot find peak of a scalar');
end
nc=(nm+1)*(nm+2)/2;  % number of columns in A matrix
if min(vz)<3
    error('Need at least 3 points in each non-singleton dimension');
end
if isempty(wz) || numel(wz)~=numel(vz) || ~all(wz==vz)
    wz=vz;
    a=ones(psz,nc);
    ix=(0:psz-1)';
    for i=1:nm
        jx=floor(ix/sz(mz(i)));
        a(:,i+nc-nm-1)=1+ix-jx*sz(mz(i));
        ix=jx;
        a(:,(i^2-i+2)/2:i*(i+1)/2)=a(:,nc-nm:i+nc-nm-1).*repmat(a(:,i+nc-nm-1),1,i);
    end
    a=(a'*a)\a';        % converts to polynomial coeficients {x^2 xy y^2 x y 1]
end

% now find the peak

c=a*z(:);   % polynomial coefficients for this data
w=zeros(nm+1,nm+1);
i=1:(nm+1)*(nm+2)/2;
j=floor((sqrt(8*i-7)-1)/2);
w(i+j.*(2*nm+1-j)/2)=c;
w=(w+w.')/2; % make it symmetrical
mr=w(1:nm,1:nm);
we=w(1:nm,nm+1);
y=-(mr\we);
v=y'*we+w(nm+1,nm+1);  % value at peak

% insert singleton dimensions into outputs

x=zeros(dx,1);
x(mz)=y;
m=zeros(dx+1,dx+1);
mz(nm+1)=dx+1;
m(mz,mz)=w;
if nargout>2
    ev=eig(mr);
    t=all(ev>0)-all(ev<0);
end
if nargout>4
    ze=zeros(sz);
    scp=cumprod([1 sz(1:end-1)]);
    ivec=fix(repmat((0:psz-1)',1,dz)./repmat(scp,psz,1));
    xe=[1+ivec-repmat(sz,psz,1).*fix(ivec./repmat(sz,psz,1)) ones(psz,1)];
    ze=reshape(sum((xe*m).*xe,2),sz);
end

if ~nargout && nm<=2
    % plot the data
    desc={'Maximum','Saddle Point','Minimum'};
    if nargout<=2
        ev=eig(mr);
        t=all(ev>0)-all(ev<0);
    end
    if nm==1
        xax=linspace(1,psz,100);
        plot(xax,c(1)*xax.^2+c(2)*xax+c(3),'-r',1:psz,z(:),'ob',x,v,'^k');
        set(gca,'xlim',[0.9 psz+0.1]);
        ylabel('z');
        xlabel(sprintf('x%d',mz(1)));
        title(sprintf('\\Delta = %s: z(%.2g) = %.2g',desc{t+2},y(1),v));
    else
        ngr=17;
        xax=repmat(linspace(1,vz(1),ngr)',1,ngr);
        yax=repmat(linspace(1,vz(2),ngr),ngr,1);
        zq=(c(1)*xax+c(2)*yax+c(4)).*xax+(c(3)*yax+c(5)).*yax+c(6);
        hold off
        mesh(xax,yax,zq,'EdgeColor','r');
        hold on
        plot3(repmat((1:vz(1))',1,vz(2)),repmat(1:vz(2),vz(1),1),reshape(z,vz),'ob',y(1),y(2),v,'^k');
        hold off
        set(gca,'xlim',[0.9 vz(1)+0.1],'ylim',[0.9 vz(2)+0.1]);
        xlabel(sprintf('x%d',mz(1)));
        ylabel(sprintf('x%d',mz(2)));
        zlabel('z');
        title(sprintf('\\Delta = %s: z(%.2g,%.2g) = %.2g',desc{t+2},y(1),y(2),v));
    end
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function planC = createStapleStruct(structNumV,probCutoff,structureName,planC)
% function planC = createStapleStruct(structNumV,probCutoff,structureName,planC)
%
% Function to create STAPLE agreement structure at the passed probability
% probCutoff.
%
% APA, 11/13/2017

if ~exist('planC','var')
    global planC
end
indexS = planC{end};

scanNum = getStructureAssociatedScan(structNumV(1),planC);
numObs = length(structNumV);
siz = getUniformScanSize(planC{indexS.scan}(scanNum));
rateMatM = false(prod(siz),numObs);
for i=1:numObs
    mask3M = getUniformStr(structNumV(i),planC);
    rateMatM(:,i) = mask3M(:);
end
indV = sum(rateMatM,2) > 0;
rateMatM = rateMatM(indV,:);

iterlim=100;
senstart=0.9999*ones(1,numObs);
specstart=0.9999*ones(1,numObs);
%[stapleV, sen, spec, Sall] = staple(rateMatM,iterlim, single(senstart), single(specstart));
[stapleV, sen, spec, Sall] = gpuStaple(rateMatM,iterlim, single(senstart), single(specstart));

staple3M = zeros(siz);
staple3M(indV) = stapleV;
newStrMask3M = staple3M >= probCutoff;

planC = maskToCERRStructure(newStrMask3M, 1, scanNum, structureName, planC);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function out=locfit_all(varargin)

% Smoothing noisy data using Local Regression and Likelihood.
%
% This is a combination of the locfit and predict functions
%

% Minimal input validation
if nargin < 1
   error( 'At least one input argument required' );
end

predict_args = {};

locfit_args = varargin{1};

if nargin==2
predict_args = varargin{2};
end;

fit = locfit( locfit_args{:} );

predict_out = predict( fit, predict_args{:} );

out = {fit predict_out};

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculating average precision using 11 point
% averaging.
%
% We first linearly interpolate the PR curves,
% which turns out to be more robust to the
% number of points sampled on the PR cureve.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ap = calcAP(rec, prec, interval)
if nargin < 3
    interval = 0:0.1:1;
end

% linear interpolation
[rec,ii] = sort(rec);
prec = prec(ii);
[rec,ii] = unique(rec);
prec = prec(ii);
Rq = 0:0.01:1;
Pq = interp1(rec,prec,Rq);
Pq(isnan(Pq)) = 0;
prec = Pq;
rec = Rq;

ap=0;
for t=interval
    p=max(prec(rec>=t));
    if isempty(p)
        p=0;
    end
    ap=ap+p/numel(interval);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [groupedImages] = groupImagesByColor(inputImages,groupVector,varargin)
	% Groups images by color based on input grouping vector.
	% Biafra Ahanonu
	% started: 2014.01.03 [19:13:01]
	% inputs
		% inputImages - [x y nImages]
		% groupVector - vector with same number of elements as images in inputImages
	% outputs
		%

	[groupedImages] = ciapkg.image.groupImagesByColor(inputImages,groupVector,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


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

%%
clear; close all; clc;

%% Syntax
% |C=doubleContraction(A,B,subContract);|

%% Description
% This function performs a double contraction of two tensors |A| and |B| to
% produce the output tensor |C|. The tensors are either both 3x3 second
% order tensors or |A| is a 4th order 3x3x3x3 tensor. The optional 3rd
% input |subContract| defines the indices to contract over for 4th order
% tensors e.g. |subContract=[1 2]| would contract over the first 2 indices.
%

%% Examples
%

%% Double contraction of a 2nd-order and a 2nd-order tensor

%%
% Creating an example tensor. Here a tensor |A| with all zeros except for
% the 3,3 component is created. Next a tensor |B| which, through double
% contraction, can "pick out" this scalar value is created. Finally the
% double contraction is performed illustrating this process.

%Create example tensor A
A=zeros(3,3); A(3,3)=pi; %Zeros with Azz=A33=pi
R=euler2DCM([0.25*pi 0.25*pi 0]); %Rotation matrix
A=R*A*R' %Rotated version of A

% Create 'structure' tensor example B
b=[0 0 1]*R.'; %Rotated Z vector
B=b.'*b %Rotated "structure" tensor for vector b

%%
% Double contraction using |doubleContraction|
s=doubleContraction(A,B)

%% Double contraction of a 4th-order and a 2nd-order tensor

%%
% Example for Hooke's law of linear elasticity

%Creating the 4th order linear elastic stiffness tensor
I=eye(3,3); %The 2nd order identity tensor
II1=dyadicProduct(I,I,1); %4th order base tensor 1
II3=dyadicProduct(I,I,3); %4th order base tensor 3

%Parameters for Hooke's law
mu=1; %The shear modulus
lambda=5; %The lambda lame parameter

%Compose stiffness tensor
C=lambda.*II1+2.*mu.*II3;

%Create strain tensor
E=rand(3,3); E=E*E'

%%
% Double contraction using |doubleContraction|
subContract=[3 4];
S=doubleContraction(C,E,subContract)

%% Using symbolic variables

%%
% Example for Hooke's law of linear elasticity

%Creating the 4th order linear elastic stiffness tensor
I=eye(3,3); %The 2nd order identity tensor
II1=dyadicProduct(I,I,1); %4th order base tensor 1
II3=dyadicProduct(I,I,3); %4th order base tensor 3

%Parameters for Hooke's law
try
    syms mu lambda
catch
    mu=1; %The shear modulus
    lambda=5; %The lambda lame parameter
end
%Compose stiffness tensor
C=lambda.*II1+2.*mu.*II3;

%Display Voigt mapped version
c=voigtMap(C)

%Create strain tensor
syms E1_1 E2_2 E3_3 E1_2 E1_3 E2_3;
E=[E1_1 E1_2 E1_3;...
   E1_2 E2_2 E2_3;...
   E1_3 E2_3 E3_3]

%%
% Double contraction using |doubleContraction|

subContract=[3 4];
S=doubleContraction(C,E,subContract)

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function rot3daxes(h,as,l,ang,axSym,killLines)
% rot3daxes preforms a simple rotation of an axis set created by make3daxes.
%
% rot3daxes(h,as,l,ang) rotates the surface objects with
% handles in array (as) in figure (h) by an angle (ang) about the axis
% defined by array (l).  The handles are assumed to be of surface objects
% generated by the function make3daxes and the angle is assumed to be in
% degrees.
%
% rot3daxes(...,axSym) plots the previous axis positions with linesytle
% axSym.
%
% rot3daxes(...,killLines) removes previous line objects from the figure,
% when killLines is set to true.
%
% Note: This function is intended to be called by animEuler.m
%
% Examples:
%     %rotate axes 45 degrees about z axis.
%     rot3daxes(h,as,[0,0,1],45,'--')

% Written by Dmitry Savransky 29 April 2009

if ~exist('axSym','var') || isempty(axSym)
    axSym = '--';
end
if ~exist('killLines','var')
    killLines = false;
end
%set figure and remove any previous line objects if requested
figure(h);

if killLines
    oldObjs = findobj(h,'Type','line');
    if ~isempty(oldObjs), delete(oldObjs); end
end

%plot orig axes
hold on
col = ['r','g','b'];
for j=1:3
    a = get(as(j),{'Xdata','Ydata','Zdata'});
    plot3(mean(a{1},2),mean(a{2},2),mean(a{3},2),[col(j),axSym],'LineWidth',2);
end

for j=1:20
    for i=1:3
        rotate(as(i),l,ang/20,[0,0,0])
    end
    pause(0.1);
end
% axis equal;
% axis([-1 1 -1 1 -1 1])
% view(45,15)
% grid on

author: Sable

repo: mcbench-benchmarks


text:
function subRoPS =  subRoPSFunc(projNeighbor,binSize)
%get the distribution matrix
neighbNum = length(projNeighbor);
distrMatrix = zeros(binSize,binSize);
minX = min(projNeighbor(:,1));
stepX = (max(projNeighbor(:,1))-minX)/binSize;
minY = min(projNeighbor(:,2));
stepY = (max(projNeighbor(:,2))-minY)/binSize;

if stepX==0 || stepY==0
    subRoPS = [0,0,0,0,0];
    return;
end

for k=1:neighbNum
    idxX = ceil((projNeighbor(k,1) - minX)/stepX);
    idxY = ceil((projNeighbor(k,2) - minY)/stepY);
    if idxX>binSize      idxX = binSize;  end
    if idxX<1                idxX = 1;            end
    if idxY>binSize     idxY = binSize;   end
    if idxY<1                idxY = 1;            end
    distrMatrix(idxX,idxY) = distrMatrix(idxX,idxY)+1;
end
distrMatrix = distrMatrix/neighbNum;%normalization
%calculate the moment of this distribution matrix
meanX = 0;
meanY = 0;
pde = 0;
for idxX = 1:binSize
    for idxY = 1:binSize
        meanX = meanX+idxX*distrMatrix(idxX,idxY);
        meanY = meanY+idxY*distrMatrix(idxX,idxY);
        if distrMatrix(idxX,idxY)>0
            pde = pde - distrMatrix(idxX,idxY)*log2(distrMatrix(idxX,idxY));
        end
    end
end
u11 = 0;
u21 = 0;
u12 = 0;
u22 = 0;

for idxX = 1:binSize
    for idxY = 1:binSize
        u11 = u11+(idxX-meanX)*(idxY-meanY)*distrMatrix(idxX,idxY);
        u21 = u21+(idxX-meanX)^2*(idxY-meanY)*distrMatrix(idxX,idxY);
        u12 = u12+(idxX-meanX)*(idxY-meanY)^2*distrMatrix(idxX,idxY);
        u22 = u22+(idxX-meanX)^2*(idxY-meanY)^2*distrMatrix(idxX,idxY);
    end
end
 subRoPS = [u11,u21,u12,u22,pde];

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
% io_loadspec_niimrs.m
% Georg Oeltzschner, Johns Hopkins University 2021
%
% USAGE:
% out = io_loadspec_niimrs(filename);
%
% DESCRIPTION:
% Reads in MRS data stored according to the NIfTI MRS format.
% See the specification under
% https://docs.google.com/document/d/1tC4ugzGUPLoqHRGrWvOcGCuCh_Dogx_uu0cxKub0EsM/edit
%
% io_loadspec_niimrs outputs the data in structure format, with fields
% corresponding to time scale, fids, frequency scale, spectra, and header
% fields containing information about the acquisition. The resulting matlab
% structure can be operated on by the other functions in this MRS toolbox.
%
% This function is currently work-in-progress and is being tested on more
% and more datasets. Please contact the FID-A developers for help with a
% particular NIfTI-MRS file that you might encounter problems with when
% using this function.
%
% Currently, this function is limited to single-voxel MRS data. It is
% planned to develop it for full compatibility with 2D and 3D multi-voxel
% and spectroscopic imaging data.
%
% DEPENDENCIES:
% This function requires the dcm2nii toolbox (Xiangrui Li) to be on the
% MATLAB path
% https://github.com/xiangruili/dicm2nii
%
% INPUTS:
% filename   = filename of NIfTI MRS file (*.nii) to load.
%
% OUTPUTS:
% out        = Input dataset in FID-A structure format.

function out = io_loadspec_niimrs(filename)

% Read in the data using the dicm2nii toolbox
% (https://github.com/xiangruili/dicm2nii)
try
    nii = nii_tool('load', filename);
catch ME
    switch ME.identifier
        case 'MATLAB:UndefinedFunction'
            error(['Cannot find the function ''nii_tool.m''.' ...
                ' Please ensure that you have downloaded the required', ...
                ' dcm2nii toolbox (https://github.com/xiangruili/dicm2nii)', ...
                ' and added it to your MATLAB path.']);
        otherwise
            rethrow(ME);
    end
end

% Extract the header and header extensions
hdr = nii.hdr;
hdr_ext = jsondecode(nii.ext.edata_decoded);

% Extract the time-domain data
fids = nii.img;

% Extract spectrometer frequency and dwell time
f0 = hdr_ext.SpectrometerFrequency;
dt = hdr.pixdim(5);
sw = 1/dt;

% Specify dimensions
% In NIfTI MRS, the three spatial dimensions and the time dimension occupy
% fixed indices in the (maximum) 7-D array
dims.x = 1;
dims.y = 2;
dims.z = 3;
dims.t = 4;

% There are some pre-defined dimension names according to the FID-A
% convention. These dimensions may or may not be stored in the NIfTI MRS
% header, so we'll initialize them as 0.
dims.coils = 0;
dims.averages = 0;
dims.subSpecs = 0;
dims.extras = 0;

% The NIfTI MRS standard reserves the remaining 3 dimensions, which are
% then explicitly specified in the JSON header extension fields dim_5,
% dim_6 and dim_7.
if isfield(hdr_ext, 'dim_5')
    dim_number = 5;
    switch hdr_ext.dim_5
        case 'DIM_COIL'
            dims.coils      = dim_number;
        case 'DIM_DYN'
            dims.averages   = dim_number;
        case 'DIM_INDIRECT_0'
            dims.extras     = dim_number;
        case 'DIM_INDIRECT_1'
            dims.extras     = dim_number;
        case 'DIM_INDIRECT_2'
            dims.extras     = dim_number;
        case 'DIM_PHASE_CYCLE'
            dims.extras     = dim_number;
        case 'DIM_EDIT'
            dims.subSpecs   = dim_number;
        case 'DIM_MEAS'
            dims.extras     = dim_number;
        case 'DIM_USER_0'
            dims.extras     = dim_number;
        case 'DIM_USER_1'
            dims.extras     = dim_number;
        case 'DIM_USER_2'
            dims.extras     = dim_number;
        case 'DIM_ISIS'
            dims.subSpecs   = dim_number;
        otherwise
            error('Unknown dimension value specified in dim_5: %s', hdr_ext.dim_5);
    end

end
if isfield(hdr_ext, 'dim_6')
        dim_number = 6;
    switch hdr_ext.dim_6
        case 'DIM_COIL'
            dims.coils      = dim_number;
        case 'DIM_DYN'
            dims.averages   = dim_number;
        case 'DIM_INDIRECT_0'
            dims.extras     = dim_number;
        case 'DIM_INDIRECT_1'
            dims.extras     = dim_number;
        case 'DIM_INDIRECT_2'
            dims.extras     = dim_number;
        case 'DIM_PHASE_CYCLE'
            dims.extras     = dim_number;
        case 'DIM_EDIT'
            dims.subSpecs   = dim_number;
        case 'DIM_MEAS'
            dims.extras     = dim_number;
        case 'DIM_USER_0'
            dims.extras     = dim_number;
        case 'DIM_USER_1'
            dims.extras     = dim_number;
        case 'DIM_USER_2'
            dims.extras     = dim_number;
        case 'DIM_ISIS'
            dims.subSpecs   = dim_number;
        otherwise
            error('Unknown dimension value specified in dim_6: %s', hdr_ext.dim_6);
    end
end
if isfield(hdr_ext, 'dim_7')
        dim_number = 7;
    switch hdr_ext.dim_7
        case 'DIM_COIL'
            dims.coils      = dim_number;
        case 'DIM_DYN'
            dims.averages   = dim_number;
        case 'DIM_INDIRECT_0'
            dims.extras     = dim_number;
        case 'DIM_INDIRECT_1'
            dims.extras     = dim_number;
        case 'DIM_INDIRECT_2'
            dims.extras     = dim_number;
        case 'DIM_PHASE_CYCLE'
            dims.extras     = dim_number;
        case 'DIM_EDIT'
            dims.subSpecs   = dim_number;
        case 'DIM_MEAS'
            dims.extras     = dim_number;
        case 'DIM_USER_0'
            dims.extras     = dim_number;
        case 'DIM_USER_1'
            dims.extras     = dim_number;
        case 'DIM_USER_2'
            dims.extras     = dim_number;
        case 'DIM_ISIS'
            dims.subSpecs   = dim_number;
        otherwise
            error('Unknown dimension value specified in dim_7: %s', hdr_ext.dim_7);
    end
end

% Parse the NIfTI hdr.dim field:
allDims = hdr.dim(2:end); % all dimensions (including singletons)

% Find the number of points
nPts = allDims(dims.t);

% Find the number of averages.  'averages' will specify the current number
% of averages in the dataset as it is processed, which may be subject to
% change.  'rawAverages' will specify the original number of acquired
% averages in the dataset, which is unchangeable.
if dims.subSpecs ~= 0
    if dims.averages ~= 0
        averages = allDims(dims.averages)*allDims(dims.subSpecs);
        rawAverages = averages;
    else
        averages = allDims(dims.subSpecs);
        rawAverages = 1;
    end
else
    if dims.averages ~= 0
        averages = allDims(dims.averages);
        rawAverages = averages;
    else
        averages = 1;
        rawAverages = 1;
    end
end

% FIND THE NUMBER OF SUBSPECS
% 'subspecs' will specify the current number of subspectra in the dataset
% as it is processed, which may be subject to change. 'rawSubspecs' will
% specify the original number of acquired  subspectra in the dataset, which
% is unchangeable.
if dims.subSpecs ~=0
    subspecs = allDims(dims.subSpecs);
    rawSubspecs = subspecs;
else
    subspecs = 1;
    rawSubspecs = subspecs;
end

% ORDERING THE DATA AND DIMENSIONS
% The FID-A array ordering conventions differ from the NIfTI MRS
% convention. Most importantly, FID-A is primarily tailored towards
% single-voxel data. We will start designing this function towards this
% purpose, and later adapt the formalism proposed in the csi_mod branch of
% the FID-A GitHub repository.

if allDims(1)*allDims(2)*allDims(3) == 1 % x=y=z=1
    dims.x = 0;
    dims.y = 0;
    dims.z = 0;
    fids = squeeze(fids);

    %Now that we've indexed the dimensions of the data array, we now need to
    %permute it so that the order of the dimensions is standardized:  we want
    %the order to be as follows:
    %   1) time domain data.
    %   2) coils.
    %   3) averages.
    %   4) subSpecs.
    %   5) extras.

    % Adjust dimension indices for the fact that we have collapsed the
    % three spatial dimensions (which we don't need for SVS data)
    sqzDims = {};
    dimsFieldNames = fieldnames(dims);
    for rr = 1:length(dimsFieldNames)
        if dims.(dimsFieldNames{rr}) ~= 0
            % Subtract 3 (x, y, z) from the dimension indices
            dims.(dimsFieldNames{rr}) = dims.(dimsFieldNames{rr}) - 3;
            sqzDims{end+1} = dimsFieldNames{rr};
        end
    end

    if length(sqzDims)==5
        fids=permute(fids,[dims.t dims.coils dims.averages dims.subSpecs dims.extras]);
        dims.t=1;dims.coils=2;dims.averages=3;dims.subSpecs=4;dims.extras=5;
    elseif length(sqzDims)==4
        if dims.extras==0
            fids=permute(fids,[dims.t dims.coils dims.averages dims.subSpecs]);
            dims.t=1;dims.coils=2;dims.averages=3;dims.subSpecs=4;dims.extras=0;
        elseif dims.subSpecs==0
            fids=permute(fids,[dims.t dims.coils dims.averages dims.extras]);
            dims.t=1;dims.coils=2;dims.averages=3;dims.subSpecs=0;dims.extras=4;
        elseif dims.averages==0
            fids=permute(fids,[dims.t dims.coils dims.subSpecs dims.extras]);
            dims.t=1;dims.coils=2;dims;averages=0;dims.subSpecs=3;dims.extras=4;
        elseif dims.coils==0
            fids=permute(fids,[dims.t dims.averages dims.subSpecs dims.extras]);
            dims.t=1;dims.coils=0;dims.averages=2;dims.subSpecs=3;dims.extras=4;
        end
    elseif length(sqzDims)==3
        if dims.extras==0 && dims.subSpecs==0
            fids=permute(fids,[dims.t dims.coils dims.averages]);
            dims.t=1;dims.coils=2;dims.averages=3;dims.subSpecs=0;dims.extras=0;
        elseif dims.extras==0 && dims.averages==0
            fids=permute(fids,[dims.t dims.coils dims.subSpecs]);
            dims.t=1;dims.coils=2;dims.averages=0;dims.subSpecs=3;dims.extras=0;
        elseif dims.extras==0 && dims.coils==0
            fids=permute(fids,[dims.t dims.averages dims.subSpecs]);
            dims.t=1;dims.coils=0;dims.averages=2;dims.subSpecs=3;dims.extras=0;
        end
    elseif length(sqzDims)==2
        if dims.extras==0 && dims.subSpecs==0 && dims.averages==0
            fids=permute(fids,[dims.t dims.coils]);
            dims.t=1;dims.coils=2;dims.averages=0;dims.subSpecs=0;dims.extras=0;
        elseif dims.extras==0 && dims.subSpecs==0 && dims.coils==0
            fids=permute(fids,[dims.t dims.averages]);
            dims.t=1;dims.coils=0;dims.averages=2;dims.subSpecs=0;dims.extras=0;
        elseif dims.extras==0 && dims.averages==0 && dims.coils==0
            fids=permute(fids,[dims.t dims.subSpecs]);
            dims.t=1;dims.coils=0;dims.averages=0;dims.subSpecs=2;dims.extras=0;
        end
    elseif length(sqzDims)==1
        dims.t=1;dims.coils=0;dims.averages=0;dims.subSpecs=0;dims.extras=0;
    end

    %Now get the size of the data array:
    sz=size(fids);

    %Compared to NIfTI MRS, FID-A needs the conjugate
    fids = conj(fids);

    %Now take fft of time domain to get fid:
    specs=fftshift(ifft(fids,[],dims.t),dims.t);

end

% Fill in additional FID-A format variables
% Nucleus (new field)
out.nucleus = hdr_ext.ResonantNucleus;
% Calculate B0 from spectrometer frequency depending on nucleus
% Gamma from Wikipedia article "Gyromagnetic ratio" (3 signif. digits)
for rr = 1:length(out.nucleus)
    switch out.nucleus{rr}
        case '1H'
            gamma = 42.577;
        case '2H'
            gamma = 6.536;
        case '3HE'
            gamma = -32.434;
        case '7LI'
            gamma = 16.546;
        case '13C'
            gamma = 10.708;
        case '19F'
            gamma = 40.052;
        case '23NA'
            gamma = 11.262;
        case '31P'
            gamma = 17.235;
        case '129XE'
            gamma = -11.777;
    end
    Bo(rr) = f0(rr) ./ gamma;
end

% Calculate t and ppm arrays using the calculated parameters:
f   =[(-sw/2) + (sw/(2*nPts)) : sw/(nPts) : (sw/2) - (sw/(2*nPts))];
ppm = -f / (Bo(1)*42.577);
ppm = ppm + 4.65;
t   = [0 : dt : (nPts-1)*dt];

% MANDATORY FIELDS
% Data & dimensions
out.fids = fids;
out.specs = specs;
out.sz = sz;
out.dims = dims;
out.Bo = Bo;
out.averages = averages;
out.rawAverages = rawAverages;
out.subspecs = subspecs;
out.rawSubspecs = rawSubspecs;

% Echo/repetition time
out.te = hdr_ext.EchoTime;
out.tr = hdr_ext.RepetitionTime;

% time and frequency axis
out.t   = t;
out.ppm = ppm;

% Dwell time & spectral width & field strength
out.spectralwidth = sw;
out.dwelltime = dt;
out.txfrq  = f0 * 10^6;
out.date = '';

% NIfTI-MRS-SPECIFIC FIELDS
% Save the NIfTI header
out.nii_mrs.hdr = hdr;
% Save the header extension
out.nii_mrs.hdr_ext = hdr_ext;
if isfield(hdr_ext, 'SequenceName')
    out.seq = hdr_ext.SequenceName;
end

%FILLING IN THE FLAGS
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
if out.dims.averages==0
    out.flags.averaged=1;
else
    out.flags.averaged=0;
end
if out.dims.coils==0
    out.flags.addedrcvrs=1;
else
    out.flags.addedrcvrs=0;
end
out.flags.subtracted=0;
out.flags.writtentotext=0;
out.flags.downsampled=0;
if out.dims.subSpecs==0
    out.flags.isFourSteps=0;
else
    out.flags.isFourSteps=(out.sz(out.dims.subSpecs)==4);
end

end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function ROTZ = ROTZ(theta)

ROTZ = [cos(theta) -sin(theta) 0;
        sin(theta) cos(theta) 0;
        0 0 1];
end

author: yorgoon

repo: minimum-snap-geometric-control


text:
classdef prtPreProcFunction < prtPreProc
    % prtPreProcFunction   Applies a function to observations
    %
    %   FUN = prtPreProcFunction creates a pre processing object that
    %   applies a specified function to the observations.
    %
    %   The function is set as a function handle in the property
    %   "transformationFunction" The default is @(x)(x), which does
    %   nothing.
    %
    %   A prtPreProcFunction object also inherits all properties and
    %   methods from the prtAction class
    %
    %   Example:
    %
    %   dataSet = prtDataGenIris;       % Load a data set.
    %   dataSet = dataSet.retainFeatures(1:2);
    %   fun = prtPreProcFunction('transformationFunction',@(x)x.^2);
    %   fun = fun.train(dataSet);       % All prtAction's must be trained
    %   dataSetNew = fun.run(dataSet);  % Normalize the data
    %
    %   % Plot
    %   subplot(2,1,1); plot(dataSet);
    %   title('Original Dataset')
    %   subplot(2,1,2); plot(dataSetNew);
    %   title('x.^2 Dataset')
    %
    %   See Also: prtPreProc, prtPreProcPca, prtPreProcPls,
    %   prtPreProcHistEq, prtPreProcZeroMeanColumns, prtPreProcLda,
    %   prtPreProcZeroMeanRows, prtPreProcLogDisc, prtPreProcZmuv,
    %   prtPreProcMinMaxRows

    properties (SetAccess=private)
        name = 'Function'
        nameAbbreviation = 'FUN'
    end
    properties
        transformationFunction = @(x)x;
    end

    properties (Hidden)
        operateOnMatrix = false; % Set to true for faster operation, but be careful
        operateOnDataStream = false;
    end

    methods
        function self = prtPreProcFunction(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
    end

    methods (Access=protected,Hidden=true)

        function self = trainAction(self,ds) %#ok<INUSD>
            % Nothing to do here
        end

        function ds = runAction(self,ds)
            if self.operateOnDataStream
                ds = self.transformationFunction(ds);
            else
               ds.X = runActionFast(self,ds.X);
            end
        end

        function x = runActionFast(self,x)
            if self.operateOnDataStream
                error('Can''t use runActionFast with operateOnDataStream = true');
            end
            if self.operateOnMatrix
                x = feval(self.transformationFunction,x);

            else % operateOnRows
                xOut = repmat(self.transformationFunction(x(1,:)),size(x,1),1);
                for i = 2:size(x,1)
                    xOut(i,:) = self.transformationFunction(x(i,:));
                end
                x = xOut;

             end
        end
    end
    methods
        function self = set.transformationFunction(self,val)
            if isa(val,'char')
                % convert strings to function handles
                % might be a good idea to do this to avoid other bundling
                % the current workspace into the function used
                val = str2func(val);
            end
            assert(isa(val,'function_handle'), 'transformationFunction must be a function handle that accepts at least one input argument.');
            self.transformationFunction = val;
        end
    end

    methods (Hidden = true)

        function yOut = runStream(self,vector)
            % yOut = runStream(self,vector)
            yOut = self.transformationFunction(vector);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [idx, frame] = registration_target(frames, useGPU)

%% Parameters
[ly, lx, nFrames] = size(frames);
data = single(frames);

if nargin < 2
  useGPU = false;
end

%% Prepare common arrays
% Taper mask
[ys, xs] = ndgrid(1:ly, 1:lx);
ys = abs(ys - mean(ys(:)));
xs = abs(xs - mean(xs(:)));
mY      = max(ys(:)) - 4;
mX      = max(xs(:)) - 4;
slope   = 1.2; % was 2

maskMul = single(1./(1 + exp((ys - mY)/slope)) ./(1 + exp((xs - mX)/slope)));
maskOffset = mean(data(:,:,1))*(1 - maskMul);
% Smoothing filter in frequency domain
sigma = 0.76; %.76 with mask
hgx = exp(-(((0:lx-1) - fix(lx/2))/sigma).^2);
hgy = exp(-(((0:ly-1) - fix(ly/2))/sigma).^2);
hg = hgy'*hgx;
fhg = real(fftn(ifftshift(single(hg/sum(hg(:))))));

eps0 = single(1e-20);
if useGPU
  data = gpuArray(data);
%   peakCorrMatrix = zeros(nFrames, nFrames, 'single', 'gpuArray');
% else
%   peakCorrMatrix = zeros(nFrames, nFrames, 'single');
end

corrRanges = zeros(nFrames, 3, 'single');

[a, b] = ndgrid(1:nFrames, 1:nFrames);
data = fft2(bsxfun(@plus, maskOffset, bsxfun(@times, maskMul, data)));
data = data./(eps0 + abs(data));
cdata = conj(data);
for fi = 1:nFrames
  % compute correlation map of frame fi with every other frame
  cmap = bsxfun(@times, data(:,:,a(:,fi)).*cdata(:,:,b(:,fi)), fhg);
  cmap = real(ifft2(cmap));
  % find peak of each correlation map
  peaks = gather(max(reshape(cmap, ly*lx, nFrames), [], 1));
  % compute first quartile of peak correlations, ignoring result
  % from correlation of this frame with itself.
  corrRanges(fi,:) = prctile(peaks((1:nFrames) ~= fi), [25 50 75]);
end
% find the frame with the best correlations with other frames
[bestCorr, idx] = max(corrRanges(:,1));
% plot(corrRanges)
frame = frames(:,:,idx);
% peakCorrMatrix = gather(peakCorrMatrix);
% [mxi, idx] = max(prctile(peakCorrMatrix, 25, 1))
% [mnx, mni] = min(prctile(peakCorrMatrix, 25, 2));
% figure, imagesc(peakCorrMatrix);
% figure, plot(1:nFrames, peakCorrMatrix(mxi,:), 1:nFrames, peakCorrMatrix(mni,:))
end

author: cortex-lab


text:
function fig = visualize_apollo_estimation(problem,apollopath,info,showfig)
if nargin < 4
    showfig = false;
end

markerSize = 6;
markerSizeest = 6;
color_detected_keypoints = 'cyan';
color_est_keypoints      = 'yellow';

if size(problem.intrinsicsorg,1) == 1 || size(problem.intrinsicsorg,2) == 1
    problem.intrinsicsorg       = vec_K_to_mat_K(problem.intrinsicsorg);
end
fullshapefname              = sprintf('%s/cad_db_values.mat',apollopath);
load(fullshapefname);
edgesfname                  = sprintf('%s/edges.mat',apollopath);
load(edgesfname);

shapes                      = permute(cad_db_values,[2,3,1]);

if problem.K == 20
    shapeids                    = 1:4:77;
elseif problem.K == 10
    shapeids                    = 1:10;
elseif problem.K == 5
    shapeids                    = 1:5;
else
    error('K can only be 10 or 20.')
end

shapes                      = shapes(:,:,shapeids);

R_est                       = info.R_est_org;
t_est                       = info.t_est_org;
shape_est                   = combine_shapes(shapes,info.c_est_org);
shape_est                   = R_est * shape_est + t_est;

K                           = problem.intrinsicsorg;
kpts2D                      = problem.kpts2D(1:2,:);
shape_homo                  = shape_est ./ shape_est(3,:);

kpts_est                    = K * shape_homo;
kpts_est                    = kpts_est(1:2,:);

imgname                     = sprintf('%s/apollo-train-set/images/%s.jpg',...
                                      apollopath,problem.imgname);
img                         = imread(imgname);

if showfig
    fig = figure;
else
    fig = figure('visible','off');
end
imshow(img,'Border','tight'); hold on;
% plot detected keypoints using squares
plot(kpts2D(1,:),kpts2D(2,:),'s',...
    'markersize',markerSize,...
    'MarkerFaceColor',color_detected_keypoints,...
    'MarkerEdgeColor',color_detected_keypoints);

hold on
plot(kpts_est(1,:),kpts_est(2,:),'o',...
        'markersize',markerSizeest,...
        'MarkerFaceColor',color_est_keypoints,...
        'MarkerEdgeColor',color_detected_keypoints);
hold on
nrEdges  = size(edges,1);
for i = 1:nrEdges
    linei = plot(kpts_est(1,edges(i,:)),kpts_est(2,edges(i,:)),'y-','lineWidth',2);
    linei.Color(4) = 0.6;
end

end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

%-----------------------------------------------------------------------
% This script demonstrates the use of function LINPROG on the basis of
% the C-VaR portfolio optimization problen given in:
% Uryasev, Rockafellar: Optimization of Conditional Value-at-Risk(1999)

clear
close all
clc
warning 'off'

tic

% load matrix data
% of annualized linear stock returns
load 'LinRetMat'

% M: number of samples
% N: number of assets
[M,N] = size(data);
% sample mean
mu = mean(data)';

% confidence level
beta = 0.99;

% number portfolios
NumPorts = 50;
% target returns
Rstar = min(mu):(max(mu)-min(mu))/(NumPorts-1):max(mu);
% objective function vector
f = objfCVaR(beta, M, N);
% constraints matrices
[A, b, Aeq, beq, lb, ub] = constCVaR(data, mu, Rstar(1));

CVaR    = zeros(NumPorts,1);
weights = zeros(N,NumPorts);

for i = 1:NumPorts
    beq(1) = -Rstar(i);
    [optimvar, CVaR(i)] = linprog(f, A, b, Aeq, beq, lb, ub);
    weights(:,i) = optimvar(1:N);
end

toc
warning 'on'

% Create figure
figure1 = figure('Color',[1 1 1]);
colormap('gray');

% plot efficient frontier
plot(CVaR,Rstar,'k','LineWidth',1.5)
legend('Mean-CVaR Frontier','Location','NorthWest')
set(gca,'xlim',[0 CVaR(end)],'ylim',[0 Rstar(end)])
title('Mean-CVaR Frontier','FontSize',17)
xlabel('$\beta$-CVaR','Interpreter','latex','FontName','Helvetia','FontSize',16)
ylabel('$\mu$','Interpreter','latex','FontName','Helvetia','FontSize',16,'Rotation',0)
ticks = get(gca,'YTick');
set(gca,'YTickLabel',[num2str(ticks'*100),repmat('%',length(ticks),1)])
ticks = get(gca,'XTick');
set(gca,'XTickLabel',[num2str(ticks'*100),repmat('%',length(ticks),1)])


text:
function [W_sorted, ind_sort, modularity, com_id, com_sizes] = gsp_sort_nodes(W, recursive, self_loops)
%GSP_SORT_NODES Sort nodes using the Louvain clustering method
%   Usage:  [Wout, ind, modularity, com_id, com_sizes] = gsp_sort_nodes(W)
%           [Wout, ind, modularity, com_id, com_sizes] = gsp_sort_nodes(W, recursive)
%           [Wout, ind, modularity, com_id, com_sizes] = gsp_sort_nodes(W, recursive, self_loops)
%
%   Inputs:
%         W         : Weigthed adjacency matrix. Can contain self-loops
%         recursive : Use recursive computation (default: 1)
%         self_loops: Take into account self loops (default: 0)
%
%   Outputs:
%         Wout      : Sorted weighted adjacency matrix
%         ind       : Indices used for sorting:       Wout = W(ind, ind)
%         modularity: Modularity of final clustering used for sorting
%         com_id    : Community ID of each node in the initial W
%         com_sizes : Number of nodes in each community
%
%   gsp_sort_nodes(W) sorts the nodes of the weighted adjacency matrix $W$
%   by clustering them according to the "Louvain" method.
%
%   By default the clustering is recursive and the level of recursion that
%   gives the maximum modularity is kept. The nodes of the graph are sorted
%   according to the clusters given.
%
%   The clustering method is using code given free online by Antoine
%   Scherrer available here:
%   https://perso.uclouvain.be/vincent.blondel/research/louvain.html
%
%   Example:::
%
%         G = gsp_nn_graph([[randn(100, 1)-1; (randn(100, 1)+7)], ...
%           [randn(50, 1); (randn(100,1)-4); (randn(50, 1)*2+8)]]);
%         [Wout, ind, modularity, com_id, com_sizes] = gsp_sort_nodes(G.W);
%         figure; imagesc(Wout); title('sorted adjacency matrix');
%         figure; gsp_plot_signal(G, ind);
%
%
%
%   See also: cluster_jl
%
%   References: blondel2008fast
%

% Author: Vassilis Kalofolias
% Date: June 2016

if nargin < 2
    % Use recursive clustering and pick the level with the highest modularity
    recursive = 1;
end
if nargin < 3
    % Don't use self loops by default
    self_loops = 0;
end
debug = 0;
verbose = 0;

COMTY = cluster_jl(W, recursive, self_loops, debug, verbose);

% Keep as best the splitting with highest modularity
ind_best_splitting = argmax(COMTY.MOD);

% Which node belongs to which community
com_id = COMTY.COM{ind_best_splitting};

% Sort according to the communities
[~, ind_sort] = sort(com_id);

modularity = COMTY.MOD(ind_best_splitting);
com_sizes = COMTY.SIZE{ind_best_splitting};

W_sorted = W(ind_sort, ind_sort);

repo: gspbox


text:
%% Machine Learning Online Class
%  Exercise 6 | Spam Classification with SVMs
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     gaussianKernel.m
%     dataset3Params.m
%     processEmail.m
%     emailFeatures.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% ==================== Part 1: Email Preprocessing ====================
%  To use an SVM to classify emails into Spam v.s. Non-Spam, you first need
%  to convert each email into a vector of features. In this part, you will
%  implement the preprocessing steps for each email. You should
%  complete the code in processEmail.m to produce a word indices vector
%  for a given email.

fprintf('\nPreprocessing sample email (emailSample1.txt)\n');

% Extract Features
file_contents = readFile('emailSample1.txt');
word_indices  = processEmail(file_contents);

% Print Stats
fprintf('Word Indices: \n');
fprintf(' %d', word_indices);
fprintf('\n\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ==================== Part 2: Feature Extraction ====================
%  Now, you will convert each email into a vector of features in R^n.
%  You should complete the code in emailFeatures.m to produce a feature
%  vector for a given email.

fprintf('\nExtracting features from sample email (emailSample1.txt)\n');

% Extract Features
file_contents = readFile('emailSample1.txt');
word_indices  = processEmail(file_contents);
features      = emailFeatures(word_indices);

% Print Stats
fprintf('Length of feature vector: %d\n', length(features));
fprintf('Number of non-zero entries: %d\n', sum(features > 0));

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 3: Train Linear SVM for Spam Classification ========
%  In this section, you will train a linear classifier to determine if an
%  email is Spam or Not-Spam.

% Load the Spam Email dataset
% You will have X, y in your environment
load('spamTrain.mat');

fprintf('\nTraining Linear SVM (Spam Classification)\n')
fprintf('(this may take 1 to 2 minutes) ...\n')

C = 0.1;
model = svmTrain(X, y, C, @linearKernel);

p = svmPredict(model, X);

fprintf('Training Accuracy: %f\n', mean(double(p == y)) * 100);

%% =================== Part 4: Test Spam Classification ================
%  After training the classifier, we can evaluate it on a test set. We have
%  included a test set in spamTest.mat

% Load the test dataset
% You will have Xtest, ytest in your environment
load('spamTest.mat');

fprintf('\nEvaluating the trained Linear SVM on a test set ...\n')

p = svmPredict(model, Xtest);

fprintf('Test Accuracy: %f\n', mean(double(p == ytest)) * 100);
pause;

%% ================= Part 5: Top Predictors of Spam ====================
%  Since the model we are training is a linear SVM, we can inspect the
%  weights learned by the model to understand better how it is determining
%  whether an email is spam or not. The following code finds the words with
%  the highest weights in the classifier. Informally, the classifier
%  'thinks' that these words are the most likely indicators of spam.
%

% Sort the weights and obtin the vocabulary list
[weight, idx] = sort(model.w, 'descend');
vocabList = getVocabList();

fprintf('\nTop predictors of spam: \n');
for i = 1:15
    fprintf(' %-15s (%f) \n', vocabList{idx(i)}, weight(i));
end

fprintf('\n\n');
fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% =================== Part 6: Try Your Own Emails =====================
%  Now that you've trained the spam classifier, you can use it on your own
%  emails! In the starter code, we have included spamSample1.txt,
%  spamSample2.txt, emailSample1.txt and emailSample2.txt as examples.
%  The following code reads in one of these emails and then uses your
%  learned SVM classifier to determine whether the email is Spam or
%  Not Spam

% Set the file to be read in (change this to spamSample2.txt,
% emailSample1.txt or emailSample2.txt to see different predictions on
% different emails types). Try your own emails as well!
filename = 'spamSample1.txt';

% Read and predict
file_contents = readFile(filename);
word_indices  = processEmail(file_contents);
x             = emailFeatures(word_indices);
p = svmPredict(model, x);

fprintf('\nProcessed %s\n\nSpam Classification: %d\n', filename, p);
fprintf('(1 indicates spam, 0 indicates not spam)\n\n');

author: Ayatans

repo: Machine-Learning-homework


text:
function success = test_mp_opt_model(verbose, exit_on_fail)
%TEST_MP_OPT_MODEL  Run all MP-Opt-Model tests.
%   TEST_MP_OPT_MODEL
%   TEST_MP_OPT_MODEL(VERBOSE)
%   TEST_MP_OPT_MODEL(VERBOSE, EXIT_ON_FAIL)
%   SUCCESS = TEST_MP_OPT_MODEL(...)
%
%   Runs all of the MP-Opt-Model tests. If VERBOSE is true (false by default),
%   it prints the details of the individual tests. If EXIT_ON_FAIL is true
%   (false by default), it will exit MATLAB or Octave with a status of 1
%   unless T_RUN_TESTS returns ALL_OK.
%
%   See also T_RUN_TESTS.

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

if nargin < 2
    exit_on_fail = 0;
    if nargin < 1
        verbose = 0;
    end
end

tests = {};

%% MP-Opt-Model tests
tests{end+1} = 't_have_fcn';
tests{end+1} = 't_nested_struct_copy';
tests{end+1} = 't_nleqs_master';
tests{end+1} = 't_pnes_master';
tests{end+1} = 't_qps_master';
tests{end+1} = 't_miqps_master';
tests{end+1} = 't_nlps_master';
tests{end+1} = 't_opt_model';
tests{end+1} = 't_om_solve_leqs';
tests{end+1} = 't_om_solve_nleqs';
tests{end+1} = 't_om_solve_pne';
tests{end+1} = 't_om_solve_qps';
tests{end+1} = 't_om_solve_miqps';
tests{end+1} = 't_om_solve_nlps';

%% run the tests
all_ok = t_run_tests( tests, verbose );

%% handle success/failure
if exit_on_fail && ~all_ok
    exit(1);
end
if nargout
    success = all_ok;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [ y, symm ] = cvx_s_upper_hankel( m, n, symm )

% CVX_S_UPPER_HANKEL Upper Hankel matrices.

c  = 0 : n - 1;
c  = c( ones( 1, m ), : );
r  = ( 0 : m - 1 )';
r  = r( :, ones( 1, n ) );
v  = abs( r + c ) + 1;
temp = v <= min( m, n );
y = sparse( v( temp ), r( temp ) + m * c( temp ) + 1, 1, min( m, n ), m * n );
symm = false;

% 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 h=som_plotplane(varargin)

%SOM_PLOTPLANE  Visualize the map prototype vectors as line graphs
%
% h=som_plotplane(lattice, msize, data, [color], [scaling], [pos])
% h=som_plotplane(topol, data, [color], [scaling], [pos])
%
%  som_plotplane('hexa',[5 5], rand(25,4), jet(25))
%  som_plotplane(sM, sM.codebook)
%
% Input and output arguments ([]'s are optional)
%  lattice   (string) grid 'hexa' or 'rect'
%  msize     (vector) size 1x2, defines the grid size
%            (matrix) size Mx2, defines explicit coordinates: in
%             this case the first argument does not matter
%  topol     (struct) map or topology struct
%  data      (matrix) Mxd matrix, M=prod(msize)
%  [color]   (matrix) size Mx3, gives an individual color for each graph
%            (string) ColorSpec gives the same color for each
%             graph, default is 'k' (black)
%  [scaling] (string) 'on' or 'off', default is 'on'
%  [pos]     (vector) 1x2 vector that determines translation.
%             Default is no translation.
%
%  h         (vector) the object handles for the LINE objects
%
%  If scaling is set on, the data will be linearly scaled in each
%  unit so that min and max values span from lower to upper edge
%  in each unit. If scaling is 'off', the proper scaling is left to
%  the user: values in range [-.5,.5] will be plotted within the limits of the
%  unit while values exceeding this range will be out of the unit.
%  Axis are set as in SOM_CPLANE.
%
% For more help, try 'type som_plotplane' or check out online documentation.
% See also SOM_PLANE, SOM_PIEPLANE, SOM_BARPLANE

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_plotplane
%
% PURPOSE
%
% Visualizes the map prototype vectors as line graph
%
% SYNTAX
%
%  h = som_plotplane(topol, data)
%  h = som_plotplane(lattice, msize, data)
%  h = som_plotplane(..., color)
%  h = som_plotplane(..., color, scaling)
%  h = som_plotplane(..., color, scaling, pos)
%
% DESCRIPTION
%
%  Visualizes the map prototype vectors as line graph
%
% KNOWN BUGS
%
%  It is not possible to specify explicit coordinates for map
%  consistig of just one unit as then the msize is interpreted as
%  map size.
%
% FEATURES
%
%  - the colors are fixed: changing colormap in the figure (see
%    COLORMAP) will not affect the coloring of the plots
%
% REQUIRED INPUT ARGUMENTS
%
% lattice  The basic topology
%
%   (string) 'hexa' or 'rect' positions the plots according to hexagonal or
%            rectangular map lattice.
%
% msize    The size of the map grid
%
%   (vector) [n1 n2] vector defines the map size (height n1 units, width n2
%            units, total M=n1 x n2 units). The units will be placed to their
%            topological locations in order to form a uniform hexagonal or
%            rectangular grid.
%   (matrix) Mx2 matrix defines arbitary coordinates for the M units.
%            In this case the argument 'lattice' has no effect.
%
% topol    Topology of the map grid
%
%   (struct) map or topology struct from which the topology is taken
%
% data     The data to be visualized
%
%   (matrix) Mxd matrix of data vectors.
%
% OPTIONAL INPUT ARGUMENTS
%
% If unspecified or given empty values ('' or []), default values
% will be used for optional input arguments.
%
% color    The color of the plots
%
%    (string) Matlab's ColorSpec (see help plot) string gives the same color
%             for each line.
%
%    (matrix) Mx3 matrix assigns an RGB color determined by the Nth row of
%             the matrix to the Nth plot.
%
%    (vector) 1x3 RGB vector gives the same color for each line.
%
% scaling  The data scaling mode
%
%    (string) 'on or 'off': if scaling is set on, the data will be
%             linearly scaled in each unit so that min and max values span from
%             lower to upper edge in each unit. If scaling is 'off', the proper
%             scaling is left to the user: values in range [-.5,.5] will be plotted
%             within the limits of the unit while values exceeding this
%             range will be out of the unit.
%
% pos      Position of the origin
%
%    (vector) This is meant for drawing the plane in arbitary location in a
%             figure. Note the operation: if this argument is given, the
%             axis limits setting part in the routine is skipped and the limits
%             setting will be left to be done by MATLAB's
%             defaults. By default no translation is done.
%
% OUTPUT ARGUMENTS
%
%  h (scalar)  Handle to the created patch object
%
% OBJECT TAG
%
% Object property 'Tag' is set to 'planePlot'.
%
% EXAMPLES
%
% %%% Create the data and make a map
%
% data=rand(1000,20); map=som_make(data);
%
% %%% Create a 'gray' colormap that has 64 levels
%
% color_map=gray(64);
%
% %%% Draw plots using red color
%
% som_plotplane(map, map.codebook,'r');
%
% %%% Calculate hits on the map and calculate colors so that
%     black = min. number of hits and white = max. number of hits
%
% hit=som_hits(map,data); color=som_normcolor(hit(:),color_map);
%
% %%% Draw plots again. Now the gray level indicates the number of hits to
%     each node
%
% som_plotplane(map, map.codebook, color);
%
% SEE ALSO
%
% som_cplane     Visualize a 2D component plane, u-matrix or color plane
% som_barplane   Visualize the map prototype vectors as bar diagrams.
% som_pieplane   Visualize the map prototype vectors as pie charts

% Copyright (c) 1999-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta Johan 160799 juuso 151199 070600

%%% Init & Check arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nargin, lattice, msize, data, color, scaling, pos] = vis_planeGetArgs(varargin{:});
error(nargchk(3, 5, nargin));  % check no. of input args is correct

s=0.8; % size of plot

if nargin < 6 || isempty(pos)
  pos=NaN;
end

if nargin < 5 || isempty(scaling)
  scaling='on';
elseif ~vis_valuetype(scaling,{'string'}) || ...
      ~any(strcmp(scaling,{'on','off'})),
  error('Scaling should be string ''on'' or ''off''.');
end

l=size(data,2);

if ~isnumeric(msize) || ndims(msize) ~= 2 || size(msize,2)~=2,
  error('msize has to be 1x2 grid size vector or a Nx2 coordinate matrix.');
elseif size(msize,1) == 1,
   xdim=msize(2);
   ydim=msize(1);
   N=xdim*ydim;
   y=repmat(repmat([1:ydim]',xdim,1),1,l);
   x=reshape(repmat([1:xdim],ydim*l,1),l,N)';
else
   x=repmat(msize(:,1),1,l);y=repmat(msize(:,2),1,l);
   N=size(msize,1);
   lattice='rect';
   if isnan(pos),
      pos=[0 0];
   end
end

switch lattice
case {'hexa', 'rect'}
otherwise
  error(['Lattice' lattice ' not implemented!']);
end

if ~isnumeric(data) || size(data,1) ~= N
  error('Data matrix is invalid or has wrong size!');
end

if nargin < 4 || isempty(color),
  color='k';
elseif vis_valuetype(color, {'colorstyle',[N 3]}),
  if ischar(color) && strcmp(color,'none'),
    error('Colorstyle ''none'' not allowed in som_plotplane.');
  end
elseif vis_valuetype(color,{'1x3rgb'})
elseif ~vis_valuetype(color,{'nx3rgb',[N 3]},'all'),
  error('The color matrix has wrong size or contains invalid RGB values or colorstyle.');
end

[linesx, linesy]=vis_line(data,scaling);

%%%% Action %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Making the lattice.
% Command view([0 90]) shows the map in 2D properly oriented

switch lattice
case 'hexa'
  t=find(rem(y(:,1),2));                  % move even rows by .5
  x(t,:)=x(t,:)-.5;
  x=(x./s+linesx).*s+.5; y=(y./s+linesy).*s;      % scale with s
case 'rect'
  x=(x./s+linesx).*s; y=(y./s+linesy).*s;         % scale with s
end

%% Draw the map! ...

h_=plot(x',y');

if size(color,1) == 1
  set(h_,'Color',color);
else
  for i=1:N,
    set(h_(i,:),'Color',color(i,:));
  end
end

if ~isnan(pos)
  x=x+pos(1);y=y+pos(2);                    % move upper left corner
end                                         % to pos(1),pos(2)

%% Set axes properties

ax=gca;
vis_PlaneAxisProperties(ax, lattice, msize, pos);

%%% Build output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

set(h_,'Tag','planePlot');                % tag the lineobject

if nargout>0, h=h_; end                   % Set h only,
                                          % if there really is output

%% Subfuntion %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,y]=vis_line(data, scaling)

s=size(data);
% normalization between [0,1] if scaling is on
if strcmp(scaling,'on')
  mn=repmat(min(data,[],2),1,s(2));
  mx=repmat(max(data,[],2),1,s(2));
  y=-((data-mn)./(mx-mn))+.5;
else                       % -sign is beacuse we do axis ij
  y=-data;
end

x=repmat(linspace(-.5, .5, size(data,2)), size(data,1),1);

author: ilarinieminen

repo: SOM-Toolbox


text:
function c = history(d)

%tstoolbox/@description/history
%   description/history
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

c = char(d.history);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function lr = LikelihoodR(ys,ar0,ma0,var0,ar1,ma1,var1)

%LIKELIHOODR Likelihood ratio between estimated ARMA models
%   lr = likelihoodR(y,ar0,ma0,var0,ar1,ma1,var1)
%   is the Likelihood ratio of the ARMA model ar0, ma0, var0
%   with respect to the ARMA model ar1, ma1, var1.
%
%   y can also be a matrix containing several segments of equal
%   length (segments in colums). The segments are considered to be
%   independent.
%
%   See also: LIKELIHOODR_MOD, KLINDEX_HAT.

% S. de Waele, March 2003.

lr = sum(KLIndex_hat(ys,ar0,ma0,var0)) - sum(KLIndex_hat(ys,ar1,ma1,var1));


text:
filename = 'BulkCaseFine';%'CantileverSquareNewFine';%'LshapeTriFine';%'CantileverSquareNew';'ArchTriFine';%'CantileverSquare';%'ArchTriFine';%'CantileverSquare';%'Lshape';%'LshapeTriFine';%'CantileverSquareSmall';%'';%'Lshape';'CantileverSquare';'Lshape';%'CantileverSquare';%'LshapeFine';%'Bridge_quad_coarse';%'Arch_quad_coarse';%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';%'Bridge';%'CantileverSquareSmall';
ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'P1';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-5;
constr_initial = 1e-5;

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 32;

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
%vademecumFileName = 'SuperEllipseQOptAnalytic';
%
% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';

kfrac = 2;
nsteps = 100;

plotting = false;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 1000;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function exact = p04_exact ( )

%*****************************************************************************80
%
%% P04_EXACT returns the exact integral for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = ( sqrt ( 32.0 ) / 3.0 ) * ( sqrt ( 27.0 ) - sqrt ( 8.0 ) - 1.0 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [outdata,outstate] = asr_process(data,srate,state,windowlen,lookahead,stepsize,maxdims,chunklen,usegpu)
% Processing function for the Artifact Subspace Reconstruction (ASR) method.
% [Data,State] = asr_process(Data,SamplingRate,State,WindowLength,LookAhead,StepSize,MaxDimensions,ChunkLength,UseGPU)
%
% This function is used to clean multi-channel signal using the ASR method. The required inputs are
% the data matrix, the sampling rate of the data, and the filter state (as initialized by
% asr_calibrate). If the data is used on successive chunks of data, the output state of the previous
% call to asr_process should be passed in.
%
% In:
%   Data : Chunk of data to process [#channels x #samples]. This is a chunk of data, assumed to be
%          a continuation of the data that was passed in during the last call to asr_process (if
%          any). The data should be *zero-mean* (e.g., high-pass filtered the same way as for
%          asr_calibrate).
%
%   SamplingRate : sampling rate of the data in Hz (e.g., 250.0)
%
%   State : initial filter state (determined by asr_calibrate or from previous call to asr_process)
%
%   WindowLength : Length of the statistcs window, in seconds (e.g., 0.5). This should not be much
%                  longer than the time scale over which artifacts persist, but the number of samples
%                  in the window should not be smaller than 1.5x the number of channels. Default: 0.5
%
%   LookAhead : Amount of look-ahead that the algorithm should use. Since the processing is causal,
%               the output signal will be delayed by this amount. This value is in seconds and should
%               be between 0 (no lookahead) and WindowLength/2 (optimal lookahead). The recommended
%               value is WindowLength/2. Default: WindowLength/2
%
%   StepSize : The statistics will be updated every this many samples. The larger this is, the faster
%              the algorithm will be. The value must not be larger than WindowLength*SamplingRate.
%              The minimum value is 1 (update for every sample) while a good value is 1/3 of a second.
%              Note that an update is always performed also on the first and last sample of the data
%              chunk. Default: 32
%
%   MaxDimensions : Maximum dimensionality of artifacts to remove. Up to this many dimensions (or up
%                   to this fraction of dimensions) can be removed for a given data segment. If the
%                   algorithm needs to tolerate extreme artifacts a higher value than the default
%                   may be used (the maximum fraction is 1.0). Default 0.66
%
%   ChunkLength : The length of the chunks to process, in samples. Larger chunks require more
%                 memory. Default: 50000
%
%   UseGPU : Whether to run on the GPU. This makes sense for offline processing if you have a a card
%            with enough memory and good double-precision performance (e.g., NVIDIA GTX Titan or
%            K20). Note that for this to work you need to have the Parallel Computing toolbox.
%            Default: false
%
% Out:
%   Data : cleaned data chunk (same length as input but delayed by LookAhead samples)
%
%   State : final filter state (can be passed in for subsequent calls)
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2012-08-31

% UC Copyright Notice
% This software is Copyright (C) 2013 The Regents of the University of California. All Rights Reserved.
%
% Permission to copy, modify, and distribute this software and its documentation for educational,
% research and non-profit purposes, without fee, and without a written agreement is hereby granted,
% provided that the above copyright notice, this paragraph and the following three paragraphs appear
% in all copies.
%
% Permission to make commercial use of this software may be obtained by contacting:
% Technology Transfer Office
% 9500 Gilman Drive, Mail Code 0910
% University of California
% La Jolla, CA 92093-0910
% (858) 534-5815
% invent@ucsd.edu
%
% This software program and documentation are copyrighted by The Regents of the University of
% California. The software program and documentation are supplied "as is", without any accompanying
% services from The Regents. The Regents does not warrant that the operation of the program will be
% uninterrupted or error-free. The end-user understands that the program was developed for research
% purposes and is advised not to rely exclusively on the program for any reason.
%
% IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
% SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF
% THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE. THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
% INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
% CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
% MODIFICATIONS.

if nargin < 4 || isempty(windowlen)
    windowlen = 0.5; end
windowlen = max(windowlen,1.5*size(data,1)/srate);
if nargin < 5 || isempty(lookahead)
    lookahead = windowlen/2; end
if nargin < 6 || isempty(stepsize)
    stepsize = 32; end
if nargin < 7 || isempty(maxdims)
    maxdims = 0.66; end
if nargin < 9 || isempty(usegpu)
    usegpu = false; end
if nargin < 8 || isempty(chunklen)
    chunklen = 50000; end
if maxdims < 1
    maxdims = round(size(data,1)*maxdims); end
if isempty(data)
    outdata = data; outstate = state; return; end

[C,S] = size(data);
N = round(windowlen*srate);
P = round(lookahead*srate);
[T,M,A,B] = deal(state.T,state.M,state.A,state.B);

% initialize prior filter state by extrapolating available data into the past (if necessary)
if isempty(state.carry)
    state.carry = repmat(2*data(:,1),1,P) - data(:,1+mod(((P+1):-1:2)-1,S)); end

data = [state.carry data];
data(~isfinite(data(:))) = 0;

% split up the total sample range into k chunks
splits = ceil(S/chunklen);
if splits > 1
    fprintf('Now cleaning data in %i blocks',splits); end

for i=1:splits
    range = 1+floor((i-1)*S/splits) : min(S,floor(i*S/splits));
    if ~isempty(range)
        % get spectrally shaped data X for statistics computation (range shifted by lookahead)
        X = double(data(:,range+P));
        if isempty(state.SOS)
            [X,state.iir] = filter(B,A,double(X),state.iir,2);
        else
            for s = 1:size(state.SOS,1)
                [X,state.Zi{s}] = filter(state.SOS(s,1:3),state.SOS(s,4:6),double(X),state.Zi{s},2); end
            X = X*state.G;
        end

        % move it to the GPU if applicable
        if usegpu && length(range) > 1000
            try X = gpuArray(X); catch,end; end
        % compute running mean covariance (assuming a zero-mean signal)
        [Xcov,state.cov] = moving_average(N,reshape(bsxfun(@times,reshape(X,1,C,[]),reshape(X,C,1,[])),C*C,[]),state.cov);
        % extract the subset of time points at which we intend to update
        update_at = min(stepsize:stepsize:(size(Xcov,2)+stepsize-1),size(Xcov,2));
        % if there is no previous R (from the end of the last chunk), we estimate it right at the first sample
        if isempty(state.last_R)
            update_at = [1 update_at];  %#ok<AGROW>
            state.last_R = eye(C);
        end
        Xcov = reshape(Xcov(:,update_at),C,C,[]);
        if usegpu
            Xcov = gather(Xcov); end
        % do the reconstruction in intervals of length stepsize (or shorter if at the end of a chunk)
        last_n = 0;
        for j=1:length(update_at)
            % do a PCA to find potential artifact components
            [V,D] = eig(Xcov(:,:,j));
            [D,order] = sort(reshape(diag(D),1,C)); V = V(:,order);
            % determine which components to keep (variance below directional threshold or not admissible for rejection)
            keep = D<sum((T*V).^2) | (1:C)<(C-maxdims);
            trivial = all(keep);
            % update the reconstruction matrix R (reconstruct artifact components using the mixing matrix)
            if ~trivial
                R = real(M*pinv(bsxfun(@times,keep',V'*M))*V');
            else
                R = eye(C);
            end
            % apply the reconstruction to intermediate samples (using raised-cosine blending)
            n = update_at(j);
            if ~trivial || ~state.last_trivial
                subrange = range((last_n+1):n);
                blend = (1-cos(pi*(1:(n-last_n))/(n-last_n)))/2;
                data(:,subrange) = bsxfun(@times,blend,R*data(:,subrange)) + bsxfun(@times,1-blend,state.last_R*data(:,subrange));
            end
            [last_n,state.last_R,state.last_trivial] = deal(n,R,trivial);
        end
    end
    if splits > 1
        fprintf('.'); end
end
if splits > 1
    fprintf('\n'); end

% carry the look-ahead portion of the data over to the state (for successive calls)
state.carry = [state.carry data(:,(end-P+1):end)];
state.carry = state.carry(:,(end-P+1):end);

% finalize outputs
outdata = data(:,1:(end-P));
if usegpu
    state.iir = gather(state.iir);
    state.cov = gather(state.cov);
end
outstate = state;

function [X,Zf] = moving_average(N,X,Zi)
% Run a moving-average filter along the second dimension of the data.
% [X,Zf] = moving_average(N,X,Zi)
%
% In:
%   N : filter length in samples
%   X : data matrix [#Channels x #Samples]
%   Zi : initial filter conditions (default: [])
%
% Out:
%   X : the filtered data
%   Zf : final filter conditions
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2012-01-10

if nargin <= 2 || isempty(Zi)
    Zi = zeros(size(X,1),N); end

% pre-pend initial state & get dimensions
Y = [Zi X]; M = size(Y,2);
% get alternating index vector (for additions & subtractions)
I = [1:M-N; 1+N:M];
% get sign vector (also alternating, and includes the scaling)
S = [-ones(1,M-N); ones(1,M-N)]/N;
% run moving average
X = cumsum(bsxfun(@times,Y(:,I(:)),S(:)'),2);
% read out result
X = X(:,2:2:end);

if nargout > 1
    Zf = [-(X(:,end)*N-Y(:,end-N+1)) Y(:,end-N+2:end)]; end

author: PatternRecognition

repo: OpenBMI


text:
function show(I,M)
%show - show content aware image resizing
%
% show(I,M) opens a figure window to allow interactive resizing of I
% in one direction based on the seam removal map M.
%
% See also removalMap.

% figure out if seams are horizontal or vertical
[h,w] = size(M);
horiz = all(sum(M,1)==sum(1:h))
vert = all(sum(M,2)==sum(1:w))
if ~xor(horiz,vert), error('seams not horizontal or vertical?!'); end

fh = figure; clf;
set(fh,'DoubleBuffer','on');
set(fh,'WindowButtonDownFcn',@mouseDown);
set(fh,'WindowButtonMotionFcn',@mouseMove);
set(fh,'WindowButtonUpFcn',@mouseUp);
set(fh,'WindowScrollWheelFcn',@mouseWheel);
% set(fh,'Interruptible','off');
% set(fh,'BusyAction','queue');
if horiz,
    ih = imshow(cat(1,I,I));
    title('horizontal seams');
else
    ih = imshow(cat(2,I,I));
    title('vertical seams');
end
set(ih,'CData',I);

if ~horiz,
    I = permute(I,[2,1,3]);
    M = M';
end

[h,w] = size(M);

n = 0; % how many seams we are removing

    % remove an additional dn seams
    function move(dn)
        n = n + dn;
        if n >= 0, % shrink image
            if horiz,
                set(ih,'CData',shrink(I,M,n));
            else
                set(ih,'CData',permute(shrink(I,M,n),[2,1,3]));
            end
        else % grow image
            if horiz,
                set(ih,'CData',expand(I,M,-n));
            else
                set(ih,'CData',permute(expand(I,M,-n),[2,1,3]));
            end
        end
        drawnow;
    end

p0 = []; % location of mouse down event
drag = false; % is the mouse being dragged?
moving = false; % are we moving the image?

    function mouseDown(src,evt)
        p0 = get(fh,'CurrentPoint');
        drag = true;
    end

    function mouseMove(src,evt)
        if ~drag, return; end
        if moving, return; end % serialize execution of this function
        moving = true;
        p = get(fh,'CurrentPoint');
        if horiz,
            move( p(2) - p0(2) );
        else
            move( p0(1) - p(1) );
        end
        p0 = p;
        moving = false;
    end

    function mouseWheel(src,evt)
        move( -evt.VerticalScrollCount * evt.VerticalScrollAmount );
    end

    function mouseUp(src,evt)
        drag = false;
    end

end


text:
function varargout = lu(varargin)
%LU   LU factorization of a SPHEREFUN.
%
%   This is not allowed and returns an error.  This function exists so that
%   the error message is meaningful to a SPHEREFUN user.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

error('CHEBFUN:SPHEREFUN:CUMPROD:notSupported',...
        'LU factorization of a SPHEREFUN is not supported.');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function model = mvuOptimise(model, display, iters)

% MVUOPTIMISE Optimise an MVU model.
% FORMAT
% DESC optimises a maximum variance unfolding model.
% ARG model : the model to be optimised.
% RETURN model : the optimised model.
%
% SEEALSO : mvuCreate, modelOptimise
%
% COPYRIGHT : Neil D. Lawrence, 2009

% MLTOOLS

if(any(any(isnan(model.Y))))
  error('Cannot run MVU when missing data is present.');
end

[X, details] = mvu(distance(model.Y'), model.k, 'solver', model.solver);

model.X = X(1:1:model.q,:)';
model.lambda = details.D/sum(details.D);

function D = distance(Y)

  D = sqrt(dist2(Y', Y'));
return

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function DCM = spm_dcm_erp_data(DCM,ERP)
% prepares structures for forward model(EEG, MEG and LFP)
% FORMAT DCM = spm_dcm_erp_data(DCM,ERP)
% DCM  - DCM structure
% ERP  - switch to average over trials (default)
%
% requires
%
%    DCM.xY.Dfile        - data file
%    DCM.options.trials  - trial codes
%    DCM.options.Tdcm    - Peri-stimulus time window
%    DCM.options.D       - Down-sampling
%    DCM.options.han     - Hanning
%    DCM.options.h       - Order of (DCT) detrending
%
% sets
%    DCM.xY.modality - 'MEG','EEG' or 'LFP'
%    DCM.xY.Time     - Time [ms] data
%    DCM.xY.pst      - Time [ms] of down-sampled data
%    DCM.xY.dt       - sampling in seconds (s)
%    DCM.xY.y        - cell array of trial-specific response {[ns x nc]}
%    DCM.xY.It       - Indices of (ns) time bins
%    DCM.xY.Ic       - Indices of (nc) good channels
%    DCM.xY.name     - names of (nc) channels
%    DCM.xY.scale    - scalefactor applied to raw data
%    DCM.xY.coor2D   - 2D coordinates for plotting
%    DCM.xY.X0       - (DCT) confounds
%    DCM.xY.R        - Residual forming matrix (with hanning)
%    DCM.xY.Hz       - Frequency bins (for Wavelet transform)
%    DCM.options.h
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_erp_data.m 7749 2019-12-05 17:05:46Z guillaume $

% Set defaults and Get D filename
%--------------------------------------------------------------------------
try
    Dfile = DCM.xY.Dfile;
catch
    errordlg('Please specify data and trials');
    error('')
end

% options
%==========================================================================

% order of drift terms
%--------------------------------------------------------------------------
if nargin < 2, ERP = 1;                         end
try, han  = DCM.options.han; catch, han   = 0;  end
try, h    = DCM.options.h;   catch, h     = 1;  end
try
    DT    = DCM.options.D;
catch
    errordlg('Please specify down sampling');
    error('')
end
try
    trial = DCM.options.trials;
catch
    errordlg('please specify trials');
    error('')
end

% load D
%--------------------------------------------------------------------------
try
    D = spm_eeg_load(Dfile);
catch
    try
        [dum,f]      = fileparts(Dfile);
        D            = spm_eeg_load(f);
        DCM.xY.Dfile = fullfile(pwd,f);
    catch
        try
            [f,p]        = uigetfile('*.mat','please select data file');
            name         = fullfile(p,f);
            D            = spm_eeg_load(name);
            DCM.xY.Dfile = fullfile(name);
        catch
            warndlg([Dfile ' could not be found'])
            return
        end
    end
end

% get time-frequency data if appropriate
%--------------------------------------------------------------------------
if isequal(D.transformtype, 'TF')
    DCM = spm_dcm_ind_data(DCM);
    return;
end

% indices of EEG channel (excluding bad channels) and peristimulus times
%--------------------------------------------------------------------------
if ~isfield(DCM.xY, 'modality')
    [mod, list] = modality(D, 0, 1);
    if isequal(mod, 'Multimodal')
        qstr    = 'Only one modality can be modelled at a time';
        if numel(list) < 4

            % Pretty dialog box
            %--------------------------------------------------------------
            options = [];
            options.Default = list{1};
            options.Interpreter = 'none';
            DCM.xY.modality = questdlg(qstr, 'Select modality', list{:}, options);

        else

            % can accommodate more buttons
            %--------------------------------------------------------------
            ind = menu(qstr, list);
            DCM.xY.modality = list{ind};
        end
    else
        DCM.xY.modality = mod;
    end
end

% good channels
%--------------------------------------------------------------------------
Ic = D.indchantype(DCM.xY.modality,'GOOD');
if isempty(Ic)
    warndlg('No good channels in these data');
    return
end

Nc            = length(Ic);               % number of channels
DCM.xY.name   = D.chanlabels(Ic);         % channel names
DCM.xY.Ic     = Ic;                       % channel indices
DCM.xY.Time   = time(D, [], 'ms');        % PST (ms)
DCM.xY.dt     = 1/D.fsample;              % time bins
DCM.xY.coor2D = D.coor2D(Ic);             % coordinates (topographic)

% time window
%--------------------------------------------------------------------------
try

    % time window and bins for modelling
    %----------------------------------------------------------------------
    T1          = DCM.options.Tdcm(1);
    T2          = DCM.options.Tdcm(2);
    [dum, T1]   = min(abs(DCM.xY.Time - T1));
    [dum, T2]   = min(abs(DCM.xY.Time - T2));

    % Time [ms] of down-sampled data
    %----------------------------------------------------------------------
    It          = (T1:DT:T2)';
    Ns          = length(It);                % number of samples
    DCM.xY.pst  = DCM.xY.Time(It);           % Down-sampled PST
    DCM.xY.dt   = DT/D.fsample;              % sampling in seconds
    DCM.xY.It   = It;                        % Indices of time bins

catch
    errordlg('Please specify time window');
    error('')
end

% confounds - DCT:
%--------------------------------------------------------------------------
if h == 0
    X0 = sparse(Ns,1);
else
    X0 = spm_dctmtx(Ns,h);
end
R      = speye(Ns) - X0*X0';

% hanning (omit second residualization for very long time-series)
%--------------------------------------------------------------------------
if han
    if Ns < 2048
        R = R*sparse(diag(spm_hanning(Ns)))*R;
    else
        R = sparse(diag(spm_hanning(Ns)))*R;
    end
end

% get trial averages - ERP
%--------------------------------------------------------------------------
cond  = D.condlist;
for i = 1:length(trial)

    % trial indices
    %----------------------------------------------------------------------
    c            = D.indtrial(cond(trial(i)), 'GOOD');
    Nt           = length(c);
    DCM.xY.nt(i) = Nt;

    % ERP
    %----------------------------------------------------------------------
    if ERP
        Y     = zeros(Ns,Nc);
        for j = 1:Nt
            Y = Y + R*D(Ic,It,c(j))';
        end
        DCM.xY.y{i} = Y/Nt;

    % all trials
    %----------------------------------------------------------------------
    else
        Y     = zeros(Ns,Nc,Nt);
        for j = 1:Nt
            Y(:,:,j) = R*D(Ic,It,c(j))';
        end
        DCM.xY.y{i} = Y;
    end

end

% condition units of measurement
%--------------------------------------------------------------------------
DCM.xY.code  = cond(trial);
DCM.xY.scale = 1;
DCM.xY.X0    = X0;
DCM.xY.R     = R;

author: spm