text:
function view(s, fontsize, fhandle)

%tstoolbox/@signal/view
%   Syntax:
%     * view(signal) (fontsize=12)
%     * view(signal, fontsize)
%     * view(signal, fontsize, figurehandle)
%
%   Signal viewer that decides from the signal's attributes which kind of
%   plot to produce, using the signal's plothint entry to get a hint which
%   kind of plot to produce
%   Possible plothints are:
%     * 'graph'
%     * 'bar'
%     * 'surrbar'
%     * 'surrerrorbar'
%     * 'points'
%     * 'xyplot'
%     * 'xypoints'
%     * 'scatter'
%     * '3dcurve'
%     * '3dpoints'
%     * 'spectrogram'
%     * 'image'
%     * 'multigraph'
%     * 'multipoints'
%     * 'subplotgraph'
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,3);

if nargin < 2
	fontsize = 13;
end

nd = ndim(s);		% Anzahl der Dimensionen mit Laenge groesser eins
dlen = dlens(s);

% check which type of plot we shall produce
pmode = plothint(s); % take this as hint

switch nd			% but check if the desired way of plotting is possible
	case 1
		default = 'graph';
		switch pmode
			case {'graph', 'bar', 'points','surrbar'}
			otherwise
				pmode = default;
		end
    case 2
		default = 'image';
		switch pmode
			case {'xyplot','xypoints', 'scatter'}
				if dlen(2)~=2, pmode = default; end
			case {'3dcurve', '3dpoints'}
				if dlen(2)~=3, pmode = default; end
			case {'spectrogram', 'image', 'multigraph', 'multipoints', 'subplotgraph','surrerrorbar'}
			otherwise
				pmode = default;
		end
    otherwise
	    error('Plotting for this kind of data not yet supported');
		return
end

% set labels, titles etc.
titel = label(s);

% do the plotting
switch pmode
	case 'surrerrorbar'
		x = spacing(s, 1);
		dat = data(s);
		plot(x,dat(:,1));
		hold on
		errorbar(x,dat(:,2),dat(:,3),dat(:,3),'r');
		hold off
		set(gca, 'FontSize', fontsize)
		xlabel(lab(getaxis(s,1)));
		yn = yname(s);
		if ~isempty(yn)
		  ylabel(yn);
		else
			ylabel(label(yunit(s)));
		end
		zoom on
		axis tight;

	case	'graph'
		x = spacing(s, 1);
		dat = data(s);
		if ~isreal(dat)
			phandle = plot(x, real(dat), x, imag(dat));
		else
	    	phandle = plot(x, dat);
		end
	    set(gca, 'FontSize', fontsize)
	    xlabel(lab(getaxis(s,1)));
		yn = yname(s);
		if ~isempty(yn)
	    	ylabel(yn);
	    else
			ylabel(label(yunit(s)));
		end
		zoom on
		axis tight;
	case	'points'
		x = spacing(s, 1);
		dat = data(s);
		if ~isreal(dat)
			phandle = plot(x, real(dat), x, imag(dat), '.');
		else
	    	phandle = plot(x, dat, '.');
		end
	    set(gca, 'FontSize', fontsize)
	    %title(titel);
	    xlabel(lab(getaxis(s,1)));
		yn = yname(s);
		if ~isempty(yn)
	    	ylabel(yn);
	    else
			ylabel(label(yunit(s)));
		end
		zoom on
		axis tight;
	case	'multigraph'
		x = spacing(s, 1);
		dat = data(s);
		phandle = plot(x, dat);
	    set(gca, 'FontSize', fontsize)
	    %title(titel);
	    xlabel(lab(getaxis(s,1)));
		yn = yname(s);
		if ~isempty(yn)
	    	ylabel(yn);
	    else
			ylabel(label(yunit(s)));
		end
		zoom on
		axis tight;
	case	'subplotgraph'
		x = spacing(s, 1);
		N = dlens(s,2); 	% number of graphs = number of subplots
		dat = data(s);
		subplot(N,1,1);
		for i=1:N
			subplot(N,1,i);
			phandle = plot(x, dat(:,i));
	    	set(gca, 'FontSize', ceil(fontsize/sqrt(N)));
	    	xlabel(lab(getaxis(s,1)));
			yn = yname(s);
			if ~isempty(yn)
	    		ylabel(yn);
	    	else
				ylabel(label(yunit(s)));
			end
			zoom on
			axis tight;
		end
	case	'multipoints'
		x = spacing(s, 1);
		dat = data(s);
		phandle = plot(x, dat, '.');
	    set(gca, 'FontSize', fontsize)
	    %title(titel);
	    xlabel(lab(getaxis(s,1)));
		yn = yname(s);
		if ~isempty(yn)
	    	ylabel(yn);
	    else
			ylabel(label(yunit(s)));
		end
		zoom on
		axis tight;

	case 'surrbar'
	        s1=histo(cut(s,1,2,length(data(s))),100);
		x = spacing(s1, 1);
		dat = data(s1);
		bar(x,dat);
		hold on
		plot([data(cut(s,1,1,1)) data(cut(s,1,1,1))],[0 max(dat)],'r','LineWidth',3);
	%		bar(data(cut(s,1,1,1)),1,(max(x)-min(x))/100, 'r');
		hold off
	    set(gca, 'FontSize', fontsize)
	    %title(titel);
	    ylabel(lab(getaxis(s,1)));
		yn = yname(s);
		if ~isempty(yn)
	    	xlabel(yn);
	    else
			ylabel(label(yunit(s)));
		end
		axis tight;
	case 'bar'
		x = spacing(s, 1);
		dat = data(s);
		bar(x,dat);
	    set(gca, 'FontSize', fontsize)
	    %title(titel);
	    xlabel(lab(getaxis(s,1)));
		yn = yname(s);
		if ~isempty(yn)
	    	ylabel(yn);
	    else
			ylabel(label(yunit(s)));
		end
		axis tight;
	case	'xyplot'
		set(gcf, 'Renderer', 'zbuffer')
		tmp = data(s);
		phandle = plot(tmp(:,1),tmp(:,2));
		set(gca, 'FontSize', fontsize);
		%title(titel);
		ylabel(yname(s));
		xlabel(yname(s));
		zoom on
	case	'xypoints'
		set(gcf, 'Renderer', 'zbuffer')
		tmp = data(s);
		phandle = plot(tmp(:,1),tmp(:,2), '.');
		set(gca, 'FontSize', fontsize);
		%title(titel);
		ylabel(yname(s));
		xlabel(yname(s));
		zoom on
	case	'scatter'
		set(gcf, 'Renderer', 'zbuffer')
		tmp = data(s);
		phandle = scatter(tmp(:,1),tmp(:,2));
		set(gca, 'FontSize', fontsize);
		%title(titel);
		ylabel(yname(s));
		xlabel(yname(s));
		zoom on
	case	'3dcurve'
		set(gcf, 'Renderer', 'zbuffer')
		tmp = data(s);
		phandle = plot3(tmp(:,1),tmp(:,2),tmp(:,3));
		set(gca, 'FontSize', fontsize);
		rotate3d on
		set(gca, 'Projection' , 'perspective', 'LineWidth' , [1]);
		grid on
	case	'3dpoints'
		set(gcf, 'Renderer', 'zbuffer')
		tmp = data(s);
		phandle = plot3(tmp(:,1),tmp(:,2),tmp(:,3), '.');
		set(gca, 'FontSize', fontsize);
		rotate3d on
		set(gca, 'Projection' , 'perspective');
		grid on
	case    'image'
		set(gcf, 'Renderer', 'zbuffer')
		v1 = spacing(s,1);
		v2 = spacing(s,2);
		phandle = imagesc(v2, v1, data(s));
		colormap('default')
		set(gca, 'FontSize', fontsize);
		%title(titel);
		set(gca, 'YDir', 'normal');
		ylabel(lab(getaxis(s,1)));
		xlabel(lab(getaxis(s,2)));
		axis tight
		zoom on
	case 	'spectrogram'
		set(gcf, 'Renderer', 'zbuffer')
	    v1 = spacing(s,1);
		v2 = spacing(s,2);
		phandle = pcolor(v2, v1, data(s));
		shading interp
		colormap('default')
		set(gca, 'FontSize', fontsize);
		%title(titel);
		ylabel(lab(getaxis(s,1)));
		xlabel(lab(getaxis(s,2)));
		axis tight
		zoom on
end 					% switch

title(titel);
set(gcf, 'PaperPositionMode', 'auto');	% Plotten immer im Aspect Ratio des wirklichen Bildes

% end of main routine

% helper functions

function r = lab(ax)
if ~isempty(label(ax))
	r = [name(ax) ' (' label(ax)  ')'];
else
	r = name(ax);
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [Theta, k] = estimate_foo(X, s, parameters)

% [Theta k] = estimate_foo(X, s)
%
% DESC:
% Template for the estimation function to be used inside RANSAC
%
% INPUT:
% X                 = 2D point correspondences
% s                 = indices of the points used to estimate the
%                     parameter vector. If empty all the points
%                     are used
% parameters        = a structure of parameters to be used by the function
%
% OUTPUT:
% Theta             = estimated parameter vector Theta = H(:). If
%                     the estimation is not successful return an
%                     empty vector. i.e. Theta = [];
% k                 = dimension of the minimal subset

% here we define the size of the MSS
k = ;

% check if the input parameters are valid
if (nargin == 0) || isempty(X)
    Theta = [];
    return;
end;

% select the points to estimate the parameter vector
if (nargin >= 2) && ~isempty(s)
    X = X(:, s);
end;

% check if we have enough points
N = size(X, 2);
if (N < k)
    error('estimate_foo:inputError', ...
        'At least k point correspondences are required');
end;

% call the estimation function foo
Theta = foo(X, parameters);

% here you may want to check if the results are meaningful.
% If not return an empty vector

return;

author: RANSAC

repo: RANSAC-Toolbox


text:
function F=FPolyKal(T,xDim,order)
%%FPOLYKAL Get the state transition matrix for a linear dynamic model of
%          the given polynomial order (number of derivatives of position
%          included) and number of dimensions (generally 3 for 3D motion).
%          order=1 means constant velocity; order=2 means constant
%          acceleration, etc. The state is ordered in terms of position,
%          velocity, acceleration, etc. This transition matrix is used with
%          discretized continuous-time models as well as with direct
%          discrete-time models.
%
%INPUTS: T The time-duration of the propagation interval.
%     xDim The dimensionality of the target state.
%          xDim is ((order+1)*numDim)X1 dimensional, where numDim is the
%          number of position dimensions of space (e.g. 2D or 3D).
%    order The order >=0 of the filter. If order=1, then it is constant
%          velocity, 2 means constant acceleration, 3 means constant jerk,
%          etc.
%
%OUTPUTS: F The state transition matrix under a linear dynamic model of the
%           given order with motion in numDim dimensions where the state is
%           stacked [position;velocity;acceleration;etc] where the number
%           of derivatives of position depends on the order given. Order=0
%           means just position.
%
%An expression for obtaining the state transition matrix of a
%polynomial-time Kalman filter is given in terms of a Taylor expansion of e
%raised to a matrix in Chapter 4 of [1]. By examining how the state
%transition matrices evolve as the order goes up, (each unit increase in
%order means one more term in the Taylor series expansion remains) one can
%obtain a general rule for an arbitrary dimensionality. For numDim=1, the
%first row is T raised to powers from 0 to order and those terms are
%divided by the factorial of their exponent. For example, T^n/factorial(n)
%for the nth column in the first row. The other rows are just the same as
%the first, but right-shifted (and zero padded).
%
%The state for the 1D case is assumed to be
%[position;velocity;acceleration;etc] --so in the order of increasing
%derivatives. in the multidimensional case, the same order is preserved, so
%position becomes numDim-dimensional as do velocity, acceleration, etc.
%This means that the values for the 1D case just get repeated.
%
%When using discretized continuous-time models, this transition matrix
%generally goes with a covariance matrix generated by QPolyKal. When using
%direct discrete-time models, it generally goes with a covariance matrix
%generated by QPolyKalDirectDisc.
%
%REFERENCES:
%[1] P. Zarchan and H. Musoff, Fundamentals of Kalman Filtering: A
%    Practical Approach, 3rd ed. Reston, VA: American Institute of
%    Aeronautics and  Astronautics, Inc., 2009.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(~isscalar(xDim)&&isvector(xDim))
    %Backwards compatibility for old version.
    xDim = length(xDim);
end
numDim=xDim/(order+1);

%First, create the matrix for 1D motion
numEl=order+1;
F1=zeros(numEl,numEl);

degList=(0:order)';
F1(1,:)=(T.^degList)./factorial(degList);
for curRow=2:numEl
    F1(curRow,curRow:numEl)=F1(1,1:(numEl-curRow+1));
end

%Now, the elements just get spread across identity matrices that are
%numDim dimensional to form a transition matrix of the desired
%dimensionality. This is done usine a Kronecker product.
F=kron(F1,eye(numDim));

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 response = get_net_attr( netfile, maxsize )
    t = caffe.get_net(netfile, 'test');
    for i = maxsize:-1:1
        input = {zeros(i,i,3,1,'single')};
        t.reshape_as_input(input);
        t.reshape();
        response(i) = size(t.blobs(t.outputs{1}).get_data(), 1);
        if response(i) < 1
            break;
        end
    end
    caffe.reset_all;
end

author: liuyuisanai

repo: RSA-for-object-detection


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

for i = 1:size(X, 1)
    for j = 1:K
        distance = norm(X(i,:)-centroids(j,:), 2).^2;
        if j == 1 || distance < min_distance
            min_distance = distance;
            idx(i) = j;
        end
    end
end

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

end

author: imLogM


text:
function Zpv = make_zpv(pv,nb,nl,f,Zb,Yd)
%MAKE_ZPV  Calculates loop impedances for all PV buses.
%
%   Zpv = make_zpv(pv,nb,nl,f,Zb,Yd)
%
%   Loop impedance of a PV bus is defined as impedance of the path between
%   the bus and the slack bus. The mutual impedance between two PV buses is
%   the impedance of the joint part of the two path going from each of the
%   PV buses to the slack bus. The impedances are calculated as bus voltages
%   in cases when at one of the PV buses we inject current of -1 A. All
%   voltages are calculated with the backward-forward sweep method. The
%   input variables are the vector of indicies with "from" buses for each
%   branch, the vector of branch impedances and indicies of PV buses.
%
%   See also CALC_V_PQ_SUM.

npv = length(pv);
Zpv = zeros(npv);
Ye = Yd;
D = zeros(nl,1);
for k = nl:-1:2
    D(k) = 1 / (1 + Zb(k)*Ye(k));
    i = f(k);
    Ye(i) = Ye(i) + D(k)*Ye(k);
end
for ipv = 1:npv
    V = zeros(nb,1);
    Je = zeros(nb,1);
    Je(pv(ipv)) = -1;
    % backward sweep
    for k = nl:-1:2
        i = f(k);
        Je(i) = Je(i) + D(k)*Je(k);
    end
    % forward sweep
    for k = 2:nl
        i = f(k);
        V(k) = D(k) * (V(i)-Zb(k)*Je(k));
    end
    Zpv(:,ipv) = V(pv);
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function write_hea(recordName, fs, datapoints, annotator, gain, offset,unit)

% write_hea(recordName, fs, datapoints, annotator, gain, offset,unit)
%
% ORIGINAL SOURCE AND AUTHORS:
%       This script written by Qiao Li
% COPYRIGHT (C) 2016
% LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information

numsig = 1;
%fs = 125;
%datapoints = 75000;
%recordName = '100';
filename = strcat(recordName, '.', annotator);
unit = 'mV';

fileID = fopen(strcat(recordName, '.hea'),'w');
fprintf(fileID,'%s %d %d %d\n', recordName, numsig, fs, datapoints);
for i = 1:numsig
    fprintf(fileID,'%s 16+24 %d/%s 12\n', filename, gain,unit);
end
fprintf(fileID,'#Creator: HRV_toolbox write_hea.m');
fclose(fileID);

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function varargout = sin(varargin)
%SIN   Sine of a DISKFUN.
%
% See also DISKFUN/SINH and DISKFUN/COS.

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

[varargout{1:nargout}] = sin@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ booleanout ] = mightbexml( filename )
%MIGHTBEXML Test whether a file might be XML or not.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

fid=fopen(filename);
current_char=' ';
while ~feof(fid)&&isspace(current_char)
  current_char=fread(fid,1,'uint8=>char')';
end
fclose(fid);
booleanout=current_char=='<'; % XML file must start with '<' symbol
if booleanout % The first test is extremely quick and filters out most non-XML files
    try % This test is slower, but verifies XML format
        xmlread(filename); % Invalid XML will cause error
    catch % An error will be caught but will still produce red text in command window.  Annoying...
        booleanout=false;
    end
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function hmm = obsinference(data,T,Gamma,hmm,residuals,XX)

N = length(T);
K = length(hmm.state);
order = hmm.train.maxorder;

if ~isstruct(data)
    data = struct('X',data);
    data.C = NaN(size(data.X,1)-order*length(T),K);
end

% Compute residuals
if nargin < 5 || isempty(residuals)
    ndim = size(data.X,2);
    if ~isfield(hmm.train,'Sind')
        orders = formorders(hmm.train.order,hmm.train.orderoffset,hmm.train.timelag,hmm.train.exptimelag);
        hmm.train.Sind = formindexes(orders,hmm.train.S) == 1;
    end
    if ~hmm.train.zeromean, hmm.train.Sind = [true(1,ndim); hmm.train.Sind]; end
    residuals =  getresiduals(data.X,T,hmm.train.S,hmm.train.maxorder,hmm.train.order,...
        hmm.train.orderoffset,hmm.train.timelag,hmm.train.exptimelag,hmm.train.zeromean);
end

% Compute XX, XXGXX
if nargin < 6 || isempty(XX)
    setxx;
end

hmm = obsupdate(Gamma,hmm,residuals,XX,XXGXX);

end

author: OHBA-analysis

repo: HMM-MAR


text:
Dataflow {
        // This is a NVDLA-like dataflow
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}

author: maestro-project

repo: maestro


text:
function Script_CECLM_menpo()

addpath(genpath('../'));

% Replace this with the location of the 300W data location
if(exist([getenv('USERPROFILE') '/Dropbox/AAM/test data/'], 'file'))
    root_test_data = [getenv('USERPROFILE') '/Dropbox/AAM/test data/'];
else
    root_test_data = 'F:\Dropbox\AAM\test data/';
end
[images, detections, labels] = Collect_wild_imgs(root_test_data);

%% loading the CE-CLM model and parameters
[patches, pdm, clmParams, early_term_params] = Load_CECLM_menpo();
% Use the multi-hypothesis model, as bounding box tells nothing about
% orientation
views = [0,0,0; 0,-30,0; 0,30,0; 0,0,30; 0,0,-30;];
views = views * pi/180;

%% Setup recording

num_points = numel(pdm.M)/3;

shapes_all = zeros(size(labels,2),size(labels,3), size(labels,1));
labels_all = zeros(size(labels,2),size(labels,3), size(labels,1));
lhoods = zeros(numel(images),1);
all_lmark_lhoods = zeros(num_points, numel(images));
all_views_used = zeros(numel(images),1);

% Change if you want to visualize the outputs
verbose = true;
output_img = false;

if(output_img)
    output_root = './ceclm_menpo_out/';
    if(~exist(output_root, 'dir'))
        mkdir(output_root);
    end
end
if(verbose)
    f = figure;
end

%% Fitting the model to the provided images

tic
for i=1:numel(images)

    image = imread(images(i).img);
    image_orig = image;

    if(size(image,3) == 3)
        image = rgb2gray(image);
    end

    bbox = detections(i,:);

    % have a multi-view version
    [shape,~,~,lhood,lmark_lhood,view_used] =...
        Fitting_from_bb_multi_hyp(image, [], bbox, pdm, patches, clmParams, views, early_term_params);

    all_lmark_lhoods(:,i) = lmark_lhood;
    all_views_used(i) = view_used;

    shapes_all(:,:,i) = shape;
    labels_all(:,:,i) = labels(i,:,:);

    if(mod(i, 200)==0)
        fprintf('%d done\n', i );
    end

    lhoods(i) = lhood;

    if(output_img)
        v_points = sum(squeeze(labels(i,:,:)),2) > 0;
        DrawFaceOnImg(image_orig, shape, sprintf('%s/%s%d.jpg', output_root, 'fit', i), bbox, v_points);
    end

    if(verbose)
        v_points = sum(squeeze(labels(i,:,:)),2) > 0;
        DrawFaceOnFig(image_orig, shape, bbox, v_points);
    end
end
toc

experiment.errors_normed = compute_error(labels_all, shapes_all + 1.0);
experiment.lhoods = lhoods;
experiment.shapes = shapes_all;
experiment.labels = labels_all;
experiment.all_lmark_lhoods = all_lmark_lhoods;
experiment.all_views_used = all_views_used;

fprintf('Done: mean normed error %.3f median normed error %.4f\n', ...
    mean(experiment.errors_normed), median(experiment.errors_normed));

%%
output_results = 'results/results_ceclm_menpo.mat';
save(output_results, 'experiment');

end

author: TadasBaltrusaitis

repo: OpenFace


text:
%  Figure 7.72      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% LTR Example
%
clf;
F=[0 1; 0 0]; % define system matrices
G=[0;1];
H=[1 0];
J=[0];
sys0=ss(F,G,H,J); % form system
H1=[1 0];
sys=ss(F,G,H1,J);
w=logspace(-1,3,1000); % form frequency vector
rho=1;
Q=rho*H1'*H1;          % define weights
r=1;
[K]=lqr(F,G,Q,r)       % compute LQR gain
sys1=ss(F,G,K,0);      % form loop gain
[maggk1,phasgk1,w]=bode(sys1,w); % Bode plot of loop gain
%pause;
rv=1;
rho=[1 10 100];
% LTR
for i=1:3,
    gam=rho(i)*G;
    Q1=gam'*gam;
    [L]=lqe(F,gam,H,Q1,rv)
    aa=F-G*K-L*H;
    bb=L;
    cc=K;
    dd=0;
    sysk=ss(aa,bb,cc,dd);
    sysgk=series(sys0,sysk);
    [maggk,phsgk,w]=bode(sysgk,w);
    [gm,phm,wcg,wcp]=margin(maggk,phsgk,w)
    subplot(2,1,1);
    loglog(w,[maggk1(:) maggk(:)]);
    xlabel('\omega (rad/sec)');
    ylabel('Magnitude');
    title('Fig. 7.72: loop gain');
    axis([0.1 1000 1e-8 1e3]);
    text(100,10e-8,'q=1');
    text(100,10e-6,'q=10');
    line([0.1 1000],[1 1]);
    line([10.5356 10.5356],[1 0.1]);
    text(11,0.5,'GM');
    text(100,10e-4,'q=100');
    text(300,1e-2,'LQR');
    grid on;
    hold on;
    subplot(2,1,2);
    semilogx(w,[phasgk1(:) phsgk(:)]);
    xlabel('\omega (rad/sec)');
    ylabel('Phase (deg)');
    text(10,-240,'q=1');
    text(10,-170,'q=10');
    text(10,-120,'q=100');
    grid;
    hold on;
 end;
axis([0.1 1000 -270 -90]);
set(gca,'YTick',[-270 -240 -210 -180 -150 -120 -90]);
set(gca,'YTickLabel',[-270 -240 -210 -180 -150 -120 -90]);
line([0.1 1000],[-180 -180]);
line([1.4124 1.4124],[-180+55.07 -180]);
text(2,-150,'PM');
text(1.1,-100,'LQR');

author: Sable

repo: mcbench-benchmarks


text:
function test_suite=test_oddeven_partitioner
% tests for cosmo_oddeven_partitioner
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_oddeven_partitioner_basics
    warning_state=cosmo_warning();
    cleaner=onCleanup(@()cosmo_warning(warning_state));
    cosmo_warning('off');

    for nchunks=[2 6 8]
        ds=cosmo_synthetic_dataset('ntargets',3,'nchunks',nchunks);

        nsamples=size(ds.samples,1);
        rp=randperm(nsamples);
        ds=cosmo_slice(ds,[rp rp]);

        unq_chunks=unique(ds.sa.chunks);
        msk_odd=cosmo_match(ds.sa.chunks,unq_chunks(1:2:end));

        idx_odd=find(msk_odd);
        idx_even=find(~msk_odd);

        fp=struct();
        fp.train_indices={idx_odd, idx_even};
        fp.test_indices={idx_even, idx_odd};
        assert_partitions_equal(fp,cosmo_oddeven_partitioner(ds,'full'));
        %c=ds.sa.chunks;
        %assert_partitions_equal(fp,cosmo_oddeven_partitioner(c,'full'));

        hp=struct();
        hp.train_indices={idx_odd};
        hp.test_indices={idx_even};
        assert_partitions_equal(hp,cosmo_oddeven_partitioner(ds,'half'));
        %c=ds.sa.chunks;
        %assert_partitions_equal(hp,cosmo_oddeven_partitioner(c,'half'));
    end

function test_oddeven_partitioner_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                    cosmo_oddeven_partitioner(varargin{:}),'');
    ds=struct();
    aet(ds);
    ds.samples=zeros(3,4);
    aet(ds);
    ds=cosmo_synthetic_dataset('nchunks',1);
    aet(ds);

    ds=cosmo_synthetic_dataset('nchunks',2);
    aet(ds,'foo');

function assert_partitions_equal(p,q)
    expected_fieldnames={'train_indices';'test_indices'};
    assertEqual(sort(fieldnames(p)),sort(expected_fieldnames));
    assertEqual(sort(fieldnames(p)),sort(fieldnames(q)));

    assert_cell_same_elements(p.train_indices,q.train_indices);
    assert_cell_same_elements(p.test_indices,q.test_indices);

function assert_cell_same_elements(p,q)
    assertEqual(size(p),size(q));
    n=numel(p);
    for k=1:n
        assertEqual(sort(p{k}),sort(q{k}));
    end

function test_unbalanced_partitions()
    warning_state=cosmo_warning();
    cleaner=onCleanup(@()cosmo_warning(warning_state));
    cosmo_warning('off');

    ds=struct();
    for ntargets=2:5
        for nchunks=2:5
            nsamples=ntargets*nchunks;
            ds.samples=randn(nsamples,1);
            ds.sa.chunks=repmat(1:nchunks,1,ntargets)';
            ds.sa.targets=ceil((1:nsamples)'/nchunks);

            idxs=cosmo_randperm(numel(ds.sa.chunks));
            ds=cosmo_slice(ds,idxs);

            % should be ok
            cosmo_oddeven_partitioner(ds);

            idx=ceil(rand()*ntargets*nchunks);
            ds.sa.chunks(idx)=ds.sa.chunks(idx)+1;
            assertExceptionThrown(@() cosmo_oddeven_partitioner(ds),'*');
        end
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function X = sum(varargin)
% SUM (overloaded)

Y = varargin{1};
X = Y;
%X.basis = [];
if nargin == 2 && isequal(varargin{2},length(X.dim))
    % smash slices
    X.basis = kron(ones(1,X.dim(end)),speye(prod(X.dim(1:end-1))))*Y.basis;
    X.dim = X.dim(1:end-1);
else
    if nargin == 1
        index = find(X.dim~=1);
        if isempty(index)
            index = length(X.dim);
        else
            index = index(1);
        end
    else
        index = varargin{2};
        if isequal(index, 'all')
            index = 1:length(X.dim);
        end
        if length(index) > 1
            for i = 1:length(index)
                X = sum(X, index(i));
            end
            return;
        end
    end
    if index > length(X.dim)
        return
    end
    % Permute to the case which we can do fast
    i = 1:length(X.dim);
    p = circshift(i',length(X.dim)-(index))';
    X = permute(X,p);
    % Permute might have squeezed to 2 array
    if length(size(X))~=length(size(Y))
        % Expand back last dimension
        if length(size(X))==2
            X = ndsdpvar(X);
        end
        X.dim = [X.dim ones(1,length(size(Y))-length(size(X)))];
    end
    X = sum(X,length(X.dim));
    if isa(X,'sdpvar')
        X = ndsdpvar(X);
        X.dim = [X.dim ones(1,length(size(varargin{1}))-length(size(X)))];
    else
        X.dim = [X.dim ones(1,length(size(varargin{1}))-length(size(X)))];
    end
    p = circshift((1:length(X.dim))',-(length(Y.dim)-(index)))';
    X = permute(X,p);
   % X.dim = Y.dim;
   % X.dim(index) = 1;
end
if isa(X,'ndsdpvar')
    X.conicinfo = [0 0];
end
X = clean(X);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [fFactorHi, fStdHi, fResHi, fPerHi, fFactorLow, fStdLow, fResLow, fPerLow] = plot_ratefac(mCat1, mCat2, fPeriod1, fPeriod2, fMc1, fMc2)
% [fFactorHi, fStdHi, fResHi, fPerHi, fFactorLow, fStdLow, fResLow, fPerLow] = plot_ratefac(mCat1, mCat2, fPeriod1, fPeriod2, fMc1, fMc2)
% ---------------------------------------------------------------------------------------------------------------------------------------
% Determine rate factor between amount of events between two time periods, differentiated for EQ with M<Mc and M>=Mc.
% Mc = min(fMc1, fMc2);
% N(Per2) = fFac * N (Per1)
% fFac is determined as mean of the factors in the magnitude bins, using a binning of 0.1.A factor for
% magnitudes M<Mc and M>=Mc is determined to explore a difference.Magnitude bins with no earthquakes
% are not used.
%
% Incoming variables:
% mCat1 : EQ catalog period 1
% mCat2 : EQ catalog period 2
% fPeriod1 : Time length period 1 /[dec. year]
% fPeriod2 : Time length period 2 /[dec. year]
% fMc1 : Magnitude of completeness period 1
% fMc2 : Magnitude of completeness period 2
%
% Outgoing variable:
% fFactorHi : Mean rate factor for EQ with M>=Mc
% fStdLow   : Standard deviation of fFactorHi
% fResLow   : Residuum for model application
% fPerHi    : Total percentage of change for EQ with M>=Mc
% fFactorLow : Mean rate factor for EQ with M<Mc
% fStdHi   : Standard deviation offFactorLow
% fResHi   : Residuum for model application
% fPerLow  : Total percentage of change for EQ with M<Mc
%
% J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 26.09.02

% Track of changes:

% Initialization
vMaxMag = [max(mCat1(:,6)) max(mCat2(:,6))];
fMaxMag =  max(vMaxMag);
fMc = min([fMc1 fMc2]);

% Binning
%% EQ with M < Mc
vSel1 = (mCat1(:,6) < fMc);
vSel2 = (mCat2(:,6) < fMc);
mCat1Mc = mCat1(vSel1,:);
mCat2Mc = mCat2(vSel2,:);
[vMag1,vBin1]=histogram(mCat1Mc(:,6),0:0.1:fMaxMag);
[vMag1sort,vBin1sort]=sort(vMag1);
[vMag2,vBin2]=histogram(mCat2Mc(:,6),0:0.1:fMaxMag);
[vMag2sort,vBin2sort]=sort(vMag2);

%% EQ with M >= Mc
mCat1McHi = mCat1(~vSel1,:);
mCat2McHi = mCat2(~vSel2,:);
[vMag1Hi,vBin1Hi]=histogram(mCat1McHi(:,6),0:0.1:fMaxMag);
[vMag1Hisort,vBin1Hisort]=sort(vMag1Hi);
[vMag2Hi,vBin2Hi]=histogram(mCat2McHi(:,6),0:0.1:fMaxMag);
[vMag2Hisort,vBin2Hisort]=sort(vMag2Hi);

%% Overall percentage of change of EQs
fPerLow = length(mCat2Mc(:,1))/length(mCat1Mc(:,1));
fPerHi = length(mCat2McHi(:,1))/length(mCat1McHi(:,1));
fPerLow = 100*fPerLow-100;
fPerHi = 100*fPerHi-100;

% Estimate factor between rates
% Factor for values below Mc of Period 1
vMag1 = vMag1/fPeriod1; % Time normalization
vMag2 = vMag2/fPeriod2; % Time normalization
fRatioLow = vMag2./vMag1;
%vSelRat = 1 - (isnan(fRatioLow) + isinf(fRatioLow));
vSelRatLow = (~isnan(fRatioLow) & ~isinf(fRatioLow));
fFactorLow = mean(fRatioLow(vSelRatLow));
fStdLow = std(fRatioLow(vSelRatLow));
vMagModelow = vMag1*fFactorLow;
fResLow = sqrt(sum((vMag2-(vMagModelow)).^2)/length(vMag2)); % Residuum

% Factor for values above or equal Mc of background
vMag1Hi = vMag1Hi/fPeriod1; % Time normalization
vMag2Hi = vMag2Hi/fPeriod2; % Time normalization
fRatioHi = vMag2Hi./vMag1Hi;
%vSelRatHi = 1 - (isnan(fRatioHi) + isinf(fRatioHi));
vSelRatHi = (~isnan(fRatioHi) & ~isinf(fRatioHi));
fFactorHi = mean(fRatioHi(vSelRatHi));
fStdHi = std(fRatioHi(vSelRatHi));
% vSelStd = (fRatioHi <= (fFactorHi+fStdHi) & fRatioHi >= (fFactorHi-fStdHi));
% fFactorHi = mean(fRatioHi(vSelStd))
vMagModel = vMag1Hi*fFactorHi;
fResHi = sqrt(sum((vMag2Hi-(vMagModel)).^2)/length(vMag2Hi)); % Residuum

%% Figures
figure_w_normalized_uicontrolunits(300);
subplot(3,2,1);
bar(vBin1,vMag1)
ylabel('Period 1');
subplot(3,2,3);
bar(vBin2,vMag2)
ylabel('Period 2');
subplot(3,2,5);
bar(vBin2,vMag2-vMag1)
ylabel('Period 2 - Period 1');
subplot(3,2,2);
bar(vBin1Hi,vMag1Hi)
ylabel('Period 1');
subplot(3,2,4);
bar(vBin2Hi,vMag2Hi)
ylabel('Period 2');
subplot(3,2,6);
bar(vBin2Hi,vMag2Hi-vMag1Hi)
ylabel('Period2 - Period 1');

figure_w_normalized_uicontrolunits(400);
subplot(2,1,1);
bar(vBin1,vMag1*fFactorLow-vMag2);
ylabel('Period2 - Period 1');
subplot(2,1,2);
bar(vBin1,vMag1Hi*fFactorHi-vMag2Hi);
ylabel('Period2 - Period 1');
return

author: CelsoReyes


text:
function D = diagonalMatrixFromVector(d)

    % DIAGONALMATRIXFROMVECTOR coverts a vector into a diagonal matrix.
    %
    % FORMAT: D = diagonalMatrixFromVector(d)
    %
    % INPUT:  - d = [n * 1] vector;
    %
    % OUTPUT: - D = [n * n] diagonal matrix.
    %
    % Authors: Daniele Pucci, Marie Charbonneau, Gabriele Nava
    %
    %          all authors are with the Italian Istitute of Technology (IIT)
    %          email: name.surname@iit.it
    %
    % Genoa, Dec 2017
    %

    %% --- Initialization ---

    D = diag(d);
end

author: robotology

repo: whole-body-controllers


text:
function tf = verifyEFluxExpressionStruct(model,expression)
% Verify the expression struct structure for EFlux
%
% USAGE:
%    tf = verifyEFluxExpressionStruct(model,expression)
%
% INPUTS:
%    model:         The COBRA model struct for the checked expression struct.
%    expression:    The expression struct (fields: value, target,
%                   preprocessed)
%
% OUTPUT:
%    tf:            Whether this struct is valid or not.

tf = true;
if isfield(expression,'preprocessed') && ~expression.preprocessed
    % checkoptional field preprocessed
    % This leads to -1 for unassociated genes
    if ~isempty(setdiff(model.genes,expression.target))
        error('The following genes are lacking assignments:\n%s\nAll genes need to be assigned in order to use eFlux. You can remove genes by using the removeGenesFromModel function',strjoin(setdiff(model.genes,expression.target),'\n'));
    end
end
% check sizes of the value/target fields. They must be of equal size.
if size(expression.value,1) ~= numel(expression.target)
    error('The number of values in the expression struct does not fit the number of  targets in the expression struct');
end

end

author: opencobra

repo: cobratoolbox


text:
%  Figure 10.08      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% fig10_08.m is a script to generate Fig. 10.8  the
% frequency response of the satellite with low-gain PD compensation

np =[0.0360    0.9100];

dp =[1.0000    0.0396    1.0010    0.0000    0.0000];

nc2=0.001*[30 1];

nol2=conv(nc2,np);
dol2=dp;

hold off ; clf
w=logspace(-2,.2);
w(46)=1;
[magol2, phol2]= bode(nol2,dol2,w);
subplot(211) ; loglog(w,magol2); grid; hold on;
xlabel('\omega (rad/sec)');
ylabel('Magnitude |D_2(s)G(s)|');
loglog(w,ones(size(magol2)),'g');
title('Fig. 10.8 Frequency response of low-gain satellite PD design')
phol2a=[phol2, -180*ones(size(phol2))];
subplot(212);  semilogx(w, phol2a); grid; hold on;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
%Bode grid
bodegrid;


text:
function [newVectors, meanValue] = remmean(vectors);
%REMMEAN - remove the mean from vectors
%
% [newVectors, meanValue] = remmean(vectors);
%
% Removes the mean of row vectors.
% Returns the new vectors and the mean.
%
% This function is needed by FASTICA and FASTICAG

% @(#)$Id: remmean.m,v 1.2 2003/04/05 14:23:58 jarmo Exp $

newVectors = zeros (size (vectors));
meanValue = mean (vectors')';
newVectors = vectors - meanValue * ones (1,size (vectors, 2));

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


function [ num ] = salientNumCnt( filesDir, totalNum )
%SALIENTNUMCNT Summary of this function goes here
%   Detailed explanation goes here

    allPoints = zeros(1,totalNum);

    fileExt = '*.txt';

    files = dir(fullfile(filesDir,fileExt));

    for i=0:length(files)-1
        fileCnt = num2str(i);
        fileName = [filesDir, fileCnt, '.txt'];

        file_t = fopen(fileName);
        pointID = fscanf(file_t, '%d');
        fclose(file_t);

        if length(pointID) > 0
            allPoints(sub2ind(size(allPoints), ones(1, length(pointID)), (pointID)')) = 1;
        end

    end

    disp('salient predicted:');
    disp(sum(allPoints));
    disp(sum(allPoints) / length(allPoints));

    num = sum(allPoints);

end


text:
function [yd2] = acre2yd2(acre)
% Convert area from acres to square yards.
% Chad A. Greene 2012
yd2 = acre*4840 ;


text:
function db=cosmo_meeg_layout_collection()
% return supported MEEG channel layouts
%
% db=cosmo_meeg_layout_collection()
%
% Output:
%     db            Nx1 cell for N layouts. Each elements has the following
%                   fields:
%                   .pos      Mx2 x and y coordinates for M channels
%                   .width    Mx1 width of channel
%                   .height   Mx1 height of channel
%                   .label    Mx1 cell string with channel labels
%                   .name     string with filename of layout
%
% Example:
%     % (This example requires FieldTrip)
%     cosmo_skip_test_if_no_external('fieldtrip');
%     %
%     % get all layouts
%     layouts=cosmo_meeg_layout_collection();
%     %
%     % find index of neuromag306cmb layout
%     names=cellfun(@(x)x.name,layouts,'UniformOutput',false);
%     i=find(cosmo_match(names,'neuromag306cmb.lay'));
%     %
%     % show neuromag306 combined planar layout (with 102 channels)
%     pl_layout=layouts{i};
%     cosmo_disp(pl_layout.label)
%     %|| { 'MEG0112+0113'
%     %||   'MEG0122+0123'
%     %||   'MEG0132+0133'
%     %||         :
%     %||   'MEG2642+2643'
%     %||   'COMNT'
%     %||   'SCALE'        }@104x1
%     cosmo_disp([pl_layout.pos pl_layout.width pl_layout.height])
%     %|| [ -0.408     0.273    0.0717    0.0791
%     %||   -0.328     0.306    0.0717    0.0791
%     %||   -0.377     0.179    0.0717    0.0791
%     %||      :         :         :         :
%     %||    0.373    -0.104    0.0717    0.0791
%     %||   -0.547      -0.5    0.0717    0.0791
%     %||    0.547      -0.5    0.0717    0.0791 ]@104x4
%     pl_layout.name
%     %|| 'neuromag306cmb.lay'
%
% Note:
%   - this function requires FieldTrip, as it uses its collection of
%     layouts
%   - the output from this function is similar to FieldTrip's
%     ft_prepare_layout.
%   - this function caches previously read layouts, for optimization
%     reasons. run "clear functions" to reset the cache.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    % cache the layouts, so that after the first call layouts don't have to
    % be read from disc
    persistent cached_layout_db;

    if isempty(cached_layout_db)
        cached_layout_db=read_all_layouts();
    end
    db=cached_layout_db;

function layouts=read_all_layouts()
    cosmo_check_external('fieldtrip');

    % see where fieldtrip is
    ft_dir=fileparts(which('ft_defaults'));
    lay_dir=fullfile(ft_dir,'template','layout');

    % get layout filenames
    lay_fns=dir(fullfile(lay_dir,'*.lay'));
    nlay=numel(lay_fns);

    if nlay==0, error('No layout founds in %s', lay_dir); end

    layout_names={lay_fns.name}';

    % allocate space for output
    layouts=cell(nlay,1);

    % read each layout
    for k=1:nlay
        lay_name=layout_names{k};
        lay_fn=fullfile(lay_dir,lay_name);
        lay=cosmo_meeg_read_layout(lay_fn);
        lay.name=lay_name;
        layouts{k}=lay;
    end

    layouts=apply_fixers(layouts);

function layouts=apply_fixers(layouts)
    fixers={@fix_elec1020_old_fieldtrip,...
            @fix_make_ft_compatible};
    for j=1:numel(fixers)
        fixer=fixers{j};
        layouts=fixer(layouts);
    end

function db=fix_elec1020_old_fieldtrip(db)
    % old fieldtrip had channel locations differently;
    % fix that there so that the unit tests pass
    i=find_layout(db,'elec1020.lay');
    lay=db{i};
    if isequal(lay.pos(1,:),[-0.308949 0.951110]) && ...
                isequal(unique(lay.width),.75)
        suffix=cellfun(@(x)x(end),lay.label)';

        suffix1=suffix=='1';
        suffix2=suffix=='2';

        assert(isequal(lay.label(suffix1),{'Fp1','O1'}'));
        assert(isequal(lay.label(suffix2),{'Fp2','O2'}'));

        lay.pos(suffix1,:)=[-.38 .89111; -.38 -.89111];
        lay.pos(suffix2,:)=[.38  .891004; .38 -.891004];
        lay.width(:)=.35;
        lay.height(:)=.25;
        db{i}=lay;
    end

function i=find_layout(db,name)
    names=cellfun(@(x)x.name,db,'UniformOutput',false);
    i=find(cosmo_match(names,name));
    assert(numel(i)==1);

function layouts=fix_make_ft_compatible(layouts)
    for j=1:numel(layouts)
        layouts{j}=make_single_ft_compatible(layouts{j});
    end

function layout=make_single_ft_compatible(layout)
    processors={@layout_ft_add_outline,...
                @layout_ft_scale,...
                @layout_ft_add_misc_labels};

    for k=1:numel(processors)
        processor=processors{k};
        layout=processor(layout);
    end

function layout=layout_ft_add_outline(layout)
    % add outline and mask to the layout
    alpha=(0:.01:1)'*(2*pi);
    head=[cos(alpha) sin(alpha)]*.5;
    nose=[0.09  0.496;
          0     0.575;
          -0.09 0.496];
    right_ear=[ 0.4970    0.0555
                0.5100    0.0775
                0.5180    0.0783
                0.5299    0.0746
                0.5419    0.0555
                0.5400   -0.0055
                0.5470   -0.0932
                0.5320   -0.1313
                0.5100   -0.1384
                0.4890   -0.1199 ];
    left_ear=bsxfun(@times,[-1 1],right_ear);

    layout.outline={head, nose, right_ear, left_ear};
    layout.mask={head};

function layout=layout_ft_scale(layout)
    % do not consider the COMNT and SCALE channels
    chan_msk=~cosmo_match(layout.label,{'COMNT','SCALE'});

    chan_pos=layout.pos(chan_msk,:);
    extent=(max(chan_pos)-min(chan_pos));

    % put everything in circle around origin with radius .45
    layout.pos=(bsxfun(@rdivide,bsxfun(@minus,...
                    layout.pos,min(chan_pos)),extent))*.9-.45;
    layout.width=layout.width/extent(1);
    layout.height=layout.height/extent(2);

function stacked_layout=stack_channels(layouts)
    % helper function to stack the channels
    % if any element is empty it is just ignored
    keys={'pos','width','height','label'};
    nlayout=numel(layouts);
    nkeys=numel(keys);

    stacked_layout=layouts{1};
    for k=1:nkeys
        key=keys{k};
        values_cell=cell(nlayout,1);
        keep=false(nlayout,1);
        for j=1:nlayout
            layout=layouts{j};
            if ~isempty(layout)
                values_cell{j}=layout.(key);
                keep(j)=true;
            end
        end
        stacked_layout.(key)=cat(1,values_cell{keep});
    end

function layout=layout_ft_add_misc_labels(layout)
    comnt_layout=get_single_label(layout,'COMNT',[-.45 -.45]);
    scale_layout=get_single_label(layout,'SCALE',[ .45 -.45]);
    layout=stack_channels({layout,comnt_layout,scale_layout});

function y=get_single_label(x,label,pos)
    % return a layout
    if cosmo_match({label},x.label);
        y=[];
    else
        y.width=mean(x.width,1);
        y.height=mean(x.height,1);
        y.label={label};
        y.pos=pos;
    end


text:
function obj = configure(obj, config, base, load_path, varargin)
    % configures the robot manipulator objects from the configuration
    % (struct or URDF/YAML) file
    %
    % Parameters:
    % config: the configuration of the joints and links @type char
    % base: the configuration of the floating base
    % coordinates  @type struct

    parser = inputParser;
    addOptional(parser,'removeFixedJoints',false);
    addOptional(parser,'roundN',6);
    parse(parser, varargin{:});
    parser_results = parser.Results;

    eval_math('Needs["RobotManipulator`"];');

    % configure base joint for the model

    if ischar(base)
        base_dofs = get_base_dofs(base);
    elseif isstruct(base)
        fields = {'Name','Type','Offset','R','Axis','Parent','Child'};
        assert(all(isfield(base,fields)),...
            'The base dof structure must have the following fields:\n %s',implode(fields,', '));
        base_dofs = base;
    elseif isempty(base)
        base_dofs = base;
    else
        error(['The second argument must be either a character vector of ',...
            'the type of the floating base coordinates, or a structure array of base DOFs.'])
    end

    if ischar(config) % configuration file name
        % parse the model configuration file
        [folder_name,file_name,file_ext] = fileparts(config);

        % check file type
        if isempty(file_ext) || strcmp(file_ext,'.urdf')
            file_ext = '.urdf'; % default file type is URDF
            config_file = GetFullPath(fullfile(folder_name, [file_name,file_ext]));

            % extract the absolute full file path of the input file
            obj.ConfigFile = config_file;

            % load model from the URDF file
            [name, links, joints, transmissions] = ros_load_urdf(config_file,parser_results.roundN);
            obj.Name = name;
        else
            error('Invalid configuration file type detected: %s.\n',file_ext(2:end));
        end

    elseif isstruct(config)
        obj.Name = config.name;

        fields = {'Name','Mass','Offset','R','Inertia'};
        assert(all(isfield(config.links,fields)),...
            'The links structure must have the following fields:\n %s',implode(fields,', '));
        links = config.links;

        fields = {'Name','Type','Offset','R','Axis','Parent','Child'};
        assert(all(isfield(config.joints,fields)),...
            'The joints structure must have the following fields:\n %s',implode(fields,', '));
        joints = config.joints;
        if isfield(config,'transmissions')
            fields = {'Joint','MechanicalReduction'};
            assert(all(isfield(config.transmissions,fields)),...
                'The transmissions structure must have the following fields:\n %s',implode(fields,', '));
            transmissions = config.transmissions;
        else
            transmissions = [];
        end
    end

    % remove links with zero mass and inertia
    link_indices_to_remove = [];
    for i=1:numel(links)
        link = links(i);
        if link.Mass == 0 && all(all(link.Inertia == zeros(3)))
            link_indices_to_remove = [link_indices_to_remove,i];
        end
    end
    if ~isempty(link_indices_to_remove)
        links(link_indices_to_remove) = [];
    end

    % remove joints with no child links (except base coordinates)
    joint_indices_to_remove = [];
    for i=1:numel(joints)
        joint = joints(i);
        child_index = str_index(joint.Child,{links.Name});
        if isempty(child_index)
            joint_indices_to_remove = [joint_indices_to_remove,i];
        end
    end
    if ~isempty(joint_indices_to_remove)
        joints(joint_indices_to_remove) = [];
    end

    if ~isempty(base_dofs)
        % non-fixed base coordiantes
        base_dofs(end).Child = findBaseLink(obj,joints);
    else
        % fixed base coordinates, we remove the base link (which is not
        % moving anyway) from the list of rigid links
        base_link = findBaseLink(obj,joints);
        % find the index of the base link
        base_link_idx = str_index(base_link,{links.Name});
        % remove the base link
        links(base_link_idx) = [];
    end

    %% Remove fixed joints
    if parser_results.removeFixedJoints

        joints_to_remove = [];
        links_to_remove = [];
        for i = 1:length(joints)
            if strcmp(joints(i).Type, 'fixed')
                joint_to_remove = joints(i);
                R_jot = Rot(joint_to_remove.R); % the rotation matrix of the current joint w.r.t. the parent joint
                p_jot = joint_to_remove.Offset'; % the offset of the current joint in the parent joint frame
                parent_index = find(strcmp({links.Name}, joint_to_remove.Parent));
                child_index = find(strcmp({links.Name}, joint_to_remove.Child));
                % Add mass together
                M = links(parent_index).Mass + links(child_index).Mass;
                % Recompute center of mass
                p_parent = links(parent_index).Offset';  % the offset of parent link CoM in the parent joint
                p_child  = R_jot*links(child_index).Offset' + p_jot; % the offset of child link CoM in the parent joint
                p_com_new = (1.0/M) * (links(parent_index).Mass * p_parent ...
                    +  links(child_index).Mass * p_child);
                % Recompute Inertia matrix
                % Compute inertias about new COM and add together

                R_parent = Rot(links(parent_index).R);    % the rotation matrix of the parent link CoM in the parent joint frame
                R_child  = R_jot*Rot(links(child_index).R); % the rotation matrix of the child link CoM in the parent joint frame

                R1 = (p_com_new - p_parent);
                R2 = (p_com_new - p_child);
                J1 = R_parent*links(parent_index).Inertia*R_parent' + links(parent_index).Mass * (dot(R1,R1)*eye(3) - (R1*R1'));
                J2 = R_child*links(child_index).Inertia*R_child' + links(child_index).Mass * (dot(R2,R2)*eye(3) - (R2*R2'));

                % Update parent link
                links(parent_index).Mass = M;
                links(parent_index).Offset = p_com_new';
                links(parent_index).Inertia = J1 + J2;
                links(parent_index).R = [0,0,0];
                % Find and remove the fixed joint and all children links
                %                 fixed_joint_children_indices = find(strcmp({links.Name}, joints(i).Child));
                %                 for j = 1:length(fixed_joint_children_indices)
                child_link = links(child_index);
                joints_with_fixed_parent = find(strcmp({joints.Parent}, child_link.Name));
                for k = 1:length(joints_with_fixed_parent)
                    joints(joints_with_fixed_parent(k)).Parent = links(parent_index).Name;
                    H01 = [R_jot, p_jot; zeros(1,3), 1];
                    H12 = [Rot(joints(joints_with_fixed_parent(k)).R), joints(joints_with_fixed_parent(k)).Offset'; zeros(1,3), 1];
                    H = H01*H12;
                    joints(joints_with_fixed_parent(k)).Offset = H(1:3, end)';
                    joints(joints_with_fixed_parent(k)).R = H(1:3, 1:3);
                end
                %                 end
                joints_to_remove = horzcat(joints_to_remove, i); %#ok<*AGROW>
                links_to_remove = horzcat(links_to_remove, child_index);
            end
        end
        joints(joints_to_remove) = [];
        links(links_to_remove) = [];
    end

    %%
    dofs = [base_dofs,joints];
    ndof = length(dofs);
    % Add state variables
    x = SymVariable('x',[ndof,1],{dofs.Name});
    dx = SymVariable('dx',[ndof,1],{dofs.Name});
    ddx = SymVariable('ddx',[ndof,1],{dofs.Name});

    obj.addState(x,dx,ddx);

    if ~isempty(transmissions)
        % Add joint actuator inputs
        mechanicalReduction = zeros(1,numel(dofs));
        for i=1:numel(dofs)
            if ~strcmp(dofs(i).Type,'fixed')
                idx = str_index({transmissions.Joint},dofs(i).Name);
                if ~isempty(idx)
                    mechanicalReduction(i) = transmissions(idx).MechanicalReduction;
                    dofs(i).Actuator.Ratio = mechanicalReduction(i);
                    if isfield(transmissions(idx), 'Inertia')
                        dofs(i).Actuator.Inertia = transmissions(idx).Inertia;
                    else
                        dofs(i).Actuator.Inertia = 0;
                    end
                end
            end
        end

        actuated_joints = dofs(mechanicalReduction~=0);
        nact = length(actuated_joints);

        gf = zeros(obj.numState, nact);
        gf(mechanicalReduction~=0,:) = diag(mechanicalReduction(mechanicalReduction~=0));

        u = SymVariable('u',[nact,1],{actuated_joints.Name});
        if isempty(load_path)
            obj.addInput('Control','u',u,gf);
        else
            obj.addInput('Control','u',u,gf,'LoadPath',load_path);
        end
    end

    %% Configure and setup the forward kinematics of the robots
    obj = configureKinematics(obj, dofs, links);

    %% add fixed joints as holonomic constraints
    % find the fixed joint indices
    fixed_joint_indices = str_index('fixed',{obj.Joints.Type});
    % enforce holonomic constraints
    if ~isempty(fixed_joint_indices)
        if isempty(load_path)
            obj = addFixedJoint(obj,fixed_joint_indices);
        else
            obj = addFixedJoint(obj,fixed_joint_indices,load_path);
        end
    end

    function R = Rot(q)
        Rz = @(th) [cos(th), -sin(th), 0; sin(th), cos(th), 0; 0,0,1];
        Ry = @(th) [cos(th), 0, sin(th); 0, 1, 0; -sin(th), 0, cos(th)];
        Rx = @(th) [1,0,0; 0, cos(th), -sin(th); 0, sin(th), cos(th)];
        if isvector(q)
            R = Rz(q(3))*Ry(q(2))*Rx(q(1));
        else
            R = q;
        end
    end
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
%compute the ycoordinate of the head in the central coordinate system

function [data,units]=compute_ycentralhead_mm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
ycentralhead_mm=cell(1,numlarvae);

for i=1:numlarvae
    larva=larvae(i);
    ycentralhead_mm{1,i}=trx(larva).yhead_mm-trx(larva).ycentral_mm;
end

units=parseunits('mm');
data=ycentralhead_mm;

author: kristinbranson

repo: JAABA


text:
function tSeries = rawTSeries(view,scanNum,ROIcoords,detrendOption)
% tSeries = rawTSeries(view,scanNum,ROIcoords,detrendOption)
% PURPOSE:
% Returns the tSeries from all the voxels in the given ROIcoods
% INPUTS:
% scanNum: scan number
% ROIcoords can be either:
%    3xN array of (y,x,z) coords
%    cell array of ROIcoords arrays
% If detrendOption is true (the default) then tSeries are detrended using
% the current settings for that scan prior to being returned. Otherwise the
% percentage tSeries numbers are returned.
% RETURNS: cell array of tSeries matrices. Or a single tSeries matrix
% if only 1 ROI was requested.
% ARW 2004-03-21
% Based on meanTSeries

% Make sure the input is in the right format
if ~exist('ROIcoords','var')
    disp('You must pass in a valid set of ROIcoords');
    return;
end
if ~exist('detrendOption','var')
    detrendOption=1;
end

ROIcoords = ROIcoords2cellArray(view,ROIcoords);

% Find the slice indices for this collection of ROIs
sliceInds = getSlicesROI(view,ROIcoords);

nROIs = length(ROIcoords);

tSeries = cell(1,nROIs);

nFrames = numFrames(view,scanNum);
detrend = detrendFlag(view,scanNum);
smoothFrames = detrendFrames(view,scanNum);

% Take first pass through ROIs to see which slices to load
switch view.viewType
    case {'Inplane' 'Flat'}

        sliceInds = [];
        for r=1:nROIs
            if isempty(ROIcoords{r})
                disp(['MeanTSeries ignoring empty ROI in this slice'])
            else
                sliceInds = [sliceInds, ROIcoords{r}(3,:)];
            end
        end

        sliceInds = unique(sliceInds);
    case {'Gray'}
        sliceInds = 1;
    otherwise
        myErrorDlg('meanTSeries: Only for Inplane, Gray, or Flat views.');
end

% Loop through slices
for iSlice = 1:length(sliceInds);

    slice = sliceInds(iSlice);

    % Load tSeries & divide by mean, but don't detrend yet.
    % Otherwise, detrending the entire tSeries is much slower. DJH
    view = percentTSeries(view,scanNum,slice,0);

    for r=1:nROIs
        % Extract time-series
        subtSeries = getTSeriesROI(view,ROIcoords{r});

        if ~isempty(subtSeries)

            % Detrend now (faster to do it now after extracting subtSeries for a small subset of the voxels)
            if (detrendOption)
                subtSeries = detrendTSeries(subtSeries,detrend,smoothFrames);
            end

            % Stack 'em up
            if isempty(tSeries{r})
                tSeries{r} = [subtSeries];
                numPts{r} = size(subtSeries,2);
            else
                tSeries{r} = [tSeries{r},subtSeries];
                numPts{r} = numPts{r} + size(subtSeries,2);

            end % End check on empty tSeries
        end % end check on empty subTSeries
    end % Next ROI
end % Next slice
% % Final pass through ROIs to turn sum into mean
%
% for r=1:nROIs
%     if isempty(numPts{r})
%         tSeries{r} = zeros(nFrames,1);
%     else
%
%         tSeries{r} = tSeries{r} / numPts{r};
%     end
% end

if (nROIs==1) % Don't mess around with cell arrays if we were only asked for a single ROI
    tSeries = tSeries{1};
end

% Clean up (because didn't detrend the whole tSeries properly)
view.tSeries=[];
view.tSeriesScan=NaN;
view.tSeriesSlice=NaN;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [Theta, k] = estimate_affine(X, s)

% [Theta k] = estimate_affine(X, s)
%
% DESC:
% estimate the parameters of an affine trasformation via least squares.
%
% X2 = A * X1 + b
%
% where Theta = [a11; a21; a12; a22; b1; b2];

%
% VERSION:
% 1.0.0
%
% INPUT:
% X                 = 2D point correspondences
% s                 = indices of the points used to estimate the parameter
%                     vector. If empty all the points are used
%
% OUTPUT:
% Theta             = estimated parameter vector
% k                 = dimension of the minimal subset

% AUTHOR:
% Marco Zuliani, email: marco.zuliani@gmail.com
% Copyright (C) 2009 by Marco Zuliani
%
% LICENSE:
% This toolbox is distributed under the terms of the GNU GPL.
% Please refer to the files COPYING.txt for more information.

% HISTORY:
% 1.0.0             = 08/27/08 - initial version

% cardinality of the MSS
k = 3;

if (nargin == 0) || isempty(X)
    Theta = [];
    return;
end;

if (nargin == 2) && ~isempty(s)
    X = X(:, s);
end;

% check if we have enough points
N = size(X, 2);
if (N < k)
    error('estimate_RST:inputError', ...
        'At least 2 point correspondences are required');
end;

[H A b]= affineLS(X(1:2, :), X(3:4, :));
Theta = [A(:); b(:)];

return;

author: RANSAC

repo: RANSAC-Toolbox


text:
% SYNTAX:
%    flag = flagExpand(flag, expand_size)
%
% DESCRIPTION:
%    expand the flag array
%    if flag is a matrix this flagging expansion will work column by column
%
% INPUT:
%   flag          [n_obs x n_arrays]
%   expand_size   n_epochs with flags to activate at the border of a flagged interval

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

function flag = flagExpand(flag, expand_size)
    % compute a moving window median to filter the data in input

    for c = 1 : size(flag, 2)
        flag(:, c) = conv(single(flag(:, c)), ones(2 * expand_size + 1, 1)', 'same') > 0;
    end
end

author: goGPS-Project


text:
%% Copyright (C) 2016-2017 Lagu
%% Copyright (C) 2017, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym ellipticF (@var{phi}, @var{m})
%% Incomplete elliptic integral of the first kind.
%%
%% The incomplete elliptic integral of the first kind with
%% amplitude @var{phi} and parameter @var{m} is given by:
%% @example
%% @group
%% syms phi m
%% ellipticF (phi, m)
%%   @result{} ans = (sym) F(φ│m)
%% @end group
%%
%% @group
%% rewrite (ans, 'Integral')         % doctest: +SKIP
%%   @result{} ans = (sym)
%%       φ
%%       ⌠
%%       ⎮          1
%%       ⎮ ──────────────────── dα
%%       ⎮    _________________
%%       ⎮   ╱        2
%%       ⎮ ╲╱  - m⋅sin (α) + 1
%%       ⌡
%%       0
%% @end group
%% @end example
%%
%% Example:
%% @example
%% @group
%% vpa (ellipticF (sym (1), sym (-1)))
%%   @result{} (sym) 0.89639378946289458637047451642060
%% @end group
%% @end example
%%
%% For the complete elliptic integral (of the first kind), @pxref{@@sym/ellipticK}.
%%
%% @strong{Note:}
%% this function (and other elliptic integrals in the Symbolic package)
%% follow the Abramowitz and Stegun convention for the ``parameter''
%% @iftex
%% @math{m}.
%% @end iftex
%% @ifnottex
%% @var{m}.
%% @end ifnottex
%% Other sources and software may use different conventions, such as
%% @iftex
%% the ``elliptic modulus'' @math{k}
%% or the ``modular angle'' @math{\alpha},
%% related by @math{m = k^2 = \sin^2(\alpha)}.
%% @end iftex
%% @ifnottex
%% the ``elliptic modulus'' k
%% or the ``modular angle'' α,
%% related by @code{@var{m} = k^2 = sin^2(α)}.
%% @end ifnottex
%% They may define these functions in terms of the sine of the amplitude
%% @iftex
%% @math{\sin(\phi)}.
%% @end iftex
%% @ifnottex
%% @code{sin(@var{phi})}.
%% @end ifnottex
%% For example, Maple uses the elliptic modulus and the sine of the amplitude.
%%
%% @seealso{@@sym/ellipticK, @@sym/ellipticE, @@sym/ellipticPi}
%% @end defmethod

function y = ellipticF (phi, m)

  if (nargin ~= 2)
    print_usage ();
  end

  % y = ellipticPi (0, phi, m);
  y = elementwise_op ('elliptic_f', sym (phi), sym (m));

end

%!error ellipticF (sym(1))
%!error ellipticF (sym(1), 2, 3)

%!assert (double (ellipticF (sym (pi)/3, sym (-105)/10)), 0.6184459461, 10e-11)
%!assert (double (ellipticF (sym (pi)/4, sym (-pi))), 0.6485970495, 10e-11)
%!assert (double (ellipticF (sym (1), sym (-1))), 0.8963937895, 10e-11)
%!assert (double (ellipticF (sym (pi)/6, sym (0))), 0.5235987756, 10e-11)

%!test
%! % compare to Maple
%! us = vpa (ellipticF (sym(11)/10, sym(9)/4), 40);
%! % > evalf(EllipticF(sin(11/10), sqrt(9/4)), 40);
%! maple = vpa ('1.206444996991058996424988192917728014427', 40) - ...
%!         vpa ('0.8157358125823472313001683083685348517476j', 40);
%! assert (abs (double (maple - us)), 0, 1e-39)

repo: octsympy


text:
function sexy(s)
%SEXY   Display a symbolic expression in human readable form.
% SEXY(S) displays the symbolic expression S in a small figure window,
% using standard mathematical notation.
%
% Examples:
%   syms x t positive
%   f=taylor(cos(x));
%   sexy(f)
%   f=int(exp(-t)*t^(x-1),t,0,inf);
%   sexy(f)
%
% Required toolboxes: Symbolic Math
%
% See also ...\symbolic\pretty.

% Author: Naor Movshovitz
% May 2006

if ~isa(s,'sym')
    error('Argument must be of class ''sym''.')
end

S=['$',latex(s),'$'];
S=strrep(S,'&','& \quad');
S=strrep(S,'{\it','\mathrm{');
h=msgbox(S,'Sexy');
h1=get(h,'children');
h2=h1(1);
h3=get(h2,'children');
if isempty(h3)
    h2=h1(2); h3=get(h2,'children');
end
set(h3,'visible','off')
set(h3,'interpreter','latex')
set(h3,'string',S)
set(h3,'fontsize',12)
w=get(h3,'extent');
W=get(h,'position');
W(3)=max(w(3)+10,125);
W(4)=w(4)+40;
set(h,'position',W)
h4=h1(2);
if ~strcmp(get(h4,'tag'),'OKButton'), h4=h1(1); end
o=get(h4,'position');
o(1)=(W(3)-o(3))/2;
set(h4,'position',o)
set(h3,'visible','on')


text:
function test_bug2185

% WALLTIME 00:20:00
% MEM 6gb
% DEPENDENCY ft_sourcegrandaverage ft_selectdata ft_selectdata_new ft_datatype_source ft_math

%% load the end-user provided test data

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/bug2185.mat');
load(filename);

% The content of the file is a cell-array that looks like this
%
% source_timelock_stim{1}
% ans =
%        time: [1x1500 double]
%         pos: [8196x3 double]
%      inside: [8196x1 double]
%     outside: [1x0 double]
%      method: 'average'
%         avg: [1x1 struct]
%         cfg: [1x1 struct]
% source_timelock_stim{1}.avg
% ans =
%          mom: {8196x1 cell}
%          pow: [8196x1500 double]
%     noisecov: {8196x1 cell}

cfg = [];
cfg.parameter = 'pow';
cfg.keepindividual = 'no';
output = ft_sourcegrandaverage(cfg, source_timelock_stim{:});
assert(isfield(output, 'pow'), 'missing output field');
assert(isfield(output, 'time'), 'missing output field');

cfg = [];
cfg.parameter = 'pow';
cfg.operation = 'multiply';
cfg.scalar = 0;
output = ft_math(cfg, output);
assert(all(output.pow(:)==0), 'power should be zero');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% intermezzo, do some interpolation and plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

source2d = output;
source2d.pow(:,:) = 0;
source2d.pow(:,1) = source2d.pos(:,3); % replace the values by something easier to visualize

% create a regular 3d-grid source structure that encompasses the cortical sheet
minpos = floor(min(source2d.pos));
maxpos = ceil(max(source2d.pos));
xgrid = minpos(1):0.5:maxpos(1);
ygrid = minpos(2):0.5:maxpos(2);
zgrid = minpos(3):0.5:maxpos(3);
dim = [length(xgrid) length(ygrid) length(zgrid)];
[X, Y, Z] = ndgrid(xgrid, ygrid, zgrid);
source3d.pos = [X(:) Y(:) Z(:)];
source3d.dim = dim;
% convert the regular 3d-grid source into a volume
volume3d = ft_checkdata(source3d, 'datatype', 'volume');

cfg = [];
cfg.parameter = 'pow';
source3d = ft_sourceinterpolate(cfg, source2d, volume3d);

cfg = [];
cfg.avgovertime = 'yes';
cfg.keeptimedim = 'no';
cfg.parameter = 'pow';
source3davg = ft_selectdata(cfg, source3d);

cfg = [];
cfg.funparameter = 'pow';
ft_sourceplot(cfg, source3davg);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of intermezzo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cfg = [];
cfg.parameter = 'pow';
cfg.keepindividual = 'yes';
output = ft_sourcegrandaverage(cfg, source_timelock_stim{:});
assert( isfield(output, 'pow'), 'missing output field');
assert(isfield(output, 'time'), 'missing output field');

cfg = [];
cfg.parameter = 'pow';
cfg.operation = 'multiply';
cfg.scalar = 0;
output = ft_math(cfg, output);
assert(all(output.pow(:)==0), 'power should be zero');

cfg = [];
cfg.parameter = 'mom';
cfg.keepindividual = 'no';
output = ft_sourcegrandaverage(cfg, source_timelock_stim{:});
assert( isfield(output, 'mom'), 'missing output field');
assert( isfield(output, 'time'), 'missing output field');

cfg = [];
cfg.parameter = 'mom';
cfg.operation = 'multiply';
cfg.scalar = 0;
output = ft_math(cfg, output);
assert(all(output.mom{1}(:)==0), 'moment should be zero');

cfg = [];
cfg.parameter = 'mom';
cfg.keepindividual = 'yes';
output = ft_sourcegrandaverage(cfg, source_timelock_stim{:});
assert( isfield(output, 'mom'), 'missing output field');
assert( isfield(output, 'time'), 'missing output field');

cfg = [];
cfg.parameter = 'mom';
cfg.operation = 'multiply';
cfg.scalar = 0;
output = ft_math(cfg, output);
assert(all(output.mom{1}(:)==0), 'moment should be zero');

cfg = [];
cfg.parameter = 'noisecov';
cfg.keepindividual = 'yes';
output = ft_sourcegrandaverage(cfg, source_timelock_stim{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [forecastmatrix]=maforecastsim(data_endo_a,data_exo_a,data_exo_p,B,Delta,p,n,m,horizon)

% function [forecastmatrix]=maforecastsim(data_endo_a,data_exo_a,data_exo_p,B,Delta,p,n,m,horizon)
% computes the matrix of unconditional forecasts (absent shocks)
% inputs:  - matrix 'data_endo_a': the matrix storing the pre-forecast endogenous data
%          - matrix 'data_exo_a': the matrix storing the pre-forecast exogenous data
%          - matrix 'data_exo_p': the matrix storing the predicted exogenous data
%          - matrix 'B': the matrix containing the VAR coefficients defined in (3.5.10)
%          - matrix 'Delta': the matrix containing the coefficients on exogenous data, defined in (3.5.10)
%          - integer 'p': the number of lags in the model
%          - integer 'n': the number of endogenous variables in the model
%          - integer 'm': the number of exogenous variables in the model
%          - integer 'horizon': the number of periods for which forecasts have to be produced
% outputs: - matrix 'forecastmatrix': the matrix storing the unconditional forecast values

% this function uses the chain rule of forecasts, defined p 38 of the user's guide, adapted for the mean_adjusted BVAR model

% compute the matrix Y
Y=data_endo_a(end-p+1:end,:);

% compute the matrix exo
exo=data_exo_a(end-p+1:end,:);

% repeat the process for periods T+1 to T+h

for jj=1:horizon

% concatenate the transpose of predicted exogenous to the top of actual exogenous
exo=[exo;data_exo_p(jj,:)];

% use the function lagx on Y to obtain the matrix X; retain only the last row
X=bear.lagx(Y,p-1);
X=X(end,:);

% use the function lagx on exo to obtain the matrix Z; retain only the last row
Z=bear.lagx(exo,p);
Z=Z(end,:);
Z(:,m+1:end)=-Z(:,m+1:end);

% obtain predicted value for T+jj
yp=X*B+Z*Delta;

% concatenate the transpose of yp to the top of Y
Y=[Y;yp];

% repeat until values are obtained for T+h
end

% record the values in the matrix forecastmatrix
forecastmatrix=Y(p+1:end,:);

author: european-central-bank

repo: BEAR-toolbox


text:
function varargout = chebcoeffs2(f)
%CHEBCOEFFS2    Bivariate expansion coefficients
%   X = CHEBCOEFFS2(F) returns the matrix of bivariate coefficients such that
%       F= sum_{i=0}^{n-1} ( sum_{j=0}^{n-1} X(i+1,j+1) T_i(y) T_j(x) ).
%
%   [A, D, B] = CHEBCOEFFS2( f ) returns the same coefficients keeping them in
%   low rank form, i.e., X = A * D * B'.
%
% See also PLOTCOEFFS2, CHEBCOEFFS.

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

if ( isempty(f) )
    varargout = { [ ] };
    return
end

if ( iszero(f) )
    varargout = { 0 } ;
    return
end

% Get the low rank representation for f:
[cols, d, rows] = cdr(f);

% Get the coeffs of the rows and the columns:
cols_coeffs = chebcoeffs(cols);
rows_coeffs = chebcoeffs(rows);

if ( nargout <= 1 )
    % Return the matrix of coefficients
    varargout = { cols_coeffs * d * rows_coeffs.' };
elseif ( nargout <= 3 )
    varargout = { cols_coeffs, d, rows_coeffs };
else
    % Two output variables are not allowed.
    error('CHEBFUN:SEPARABLEAPPROX:chebcoeffs2:outputs', ...
        'Incorrect number of outputs.');
end

end

author: chebfun

repo: chebfun


text:
%% Boundary Curvature
%
%%
% The curvature of a curve is defined by fitting localy a circle and taking
% one over its radius. Hence, a straight line will have curvature 0 and a
% circle of radius $2$ will have constant curvature $1/2$ everywhere.
% Hence, the unit of the curvature computed in MTEX is one over the unit of
% the EBSD coordinates which is usually $\mu m$. Let us demonstrate
% boundary curvature use some artifical grain shapes

% import the artificial grain shapes
mtexdata testgrains silent

% select and smooth a few interesting grains
grains = smooth(grains('id',[2 3 9 11 15 16 18 23 31 33 38 40]),10);

plotx2east

%%
% Therefore, we first extract all boundary segments and colorize them
% according to their curvature.

% extract boundary segments
gB = grains.boundary;

% plot some dark background
plot(gB,'linewidth',10,'micronbar','off');

% colorize boundaries by curvature
hold on
plot(gB,gB.curvature,'linewidth',6);
hold off

% set a specific colormap
mtexColorMap('blue2red')
setColorRange(0.25*[-1,1])
mtexColorbar

%%
% Note that all the curvature values are positive. This allways happens if
% the curvature is computed for multiple grains at one. If we consider
% single grains and compute the curvature of single grain boundaries the
% sign of the values indicates whether the grain boundary is convex for
% concav with respect to the specifc grain.

for k = 1:length(grains)

  gB = grains(k).boundary;
  plot(gB,'linewidth',10,'micronbar','off');
  hold on
  plot(gB,gB.curvature,'linewidth',6);

end
hold off

mtexColorMap('blue2red')
setColorRange(0.25*[-1,1])
drawNow(gcm,'figSize',getMTEXpref('figSize'))

%%
% Finally we illustrate the ussage of the <grainBoundary.curvature.html
% curvature> command at a real EBSD map.

% import data and reconstruct grains
mtexdata titanium silent
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'));
grains = smooth(grains,5);

% plot an ipf map
plot(ebsd('indexed'),ebsd('indexed').orientations)

hold on

% plot grain boundaries
plot(grains.boundary,'linewidth',4)

% colorize the grain boundaries of grain 42 according to curvature
plot(grains(42).boundary,grains(42).boundary.curvature(5),'linewidth',6)
hold off
mtexColorMap('blue2red')
setColorRange(0.1*[-1,1])
mtexColorbar

author: mtex-toolbox

repo: mtex

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


text:
function out = cfg_run_file_split(job)

% Split a set of files according to subset indices.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_run_file_split.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

nosel = true(1,numel(job.files));
for k = 1:numel(job.index)
    idx = job.index{k}(job.index{k}<=numel(job.files));
    nosel(idx) = false;
    out{k} = job.files(idx);
end;
out{k+1} = job.files(nosel);

author: spm


text:
%% Ball & Beam Demo
%
% The ball & beam demo was created as a lecture and homework demonstration
% of basic modes of proportional and proportional-derivative control. It
% was originally written in 2001, and most recently tested in 2011 using a current
% version of Matlab. These notes describe basic concepts incorporated in
% the demo and possible homework exercises.

%% Operation
%
% The ballbeam.m file is a function that can be run with no arguments. To
% start the demo simply enter |ballbeam| at the Matlab command prompt.

ballbeam

%% User Interaction
%
% The demo starts up in manual mode. The vertical slider moves the end of
% the beam up and down. The bottom slider adjusts the setpoint indicated by
% the red marker. Push the |Run| button to start the simulation. The ball
% will roll back and forth on the beam as you move the beam position up and
% down. The simulation is stopped by pushing the |Stop| button or when the
% ball rolls off either end of the beam. A summary of the simulation
% results is shown in the plot.
%
% *Exercise*: Move the setpoint to a new place, press |Run|,
% and manipulate the beam position to settle the ball over the new
% setpoint.

%% Plotting Simulation Data
%
% Simulation data is stored in a file |ballbeam.mat| after each run. The
% data can be loaded and plotted as follows:

load ballbeam
clf
subplot(211);
plot(tdata,xdata,'b',tdata,xspdata,'r');
legend('Ball Position','Setpoint','Location','NW');
xlabel('Time (Sec.)');
subplot(212);
plot(tdata,udata);
legend('Beam Position','Location','NW');
xlabel('Time (Sec.)');

%% Dynamics
%
% The position of the ball is governed by a second order differential
% equation of the form
%
% $$ \frac{d^2x}{dt^2} = K_{beam} u$$
%
% *Exercise*: Design and perform a simulation experiment to
% estimate the parameter $K_{beam}$.

%% Proportional Control
%
% As the name suggests, under proportional control the beam position is set
% in proportion to the difference between the ball position and desired
% setpoint.
%
% $$ u = K_p(x_{sp}-x) $$
%
% Under the |Mode| menu, select proportional control (P) and use the slider
% controls to adjust $K_p$. Adjust the setpoint, then push |Run|. You can
% adjust the setpoint and control gain during the course of a |Run|.
%
% *Exercise*: What type of behavior do you observe?  Substituting the the
% expression for the control action into the dynamical model leads to a
% second order differential equation.  What is the solution to that
% equation?
%
% *Exercise*: Use the results of the above exercise to estimate the value
% of $K_{beam}$.

%% Proportional-Derivative Control
%
% Proportional control alone leads to an oscillatory response.
% Proportional-Derivative (PD) control adds
% additional damping through a derivative term in the control law
%
% $$ u = K_p(x_{sp}-x) - K_p\tau_d\frac{dx}{dt} $$
%
% where $\tau_d$ is the derivative 'time-constant'. Under the |Mode| menu,
% select proportional-derivative control (PD) and adjust the parameters
% with the associated sliders.
%
% *Exercise*: Set $K_p = 0.1$. How large does $\tau_d$ need to be to
% completely suppress overshoot of the setpoint?
%
% *Exercise*: Explore the response of the controlled system to a step
% change in the setpoint. Would you prefer an underdamped, critcally
% damped, or overdamped control system? Explain.
%

author: Sable


text:
function [X]=tt_minres_selfprec(A, tol, eps, maxit, prec_type)
%Computation of the approximate TT-matrix inverse using self-prec method
%   X=TT_MINRES_SELFPREC(A,TOL,EPS,MAXIT,PREC_TYPE) Computes the
%   approximate matrix inverse of the TT-matrix A (given in TT1.0 format)
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

matvec='mm+compr';
max_rank=10;

d=max(size(A));
ns=tt_size(A);

I=tt_eye(ns, d); I=core(I);
if (strcmp(prec_type, 'left'))
    X=tt_scal(tt_transp(A), 1e-15);
else
    X=tt_scal(I, 1e-150);
end;

normf=sqrt(tt_dot(tt_mat_to_vec(I), tt_mat_to_vec(I)));
err=1;
err_old=2;
sp=0;

for iout=1:maxit
    if (strcmp(matvec, 'mmk2'))
        resid=tt_mmk2(A, X, eps, max_rank);
    else
        resid=tt_mm(A,X);
        resid=tt_mat_compr(resid, eps,max_rank*2);
    end;
    resid=ttm_add(I, tt_scal(resid, -1));

    if (iout>1)||(strcmp(prec_type, 'left'))
        if (strcmp(matvec, 'mmk2'))
            Xresid=tt_mmk2(X, resid, eps,max_rank);
        else
            Xresid=tt_mm(X,resid);
            Xresid=tt_mat_compr(Xresid, eps,max_rank*2);
        end;
    else
        Xresid=resid;
    end;

    max_xrrank=max(tt_ranks(tt_mat_to_vec(Xresid)));

    if (strcmp(prec_type, 'left')) % we use left preconditioner (igmres)
        if (strcmp(matvec, 'mmk2'))
            w=tt_mmk2(A, Xresid, eps,max_rank);
            w=tt_mmk2(X, w, eps,max_rank);
        else
            w=tt_mm(A,Xresid);
            w=tt_mat_compr(w, eps,max_rank);
            w=tt_mm(X,w);
%             w=tt_mat_compr(w, eps);
        end;

        w = tt_mat_to_vec(w);
        beta=sqrt(tt_dot(tt_mat_to_vec(Xresid),tt_mat_to_vec(Xresid)));
        wr = tt_dot(w, tt_mat_to_vec(Xresid));
        ww = tt_dot(w,w);

        H=zeros(2,1);
        H(1,1)=wr/(beta^2);
        H(2,1)=sqrt(ww/(beta^2)-wr^2/(beta^4));
        rhs=H(1,1)*beta;

        y = rhs/(H'*H);
        err = err*norm(H*y-[beta; 0])/beta;
        y=y/beta;

%         y = tt_dot(w,w);
%         if (y~=0)
%             y = tt_dot(w, tt_mat_to_vec(Xresid))/y;
%         else
%             y=1;
%         end;
    else % Use right preconditioner and fgmres
        if (strcmp(matvec, 'mmk2'))
            w=tt_mmk2(A, Xresid, eps,max_rank);
        else
            w=tt_mm(A,Xresid);
%             w=tt_mat_compr(w, eps);
        end;

        w = tt_mat_to_vec(w);

        beta=sqrt(tt_dot(tt_mat_to_vec(resid),tt_mat_to_vec(resid)));
        wr = tt_dot(w, tt_mat_to_vec(resid));
        ww = tt_dot(w,w);

        H=zeros(2,1);
        H(1,1)=wr/(beta^2);
        H(2,1)=sqrt(ww/(beta^2)-wr^2/(beta^4));
        rhs=H(1,1)*beta;

        y = rhs/(H'*H);
        err = norm(H*y-[beta; 0])/normf;
        y=y/beta;

%         y2=ww;
%         if (y~=0)
%             y2 = wr/y2;
%         else
%             y2=1;
%         end;
    end;

    max_wrank=max(tt_ranks(w));

    Xresid=tt_mat_to_vec(Xresid);

    X=tt_mat_to_vec(X);
    if (err<err_old)
        X = tt_axpy(1, X, y, Xresid, eps, max_rank);
        err_old=err;
    else
        sp=sp+1;
    end;
    max_xrank=max(tt_ranks(X));

    X = tt_vec_to_mat(X, ns, ns);

%     if (strcmp(matvec, 'mmk2'))
%         resid=tt_mmk2(A, X, eps, max_rank);
%     else
%         resid=tt_mm(A, X);
% %         resid=tt_mat_compr(resid, eps);
%     end;
%     resid=ttm_add(I, tt_scal(resid, -1));
%     resid=tt_mat_to_vec(resid);
%
%     err_real=sqrt(tt_dot(resid,resid))/normf;
    err_real=0;
    fprintf('iter [%d], res_real=%3.3e, x_rank=%d, Xr_rank=%d, w_rank=%d, err=%3.3e\n', iout, err_real, max_xrank, max_xrrank, max_wrank, err);
    if (err<tol)||(sp>0)
        break;
    end;

end;

end

author: oseledets

repo: TT-Toolbox


text:
function polygon_integrals_test ( )

%*****************************************************************************80
%
%% POLYGON_INTEGRALS_TEST tests the POLYGON_INTEGRALS library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 October 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( '\n' );
  fprintf ( 'POLYGON_INTEGRALS_TEST:\n' );
  fprintf ( '  MATLAB version\n' );
  fprintf ( '  Test the POLYGON_INTEGRALS library.\n' );

  polygon_integrals_test01 ( );
%
%  Terminate.
%
  fprintf ( '\n' );
  fprintf ( 'POLYGON_INTEGRALS_TEST:\n' );
  fprintf ( '  Normal end of execution.\n' );
  fprintf ( '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function results = run_HOGLR_AFFINE(seq, res_path, bSaveImage)

    project_dir = fileparts(fileparts(fileparts(mfilename('fullpath'))));
    addpath(genpath(fullfile(project_dir, 'trackers', 'HOG_LR')));
    addpath(genpath(fullfile(project_dir, 'trackers', 'HOG_LR', 'FeatureExtractor')));
    addpath(genpath(fullfile(project_dir, 'trackers', 'HOG_LR', 'MotionModel')));
    addpath(genpath(fullfile(project_dir, 'trackers', 'HOG_LR', 'ObservationModel')));
    addpath(genpath(fullfile(project_dir, 'trackers', 'HOG_LR', 'sampler')));
    addpath(genpath(fullfile(project_dir, 'trackers', 'HOG_LR', 'UpdateJudger')));
    addpath(genpath(fullfile(project_dir, 'trackers', 'HOG_LR', 'Utils')));

    config_AFFINE;
    rng(0);
    seq.opt = opt;
    rect=seq.init_rect;
    p = [rect(1)+rect(3)/2, rect(2)+rect(4)/2, rect(3), rect(4)]; % Center x, Center y, height, width
    frame = imread(seq.s_frames{1});
    if size(frame,3)==1
        frame = repmat(frame,[1,1,3]);
    end
    frame = rgb2gray(frame);
    if (seq.opt.useNormalSize)
        scaleHeight = size(frame, 1) / seq.opt.normalHeight;
        scaleWidth = size(frame, 2) / seq.opt.normalWidth;
        p(1) = p(1) / scaleWidth;
        p(3) = p(3) / scaleWidth;
        p(2) = p(2) / scaleHeight;
        p(4) = p(4) / scaleHeight;
    end

    duration = 0;
    tic;
    reportRes = [];

    %% compute feature points
    cvfeat.detector = cv.FeatureDetector('SURF');
    cvfeat.extractor = cv.DescriptorExtractor('SURF');
    cvfeat.matcher = cv.DescriptorMatcher('FlannBased');

    for f = 1:length(seq.s_frames)

        frame = imread(seq.s_frames{f});
        if size(frame,3)==1
            frame = repmat(frame,[1,1,3]);
        end

        if (seq.opt.useNormalSize)
%           frame = imresize(frame, [seq.opt.normalHeight, seq.opt.normalWidth]);
            frame = mexResize(frame, [seq.opt.normalHeight, seq.opt.normalWidth], 'auto');
        end

        %% added code for camera motion estimation
        % resize_im = mexResize(frame, [opt.feature.normalHeight, opt.feature.normalWidth]);
        if f ~= 1
            s_pos = [p(2), p(1)];
            s_size = [p(4), p(3)];
            [cvfeat, H] = find_affine(cvfeat, frame, s_pos, s_size);
        else
            cvfeat.last_keypoints = cvfeat.detector.detect(frame);
            cvfeat.last_descriptors = cvfeat.extractor.compute(frame, cvfeat.last_keypoints);
        end

        frame = im2double(frame);

        if (f ~= 1)
            % first apply the H transformation
            tmpl(:, 1:2) = project_t(tmpl(:, 1:2), H);
            % in case the transformed bounding box get out the image
            % boundary, manually bound them by the size of the bounding box
            tmpl(:, 1) = max(tmpl(:, 1), -tmpl(:, 3)/2);
            tmpl(:, 2) = max(tmpl(:, 2), -tmpl(:, 4)/2);
            tmpl(:, 1) = min(tmpl(:, 1), size(frame, 2)+tmpl(:, 3)/2);
            tmpl(:, 2) = min(tmpl(:, 2), size(frame, 1)+tmpl(:, 4)/2);
            tmpl    = globalParam.MotionModel(tmpl, prob, seq.opt);
            [feat, seq.opt] = globalParam.FeatureExtractor(frame, tmpl, seq.opt);
            prob    = globalParam.ObservationModelTest(feat, model);

            [maxProb, maxIdx] = max(prob);
            p = tmpl(maxIdx, :);
            model.lastOutput = p;
            model.lastProb = maxProb;
        else
            tmpl = globalParam.MotionModel(p, 1, seq.opt);
            prob = ones(1, size(tmpl, 1));
        end

        if (f == 1)
            tmplPos = globalParam.PosSampler(p, seq.opt);
            tmplNeg = globalParam.NegSampler(p, seq.opt);
            [dataPos, seq.opt] = globalParam.FeatureExtractor(frame, tmplPos, seq.opt);
            [dataNeg, seq.opt] = globalParam.FeatureExtractor(frame, tmplNeg, seq.opt);
            model   = globalParam.ObservationModelTrain(dataPos, dataNeg, seq.opt);
            if (seq.opt.useFirstFrame)
                assert(~strcmp(func2str(globalParam.ObservationModelTrain), 'SOSVMTrain'), ...
                    'SOSVM does not support useFirstFrame option!!');
                dataPosFirstFrame = dataPos;
            end
        else
            if (globalParam.ConfidenceJudger(model, seq.opt))
                tmplPos = globalParam.PosSampler(p, seq.opt);
                tmplNeg = globalParam.NegSampler(p, seq.opt);
                [dataPos, seq.opt] = globalParam.FeatureExtractor(frame, tmplPos, seq.opt);
                [dataNeg, seq.opt] = globalParam.FeatureExtractor(frame, tmplNeg, seq.opt);
%                 disp(f);
                if (seq.opt.useFirstFrame)
                    dataPos.feat = [dataPosFirstFrame.feat, dataPos.feat];
                    dataPos.tmpl = [zeros(size(dataPosFirstFrame.tmpl)); dataPos.tmpl];
                end
                model   = globalParam.ObservationModelTrain(dataPos, dataNeg, seq.opt, model);
            end
        end
%
        if opt.visualization
            figure(1),imagesc(frame);
    %         pause(0.1);
            imshow(frame);
            rectangle('position', [p(1) - p(3) / 2, p(2) - p(4) / 2, p(3), p(4)], ...
                'EdgeColor','r', 'LineWidth',2);
            drawnow;
        end
        if (seq.opt.useNormalSize)
            p(1) = p(1) * scaleWidth;
            p(3) = p(3) * scaleWidth;
            p(2) = p(2) * scaleHeight;
            p(4) = p(4) * scaleHeight;
        end
        rect = [p(1) - p(3) / 2, p(2) - p(4) / 2, p(3), p(4)];
        reportRes = [reportRes; round(rect)];
        fprintf('frame %d\n', f);
    end

    duration = duration + toc;
    fprintf('%d frames took %.3f seconds : %.3fps\n',f,duration,f/duration);
    results.res=reportRes;
    results.type='rect';
    results.fps = f/duration;

end

author: flyers

repo: drone-tracking