text:
function varargout = dtmf_story(varargin)
% DTMF_STORY M-file for dtmf_story.fig
%      DTMF_STORY, by itself, creates a new DTMF_STORY or raises the existing
%      singleton*.
%
%      H = DTMF_STORY returns the handle to a new DTMF_STORY or the handle to
%      the existing singleton*.
%
%      DTMF_STORY('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in DTMF_STORY.M with the given input arguments.
%
%      DTMF_STORY('Property','Value',...) creates a new DTMF_STORY or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before dtmf_story_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to dtmf_story_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help dtmf_story

% Last Modified by GUIDE v2.5 16-Jul-2010 19:25:03

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @dtmf_story_OpeningFcn, ...
                   'gui_OutputFcn',  @dtmf_story_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before dtmf_story is made visible.
function dtmf_story_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to dtmf_story (see VARARGIN)

% Choose default command line output for dtmf_story
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes dtmf_story wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = dtmf_story_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes during object creation, after setting all properties.
function axess_CreateFcn(hObject, eventdata, handles)
imshow(imread('toni.jpg'))
% hObject    handle to axess (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axess


text:
function planC = reRasterAndUniformize(planC,scanNumV)
%reRasterAndUniformize.m
%script to delete all existing rasterSegments and re-uniformize the plan
%
%APA, 12/22/2006
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

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

if ~exist('scanNumV','var')
    scanNumV = 1:length(indexS.scan);
end

for scanNum = scanNumV
    planC{indexS.scan}(scanNum).uniformScanInfo         = [];
    planC{indexS.scan}(scanNum).scanArraySuperior       = [];
    planC{indexS.scan}(scanNum).scanArrayInferior       = [];
    planC{indexS.structureArray}(scanNum).indicesArray  = [];
    planC{indexS.structureArray}(scanNum).bitsArray     = [];
    planC{indexS.structureArrayMore}(scanNum).indicesArray  = [];
    planC{indexS.structureArrayMore}(scanNum).bitsArray     = [];
end

%Read CERROptions.json to get import flags
pathStr = getCERRPath;
optName = fullfile(pathStr,'CERROptions.json');
optS    = opts4Exe(optName);
planC{indexS.CERROptions} = optS;

numStructs = length(planC{indexS.structures});
assocScanV = getStructureAssociatedScan(1:numStructs,planC);
strIndV = find(ismember(assocScanV,scanNumV));
%re-generate raster segments
if ~isempty(strIndV)
    [planC{indexS.structures}(strIndV).rasterized] = deal(0);
    planC = getRasterSegs(planC,strIndV);
end

%uniformize
planC = setUniformizedData(planC,[],scanNumV);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%ENCODER PORTION

%here,  fs = sampling frequency
%       aCoeff = LP coefficients
%       pitch = pitch periods
%       v = voiced or unvoiced decision bit
%       g = gain of frames

function [aCoeff, pitch_plot, voiced, gain] = f_ENCODER(x, fs, M);

if (nargin<3), M = 10; end   %prediction order=10;

%INITIALIZATION;
b=1;        %index no. of starting data point of current frame
fsize = 30e-3;    %frame size
frame_length = round(fs .* fsize);   %=number data points in each framesize
                                %of "x"
N= frame_length - 1;        %N+1 = frame length = number of data points in
                            %each framesize

%VOICED/UNVOICED and PITCH;     [independent of frame segmentation]
[voiced, pitch_plot] = f_VOICED (x, fs, fsize);

%FRAME SEGMENTATION for aCoeff and GAIN;
for b=1 : frame_length : (length(x) - frame_length),
    y1=x(b:b+N);     %"b+N" denotes the end point of current frame.
                %"y" denotes an array of the data points of the current
                %frame
    y = filter([1 -.9378], 1, y1);  %pre-emphasis filtering

    %aCoeff [LEVINSON-DURBIN METHOD];
    [a, tcount_of_aCoeff, e] = func_lev_durb (y, M); %e=error signal from lev-durb proc
    aCoeff(b: (b + tcount_of_aCoeff - 1)) = a;  %aCoeff is array of "a" for whole "x"

    %GAIN;
        pitch_plot_b = pitch_plot(b); %pitch period
        voiced_b = voiced(b);
    gain(b) = f_GAIN (e, voiced_b, pitch_plot_b);
end

author: Sable

repo: mcbench-benchmarks


text:
function [Pn] = lomb2(y, t, f,flagNorm)
%
%  [Pn] = lomb(t, y, f)
%
%  Uses Lomb's method to compute normalized
%  periodogram values "Pn" as a function of
%  supplied vector of frequencies "f" for
%  input vectors "t" (time) and "y" (observations).
%  Also returned is probability "Prob" of same
%  length as Pn (and f) that the null hypothesis
%  is valid. If f is not supplied it assumes
%  f =  [1/1024 : 1/1024 : 0.5/min(diff(t))];
%  x and y must be the same length.

% See also:
%
% [1] N.R. Lomb, ``Least-squares frequency analysis of
% unequally spaced data,'' Astrophysical and Space Science,
% (39) pp. 447--462, 1976.   ... and
%
% [2] J.~D. Scargle, ``Studies in astronomical time series analysis.
% II. Statistical aspects of spectral analysis of unevenly spaced data,''
% Astrophysical Journal, vol. 263, pp. 835--853, 1982.
%
% [3] T. Thong, "Lomb-Welch Periodogram for Non-uniform Sampling",
% Proceedings for the 26th anual international conference of the IEEE EMBS,
% Sept 1-5, 2004.

if nargin<4
    flagNorm=true;
end

if nargin<3
 upper_freq = 0.5/min(diff(t));
 f =  [1/1024 : 1/1024 : upper_freq];
end

if nargin < 2
 fprintf('assuming y=diff(t)\n');
 y=diff(t); % RR tachogram?
 t=t(1:length(y)); % shorter by one time stamp now
end

% check inputs
if length(t) ~= length(y);
 error('t and y not same length');
 exit;
end;

% subtract mean, compute variance, initialize Pn
z = y - mean(y);
var = std(y)^2;
N=length(f);
Pn=zeros(size(f));

%	now do main loop for all frequencies
for i=1:length(f)
    w=2*pi*f(i);
    if w > 0
       twt = 2*w*t;
       tau = atan2(sum(sin(twt)),sum(cos(twt)))/2/w;
       wtmt = w*(t - tau);
       Pn(i) = (sum(z.*cos(wtmt)).^2)/sum(cos(wtmt).^2) + ...
		(sum(z.*sin(wtmt)).^2)/sum(sin(wtmt).^2);
     else
	Pn(i) = (sum(z.*t).^2)/sum(t.^2);
     end
end

if flagNorm %normalize by variance
    Pn=Pn./(2*var);
else % return denormalized spectrum (see T. Thong)
    Pn=Pn./length(y);
end

author: jramshur

repo: HRVAS

save_path: github-repos/MATLAB/jramshur-HRVAS


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

%%
clear; close all; clc;

%% Syntax
% |[varargout]=meshBoundary(varargin);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef RWMOP9 < PROBLEM
% <multi> <real> <constrained>
% Four bar plane truss

%------------------------------- Reference --------------------------------
% A. Kumar, G. Wu, M. Ali, Q. Luo, R. Mallipeddi, P. Suganthan, and S. Das,
% A benchmark-suite of real-world constrained multi-objective optimization
% problems and some baseline results, Swarm and Evolutionary Computation,
% 2021, 67: 100961.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Initialization
        function Setting(obj)
            obj.M        = 2;
            obj.D        = 4;
            obj.lower    = [1,sqrt(2),sqrt(2),1];
            obj.upper    = [3,3,3,3];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x  = varargin{1};
            F  = 10; E = 2e5; L = 200; sig = 10;
            x1 = x(:,1);
            x2 = x(:,2);
            x3 = x(:,3);
            x4 = x(:,4);
            % Objective function
            f(:,1) = L.*(2*x1+sqrt(2).*x2+sqrt(2).*x3+x4);
            f(:,2) = F.*L./E.*(2./x1+2.*sqrt(2)./x2-2.*sqrt(2)./x3+2./x4);
            % Constraints
            g = zeros(size(x,1),1);
            Population = SOLUTION(varargin{1},f,g,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
         %% Generate a point for hypervolume calculation
        function R = GetOptimum(obj,~)
            R = [3.0485281e+03   4.0000000e-02];
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [params, names] = ngaussNoiseExtractParam(noise)

% NGAUSSNOISEEXTRACTPARAM Extract parameters from the NGAUSS noise structure.
% FORMAT
% DESC extracts parameters from the noiseless Gaussian
% noise structure into a vector of parameters for optimisation.
% ARG noise : the noise structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the noise. If
% the field 'transforms' is not empty in the noise structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC extracts parameters and parameter names from the noiseless Gaussian
% noise structure.
% ARG noise : the noise structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the noise. If
% the field 'transforms' is not empty in the noise structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings containing names for each
% parameter.
%
% SEEALSO ngaussNoiseParamInit, ngaussNoiseExpandParam, noiseExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005
%
% NOISE

params = [noise.bias];

if nargout > 1
  for i = 1:noise.numProcess
    names{i} = ['bias ' num2str(i)];
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [out] = prep_normalize(dat,varargin)
% prep_normalize (Pre-processing procedure):
%
% Description:
%     This function normalizes the feature vection (Input should be two-dimensional data)
%
% Example:
% [out] = prep_normalize(dat,{'Method','std'});
%
% Input:
%     dat - Feature vector, struct or data itself
% Options:
%     Method - 'std' (default)
%              'max'
%     Type - 'trial'   : Normalizing the data for each trial
%            'feature' : ... for each feature dimension (default)
%
% Seon Min Kim, 05-2016
% seonmin5055@gmail.com

% Options of normalizing methods should be added later

opt = opt_cellToStruct(varargin{:});
if ~isfield(opt,'Method')
    opt.Method = 'std';
end
if ~isfield(opt,'Type')
    opt.Type = 'feature';
end

if isstruct(dat)
    if ~isfield(dat,'x')
        warning('OpenBMI: Data structure must have a field named ''x''');return
    end
    fv=dat.x;
elseif isnumeric(dat)
    fv=dat;
else
    warning('OpenBMI: Check for format of the data');return
end
z = size(fv);
if length(z)~=2
    warning('OpenBMI: Data should be in a form of "feature x trial"');return
end

switch opt.Method
    case 'std'
        if strcmp(opt.Type,'feature')
                c = repmat(std(fv,0,2),[1,z(2)]);
        elseif strcmp(opt.Type,'trial')
                c = repmat(std(fv,0,1),[z(1),1]);
        else
            warning('OpenBMI: Check for the ''Type'' option. It is either ''feature'' or ''trial''');return
        end
        fv = fv.*(1./c);
    case 'max'
        disp('OpenBMI: It might be updated soon')
    otherwise
        warning('OpenBMI: Unknown nomalizing method. It might be updated soon')
end

if isstruct(dat)
    out = rmfield(dat,'x');
    out.x = fv;
elseif isnumeric(dat)
    out = fv;
end

author: PatternRecognition

repo: OpenBMI


text:
function g = linard2KernGradient(kern, x, varargin)

% LINARD2KERNGRADIENT Gradient of LINARD2 kernel's parameters.
% FORMAT
% DESC computes the gradient of functions with respect to the
% automatic relevance determination linear
% kernel's parameters. As well as the kernel structure and the
% input positions, the user provides a matrix PARTIAL which gives
% the partial derivatives of the function with respect to the
% relevant elements of the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x : the input locations for which the gradients are being
% computed.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The argument takes
% the form of a square matrix of dimension  numData, where numData is
% the number of rows in X.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters. The ordering of the vector should match
% that provided by the function kernExtractParam.
%
% FORMAT
% DESC computes the derivatives as above, but input locations are
% now provided in two matrices associated with rows and columns of
% the kernel matrix.
% ARG kern : the kernel structure for which the gradients are being
% computed.
% ARG x1 : the input locations associated with the rows of the
% kernel matrix.
% ARG x2 : the input locations associated with the columns of the
% kernel matrix.
% ARG partial : matrix of partial derivatives of the function of
% interest with respect to the kernel matrix. The matrix should
% have the same number of rows as X1 and the same number of columns
% as X2 has rows.
% RETURN g : gradients of the function of interest with respect to
% the kernel parameters.
%
% SEEALSO linard2KernParamInit, kernGradient, linard2KernDiagGradient, kernGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2009
%
% COPYRIGHT : Michalis K. Titsias, 2009

% KERN

  g = zeros(1, size(x, 2));
  if nargin < 4
    [k, sk] = linard2KernCompute(kern, x);
  else
    [k, sk] = linard2KernCompute(kern, x, varargin{1});
  end
  if nargin < 4
    for i = 1:size(x, 2)
      g(i) =  x(:, i)'*varargin{end}*x(:, i);
    end
  else
    for i = 1:size(x, 2)
      g(i) =  x(:, i)'*varargin{end}*varargin{1}(:, i);
    end
  end
end


text:
function [varargout] = meanDose(planC, varargin)
%Calculate a meanDose metric.
%Stand alone metric :       meanDose(planC, structNum, doseNum, 'Percent' or 'Absolute');
%Request m object   :       meanDose(planC, 'getnewmetric');
%Evaluate m object  :       meanDose(planC, m, 'evaluate');
%LM:  4 Jan 06, JOD: changed interface.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};
optS = planC{indexS.CERROptions};

if(nargin > 1)
    call = varargin{end};
else
    warning('Incorrect Usage, try: meanDose(planC, structNum, doseNum, ''Percent'' or ''Absolute'');');
    return;
end

switch upper(call)
    case 'GETNEWMETRIC'
        m.name = 'meanDose';
        m.valueV = [];
        m.description = 'Returns the mean dose in specified structure';
        m.functionName = @meanDose;
        m.note = '';
        m.params(1).name = 'Structure';
        m.params(1).type = 'DropDown';
        m.params(1).list = {planC{indexS.structures}.structureName};
        m.params(1).value = 1;
        m.params(2).name = 'IsTarget?';
        m.params(2).type = 'DropDown';
        m.params(2).list ={'Target', 'Not Target'};
        m.params(2).value = 2;
        m.units = [];
        m.range = [0 inf];
        m.doseSets = [1];
        varargout = {planC,m};
    return;

    case 'EVALUATE'
        m = varargin{1};
        structName = planC{indexS.structures}(m.params(1).value).structureName;
        m.note = structName;
        m.valueV = [];

        maxD = -inf;
		for i=1:length(m.doseSets)
            [planC, doseBinsV, volsHistV] = getDVHMatrix(planC, m.params(1).value, m.doseSets(i));
            ans = calc_meanDose(doseBinsV, volsHistV, m.params(2).value);
            m.valueV = [m.valueV ans];
            doseArray = getDoseArray(planC{indexS.dose}(m.doseSets(i)));
            maxD = max(maxD, max(doseArray(:)));
        end
        if m.params(2).value == 1
            %Is target, high dose is best.
            m.range = [0 max(m.valueV)];
        elseif m.params(2).value == 2
            m.range = [maxD min(m.valueV)];
        end
        m.units = getDoseUnitsStr(i,planC);
        varargout = {planC,m};
        return;

    otherwise  %No specific call has been made, assume user just wants raw metric.
        if(nargin > 3)
            structNum = varargin{1};
            doseNum = varargin{2};
            volumeTypeString = varargin{3};
            if(strcmp(upper(volumeTypeString), 'ABSOLUTE'))
                volumeType = 2;
            else
                volumeType = 1;
            end
            [planC, doseBinsV, volsHistV] = getDVHMatrix(planC, structNum, doseNum);
            varargout = {calc_meanDose(doseBinsV, volsHistV, volumeType)};
        else
            error('Not enough parameters, try: meanDose(planC, structNum, doseNum, ''Percent'' or ''Absolute'');')
            return;
        end
    end


text:
function [J_c, J_c_dot_nu, holonomicConstraintActivation] = ComputeContactJacobianWithFeetContactFCN(J_L, J_R, J_dot_nu_L, J_dot_nu_R, feet_contact_status)

NDOF = length(J_L(1,7:end));

J_c        = zeros(12,6+NDOF);
J_c_dot_nu = zeros(12,1);
holonomicConstraintActivation = 0.0;
% if left foot in contact
if(feet_contact_status(1))
    J_c (1:6,:)     = J_L;
    J_c_dot_nu(1:6) = J_dot_nu_L;
    holonomicConstraintActivation= 1.0;
end

% if right foot in contact
if(feet_contact_status(2))
    J_c(7:end,:)      = J_R;
    J_c_dot_nu(7:end) = J_dot_nu_R;
    holonomicConstraintActivation = 1.0
end

end

author: robotology

repo: whole-body-controllers


text:
function [fClusterDeg] = calc_ClusterDeg(mCatalog, vCluster)
    % determine degree of clustering
    %
    % [fClusterDeg] = calc_ClusterDeg(mCatalog, vCluster);
    %---------------------------------------------------
    % determine degree of clustering
    %
    % Incoming variables
    % mCatalog : EQ catalog in ZMAP format
    % vCluster : Vector of cluster numbers
    %
    % Outgiong variables:
    % fClusterDeg : Percentage of clustering
    %
    % J. Woessner, woessner@seismo.ifg.ethz.ch
    % updated: 19.08.02

    report_this_filefun();

    vSel = (vCluster(:,1) > 0);
    mCatalogDecl = mCatalog.subset(vSel);
    if isempty(mCatalogDecl) % This means no events in cluster!
        fClusterDeg = NaN;
    else
        fClusterDeg = length(mCatalogDecl(:,1))/mCatalog.Count;
    end
end

author: CelsoReyes


text:
function test_bug1651

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_prepare_mesh ft_datatype_segmentation

readfromdisk = true;

if readfromdisk
  % the segmentation takes quite some time
  % furthermore we want to ensure that we have the old-stype segmentation
  load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1651.mat'));
else
  % this is the original code to create the segmentation
  % create a 3 layered segmentation
  mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.mri'));

  % this speeds up the subsequent stuff
  cfg = [];
  cfg.downsample = 2;
  mri = ft_volumedownsample(cfg, mri);

  cfg = [];
  cfg.output = {'brain', 'scalp', 'skull'};
  seg2 = ft_volumesegment(cfg, mri);
end

% the following results in the error
cfg = [];
cfg.tissue = {'brain', 'skull', 'scalp'};
cfg.numvertices = [3000 2000 1000];
bnd = ft_prepare_mesh(cfg, seg2);

figure
ft_plot_mesh(bnd)
alpha 0.1

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ rhs, ios ] = dlap_rhs_read ( iunit, n )

%*****************************************************************************80
%
%% DLAP_RHS_READ reads the right hand side vector from a DLAP file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Mark Seager,
%    A SLAP for the Masses,
%    Lawrence Livermore National Laboratory,
%    Technical Report UCRL-100267, December 1988.
%
%  Parameters:
%
%    Input, integer IUNIT, the I/O device to be read.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real RHS(N), the right hand side vector.
%
%    Output, integer IOS, the I/O status variable, which is 0 if
%    no I/O error occurred.
%
  ios = 0
  rhs = zeros(n,1);

  for i = 1 : n

    line = fgets ( iunit );

    if ( line == -1 )
      ios = 1;
      break
    end

    [ rhs(i), count ] = sscanf ( line, '%f' );

    if ( count ~= 1 )
      ios = 2;
      break
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function I = quad2d(f, a, b, c, d, varargin)
%QUAD2D  Complete definite integral of SEPARABLEAPPROX.
%   I = QUAD2D( F ), returns the definite integral of a SEPARABLEAPPROX integrated
%   over its domain of definition.
%
%   I = QUAD2D(F, a, b, c, d), returns the definite integral of a SEPARABLEAPPROX.
%   Integrated over the domangle [a b] x [c d].
%
% See also INTEGRAL2, SUM2, INTEGRAL.

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

% Is [a b c d]  subset of the domain of f ?
dom = f.domain;
if ( ( a < dom(1) ) || ( b > dom(2) ) || ( c < dom(3) ) || ( d > dom(4) ) )
    error('CHEBFUN:SEPARABLEAPPROX:quad2d:domain', ...
        'Can only integrate within the SEPARABLEAPPROX''s domain');
end

% Form a new SEPARABLEAPPROX and integrate. (This can be faster.)
I = integral2( restrict( f, [a, b, c, d] ) ); % Any extra arguments are ignored.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% M=2 scattering + freq scatt, mult Q1, mult T, cv parameters

run_name = 'DSS_Table2_TIMIT_m2_freq_multQ1_multT';

src = phone_src('/path/to/timit');

[train_set,test_set,valid_set] = phone_partition(src);

N = 2^13;
T_s = 2560;

filt1_opt.filter_type = {'gabor_1d','morlet_1d'};
filt1_opt.Q = [8 1];
filt1_opt.J = T_to_J(512,filt1_opt);

sc1_opt.M = 2;

ffilt1_opt.filter_type = 'morlet_1d';
ffilt1_opt.J = 6;

fsc1_opt.M = 1;

Wop1 = wavelet_factory_1d(N, filt1_opt, sc1_opt);
fWop1 = wavelet_factory_1d(64, ffilt1_opt, fsc1_opt);

scatt_fun1 = @(x)(log_scat(renorm_scat(scat(x,Wop1))));
fscatt_fun1 = @(x)(func_output(@scat_freq,2,scatt_fun1(x),fWop1));
format_fun1 = @(x)(format_scat(fscatt_fun1(x)));

filt2_opt = filt1_opt;
filt2_opt.Q = [1 1];
filt2_opt.J = T_to_J(512,filt2_opt);

sc2_opt = sc1_opt;

ffilt2_opt = ffilt1_opt;
ffilt2_opt.J = 4;

fsc2_opt = fsc1_opt;

Wop2 = wavelet_factory_1d(N, filt2_opt, sc2_opt);
fWop2 = wavelet_factory_1d(16, ffilt2_opt, fsc2_opt);

scatt_fun2 = @(x)(log_scat(renorm_scat(scat(x,Wop2))));
fscatt_fun2 = @(x)(func_output(@scat_freq,2,scatt_fun2(x),fWop2));
format_fun2 = @(x)(format_scat(fscatt_fun2(x)));

filt3_opt = filt1_opt;
filt3_opt.J = T_to_J(2*512,filt3_opt);

sc3_opt = sc1_opt;

Wop3 = wavelet_factory_1d(N, filt3_opt, sc3_opt);

scatt_fun3 = @(x)(log_scat(renorm_scat(scat(x,Wop3))));
fscatt_fun3 = @(x)(func_output(@scat_freq,2,scatt_fun3(x),fWop1));
format_fun3 = @(x)(format_scat(fscatt_fun3(x)));

filt4_opt = filt2_opt;
filt4_opt.J = T_to_J(2*512,filt4_opt);

sc4_opt = sc2_opt;

Wop4 = wavelet_factory_1d(N, filt4_opt, sc4_opt);

scatt_fun4 = @(x)(log_scat(renorm_scat(scat(x,Wop4))));
fscatt_fun4 = @(x)(func_output(@scat_freq,2,scatt_fun4(x),fWop2));
format_fun4 = @(x)(format_scat(fscatt_fun4(x)));

filt5_opt = filt1_opt;
filt5_opt.J = T_to_J(4*512,filt5_opt);

sc5_opt = sc1_opt;

Wop5 = wavelet_factory_1d(N, filt5_opt, sc5_opt);

scatt_fun5 = @(x)(log_scat(renorm_scat(scat(x,Wop5))));
fscatt_fun5 = @(x)(func_output(@scat_freq,2,scatt_fun5(x),fWop1));
format_fun5 = @(x)(format_scat(fscatt_fun5(x)));

filt6_opt = filt2_opt;
filt6_opt.J = T_to_J(4*512,filt6_opt);

sc6_opt = sc2_opt;

Wop6 = wavelet_factory_1d(N, filt6_opt, sc6_opt);

scatt_fun6 = @(x)(log_scat(renorm_scat(scat(x,Wop6))));
fscatt_fun6 = @(x)(func_output(@scat_freq,2,scatt_fun6(x),fWop2));
format_fun6 = @(x)(format_scat(fscatt_fun6(x)));

duration_fun = @(x,obj)(32*duration_feature(x,obj));

features = {format_fun1, format_fun2, format_fun3, ...
	format_fun4, format_fun5, format_fun6, duration_fun};

for k = 1:length(features)
	fprintf('testing feature #%d...',k);
	tic;
	if nargin(features{k}) == 1
		sz = size(features{k}(randn(N,1)));
	else
		sz = size(features{k}(randn(N,1),struct('u1',1,'u2',N)));
	end
	aa = toc;
	fprintf('OK (%.2fs) (size [%d,%d])\n',aa,sz(1),sz(2));
end

database_opt.input_sz = N;
database_opt.output_sz = T_s;
database_opt.obj_normalize = 2;
database_opt.collapse = 1;

db = prepare_database(src,features,database_opt);
db.features = single(db.features);
db = svm_calc_kernel(db,'gaussian','triangle',[db.indices{train_set}]);

optt.kernel_type = 'gaussian';
optt.gamma = 2.^[-14:2:-10];
optt.C = 2.^[2:2:6];
optt.search_depth = 2;
optt.full_test_kernel = 1;

[dev_err_grid,C_grid,gamma_grid] = ...
	svm_adaptive_param_search(db,train_set,valid_set,optt);

[dev_err,ind] = min(dev_err_grid{end});
C = C_grid{end}(ind);
gamma = gamma_grid{end}(ind);

optt1 = optt;
optt1.C = C;
optt1.gamma = gamma;

model = svm_train(db,train_set,optt1);
labels = svm_test(db,model,test_set);
err = classif_err(labels,test_set,db.src);

save([run_name '.mat'],'labels','err','C','gamma');

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [] = SetLEDsRoomba(serPort, LED,Color, Intensity)
%[] = SetLEDsRoomba(serPort, LED,Color, Intensity)
% Manipulates LEDS
% LED = 0;both off LED=1;Play on LED=2;Advance on LED=3;Both on
% Color determines which color(Red/Green) that the Power LED will
% illuminate as, from 0-100%
% 0 is pure green, 100 is pure red.
% Intensity determines how bright the Power LED appears from 1-100%

% By; Joel Esposito, US Naval Academy, 2011
try

%Flush Buffer
N = serPort.BytesAvailable();
while(N~=0)
fread(serPort,N);
N = serPort.BytesAvailable();
end

warning off
global td
%fwrite(serPort, [139, 2,0, 255])
if(LED==0)
    LED=bin2dec('00000000');
end
if(LED==2)
    LED=bin2dec('00001000');
end
if(LED==1)
    LED=bin2dec('00000010');
end

if (LED==3)
    LED=bin2dec('00001010');
end

L= LED;
aColor= (Color/100)*255;
aIntensity = (Intensity/100)*255;

fwrite(serPort, [139]);  fwrite(serPort,L); fwrite(serPort,aColor);fwrite(serPort,aIntensity);
disp('LEDs Changing')
pause(td)
catch
    disp('WARNING:  function did not terminate correctly.  Output may be unreliable.')
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [ x, w ] = legendre_set_x2_01 ( n )

%*****************************************************************************80
%
%% LEGENDRE_SET_X2_01 sets a Gauss-Legendre rule for X**2 * F(X) on [0,1].
%
%  Discussion:
%
%    The integral:
%
%      Integral ( 0 <= X <= 1 ) X * X * F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( X(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Abramowitz, Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964, page 921.
%
%  Parameters:
%
%    Input, integer N, the order.
%    N must be between 1 and 8.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%
  x = zeros ( n, 1 );
  w = zeros ( n, 1 );

  if ( n == 1 )

    x(1) =   0.75;

    w(1) = 1.0 / 3.0;

  elseif ( n == 2 )

    x(1) = 0.4558481560;
    x(2) = 0.8774851773;

    w(1) = 0.1007858821;
    w(2) = 0.2325474513;

  elseif ( n == 3 )

    x(1) = 0.2949977901;
    x(2) = 0.6529962340;
    x(3) = 0.9270059759;

    w(1) = 0.0299507030;
    w(2) = 0.1462462693;
    w(3) = 0.1571363611;

  elseif ( n == 4 )

    x(1) = 0.2041485821;
    x(2) = 0.4829527049;
    x(3) = 0.7613992624;
    x(4) = 0.9514994506;

    w(1) = 0.0103522408;
    w(2) = 0.0686338872;
    w(3) = 0.1434587898;
    w(4) = 0.1108884156;

  elseif ( n == 5 )

    x(1) = 0.1489457871;
    x(2) = 0.3656665274;
    x(3) = 0.6101136129;
    x(4) = 0.8265196792;
    x(5) = 0.9654210601;

    w(1) = 0.0041138252;
    w(2) = 0.0320556007;
    w(3) = 0.0892001612;
    w(4) = 0.1261989619;
    w(5) = 0.0817647843;

  elseif ( n == 6 )

    x(1) = 0.1131943838;
    x(2) = 0.2843188727;
    x(3) = 0.4909635868;
    x(4) = 0.6975630820;
    x(5) = 0.8684360583;
    x(6) = 0.9740954449;

    w(1) = 0.0018310758;
    w(2) = 0.0157202972;
    w(3) = 0.0512895711;
    w(4) = 0.0945771867;
    w(5) = 0.1073764997;
    w(6) = 0.0625387027;

  elseif ( n == 7 )

    x(1) = 0.0888168334;
    x(2) = 0.2264827534;
    x(3) = 0.3999784867;
    x(4) = 0.5859978554;
    x(5) = 0.7594458740;
    x(6) = 0.8969109709;
    x(7) = 0.9798672262;

    w(1) = 0.0008926880;
    w(2) = 0.0081629256;
    w(3) = 0.0294222113;
    w(4) = 0.0631463787;
    w(5) = 0.0917338033;
    w(6) = 0.0906988246;
    w(7) = 0.0492765018;

  elseif ( n == 8 )

    x(1) = 0.0714910350;
    x(2) = 0.1842282964;
    x(3) = 0.3304477282;
    x(4) = 0.4944029218;
    x(5) = 0.6583480085;
    x(6) = 0.8045248315;
    x(7) = 0.9170993825;
    x(8) = 0.9839022404;

    w(1) = 0.0004685178;
    w(2) = 0.0044745217;
    w(3) = 0.0172468638;
    w(4) = 0.0408144264;
    w(5) = 0.0684471834;
    w(6) = 0.0852847692;
    w(7) = 0.0768180933;
    w(8) = 0.0397789578;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEGENDRE_SET_X2_01 - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of N = %d\n', n );
    error ( 'LEGENDRE_SET_X2_01 - Fatal error!' );

  end

  return
end


text:
function han = make_surface_figure(varargin)

% This is a function to create a figure with surface-brain plots from
% gifti-files (.gii) as available from HCP. The returned handles are passed
% to the function plot_surface_map, which overlays functional data on the
% brain surfaces. plot_surface_map calls this function, so you don't need
% to call it directly.
%
% You need the gifti-toolbox installed
% (http://www.artefact.tk/software/matlab/gifti/, https://github.com/gllmflndn/gifti)
% and HCP S1200 release surface files on your MATLAB path
% (e.g., S1200.L.inflated_MSMAll.32k_fs_LR.surf.gii,
%        S1200.R.inflated_MSMAll.32k_fs_LR.surf.gii), unless you specify
%  your own surface files (see 'surfacefile' input options below).
%
%
%
% :Usage:
% ::
%     h = make_surface_figure(varargin);
%
% See plot_surface_map for plotting functional data.
%
%
%
% :Optional inputs:
%
%   **'surface'**
%       followed by the name of the HCP surface type to use. Will be
%       ignored if you specify 'surfacefile' (see below).
%       options are: inflated [default], midthickness, pial, very_inflated
%
%   **'surfacefiles'**
%       followed by 2x1 cell array with the fullpath to your individually
%       selected surface files, e.g. {path_to_left_hem; path_to_right_hem}.
%       If specified, the 'surface' option above will be ignored. If not
%       specified, the function looks for the HCP files specified with the
%       'surface' option on the matlab path.
%
%   **'facecol'**
%       followed by a [1 x 3] RGB vector specifying the color of the brain
%       surfaces. Default is [0.6 0.6 0.6].
%
%   **'facealpha'**
%       followed by a scalar value for the datamap facealpha value
%       (default=1).
%
%   **'bgcolor'**
%       followed by a [1 x 3] RGB vector setting the figure background color
%
%   **'newfig'**
%       will create a new figure instead of overwriting any existing
%       surface figure (default).
%
% Examples
% see plot_surface_map
%
% Example of GIFTI plotting:
% f = which('S1200.L.inflated_MSMAll.32k_fs_LR.surf.gii')
% s = gifti(f);
% figure; handles = plot(s);
% set(handles, 'FaceColor', [.5 .5 .5], 'FaceAlpha', .4)
%
% ..
%     Author and copyright information:
%     -------------------------------------------------------------------------
%     Copyright (C) 2018 Stephan Geuter
%
%     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/>.
% ..
%

%
% ..
%    Programmers' notes:
%    List dates and changes here, and author of changes
%
%   1/15/2018 - created
%   Stephan Geuter, sgeuter@jhmi.edu
%
%   2/21/2018 - updated surface file options
%   Fred S Barrett & Stephan Geuter
%

newfig = 0;
surftype = 'inflated';
facecol  = [0.6 0.6 0.6];
bgcol    = [1 1 1];
txtcol   = [0 0 0];
facealpha = 1;

if numel(varargin)>0

    % loop input arguments
    for j=1:numel(varargin)

        switch varargin{j}
            case {'surface'}, surftype = varargin{j+1}; varargin{j+1} = '';
            % valid arguments: inflated [default], midthickness, pial, very_inflated

             case {'surfacefile','surfacefiles'}
                 fileL = varargin{j+1}{1};
                 fileR = varargin{j+1}{2};
                 varargin{j+1} = '';
            % if not provided, will default to S1200.[L/R].[surface]_MSMAll.32k_fs_LR.surf.gii

            case {'facecol'}, facecol = varargin{j+1}; varargin{j+1} = '';
            % color for surfaces. default = [0.6 0.6 0.6]

            case {'facealpha'}, facealpha = varargin{j+1}; varargin{j+1} = '';
            % facealpha, default = 1

            case {'bgcolor'}, bgcol = varargin{j+1}; varargin{j+1} = '';
                            if mean(bgcol)<0.5, txtcol = [1 1 1]; end

            case {'newfig','newfigure'}, newfig = 1;
            % make a new figure to plot, default = 0

            otherwise
        end
    end
end

if exist('gifti','file')==0
    warning('''gifti'' function not found on path. Trying to add.');
    addpath(genpath('/Users/sgeuter/Documents/MATLAB/toolbox/gifti'));
end

% check for surface files
if ~exist('fileL','var'), fileL = which(sprintf('S1200.L.%s_MSMAll.32k_fs_LR.surf.gii',surftype)); end
if ~exist('fileR','var'), fileR = which(sprintf('S1200.R.%s_MSMAll.32k_fs_LR.surf.gii',surftype)); end

if isempty(fileL) || isempty(fileR)
    error('Surface files not found on MATLAB path. Check for: \n%s\n%s',...
        sprintf('S1200.L.%s_MSMAll.32k_fs_LR.surf.gii',surftype),...
        sprintf('S1200.R.%s_MSMAll.32k_fs_LR.surf.gii',surftype));
else
    % load surfaces
    hemL = gifti(fileL);
    hemR = gifti(fileR);
end

%% start drawing

f = findobj('Type','figure','-and','Tag','SurfaceFig');
if newfig || isempty(f)
    f = figure('Tag','SurfaceFig');
end
figure(f(1)); clf; set(f,'color',bgcol,'DefaultTextColor',txtcol,'position',[1 600 780 350]);

%%% left lateral
ax(1)=axes('position',[.025 .31 .39 .66]);
p(1) = patch('Faces',hemL.faces,'Vertices',hemL.vertices,'FaceColor',facecol, ...
        'EdgeColor','none','SpecularStrength',.2,'FaceAlpha',facealpha,'SpecularExponent',200);

view(-90,0);
axis off; axis image; material dull;

h(1) = light;
lightangle(h(1),-90,5); % lateral light
lighting gouraud;

%%% left medial
ax(2)=axes('position',[.245 .02 .245 .4]);
p(2) = patch('Faces',hemL.faces,'Vertices',hemL.vertices,'FaceColor',facecol, ...
        'EdgeColor','none','SpecularStrength',.2,'FaceAlpha',facealpha,'SpecularExponent',200);

view(90,0);
axis off; axis image; material dull;

h(2) = light;
lightangle(h(2),90,5); % medial light
lighting gouraud;

%%% right lateral
ax(3)=axes('position',[.585 .31 .39 .66]);
p(3) = patch('Faces',hemR.faces,'Vertices',hemR.vertices,'FaceColor',facecol, ...
        'EdgeColor','none','SpecularStrength',.2,'FaceAlpha',facealpha,'SpecularExponent',200);

view(90,0);
axis off; axis image; material dull;

h(3) = light;
lightangle(h(3),90,5); % lateral light
lighting gouraud;

%%% right medial
ax(4)=axes('position',[.51 .02 .245 .4]);
p(4) = patch('Faces',hemR.faces,'Vertices',hemR.vertices,'FaceColor',facecol, ...
        'EdgeColor','none','SpecularStrength',.2,'FaceAlpha',facealpha,'SpecularExponent',200);

view(-90,0);
axis off; axis image; material dull;

h(4) = light;
lightangle(h(4),-90,5); % medial light
lighting gouraud;

%%%
drawnow;
figure(f);

% output
han.fig = f;
han.ax = ax;
han.obj = p;
han.light = h;
han.label = {'left lateral','left medial','right lateral','right medial'};
han.map = {};

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function plane = fitPlane(points)
%FITPLANE  Fit a 3D plane to a set of points.
%
%   PLANE = fitPlane(POINTS)
%
%   Example
%     pts = randn(300, 3);
%     pts = transformPoint3d(pts, createScaling3d([6 4 2]));
%     pts = transformPoint3d(pts, createRotationOx(pi/6));
%     pts = transformPoint3d(pts, createRotationOy(pi/4));
%     pts = transformPoint3d(pts, createRotationOz(pi/3));
%     pts = transformPoint3d(pts, createTranslation3d([5 4 3]));
%     elli = equivalentEllipsoid(pts);
%     figure; drawPoint3d(pts); axis equal;
%     hold on; drawEllipsoid(elli, ...
%         'drawEllipses', true, 'EllipseColor', 'b', 'EllipseWidth', 3);
%     plane = fitPlane(pts);
%     drawPlane3d(plane, 'm');
%
%   See also
%     planes3d, equivalentEllipsoid, fitLine3d
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2012-11-11, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

% number of points
n = size(points, 1);

% compute centroid
center = mean(points);

% compute the covariance matrix
covPts = cov(points)/n;

% perform a principal component analysis with 2 variables,
% to extract inertia axes
[U, S] = svd(covPts);

% sort axes from greater to lower
[dummy, ind] = sort(diag(S), 'descend'); %#ok<ASGLU>

% format U to ensure first axis points to positive x direction
U = U(ind, :);
if U(1,1) < 0
    U = -U;
    % keep matrix determinant positive
    U(:,3) = -U(:,3);
end

plane = [center U(:,1)' U(:,2)'];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function RightStance2Constraints(nlp, bounds, varargin)

    ip = inputParser;
    ip.addParameter('LoadPath',[],@ischar);
    ip.parse(varargin{:});

    domain = nlp.Plant;

    %% virtual constraints
    opt.constraint.virtual_constraints(nlp, bounds, ip.Results.LoadPath);
    %% foot clearance
    [left_foot_frame] = sys.frames.LeftFoot(domain);
    opt.constraint.foot_clearance(nlp, bounds, left_foot_frame);

    %% swing toe position
    %     trans_opt.constraint.step_distance(nlp, bounds);

    %% swing foot velocity
    opt.constraint.impact_velocity(nlp, bounds, left_foot_frame);

    %% feet distance
    %     opt.constraint.feet_distance(nlp, bounds);

    opt.constraint.yaw_start(nlp, bounds);

    opt.constraint.knee_angle(nlp, bounds);

    %     opt.constraint.average_velocity(nlp, bounds);

    % trans_opt.constraint.periodicity(nlp, floor(nlp.NumNode/2)+1, bounds);
end

author: ayonga

repo: frost-dev

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


text:
function plotSolution(route,model)
city = model.city;
veh = model.veh;
x0 = model.x0;
y0 = model.y0;
x = model.x;
y = model.y;

oneRoute = find(route>city);
From = [0 oneRoute]+1;
To = [oneRoute city+veh]-1;
routes = cell(veh,1);
subplot(1,2,1);

for i = 1:veh
       routes{i} = route(From(i):To(i));
end

colors=hsv(veh);

for r = 1:veh
    X = [x0 x(routes{r}) x0];
    Y = [y0 y(routes{r}) y0];
    color = 0.8*colors(r,:);

     plot(X,Y,'-o',...
            'Color',color,...
            'LineWidth',2,...
            'MarkerSize',10,...
            'MarkerFaceColor','white');
     hold on;
end

plot(x0,y0,'ks',...
        'LineWidth',2,...
        'MarkerSize',18,...
        'MarkerFaceColor','yellow');

  hold off;
  grid on;
  axis equal;

end

author: lzane

repo: VRP-using-SA-with-Matlab


text:
function test_suite=test_testcase_class()
    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_assert_forwarding(t)
    t.assertElementsAlmostEqual(1, 1+eps);
    t.assertEqual(1, 1);
    t.assertError(...
        @()assertEqual([1],'a'), 'assertEqual:classNotEqual');
    t.assertExceptionThrown(...
        @()assertEqual([1],'a'), 'assertEqual:classNotEqual');
    t.assertFalse(true == false);
    t.assertLessThan(0, 1);
    t.assertGreaterThan(1, 0);
    t.assertTrue(false == false);
    t.assertVectorsAlmostEqual(ones(2,1), ones(2,1)+eps);
    t.assertWarning(@()warning('moxunit:warning','msg'),...
                    'moxunit:warning');

function test_undefined_forward(t)
    assertExceptionThrown(@()t.assertDoesNotExist(1, 2),...
                          'moxunit:undefinedMethod');

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function Cdu = computeCdeltau(nu)
    Cu_tmp1 = cell(nu,nu);
    Cu_tmp2 = cell(nu,nu);

    for i=1:nu
       for j=1:nu
          if i==j
              Cu_tmp1{i,j} = eye(2);
              Cu_tmp2{i,j} = -eye(2);
          else
              Cu_tmp1{i,j} = zeros(2,2);
              Cu_tmp2{i,j} = zeros(2,2);
          end
       end
    end

    Cdu = [cell2mat(Cu_tmp1);cell2mat(Cu_tmp2)];
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
function grains = rotate(grains,rot,varargin)
% rotate grains
%
% Syntax
%
%   % rotate the grains about the z-axis by 90*degree
%   grains = rotate(grains,10*degree)
%
%   % rotate about the x-axis
%   grains = rotate(grains,rotation.byAxisAngle(xvector,180*degree))
%
%   % rotate about a specific point
%   grains = rotate(grains,180*degree,'center',[0,0])
%
% Input
%  grains- @grain2d
%  angle - double
%  rot   - @rotation
%
% Options
%  center - [x,y] center of rotation, default is (0,0)
%
% Flags
%  keepXY    - rotate only the orientation data, i.e. the Euler angles
%  keepEuler - rotate only the spatial data, i.e., the x,y, and z values
%
% Output
%  grains - @grain2d

if isa(rot,'double'), rot = rotation.byAxisAngle(vector3d.Z,rot); end

% rotate the orientations
if ~check_option(varargin,'keepEuler')
  grains.prop.meanRotation = quaternion(rot .* grains.prop.meanRotation);
end

% rotate the spatial data
if ~check_option(varargin,'keepXY')

  center = get_option(varargin,'center',[0,0]);

  grains = grains - center;

  % store coordinates as vector3d
  V = vector3d(grains.V(:,1),grains.V(:,2),0);

  % rotate vertices
  V = rot * V;

  % store back
  grains.V = [V.x(:),V.y(:)];

  grains = grains + center;

end

author: mtex-toolbox

repo: mtex

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


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/ra_mmv_phase_transition_snr_30db_s_16.mat';
N = 256;
S = 16;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
pta.SNR = 30;

% options for CoSaMP MMV solver
solver_options.RankAwareResidual = true;
P = 2;

% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K, P, solver_options).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);

author: indigits

repo: sparse-plex


text:
classdef TestGradientDeriche
    %TestGradientDeriche

    methods (Static)
        function test_1
            img = cv.imread(fullfile(mexopencv.root(),'test','cat.jpg'));
            for i=1:2
                if i==2
                    img = cv.cvtColor(img, 'RGB2GRAY');
                end

                Gx = cv.GradientDeriche(img, 'X');
                Gy = cv.GradientDeriche(img, 'Y');
                G = hypot(Gx,Gy);

                validateattributes(Gx, {'numeric'}, {'real', 'size',size(img)});
                validateattributes(Gy, {'numeric'}, {'real', 'size',size(img)});
            end
        end

        function test_error_argnum
            try
                cv.GradientDeriche();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%  INTERNAL FUNCTION: find unique permutations arising from tensors of sums
%
%  ::
%
%    c=find_combos(acell,bcell)
%
%  Args:
%
%     - **acell** [2x1 or 1x2 cell]: where the first element is a number or a
%       char and the second element is an integer describing the number of times
%       the element in the first cell occurs
%     - **bcell** [2x1 or 1x2 cell]: where the first element is a number or a
%       char and the second element is an integer describing the number of times
%       the element in the first cell occurs
%
%  Returns:
%     :
%
%     - **c** [n x k matrix]: char or doubles depending on the type of inputs
%
%  Note:
%
%     - using the Pascal's triangle, we know that
%       (a+b)^5=a^5+5a^4*b+10a^3*b^2+10a^2*b^3+5a*b^4+b^5. Suppose now that a and
%       b are matrices and let's focus on one term, say 10a^3*b^2. We know that
%       we will have 10 combinations of kronecker products in which a appears 3
%       times and b appears twice. But then how to find them? This routine is
%       designed to solve that problem efficiently.
%
%  Example:
%
%     - suppose we want to find all permutations in which u appears 4 times and
%       h appears twice. c=find_combos({'u',4},{'h',2}),size(c)
%     - suppose we want to find all permutations in which 1 appears 4 times and
%       2 appears three times. c=find_combos({1,4},{2,3})
%
%

author: jmaih


text:
function calc_dem_offset()
%CALC_DEM_OFFSET calculate offset of dem
%
%   Andy Hooper, Dec 2007, modified after a script by Judicael Decriem
%
%   ======================================================================
%   06/2009 AH: Large variables cleared from memory
%   10/2009 AH: More large variables cleared from memory
%   ======================================================================

im=combine_amp_dem;
im1=im(:,:,3);
im2=im(:,:,1);
clear im

[s1,s2] = size(im1);
im1 = single(im1);
im2 = single(im2);

a = fft2(im1-mean(mean(im1)));
clear im1
b = ifft2(im2-mean(mean(im2)));
clear im2
c = a.*b;
clear a b
c = ifft2(c);

[C,I1] = max(abs(c));
[C,I2] = max(max(abs(c)));

clear c

if(I2 > s2/2)
   dl = s2-I2+1;
else
   dl = 1-I2;
end
if(I1(I2) > s1/2)
    ds = s1-I1(I2)+1;
else
    ds = 1-I1(I2);
end

fprintf('DEM down  = %d\n',ds)
fprintf('DEM right = %d\n',dl)
figure
plot_amp_dem(ds,dl);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function H=channel2tf(h,delays,freqs);

% H=channel2tf(h,delays,freqs);
%
% Convert a channel specified as an impulse response into a frequency
% response measured at the frequencies specified in 'freqs'.
%
% The expected format of 'h' is n_ue_ant*n_bs_ant*n_delays.
% The format of H is n_ue_ant*n_bs_ant*n_freqs

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% 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/>.

error(nargchk(3,3,nargin));
error(chk_param(h,'H','numeric'));
error(chk_param(delays,'delays','real','vector','horizontal'));
error(chk_param(freqs,'freqs','vector','real'));

n_ant_ue=size(h,1);
n_ant_bs=size(h,2);

H=zeros(n_ant_ue,n_ant_bs,length(freqs));
for k=1:length(delays)
  H(:,:,:)=H(:,:,:)+reshape(kron(exp(-j*2*pi*freqs*delays(k)),h(:,:,k)),n_ant_ue,n_ant_bs,length(freqs));
end

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
%COLNORM Column-wise norm of a matrix
%
% CN = COLNORM(A) is a vector (1xM) comprising the Euclidean norm of each column of the
% matrix A (NxM).
%
% See also norm.

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math
function n = colnorm(a)

	n = sqrt(sum(a.^2));

repo: spatialmath-matlab


text:
function wavelet_test05 ( )

%*****************************************************************************80
%
%% WAVELET_TEST05 tests DAUB10_TRANSFORM and DAUB10_TRANSFORM_INVERSE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WAVELET_TEST05\n' );
  fprintf ( 1, '  DAUB10_TRANSFORM computes the DAUB10 transform of a vector.\n' );
  fprintf ( 1, '  DAUB10_TRANSFORM_INVERSE inverts it.\n' );
%
%  Random data.
%
  n = 16;
  seed = 123456789;
  [ u, seed ] = r8vec_uniform_01 ( n, seed );

  v = daub10_transform ( n, u );

  w = daub10_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)       D10(U)(i)  D10inv(D10(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Constant signal.
%
  n = 8;
  u(1:n) = 1.0;

  v = daub10_transform ( n, u );

  w = daub10_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D10(U)(i)  D10inv(D10(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Linear signal.
%
  n = 16;
  a_first = 1.0;
  a_last = n;
  u = linspace ( a_first, a_last, n );

  v = daub10_transform ( n, u );

  w = daub10_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D10(U)(i)  D10inv(D10(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end
%
%  Quadratic data.
%
  n = 8;
  u(1) = 25.0;
  u(2) = 16.0;
  u(3) = 9.0;
  u(4) = 4.0;
  u(5) = 1.0;
  u(6) = 0.0;
  u(7) = 1.0;
  u(8) = 4.0;

  v = daub10_transform ( n, u );

  w = daub10_transform_inverse ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   i      U(i)        D10(U)(i)  D10inv(D10(U))(i)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %2d  %10.4f  %10.4f  %10.4f\n', i, u(i), v(i), w(i) );
  end

  return
end


text:
function [u, v, p, T, s] = twod_readcad5out(file,e_conn,n_sens)
%-----------------------------------------------------------------------
%  twod_readcad5out.m - Reads in solution data from a *.cad5out file
%                       generated with cadyf.
%
%                     This code is quite problem specific, since the
%                     *.cad5out file is setup to solve the Navier-Stokes
%                     equations.
%
%                     ! NOTE: Matlab has difficulty with D+00 double
%                     precision notation, globally change to E+00
%
%  Copyright (c) 2002, Jeff Borggaard, Virginia Tech
%  Version: 1.0a
%
%  Usage:     [u, v, p, T] = twod_readcad5out(file,e_conn)
%          or
%             [u, v, p   ] = twod_readcad5out(file,e_conn)
%          or
%             [u, v, p, T, s] = twod_readcad5out(file,e_conn,n_sens)
%          or
%             [u, v, p,    s] = twod_readcad5out(file,e_conn,n_sens)
%
%  Variables:     file
%                        Name of the *.cad5out file ('string').
%                 e_conn
%                        Element connectivity (required to calculate
%                        nodal values of pressure)
%                 n_sens
%                        Number of sensitivity variables
%
%                 u,v (,T)   (T is optional)
%                        Nodal information
%                 p
%                        values of pressure nodes on each linear element
%-----------------------------------------------------------------------

[nelem,nel_dof] = size(e_conn);

if ( nargin==2 )
  n_sens = 0;
end

if ( n_sens==0 )
  s = 0;
end

thrmiq = 1;   % 1 if thermal problem, 0 otherwise

iturb = 0;   % if we do turbulent problems, we need to figure out
             % when we have this case (w and w/o thrmiq)

%-----------------------------------------------------------------------
fid = fopen(file);

% skip over beginning stuff
for i=1:137
  line = fgetl(fid);
end

line = fgetl(fid);
while( line(10:15) ~= 'number')
  line = fgetl(fid);
  line = strcat(line,'1                     1');
end

numnp = str2num(line(80:end));

line = fgetl(fid);
line = fgetl(fid);

nelem = str2num(line(80:end));

while ( ~strcmp(line(2:16),'v e l o c i t y') )
  line = fgetl(fid);
  line = strcat(line,'1                     1');
end

for i=1:4
  line = fgetl(fid);
end

% we only assume twod or axi for now
if ( thrmiq )
  for i=1:numnp
    A = fscanf(fid,'%g',8);
    node    = A(1);
    x(node) = A(2);
    y(node) = A(3);
    u(node) = A(5);
    v(node) = A(6);
    T(node) = A(8);
  end
else
  for i=1:numnp
    A = fscanf(fid,'%g',7);
    node = A(1);
    x(i) = A(2);
    y(i) = A(3);
    u(i) = A(5);
    v(i) = A(6);
  end
    T       = 0; % no reason to get the size right
end

for i=1:6
  line = fgetl(fid);
end

% assume ipres = 0
for i=1:nelem
  A = fscanf(fid,'%g',7);
  node = A(1);
  press(i,:) = A(5:7)';
end

% need to convert press (centroid + derivatives) to values at nodes
for i=1:nelem
  nodes_local = e_conn(i,1:3);
  x_local     = x(nodes_local);
  y_local     = y(nodes_local);

  xc          = sum(x_local)/3;
  yc          = sum(y_local)/3;

  p(i,:)      = press(i,1) + press(i,2)*(x_local-xc) + press(i,3)*(y_local-yc);
end

if (iturb)
  for i=1:6
    line = fgetl(fid);
  end

  for i=1:numnp
    A = fscanf(fid,'%g',6);
    k(i) = A(5);
    e(i) = A(6);
  end
end

for np=1:n_sens

  while ( ~strcmp(line(1:20),'SENSITIVITY SOLUTION') )
    line = fgetl(fid);
    line = strcat(line,'1                     1');
  end

  for i=1:9
    line = fgetl(fid);
  end

  % again, assuming either twod or axi
  if (thrmiq)
    for i=1:numnp
      A = fscanf(fid,'%g',4);
      s(np).u(i) = A(2);
      s(np).v(i) = A(3);
      s(np).T(i) = A(4);
    end
  else
    for i=1:numnp
      A = fscanf(fid,'%g',3);
      s(np).u(i) = A(2);
      s(np).v(i) = A(3);
    end
  end

  for i=1:6
    line = fgetl(fid);
  end

  % again, assuming discontinuous pressure
  for i=1:nelem
    A = fscanf(fid,'%g',4);
    node = A(1);
    s_press(i,:) = A(2:4)';
  end

  % need to convert s_press (centroid + derivatives) to values at nodes
  for i=1:nelem
    nodes_local = e_conn(i,1:3);
    x_local     = x(nodes_local);
    y_local     = y(nodes_local);

    xc          = sum(x_local)/3;
    yc          = sum(y_local)/3;

    s(np).p(i,:) = s_press(i,1) + s_press(i,2)*(x_local-xc) + ...
                                  s_press(i,3)*(y_local-yc);
  end

  if (iturb)
    % read in turbulence sensitivities
    for i=1:6
      line = fgetl(fid);
    end

    for i=1:numnp
      A = fscanf(fid,'%g',3);
      s(np).k(i) = A(2);
      s(np).e(i) = A(3);
    end
  end

end


text:
function roi = dtiRoiBuildFromText(fname)
%
% roi = dtiRoiBuildFromText(fname)
%
%
%
% 2007.08.03 RFD wrote it.
%

if(~exist('fname','var')||isempty(fname))
  [f,p] = uigetfile({'*.txt','text-file';'*.*','All files'},'Select a text file of coordinates...');
  if(isnumeric(f))
	disp('User canceled.'); return;
  end
  fname = fullfile(p,f);
end

[p,roiName,e] = fileparts(fname);
coordSpace = 'MNI';
color = 'm';
coords = dlmread(fname);
if(size(coords,2)~=3)
  coords = coords';
end
if(size(coords,2)~=3)
  error('text file must be arranged as three columns.');
end
roi = dtiNewRoi(roiName, color, coords);

outFname = fullfile(p,[roiName '.mat']);
disp(['Saving mrDiffusion roi to ' outFname '...']);
dtiWriteRoi(roi, outFname, [], coordSpace);

if(nargout==0)
    clear roi;
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function f = f11_f0 ( n, x, y )

%*****************************************************************************80
%
%% F11_F0 returns the value of function 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real F(N,1), the function values.
%
  f(1:n,1) = x(1:n,1) .* ( y(1:n,1) + 1.0 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [data,units] = compute_dlength_inmost_wing(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  dlengthl = diff(trx(fly).wing_lengthl_mm,1,2);
  dlengthr = diff(trx(fly).wing_lengthr_mm,1,2);
  data{i} = dlengthl;
  idx = trx(fly).wing_angler(1:end-1) <= -trx(fly).wing_anglel(1:end-1);
  data{i}(idx) = dlengthr(idx);

  data{i} = data{i} ./ trx(fly).dt;
end
units = parseunits('mm/s');

author: kristinbranson

repo: JAABA


text:
%% brute_force_tune
% Code to test the performance of various tuning parameters
% Works sorta like RANSAC I guess?
% Adam Werries 2016, see Apache 2.0 license.

% Specify range
epoch = linspace(0.01,0.1,200);
num_items = length(epoch);
rms_error_filter = Inf*ones(1,num_items);
max_error_filter = Inf*ones(1,num_items);
parfor i = 1:num_items
    fprintf('Iteration: %d, Epoch: %08.7f\n', i, epoch(i));
    temp_filter_time = 0:epoch(i):imu_time(length(imu_time));
    [out_profile,out_IMU_bias_est,out_KF_SD] = ...
        Loosely_coupled_INS_GNSS(init_cond, temp_filter_time, epoch(i), lla, gps, imu, LC_KF_config, est_IMU_bias);
    xyz = out_profile(:,2:4);
    if ~any(any(isnan(xyz))) && ~any(any(isinf(xyz)))
        llh = ecef2lla(xyz);
        [x,y] = deg2utm(llh(:,1),llh(:,2));
        x = x-min_x;
        y = y-min_y;
        ground_truth_full = zeros(length(temp_filter_time), 9);
        ground_truth_full(:,1) = (interp1(gt_t,gt_xyz(:,2),temp_filter_time)-min_x)';
        ground_truth_full(:,2) = (interp1(gt_t,gt_xyz(:,1),temp_filter_time)-min_y)';
        distance = ((ground_truth_full(:,1)-x).^2 + (ground_truth_full(:,2)-y).^2).^0.5;
        rms_error_filter(i) = rms(distance);
        max_error_filter(i) = max(distance);
    end
end

[minmax, i] = min(max_error_filter);
fprintf('\nBest max: %08.7f, rms is %08.7f\n', minmax, rms_error_filter(i));
fprintf('Best iteration for max: %d, epoch: %08.7f\n', i, epoch(i));
[minrms, i] = min(rms_error_filter);
fprintf('Best rms: %08.7f, max is %08.7f\n', minrms, max_error_filter(i));
fprintf('Best iteration for rms: %d, epoch: %08.7f\n', i, epoch(i));

author: awerries

repo: kalman-localization


text:
function f=mandrill()
f  = imreadBW('mandrill.png');
 end

author: scatnet

repo: scatnet


text:
function qNormalised = normaliseImagePoints(q,K,kc)
%function to normalise a set of image points according to the perspective
%intrinsic calibration matrix K and distortion parameters kc. This is the model
%from Oulu university and described clearly in Bouguet's Matlab camera calibration toolbox
% (http://www.vision.caltech.edu/bouguetj/calib_doc/). This is done using
% code from Bouguet's toolbox.
%
%inputs:
%q: a 2*N matrix of un-normalised points in an image.
%
%K: a 3x3 perspective intrinsic calibration matrix
%
%kc: a 5x1 vector of radial and tangential distortion coefficients.
%
%outputs:
%qNormalised: the normalised points
%

Kinv = inv(K);
q(3,:) = 1;
qNormalised = Kinv*q;
qNormalised(1,:) = qNormalised(1,:)./qNormalised(3,:);
qNormalised(2,:) = qNormalised(2,:)./qNormalised(3,:);
qNormalised = qNormalised(1:2,:);

%now compensate for radial and tangential distortion:

if norm(kc)>eps
    k1 = kc(1);
    k2 = kc(2);
    k3 = kc(5);
    p1 = kc(3);
    p2 = kc(4);

    xd = qNormalised;
    x = xd; 				% initial guess

    for kk=1:20,

        r_2 = sum(x.^2);
        k_radial =  1 + k1 * r_2 + k2 * r_2.^2 + k3 * r_2.^3;
        delta_x = [2*p1*x(1,:).*x(2,:) + p2*(r_2 + 2*x(1,:).^2);
            p1 * (r_2 + 2*x(2,:).^2)+2*p2*x(1,:).*x(2,:)];
        x = (xd - delta_x)./(ones(2,1)*k_radial);

    end;
    qNormalised = x;
end

author: tobycollins

repo: IPPE

save_path: github-repos/MATLAB/tobycollins-IPPE


text:
function test_datatype_source

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY

% this function defines a bunch of 'ideal' source structures

% create a set of 3D grid positions
[x,y,z] = ndgrid(-2:5,-3:3,-1:2);
pos     = [x(:) y(:) z(:)];
npos    = size(pos,1);
dim     = size(x);
inside  = false(dim);
inside(2:end-1,2:end-1,:) = true;
outside = find(~inside);
inside  = find(inside);

clear x y z

% create a triangulated mesh
[pnt,tri] = mesh_sphere(162);
npnt      = size(pnt,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the definitions start here

% 3D grid average time domain
source1     = [];
source1.pos = pos;
source1.dim = dim;
source1.time = (-10:89)./100;
source1.inside  = inside(:)';
source1.outside = outside(:)';
source1.mom     = cell(1,npos);
for k = 1:numel(inside)
  source1.mom{inside(k)} = randn(3,100);
end
source1.momdimord = '{pos}_dir_time';
source1.pow       = source1.mom;
for k = 1:numel(inside)
  source1.pow{inside(k)} = sum(abs(source1.mom{inside(k)}).^2);
end
source1.powdimord = '{pos}_time';

% 3D grid multiple trials time domain
source2 = [];
source2.pos = pos;
source2.dim = dim;
source2.time = (-10:89)./100;
source2.inside  = inside(:)';
source2.outside = outside(:)';
source2.mom     = cell(1,npos);
for k = 1:numel(inside)
  source2.mom{inside(k)} = randn(20,3,100);
end
source2.momdimord = '{pos}_rpt_dir_time';
source2.pow       = source2.mom;
for k = 1:numel(inside)
  source2.pow{inside(k)} = squeeze(sum(abs(source2.mom{inside(k)}).^2,2));
end
source2.powdimord = '{pos}_rpt_time';

% 3D grid multiple trials time domain fixed orientation
source3 = [];
source3.pos = pos;
source3.ori = nan(npos,3); %FIXME or should this be 3xnpos?
source3.ori(inside,:) = randn(numel(inside),3);
for k = 1:numel(inside)
  source3.ori(inside(k),:) = source3.ori(inside(k),:)./norm(source3.ori(inside(k),:));
end
source3.dim = dim;
source3.time = (-10:89)./100;
source3.inside  = inside(:)';
source3.outside = outside(:)';
source3.mom     = cell(1,npos);
for k = 1:numel(inside)
  source3.mom{inside(k)} = randn(20,100);
end
source3.momdimord = '{pos}_rpt_time';
source3.pow       = source2.mom;
for k = 1:numel(inside)
  source3.pow{inside(k)} = sum(abs(source3.mom{inside(k)}).^2,2);
end
source3.powdimord = '{pos}_rpt_time';

% cortical sheet average time domain
source4 = [];
source4.pos = pnt;
source4.tri = tri;
source4.time = (-10:89)./100;
source4.inside  = 1:npnt;
source4.outside = [];
source4.mom     = cell(1,npnt);
inside = source4.inside;
for k = 1:numel(inside)
  source4.mom{inside(k)} = randn(3,100);
end
source4.momdimord = '{pos}_time';
source4.pow       = zeros(npnt,100);
for k = 1:numel(inside)
  source4.pow(k,:) = sum(abs(source4.mom{inside(k)}).^2);
end
source4.powdimord = 'pos_time';

% cortical sheet multiple trials
source5 = [];
source5.pos = pnt;
source5.tri = tri;
source5.time = (-10:89)./100;
source5.inside  = 1:npnt;
source5.outside = [];
source5.mom     = cell(1,npnt);
inside = source5.inside;
for k = 1:numel(inside)
  source5.mom{inside(k)} = randn(20,3,100);
end
source5.momdimord = '{pos}_rpt_dir_time';
source5.pow       = zeros(npnt,20,100);
for k = 1:numel(inside)
  source5.pow(k,:,:) = squeeze(sum(abs(source5.mom{inside(k)}).^2,2));
end
source5.powdimord = 'pos_rpt_time';

% cortical sheet multiple trials with fixed orientation
source6 = [];
source6.pos = pnt;
source6.tri = tri;
source6.ori = nan(npos,3); %FIXME or should this be 3xnpos?
source6.ori(inside,:) = randn(numel(inside),3);
for k = 1:numel(inside)
  source6.ori(inside(k),:) = source6.ori(inside(k),:)./norm(source6.ori(inside(k),:));
end
source6.time = (-10:89)./100;
source6.inside  = 1:npnt;
source6.outside = [];
source6.mom     = cell(1,npnt);
inside = source6.inside;
for k = 1:numel(inside)
  source6.mom{inside(k)} = randn(20,100);
end
source6.momdimord = '{pos}_rpt_time';
source6.pow       = zeros(npnt,20,100);
for k = 1:numel(inside)
  source6.pow(k,:,:) = abs(source6.mom{inside(k)}).^2;
end
source6.powdimord = 'pos_rpt_time';

% 3D grid average freq domain
source7      = [];
source7.pos  = pos;
source7.dim  = dim;
source7.freq = 10;
source7.inside  = inside(:)';
source7.outside = outside(:)';
source7.csd     = cell(1,npos);
for k = 1:numel(inside)
  tmp = randn(3,1)+randn(3,1).*1i;
  source7.csd{inside(k)} = tmp*tmp';
end
source7.csddimord = '{pos}_dir_dir_freq';
source7.pow       = zeros(npos,1);
for k = 1:numel(inside)
  source7.pow(inside(k)) = trace(source7.csd{inside(k)});
end
source7.powdimord = 'pos_freq';

author: fieldtrip

repo: fieldtrip


text:
%  Figure 7.64      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 7. 64
% Robust Servo Example
%% Example 7.36; Figure 7.64
clf;
clear all;
F=[0 1; 0 -1];
G=[0;1];
H=[1 0];
J=[0];
omega=1;
a=[0 1 0 0;-omega*omega 0 1 0;0 0 0 1;0 0 0 -1];
b=[0;0;G];

% desired closed-loop poles
pc=[-1+sqrt(3)*j;-1-sqrt(3)*j;-sqrt(3)+j;-sqrt(3)-j];
k=acker(a,b,pc);

% form controller matrices
k1=k(:,1:2);
ko=k(:,3:4);
ac=[0 1;-omega*omega 0];
bc=-[k(2);k(1)];;
cc=[1 0];
dc=[0];

% closed-loop system
acl=[F-G*ko G*cc;bc*H ac];
bcl=[0;0;-bc];
ccl=[H 0 0];
dcl=[0];
syscl=ss(acl,bcl,ccl,dcl);
pole(syscl)
tzero(syscl)
% blocking zeros
dcle=-1;
syse=ss(acl,bcl,ccl,dcle);
tzero(syse)
% Closed-loop response
t=0:.01:25;
r=sin(t);
y=lsim(syscl,r,t);
figure(1)
plot(t,y,t,r);
text(0.5,0.8,'r');
text(1.5,.4,'y');
xlabel('Time (sec)');
ylabel('Reference, output');
title('Fig. 7.64 (a): Tracking response');
nicegrid;

%Control effort
cclu=[-ko cc];
dclu=[0];
sysclu=ss(acl,bcl,cclu,dclu);
u=lsim(sysclu,r,t);
figure(2)
plot(t,u);
xlabel('Time (sec)');
ylabel('Control, u');
title('Fig. 7.64 (b): Control effort');
nicegrid

%error signal
ccle=[-H 0 0];
dcle=[1];
syscle=ss(acl,bcl,ccle,dcle);
e=lsim(syscle,r,t);
figure(3)
plot(t,e);
xlabel('Time (sec)');
ylabel('Error signal, e');
title('Fig. 7.64 (c): Tracking error');
nicegrid;

author: Sable

repo: mcbench-benchmarks


text:
function [C, sigma] = dataset3Params(X, y, Xval, yval)
%EX6PARAMS returns your choice of C and sigma for Part 3 of the exercise
%where you select the optimal (C, sigma) learning parameters to use for SVM
%with RBF kernel
%   [C, sigma] = EX6PARAMS(X, y, Xval, yval) returns your choice of C and
%   sigma. You should complete this function to return the optimal C and
%   sigma based on a cross-validation set.
%

% You need to return the following variables correctly.
C = 1;
sigma = 0.3;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the optimal C and sigma
%               learning parameters found using the cross validation set.
%               You can use svmPredict to predict the labels on the cross
%               validation set. For example,
%                   predictions = svmPredict(model, Xval);
%               will return the predictions on the cross validation set.
%
%  Note: You can compute the prediction error using
%        mean(double(predictions ~= yval))
%

Carr = [0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30];
sigmaArr = [0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30];

% bestC = 0;
% bestSigma = 0;

model= svmTrain(X, y, C, @(x1, x2) gaussianKernel(x1, x2, sigma));
predictions = svmPredict(model, Xval);
minErr = mean(double(predictions~=yval));

for c = 1:length(Carr)
    for s = 1:length(sigmaArr)
        model= svmTrain(X, y, Carr(c), @(x1, x2) gaussianKernel(x1, x2, sigmaArr(s)));
        predictions = svmPredict(model, Xval);
        err = mean(double(predictions~=yval));
        if(err<=minErr)
            C = Carr(c);
            sigma = sigmaArr(s);
            minErr = err;
        end;
    end;
end;

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

end

repo: machine-learning-coursera


text:
function [BE,I]=vms_sp_read(filename,pl);
% function [BE,I]=vms_sp_read(filename,pl);
%
% function for reading Kratos spectral VMS files into two vectors
% BE - binding energy scale
% I - intensity values
% pl=1 - the spectra will be plotted
%
% example:
% [X1,C1]=vms_sp_read('s_C_1s_1.vms',0);
%
% Written by K. Artyushkova

% Kateryna Artyushkova
% Postdoctoral Scientist
% Department of Chemical and Nuclear Engineering
% The University of New Mexico
% (505) 277-0750
% kartyush@unm.edu
%
% 11/13/2003

n = textread(filename,'%f',1,'delimiter','.','headerlines',60); % number of points in spectrum
I = textread(filename,'%f',n,'delimiter','.','headerlines',63); % vector of intensity values

PhE = textread(filename,'%f',1,'delimiter','.','headerlines',27); % Photon Kinetic energy
KEstart = textread(filename,'%f',1,'delimiter','.','headerlines',47); % start Kineting energy
BEstart=PhE-KEstart; % start Binding energy

step = textread(filename,'%f',1,'delimiter','.','headerlines',48); %scanning step in eV
BEend=BEstart-n*step+step; % end Binding energy
BE=[BEstart:-step:BEend]; % vector of Binding energy
BE=BE';

%displaying image
if pl==1
plot(BE,I)
set(gca,'Xdir','reverse')
end


text:
function evalRP_group(expidxs,bUseHeadSize,bSave,bHeadSizeFromRect)

fprintf('evalRP()\n');

if (nargin < 2)
    bUseHeadSize = true;
end

if (nargin < 3)
    bSave = true;
end

if (nargin < 4)
    bHeadSizeFromRect = true;
end

if (bUseHeadSize)
    prefix = 'pck-head-rpc-';
    range = 0:0.05:0.5;
else
    prefix = 'pck-torso-rpc-';
    range = 0:0.02:0.2;
end

fontSize = 18;
figure(100); clf; hold on;
legendName = cell(0);
set(0,'DefaultAxesFontSize', fontSize)
set(0,'DefaultTextFontSize', fontSize)

[~,parts] = util_get_parts24();
nJoints = 14;
jidxsUpperBody = 7:12;

for expidx = expidxs
    fprintf('**********************************************************************************\n');
    % load experiment parameters
    p = rcnn_exp_params(expidx);

    if (isfield(p,'nms_thresh_det'))
        nms_thresh_det = p.nms_thresh_det;
    else
        nms_thresh_det = 0.5;
    end
%     % load ground truth
%     [annolist, imgidxs_gt, rectidxs_gt, rect_ignore, ~, keypointsidxs] = getAnnolist(expidx);
%     annolist_gt = annolist(imgidxs_gt);

    fname = '/home/andriluk/IMAGES/human_pose_dataset/dataset/dataset_release_candidate1/test/multPerson/h400/annolist-2-people';
    try
        load(fname,'annolist','imgidxs_gt','rectidxs_gt','rect_ignore','keypointsidxs');
        annolist_gt = annolist(imgidxs_gt);
    catch
        [annolist, imgidxs_gt, rectidxs_gt, rect_ignore, ~, keypointsidxs] = getAnnolist(expidx);
        annolist_gt = annolist(imgidxs_gt);
        save(fname,'annolist','imgidxs_gt','rectidxs_gt','rect_ignore','keypointsidxs')
    end

    mkdir_if_missing(p.latexDir);
    mkdir_if_missing(p.plotsDir);

    assert(length(p.pidxs) == nJoints);

    % load keypoints
    if (exist([p.evalTest '.mat'],'file')>0)
        fprintf('load %s\n',p.evalTest);
        load([p.evalTest '.mat']);
    else
        warning('file not found: %s\n',[p.evalTest '.mat']);
        keypointsAll = backprojectKeypointsImg(expidx,1,length(annolist_gt),false);
    end

    if (isfield(p,'bp'))
        bp = p.bp;
    else
        bp = false;
    end

    if (~bp)
        keypointsAll = keypointsAll(keypointsidxs);
    end
    assert(length(annolist_gt) == length(keypointsAll));

    if (isfield(p,'img_2_people_only'))
        nrects = cellfun(@length, rectidxs_gt(imgidxs_gt));
        annolist_gt = annolist_gt(nrects == 2);
        imgidxs_gt = imgidxs_gt(nrects == 2);
        keypointsidxs = keypointsidxs(nrects == 2);
        keypointsAll = keypointsAll(nrects == 2);
    end

    % get detections from neck detector
    if (isfield(p,'neckDetFilename'))
        fnameNeck = p.neckDetFilename;
    else
        fnameNeck = [fileparts(p.evalTest) '/boxesNeck'];
    end

    if (isfield(p,'oracle_det'))
        oracle_det = p.oracle_det;
    else
        oracle_det = false;
    end

    if (~bp)
        if (oracle_det)
            boxesNeck = struct();
            for imgidx = 1:length(annolist_gt)
                rect = annolist_gt(imgidx).annorect;
                boxes = [];
                for ridx = 1:length(rect)
                    if (isfield(rect(ridx), 'annopoints') && isfield(rect(ridx).annopoints, 'point'))
                        points = rect(ridx).annopoints.point;
                        x = [points.x];
                        y = [points.y];
                        x1 = min(x);
                        x2 = max(x);
                        y1 = min(y);
                        y2 = max(y);
                        boxes = [boxes; x1 y1 x2 y2];
                    end
                end
                boxesNeck(imgidx).det = boxes;
            end
        else
            try
                %         assert(false);
                load(fnameNeck,'boxesNeck');
            catch
                fboxes = [fileparts(p.evalTest) '/bboxAll'];
                load(fboxes,'bboxAll');
                refHeadSizeAll = zeros(length(keypointsAll),1);
                nNecks = zeros(length(keypointsAll),1);
                for imgidx = 1:length(keypointsAll)
                    keypointsNeck(imgidx).imgname = keypointsAll(imgidx).imgname;
                    rect = annolist_gt(imgidx).annorect;
                    headSize = [];
                    for ridx = 1:length(rect)
                        headSize = [headSize; util_get_head_size(rect(ridx))];
                    end
                    refDist = mean(headSize);
                    keep = nms_dist(keypointsAll(imgidx).det{9},nms_thresh_det,refDist);
                    keypointsNeck(imgidx).det = {keypointsAll(imgidx).det{9}(keep,:)};
                    refHeadSizeAll(imgidx) = refDist;
                    nNecks(imgidx) = size(keypointsNeck(imgidx).det{1},1);
                    %         keypointsAll(imgidx).det(9) = keypointsNeck(imgidx).det;
                end

                boxesNeck = struct();
                for imgidx = 1:length(keypointsAll)
                    fprintf('.');
                    img = imread(keypointsNeck(imgidx).imgname);
                    [Y,X,~] = size(img);
                    boxesNeck(imgidx).imgname = keypointsNeck(imgidx).imgname;
                    %         [val,idxs] = sort(keypointsNeck(imgidx).det{1}(:,end),'descend');
                    points = keypointsNeck(imgidx).det{1};
                    boxes = [points(:,1)-2*refHeadSizeAll(imgidx) points(:,2)-1.5*refHeadSizeAll(imgidx) ...
                        points(:,1)+2*refHeadSizeAll(imgidx) points(:,2)+7*refHeadSizeAll(imgidx)];
                    boxes(:,1) = max(1,boxes(:,1));
                    boxes(:,2) = max(1,boxes(:,2));
                    boxes(:,3) = min(X,boxes(:,3));
                    boxes(:,4) = min(Y,boxes(:,4));
                    boxesNeck(imgidx).det = boxes;
                    figure(200); clf; imagesc(imread(annolist_gt(imgidx).image.name));
                    hold on; axis equal;
                    for i=1:size(boxes,1)
                        plot(points(i,1),points(i,2),'r+');
                        rectangle('Pos',[boxes(i,1) boxes(i,2) boxes(i,3) - boxes(i,1) boxes(i,4) - boxes(i,2)],'lineWidth',5,'edgeColor','g');
                    end
                    if (~mod(imgidx, 100))
                        fprintf(' %d/%d\n',imgidx,length(keypointsAll));
                    end
                end
                save(fnameNeck,'boxesNeck');
            end
        end
    end

    if (isfield(p,'img_2_people_only'))
        for imgidx = 1:length(boxesNeck)
        end
    end

    tp_score_all = cell(length(p.pidxs),1);
    fp_score_all = cell(length(p.pidxs),1);
    totalPos_all = cell(length(p.pidxs),1);
    tp_points_all = cell(length(p.pidxs),1);
    fp_points_all = cell(length(p.pidxs),1);
    tp_idxs_all = cell(length(p.pidxs),1);
    fp_idxs_all = cell(length(p.pidxs),1);
    fp_dist_all = cell(length(p.pidxs),1);
    missing_recall_all = cell(length(p.pidxs),1);

    fprintf('compute tp/fp\n');
    if (bp)
        if (isfield(p,'marg_scores'))
            marg_scores = p.marg_scores;
        else
            marg_scores = false;
        end

        [distAll,scoresAll,pointsAll] = assignGTinference(keypointsAll,annolist_gt,p.pidxs,parts,0.5,marg_scores);
        for i = 1:length(p.pidxs)
            pidx = p.pidxs(i);
            fprintf('pidx %d\n',pidx);
            [tp_score_all{i},fp_score_all{i},totalPos_all{i},tp_points_all{i},fp_points_all{i},tp_idxs_all{i},fp_idxs_all{i},fp_dist_all{i},missing_recall_all{i}] = ...
                computeTFP(distAll{i},scoresAll{i},pointsAll{i},rectidxs_gt(imgidxs_gt),rect_ignore(imgidxs_gt));
        end
    else
        for i = 1:length(p.pidxs)
            pidx = p.pidxs(i);
            fprintf('pidx %d\n',pidx);
            % part is a joint
            assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
            jidx = parts(pidx+1).pos(1);

            for imgidx = 1:length(keypointsAll)
                [val,idxs] = sort(keypointsAll(imgidx).det{jidx+1}(:,end),'descend');
                keypointsAll(imgidx).det{jidx+1} = keypointsAll(imgidx).det{jidx+1}(idxs(1:p.top_det),:);
            end

            [dist,score,points] = getNormGTJointDistRP_group(keypointsAll,annolist_gt,jidx,boxesNeck);

            [tp_score_all{i},fp_score_all{i},totalPos_all{i},tp_points_all{i},fp_points_all{i},tp_idxs_all{i},fp_idxs_all{i},fp_dist_all{i},missing_recall_all{i}] = ...
                computeTFP(dist,score,points,rectidxs_gt(imgidxs_gt),rect_ignore(imgidxs_gt));%,rect_ignore(imgidxs_gt));
        end
    end
%         pidxs = p.pidxs;
%         save(fnameDist,'tp_score_all','fp_score_all','totalPos_all','pidxs','tp_points_all','fp_points_all','tp_idxs_all','fp_idxs_all');

%     groups = {22, 21, 12, 19, 14, 17, 16, 9, 0, 7, 2, 5, 4, 23};%, [0 9], [2 7], [4 5], [12 21], [14 19], [16 17], [22 23]
    groups = {0,2,4,5,7,9,12,14,16,17,19,21,22,23};
%     labels = {'neck','lwrist','rwrist','lelbow','relbow','lshoulder','rshoulder','lankle','rankle','lknee','rknee','lhip','rhip','tophead'};
%     labels = {'neck','lwrist','rwrist','lelbow','relbow','lshoulder','rshoulder','lankle','rankle','lknee','rknee','lhip','rhip','tophead'};
    labels = {'rankle','rknee','rhip','lhip','lknee','lankle','rwrist','relbow','rshoulder','lshoulder','lelbow','lwrist','neck','tophead'};
    apAll = zeros(length(groups),1);
    figure(100);
    for i = 1:length(groups)
        clf;
        if (length(groups{i}) == 2)
            tp_score = [tp_score_all{p.pidxs == groups{i}(1)};tp_score_all{p.pidxs == groups{i}(2)}];
            fp_score = [fp_score_all{p.pidxs == groups{i}(1)};fp_score_all{p.pidxs == groups{i}(2)}];
            totalPos = totalPos_all{p.pidxs == groups{i}(1)} + totalPos_all{p.pidxs == groups{i}(2)};
        else
            tp_score = tp_score_all{p.pidxs == groups{i}(1)};
            fp_score = fp_score_all{p.pidxs == groups{i}(1)};
            totalPos = totalPos_all{p.pidxs == groups{i}(1)};
        end

        [precision,recall,sorted_scores,sortidx,sorted_labels] = getRPC([tp_score; fp_score], [ones(length(tp_score),1);zeros(length(fp_score),1)],totalPos);
        plotRPC(precision,recall,p.colorName,'-',labels{i});
        ap = VOCap(recall,precision)*100;
        apAll(i) = ap;
        legendName = sprintf('%s, AP: %1.1f%%', p.name, ap);
        legend(legendName,'Location','NorthWest');
        fprintf('%s\n',legendName);

        if (bSave)
            print(gcf, '-dpng', [p.plotsDir '/' prefix labels{i} '-' 'expidx' num2str(expidx) '.png']);
            printpdf([p.plotsDir '/' prefix labels{i} '-' 'expidx' num2str(expidx) '.pdf']);
        end

        if (length(groups{i}) == 1)
            fp_idxs_sort = sortidx(sorted_labels == 0) - length(tp_score);
            first_idxs = fp_idxs_all{p.pidxs == groups{i}(1)}(fp_idxs_sort);
            first_scores = fp_score_all{p.pidxs == groups{i}(1)}(fp_idxs_sort);
            first_points = fp_points_all{p.pidxs == groups{i}(1)}(fp_idxs_sort,:);
            first_dist = fp_dist_all{p.pidxs == groups{i}(1)}(fp_idxs_sort);
            %         tp_idxs_sort = sortidx(sorted_labels == 1);
            %         first_idxs = tp_idxs_all{p.pidxs == groups{i}(1)}(tp_idxs_sort);
            %         first_scores = tp_score_all{p.pidxs == groups{i}(1)}(tp_idxs_sort);
            %         first_points = tp_points_all{p.pidxs == groups{i}(1)}(tp_idxs_sort,:);

%             visTFP(first_idxs(1:50),first_points,first_dist,first_scores,annolist_gt,rectidxs_gt(imgidxs_gt),parts(groups{i}(1)+1).pos(1),[fileparts(p.evalTest) '/vis/fp/pidx_' num2str(groups{i}(1))]);
            missing_recall = missing_recall_all{p.pidxs == groups{i}(1)};
            idxs = find(cellfun(@isempty,missing_recall) == 0);

            fprintf('missing recall: %1.2f%%\n',sum(cellfun(@length,missing_recall))/totalPos_all{p.pidxs == groups{i}(1)}*100);
%             visMissingRecall(missing_recall(idxs(1:50)),annolist_gt(idxs(1:50)),parts(groups{i}(1)+1).pos(1),[fileparts(p.evalTest) '/vis/mr/pidx_' num2str(groups{i}(1))]);
        end
    end
    if (bSave)
        tableFilename = [p.latexDir '/ap-expidx' num2str(expidx) '.tex'];
        apAll(15,:) = mean(apAll(7:14));
        apAll(16,:) = mean(apAll);
        [row, header] = genTable(apAll,p.name);
        fid = fopen(tableFilename,'wt');assert(fid ~= -1);
        fprintf(fid,'%s\n',row{1});fclose(fid);
        fid = fopen([p.latexDir '/ap-header.tex'],'wt');assert(fid ~= -1);
        fprintf(fid,'%s\n',header);fclose(fid);
    end
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function rescale = spm_cfg_eeg_tf_rescale
% configuration file for rescaling spectrograms
%__________________________________________________________________________
% Copyright (C) 2009-2013 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_cfg_eeg_tf_rescale.m 6825 2016-07-04 10:03:57Z vladimir $

%--------------------------------------------------------------------------
% D
%--------------------------------------------------------------------------
D        = cfg_files;
D.tag    = 'D';
D.name   = 'File Name';
D.filter = 'mat';
D.num    = [1 1];
D.help   = {'Select the M/EEG mat file.'};

%--------------------------------------------------------------------------
% Db
%--------------------------------------------------------------------------
Db        = cfg_files;
Db.tag    = 'Db';
Db.name   = 'External baseline dataset';
Db.filter = 'mat';
Db.num    = [0 1];
Db.val    = {[]};
Db.help   = {'Select the baseline M/EEG mat file. Leave empty to use the input dataset'};

%--------------------------------------------------------------------------
% timewin
%--------------------------------------------------------------------------
timewin         = cfg_entry;
timewin.tag     = 'timewin';
timewin.name    = 'Baseline time window';
timewin.help    = {'Start and stop of the baseline time window [ms].'};
timewin.strtype = 'r';
timewin.num     = [1 2];
timewin.val     = {[-Inf 0]};

%--------------------------------------------------------------------------
% pooledbaseline
%--------------------------------------------------------------------------

pooledbaseline = cfg_menu;
pooledbaseline.tag = 'pooledbaseline';
pooledbaseline.name = 'Pool baseline across trials';
pooledbaseline.labels = {'Yes', 'No'};
pooledbaseline.val = {0};
pooledbaseline.values = {1,0};
pooledbaseline.help = {'Combine baseline across trials to avoid bias',...
    'See Ciuparu and Muresan Eur J Neurosci. 43(7):861-9, 2016'};

%--------------------------------------------------------------------------
% baseline
%--------------------------------------------------------------------------
baseline         = cfg_branch;
baseline.tag     = 'baseline';
baseline.name    = 'Baseline';
baseline.help    = {'Baseline parameters.'};
baseline.val     = {timewin, pooledbaseline, Db};

%--------------------------------------------------------------------------
% method_logr
%--------------------------------------------------------------------------
method_logr      = cfg_branch;
method_logr.tag  = 'LogR';
method_logr.name = 'Log Ratio';
method_logr.val  = {baseline};
method_logr.help = {'Log Ratio.'};

%--------------------------------------------------------------------------
% method_diff
%--------------------------------------------------------------------------
method_diff      = cfg_branch;
method_diff.tag  = 'Diff';
method_diff.name = 'Difference';
method_diff.val  = {baseline};
method_diff.help = {'Difference.'};

%--------------------------------------------------------------------------
% method_rel
%--------------------------------------------------------------------------
method_rel       = cfg_branch;
method_rel.tag   = 'Rel';
method_rel.name  = 'Relative';
method_rel.val   = {baseline};
method_rel.help  = {'Relative.'};

%--------------------------------------------------------------------------
% method_zscore
%--------------------------------------------------------------------------
method_zscore       = cfg_branch;
method_zscore.tag   = 'Zscore';
method_zscore.name  = 'Zscore';
method_zscore.val   = {baseline};
method_zscore.help  = {'Z score'};

%--------------------------------------------------------------------------
% method_log
%--------------------------------------------------------------------------
method_log       = cfg_const;
method_log.tag   = 'Log';
method_log.name  = 'Log';
method_log.val   = {1};
method_log.help  = {'Log.'};

%--------------------------------------------------------------------------
% method_logeps
%--------------------------------------------------------------------------
method_logeps       = cfg_const;
method_logeps.tag   = 'Logeps';
method_logeps.name  = 'Log+eps';
method_logeps.val   = {1};
method_logeps.help  = {'Log + epsilon (to avoid boosting very low values)'};

%--------------------------------------------------------------------------
% method_sqrt
%--------------------------------------------------------------------------
method_sqrt      = cfg_const;
method_sqrt.tag  = 'Sqrt';
method_sqrt.name = 'Sqrt';
method_sqrt.val  = {1};
method_sqrt.help = {'Square Root.'};

%--------------------------------------------------------------------------
% method_sqrt
%--------------------------------------------------------------------------
method_none      = cfg_const;
method_none.tag  = 'None';
method_none.name = 'None';
method_none.val  = {1};
method_none.help = {'No rescaling - just copy. Useful for optimising batch pipelines.'};

%--------------------------------------------------------------------------
% method
%--------------------------------------------------------------------------
method        = cfg_choice;
method.tag    = 'method';
method.name   = 'Rescale method';
method.val    = {method_logr};
method.help   = {'Select the rescale method.'};
method.values = {method_logr method_diff method_rel method_zscore method_log method_logeps method_sqrt, method_none};

%--------------------------------------------------------------------------
% prefix
%--------------------------------------------------------------------------
prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the output dataset. Default prefix is ''r''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.val     = {'r'};

%--------------------------------------------------------------------------
% rescale
%--------------------------------------------------------------------------
rescale          = cfg_exbranch;
rescale.tag      = 'rescale';
rescale.name     = 'Time-frequency rescale';
rescale.val      = {D, method, prefix};
rescale.help     = {'Rescale (avg) spectrogram with nonlinear and/or difference operator.'
              'For ''Log'' and ''Sqrt'', these functions are applied to spectrogram.'
              'For ''LogR'', ''Rel'' and ''Diff'' this function computes power in the baseline.'
              'p_b and outputs:'
              '(i) p-p_b for ''Diff'''
              '(ii) 100*(p-p_b)/p_b for ''Rel'''
              '(iii) log (p/p_b) for ''LogR'''}';
rescale.prog     = @eeg_tf_rescale;
rescale.vout     = @vout_eeg_tf_rescale;
rescale.modality = {'EEG'};

%==========================================================================
function out = eeg_tf_rescale(job)
% construct the S struct
S.D         = job.D{1};
S.method    = char(fieldnames(job.method));
S.prefix    = job.prefix;

if ismember(lower(S.method), {'logr','diff', 'rel', 'zscore'})
    S.timewin        = job.method.(S.method).baseline.timewin;
    S.pooledbaseline = job.method.(S.method).baseline.pooledbaseline;
    if ~(isempty(job.method.(S.method).baseline.Db) || isequal(job.method.(S.method).baseline.Db, {''}))
        S.Db = job.method.(S.method).baseline.Db{1};
    end
end

out.D          = spm_eeg_tf_rescale(S);
out.Dfname     = {fullfile(out.D)};

%==========================================================================
function dep = vout_eeg_tf_rescale(job)
% return dependencies
dep(1)            = cfg_dep;
dep(1).sname      = 'Rescaled TF Data';
dep(1).src_output = substruct('.','D');
dep(1).tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2)            = cfg_dep;
dep(2).sname      = 'Rescaled TF Datafile';
dep(2).src_output = substruct('.','Dfname');
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm


text:
function srMethod = SRMethodConfig(srMethodName, type, estimateMotion, estimateMotionBackward, pixelSampling)

    if nargin < 5
        pixelSampling = 'center';
    end

    % Name of the super-resolution method in plain text.
    srMethod.name = srMethodName;

    srMethod.type = type;

    % Flag to enable motion estimation.
    % 'false': No motion estimation
    % 'true': Perform motion estimation (for multi-frame algorithms only)
    srMethod.estimateMotion = estimateMotion;

    % Flag to enable motion estimation in backward direction
    % 'false': No motion estimation
    % 'true: Perform motion estimation in backward direction (required for
    %        some multi-frame algorithms only)
    srMethod.estimateMotionBackward = estimateMotionBackward;

	% Pixel sampling scheme. 'center' (default) samples pixels at the center
	% point, whereas 'topLeft' samples pixels at the top-left corner.
	srMethod.pixelSampling = pixelSampling;

author: thomas-koehler

repo: SupER


text:
function test_bug2364

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY

nTrials = 20;
data = [];
data.fsample = 256;
data.trial = arrayfun(@(x) rand(10, 512), 1:nTrials, 'UniformOutput', false);
data.time = repmat({(0:length(data.trial{1})-1)/data.fsample}, [1, nTrials]);
data.label = cellfun(@num2str, num2cell(1:10), 'UniformOutput', false);

for iTrial = 1:nTrials
  data.trial{iTrial}(1,:) = nan;
end

%%

cfg = [];
cfg.method = 'wavelet';
cfg.width = 2;
cfg.foilim = [10 50];
cfg.toi = 0:0.01:2;

freq1 = ft_freqanalysis(cfg, data);
count1 = sum(~isnan(freq1.powspctrm(:)));
fprintf('number of non-nans for wavelet method: %i\n', count1)

%%

cfg = [];
cfg.method = 'mtmconvol';
cfg.taper = 'hanning';
cfg.foi = 10:50;
cfg.t_ftimwin = 2./cfg.foi;
cfg.toi = 0:0.01:2;

freq2 = ft_freqanalysis(cfg, data);
count2 = sum(~isnan(freq2.powspctrm(:)));
fprintf('number of non-nans for mtmconvol method: %i\n', count2)

%%

cfg = [];
cfg.method = 'mtmfft';
cfg.taper = 'hanning';

freq3 = ft_freqanalysis(cfg, data);
count3 = sum(~isnan(freq3.powspctrm(:)));
fprintf('number of non-nans for mtmfft method: %i\n', count3)

%%

assert(count1==140067)
assert(count2==70803)
assert(count3==2313);

%% deal with regression problem in single-channel data

data1 = [];
data1.label = {'1'};
data1.time = {(1:1000)/1000};
data1.trial = {randn(1,1000)};

data2 = [];
data2.label = {'1', '2'};
data2.time = {(1:1000)/1000};
data2.trial = {randn(2,1000)};
data2.trial{1}(1,:) = data1.trial{1};

cfg = [];
cfg.method = 'wavelet';
cfg.width = 2;
cfg.foilim = [10 50];
cfg.toi = 0:0.01:1;
tf1 = ft_freqanalysis(cfg, data1);
tf2 = ft_freqanalysis(cfg, data2);

figure; imagesc(abs(squeeze(tf1.powspctrm(1,:,:))))
figure; imagesc(abs(squeeze(tf2.powspctrm(1,:,:))))

difference = abs(tf1.powspctrm(1,:,:) - tf2.powspctrm(1,:,:));
assert(all(difference(~isnan(difference(:)))<100*eps)); % there is a small numerical difference