text:
% POP_SELECTCOMPS - Display components with button to visualize their
%                  properties and label them for rejection.
% Usage:
%       >> OUTEEG = pop_selectcomps( INEEG, compnum );
%
% Inputs:
%   INEEG    - Input dataset
%   compnum  - vector of component numbers
%
% Output:
%   OUTEEG - Output dataset with updated rejected components
%
% Note:
%   if the function POP_REJCOMP is ran prior to this function, some
%   fields of the EEG datasets will be present and the current function
%   will have some more button active to tune up the automatic rejection.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: POP_PROP, EEGLAB

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad

function [EEG, com] = pop_selectcomps( EEG, compnum, fig );

COLREJ = '[1 0.6 0.6]';
COLACC = '[0.75 1 0.75]';
PLOTPERFIG = 35;

com = '';
if nargin < 1
	help pop_selectcomps;
	return;
end;

if nargin < 2
    uilist = { { 'style' 'text' 'string' 'Components to plot:' } ...
               { 'style' 'edit' 'string'  ['1:' int2str(size(EEG.icaweights,1)) ] } ...
               {} ...
               { 'style' 'text' 'string' [ 'Note: in the next interface, click on buttons to see' char(10) ...
                                           'component properties and label them for rejection.' char(10) ...
                                           'To actually reject labelled components use menu item' char(10) ...
                                           '"Tools > Remove components" or use STUDY menus.' ] } };

    result = inputgui('uilist', uilist, 'geometry', { [1 1] 1 1 }, 'geomvert', [1 0.3 3], 'title', 'Reject comp. by map -- pop_selectcomps');
    if isempty(result), return; end
    compnum = eval( [ '[' result{1} ']' ]);

    if length(compnum) > PLOTPERFIG
        ButtonName=questdlg2(strvcat(['More than ' int2str(PLOTPERFIG) ' components so'],'this function will pop-up several windows'), ...
                             'Confirmation', 'Cancel', 'OK','OK');
        if ~isempty( strmatch(lower(ButtonName), 'cancel')), return; end
    end

end
fprintf('Drawing figure...\n');
currentfigtag = ['selcomp' num2str(rand)]; % generate a random figure tag

if length(compnum) > PLOTPERFIG
    for index = 1:PLOTPERFIG:length(compnum)
        pop_selectcomps(EEG, compnum([index:min(length(compnum),index+PLOTPERFIG-1)]));
    end

    com = [ 'pop_selectcomps(EEG, ' vararg2str(compnum) ');' ];
    return;
end

if isempty(EEG.reject.gcompreject)
	EEG.reject.gcompreject = zeros( size(EEG.icawinv,2));
end
try, icadefs;
catch,
	BACKCOLOR = [0.8 0.8 0.8];
	GUIBUTTONCOLOR   = [0.8 0.8 0.8];
end

% set up the figure
% -----------------
column =ceil(sqrt( length(compnum) ))+1;
rows = ceil(length(compnum)/column);
if ~exist('fig','var')
	figure('name', [ 'Reject components by map - pop_selectcomps() (dataset: ' EEG.setname ')'], 'tag', currentfigtag, ...
		   'numbertitle', 'off', 'color', BACKCOLOR);
	set(gcf,'MenuBar', 'none');
	pos = get(gcf,'Position');
	set(gcf,'Position', [pos(1) 20 800/7*column 600/5*rows*1.2]);
    incx = 120;
    incy = 110;
    sizewx = 100/column;
    if rows > 2
        sizewy = 90/rows;
	else
        sizewy = 80/rows;
    end
    pos = get(gca,'position'); % plot relative to current axes
	hh = gca;
	q = [pos(1) pos(2) 0 0];
	s = [pos(3) pos(4) pos(3) pos(4)]./100;
	axis off;
end

% figure rows and columns
% -----------------------
if EEG.nbchan > 64
    disp('More than 64 electrodes: electrode locations not shown');
    plotelec = 0;
else
    plotelec = 1;
end
count = 1;
for ri = compnum
	if exist('fig','var')
		button = findobj('parent', fig, 'tag', ['comp' num2str(ri)]);
		if isempty(button)
			error( 'pop_selectcomps(): figure does not contain the component button');
		end;
	else
		button = [];
	end;

	if isempty( button )
		% compute coordinates
		% -------------------
		X = mod(count-1, column)/column * incx-10;
        	Y = (rows-floor((count-1)/column))/rows * incy - sizewy*1.3;

		% plot the head
		% -------------
		if ~strcmp(get(gcf, 'tag'), currentfigtag);
		    figure(findobj('tag', currentfigtag));
		end
		ha = axes('Units','Normalized', 'Position',[X Y sizewx sizewy].*s+q);
		if plotelec
		    topoplot( EEG.icawinv(:,ri), EEG.chanlocs, 'verbose', ...
			      'off', 'chaninfo', EEG.chaninfo, 'numcontour', 8);
		else
		    topoplot( EEG.icawinv(:,ri), EEG.chanlocs, 'verbose', ...
			      'off', 'electrodes','off', 'chaninfo', EEG.chaninfo, 'numcontour', 8);
		end

		% labels
		% -------------
		if isfield(EEG.etc, 'ic_classification')
			classifiers = fieldnames(EEG.etc.ic_classification);
			if ~isempty(classifiers)
				if ~exist('classifier_name', 'var') || isempty(classifier_name)
					if any(strcmpi(classifiers, 'ICLabel'));
						classifier_name = 'ICLabel';
					else
						classifier_name = classifiers{1};
					end
				else
					classifier_name = classifiers{strcmpi(classifiers, classifier_name)};
				end
				if ri == compnum(1) && size(EEG.icawinv, 2) ...
						~= size(EEG.etc.ic_classification.(classifier_name).classifications, 1)
					warning(['The number of ICs do not match the number of IC classifications. This will result in incorrectly plotted labels. Please rerun ' classifier_name])
				end
				[prob, classind] = max(EEG.etc.ic_classification.(classifier_name).classifications(ri, :));
				t = title(sprintf('%s : %.1f%%', ...
					EEG.etc.ic_classification.(classifier_name).classes{classind}, ...
					prob*100));
				set(t, 'Position', get(t, 'Position') .* [1 -1.2 1])
			end
		end
		axis square;

		% plot the button
		% ---------------
         if ~strcmp(get(gcf, 'tag'), currentfigtag);
             figure(findobj('tag', currentfigtag));
         end
		button = uicontrol(gcf, 'Style', 'pushbutton', 'Units','Normalized', 'Position',...
                           [X Y+sizewy sizewx sizewy*0.18].*s+q, 'tag', ['comp' num2str(ri)]);
        command = sprintf('pop_prop( EEG, 0, %d, gcbo, { ''freqrange'', [1 50] });', ri);
		set( button, 'callback', command );
	end
	set( button, 'backgroundcolor', eval(fastif(EEG.reject.gcompreject(ri), COLREJ,COLACC)), 'string', int2str(ri));
	drawnow;
	count = count +1;
end

% draw the bottom button
% ----------------------
if ~exist('fig','var')
    if ~strcmp(get(gcf, 'tag'), currentfigtag);
        figure(findobj('tag', currentfigtag));
    end
	hh = uicontrol(gcf, 'Style', 'pushbutton', 'string', 'Cancel', 'Units','Normalized', 'backgroundcolor', GUIBUTTONCOLOR, ...
			'Position',[-10 -10  15 sizewy*0.25].*s+q, 'callback', 'close(gcf); fprintf(''Operation cancelled\n'')' );
	hh = uicontrol(gcf, 'Style', 'pushbutton', 'string', 'Set threhsolds', 'Units','Normalized', 'backgroundcolor', GUIBUTTONCOLOR, ...
			'Position',[10 -10  15 sizewy*0.25].*s+q, 'callback', 'pop_icathresh(EEG); pop_selectcomps( EEG, gcbf);' );
	if isempty( EEG.stats.compenta	), set(hh, 'enable', 'off'); end;
	hh = uicontrol(gcf, 'Style', 'pushbutton', 'string', 'See comp. stats', 'Units','Normalized', 'backgroundcolor', GUIBUTTONCOLOR, ...
			'Position',[30 -10  15 sizewy*0.25].*s+q, 'callback',  ' ' );
	if isempty( EEG.stats.compenta	), set(hh, 'enable', 'off'); end;
	hh = uicontrol(gcf, 'Style', 'pushbutton', 'string', 'Help', 'Units','Normalized', 'backgroundcolor', GUIBUTTONCOLOR, ...
			'Position',[70 -10  15 sizewy*0.25].*s+q, 'callback', 'pophelp(''pop_selectcomps'');' );
	command = '[ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); eegh(''[ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET);''); close(gcf)';
	hh = uicontrol(gcf, 'Style', 'pushbutton', 'string', 'OK', 'Units','Normalized', 'backgroundcolor', GUIBUTTONCOLOR, ...
			'Position',[90 -10  15 sizewy*0.25].*s+q, 'callback',  command);
			% sprintf(['eeg_global; if %d pop_rejepoch(%d, %d, find(EEG.reject.sigreject > 0), EEG.reject.elecreject, 0, 1);' ...
		    %		' end; pop_compproj(%d,%d,1); close(gcf); eeg_retrieve(%d); eeg_updatemenu; '], rejtrials, set_in, set_out, fastif(rejtrials, set_out, set_in), set_out, set_in));
end

com = [ 'pop_selectcomps(EEG, ' vararg2str(compnum) ');' ];
return;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function f = p00_f ( problem, n, x )

%*****************************************************************************80
%
%% P00_F evaluates the objective function for any problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem number.
%
%    Input, integer N, the number of arguments.
%
%    Input, real X(N), the arguments.
%
%    Output, real F(N), the function at each argument.
%
  if ( problem == 1 )
    f = p01_f ( n, x );
  elseif ( problem == 2 )
    f = p02_f ( n, x );
  elseif ( problem == 3 )
    f = p03_f ( n, x );
  elseif ( problem == 4 )
    f = p04_f ( n, x );
  elseif ( problem == 5 )
    f = p05_f ( n, x );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_F - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number PROBLEM = %d\n', problem );
    error ( 'P00_F - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [trainData,testData] = fetchStockData
% Copyright 2020 The MathWorks, Inc.
% helper function for downloading the data used for this example

% This example requires you to obtain an api-key
% To get your free api key go to:
% https://www.alphavantage.co/support/#api-key

%%%%%  Enter your api key here
apikey = '';
% Example:
% apikey = '6XRD9S9SYPOI0M96'; %note this is not a real key

%%%%%%
symbol = 'AMZN';
AMZN = fetchDaily(symbol,apikey);
symbol = 'GOOGL';
GOOGL = fetchDaily(symbol,apikey);
symbol = 'UNH';
UNH = fetchDaily(symbol,apikey);

dataC = synchronize(AMZN,GOOGL);
dataC = synchronize(dataC,UNH);

%setting dates for training and testing data sets
d1 = datetime('2004-08-19','InputFormat','yyyy-MM-dd');
d2 = datetime('2014-12-11','InputFormat','yyyy-MM-dd');
d3 = datetime('2014-12-12','InputFormat','yyyy-MM-dd');
d4 = datetime('2019-05-16','InputFormat','yyyy-MM-dd');

I1 = find(dataC.Close_Date == d1);
I2 = find(dataC.Close_Date == d2);
I3 = find(dataC.Close_Date == d3);
I4 = find(dataC.Close_Date == d4);

trainData = dataC{I1:I2,:};
trainData = array2table(trainData);
testData = dataC{I3:I4,:};
testData = array2table(testData);

author: matlab-deep-learning


text:
% Copyright 2017 Lime Microsystems Ltd.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%    http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%
% Just for fun, recover ASCII text from DPDCH with BCH(15,11) etc FEC
% Note real encoding uses Viterbi over several radio frames!  Viterbi is more efficient than block coding.
% also have not muxed DPCCH into DPDCH
function [y,evm,din,dout]=WCDMADLtxtMsgRead2( x )
	m=4;
	n=2^m-1; % also from bchpoly (1)
  ply=bchpoly(n);
  szp=size(ply);
  fecLevel=1;
  if fecLevel>szp(1)
    printf("Required FEClevel=%i too high for m=%i reduced to %i\n", fecLevel, m, szp(1));
    fecLevel=szp(1);
  end
	k=ply(fecLevel,2); % from bchpoly (2)
	t=ply(fecLevel,3); % from bchpoly (3)
  agc=1/abs(sqrt(sum(x.*conj(x))/length(x))); % agc as each channel is not a fixed level
  y=x*agc;
	% convert qpsk to bit stream
  qam=[1+i,1-i,-1+i,-1-i]/sqrt(2);
  bitsPerQAMSymbol=2;
	% convert qpsk to bit stream
%	y=WCDMAbpskDemod(x);
	evmref2=sum(abs(qam))/length(qam);
	[evm,loc]=min(abs(repmat(conj(qam'),1,length(y))-repmat(y,length(qam),1))); % decode QAM symbols
%	evm=evm./abs(qam(loc)); % some EVMs define relative to desired symbol
	evm/=evmref2; % some EVMs define relative to rms of QAM symbol table
	din=reshape(de2bi(loc-1,bitsPerQAMSymbol,2,"left-msb")',1,bitsPerQAMSymbol*length(loc)); % convert to binary stream
	bitsPerFrame=length(din);
	%convert bistream back to text
	chars6b=floor(floor(2*length(x)/n)*k/6); % DPDCH is complex x2 data
	symbols=floor(2*length(x)/n);
	framelen=symbols*n;
  y=din(1:framelen); % discarding padding
  y=matdeintrlv(y,n,symbols);% deinterleave
%	y=reshape(reshape(y,40,15)',1,600); % was 10 15 150
	% discard padding bits at end
	y=y(1:framelen);
	% convert into a nx? matrix
	l=size(y);
	y=reshape(y,l(2)/n,n); % group n bits for decoding
	y=bchdeco(y,k,t);
	y=reshape(y,1,k*l(2)/n);
	% convert back to binary, convert to 6 bit words, discarding padding
	dout=y(1:(6*chars6b));
	y=ASCII6dec(dout);
	printf('WCDMAtxtMsgDLRead[%s],BCH(%g,%g)\n', y, n, k );
end

% testing
% a=WCDMAtxtMsgDLWrite2("HelloZ",300);
% WCDMAtxtMsgDLRead2(a)

author: myriadrf


text:
function status = besa_save2SurfacePoint(file_path, file_name, ...
    channel_labels, channel_coords, fiducial_labels, fiducial_coords, ...
    headsurfacepoint_coords)
% besa_save2SurfacePoint writes the channel and fiducial coordinates into
% the file 'file_name'.
%
% Parameters:
%     [file_path]
%         Full path to the folder where the file should be saved.
%
%     [file_name]
%         The name of the file where the output should be written.
%
%     [channel_labels]
%         A cell (size: number of channels x 1) containing the channel
%         labels for the current data.
%
%     [channel_types]
%         An array (size: number of channels x 3) containing the channel
%         coordinates for the current data.
%
%     [fiducial_labels]
%         A cell (size: number of fiducials x 1) containing the fiducial
%         labels for the current data.
%
%     [fiducial_coords]
%         An array (size: number of fiducials x 3) containing the fiducials
%         coordinates for the current data.
%
%     [headsurfacepoint_coords]
%         Optional. An array (size: number of fiducials x 3) containing
%         additional head surface point coordinates for the current data.
%         Labels for head surface points will be created automatically,
%         starting with Sfh1.
%
% Return:
%     [status]
%         The status of the writing process. Still unused.
%

% Copyright (C) 2013, BESA GmbH
%
% File name: besa_save2SurfacePoint.m
%
% This file is part of MATLAB2BESA.
%
%    MATLAB2BESA 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.
%
%    MATLAB2BESA 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 MATLAB2BESA. If not, see <http://www.gnu.org/licenses/>.
%
% Author: Robert Spangler
% Created: 2013-12-22

n_FiducialLabels = size(fiducial_labels, 1);
n_FiducialCoords = size(fiducial_coords, 1);
n_ChannelLabels = size(channel_labels, 1);
n_ChannelCoords = size(channel_coords, 1);

% Get number of channels and fiducials
n_NumFiducials = n_FiducialLabels;
if n_FiducialLabels ~= n_FiducialCoords
    n_NumFiducials = min(n_FiducialLabels, n_FiducialCoords);
end;
n_NumChannels = n_ChannelLabels;
if n_ChannelLabels ~= n_ChannelCoords
    n_NumChannels = min(n_ChannelLabels, n_ChannelCoords);
end;

if exist('headsurfacepoint_coords', 'var');
    n_HSPCoords = size(headsurfacepoint_coords, 1);
else
    n_HSPCoords = 0;
end;

% Change to the directory where the data should be saved.
cd(file_path)

% Open file for writing.
fid = fopen(file_name, 'w');

% Write the fiducial labels and the coordinates to the file.
for i=1:n_NumFiducials
    fprintf(fid, '%s\t\t\t%2.6f\t%2.6f\t%2.6f\n', fiducial_labels{i}, ...
        fiducial_coords(i,1), fiducial_coords(i,2), fiducial_coords(i,3));
end;

% Write the channel labels and the coordinates to the file.
for i=1:n_NumChannels
    curChannelLabel = strrep(channel_labels{i}, 'Ele_', '');
    fprintf(fid, '%s\t\t\t%2.6f\t%2.6f\t%2.6f\n', curChannelLabel, ...
        channel_coords(i,1), channel_coords(i,2), channel_coords(i,3));
end

% Write additional head surface point labels and coordinates to the file.
for i=1:n_HSPCoords
    curHSPLabel = strcat('Sfh', num2str(i));
    fprintf(fid, '%s\t\t\t%2.6f\t%2.6f\t%2.6f\n', curHSPLabel, ...
        headsurfacepoint_coords(i,1), headsurfacepoint_coords(i,2), ...
        headsurfacepoint_coords(i,3));
end;

fclose(fid);
status = 1;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

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

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

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

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

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

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

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

    ifault = 1;
    value = 0.0;

  end

  return
end


text:
function [fileNameOut, extension] = lrsWriteRay(Q,modelName,vertexBool,param)
% Outputs a file for lrs to convert an V-representation (vertex / ray) of a
% polyhedron to a H-representation (half-space) via facet enumeration
%
% V-representation:
%     m is the number of input rows, each being a vertex, ray or line.
%     n is the number of input columns and d=n-1 is dimension of the input.
%     Each vertex is given in the form:
%
%     1   v_1   v_1  ... v_d
%
%     Each ray is given in the form:
%
%     0   r_1   r_2...   r_d
%
%     where  r_1  ...   r_d  is a point on the ray.
%
%     There must be at least one vertex in each file.
%     For bounded polyhedra there will be no rays entered.
%     The coefficients can be entered as integers or rationals in the format x/y.
%     An input line can be specified as a ray and then included in the linearity option (see below).
%
% INPUT
% Q             m x n integer matrix where each row is a variable and each column is a vertex or ray
% modelName     string giving the prefix of the *.ext file that will contain the vertex representation
%               It is assumed the file is pwd/*.ine, otherwise provide the full path.
%
% OPTIONAL INPUT
% vertexBool         n x 1 Boolean vector indicating which columns of Q are vertices
%                    By default, all columns of Q are assumed to be rays.
%
% param:             parameter structure with the following fields:
%    *.positivity:    if equals to 1, then positive orthant base
%    *.inequality:    if equals to 1, then represent as two inequalities rather than a single equality
%    *.shellScript:   if equals to 1, then lrs is run through a bash script
%    *.redund         if equals to 0, then remove redundant linear equalities

% Ronan Fleming 2021

[nMonomial,nRay]=size(Q);
if ~exist('vertexBool','var')
    fprintf('%s\n',['Assuming ' int2str(nRay) ' rays']);
    vertexBool = zeros(nRay,1);
end

if ~exist('param','var')
    param = struct();
end
if ~isfield(param,'positivity')
    param.positivity  = 0;
end
if ~isfield(param,'inequality')
    param.inequality  = 0;
end
if ~isfield(param,'shellScript')
    param.shellScript  = 0;
end
if ~isfield(param,'facetEnumeration')
    %assume vertex enumeration, unless specified that it is facet enumeration
    param.facetEnumeration  = 1;
end
if ~isfield(param,'redund')
    param.redund  = 1;
end

if param.redund ==0
    %remove all redundant halfspaces
    redundCmd = 'redund 0 0';
end

%lrs wants each row to be a ray so transpose
Q=Q';

extension = '.ext';
fileNameOut=[modelName extension];
fid=fopen(fileNameOut,'w');
fprintf(fid,'%s\n%s\n',modelName,'V-representation');

if any(~vertexBool)
    ind = find(vertexBool);
    if ~isempty(ind)
        fprintf(fid,'%s%s','linearity ',int2str(nnz(vertexBool)));
        for j = 1:length(ind)
            fprintf(fid,'%s%s',' ',int2str(ind(j)));
        end
        fprintf(fid,'\n');
    end
end

fprintf(fid,'%s\n','begin');

if ~any(vertexBool)
    fprintf(fid,'%s\n',[int2str(nRay+1) ' ' int2str(nMonomial+1) ' integer']);
    %there must be at least one vertex
    for r=1:nRay+1
        if r==1
            for d=1:nMonomial+1
                if d==1
                    fprintf(fid,'%s',[int2str(1)]);
                else
                    fprintf(fid,'%s',[' ' int2str(0)]);
                end
            end
        else
            for d=1:nMonomial+1
                if d==1
                    fprintf(fid,'%s',[int2str(0)]);
                else
                    fprintf(fid,'%s',[' ' int2str(Q(r-1,d-1))]);
                end
            end
        end
        fprintf(fid,'\n');
    end
else
    fprintf(fid,'%s\n',[int2str(nRay) ' ' int2str(nMonomial+1) ' integer']);
    for r=1:nRay
        for d=1:nMonomial+1
            if d==1
                fprintf(fid,'%s',[int2str(vertexBool(r)+0)]);
            else
                fprintf(fid,'%s',[' ' int2str(Q(r,d-1))]);
            end
        end
        fprintf(fid,'\n');
    end
end

%set of vertices specifying positive orthant
if param.positivity
    for i=1:nMonomial
        fprintf(fid,'%s',int2str(1));
        for j=1:nMonomial
            if i==j
                fprintf(fid,'%s',[' ' int2str(1)]);
            else
                fprintf(fid,'%s',[' ' int2str(0)]);
            end
        end
        fprintf(fid,'\n');
    end
end

fprintf(fid,'%s\n','end');
if param.redund==0
    %remove redundant halfspaces
    fprintf(fid,'%s\n',redundCmd);
end
fclose(fid);

author: opencobra

repo: cobratoolbox


text:
function camAAr2c = projICP(camAAr2c)

    global VMap;
    global NMap;
    global XYZcamBilateral;
    global NcamBilateral;
    global K;

    global NMapCam;
    global VMapCam;
    global Vku;

    for iterID = 1:10

        %% data association

        % transform the point cloud

        VMapCam = AngleAxisRotatePoint(camAAr2c, VMap) + repmat(camAAr2c(4:6)',1,size(VMap,2));
        NMapCam = AngleAxisRotatePoint(camAAr2c, NMap);

        % projective association

        px = round(K(1,1)*(VMapCam(1,:)./VMapCam(3,:)) + K(1,3));
        py = round(K(2,2)*(VMapCam(2,:)./VMapCam(3,:)) + K(2,3));

        isValid = (1<=px & px <= 640 & 1<=py & py<= 480);
        VMapCam = VMapCam(:,isValid);
        NMapCam = NMapCam(:,isValid);
        px = px(isValid);
        py = py(isValid);

        ind = sub2ind([480 640],py,px);

        %{
        cla(subplot(6,3,12))
        validMap = zeros(480,640);
        validMap(ind) = 1;
        imagesc(validMap);
        axis equal;
        axis tight;
        title(sprintf('Iteration %d: init valid Map',iterID))
        %}

        isValid = XYZcamBilateral(640*480*2+ind)~=0;
        VMapCam = VMapCam(:,isValid);
        NMapCam = NMapCam(:,isValid);
        ind = ind(isValid);

        % outlier rejection

        diffD = (VMapCam(1,:) - XYZcamBilateral(ind)).^2 + (VMapCam(2,:) - XYZcamBilateral(640*480+ind)).^2 + (VMapCam(3,:) - XYZcamBilateral(640*480*2+ind)).^2 ;
        dotProdN = sum([NcamBilateral(ind); NcamBilateral(640*480+ind); NcamBilateral(640*480*2+ind)] .* NMapCam,1);

        subplot(6,3,17)
        validMap = zeros(480,640);
        validMap(ind) = diffD;
        imagesc(validMap);
        axis equal;
        axis tight;
        title(sprintf('Iteration %d: diffD',iterID))

        subplot(6,3,18)
        validMap = zeros(480,640);
        validMap(ind) = dotProdN;
        imagesc(validMap);
        axis equal;
        axis tight;
        title(sprintf('Iteration %d: dotProdN',iterID))

        %isValid = (diffD<0.1^2) & (dotProdN > cos(pi/3));
        isValid = (diffD<1^2);

        VMapCam = VMapCam(:,isValid);
        NMapCam = NMapCam(:,isValid);
        ind = ind(isValid);
        Vku = [XYZcamBilateral(ind); XYZcamBilateral(640*480+ind); XYZcamBilateral(640*480*2+ind)];

        subplot(6,3,15)
        validMap = zeros(480,640);
        validMap(ind) = 1;
        imagesc(validMap);
        axis equal;
        axis tight;
        title(sprintf('Iteration %d: valid Map',iterID))

        %% optimization
        % objective function
        E = sum(( Vku - VMapCam ) .* NMapCam,1);

        subplot(6,3,16)
        validMap = zeros(480,640);
        validMap(ind) = E;
        imagesc(validMap);
        axis equal;
        axis tight;
        title(sprintf('Iteration %d: Distance Map',iterID))

        fprintf('initial error: #inliers = %.2f(%d/%d), sum = %f, mean = %f, median = %f\n', sum(isValid)/length(isValid), sum(isValid), length(isValid), sum(E.^2), mean(E.^2), median(E.^2));

        %options = optimset('Display','iter', 'Algorithm','levenberg-marquardt');
        options = optimset('display','off','Algorithm','levenberg-marquardt');
        [AA_gk, resnorm, residual, exitflag, output] = lsqnonlin(@residualFunction, [0 0 0 0 0 0],[],[],options);

        camAAr2c = cameraRt2AngleAxis(concatenateCameraRt(transformCameraRt(cameraAngleAxis2Rt(AA_gk)), cameraAngleAxis2Rt(camAAr2c)));
    end

end

function residuals = residualFunction(Tgk)
    global NMapCam;
    global VMapCam;
    global Vku;

    VkuTran = AngleAxisRotatePoint(Tgk, Vku) + repmat(Tgk(4:6)',1,size(Vku,2));

    residuals = ( VkuTran - VMapCam ) .* NMapCam;
end

%{
% for visualization

figure
image2show = zeros(480,640);
image2show(ind) = distance(1,:);
imagesc(image2show); axis equal; axis tight;
title('Vm(1) - Vd(1)')

figure
image2show = zeros(480,640);
image2show(ind) = distance(2,:);
imagesc(image2show); axis equal; axis tight;
title('Vm(2) - Vd(2)')

figure
image2show = zeros(480,640);
image2show(ind) = distance(3,:);
imagesc(image2show); axis equal; axis tight;
title('Vm(3) - Vd(3)')

figure
image2show = zeros(480,640);
image2show(ind) = sum(distance.^2,1);
imagesc(image2show); axis equal; axis tight;
title('(Vm - Vd)^2')

figure
image2show = zeros(480,640);
image2show(ind) = sum(distance .* NMapCam,1).^2;
imagesc(image2show); axis equal; axis tight;
title('E')

%}

% coarse to fine

author: jianxiongxiao

repo: ProfXkit


text:
function varargout = cpower(varargin)
%CPOWER Power of SDPVAR variable with convexity knowledge
%
% CPOWER is recommended if your goal is to obtain
% a convex model, since the function CPOWER is implemented
% as a so called nonlinear operator. (For p/q ==2 you can
% however just as well use the overloaded power)
%
% t = cpower(x,p/q)
%
% For negative p/q, the operator is convex.
% For positive p/q with p>q, the operator is convex.
% For positive p/q with p<q, the operator is concave.
%
% A domain constraint x>0 is automatically added if
% p/q not is an even integer.
%
% Note, the complexity of generating the conic representation
% of these variables are O(2^L) where L typically is the
% smallest integer such that 2^L >= min(p,q)

switch class(varargin{1})

    case 'double'
        varargout{1} = power(varargin{1},varargin{2});
    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.
        X = varargin{1};
        if isreal(X)
            dim = size(X);
            X = reshape(X,prod(dim),1);
            y = [];
            for i = 1:prod(dim)
                y = [y;yalmip('define',mfilename,extsubsref(X,i),varargin{2})];
            end
            y = reshape(y,dim);
            varargout{1} = y;
        else
            error('CPOWER can only be applied to real vectors.');
        end

    case 'char' % YALMIP send 'model' when it wants the epigraph or hypograph
        if isequal(varargin{1},'graph')
            t = varargin{2}; % Second arg is the extended operator variable
            X = varargin{3}; % Third arg and above are the args user used when defining t.
            p = varargin{4};

            if p>0
                [p,q] = rat(abs(p));
                F = pospower(X,t,p,q);
                if p>q
                    convexity = 'convex';
                    monotonicity = 'increasing';
                else
                    convexity = 'concave';
                    monotonicity = 'decreasing';
                end
            else
                [p,q] = rat(abs(p));
                F = negpower(X,t,p,q);
                convexity = 'convex';
                monotonicity = 'decreasing';
            end

            varargout{1} = F;
            varargout{2} = struct('convexity',convexity,'monotonicity',monotonicity,'definiteness','positive','model','graph');
            varargout{3} = X;
        end
    otherwise
end

function F = pospower(x,t,p,q)
if p>q
    l = ceil(log2(abs(p)));
    r = 2^l-p;
    y = [ones(r,1)*x;ones(q,1)*t;ones(2^l-r-q,1)];
    F = detset(x,y);
else
    l = ceil(log2(abs(q)));
    y = [ones(p,1)*x;ones(2^l-q,1)*t;ones(q-p,1)];
    F = detset(t,y);
end

function F = negpower(x,t,p,q)
l = ceil(log2(abs(p+q)));
p = abs(p);
q = abs(q);
y = [ones(2^l-p-q,1);ones(p,1)*x;ones(q,1)*t];
F = detset(1,y);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [class_name,rotation] = modelnet_extract_classname_and_rot_from_filename(filename)

s = strsplit(filename,{'.','/'});
class_name = s{end-4};
r = strsplit(s{end-1},'_');
rotation = str2double(r(end));

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function ft_realtime_packettimer(cfg)

% FT_REALTIME_PACKETTIMER can be used to time the rate at which data can be processed
%
% Use as
%   ft_realtime_packettimer(cfg)
% with the following configuration options
%   cfg.bcifun    = processing of the data (default = @bcifun_timer)
%   cfg.npackets  = the number of packets shown in one plot (default=1000)
%                     after reaching the end
%   cfg.saveplot  = if path is specified, first plot is saved (default=[]);
%   cfg.rellim = y limits of subplot 1 (default = [-100 100])
%
% SEE ALSO:
%   FT_REALTIME_PROCESS

% TO DO:
%   jitter in het binnenhalen van de data; scatterplot!
%   triggers sturen en herhalen (loop closen)
%   tijd schatten waarin matlab nog kan processen

% Copyright (C) 2009, Marcel van Gerven
%
% $Id$

if ~isfield(cfg,'bcifun'),    cfg.bcifun = @bcifun_timer; end
if ~isfield(cfg,'npackets'),  cfg.npackets = 10^2; end
if ~isfield(cfg,'rellim'),    cfg.rellim = [-1 1]; end
if ~isfield(cfg,'saveplot'),  cfg.saveplot= []; end

close all;
f1 = figure();
set(f1,'units','normalized','outerposition',[0 0 1 1]);

% reset persistent variables
cfg.bcifun();

try
  ft_realtime_process(cfg);
catch
  fprintf('%s\n',lasterr);
  close;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BCIFUN_TIMER plots real time versus packet time
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bcifun_timer(cfg,data)

persistent startsample;
persistent olddf;

persistent idx;

if nargin == 0
  idx = []; % reset persistent
  return;
end

if isempty(idx)

  tic;

  startsample = data.endsample;

  idx       = 1;
  olddf     = 0;

  return;

end

if mod(idx,cfg.npackets) == 1

  xl = [floor(idx/cfg.npackets)*cfg.npackets+1 (floor(idx/cfg.npackets)+1)*cfg.npackets];

  hold off
  plot(xl,[0 0],'k--');
  xlim(xl);
  ylim(cfg.rellim);
  xlabel('packet number');
  ylabel('delay (sample time - real time)');
  hold on;

end

cursample = data.endsample;
curtime   = toc;

smptime = (cursample - startsample)/data.fsample;
df = smptime - curtime;
if df < 0
  plot([idx-1 idx],[olddf df],'r');
else
  plot([idx-1 idx],[olddf df],'g');
end
olddf = df;
title(sprintf('packet timing for %g s blocks in %d channels at sample frequency %d; sample time = %g, real time = %g, delay = %g',...
  data.blocksize/data.fsample,length(data.label),data.fsample,smptime,curtime,smptime - curtime));

idx = idx + 1;

drawnow;

if idx==cfg.npackets && strcmp(cfg.saveplot,'yes')
  saveas(gcf,cfg.saveplot,'jpg');
end


text:
function stable_camva(a)
  % STABLE_CAMVA Act like camva without changing zoom (aka dolly zoom)
  %
  % Inputs:
  %   a  angle input to camva
  %
  % Note: This will mess up interactive view orbitting. So set up your
  % view first...
  %
  % See also: camva
  old_a = camva;
  old_t = camtarget;
  old_p = campos;
  camva(a);
  s = tan(a/2*pi/180)/tan(old_a/2*pi/180);
  p = (old_p-old_t)/s+old_t;
  campos(p);
end

author: alecjacobson

repo: gptoolbox


text:
function out = sig_prbs(n,m)
% Pseudo-random binary signal (PRBS)
%
%% Syntax
%  out = sig_prbs(n,m)
%
%% Description
% Function generates the pseudo-random binary signal (PRBS) with the uniform
% amplitude , see R. Iserman: Prozessidentifikation, 1974, Springer.
%
% Input:
% * n ... the length of the internal register
% * m ... the number of samples per the width of the narrowest
%         impulse of PRBS
%
% Output:
% * out ... PRBS signal values
%
% See also:
% sig_prs_minmax
%
%%

registry1=0:11;
registry2=[0 2 3 4 5 6 7 8 9 10 11 15 20 40 100];
store=ones(n,1);

out(1:m,1)=store(n);

for i=2:2^n-1
    storex=store(registry1(n))+store(registry2(n));
    for j=0:n-2
        store(n-j)=store(n-j-1);
    end
    store(1)=storex;
    if store(1)==2
        store(1)=0;
    end
    for j=(i-1)*m+1:m*i
        out(j,1)=2*store(n)-1;
    end
end

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function [ o, x, w ] = cn_jac_02_xiu ( n, alpha, beta )

%*****************************************************************************80
%
%% CN_JAC_02_XIU implements the Xiu rule for region CN_JAC.
%
%  Discussion:
%
%    The rule has order
%
%      O = N + 1.
%
%    The rule has precision P = 2.
%
%    CN is the cube [-1,+1]^N with the Jacobi (beta) weight function
%
%      w(alpha,beta;x) = product ( 1 <= i <= n ) (1-x(i))^beta (1+x(i))^alpha.
%
%    with -1 < alpha, -1 < beta.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dongbin Xiu,
%    Numerical integration formulas of degree two,
%    Applied Numerical Mathematics,
%    Volume 58, 2008, pages 1515-1520.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, real ALPHA, BETA, the parameters.
%    -1.0 < ALPHA, -1.0 < BETA.
%
%    Output, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  if ( alpha <= -1.0 )
    fprint ( 1, '\n' );
    fprint ( 1, 'CN_JAC_02_XIU - Fatal error!\n' );
    fprint ( 1, '  ALPHA <= -1.0\n' );
    error ( 'CN_JAC_02_XIU - Fatal error!' );
  end

  if ( beta <= -1.0 )
    fprint ( 1, '\n' );
    fprint ( 1, 'CN_JAC_02_XIU - Fatal error!\n' );
    fprint ( 1, '  BETA <= -1.0\n' );
    error ( 'CN_JAC_02_XIU - Fatal error!' );
  end

  o = n + 1;

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  for j = 1 : o

    i = 0;
    for r = 1 : floor ( n / 2 )
      arg = 2 * r * ( j - 1 ) * pi / ( n + 1 );
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * cos ( arg );
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * sin ( arg );
    end

    if ( i < n )
      i = i + 1;
      x(i,j) = r8_mop ( j - 1 );
    end

  end

  gamma0 = ( alpha + beta + 2.0 ) / 2.0;
  delta0 = ( alpha - beta ) / 2.0;
  c1 = 2.0 * ( alpha + 1.0 ) * ( beta + 1.0 ) / ( alpha + beta + 3.0 ) ...
    / ( alpha + beta + 2.0 );

  x(1:n,1:o) = ( sqrt ( gamma0 * c1 ) * x(1:n,1:o) - delta0 ) / gamma0;

  expon = 0;
  volume_1d = c1_jac_monomial_integral ( alpha, beta, expon );
  volume = volume_1d ^ n;

  w(1:o) = volume / o;

  return
end


text:
function A = int_area(X1,Y1,X2,Y2)
% Compute area of intersection:
%
% area = int_area(X1,Y1,X2,Y2)

X1 = double(X1);
Y1 = double(Y1);
X2 = double(X2);
Y2 = double(Y2);

max_res = 100;

min_x = min([X1(:); X2(:)]);
max_x = max([X1(:); X2(:)]);
min_y = min([Y1(:); Y2(:)]);
max_y = max([Y1(:); Y2(:)]);

X1 = X1-min_x;
X2 = X2-min_x;
Y1 = Y1-min_y;
Y2 = Y2-min_y;
max_x = max_x - min_x;
max_y = max_y - min_y;

max_dim = max(max_x,max_y);

ratio = 1;
% $$$   if max_dim > max_res
% $$$     ratio = max_res/max_dim;
% $$$     X1 = X1*ratio;
% $$$     X2 = X2*ratio;
% $$$     Y1 = Y1*ratio;
% $$$     Y2 = Y2*ratio;
% $$$     max_x = round(max_x*ratio);
% $$$     max_y = round(max_y*ratio);
% $$$   end

X1 = X1+1; X2 = X2+1;
Y1 = Y1+1; Y2 = Y2+1;
max_x = ceil(max_x+1);
max_y = ceil(max_y+1);

M1 = poly2mask(X1,Y1,max_y,max_x);
M2 = poly2mask(X2,Y2,max_y,max_x);

A = sum(sum(double(M1&M2)))/ratio;

author: CSAILVision

repo: LabelMeToolbox


text:
function pdf = normal_pdf ( x, a, b )

%*****************************************************************************80
%
%% NORMAL_PDF evaluates the Normal PDF.
%
%  Discussion:
%
%    The normal PDF is also known as the Gaussian PDF.
%
%  Formula:
%
%    PDF(X)(A,B) = EXP ( - 0.5 * ( ( X - A ) / B )^2 ) / SQRT ( 2 * PI * B^2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(), the argument of the PDF.
%
%    Input, real A, a parameter of the PDF.
%    A represents the mean value of the variables.
%
%    Input, real B, a parameter of the PDF.
%    B represents the standard deviation of the variables.
%
%    Output, real PDF(), the value of the PDF.
%
  pdf = exp ( - 0.5 * ( ( x - a ) / b ).^2 )  / sqrt ( 2.0 * pi * b^2 );

  return
end


text:
function b = r83_mv ( m, n, a, x )

%*****************************************************************************80
%
%% R83_MV multiplies an R83 matrix times an R8VEC.
%
%  Discussion:
%
%    The R83 storage format is used for a tridiagonal matrix.
%    The superdiagonal is stored in entries (1,2:N), the diagonal in
%    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
%    original matrix is "collapsed" vertically into the array.
%
%  Example:
%
%    Here is how an R83 matrix of order 5 would be stored:
%
%       *  A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54  *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(3,N), the R83 matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(M), the product A * x.
%
  b = zeros ( m, 1 );

  mn = min ( m, n );

  if ( n == 1 )
    b(1,1) = a(2,1) * x(1,1);
    if ( 1 < m )
      b(2,1) = a(3,1) * x(1,1);
    end
    return
  end

  b(1,1)      = a(2,1)       * x(1,1) ...
              + a(1,2)       * x(2,1);

  b(2:mn-1,1) = a(3,1:mn-2)' .* x(1:mn-2,1) ...
              + a(2,2:mn-1)' .* x(2:mn-1,1) ...
              + a(1,3:mn)'   .* x(3:mn,1);

  b(mn,1)     = a(3,mn-1)    * x(mn-1,1) ...
              + a(2,mn)      * x(mn,1);

  if ( n < m )
    b(mn+1,1) = b(mn+1,1) + a(3,mn) * x(mn,1);
  end

  if ( m < n )
    b(mn,1) = b(mn,1) + a(1,mn+1) * x(mn+1,1);
  end

  return
end


text:
function p = parent(tree,uid)
% XMLTREE/PARENT Parent Method
% FORMAT uid = parent(tree,uid)
%
% tree   - XMLTree object
% uid    - UID of the lonely child
% p      - UID of the parent ([] if root is the child)
%__________________________________________________________________________
%
% Return the uid of the parent of a node.
%__________________________________________________________________________
% Copyright (C) 2002-2011  http://www.artefact.tk/

% Guillaume Flandin
% $Id: parent.m 4460 2011-09-05 14:52:16Z guillaume $

p = tree.tree{uid}.parent;


text:
% Two-input NAND gate sizing (GP)
% Boyd, Kim, Patil, and Horowitz, "Digital circuit optimization
% via geometric programming"
% Written for CVX by Almir Mutapcic 02/08/06
% (a figure is generated)
%
% This is an example taken directly from the paper:
%
%   Digital circuit optimization via geometrical programming
%   by Boyd, Kim, Patil, and Horowitz
%   Operations Research 53(6): 899-932, 2005.
%
% Solves the problem of choosing device widths w_i for the given
% NAND2 gate in order to achive minimum Elmore delay for different
% gate transitions, subject to limits on the device widths,
% gate area, power, and so on. The problem is a GP:
%
%   minimize   D = max( D_1, ..., D_k )  for k transitions
%       s.t.   w_min <= w <= w_max
%              A <= Amax, etc.
%
% where variables are widths w.
%
% This code is specific to the NAND2 gate shown in figure 19
% (page 926) of the paper. All the constraints and the objective
% are hard-coded for this particular circuit.

%********************************************************************
% problem data and hard-coded GP specs (evaluate all transitions)
%********************************************************************
N = 4;       % number of devices
Cload = 12;  % load capacitance
Vdd = 1.5;   % voltage

% device specs
NMOS = struct('R',0.4831, 'Cdb',0.6, 'Csb',0.6, 'Cgb',1, 'Cgs',1);
PMOS = struct('R',2*0.4831, 'Cdb',0.6, 'Csb',0.6, 'Cgb',1, 'Cgs',1);

% maximum area and power specification
Amax = 24;
wmin = 1;

% varying parameters for the tradeoff curve
Npoints = 25;
Amax = linspace(5,45,Npoints);
Dopt = [];

disp('Generating the optimal tradeoff curve...')
need_sedumi = strncmpi(cvx_solver,'sdpt',4);
if need_sedumi,
    warning('This model does not converge with SDPT3... switching to SeDuMi.');
end
for k = 1:Npoints
    fprintf(1,'  Amax = %5.2f:', Amax(k));
    cvx_begin gp quiet
        if need_sedumi,
            cvx_solver sedumi
        end

        % device width variables
        variable w(N)

        % device specs
        device(1:2) = PMOS; device(3:4) = NMOS;

        for num = 1:N
            device(num).R   = device(num).R/w(num);
            device(num).Cdb = device(num).Cdb*w(num);
            device(num).Csb = device(num).Csb*w(num);
            device(num).Cgb = device(num).Cgb*w(num);
            device(num).Cgs = device(num).Cgs*w(num);
        end

        % capacitances
        C1 = sum([device(1:3).Cdb]) + Cload;
        C2 = device(3).Csb + device(4).Cdb;

        % input capacitances
        Cin_A = sum([ device([2 3]).Cgb ]) + sum([ device([2 3]).Cgs ]);
        Cin_B = sum([ device([1 4]).Cgb ]) + sum([ device([1 4]).Cgs ]);

        % resistances
        R = [device.R]';

        % area definition
        area = sum(w);

        % delays and dissipated energies for all six possible transitions
        % transition 1 is A: 1->1, B: 1->0, Z: 0->1
        D1 = R(1)*(C1 + C2);
        E1 = (C1 + C2)*Vdd^2/2;
        % transition 2 is A: 1->0, B: 1->1, Z: 0->1
        D2 = R(2)*C1;
        E2 = C1*Vdd^2/2;
        % transition 3 is A: 1->0, B: 1->0, Z: 0->1
        % D3 = C1*R(1)*R(2)/(R(1) + R(2)); % not a posynomial
        E3 = C1*Vdd^2/2;
        % transition 4 is A: 1->1, B: 0->1, Z: 1->0
        D4 = C1*R(3) + R(4)*(C1 + C2);
        E4 = (C1 + C2)*Vdd^2/2;
        % transition 5 is A: 0->1, B: 1->1, Z: 1->0
        D5 = C1*(R(3) + R(4));
        E5 = (C1 + C2)*Vdd^2/2;
        % transition 6 is A: 0->1, B: 0->1, Z: 1->0
        D6 = C1*R(3) + R(4)*(C1 + C2);
        E6 = (C1 + C2)*Vdd^2/2;

        % objective is the worst-case delay
        minimize( max( [D1 D2 D4] ) )
        subject to
            area <= Amax(k);
            w >= wmin;
    cvx_end
    % display and store computed values
    fprintf(1,' delay = %3.2f\n',cvx_optval);
    Dopt = [Dopt cvx_optval];
end

% plot the tradeoff curve
plot(Dopt,Amax);
xlabel('Dmin'); ylabel('Amax');
disp('Optimal tradeoff curve plotted.')

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function pass = test_kron(pref)
% Test the chebfun/kron() command.

if ( nargin < 1 )
    pref = chebfunpref();
end

tol = 200*pref.cheb2Prefs.chebfun2eps;

%% Kronecker products resulting in CHEBFUN2 objects

% Rank 1 chebfun2
f = chebfun(@(x) x.^2, pref);
g = chebfun(@(y) sin(y), pref);

h1 = chebfun2(@(x,y) x.^2.*sin(y));
h2 = chebfun2(@(x,y) y.^2.*sin(x));

pass(1) = (norm(h1 - kron(f', g)) < tol);
pass(2) = (norm(h2 - kron(f, g')) < tol);
pass(3) = (norm(h1 - g*f') < tol);
pass(4) = (norm(h2 - f*g') < tol);

% Different domain in x and y
d = [-2 pi -pi 2];
f = chebfun(@(x) x.^2, d(1:2), pref);
g = chebfun(@(y) sin(y), d(3:4), pref);

h1 = chebfun2(@(x,y) x.^2.*sin(y), d);
h2 = chebfun2(@(x,y) y.^2.*sin(x), [d(3:4) d(1:2)]);

pass(5) = ( norm(h1 - kron(f', g)) < tol);
pass(6) = (norm(h2 - kron(f, g')) < tol);
pass(7) = ( norm(h1 - g*f') < tol);
pass(8) = (norm(h2 - f*g') < tol);

% Quasimatrices and rank 4 chebfun2
x = chebfun('x', [-1 1], pref);
F = [1 x x.^2 x.^4];
G = [1 cos(x) sin(x) x.^5];

h1 = chebfun2(@(x,y) 1 + x.*cos(y) + x.^2.*sin(y) + x.^4.*y.^5);
h2 = chebfun2(@(x,y) 1 + y.*cos(x) + y.^2.*sin(x) + y.^4.*x.^5);

pass(9) = (norm(h1 - kron(F', G)) < tol);
pass(10) = (norm(h2 - kron(F, G')) < tol);
pass(11) = (norm(h1 - G*F') < tol);
pass(12) = (norm(h2 - F*G') < tol);

% Different domains and quasimatrices
x = chebfun('x',[-2 1], pref);
y = chebfun('x',[-1 1], pref);
d = [-2 1 -1 1];
F = [1 x x.^2 x.^4];
G = [1 cos(y) sin(y) y.^5];
h1 = chebfun2(@(x,y) 1 + x.*cos(y) + x.^2.*sin(y) + x.^4.*y.^5,d);
h2 = chebfun2(@(x,y) 1 + y.*cos(x) + y.^2.*sin(x) + y.^4.*x.^5,[d(3:4) d(1:2)]);

pass(13) = (norm(h1 - kron(F', G)) < tol);
pass(14) = (norm(h2 - kron(F, G')) < tol);
pass(15) = (norm(h1 - G*F') < tol);
pass(16) = (norm(h2 - F*G') < tol);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
 function ob = Gtomo2(chat)
%function ob = Gtomo2(chat)
%
%	create default "parent" 2D tomography object (structure)
%	with some basic methods common to all child tomography objects
%	chat: 0 to disable debugging messages, >=1 to display them
%
%	Copyright Mar 2001, Jeff Fessler, University of Michigan

if nargin < 1, chat = 0; end

%
%	the dimensions are usually [ob.nb * ob.na, ob.nx * ob.ny]
%	unless masked, in which case [ob.nb * ob.na, sum(mask(:))]
%	or if tranposed (in which case swapped)
%
ob.dims = [0 0];

ob.nx = 0;	% input image size
ob.ny = 0;
ob.nb = 0;	% output sinogram size (bins)
ob.na = 0;	% output sinogram size (angles)

ob.apower	= 1;		% array power, becomes 2 for G.^2
ob.scale	= 1;		% global scaling factor used in some children
ob.mask		= [];		% [nx,ny] logical array
ob.chat		= chat;
ob.index1	= [];		% row indices (unless transposed) e.g. G(3:9,:)
ob.index2	= [];		% col indices (unless transposed) e.g. G(:,3:9)
				% the empty default means *all* rows/cols
ob.is_transpose = false;
ob.is_masked	= false;	% set to 1 if G(:,mask(:))
%ob.version	= 1.0;

%ob = class(ob, 'Gtomo2');	% matlab inheritence sucks

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [zPred,PzPred,otherInfo]=reducedStateMeasPred(xPred,PPred,H)
%%REDUCEDSTATEMEASPRED Perform the measurement prediction part of the
%           measurement update step of the reduced state estimator. The
%           function reducedStateMeasUpdateWithPred can be used to complete
%           the measurement update. Separating the measurement prediction
%           step from the rest of the update step can make the creation of
%           multiple measurement association hypotheses from a single
%           target prediction more efficient. The full measurement update
%           function is reducedStateUpdate.
%
%INPUTS: xPred The xDimXnumComp predicted target states.
%        PPred The xDimXxDimXnumComp predicted state covariance matrices.
%            H The zDimXxDim measurement matrix. The measurement is modeled
%              as z=H*x+noise.
%
%OUTPUTS: zPred The zDimXnumComp measurement predictions from the filter.
%        PzPred The zDimXzDimXnumComp covariance matrix associated with the
%               values in zPred.
%     otherInfo A structure containing members of intermediate results of
%               this function that can be passed to KalmanUpdateWithPred
%               when updating with a measurement.
%
%The measurement prediction step in the measurement update step of the
%reduced state estimator is the same as that in the Kalman filter. Thus,
%this function just calls KalmanMeasPred.
%
%The filter is taken from [1]. See the comments to reducedStateUpdate for
%more information.
%
%EXAMPLE:
%With this example, we demonstrate that one gets the same result using
%reducedStateUpdate versus calling reducedStateMeasPred and then
%reducedStateMeasUpdateWithPred.
% xPred=[1e3;-2e3;100;200];
% MPred=[28,   3.5,    6,  8.5;
%      3.5,    23,  8.5,   11;
%        6,   8.5,   18, 13.5;
%      8.5,    11, 13.5,   13];
% DPred=[8, 0;
%        0, 8;
%        8, 0;
%        0, 8];
% PPred=MPred+(1/2)*(DPred*DPred');
% z=1e3*[-5.498856156296510;
%        1.199241491470584];
% R=eye(2);
% H=[0, 4, 9, 8;
%    6, 3, 0, 6];
% %The update in one step.
% [xUpdate,MUpdate,DUpdate,innov,Pzz,W]=reducedStateUpdate(xPred,PPred,MPred,DPred,z,R,H);
% %The update in two steps.
% [zPred,PzPred,otherInfo]=reducedStateMeasPred(xPred,PPred,H);
% [xUpdate1,MUpdate1,DUpdate1,innov1,Pzz1,W1]=reducedStateMeasUpdateWithPred(z,R,zPred,PzPred,otherInfo,MPred,DPred);
% %One will see that the one and two step updates agree.
% max(abs([xUpdate1-xUpdate;MUpdate1(:)-MUpdate(:);DUpdate1(:)-DUpdate(:);innov1(:)-innov;Pzz1(:)-Pzz(:);W1(:)-W(:)]))
%
%REFERENCES:
%[1] P. Mookerjee and F. Reifler, "Reduced state estimator for systems with
%    parametric inputs," IEEE Transactions on Aerospace and Electronic
%    Systems, vol. 40, no. 2, pp. 446-461, Apr. 2004.
%
%June 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

[zPred,PzPred,otherInfo]=KalmanMeasPred(xPred,PPred,H);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function h = montageKPM(arg)
% montageKPM is like the built-in montage, but assumes input is MxNxK or filenames
%
% Converts patches (y,x,i) into patches(y,x,1,i)
% Also, adds a black border aroudn them

if iscell(arg)
  h= montageFilenames(arg);
else
  nr = size(arg,1); nc = size(arg,2); Npatches = size(arg,3);
  patchesColor = reshape(arg, [nr nc 1 Npatches]);
  patchesColor = patchesColor ./ max(patchesColor(:));

  if 1
    %put a black border around them for display purposes
    border = 5;
    bgColor = ones(1,1,class(patchesColor));
    patchesColorBig = bgColor*ones(nr+2*border, nc+2*border, 1, Npatches, class(patchesColor));
    %patchesColorBig = zeros(nr+2*border, nc+2*border, 1, Npatches, class(patchesColor));
    patchesColorBig(border+1:end-border, border+1:end-border, :, :) = patchesColor;
  else
    patchesColorBig = patchesColor;
  end
  montage(patchesColorBig)

end

%%%%%%%%%%%%%

function h = montageFilenames(filenames)

%[nRows, nCols, nBands, nFrames] = size(a);

% Estimate nMontageColumns and nMontageRows given the desired ratio of
% Columns to Rows to be one (square montage).
aspectRatio = 1;
nMontageCols = sqrt(aspectRatio * nRows * nFrames / nCols);

% Make sure montage rows and columns are integers. The order in the adjustment
% matters because the montage image is created horizontally across columns.
nMontageCols = ceil(nMontageCols);
nMontageRows = ceil(nFrames / nMontageCols);

% Create the montage image.
b = a(1,1); % to inherit type
b(1,1) = 0; % from a
b = repmat(b, [nMontageRows*nRows, nMontageCols*nCols, nBands, 1]);

rows = 1 : nRows;
cols = 1 : nCols;

for i = 0:nMontageRows-1
  for j = 0:nMontageCols-1,
    k = j + i * nMontageCols + 1;
    if k <= nFrames
      b(rows + i * nRows, cols + j * nCols, :) = a(:,:,:,k);
    else
      break;
    end
  end
end

if isempty(cm)
  hh = imshow(b);
else
  hh = imshow(b,cm);
end

if nargout > 0
    h = hh;
end

%--------------------------------------------------------------
%Parse Inputs Function

function [I,map] = parse_inputs(varargin)

% initialize variables
map = [];

iptchecknargin(1,2,nargin,mfilename);
iptcheckinput(varargin{1},{'uint8' 'double' 'uint16' 'logical' 'single' ...
                    'int16'},{},mfilename, 'I, BW, or RGB',1);
I = varargin{1};

if nargin==2
  if isa(I,'int16')
    eid = sprintf('Images:%s:invalidIndexedImage',mfilename);
    msg1 = 'An indexed image can be uint8, uint16, double, single, or ';
    msg2 = 'logical.';
    error(eid,'%s %s',msg1, msg2);
  end
  map = varargin{2};
  iptcheckinput(map,{'double'},{},mfilename,'MAP',1);
  if ((size(map,1) == 1) && (prod(map) == numel(I)))
    % MONTAGE(D,[M N P]) OBSOLETE
    eid = sprintf('Images:%s:obsoleteSyntax',mfilename);
    msg1 = 'MONTAGE(D,[M N P]) is an obsolete syntax.';
    msg2 = 'Use multidimensional arrays to represent multiframe images.';
    error(eid,'%s\n%s',msg1,msg2);
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [X_hat]=SD_detector(y,H,nT)
    % Input parameters
    %     y : received signal, nRx1
    %     H : Channel matrix, nRxnT
    %    nT : number of Tx antennas
    % Output parameter
    %    X_hat : estimated signal, nTx1

    global x_list;         % candidate symbols in real constellations
    global x_now;          % temporary x_vector elements
    global x_hat;          % inv(H)*y
    global x_sliced;       % sliced x_hat
    global x_pre;          % x vectors obtained in the previous stage
    global real_constellation; % real constellation
    global R;               % R in the QR decomposition
    global radius_squared; % radius^2
    global x_metric;       % ML metrics of previous stage candidates
    global len;             % nT*2
    QAM_table2 = [-3-3j, -3-j, -3+3j, -3+j, -1-3j, -1-j, -1+3j, -1+j,3-3j, ...
               3-j, 3+3j, 3+j, 1-3j, 1-j, 1+3j, 1+j]/sqrt(10); % 16-QAM
    real_constellation = [-3 -1 1 3]/sqrt(10);
    y =[real(y); imag(y)];     % y : complex vector  -> real vector
    H =[real(H)  -(imag(H)) ; imag(H)   real(H)];
    % H : complex vector  -> real vector
    len = nT*2; % complex -> real
    x_list = zeros(len,4); % 4 : real constellation length, 16-QAM
    x_now = zeros(len,1); x_hat = zeros(len,1); x_pre = zeros(len,1); x_metric = 0;
    [Q,R] = qr(H);     % nR x nT QR decomposition
    x_hat = inv(H)*y;                % zero forcing equalization
    x_sliced = QAM16_real_slicer(x_hat,len)';  % slicing
    radius_squared  = norm(R*(x_sliced-x_hat))^2;  % Radious^2
    transition = 1;
    % meaning of transition
    % 0 : radius*2, 1~len : stage number
    % len+1 : compare two vectors in terms of norm values
    % len+2 : finish
    flag = 1;
    % transition tracing 0 : stage index increases by +1
    %1 : stage index decreases by -1
    %2 : 1->len+2 or len+1->1
    while (transition<len+2)
       if transition==0    % radius_squared*2
         [flag,transition,radius_squared,x_list]= radius_control(radius_squared,transition);
        elseif transition <= len
         [flag,transition] = stage_processing(flag,transition);
        elseif transition == len+1 %
         [flag,transition] = compare_vector_norm(transition);
       end
    end
    ML = x_pre;
    for i=1:len/2
        X_hat(i) = ML(i)+j*ML(i+len/2);
    end
end

function [flag,transition] = stage_processing(flag,transition)
    % Input parameters
    %    flag : previous stage index
    %       flag = 0 : stage index decreased -> x_now empty -> new x_now
    %       flag = 1 : stage index decreased -> new x_now
    %       flag = 2 : previous stage index =len+1 ->  If R>R'? start from the first stage
    %     transition : stage number
    % Output parameters
    %     flag : stage number is calculated from flag
    %     transition : next stage number, 0 : R*2, 1: next stage, len+2: finish
    global x_list x_metric x_now x_hat real_constellation R radius_squared x_sliced;

    global x_list;
    global x_metric;
    global x_now;
    global x_hat;
    global real_constellation;
    global R;
    global radius_squared;
    global x_sliced;
    stage_index = length(R(1,:))-(transition-1);
    if flag == 2  % previous stage=len+1 : recalculate radius R'
      radius_squared  = norm(R*(x_sliced-x_hat))^2;
    end
    if flag ~= 0 % previous stage=len+1 or 0
    -> upper and lower bound calculation, x_list(stage_index,:)
        [bound_lower bound_upper] = bound(transition);
        for i =1:4    % search for a candidate in x_now(stage_index),
           % 4=size(real_constellation), 16-QAM assumed
           if bound_lower <= real_constellation(i) && real_constellation(i) <= bound_upper
             list_len = list_length(x_list(stage_index,:));
             x_list(stage_index,list_len+1) = real_constellation(i);
           end
        end
    end
    list_len = list_length(x_list(stage_index,:));
    if list_len == 0     % no candidate in x_now
      if x_metric == 0 || transition ~= 1
        % transition >=2 ? if no candidate ? decrease stage index
        flag = 0;
        transition = transition-1;
       elseif x_metric ~= 0 && transition == 1
        % above two conditions are met? ML solution found
        transition = length(R(1,:))+2;  % finish stage
      end
    else              % candidate exist in x_now ? increase stage index
      flag = 1;
      transition = transition+1;
      x_now(stage_index) = x_list(stage_index,1);
      x_list(stage_index,:) = [x_list(stage_index,[2:4]) 0];
end

function [bound_lower bound_upper]=bound(transition)
    % Input parameters
    %     R : [Q R] = qr(H)
    %     radius_squared : R^2
    %     transition : stage number
    %     x_hat : inv(H)*y
    %     x_now : slicing x_hat
    % Output parameters
    %     bound_lower : bound lower
    %     bound_upper : bound upper

    global R  radius_squared  x_now  x_hat;
    len = length(x_hat);
    temp_sqrt = radius_squared;
    temp_k=0;
    for i=1:1:transition-1
       temp_abs=0;
       for k=1:1:i
          index_1 = len-(i-1);
          index_2 = index_1+ (k-1);
          temp_k = R(index_1,index_2)*(x_now(index_2)-x_hat(index_2));
          temp_abs=temp_abs+temp_k;
       end
       temp_sqrt = temp_sqrt - abs(temp_abs)^2;
    end
    temp_sqrt = sqrt(temp_sqrt);
    temp_no_sqrt = 0;
    index_1 = len-(transition-1);
    index_2 = index_1;
    for i=1:1:transition-1
       index_2 = index_2+1;
       temp_i = R(index_1,index_2)*(x_now(index_2)-x_hat(index_2));
       temp_no_sqrt = temp_no_sqrt - temp_i;
    end
    temp_lower = -temp_sqrt + temp_no_sqrt;
    temp_upper = temp_sqrt + temp_no_sqrt;
    index = len-(transition-1);
    bound_lower = temp_lower/R(index,index) + x_hat(index);
    bound_upper = temp_upper/R(index,index) + x_hat(index);
    bound_upper = fix(bound_upper*sqrt(10))/sqrt(10);
    bound_lower = ceil(bound_lower*sqrt(10))/sqrt(10);
end

function [len]=list_length(list)
    % Input parameter
    %     list : vector type
    % Output parameter
    %     len : index number

    len = 0;
    for i=1:4
       if list(i)==0,  break;  else len = len+1;  end
    end
end

function [flag,transition,radius_squared,x_list] =radius_control(radius_squared,transition)
    % Input parameters
    %     radius_squared : current radius
    %     transition : current stage number
    % Output parameters
    %     radius_squared : doubled radius
    %     transition : next stage number
    %     flag : next stage number is calculated from flag
    global len;
    radius_squared = radius_squared*2;
    transition = transition+1;
    flag = 1;
    x_list(len,:)=zeros(1,4);
end

function [check]=vector_comparison(vector_1,vector_2)
    % check if the two vectors are the same
    % Input parameters
    %   pre_x : vector 1
    %   now_x : vector 2
    % Output parameters
    %   check : 1-> same vectors, 0-> different vectors
    check = 0;
    len1 = length(vector_1);  len2 = length(vector_2);
    if len1 ~= len2
      error('vector size is different');
    end
    for column_num = 1:len1
       if vector_1(column_num,1) == vector_2(column_num,1)
         check = check + 1;
       end
    end
    if check == len1,  check = 1;
     else    check = 0;
    end
end

function [flag,transition]=compare_vector_norm(transition)
    % stage index increased(flag = 1) : recalculate x_list(index,:)
    % stage index decreased(flag = 0) : in the previous stage, no candidate x_now in x_list
    % Input parameters
    %     flag : previous stage
    %     transition : stage number
    % Output parameters
    %    flag : next stage number is calculated from flag
    %    transition : next stage number
    global x_list x_pre x_metric x_now x_hat R radius_squared x_sliced len;
    vector_identity = vector_comparison(x_pre,x_now);
    % check if the new candidate is among the ones we found before
    if vector_identity == 1
      % if 1 ? ML solution found
      len_total = 0;
      for i=1:len  % if the vector is unique ? len_total = 0
         len_total = len_total + list_length(x_list(i,:));
      end
      if len_total == 0      % ML solution vector found
        transition = len+2; % finish
        flag = 1;
       else                      % more than one candidates
        transition = transition-1;  % go back to the previous stage
        flag =0;
      end
     else  % if 0 ? new candidate vector is different from the previous candidate vector and norm is smaller ? restart
      x_sliced_temp = x_now;
      metric_temp  = norm(R*(x_sliced_temp-x_hat))^2;
      if metric_temp <=  radius_squared
        % new candidate vector has smaller metric ? restart
        x_pre = x_now;  x_metric = metric_temp;
        x_sliced = x_now;  transition = 1;       % restart
        flag = 2;  x_list=zeros(len,4); % initialization
        x_now=zeros(len,1);  % initialization
       else % new candidate vector has a larger ML metric
        transition = transition-1;  % go back to the previous stage
        flag =0;
      end
    end
end

author: LyricYang


text:
function [ x, seed ] = sech_sample ( a, b, seed )

%*****************************************************************************80
%
%% SECH_SAMPLE samples the Hyperbolic Secant PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ cdf, seed ] = r8_uniform_01 ( seed );

  x = a + b * log ( tan ( 0.5 * pi * cdf ) );

  return
end

author: johannesgerer


function video_feat = process_box_feat_prune(box_feat, video_info, prune, lp, pyramid)

nVideo = size(video_info, 1);
video_feat = zeros(size(box_feat, 1) * (pyramid), nVideo);
for n = 1:nVideo
    feature_set = box_feat(:, video_info(n, 1):video_info(n, 2));
    feature_set = double(feature_set);

    if prune ~= 0
        %% prune video track ****************
        prune_feat = feature_set;

        sum_val = sqrt(sum(prune_feat.^2));
        for qqq = 1:size(prune_feat, 1)
            prune_feat(qqq, :) = prune_feat(qqq, :)./sum_val;
        end

        dist_prune = pdist2(prune_feat', prune_feat');
        avg_sum = sum(dist_prune, 2);
        avg = mean(avg_sum);

        feature_set = feature_set(:, avg_sum(:, 1) <= avg*prune);

    end
    %% pyramid pooling
    if pyramid ~= 1

        video_len = size(feature_set, 2);

        video_feat(1: 1024, n) = sum(feature_set.^lp, 2).^(1/lp) * (1/size(feature_set, 2)).^(1/lp);

        py_2_1 = round(video_len/4);
        py_2_2 = 2 *round(video_len/4);
        py_2_3 = 3 * round(video_len/4);

        %% pyramid 2
        if video_len <= 20

            video_feat(1025 : 2048, n) = sum(feature_set.^lp, 2).^(1/lp) * (1/size(feature_set, 2)).^(1/lp);
            video_feat(2049 : 3072, n) = sum(feature_set.^lp, 2).^(1/lp) * (1/size(feature_set, 2)).^(1/lp);
            video_feat(3073 : 4096, n) = sum(feature_set.^lp, 2).^(1/lp) * (1/size(feature_set, 2)).^(1/lp);
            video_feat(4097 : 5120, n) = sum(feature_set.^lp, 2).^(1/lp) * (1/size(feature_set, 2)).^(1/lp);

        else
            feature_set_2_1 = feature_set(:, 1: py_2_1);
            feature_set_2_2 = feature_set(:, py_2_1  : py_2_2);
            feature_set_2_3 = feature_set(:, py_2_2 : py_2_3);
            feature_set_2_4 = feature_set(:, py_2_3 : end);
            video_feat(1025 : 2048, n) = sum(feature_set_2_1.^lp, 2).^(1/lp) * (1/size(feature_set_2_1, 2)).^(1/lp);
            video_feat(2049 : 3072, n) = sum(feature_set_2_2.^lp, 2).^(1/lp) * (1/size(feature_set_2_2, 2)).^(1/lp);
            video_feat(3073 : 4096, n) = sum(feature_set_2_3.^lp, 2).^(1/lp) * (1/size(feature_set_2_3, 2)).^(1/lp);
            video_feat(4097 : 5120, n) = sum(feature_set_2_4.^lp, 2).^(1/lp) * (1/size(feature_set_2_4, 2)).^(1/lp);
        end

    else
        %% lp-norm pooling
        feature_set = (feature_set).^lp;
        video_feat(:, n) = sum(feature_set, 2).^(1/lp) * (1/size(feature_set, 2)).^(1/lp);

    end

    %% **********************
    % video_feat(:, n) = max(feature_set, [], 2); % max pooling
    %     video_feat(:, n) = mean(feature_set, 2); % avg pooling
end

%%% normalize train and test features
sum_val = sqrt(sum(video_feat.^2));
for n = 1:size(video_feat, 1)
    video_feat(n, :) = video_feat(n, :)./sum_val;
end


text:
function DEM_demo_Cornsweet
% The Cornsweet effect: This demo illustrates the inference underlying the
% Cornsweet effect or illusion. It exploits formal priors on the spatial
% contiguity of the illuminant and reflectance; where the illuminant does not
% have edges, but the reflectance can. This is implemented using a
% discrete cosine set (DCT) as the spatial basis for the illuminant and a
% (Haar) Discrete Wavelet transform (DWT) for the reflectance. Appropriate
% shrinkage priors on the (implicit) transform coefficients ensure that the
% explanation for visual input (reflectance times illuminant) assigns edges
% to the reflectance; thereby producing the Cornsweet effect.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_demo_Cornsweet.m 4851 2012-08-20 15:03:48Z karl $

% Illustrate the Cornsweet effect
%==========================================================================
spm_figure('GetWin','Figure 1');
colormap((1:255)'*[1 1 1]/255)

% basic profile
%--------------------------------------------------------------------------
nx    = 64;
YX    = [zeros(1,nx*3/8) 1:(nx/8) -(nx/8):-1 zeros(1,nx*3/8)]*8/nx;
YT    = ones(32,1);

% present with different contrasts (C)
%--------------------------------------------------------------------------
C     = [1 4 32 64];
nc    = length(C);
for i = 1:length(C)

    subplot(nc,2,2*i - 1)
    image(YT*YX*C(i) + 128)
    title('Cornsweet effect','FontSize',16)
    xlabel('eccentricity ','FontSize',12)
    axis square

    subplot(nc,2,2*i)
    plot(YX*C(i))
    title('Reflectance ','FontSize',16)
    xlabel('eccentricity ','FontSize',12)
    axis square, axis([1 nx -64 64]);

end

% Simulations
%==========================================================================

% Basis functions of illuminant and reflectance
%--------------------------------------------------------------------------
nx      = 32;                               % number of pixels
P.R     = spm_dwtmtx(nx,1,1);               % DWT for hidden causes
P.I     = spm_dctmtx(nx,3);                 % DCT for hidden causes
P.I     = P.I/diag(max(P.I));
nr      = size(P.R,2);
ni      = size(P.I,2);
W       = log2(nx./sum(P.R > 0))*4;         % scale of reflectance
W(1)    = 16;

% initial hidden states and causes
%--------------------------------------------------------------------------
x       = sparse(nr,1);
v       = sparse(ni + nr,1);
v(1)    = -1;
ii      = (1:ni);
ir      = (1:nr) + ni;

% level 1
%--------------------------------------------------------------------------
M(1).f  = inline('(v(4:end) - x)/16','x','v','P');
M(1).g  = inline('exp(P.R*x + P.I*v(1:3))','x','v','P');
M(1).pE = P;                                % The prior expectation
M(1).x  = x;                                % The prior expectation
M(1).V  = exp(6);                           % error precision (data)
M(1).W  = exp(10);                          % error precision (motion)
M(1).xP = diag(exp(W));                     % error precision (motion)

% level 2
%--------------------------------------------------------------------------
M(2).v  = v;                                % hidden causes
M(2).V  = exp(0);                           % error precision (cause)

% Create stimulus
%==========================================================================

% Create stimulus (with spatial and temporal envelopes YX and YT)
%--------------------------------------------------------------------------
M(1).E.n = 1;
N        = 64;                             % length of sequence
t        = ([1:N] - N/4)*8;                % time (ms)
YX       = [zeros(1,nx*3/8) 1:(nx/8) -(nx/8):-1 zeros(1,nx*3/8)]*8/nx;
YX       = YX'/4 + 1;
YX       = exp(P.R*pinv(P.R)*log(YX));
YT       = exp((tanh(([1:N] - N/3)/8) - 1)*4);
YT       = exp(-([1:N] - N/2).^2/(2*(N/8)^2));
Y        = YX*YT;

% invert
%==========================================================================
DEM.M    = M;
DEM.Y    = Y;
DEM      = spm_DEM(DEM);

% render true and perceived stimuli
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 2'); clf

subplot(2,2,1)
PI  = P.I*DEM.qU.v{2}(ii,:);
imagesc(exp(PI))
title('Perceived illunimant','FontSize',16)
xlabel('time (bins)','FontSize',12)
ylabel('eccentricity ','FontSize',12)
axis square

subplot(2,2,2)
PR = P.R*DEM.qU.x{1};
imagesc(exp(PR))
title('Perceived reflectance','FontSize',16)
xlabel('time (bins)','FontSize',12)
ylabel('eccentricity ','FontSize',12)
axis square

subplot(2,2,3)
imagesc(DEM.qU.v{1})
title('Predicted stimulus','FontSize',16)
xlabel('time (bins)','FontSize',12)
ylabel('eccentricity ','FontSize',12)
axis square

% first order effect: (R1) - Cornsweet
%--------------------------------------------------------------------------
C     = sparse(1,nx);
i     = nx/2 - nx/4 - 1;
C(i)  =  1;
i     = nx/2 + nx/4 + 2;
C(i)  = -1;
R1    = C*P.R;

% second order effect: (R2) - Mach bands
%--------------------------------------------------------------------------
i     = nx/2 - nx/8 - 1;
C(i)  = -1;
i     = nx/2 + nx/8 + 2;
C(i)  =  1;
R2    = C*P.R;

D1    = R1*DEM.qU.x{1};
D2    = R2*DEM.qU.x{1};
for i = 1:N
    V1(i) = R1*DEM.qU.S{i}*R1';
    V2(i) = R2*DEM.qU.S{i}*R2';
end

subplot(2,2,4)
spm_plot_ci(D1,V1,t), hold on
spm_plot_ci(D2,V2,t), hold off
title('Perceived reflectance','FontSize',16)
xlabel('time (ms)','FontSize',12)
ylabel('eccentricity ','FontSize',12)
axis square
spm_axis tight
drawnow

% Cycle over different levels of visual precision (cf contrast)
%==========================================================================
C   = -2:2:12;
for i = 1:length(C)

    % Change precision (contrast) and invert
    %----------------------------------------------------------------------
    SIM{i}        = DEM;
    SIM{i}.M(1).V = exp(C(i));
    SIM{i}        = spm_DEM(SIM{i});

    % Record conditional peripheral reflectance difference
    %----------------------------------------------------------------------
    d1(i) = R1*SIM{i}.qU.x{1}(:,N/2);
    v1(i) = R1*SIM{i}.qU.S{N/2}*R1';
    d2(i) = R2*SIM{i}.qU.x{1}(:,N/2);
    v2(i) = R2*SIM{i}.qU.S{N/2}*R2';

end

% show Cornsweet effect as a function of precision (contrast)
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 3'); clf

subplot(3,1,1)
spm_plot_ci(d1,v1,C),                hold on
plot(C,d1,'ob',C,C*0,'LineWidth',2), hold off
title('Conditional difference (Cornsweet)','FontSize',16)
xlabel('log-precision (contrast)','FontSize',12)
ylabel('reflectance difference','FontSize',12)
spm_axis tight square

subplot(3,1,2)
spm_plot_ci(d2,v2,C),                hold on
plot(C,d2,'or',C,C*0,'LineWidth',2), hold off
title('Conditional difference (Mach Band)','FontSize',16)
xlabel('log-precision (contrast)','FontSize',12)
ylabel('reflectance difference','FontSize',12)
spm_axis tight square

% and associated percepts
%--------------------------------------------------------------------------
colormap([1:255]'*[1 1 1]/255)
j     = [1 5 7];
nj    = length(j);
for i = 1:nj

    % predictions
    %----------------------------------------------------------------------
    PR    = P.R*SIM{j(i)}.qU.x{1};
    PR    = PR*256 + 128;

    subplot(3,nj,i + 2*nj)
    image(PR)
    ylabel('eccentricity','FontSize',12)
    axis square

end

% and associated prediction errors (cf ERPs)
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 4'); clf
as    = 0;
ax    = 0;
av    = 0;
for i = 1:nj

    % predictions and errors
    %----------------------------------------------------------------------
    Ev = SIM{j(i)}.qU.z{2}(ir,:);
    Es = SIM{j(i)}.qU.z{1};
    av = max(max(abs(Ev(:))),av);
    as = max(max(abs(Es(:))),as);

    subplot(3,nj,i)
    plot(t,Es,'r')
    title('Error (sensory)','FontSize',14)
    xlabel('time (ms)','FontSize',12)
    ylabel('error (state)','FontSize',12)

    subplot(3,nj,i + nj)
    plot(t,Ev,'r')
    title('Error (hidden states)','FontSize',14)
    xlabel('time (ms)','FontSize',12)
    ylabel('error (state)','FontSize',12)

end

% adjust axes
%--------------------------------------------------------------------------
for i = 1:nj

    subplot(3,nj,i)
    axis square, axis([t(1) t(end) -as*1.2 as*1.2])

    subplot(3,nj,i + nj)
    axis square, axis([t(1) t(end) -av*1.2 av*1.2])

end

% redraw inference for effect
%--------------------------------------------------------------------------
spm_figure('GetWin','DEM');
spm_DEM_qU(SIM{7}.qU)

% Fitting behavioural data
%==========================================================================
try, load cornsweet_data, catch, return, end

% save simulation results
%--------------------------------------------------------------------------
sim.Con    = C;               % simulated contrast
sim.Ecs    = d1;              % conditional expectation of Cornsweet effect
sim.Emb    = d2;              % conditional expectation of mach band effect
sim.Vcs    = v1;              % conditional dispersion of Cornsweet effect
sim.Vmb    = v2;              % conditional dispersion of mach band effect

% add emprical contrasts that were used
%--------------------------------------------------------------------------
sim.Con_cs = cornsweet.cornsweet;
sim.Con_mb = mach.machContrast;
G.sim      = sim;     % place in model

% empirical responses
%--------------------------------------------------------------------------
Y          = struct;
Y.y{1}     = cornsweet.stepMatch(:);
Y.y{2}     = mach.pSeeMach(:);
Y.Q        = spm_Ce([length(Y.y{1}) length(Y.y{2})]); % error precisions

% model parameters
%--------------------------------------------------------------------------
B.thresh = -2;                % log-contrast threshold for Mach band
B.sig    = 1;                 % parameter of contrast function
B.off    = 2;                 % parameter of contrast function
B.con    = 1/32;              % scaling of reported contrast
B.sen    = 4;                 % parameter of response probability
nb       = length(spm_vec(B));

% setup model of empirical responses
%--------------------------------------------------------------------------
G.IS = 'spm_cornsweet';       % function name f(P,M,U) - generative model
G.pE = B;                     % prior expectation of model parameters
G.pC = speye(nb,nb)/8;        % prior covariance  of model parameters
G.hE = [16; 8];               % prior expectation of log-precisions
G.hC = exp(-4);               % prior covariance of log-precisions

% invert model of empirical responses and plot
%--------------------------------------------------------------------------
Ep   = spm_nlsi_GN(G,[],Y);
spm_cornsweet(Ep,G,Y);

return

% plot behavioural results alone.
%--------------------------------------------------------------------------
subplot(2,2,1)
SE   = cornsweet.error;
errorbar(cornsweet.cornsweet,cornsweet.stepMatch,SE/2)
xlabel('empirical contrast','Fontsize',12)
ylabel('reported contrast','Fontsize',12)
title('Cornsweet','Fontsize',16)
set(gca,'XLim',[0 .15]);
axis square

subplot(2,2,2)
SE   = mach.error;
errorbar(mach.machContrast,mach.pSeeMach,SE/2)
xlabel('empirical contrast','Fontsize',12)
ylabel('report probability','Fontsize',12)
title('Mach Bands','Fontsize',16)
set(gca,'XLim',[0 .15]);
axis square

% Ancillary code for producing figures:
%==========================================================================
spm_figure('GetWin','paper'); clf

n = 128;
a = 4;
i = tanh(([1:n] - n/2)/(n/8))/a;
r = kron([1 -1],ones(1,n/2))/a;
s = exp(i + r);
c = zeros(1,n);

subplot(5,1,1), plot(i),     spm_axis tight, title('illuminant')
subplot(5,1,2), plot(r),     spm_axis tight, title('reflectant')
subplot(5,1,3), plot(s),     spm_axis tight, title('stimulus')
subplot(5,1,4), plot(i + r), spm_axis tight, title('reflectant')
subplot(5,1,5), plot(c),     spm_axis tight, title('illuminant')

% Illustrate the generative model (for 1D)
%==========================================================================
spm_figure('GetWin','paper'); clf

DEM      = spm_DEM_generate(M,64,P,{16 0},{16});

subplot(3,2,1)
PI  = P.I*DEM.pU.v{2}(ii,:);
imagesc(exp(PI))
title('Illunimant','FontSize',16)
xlabel('time (bins)','FontSize',12)
ylabel('location ','FontSize',12)
axis square

subplot(3,2,3)
PR = P.R*DEM.pU.x{1};
imagesc(exp(PR))
title('Reflectant','FontSize',16)
xlabel('time (bins)','FontSize',12)
ylabel('location ','FontSize',12)
axis square

subplot(3,2,5)
imagesc(DEM.pU.v{1})
title('Stimulus','FontSize',16)
xlabel('time (bins)','FontSize',12)
ylabel('location ','FontSize',12)
axis square

subplot(3,2,2), imagesc(P.I), axis image off
subplot(3,2,4), imagesc(P.R), axis image off

% Illustrate the generative model (for 2D)
%==========================================================================
nx    = 128;
P.R   = spm_dwtmtx(nx,2);                 % DWT for hidden causes
P.I   = spm_dctmtx(nx,3);                 % DCT for hidden causes
I     = 0;
R     = 0;
for i = 1:size(P.I,2)
    for j = 1:size(P.I,2)
        I = I + (P.I(:,i)*P.I(:,j)')*randn;
    end
end
for i = 1:size(P.R,2)
    for j = 1:size(P.R,2)
        v = sum(P.R(:,i) > 0)*sum(P.R(:,j) > 0);
        R = R + (P.R(:,i)*P.R(:,j)')*randn*sqrt(v);
    end
end

R   = R/std(R(:))/2;
I   = I/std(I(:))/2;

% Plot
%--------------------------------------------------------------------------
spm_figure('GetWin','Paper');

subplot(2,2,1)
imagesc(exp(I))
title('Illuminant ','FontSize',16)
axis square

subplot(2,2,2)
imagesc(exp(R))
title('Reflectance ','FontSize',16)
axis square

subplot(2,1,2)
imagesc(exp(I + R))
title('Image','FontSize',16)
axis square, drawnow

author: spm


text:
% SUMMARY:  convert LogGamma, LogKsi to Gamma, Ksi. Subtract max to avoid
%           underflow
% AUTHOR:   QIUQIANG KONG
% Created:  17-11-2015
% Modified: 25-11-2015 Add annotation
% -----------------------------------------------------------
% input:
%   LogGamma   cell{ ln p(zn|X) }
%   LogKsi     cell{ ln p(zn,zn-1|X) }
% output:
%   Gamma      cell{ gamma }
%   Ksi        cell{ ksi }
% ===========================================================
function [Gamma, Ksi] = UniformLogGammaKsi(LogGamma, LogKsi)
obj_num = length(LogGamma);
Q = size(LogGamma{1}, 2);
for q = 1:Q
    max_gamma_ary = zeros(1, obj_num);
    max_ksi_ary = zeros(1, obj_num);
    for r = 1:obj_num
        [Nr, Q] = size(LogGamma{r});
        max_gamma_ary(r) = max(LogGamma{r}(:,q));
        max_ksi_ary(r) = max(reshape(LogKsi{r}(:,q,:), (Nr-1)*Q, 1));
    end
    max_gamma = max(max_gamma_ary);
    max_ksi = max(max_ksi_ary);

    for r = 1:obj_num
        LogGamma{r}(:,q) = LogGamma{r}(:,q) - max_gamma;
        LogKsi{r}(:,q,:) = LogKsi{r}(:,q,:) - max_ksi;
    end
end

for r = 1:obj_num
    Gamma{r} = exp(LogGamma{r});
    Ksi{r} = exp(LogKsi{r});
end
end

author: qiuqiangkong

repo: matlab-hmm


text:
function y = ivmOut(model, x);

% IVMOUT Evaluate the output of an IVM model.
% FORMAT
% DESC evaluates the output of a given IVM model.
% ARG model : the model for which the output is being evaluated.
% ARG x : the input position for which the output is required.
% RETURN y : the output of the GP model. The function checks if
% there is a noise model associated with the GP, if there is, it is
% used, otherwise the mean of gpPosteriorMeanVar is returned.
%
% SEEALSO : ivmCreate, ivmPosteriorMeanVar, noiseOut
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2005

% IVM

if nargin < 2
  % This implies evaluate for the training data.
  mu = model.mu;
  varsigma = model.varSigma;
else
  [mu, varsigma] = ivmPosteriorMeanVar(model, x);
end

y = noiseOut(model.noise, mu, varsigma);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_example_megrealign

% MEM 6gb
% WALLTIME 00:10:00

%
%% Interpolating data from the CTF151 to the CTF275 sensor array using megrealign
%
% the two example files with the 151 and 275 channel CTF specifications are available from
% ftp://ftp.fieldtriptoolbox.org/pub/fieldtrip/example/megrealign/

grad151 = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/data/ftp/example/megrealign/ctf151.mat'));
grad275 = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/data/ftp/example/megrealign/ctf275.mat'));

vol = [];
vol.r = 12;
vol.o = [0 0 4];

% in this example the dipole is carefully positioned on the same depth as
% where the dipole layer for the interpolation will be located. The center
% of the head is at [0 0 4], the radius is 12 cm, and the dipole layer that is
% specified in ft_megrealign is 2.5 cm shifted inward from the head surface
cfg = [];
cfg.dip.pos = [0 0 13.5];  % 4 + 12 - 2.5
cfg.dip.frequency = 1;
cfg.headmodel = vol;
cfg.grad = grad151;
data151 = ft_dipolesimulation(cfg);
cfg.grad = grad275;
data275 = ft_dipolesimulation(cfg);

cfg = [];
avg151 = ft_timelockanalysis(cfg, data151);
avg275 = ft_timelockanalysis(cfg, data275);

% here I'll not only realign the 151 channel MEG data to the 275 channel and
% vice versa, but also apply the ft_megrealign function from 151 to 151. In the
% realignment the signal will be slightly distorted as a comparison of avg151_151
% to avg151 would show.
cfg = [];
cfg.inwardshift = 3;
cfg.headmodel = vol;
cfg.template{1} = grad151; avg151_151 = ft_timelockanalysis([], ft_megrealign(cfg, avg151));
cfg.template{1} = grad275; avg151_275 = ft_timelockanalysis([], ft_megrealign(cfg, avg151));
cfg.template{1} = grad151; avg275_151 = ft_timelockanalysis([], ft_megrealign(cfg, avg275));
cfg.template{1} = grad275; avg275_275 = ft_timelockanalysis([], ft_megrealign(cfg, avg275));

% plot the realigned datasets
cfg = [];
figure; ft_multiplotER(cfg, avg151_151);
title('Interpolated from 151 to 151');
orient portrait; %print -dpng -r300 avg151_151.png

%
figure; ft_multiplotER(cfg, avg151_275);
title('Interpolated from 151 to 151');
orient portrait; %print -dpng -r300 avg151_275.png

%
figure; ft_multiplotER(cfg, avg275_151);
title('Interpolated from 151 to 151');
orient portrait; %print -dpng -r300 avg275_151.png

%
figure; ft_multiplotER(cfg, avg275_275);
title('Interpolated from 151 to 151');
orient portrait; %print -dpng -r300 avg275_275.png

%
% now plot them together
cfg = [];
figure; ft_multiplotER(cfg, avg151, avg151_151, avg275_151);
orient portrait; %print -dpng -r300 compare151.png

%
figure; ft_multiplotER(cfg, avg275, avg275_275, avg151_275);
orient portrait; %print -dpng -r300 compare275.png

author: fieldtrip


text:
function b = r8pbu_mxv ( n, mu, a, x )

%*****************************************************************************80
%
%% R8PBU_MXV multiplies a R8PBU matrix times a vector.
%
%  Discussion:
%
%    The R8PBU storage format is for a symmetric positive definite band matrix.
%
%    To save storage, only the diagonal and upper triangle of A is stored,
%    in a compact diagonal format that preserves columns.
%
%    The diagonal is stored in row MU+1 of the array.
%    The first superdiagonal in row MU, columns 2 through N.
%    The second superdiagonal in row MU-1, columns 3 through N.
%    The MU-th superdiagonal in row 1, columns MU+1 through N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer MU, the number of superdiagonals in the matrix.
%    MU must be at least 0 and no more than N-1.
%
%    Input, real A(MU+1,N), the R8PBU matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(N), the result vector A * x.
%

%
%  Multiply X by the diagonal of the matrix.
%
  b(1:n) = a(mu+1,1:n) .* x(1:n);
%
%  Multiply X by the superdiagonals of the matrix.
%
  for i = mu : -1 : 1
    for j = mu+2-i : n
      ieqn = i + j - mu - 1;
      b(ieqn) = b(ieqn) + a(i,j) * x(j);
      b(j) = b(j) + a(i,j) * x(ieqn);
    end
  end

  return
end


text:
function [segmentsM] = mask2scan(maskM, optS, sliceNum)
%function [segmentsM] = mask2scan(maskM, optS, sliceNum)
%segmentsM holds the mask in 'scan segment' format:
%each row of this 2-D matrix has the elements:
%yValue, xStart, xStop, delta_x (all in AAPM coordinates).
%xStart and xStop are the points at which the filled-in segment
%begins and ends; yValue is the y-value of that row.  delta_x is the
%interval between x values.
%
%Get scan format information:
%
%Starting points for filled segments will be those after diff = 1; stopping points will
%those where diff = -1.  This also works for jutting points (one point
%in the segment).
%
%J.O.Deasy, deasy@radonc.wustl.edu.
%LM:  18 Mar 02, JOD;
%Major bug fix: 28 dec 02: Scan row was shifted up one row;
%scan column was shifted to the right one column.  Fixed prior to release
%of version 2 beta, JOD.

imageSizeV = optS.ROIImageSize;

%First convert to pixel-based coords:
xOffset = optS.xCTOffset;
yOffset = optS.yCTOffset;

zerosV = zeros(size(maskM,1),1);  %padding vector

delta_x = optS.ROIxVoxelWidth;
delta_y = optS.ROIyVoxelWidth;

mask2M = [zerosV, maskM, zerosV];

diffM = diff(mask2M');  %note: mask is rotated here!
startM = [diffM == 1];
startM = startM(:,2:end-1);
stopM  = [diffM == -1];
stopM = stopM(:,2:end-1);

[i1V, j1V] = find(startM);
[xStartV, yStartV] = mtoaapm(j1V + 1, i1V, size(maskM));
xStartV = xStartV * delta_x;
yStartV = yStartV * delta_y;

xStartV = xStartV + optS.xCTOffset;
yStartV = yStartV + optS.yCTOffset;

[i2V, j2V] = find(stopM);
[xStopV, yStopV] = mtoaapm(j2V + 1, i2V - 1 , size(maskM));
xStopV = xStopV * delta_x;
yStopV = yStopV * delta_y;

xStopV = xStopV + optS.xCTOffset;
yStopV = yStopV + optS.yCTOffset;

if any(yStartV ~= yStopV)
    error('Ooops! Problem in converting a mask to scan segments.  Check options.')
end

%Store: row y, start x, stop x, delta x, CT mask row, CT start col, CT end col.
z1V = ones(length(j1V),1) * sliceNum;
segmentsM = [yStartV(:), xStartV(:), xStopV(:), ones(length(xStopV),1) * delta_x, z1V, j1V(:) + 1, i1V(:), i2V(:) - 1];

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ iarray, more ] = vec_next ( n, ibase, iarray, more )

%*****************************************************************************80
%
%% VEC_NEXT generates all N-vectors of integers modulo a given base.
%
%  Discussion:
%
%    The items are produced one at a time.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Parameters:
%
%    Input, integer N, the size of the vectors to be used.
%
%    Input/output, integer IARRAY(N).  On each return from VECNEX,
%    IARRAY will contain entries in the range 0 to IBASE-1.
%
%    Input/output, logical MORE.  Set this variable .FALSE. before
%    the first call.  Normally, MORE will be returned .TRUE. but
%    once all the vectors have been generated, MORE will be
%    reset .FALSE. and you should stop calling the program.
%
%    Input, integer IBASE, the base to be used.  IBASE = 2 will
%    give vectors of 0's and 1's, for instance.
%
  persistent kount
  persistent last

  if ( ~ more )

    kount = 1;
    last = ibase ^ n;
    more = 1;
    iarray(1:n) = 0;

  else

    kount = kount + 1;

    if ( kount == last )
      more = 0;
    end

    iarray(n) = iarray(n) + 1;

    for i = 1 : n

      nn = n - i;

      if ( iarray(nn+1) < ibase )
        return
      end

      iarray(nn+1) = 0;

      if ( nn ~= 0 )
        iarray(nn) = iarray(nn) + 1;
      end

    end

  end

  return
end

repo: jburkardt-m