text:
function stats = calcStatistics(optObj,limit,smoothBounds)
%Calculate confidence and variance statistics of an OPTI curve fitting problem
%
%   Called By opti fitStats

%   Copyright (C) 2011-2014 Jonathan Currie (IPL)

% This function uses ideas from:
% - David M. Himmelblau. Process Analysis by Statistical Methods.
%   John Wiley & Sons, 1970
% - David I. Wilson. Advanced Control Using MATLAB, AUT University, 2014.
% - NonLinearModel.m from the Statistics Toolbox
% - confint.m from the Curve Fitting Toolbox
% - predint.m from the Curve Fitting Toolbox

% Value explanations
% http://www.ats.ucla.edu/stat/sas/output/reg.htm

if(nargin < 3 || isempty(smoothBounds)), smoothBounds = true; end
if(nargin < 2 || isempty(limit)), limit = 0.95; end

if(limit <= 0 || limit >= 1)
    error('Confidence limit must be 0 < limit < 1');
end

%Get Problem
prob = optObj.prob;
x = optObj.sol;
nparam = length(x);
ndata = length(prob.ydata); %check whether we need original data?
wlevel = optiWarnLevel(optObj.opts.warnings);
%Get Solution SSE
sse = optObj.obj;
dfe = ndata-nparam;
rmse = sqrt(sse/dfe);
if(wlevel>1)
    optiwarn('optistats:dfe','Your problem has #data <= #parameters, confidence limits and most statistics cannot be estimated');
end
%Build Return Structure
stats = struct('SSE',sse,'Rsquare',NaN,'AdjRsquare',NaN,'RMSE',NaN,'DFE',ndata-nparam,'ConfInt',NaN,'ConfBnds',struct('xdata',NaN,'bnds',NaN),'Cov',NaN,'BndIdx',NaN,...
               'Param',struct('StdError',NaN,'tStat',NaN,'pValues',NaN),'Model',struct('FStat',NaN,'pValue',NaN),'Conf',limit,'SolverStatus','');

nidx = strfind(prob.Name,'optifit:');
if(~isempty(nidx))
    stats.ModelStructure = prob.Name(9:end);
end

%Don't be overly enthusiastic
if(~isempty(strfind(optObj.info.Status,'Optimal')))
    stats.SolverStatus = 'Converged';
else
    stats.SolverStatus = optObj.info.Status;
end

%Ensure we have a NLS or DNLS
if(~any(strcmpi(prob.type,{'NLS','DNLS'})))
    error('This function can only be used on curve fitting (NLS,DNLS) problems');
end

%Get ydata and weights
ydata = prob.ydata;
weights = prob.weighting;

%Get Partial Derivative Matrix (data/parameters) [original xdata]
if(~isempty(prob.xdata))
    try
        X = prob.misc.fitGrad(x,prob.xdata);
        gmode = 2;
    catch
        X = prob.misc.fitGrad(x);
        gmode = 1;
    end
else
    try
        X = prob.misc.fitGrad(x);
    catch
        error('Unknown gradient field!');
    end
    gmode = 1;
end
%If we have weights, weight the Jacobian
Xnowts = X;
if(~isempty(weights))
    X = bsxfun(@times, weights, X);
end
%Check for active bounds, if so, remove cols of X
bidx = false(size(X,2),1);
if(~isempty(prob.lb) || ~isempty(prob.ub))
    if(any(~isinf(prob.lb)))
        bidx = bidx | abs(x-prob.lb) < sqrt(eps);
    end
    if(any(~isinf(prob.ub)))
        bidx = bidx | abs(x-prob.ub) < sqrt(eps);
    end
end
if(any(bidx))
    if(wlevel > 1)
        optiwarn('opti:conf','One or more bounds are active at the solution. The confidence interval and parameter statistics cannot be calculated for these variables.');
    end
    X = X(:,~bidx);
    %Correct dfe, rmse
    nparam      = nparam - sum(bidx);
    dfe         = ndata-nparam;
    rmse        = sqrt(sse/dfe);
    stats.DFE   = dfe;
    stats.RMSE  = rmse;
end

try
    %Solve inverse of X'*X, but only interested in diagonal
    [~,R] = qr(X,0);
    if(size(R,1)~=size(R,2))
        throw(MException('opti:conf','QR R not square'));
    end
    s = warning('off','MATLAB:singularMatrix');
    s1 = warning('off','MATLAB:nearlySingularMatrix');
    Rinv = R \ eye(length(R));
    warning(s1);
    warning(s);
catch
    error('The problem is poorly scaled and resulted in a (near) singular matrix, limits cannot be calculated');
end

%Solve Covariance Matrix
s = warning('off','MATLAB:singularMatrix');
s1 = warning('off','MATLAB:nearlySingularMatrix');
try
    L = chol(X'*X,'lower');
    stats.Cov = L'\(L\eye(length(L)))*rmse^2;
catch
    if(wlevel)
        optiwarn('opti:conf','The X''*X matrix has been found to be singular, and a general inverse will be used to (attempt to) solve the system covariance.\nResults are now suspicious at best...');
    end
    stats.Cov = inv(X'*X)*rmse^2; %#ok<MINV> %poor method - any better ideas?
end
warning(s1);
warning(s);
stats.BndIdx = bidx;

%Generate ypred for ANOVA stuff
if(nargin(prob.misc.fitFun)==2)
    sumypred = sum(prob.misc.fitFun(x,prob.xdata).^2);
else
    sumypred = sum(prob.misc.fitFun(x).^2);
end

%Solve Confidence Interval
v = sum(Rinv.^2,2) * (sse / dfe);
alpha = (1-limit)/2;
ConfInt = (-rmathlib('qt',alpha,dfe)* sqrt(v'))';
stats.ConfInt = NaN(length(bidx),1);
stats.ConfInt(~bidx) = ConfInt;

% Return X to the unweighted version
% This seems to be the accepted convention for bounds calculation
X = Xnowts;

%Solve Confidence Bounds for Plotting (note we must have system covariance for this to work!)
if(~isempty(stats.Cov))
    %If we have an ODE and simple time step (i.e. not a cell array or repeated measurements, try for a smooth curve)
    if(isfield(prob.misc,'xdata_orig') && ~isempty(prob.misc.xdata_orig) && ~iscell(prob.misc.xdata_orig))
        havRep = length(unique(prob.misc.xdata_orig)) ~= length(prob.misc.xdata_orig);
    else
        havRep = false;
    end
    % Try generate smooth bounds
    if (smoothBounds == true)
        if(~isempty(prob.ode) && ~iscell(prob.misc.xdata_orig) && length(unique(prob.xdata))==length(prob.xdata) && ~havRep)
            try
                %Lazy way to convert problem again
                xd = linspace(min(prob.xdata),max(prob.xdata),max(10*length(prob.xdata),1e2))';
                prob.ydata = ones(length(prob.ydata)/length(prob.xdata)*length(xd),1); prob.xdata = xd;
                prob = DNLS2NLS(prob,optObj.opts);
                Xb = prob.misc.fitGrad(x);
                ypred = prob.misc.fitFun(x,xd);
            catch %no luck
                xd = prob.xdata;
                Xb = X;
                if(nargin(prob.misc.fitFun)==2)
                    ypred = prob.misc.fitFun(x,xd);
                else
                    ypred = prob.misc.fitFun(x);
                end
            end
        else %algebraic system or complicated ode problem
            if(gmode==2)
                %try evaluate at many intermediate points for a smoother curve
                try
                    %Assume if xdata is a matrix, or not sorted, don't try smooth
                    if(size(prob.xdata,1) > 1 && size(prob.xdata,2) > 1 || (any(sort(prob.xdata) ~= prob.xdata) && any(sort(prob.xdata,'descend') ~= prob.xdata)))
                        error('skip');
                    end
                    xd = linspace(min(prob.xdata),max(prob.xdata),max(10*length(prob.xdata),1e2))';
                    Xb = prob.misc.fitGrad(x,xd);
                    ypred = prob.misc.fitFun(x,xd);
                catch
                    xd = prob.xdata;
                    Xb = X;
                    ypred = prob.misc.fitFun(x,prob.xdata);
                end
            else %no luck
                xd = prob.xdata;
                Xb = X;
                if(nargin(prob.misc.fitFun)==2)
                    %if function has two input arguments, try a numerical gradient
                    try
                        if(~isempty(strfind(char(prob.misc.fitFun),'odeEstim'))), error('no luck'); end %integrator most likely has set time points it is expecting
                        xd = linspace(min(prob.xdata),max(prob.xdata),max(10*length(prob.xdata),1e2))';
                        Xb = mklJac(@(x) prob.misc.fitFun(x,xd),x);
                        ypred = prob.misc.fitFun(x,xd);
                    catch
                        xd = prob.xdata;
                        Xb = X;
                        ypred = prob.misc.fitFun(x,xd);
                    end
                else
                    ypred = prob.misc.fitFun(x);
                end
            end
        end
    else
        xd = prob.xdata;
        Xb = X;
        if(nargin(prob.misc.fitFun)==2)
            ypred = prob.misc.fitFun(x,xd);
        else
            ypred = prob.misc.fitFun(x);
        end
    end

    %Ensure ypred is a column
    if(size(ypred,2) > 1), ypred = ypred'; end
    %Solve Prediction Bounds (Simultaneous Functional)
    crit = sqrt(length(x) * rmathlib('qf',limit, length(x), dfe));
    %If we have active bounds, best to calculate Jacobian at each point we want to plot, then drop columns, rather than use Covariance below
    if(any(bidx))
        try
            J = prob.misc.fitGrad(x,xd);
            E = J(:,~bidx)*Rinv;
            delta = crit * sqrt(sum(E.*E,2)) * sqrt(sse/dfe);
        catch
            delta = crit * sqrt(sum((Xb*stats.Cov) .* Xb,2));
        end
    else
        delta = crit * sqrt(sum((Xb*stats.Cov) .* Xb,2));
    end
    if(isempty(xd))
        stats.ConfBnds.xdata = (1:length(ypred))';
    else
        stats.ConfBnds.xdata = xd;
    end
    stats.ConfBnds.bnds = [ypred-delta ypred+delta];

    %Reshape based on number of curves
    if(~isempty(prob.odez0))
        if(~isempty(optObj.opts.dynamicOpts) && isfield(optObj.opts.dynamicOpts,'stateIndex') && ~isempty(optObj.opts.dynamicOpts.stateIndex))
            nc = length(optObj.opts.dynamicOpts.stateIndex);
        else
            nc = length(prob.odez0);
        end
        nd = size(ypred,1)/nc;
        if(floor(nd)==nd) %have to work out xdata associated with which point... anyone interested?
            cb = zeros(nd,nc*2);
            idx = 1;
            for i = 1:2:nc*2
                cb(:,i:i+1) = stats.ConfBnds.bnds(idx:idx+nd-1,:);
                idx = idx + nd;
            end
            stats.ConfBnds.bnds = cb;
        else
            stats.ConfBnds.bnds = NaN;
            stats.ConfBnds.xdata = NaN;
            if(wlevel)
                optiwarn('opti:confbnds','Cannot currently determine confidence prediction bounds for this type of problem, sorry!');
            end
        end
    end
    %If we had repeated points, try index out...
    if(havRep)
        %Check for repeated measurements
        try
            [~,ia] = unique(prob.misc.xdata_orig,'legacy');
        catch
            [~,ia] = unique(prob.misc.xdata_orig);
        end
        stats.ConfBnds.bnds = stats.ConfBnds.bnds(ia,:);
    end

    %Make sure in order for plotting
    [stats.ConfBnds.xdata,sbidx] = sort(stats.ConfBnds.xdata);
    stats.ConfBnds.bnds = stats.ConfBnds.bnds(sbidx,:);
end

%Solve Standard Statistics
[SST,stats.Rsquare,stats.AdjRsquare,stats.RMSE,SST0] = fitstats(ydata,weights,sse,dfe,ndata);
smodel = sumypred-sse;
stats.DF = [nparam dfe ndata]; %[model error uncorrected]
stats.SOS = [smodel sse sumypred];
stats.MS = [smodel/nparam sse/(ndata-nparam)];

%Parameter statistics require covariance too...
if(~isempty(stats.Cov))
    %Parameter Statistics
    SE = NaN(nparam,1); T = SE; P = SE;
    SE(~bidx) = sqrt(diag(stats.Cov));
    T(~bidx) = x(~bidx) ./ SE(~bidx);
    P(~bidx) = 2*rmathlib('pt',-abs(T(~bidx)),dfe);
    stats.Param.StdError = SE;
    stats.Param.tStat = T;
    stats.Param.pValues = P;
end

%Model F-Test & P Value
%Check for intercept (constant column in X)
Xmin = min(X,[],1); Xmax = max(X,[],1);
if(nparam > 1 && any(abs(Xmax-Xmin) <= sqrt(eps)))
    %Intercept Model
    DFR = nparam - 1;
    DFE = ndata - 1 - DFR;
    SSR = max(SST - sse,0);
    stats.Model.FStat = (SSR/DFR) / (sse / DFE);
    stats.Model.pValue = rmathlib('pf',1/stats.Model.FStat,DFE,DFR);
    stats.Model.Int = true;
else
    %Zero Model
    DFR = nparam;
    DFE = ndata - nparam;
    SSR = max(SST0 - sse,0);
    stats.Model.FStat = (SSR/DFR) / (sse / DFE);
	stats.Model.pValue = rmathlib('pf',1/stats.Model.FStat,DFE,DFR);
    stats.Model.Int = false;
end

function [sst,rsquare,adjrsquare,rmse,sst0] = fitstats(ydata,weights,sse,dfe,ndata)
%Calculate standard statistical measures

%Spread about mean
if(isempty(weights))
    ybar = mean(ydata);
    sst = sum((ydata - ybar).^2);
    sst0 = sum(ydata.^2);
else
    ybar = sum(ydata.*weights.^2)/sum(weights.^2);
    sst = sum(weights.^2.*(ydata - ybar).^2);
    sst0 = sum(weights.^2.*ydata.^2);
end
%R^2
if(sst~=0)
    rsquare = 1 - sse/sst;
    adjrsquare = 1 - (1-rsquare)*(ndata-1)/dfe; %1 - (sse/dfe)/(sst/ndata);
else
    rsquare = NaN;
    adjrsquare = NaN;
end
%RMSE
if(dfe > 0)
    rmse = sqrt(sse/dfe);
else
    rmse = NaN;
end


text:
%PRTOOLS_MOD Read PRTools Message Of the Day

function mod = prtools_mod

url = 'http://37steps.com/files/prtools_mod.text';

if ~usejava('jvm') & isunix
	[stat,s] = unix(['wget -q -O - ' url]);
	status = (stat == 0);
else
  if verLessThan('matlab','8.0')
    status = 0; % Cannot read if server is down
  else
    [s,status] = urlread(url,'TimeOut',5);
  end
  %[s,status] = urlread(url);
  % needs updatating. s may contain a firewall message while status = 1
end

if ~status
  mod = [];
else
  mod = s;
end

author: marianux

repo: ecg-kit

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


text:
function events = NewEvents(times,description)

%NewEvents - Create events structure.
%
%  USAGE
%
%    events = NewEvents(times,description)
%
%    times               event timestamps
%    description         (common) event description

% Copyright (C) 2008-2011 by Michaël Zugaro
%
% 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.

events.time = times;
events.description = cell(size(times));
[events.description{:}] = deal(description);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function grid = reshapeGrid(~, grid)
%RESHAPEGRID   Add the repeated endpoints to a 2D periodic grid.

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

% Note: We use periodic discretizations in 1D/2D/3D. These discretizatioins
% do not include the repeated endpoints. Before plotting data with
% SPINOPERATOR/INITIALIZEMOVIE, we add these repeated endpoints to the grid.

% Get the points:
xx = grid{1};
yy = grid{2};

% Add the endpoints:
xx = [xx, 2*xx(:,end) - xx(:,end-1)];
xx = [xx; xx(1,:)];
yy = [yy; 2*yy(end,:) - yy(end-1,:)];
yy = [yy, yy(:,1)];

% Output the new grid:
grid{1} = xx;
grid{2} = yy;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function varargout = drawPolyline3d(varargin)
%DRAWPOLYLINE3D Draw a 3D polyline specified by a list of vertex coords.
%
%   drawPolyline3d(POLY);
%   packs coordinates in a single N-by-3 array.
%
%   drawPolyline3d(PX, PY, PZ);
%   specifies coordinates in separate numeric vectors (either row or
%   columns)
%
%   drawPolyline3d(..., CLOSED);
%   Specifies if the polyline is closed or open. CLOSED can be one of:
%   - 'closed'
%   - 'open'    (the default)
%   - a boolean variable with value TRUE for closed polylines.
%
%   drawPolyline3d(..., PARAM, VALUE);
%   Specifies style options to draw the polyline, see plot for details.
%
%   H = drawPolyline3d(...);
%   also returns a handle to the list of created line objects.
%
%   Example
%     t = linspace(0, 2*pi, 100)';
%     xt = 10 * cos(t);
%     yt = 5 * sin(t);
%     zt = zeros(1,100);
%     figure; drawPolyline3d(xt, yt, zt, 'b');
%
%   See also
%   polygons3d, drawPolygon3d, fillPolygon3d
%

% ------
% Author : David Legland
% E-mail: david.legland@inra.fr
% Created: 2005-02-15
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

%% Process input arguments

% Check if axes handle is specified
if isAxisHandle(varargin{1})
    hAx = varargin{1};
    varargin(1) = [];
else
    hAx = gca;
end

% check case we want to draw several curves, stored in a cell array
var = varargin{1};
if iscell(var)
    hold on;
    h = zeros(length(var(:)), 1);
    for i = 1:length(var(:))
        h(i) = drawPolyline3d(hAx, var{i}, varargin{2:end});
    end
    if nargout > 0
        varargout = {h};
    end
    return;
end

% extract curve coordinates
if min(size(var)) == 1
    % if first argument is a vector (either row or column), then assumes
    % first argument contains x coords, second argument contains y coords
    % and third one the z coords
    px = var;
    if length(varargin) < 3
        error('geom3d:drawPolyline3d:Wrong number of arguments in drawPolyline3d');
    end
    if  isnumeric(varargin{2}) && isnumeric(varargin{3})
        py = varargin{2};
        pz = varargin{3};
        varargin(1:3) = [];
    else
        px = var(:, 1);
        py = var(:, 2);
        pz = var(:, 3);
        varargin(1) = [];
    end
else
    % all coordinates are grouped in the first argument
    px = var(:, 1);
    py = var(:, 2);
    pz = var(:, 3);
    varargin(1) = [];
end

% check if curve is closed or open (default is open)
closed = false;
if ~isempty(varargin)
    var = varargin{1};
    if islogical(var)
        % check boolean flag
        closed = var;
        varargin = varargin(2:end);

    elseif ischar(var)
        % check string indicating close or open
        if strncmpi(var, 'close', 5)
            closed = true;
            varargin = varargin(2:end);

        elseif strncmpi(var, 'open', 4)
            closed = false;
            varargin = varargin(2:end);
        end

    end
end

%% draw the curve

% for closed curve, add the first point at the end to close curve
if closed
    px = [px(:); px(1)];
    py = [py(:); py(1)];
    pz = [pz(:); pz(1)];
end

h = plot3(hAx, px, py, pz, varargin{:});

if nargout > 0
    varargout = {h};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ a_lu, info ] = r8ge_np_trf ( m, n, a )

%*****************************************************************************80
%
%% R8GE_NP_TRF computes the LU factorization of a R8GE matrix.
%
%  Discussion:
%
%    The R8GE storage format is used for a general M by N matrix.  A storage
%    space is made for each logical entry.  The two dimensional logical
%    array is mapped to a vector, in which storage is by columns.
%
%    R8GE_NP_TRF is a nonpivoting version of R8GE_TRF, and will fail if
%    a zero element is encountered along the diagonal.
%
%    The factorization has the form
%      A = L * U
%    where L is lower triangular with unit diagonal elements (lower
%    trapezoidal if N < M), and U is upper triangular (upper trapezoidal
%    if M < N).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix A.  0 <= M.
%
%    Input, integer N, the number of columns of the matrix A.  0 <= N.
%
%    Input, real A(M,N), the matrix to be factored.
%
%    Output, integer INFO.
%    = 0: successful exit
%    < 0: if INFO = -K, the K-th argument had an illegal value
%    > 0: if INFO = K, U(K,K) is exactly zero. The factorization
%         has been completed, but the factor U is exactly
%         singular, and division by zero will occur if it is used
%         to solve a system of equations.
%
%    Output, real A_LU(M,N), the factors L and U from the factorization
%    A = L*U; the unit diagonal elements of L are not stored.
%
  info = 0;
  a_lu(1:m,1:n) = a(1:m,1:n);
%
%  Test the input parameters.
%
  if ( m < 0 )
    info = -1;
    return
  elseif ( n < 0 )
    info = -2;
    return
  end

  if ( m == 0 | n == 0 )
    return
  end

  for j = 1 : min ( m, n )
%
%  Compute elements J+1:M of the J-th column.
%
    if ( a(j,j) ~= 0.0E+00 )
      a_lu(j+1:m,j) = a_lu(j+1:m,j) / a_lu(j,j);
    elseif ( info == 0 )
      info = j;
    end
%
%  Update the trailing submatrix.
%
    if ( j < min ( m, n ) )

      for ii = j+1 : m
        a_lu(ii,j+1:n) = a_lu(ii,j+1:n) - a_lu(ii,j) * a_lu(j,j+1:n);
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [file] = read_ctf_ascii(filename)

% READ_CTF_ASCII reads general data from an CTF configuration file
%
% The file should be formatted like
%    Group
%    {
%      item1 : value1a value1b value1c
%      item2 : value2a value2b value2c
%      item3 : value3a value3b value3c
%      item4 : value4a value4b value4c
%    }
%
% This fileformat structure is used in
%   params.avg
%   default.hdm
%   multiSphere.hdm
%   processing.cfg
% and maybe for other files as well.

% Copyright (C) 2003, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

fid = fopen_or_error(filename, 'r');

line = '';
while ischar(line)
  line = cleanline(fgetl(fid));
  if isempty(line) || (length(line)==1 && all(line==-1))
    continue
  end

  % the line is not empty, which means that we have encountered a chunck of information
  subline = cleanline(fgetl(fid));  % read the {
  subline = cleanline(fgetl(fid));  % read the first item
  while ~contains(subline, '}')
    if ~isempty(subline)
      [item, value] = strtok(subline, ':');
      value(1) = ' ';           % remove the :
      value  = strtrim(value);
      item   = strtrim(item);

      % the item name should be a real string, otherwise I cannot put it into the structure
      if strcmp(sprintf('%d', str2num(deblank(item))), deblank(item))
        % add something to the start of the string to distinguish it from a number
        item = ['item_' item];
      end

      % the value can be either a number or a string, and is put into the structure accordingly
      if isempty(str2num(value))
        % the value appears to be a string
        eval(sprintf('file.%s.%s = [ ''%s'' ];', line, item, value));
      else
        % the value appears to be a number or a list of numbers
        eval(sprintf('file.%s.%s = [ %s ];', line, item, value));
      end

    end
    subline = cleanline(fgetl(fid));    % read the first item
  end
end

fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function line = cleanline(line)
  if isempty(line) || (length(line)==1 && all(line==-1))
    return
  end
  comment = strfind(line, '//');
  if ~isempty(comment)
    line(min(comment):end) = ' ';
  end
  line = strtrim(line);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% emacs, this is -*-Matlab-*- mode
%
% readradfile    reads the BlueC *.rad files
%
% *.rad files contain paprameters of the radial distortion
% [K,kc] = readradfile(name)
% name ... name of the *.rad file with its full path
%
% K ... 3x3 calibration matrix
% kc ... 4x1 vector of distortion parameters
%
% $Id: readradfile.m,v 2.0 2003/06/19 12:07:16 svoboda Exp $
function [K,kc] = readradfile(name);

fid = fopen(name,'r');
if fid<0
  error(sprintf('Could not open %s. Missing rad files?',name'))
end

for i=1:3,
  for j=1:3,
	buff = fgetl(fid);
        str_end = buff(7:end);
        if str_end(end)==';'
          str_end = str_end(1:end-1);
        end
        K(i,j) = str2num(str_end);
  end
end

buff = fgetl(fid);
for i=1:4,
  buff = fgetl(fid);
  str_end = buff(7:end);
  if str_end(end)==';'
    str_end = str_end(1:end-1);
  end
  kc(i) = str2num(str_end);
end

fclose(fid);

return;

author: strawlab

repo: MultiCamSelfCal


text:
addpath([cd '\ballistics_problem_files\']);
ballist_prob;

author: Sable

repo: mcbench-benchmarks


text:
classdef (Abstract) math_model_cpf_acc < mp.math_model_cpf
%MP.MATH_MODEL_CPF_ACC  MATPOWER mathematical model for continuation power flow (CPF) problem.
%   ?
%
%   MP.MATH_MODEL_CPF_ACC ... power flow ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        %% constructor
        function obj = math_model_cpf_acc()
            obj@mp.math_model_cpf();
            obj.element_classes = { @mp.mme_bus_pf_acc, @mp.mme_gen_pf_ac, ...
                @mp.mme_branch_pf_ac, @mp.mme_load_cpf, @mp.mme_shunt_cpf };
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
classdef RWMOP27 < PROBLEM
% <multi> <real> <constrained>
% Process flow sheeting problem

%------------------------------- 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        = 3;
            obj.lower    = [0.2,-2.22554,-0.49];
            obj.upper    = [1,-1,1.49];
            obj.encoding = ones(1,obj.D);
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            x  = varargin{1};
            x1 = x(:,1); x2 = x(:,2); x3 = x(:,3);
            % Objective function
            f(:,1) = -0.7.*x3 + 0.8 + 5.*(0.5 - x1).^2;
            f(:,2) = x1 - x3;
            % Constraints
            g(:,1) = -(exp(x1 - 0.2) + x2);
            g(:,2) = x2 + 1.1.*x3 - 1;
            g(:,3) = x1 - x3 - 0.2;
            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 = [-2.4300000e-01   0.0000000e+00];
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function M=movwv2at(N,Np);
%MOVWV2AT Oscillating structure of the interferences of the WVD.
%	M=MOVWV2AT(N,Np) generates the movie frames illustrating the
%	influence of the distance between two components on the oscillating
%	structure of the interferences of the WVD.
%
%	N : number of points for the signal;
%	Np : number of snapshots (default : 9)
%	M : matrix of movie frames.
%
%	Example :
%	 M=movwv2at(128,15); movie(M,10);

%	O. Lemoine - May 1996.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if nargin<1,
 error('At least one argument required');
elseif nargin==1,
 Np=9;
end
Np=odd(Np);

M  = moviein(Np);

c  = Np/0.2;

ampl=amgauss(N/2,N/2,N/2);

for k=1:(Np+1)/2,
 sig=[ampl.*fmconst(N/2,0.25-(k-1)/c);ampl.*fmconst(N/2,0.25+(k-1)/c)];
 [tfr,t,f]=tfrwv(sig);
 Max=max(max(tfr));V=[0.1 0.3 0.5 0.7 0.9]*Max;
 contour(t,f,tfr,V);xlabel('Time'); ylabel('Frequency');
 title('Wigner-Ville distribution'); axis('xy')
 M(:,k) = getframe;
end

for k=(Np+3)/2:Np,
 M(:,k) = M(:,Np+1-k);
end

author: HeLiangHIT


text:
function eval_tracker(seqs, trackers, eval_type, name_tracker_all, tmp_mat_path, path_anno, rp_all, norm_dst)
% evaluate each tracker
num_tracker = numel(trackers);

threshold_set_overlap = 0:0.05:1;
threshold_set_error   = 0:50;
if norm_dst
    threshold_set_error = threshold_set_error / 100;
end

for i = 1:numel(seqs) % for each sequence
    s    = seqs{i};      % name of sequence

    % load GT and the absent flags
    anno        = dlmread([path_anno s '.txt']);
    absent_anno = dlmread([path_anno 'absent/' s '.txt']);

    for k = 1:num_tracker  % evaluate each tracker
        t = trackers{k};   % name of tracker

        % load tracking result
        res = dlmread([rp_all t.name '_tracking_result/' s '.txt']);
        fprintf(['evaluating ' t.name ' on ' s ' ...\n']);

        success_num_overlap = zeros(1, numel(threshold_set_overlap));
        success_num_err     = zeros(1, numel(threshold_set_error));

        if isempty(res)
            break;
        end

        [err_coverage, err_center] = calc_seq_err_robust(res, anno, absent_anno, norm_dst);

        for t_idx = 1:numel(threshold_set_overlap)
            success_num_overlap(1, t_idx) = sum(err_coverage > threshold_set_overlap(t_idx));
        end

        for t_idx = 1:length(threshold_set_error)
            success_num_err(1, t_idx) = sum(err_center <= threshold_set_error(t_idx));
        end

        len_all = size(anno, 1);  % number of frames in the sequence

        ave_success_rate_plot(k, i, :)     = success_num_overlap/(len_all + eps);
        ave_success_rate_plot_err(k, i, :) = success_num_err/(len_all + eps);
    end
end

% save results
if ~exist(tmp_mat_path, 'dir')
    mkdir(tmp_mat_path);
end

dataName1 = [tmp_mat_path 'aveSuccessRatePlot_' num2str(num_tracker) 'alg_overlap_' eval_type '.mat'];
save(dataName1, 'ave_success_rate_plot', 'name_tracker_all');

dataName2 = [tmp_mat_path 'aveSuccessRatePlot_' num2str(num_tracker) 'alg_error_' eval_type '.mat'];
ave_success_rate_plot = ave_success_rate_plot_err;
save(dataName2, 'ave_success_rate_plot', 'name_tracker_all');

end

author: HengLan


text:
%SNR
function Sn=SNR(Im_original,Im_modified)

if (size(Im_original)~=size(Im_modified))
    error ('error:image sizes do not agree')
end

else
    A=double(Im_original);
    B=double(Im_modified);
end

[m,n]=size(A);
sumaI=0;
sumaDif=0;
for u=1:m
    for v=1:n
        sumaI=sumaI+A(u,v)^2;
        sumaDif=sumaDif+(A(u,v)-B(u,v))^2;
    end
end

if (sumaDif==0)
    sumaDif=1;
end

Sn=sumaI/sumaDif;
Sn=10*log10(Sn);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
classdef ncs_wc_reader < handle
	properties
        sr
        max_segments
        raw_filename
        %channel
        opened_file
        recmin
        recmax
        num_scale_factor
        TimeStamps
        dt
    end
	methods
        function obj = ncs_wc_reader(par, raw_filename)
            obj.raw_filename = raw_filename;
            obj.opened_file = fopen(obj.raw_filename, 'r', 'l');
            fseek(obj.opened_file,16384,'bof');                       %Skip Header, put pointer to the first record
            TimeStamps = fread(obj.opened_file,inf,'int64',(4+4+4+2*512));          %Read all TimeStamps, in us
            fseek(obj.opened_file, 16384+8+4+4+4, 'bof');                             %put pointer to the beginning of data

            % conversion to datapoints
            tsdiff = diff(TimeStamps);
            obj.dt = median(tsdiff);
            min_dt = min(diff(TimeStamps));
            if min_dt<=0 %check for corrupt TimeStamps that would violate the monotonic increase
                warning('corrupt TimeStamps - attempting correction')
                ind = find(tsdiff <= 0); %corrupt TimeStamps are usually too low
                TimeStamps(ind+1) = TimeStamps(ind) + obj.dt;
                clear tsdiff;
                min_dt = min(diff(TimeStamps));
                if min_dt <= 0
                    error('Automated correction of corrupt TimeStamps failed - try manually!');
                end
            end

            %remove the following if statement to keep unix time in new
            %versions
            if TimeStamps(1) > 1e+15 %reset time to the beggining of the file
                TimeStamps = TimeStamps - TimeStamps(1);
            end
            obj.sr = 512*1e6/obj.dt;            % sampling rate (in Hz).
            time0 = TimeStamps(1);

            tmin = time0 + par.tmin*1e6; %min time to read (in micro-sec)
            index_tinitial = find(tmin > TimeStamps);
            if isempty(index_tinitial) == 1
                index_tinitial = 0;
            else
                index_tinitial = index_tinitial(end)-1;
            end

            if strcmp(par.tmax,'all')
                index_tfinal = length(TimeStamps);
                tmax = TimeStamps(end);
            else
                tmax = time0 + par.tmax*1e6;                   %max time to read (in micro-sec)
                index_tfinal = find(tmax < TimeStamps);
                if isempty(index_tfinal) ==1;
                    index_tfinal = length(TimeStamps);
                else
                    index_tfinal = index_tfinal(1);
                end
            end
            fseek(obj.opened_file,16384+8+4+4+4+index_tinitial*(8+4+4+4+2*512),'bof');

            lts = index_tfinal - index_tinitial;
            %Segments the data in par.segments pieces
            obj.max_segments = ceil((tmax - tmin)/ ...
                    (par.segments_length * 1e6 * 60));         %number of segments in which data is cutted
            segmentLength = floor (lts/obj.max_segments);

            tmin = 1 : segmentLength :lts;
            tmin = tmin(1:obj.max_segments);
            tmax = tmin - 1;
            tmax = tmax (2:end);
            tmax = [tmax, lts];
            obj.recmax = tmax;
            obj.recmin = tmin;
            obj.TimeStamps = TimeStamps;
            header = textread(obj.raw_filename,'%s',50);
			obj.num_scale_factor = 1e6 * str2num(header{find(strcmp('-ADBitVolts', header))+1});
			%for older versions:
            %scale_factor = textread(obj.raw_filename,'%s',43);
			%if(str2num(scale_factor{41})*1e6 > 0.5)
            %    obj.num_scale_factor = 1e6 * str2num(scale_factor{43}); %for the new CSC format
            %else
            %    obj.num_scale_factor = 1e6 * str2num(scale_factor{41}); %for the old CSC format
            %end

        end

        function [sr,max_segments,with_raw,with_spikes] = get_info(obj)
        	sr = obj.sr;
            max_segments = obj.max_segments;
            with_raw = true;
            with_spikes = false;
        end

        function index_ts = index2ts(obj,index,i)
            index_ts = (obj.TimeStamps(obj.recmin(i)+floor(index/512))'+mod(index,512)/512*obj.dt)/1000;
        end

        function x = get_segment(obj,i)

            Samples = fread(obj.opened_file,512*(obj.recmax(i)- ...
                obj.recmin(i)+1),'512*int16=>int16',8+4+4+4); %the recmax are in timestamps indexs

            x = double(Samples(:))' *  obj.num_scale_factor;

            if i == obj.max_segments
                fclose(obj.opened_file);
            end
        end
    end
end

author: csn-le


text:
function paramval = fl_inputparser(parseobj,paramname,default,validator)
% if default == [], then parameter is required

%% if no input return default
if isempty(parseobj)
    paramval = default;
    return
end

%% find parameter value
if isstruct(parseobj{1}) %if input is given as a structure

    if isfield(parseobj{1},paramname)
        paramval = parseobj{1}.(paramname);
    end

elseif iscell(parseobj) %if input is given as name/value pairs

    %check whether input is in pairs
    if mod(length(parseobj),2)
        error('Optional input parameters should come in name/value pairs.');
    end

    %separate parameters and values
    paramnames = parseobj(1:2:end);
    paramvalues = parseobj(2:2:end);

    %check if parameter name is invalid
    if isempty(paramname)
        for i = 1:length(paramnames)
            if ischar(paramnames{i}) & ~nnz(ismember(validator,paramnames{i})) %if parameter name not valid
                error(['Input ''' paramnames{i} ''' not valid.']);
            end
        end
        paramval = [];
        return
    end

    ndx = find(strcmp(paramnames,paramname));

    if length(ndx)>1
        error('Optional input parameters should only be defined once; multiple definitions detected.');
    end
    if ndx
        paramval = paramvalues{ndx};
    end

end

%% check if paramval was not defined
if ~exist('paramval')
    if ~exist('default') | isempty('default') %check if parameter is required
        error(['Input parameter '''  paramname ''' is required.']);
    else
        paramval = default;
    end
end

%% validate input
if exist('validator') %validator is optional
    if isa(validator,'function_handle') %if validator is a function handle
        if ~validator(paramval) %if not valid
            error(['Optional input parameter '''  paramname ''' has invalid value.']);
        end
    elseif iscell(validator) %directly provide cell; makes it easy to check strings

        paramval = validatestring(paramval,validator);

    end
end

author: brainstorm-tools


text:
%%%  PosControl_Sim
clear
path('./icon/',path);
Init;

% constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
% throttle when UAV is hovering
THR_HOVER = 0.609;

%% control parameter
% attitude PID parameters
Kp_PITCH_ANGLE = 6.5;
Kp_PITCH_AngleRate = 0.1;
Ki_PITCH_AngleRate = 0.02;
Kd_PITCH_AngleRate = 0.001;
Kp_ROLL_ANGLE = 6.5;
Kp_ROLL_AngleRate = 0.1;
Ki_ROLL_AngleRate = 0.02;
Kd_ROLL_AngleRate = 0.001;
Kp_YAW_AngleRate = 0.5;
Ki_YAW_AngleRate = 0.01;
Kd_YAW_AngleRate = 0.00;
% position PID parameters
Kpxp = 1.0;
Kpyp = 1.0;
Kpzp = 4.0;
Kvxp = 2.5; Kvxi = 0.4; Kvxd = 0.01;
Kvyp = 2.5; Kvyi = 0.4; Kvyd = 0.01;
Kvzp = 0.45; Kvzi = 0.01; Kvzd = 0.005;
% integral saturation
Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
Saturation_I_ah = 3.43;
Saturation_I_az = 5;

% max control angle, default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
% max control angle rate, rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
% max control speed, m/s
MAX_CONTROL_VELOCITY_XY = 5;
MAX_CONTROL_VELOCITY_Z = 3;
% throttle amplitude
MAX_MAN_THR = 0.9;
MIN_MAN_THR = 0.05;
%% run simulink model
PosControl_Sim

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function ref = physical_to_reference_t3 ( t, n, phy )

%*****************************************************************************80
%
%% PHYSICAL_TO_REFERENCE_T3 maps a physical point to a reference point.
%
%  Discussion:
%
%    Given the vertices of an order 3 physical triangle and a point
%    (X,Y) in the physical triangle, the routine computes the value
%    of the corresponding image point (XSI,ETA) in reference space.
%
%    Note that this routine may also be appropriate for an order 6
%    triangle, if the mapping between reference and physical space
%    is linear.  This implies, in particular, that the sides of the
%    image triangle are straight and that the "midside" nodes in the
%    physical triangle are halfway along the sides of
%    the physical triangle.
%
%  Reference Element T3:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  |  \
%    |  |   \
%    |  |    \
%    0  1-----2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the X and Y coordinates
%    of the vertices.  The vertices are assumed to be the images of
%    (0,0), (1,0) and (0,1) respectively.
%
%    Input, integer N, the number of points to transform.
%
%    Input, real PHY(2,N), the coordinates of points in the physical space.
%
%    Output, real REF(2,N), the coordinates of the corresponding
%    points in the reference space.
%
  ref(1,1:n) = ( ( t(2,3) - t(2,1) ) * ( phy(1,1:n) - t(1,1) )   ...
               - ( t(1,3) - t(1,1) ) * ( phy(2,1:n) - t(2,1) ) ) ...
             / ( ( t(2,3) - t(2,1) ) * ( t(1,2)   - t(1,1) )   ...
               - ( t(1,3) - t(1,1) ) * ( t(2,2)   - t(2,1) ) );

  ref(2,1:n) = ( ( t(1,2) - t(1,1) ) * ( phy(2,1:n) - t(2,1) )   ...
               - ( t(2,2) - t(2,1) ) * ( phy(1,1:n) - t(1,1) ) ) ...
             / ( ( t(2,3) - t(2,1) ) * ( t(1,2)   - t(1,1) )   ...
               - ( t(1,3) - t(1,1) ) * ( t(2,2)   - t(2,1) ) );

  return
end


text:
warning('off','all');

tic;
image = double(imread('forest.jpg'))/255;

image = imresize(image, 0.1);

result = dehaze(image, 0.95, 15);
toc;

figure, imshow(image)
figure, imshow(result)

warning('on','all');

author: sjtrny

repo: Dark-Channel-Haze-Removal


text:
function wulff()
    % wulff -- Program for plotting a Wulff net
    % to plot points, first calculate theta = pi*(90-azimuth)/180
    % then rho = tan(pi*(90-dip)/360), and finally the components
    % xp = rho*cos(theta) and yp = rho*cos(theta)
    % turned into function by Celso G Reyes 2017

    N = 50;
    cx = cos(0:pi/N:2*pi);                           % points on circle
    cy = sin(0:pi/N:2*pi);
    xh = [-1 1];                                     % horizontal axis
    yh = [0 0];
    xv = [0 0];                                      % vertical axis
    yv = [-1 1];
    axis([-1 1 -1 1]);
    axis('square');
    plot(xh,yh,'-k',xv,yv,'-k');                     %plot green axes
    axis off;
    set(gca,'NextPlot','add');
    plot(cx,cy,'-k');                                %plot white circle
    psi = 0:pi/N:pi;
    for i = 1:8                                      %plot great circles
        rdip = i*(pi/18);                             %at 10 deg intervals
        radip = atan(tan(rdip)*sin(psi));
        rproj = tan((pi/2 - radip)/2);
        x1 = rproj .* sin(psi);
        x2 = rproj .* (-sin(psi));
        y = rproj .* cos(psi);
        plot(x1,y,':k',x2,y,':k');
    end
    for i = 1:8                                     %plot small circles
        alpha = i*(pi/18);
        xlim = sin(alpha);
        % ylim = cos(alpha);
        x = -xlim:0.01:xlim;
        d = 1/cos(alpha);
        rd = d*sin(alpha);
        y0 = sqrt(rd*rd - (x .* x));
        y1 = d - y0;
        y2 = - d + y0;
        plot(x,y1,':k',x,y2,':k');
    end
    axis('square');
    set(gcf,'color','w');
end

author: CelsoReyes


text:
classdef ParadigmPAL < ParadigmBaseSimplified
    % Experimental implementation of the Pattern Alignment Learning (PAL) method.
    %
    % This paradigm is not yet completed, please move on! :-)
    %
    % Name:
    %   Wave Alignment, work in progress
    %
    %                            Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
    %                            2011-08-28

    methods

        function defaults = preprocessing_defaults(self)
            % define the default pre-processing parameters of this paradigm
            defaults = {'SpectralSelection',[0.1 35],'EpochExtraction',[-0.5 1.5], 'Resampling',100};
        end

        function model = calibrate_simple(self,varargin)
            % Override this function to implement your calibration code
            % Model = calibrate_simple(Signal,Arguments...)
            %
            % In:
            %   Signal : a single continuous EEGLAB data set (annotated with target markers; see
            %            set_targetmarkers for more info)
            %
            %   Arguments : further optional user arguments
            %
            % Out:
            %   Model : a model struct with a mandatory field .filter_graph and arbitrary other content
            %           * The .filter_graph filed is a 1x1 cell array that contains the desciption of
            %             filter steps that is to be applied to the data, and is usually passed as either
            %             the .tracking.online_expression field of the processed Signal or the processed
            %             Signal itself. If no signal processing is performed by this paradigm, the raw
            %             signal may be passed.
            %
            %           * May have optional fields .prediction_function, .prediction_window and
            %             .prediction_channels - though these are generally auto-deduced
            %
            % Notes:

            args = arg_define(varargin, ...
                arg_norep('signal'), ...
                arg_sub({'flt','SignalProcessing'}, self.preprocessing_defaults(), @flt_pipeline, 'Signal processing stages. These parameters control filter stages that run on the signal level; they can be enabled, disabled and configured for the given paradigm. The prediction operates on the outputs of this stage.'), ...
                arg({'lambdas','Lambdas'},[2.^(8:-0.125:-5)],[0 Inf],'Regularization parameter. Larger means stronger regularization, 1.0 is the Bayesian ad hoc choice.'), ...
                arg({'duration','WaveDuration'},0.8,[0 Inf],'Duration of Waveform. This is the length of the waveform that will be searched in the epoch.'), ...
                arg({'stepping','Stepping'},3,uint32([1 100]),'Relative jitter stepping. This is the precision (in samples) at which relative jitters will be considered.'), ...
                arg({'admm_iters','ADMMIterations'},30,uint32([1 1000]),'Max Iterations. Number of iterations for the overall ADMM solver.'), ...
                arg({'inner_iters','InnerIterations'},10,uint32([1 1000]),'Max Inner Iterations. Number of iterations for the inner ADMM solver that solves the sparsity vs trace norm.'), ...
                arg({'cg_iters','CgIterations'},30,uint32([1 1000]),'Max CG Iterations. Number of iterations for the inner conjugate-gradient solver.'), ...
                arg({'lanczos_iters','LanczosIterations'},10,uint32([1 1000]),'Max Lanczos Iterations. Number of iterations for the inner Lanczos solver.'), ...
                arg({'abs_tol','AbsoluteTolerance'},1e-6,[],'Absolute tolerance.'), ...
                arg({'rel_tol','RelativeTolerance'},1e-3,[],'Relative tolerance.'), ...
                arg({'alpha','OverRelaxation'},1,[1 1.8],'Over-relaxation parameter. Some references suggest that values between 1.5 and 1.8 can improve convergence.'), ...
                arg({'rho','AugmentedLagrangian'},1.0,[0 Inf],'Augmented Lagrangian parameter.'), ...
                arg({'rho_update','RhoUpdate'},true,[],'Update Rho. Whether to update rho dynamically according to 3.4.1 in [1]'), ...
                arg({'rho_cutoff','RhoUpdateThreshold'},10.0,[0 Inf],'Rho update threshold.'), ...
                arg({'rho_incr','RhoUpdateIncr'},2.0,[1 Inf],'Rho update increment factor.'), ...
                arg({'rho_decr','RhoUpdateDecr'},2.0,[1 Inf],'Rho update decrement factor.'), ...
                arg({'kappa','SparsityTraceTradeoff'},0.75,[0 1],'Sparsity (=1) vs. Trace norm (=0) tradeoff. Between 0 and 1, similar to the elastic net criterion.'), ...
                arg({'patterns','NumPatterns'},20,uint32([1 1000]),'Number of weight patterns. This is the maximum number of (possibly distinct) weight patterns to recover.'), ...
                arg({'verbose','Verbose'},'really',[],'Verbose output.'), ...
                arg({'scaling','Scaling'}, 'std', {'none','center','std','minmax','whiten'}, 'Pre-scaling of the data. For the regulariation to work best, the features should either be naturally scaled well, or be artificially scaled.'), ...
                arg({'fast','Fast'},true,[],'Fast mode.'), ...
                arg({'arg_dialogsel','ConfigLayout'},self.dialog_layout_defaults(),[],'Parameters displayed in the config dialog. Cell array of parameter names to display (dot-notation allowed); blanks are translated into empty rows in the dialog. Referring to a structure argument lists all parameters of that struture, except if it is a switchable structure - in this case, a pulldown menu with switch options is displayed.','type','cellstr','shape','row'));

            % first pre-process the data (symbolically)
            signal = flt_pipeline('signal',args.signal, args.flt); %#ok<*NODEF>

            % evaluate this in an optimized fashion
            signal = exp_eval_optimized(signal);

            % get trial labels
            labels = set_gettarget(signal);

            % build meta-data
            model.args = rmfield(args,'signal');
            model.classes = unique(labels,'rows');
            model.tracking.filter_graph = signal;

            if length(model.classes) > 2
                error('Currently, only two classes are supported here.'); end

            % extract time-shifted windows from each trial...
            duration = round(args.duration * signal.srate);
            offsets = 1 : args.stepping : (signal.pnts - duration);
            trials = cell(1,length(offsets));

            % get the entire epoch and figure out an appropriate rescaling
            alldata = reshape(permute(signal.data,[2 1 3]),[],signal.trials);
            sc_info = hlp_findscaling(alldata,args.scaling);

            % for each time shift...
            for o = offsets
                % determine the window to cut out of the data
                cutwindow = o + (1:duration);
                % extract from all trials
                data = signal.data(:,cutwindow,:);
                % vectorize into #trials x #features, and rescale
                data = permute(data,[2 1 3]);
                trials{o} = hlp_applyscaling(reshape(data,[],signal.trials),sc_info)';
            end

            T = signal.trials;
            fprintf('#trials=%.0f; #shifts=%.0f; #features=%.0f; #channels=%.0f; #windows=%.0f\n',T,length(offsets),signal.nbchan*length(cutwindow),signal.nbchan,length(cutwindow));

            % vertically concatenate all shifts (into #shifts*#trials x #features)
            A = vertcat(trials{:});

            % remap the labels to -1 / +1
            labels(labels == model.classes(1)) = -1;
            labels(labels == model.classes(2)) = +1;

            % and generate shift-replicated target matrix b
            b = repmat(labels,1,length(offsets));
            b = b(:);

            C0 = -[b A];   % (bias is first elem, followed by rest...)
            ccp = C0';

            % weight vectors
            [m,n] = size(A);
            Wz = zeros(n+1,m);          % consensus variable
            Wx = zeros(n+1,m);          % data weights variable, FxN; n+1 accounts for the added bias (this is actually just the initial conditions vector...)
            Wu = zeros(n+1,m);          % scaled dual parameter u for Wx

            alpha = args.alpha;         % over-relaxation parameter (if 1, OR turned off)
            rho = args.rho;             % augmented Lagrangian parameter; probably determines sort of the tightness of the constraints...
            kappa = args.kappa;         % blend factor between trialwise sparsity and trialwise agreement
            K = args.patterns;          % the maximum number of SV's to recover

            rescale = 1;                  % scale of the left-hand side term (could be m/length(offsets)
            for mu = args.lambdas       % regularization parameter (should be decreasing)
                fprintf('Now using lambda = %.3f\n',mu);
                if args.verbose
                    fprintf('%3s\t%10s\t%10s\t%10s\t%10s\t%10s\n', 'iter', ...
                        'r norm', 'eps pri', 's norm', 'eps dual', 'objective');
                end

                for k = 1:args.admm_iters
                    Wzold = Wz;
                    fprintf('Outer step #%.0f\n',k);

                    %% x-update
                    Wx = ParadigmPAL.minimize_batch(Wx, @ParadigmPAL.x_objective_batch, args.cg_iters, Wz, Wu, ccp, rho);
                    Wx_hat = alpha*Wx + (1-alpha)*Wzold;  % (optional over-relaxation)

                    %% z-update
                    Wz = Wx_hat + Wu;
                    [Wz,fuuu,U,V] = ParadigmPAL.prox_sum(Wz,rescale,mu,kappa,rho,K,T,m,args);

                    %% dual parameter update
                    fprintf('  updating u''s\n');
                    Wu = Wu + (Wx_hat - Wz);

                    %% diagnostics, reporting, termination checks
                    fprintf('  diagnostics...\n');
                    history.objval(k)  = 0; %ParadigmPAL.full_objective(ccp, mu, Wx, Wz);
                    history.r_norm(k)  = norm(Wx(:) - Wz(:));
                    history.s_norm(k)  = norm(rho*(Wz(:) - Wzold(:)));
                    history.eps_pri(k) = args.abs_tol + args.rel_tol*max(norm(Wx(:)), norm(Wz(:)));
                    history.eps_dual(k)= args.abs_tol + args.rel_tol*norm(rho*Wu(:));

                    if args.verbose
                        fprintf('%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n', k, ...
                            history.r_norm(k), history.eps_pri(k), ...
                            history.s_norm(k), history.eps_dual(k), history.objval(k));
                        if strcmp(args.verbose,'really') && nnz(fuuu)
                            f=figure; set(f,'Position',get(f,'Position')-[0 500 0 0]);
                            imagesc(Wz); ylabel('features'); xlabel('shifted trials'); title(sprintf('ADMM iter #%.0f, lambda=%.3f',k,mu));
                            figure;
                            % weight patterns
                            N = nnz(fuuu);
                            cols = ceil(sqrt(N));
                            rows = ceil(N/cols);
                            for p=1:N
                                subplot(cols,rows,p);
                                plot((0:round(args.duration * signal.srate)-1)/signal.srate,reshape(U(2:end,p),length(cutwindow),signal.nbchan)); title(sprintf('weight array #%.0f',p)); xlabel('time (s)'); ylabel('weight');
                            end
                            % topoplots...
                            f=figure; set(f,'Position',get(f,'Position')-[550 0 0 0]);
                            N = nnz(fuuu);
                            cols = ceil(sqrt(N));
                            rows = ceil(N/cols);
                            for p=1:N
                                subplot(cols,rows,p);
                                mat = reshape(U(2:end,p),length(cutwindow),signal.nbchan);
                                mat_inner = mat(round(end*1/4):round(end*3/4),:);
                                [dummy,peakidx] = max(mean(abs(mat_inner),2)); % #ok<NASGU>
                                topoplot(mat_inner(peakidx,:),signal.chanlocs);
                            end
                            drawnow;
                        end
                    end

                    if history.r_norm(k) < history.eps_pri(k) && history.s_norm(k) < history.eps_dual(k)
                        break; end

                    %% update of rho
                    if args.rho_update
                        if history.r_norm(k) > args.rho_cutoff * history.s_norm(k)
                            disp('  incrementing rho.');
                            rho = rho .* args.rho_incr;
                            Wu = Wu ./ args.rho_incr;
                        elseif history.s_norm(k) > args.rho_cutoff * history.r_norm(k)
                            disp('  decrementing rho.');
                            rho = rho ./ args.rho_incr;
                            Wu = Wu .* args.rho_incr;
                        end
                    end
                end
            end
            1
            %% plotting...

%             % full weight matrix
%             figure;imagesc(Wz); title('full weight matrix'); xlabel('features'); ylabel('shifted trials');
%
%             % shift matrix (?)
%             figure; imagesc(reshape(V(:,1),length(offsets),[])); title('shift matrix'); xlabel('trials'); ylabel('shifts');
        end

        function outputs = predict_simple(self,signal,model)
            % Override this function to implement your prediction code
            % Outputs = predict_simple(Sginal,Model)
            %
            % In:
            %   Signal : a signal pre-processed according to the model's filter graph
            %
            %   Model : a predictive model as created by your calibrate_simple() function
            %
            % Out:
            %   Outputs : a prediction/estimate for the most recent time point in the data (or one for
            %             every epoch if the signal is epoched); see ml_predict for the allowed formats

            error('This paradigm implements no predict() function.');
        end

        function visualize(self,model)
            % Optionally override this function to implement your visualization code
            % visualize(Model)
            %
            % In:
            %   Model: a model as created by your calibrate() function;
            %          a plot or GUI will be produced to inspect the model
            %

            error('This paradigm implements no visualize() function.');
        end

        % --- implementation ----

        function model = calibrate(self,varargin)
            % Implementation of the calibrate() function (see ParadigmBase)
            % Model = calibrate(Collection,GoalIdentifier,Arguments...)
            %
            % This is essentially a wrapper around calibrate_simple() which peels off any additional
            % data structure (for multiple collections and stream bundles) before passing it on.
            %
            % In:
            %   Collection : a collection (cell array) of stream bundles
            %
            %   GoalIdentifier: the goal-identifier
            %
            %   Arguments... : further optional user arguments
            %
            % Out:
            %   Model : The returned model struct

            % get the collection argument and extract the signal argument
            collection = arg_extract(varargin,{'collection','Collection'},[],{});
            if ~isempty(collection)
                if length(collection) > 1
                    error('This paradigm does not support collections of more than one recording.'); end
                if length(collection{1}.streams) > 1
                    error('This paradigm does not support more than one stream in parallel.'); end
                signal = collection{1}.streams{1};
            else
                signal = {};
            end

            % call calibrate_simple with the signal passed in as additional argument
            model = self.calibrate_simple('signal',signal,varargin{:});
        end

        function outputs = predict(self,bundle,model)
            % Override this function to implement your prediction code
            % Outputs = predict(Bundle,Model)
            %
            % This is a wrapper around the predict_simple() function which peels off the stream bundle
            % representation and just passes on the contained signal.
            %
            % In:
            %   Bundle : stream bundle preprocessed by the model's filter graph
            %
            %   Model : a predictive model
            %
            % Out:
            %   Outputs : the outputs of calibrate_simple()

            if length(bundle.streams) > 1
                error('This paradigm does not support more than one stream in parallel.'); end
            outputs = self.predict_simple(bundle.streams{1},model);
        end

    end

    methods(Static)
        function obj = full_objective(ccp, mu, Wx, Wz)
            % this is the overall objective function of alignment learning
            obj = sum(log(1 + exp(sum(ccp.*Wx)))) + mu*(norm_nuc(Wz) + sum(norms(Wz))) * size(ccp,2);
        end

        function [val,grad] = x_objective(x, C, z, u, rho)
            % this is the objective function for the decoupled x criterion (basically l2 logreg)
            ecx = exp(C*x);
            val = sum(log(1 + ecx)) + (rho/2)*norm(x - z + u).^2;
            grad = C'*(ecx./(1 + ecx)) + rho*(x - z + u);
        end

        function [vals,grads] = x_objective_batch(Wx,Wz,Wu,ccp,rho)
            % this is the batch version of the above (for all samples & variables in parallel)
            ecx = exp(sum(ccp.*Wx));
            d = Wx - Wz + Wu;
            vals = log(1 + ecx) + (rho/2)*sum(d.*d);
            grads = bsxfun(@times,ccp,(ecx./(1 + ecx))) + rho*d;
        end

        % l1/l2 shrinkage operator
        function Wz = shrinkage(Wa, kappa)
            Wz = pos(1 - kappa/norm(Wa))*Wa;
        end

        % l1 shrinkage operator
        function z = shrinkage_l1(a, kappa)
            z = max(0, a-kappa) - max(0, -a-kappa);
        end

        % trace norm proximal operator
        function [W,fuuu,U,V] = prox_tr(W,rescale,mu,kappa,rho,K,args)
            fprintf('  trace norming r\n');
            [U,S,V] = pca(W,K,args.lanczos_iters);
            fuuu = ParadigmPAL.shrinkage_l1(diag(S),(rescale*mu*(1-kappa))/rho); nn = length(fuuu); % we shrink the SV spectrum using prox operator
            Ws = sparse(1:nn,1:nn,fuuu,size(S,1),size(S,2));
            % clamp minority of V's members at 0
            smv = sign(median(V));
            V(:,smv==-1) = min(0,V(:,smv==-1));
            V(:,smv==1) = max(0,V(:,smv==1));
            W = U*(Ws*V');

            if strcmp(args.verbose,'ultra')
                f=figure; set(f,'Position',get(f,'Position')-[0 500 0 0]);
                imagesc(W); ylabel('features'); xlabel('shifted trials'); title(sprintf('lambda=%.3f',mu));
                figure;
                % weight patterns
                N = nnz(spec);
                cols = ceil(sqrt(N));
                rows = ceil(N/cols);
                signal=args.signal;
                for p=1:N
                    subplot(cols,rows,p);
                    plot((0:(numel(U(2:end,p))/signal.nbchan)-1)/signal.srate,reshape(U(2:end,p),[],signal.nbchan)); title(sprintf('weight array #%.0f',p)); xlabel('time (s)'); ylabel('weight');
                end
                drawnow;
            end

        end

        % group-sparse proximal operator
        function W = prox_gl(W,rescale,mu,kappa,rho,T,m)
            fprintf('  sparsifying z\n');
            for i = 1:m
                W(2:end,i) = ParadigmPAL.shrinkage(W(2:end,i), (rescale*mu*kappa/10)/rho); end
        end

        % summed proximal operator (using some splitting scheme, see Bauscke, 2008)
        function [x,fuuu,U,V] = prox_sum(r,rescale,mu,kappa,rho,K,T,m,args)
            y = r;
            p = zeros(size(y));
            q = zeros(size(y));
            for n=1:args.inner_iters
                [x,fuuu,U,V] = ParadigmPAL.prox_tr(y+q,rescale,mu,kappa,rho,K,args);
                q = y+q-x;
                y = ParadigmPAL.prox_gl(x+p,rescale,mu,kappa,rho,T,m);
                p = x+p-y;
                if norm(x(:) - y(:)) < args.rel_tol*max(norm(x(:)),norm(y(:)))
                    break; end
            end
        end

        function X = minimize_batch(X, f, len, Wz,Wu,ccp,rho)
            % Copyright (C) 2001 - 2010 by Carl Edward Rasmussen, 2010-01-03

            INT = 0.1;    % don't reevaluate within 0.1 of the limit of the current bracket         % const
            EXT = 3.0;                  % extrapolate maximum 3 times the current step-size         % const
            MAX = 20;                         % max 20 function evaluations per line search         % const
            RATIO = 10;                                       % maximum allowed slope ratio         % const
            SIG = 0.1; RHO = SIG/2; % SIG and RHO are the constants controlling the Wolfe-          % const
            if max(size(len)) == 2, red=len(2); len=len(1); else red=1; end                         % scalar

            [f0 df0] = feval(f,X,Wz,Wu,ccp,rho);          % get function value and gradient
            N = length(f0);

            [A,B,d1,d2,d3,d4,f1,f2,f3,f4,x1,x2,x4] = deal(zeros(1,N)); % initialize a few variables...
            df3 = zeros(size(X));
            i = 0;                                            % zero the run length counter         % scalar
            ls_failed = false(1,N);                             % no previous line search has failed% boolean vector
            i = i + (len<0);                                            % count epochs?!            % scalar
            s = -df0;                                 % initial search direction (steepest)         % column MATRIX
            d0 = -sum(s.*s);                                                    % and slope         % row vector
            x3 = red./(1-d0);                                  % initial step is red/(|s|+1)        % row vector

            keep_going = true(1,N);                                                                 % mask for the outermost loop updates...
            while any(keep_going) && i < abs(len)                                      % while not finished            % scalar
                i = i + (len>0);                                      % count iterations?!          % scalar
                fprintf('  Cg step #%.0f\n',i);

                X0 = X;                                                                             % column MATRIX
                F0 = f0;                                                                            % row vector
                dF0 = df0;                   % make a copy of current values                        % column MATRIX
                if len>0, M = MAX; else M = min(MAX, -len-i); end                                   % scalar

                om = keep_going;                                                                    % mask for outer-loop updates
                while 1                             % keep extrapolating as long as necessary
                    x2(om) = 0;                                                                     % row vector
                    f2(om) = f0(om);                                                                % row vector
                    d2(om) = d0(om);                                                                % row vector
                    f3(om) = f0(om);                                                                % row vector
                    df3(:,om) = df0(:,om);                                                          % column matrix
                    ok(om) = false;                                                                 % boolean vector
                    while ~all(ok(om)) && M > 0                                                     % mask check
                        M = M - 1; i = i + (len<0);                         % count epochs?!        % scalar
                        nokom = ~ok&om;
                        [f3(nokom) df3(:,nokom)] = feval(f, X(:,nokom) + bsxfun(@times,x3(nokom),s(:,nokom)), Wz(:,nokom), Wu(:,nokom), ccp(:,nokom), rho); % masked array op
                        ok(om) = ok(om) | ~(isnan(f3(om)) | isinf(f3(om)) | any(isnan(df3(om))+isinf(df3(om))));            % mask update
                        x3(~ok&om) = (x2(~ok&om)+x3(~ok&om))/2;         % bisect and try again      % maked update
                    end
                    mask = f3<F0;                                       % keep best values          % mask generation
                    X0(:,mask&om) = X(:,mask&om) + bsxfun(@times,x3(mask&om),s(:,mask&om));         % masked updates
                    F0(mask&om) = f3(mask&om);
                    dF0(:,mask&om) = df3(:,mask&om);

                    d3(om) = sum(df3(:,om).*s(:,om));                                    % new slope% row vector

                    om(d3 > SIG*d0 | f3 > f0 + x3.*d0*RHO) = false;             % update outer mask
                    if ~any(om) || M == 0  % are we done extrapolating?
                        break; end

                    x1(om) = x2(om); f1(om) = f2(om); d1(om) = d2(om); % move point 2 to point 1    % row vectors
                    x2(om) = x3(om); f2(om) = f3(om); d2(om) = d3(om); % move point 3 to point 2    % row vectors
                    A(om) = 6*(f1(om)-f2(om))+3*(d2(om)+d1(om)).*(x2(om)-x1(om));% make cubic extrap.% row vectors
                    B(om) = 3*(f2(om)-f1(om))-(2*d1(om)+d2(om)).*(x2(om)-x1(om));
                    x3(om) = x1(om)-d1(om).*(x2(om)-x1(om)).^2./(B(om)+sqrt(B(om).*B(om)-A(om).*d1(om).*(x2(om)-x1(om)))); % num. error possible, ok!

                    % implement fixups: num prob | wrong sign | new point beyond extrapolation limit?
                    issue_mask = ~isreal(x3) | isnan(x3) | isinf(x3) | x3 < 0 | x3 > x2*EXT;
                    x3(om & issue_mask) = x2(om & issue_mask)*EXT;
                    % new point too close to previous point?
                    issue_mask = ~issue_mask & (x3 < x2+INT*(x2-x1));
                    x3(om & issue_mask) = x2(om & issue_mask) + INT*(x2(om & issue_mask) - x1(om & issue_mask));
                end                                                       % end extrapolation

                while 1
                    om = keep_going & (abs(d3) > -SIG*d0 | f3 > f0+x3.*d0*RHO);    % update outer mask
                    if ~(any(om) && M > 0)
                        break; end
                    tog = d3 > 0 | f3 > f0+x3.*d0*RHO;           % subinterval toggle
                    x4(om&tog) = x3(om&tog);
                    f4(om&tog) = f3(om&tog);
                    d4(om&tog) = d3(om&tog);
                    x2(om&~tog) = x3(om&~tog);
                    f2(om&~tog) = f3(om&~tog);
                    d2(om&~tog) = d3(om&~tog);
                    tog=f4>f0;                                   % interpolation toggle
                    x3(om&tog) = x2(om&tog) - (0.5*d2(om&tog).*(x4(om&tog)-x2(om&tog)).^2) ./ (f4(om&tog)-f2(om&tog)-d2(om&tog).*(x4(om&tog)-x2(om&tog)));  % quadratic interpolation
                    A(om&~tog) = 6*(f2(om&~tog)-f4(om&~tog))./(x4(om&~tog)-x2(om&~tog)) + 3*(d4(om&~tog)+d2(om&~tog));                    % cubic interpolation
                    B(om&~tog) = 3*(f4(om&~tog)-f2(om&~tog)) - (2*d2(om&~tog)+d4(om&~tog)).*(x4(om&~tog)-x2(om&~tog));
                    x3(om&~tog) = x2(om&~tog) + (sqrt(B(om&~tog).*B(om&~tog) - A(om&~tog).*d2(om&~tog).*(x4(om&~tog)-x2(om&~tog)).^2)-B(om&~tog))./A(om&~tog);  % num. error possible, ok!
                    numprob = isnan(x3) | isinf(x3);             % if we had a numerical problem then bisect
                    x3(om&numprob) = (x2(om&numprob)+x4(om&numprob))/2;

                    x3(om) = max(min(x3(om), x4(om)-INT*(x4(om)-x2(om))),x2(om)+INT*(x4(om)-x2(om)));  % don't accept too close
                    [f3(om) df3(:,om)] = feval(f, X(:,om) + bsxfun(@times,x3(om),s(:,om)), Wz(:,om), Wu(:,om), ccp(:,om), rho); % masked array op
                    mask = om & f3<F0;                                      % keep best values
                    X0(:,mask) = X(:,mask) + bsxfun(@times,x3(mask),s(:,mask));
                    F0(mask) = f3(mask);
                    dF0(:,mask) = df3(:,mask);
                    M = M - 1; i = i + (len<0);                             % count epochs?!
                    d3(om) = sum(df3(:,om).*s(:,om));                                    % new slope% row vector
                end

                % --- we are here ---

                okay = abs(d3) < -SIG*d0 & f3 < f0+x3.*d0*RHO;

                upd = okay & keep_going;
                X(:,upd) = X(:,upd) + bsxfun(@times,x3(upd),s(:,upd));
                f0(upd) = f3(upd);
                s(:,upd) = bsxfun(@times,sum(df3(:,upd).*df3(:,upd)) - sum(df0(:,upd).*df3(:,upd)) ./ sum(df0(:,upd).*df0(:,upd)),s(:,upd)) - df3(:,upd);
                df0(:,upd) = df3(:,upd);                                % swap derivatives
                d3(upd) = d0(upd);
                d0(upd) = sum(df0(:,upd).*s(:,upd));
                mask = upd & d0 > 0;                         % otherwise use steepest direction
                s(:,mask) = -df0(:,mask);
                d0(mask) = -sum(s(:,mask).*s(:,mask));
                x3(upd) = x3(upd) .* min(RATIO, d3(upd)./(d0(upd)-realmin));  % slope ratio but max RATIO
                ls_failed(upd) = false;                              % this line search did not fail

                upd = ~okay & keep_going;
                X(:,upd) = X0(:,upd);
                f0(upd) = F0(upd);
                df0(:,upd) = dF0(:,upd);

                keep_going(upd) = keep_going(upd) & ~ls_failed(upd); % line search failed twice in a row
                                                                     % so we give up on these guys
                upd = upd & keep_going;
                if i <= abs(len)
                    s(:,upd) = -df0(:,upd);
                    d0(:,upd) = -sum(s(:,upd).*s(:,upd));                             % try steepest
                    x3(upd) = 1./(1-d0(upd));
                    ls_failed(upd) = true;
                end
            end
        end

        function layout = dialog_layout_defaults(self)
            layout = {'SignalProcessing.Resampling.SamplingRate','SignalProcessing.SpectralSelection.FrequencySpecification', ...
                'SignalProcessing.EpochExtraction','','Lambdas','WaveDuration','Stepping','ADMMIterations','CgIterations',...
                'LanczosIterations','NumPatterns','SparsityTraceTradeoff'};
        end

    end
end

% (turn off a few editor warnings because some actual implementations are missing in this file)
%#ok<*INUSD,*STOUT,*MANU>

repo: REST


text:
function histmat  = hist2D(x, y, xedges, yedges)
% function histmat  = hist2D(x, y, xedges, yedges)
%
% Extract 2D histogram data containing the number of events
% of [x , y] pairs that fall in each bin of the grid defined by
% xedges and yedges. The edges are vectors with monotonically
% non-decreasing values.
%
%EXAMPLE
%
% events = 1000000;
% x1 = sqrt(0.05)*randn(events,1)-0.5; x2 = sqrt(0.05)*randn(events,1)+0.5;
% y1 = sqrt(0.05)*randn(events,1)+0.5; y2 = sqrt(0.05)*randn(events,1)-0.5;
% x= [x1;x2]; y = [y1;y2];
%
%For linearly spaced edges:
% xedges = linspace(-1,1,64); yedges = linspace(-1,1,64);
% histmat = hist2(x, y, xedges, yedges);
% figure; pcolor(xedges,yedges,histmat'); colorbar ; axis square tight ;
%
%For nonlinearly spaced edges:
% xedges_ = logspace(0,log10(3),64)-2; yedges_ = linspace(-1,1,64);
% histmat_ = hist2(x, y, xedges_, yedges_);
% figure; pcolor(xedges_,yedges_,histmat_'); colorbar ; axis square tight ;

% University of Debrecen, PET Center/Laszlo Balkay 2006
% email: balkay@pet.dote.hu

if nargin ~= 4
    error ('The four input arguments are required!');
    return;
end
if any(size(x) ~= size(y))
    error ('The size of the two first input vectors should be same!');
    return;
end

[xn, xbin] = histc(x,xedges);
[yn, ybin] = histc(y,yedges);

%xbin, ybin zero for out of range values
% (see the help of histc) force this event to the
% first bins
xbin(find(xbin == 0)) = inf;
ybin(find(ybin == 0)) = inf;

xnbin = length(xedges);
ynbin = length(yedges);

if xnbin >= ynbin
    xy = ybin*(xnbin) + xbin;
      indexshift =  xnbin;
else
    xy = xbin*(ynbin) + ybin;
      indexshift =  ynbin;
end

%[xyuni, m, n] = unique(xy);
xyuni = unique(xy);
xyuni(end) = [];
hstres = histc(xy,xyuni);
clear xy;

histmat = zeros(ynbin,xnbin);
histmat(xyuni-indexshift) = hstres;

author: buzsakilab


text:
function [f,bandwidth] = calcDensity(x,varargin)
% kernel density estimation from real valued data
%
% Syntax
%   f = calcDensity(x)
%   f = calcDensity(x,'range',[xmin;xmax])
%
%   f = calcDensity([x,y],'range',[xmin,ymin;xmax,ymax])
%   f = calcDensity([x,y,z],'range',[xmin,ymin,zmin;xmax,ymax,zmax])
%
% Input
%  x,y,z - random samples as n x 1 vectors
%
% Output
%  f - density as <https://de.mathworks.com/help/matlab/ref/griddedinterpolant.html griddedinterpolant>
%
% See also
% vector3d/calcDensity orientation/calcDensity

if check_option(varargin,'periodic')

  f = calcS1Density(x,varargin{:});
  return
end

range = get_option(varargin,'range',[min(x);max(x)]);
varargin = delete_option(varargin,'range',1);

% the one dimensional case
if length(x) == numel(x)

  [bandwidth,density,grid] = kde(x,2^14,range(1),range(2),varargin{:});
  grid = {grid};

else % the multidimensional case

  dim = size(x,2);
  N = round(1000000^(1/dim));
  for d = 1:dim
    gs{d} = linspace(range(1,d),range(2,d),N);
  end

  [gs{:}] = ndgrid(gs{:});
  grid = gs;

  gs = cellfun(@(y) y(:),gs,'UniformOutput',false);
  density = reshape(kdeN(x,[gs{:}]),size(grid{1}));

end

f = griddedInterpolant(grid{:},density,'spline');

author: mtex-toolbox

repo: mtex

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


text:
function d2 = dtiXformTensors(d, x)
% dt6_new = dtiXformTensors(dt6, xform)
%
% Efficiently applies the 3x3 transform to a tensor volume in dt6 format.
% dt6 is a XxYxZx3x3 array.
%
% This is just a matlab-efficient way to do xform*T*xform' on a
% bunch of tensors (T) in dt6 format, a 3d array of the 6 unique
% tensor elements in [Dxx, Dyy, Dzz, Dxy, Dxz, Dyz] order.
%
%
%
% HISTORY:
% 2003.12.15 RFD & AHS Wrote it.
% 2006.08.03 RFD: rewote it to be much more memory efficient. It's
% also now ~2.5x faster.

if(size(x,2)>3)
  x = x(1:3,1:3);
end

d2 = zeros(size(d));
d2(:,:,:,1) = x(1)*(x(1).*d(:,:,:,1) + x(4).*d(:,:,:,4) + x(7).*d(:,:,:,5)) ...
            + x(4)*(x(1).*d(:,:,:,4) + x(4).*d(:,:,:,2) + x(7).*d(:,:,:,6)) ...
	    + x(7)*(x(1).*d(:,:,:,5) + x(4).*d(:,:,:,6) + x(7).*d(:,:,:,3));

d2(:,:,:,2) = x(2)*(x(2).*d(:,:,:,1) + x(5).*d(:,:,:,4) + x(8).*d(:,:,:,5)) ...
            + x(5)*(x(2).*d(:,:,:,4) + x(5).*d(:,:,:,2) + x(8).*d(:,:,:,6)) ...
	    + x(8)*(x(2).*d(:,:,:,5) + x(5).*d(:,:,:,6) + x(8).*d(:,:,:,3));

d2(:,:,:,3) = x(3)*(x(3).*d(:,:,:,1) + x(6).*d(:,:,:,4) + x(9).*d(:,:,:,5)) ...
            + x(6)*(x(3).*d(:,:,:,4) + x(6).*d(:,:,:,2) + x(9).*d(:,:,:,6)) ...
	    + x(9)*(x(3).*d(:,:,:,5) + x(6).*d(:,:,:,6) + x(9).*d(:,:,:,3));

d2(:,:,:,4) = x(2)*(x(1).*d(:,:,:,1) + x(4).*d(:,:,:,4) + x(7).*d(:,:,:,5)) ...
            + x(5)*(x(1).*d(:,:,:,4) + x(4).*d(:,:,:,2) + x(7).*d(:,:,:,6)) ...
	    + x(8)*(x(1).*d(:,:,:,5) + x(4).*d(:,:,:,6) + x(7).*d(:,:,:,3));

d2(:,:,:,5) = x(3)*(x(1).*d(:,:,:,1) + x(4).*d(:,:,:,4) + x(7).*d(:,:,:,5)) ...
            + x(6)*(x(1).*d(:,:,:,4) + x(4).*d(:,:,:,2) + x(7).*d(:,:,:,6)) ...
	    + x(9)*(x(1).*d(:,:,:,5) + x(4).*d(:,:,:,6) + x(7).*d(:,:,:,3));

d2(:,:,:,6) = x(3)*(x(2).*d(:,:,:,1) + x(5).*d(:,:,:,4) + x(8).*d(:,:,:,5)) ...
            + x(6)*(x(2).*d(:,:,:,4) + x(5).*d(:,:,:,2) + x(8).*d(:,:,:,6)) ...
	    + x(9)*(x(2).*d(:,:,:,5) + x(5).*d(:,:,:,6) + x(8).*d(:,:,:,3));

return;

% OLD CODE (pre 2006.08.03 rewrite):

temp = zeros([size(d,1), size(d,2), size(d,3), 3, 3]);
temp(:,:,:,1,1) = x(1,1).*d(:,:,:,1) + x(1,2).*d(:,:,:,4) + x(1,3).*d(:,:,:,5);
temp(:,:,:,1,2) = x(1,1).*d(:,:,:,4) + x(1,2).*d(:,:,:,2) + x(1,3).*d(:,:,:,6);
temp(:,:,:,1,3) = x(1,1).*d(:,:,:,5) + x(1,2).*d(:,:,:,6) + x(1,3).*d(:,:,:,3);
temp(:,:,:,2,1) = x(2,1).*d(:,:,:,1) + x(2,2).*d(:,:,:,4) + x(2,3).*d(:,:,:,5);
temp(:,:,:,2,2) = x(2,1).*d(:,:,:,4) + x(2,2).*d(:,:,:,2) + x(2,3).*d(:,:,:,6);
temp(:,:,:,2,3) = x(2,1).*d(:,:,:,5) + x(2,2).*d(:,:,:,6) + x(2,3).*d(:,:,:,3);
temp(:,:,:,3,1) = x(3,1).*d(:,:,:,1) + x(3,2).*d(:,:,:,4) + x(3,3).*d(:,:,:,5);
temp(:,:,:,3,2) = x(3,1).*d(:,:,:,4) + x(3,2).*d(:,:,:,2) + x(3,3).*d(:,:,:,6);
temp(:,:,:,3,3) = x(3,1).*d(:,:,:,5) + x(3,2).*d(:,:,:,6) + x(3,3).*d(:,:,:,3);

temp = dtiXformVectors(temp, x', 'post');
d2 = temp(:,:,:,[1 5 9 4 7 8]);

% This is equivalent to:
% dt6_new = zeros(size(dt6));
% for (x=1:size(dt6,1)),
%   for (y=1:size(dt6,2)),
%       for (z=1:size(dt6,3)),
%           t = dti6to33(dt6_new(x,y,z,:));
%           t = xform * t * xform';
%           dt6_new(x,y,z,:) = dti33to6(t);
%       end
%   end
% end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function eventLog = ma_executeDVManeuver_dv_inertial_spherical(dvVectInertialSpherical, thruster, initialState, eventNum)
%ma_executeDVManeuver_dv_inertial_spherical Summary of this function goes here
%   Detailed explanation goes here

    azimuth = dvVectInertialSpherical(1);
    elevation = dvVectInertialSpherical(2);
    r = dvVectInertialSpherical(3);

    [x,y,z] = sph2cart(azimuth,elevation,r);
    dvVectInertial = [x,y,z];

    eventLog = ma_executeDVManeuver_dv_inertial(dvVectInertial, thruster, initialState, eventNum);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ vX ] = ProjectSimplexExact( vY, ballRadius )
% ----------------------------------------------------------------------------------------------- %
% [ vX ] = ProjectSimplexExact( vY, ballRadius )
%   Solving the Orthogonal Projection Problem of the input vector onto the
%   Simplex Ball using Dual Function and exact solution by solving linear
%   equation.
% Input:
%   - vY            -   Input Vector.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - ballRadius    -   Ball Radius.
%                       Sets the Radius of the Simplex Ball. For Unit
%                       Simplex set to 1.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: (0, inf).
% Output:
%   - vX            -   Output Vector.
%                       The projection of the Input Vector onto the Simplex
%                       Ball.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  https://math.stackexchange.com/questions/2327504.
% Remarks:
%   1.  The solver finds 2 points which one is positive and the other is
%       negative. Then, since the objective function is linear, finds the
%       exact point where the linear function has value of zero.
% TODO:
%   1.  U.
% Release Notes:
%   -   1.0.000     14/04/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

numElements = size(vY, 1);

if((abs(sum(vY) - ballRadius) < 1e-9) && all(vY >= 0))
    % The input is already within the Simplex.
    vX = vY;
    return;
end

vZ = sort(vY, 'ascend');

vParamMu = [vZ(1) - ballRadius; vZ; vZ(numElements) + ballRadius];
hObjFun = @(paramMu) sum( max(vY - paramMu, 0) ) - ballRadius;

vObjVal = zeros(numElements + 2, 1);
for ii = 1:(numElements + 2)
	vObjVal(ii) = hObjFun(vParamMu(ii));
end

if(any(vObjVal == 0))
    paramMu = vParamMu(vObjVal == 0);
else
    % Working on when an Affine Function have the value zero
    valX1Idx = find(vObjVal > 0, 1, 'last');
    valX2Idx = find(vObjVal < 0, 1, 'first');

    valX1 = vParamMu(valX1Idx);
    valX2 = vParamMu(valX2Idx);
    valY1 = vObjVal(valX1Idx);
    valY2 = vObjVal(valX2Idx);

    paramA = (valY2 - valY1) / (valX2 - valX1);
    paramB = valY1 - (paramA * valX1);
    paramMu = -paramB / paramA;
end

vX = max(vY - paramMu, 0);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
classdef BoundayConditionsPlotter_2D < BoundayConditionsPlotter_Abstract

    methods (Access = public)

        function obj = BoundayConditionsPlotter_2D(axes,mesh,bc)
            obj.init(axes,mesh,bc);
        end

        function plotDirichlet(obj)
            if ~isempty(obj.iD)
                plot(obj.axes,obj.mesh.coord(obj.iD,1),obj.mesh.coord(obj.iD,2),'>','Color',obj.colorD,'MarkerSize',4,'MarkerFaceColor',obj.colorD)
                quiver(obj.axes,obj.mesh.coord(obj.iD,1),obj.mesh.coord(obj.iD,2),obj.vD(:,1),obj.vD(:,2),'Color',obj.colorD,'AutoScaleFactor',obj.scaleD,'LineWidth',obj.lineWidth,'MaxHeadSize',obj.maxHeadSize);
            end
        end

        function plotNeumann(obj)
            if ~isempty(obj.iN)
                plot(obj.axes,obj.mesh.coord(obj.iN,1),obj.mesh.coord(obj.iN,2),'o','MarkerFaceColor',obj.colorN,'Color',obj.colorN)
                quiver(obj.axes,obj.mesh.coord(obj.iN,1),obj.mesh.coord(obj.iN,2),obj.vN(:,1),obj.vN(:,2),'Color',obj.colorN,'AutoScaleFactor',obj.scaleN,'LineWidth',obj.lineWidth,'MaxHeadSize',obj.maxHeadSize);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%create_nullspace Create null-space for scene with perspective camera.
%
%  Parameters:
%    depths .. matrix of zeros and ones meaning whether corresponding element
%              in M is scaled (i.e. multiplied by known perspective depth)

function [nullspace, result] = create_nullspace(M, depths, central, opt)

if nargin < 4, opt.trial_coef = 1;
               opt.threshold  = .01; end

I = ~isnan(M(1:3:end,:)); [m n] = size(I); show_mod = 10; use_maxtuples = 0;
if opt.verbose, fprintf(1, 'Used 4-tuples (.=%d): ', show_mod); tic; end

if central,  cols_scaled(1:n) = 0; cols_scaled(find(I(central,:) > 0)) = 1;
else, cols_scaled = []; end

num_trials = round(opt.trial_coef * n);

if opt.verbose, fprintf(1,'(Allocating memory...'); end
nullspace(size(M,1),num_trials) = 0; % Memory allocation: at least this
                                     % because at least one column is
                                     % added per one 4-tuple.
width                           = 0;
if opt.verbose, fprintf(1, ')'); end

tenth = .1;  % because of the first %
result.used = 0; result.failed = 0;
for i = 1:num_trials
  % choose a 4/max-tuple
  cols = 1:n;
  rows = 1:m;

  cols_chosen = []; t=1; failed = 0;
  if central,
    scaled_ensured = 0;
  else
    scaled_ensured = 1;   % trial version: no scale controling when cutting
  end
  for t = 1:4
    % choose one column, cut useless parts etc.
    [c, cols] = random_element(cols);
    cols_chosen = [cols_chosen c];

    % check just added column
    rows = intersect(rows, find(I(:,c) > 0));

    if t < 4,
      [rows, cols, scaled_ensured] = cut_useless(I, cols_scaled, ...
                                cols_chosen, rows, cols, 4-t, scaled_ensured);
    end

    if isempty(rows), failed = 1; break; end
  end

  if ~failed
    % use the 4/max-tuple
    d = depths(rows,cols_chosen);
    % see ``debug code'' in the comment lower

    rowsbig   = k2i(rows);
    submatrix=[]; for j=1:length(cols_chosen) % 4,
      submatrix=[ submatrix ...
                  spread_depths_col(M(rowsbig,cols_chosen(j)), d(:,j)) ]; end
    debug=1; if debug, if size(submatrix, 1)<=size(submatrix,2) && opt.verbose
        fprintf(1,'-'); end;end
    subnull = nulleps(submatrix,opt.threshold); %svd(submatrix)
    if size(subnull,2)>0  &&  ( use_maxtuples || ...
       size(submatrix,1) == size(submatrix,2) + size(subnull,2))
      nulltemp            = zeros(size(M,1),size(subnull,2));
      nulltemp(rowsbig,:) = subnull; % * (length(rows)/m); % weighting
      if width+size(nulltemp,2) > size(nullspace,2) % Memory allocation:
        if opt.verbose, fprintf(1,'(Allocating memory...)'); end
        mean_added = width/i;
        nullspace(size(M,1), size(nullspace,2) ...
                  + round(mean_added * (num_trials-i))) = 0;
      end
      nullspace(:, width+1 : width+size(nulltemp,2)) = nulltemp;
      width                                          = width +size(nulltemp,2);
      result.used         = result.used +1;
      if mod(result.used, show_mod)==0 && opt.verbose, fprintf(1,'.'); end
    end
  else
    result.failed = result.failed +1;
  end

  if i/num_trials > .1*tenth
    if opt.verbose, fprintf(1,'%d%%', tenth*10); end
    if tenth < 1, tenth=0; end
    tenth = tenth + 1;
  end
end

if size(nullspace,2) > width,   % cut off unused allocated memory
  if opt.verbose, fprintf(1,'(Cutting unused memory...'); end
  nullspace = nullspace(:,1:width);
  if opt.verbose, fprintf(1, ')'); end
end
if opt.verbose, fprintf(1, ['(' num2str(toc) ' sec)\n']); end
result.tried = i;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% debug code
%    if size(d, 1) == 1, bla; end
%    if sum(d(1,:)==0) > 2, blabla; end
%    if use_maxtuples | ...
%          size(d,1)*3 > size(d,2) + sum(d(:)==0)... % i.e. size(submatrix,1)>
%          ...                                       % size(submatrix,2)
%          - sum(cols_scaled(cols_chosen)==0) % i.e. subtract 1 per each
%                                             % column full of zeros
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [rows, cols, scaled_ensured] = cut_useless( ... %)
    I, cols_scaled, ... % this is always same
    cols_chosen, rows, cols, demanded, scaled_ensured)

if ~scaled_ensured
  % check scaled columns
  if length(rows) == 2, demanded_scaled = 3; demanded_rows = 2;
  else                  demanded_scaled = 2; demanded_rows = 3; end
    cols_scaled_chosen = sum(cols_scaled(cols_chosen) > 0);

  % if no unscaled are allowed, they must be all cut
  if demanded == demanded_scaled - cols_scaled_chosen,
    cols           = intersect(cols, find(cols_scaled > 0));
    scaled_ensured = 1;
  end
else demanded_rows = 2; end

% check columns
cols = cols(find(sum(I(rows,cols)) >= demanded_rows));

% check rows
rows = rows(find(sum(I(rows,cols)') >= demanded));

function [element, rest] = random_element(set)
% Take a random element out of a set.
element = set(random_int(1, length(set)));
rest    = setdiff(set, element);

function y = random_int(from,to)
y = floor(from + (1 + to - from)*rand);

function [N,s] = nulleps(M,tol)
% Find the nullspace of M.  This is the regular nullspace, augmented by
% vectors that aren't really in the nullspace, but have low singular
% values associated with them.  tol is the threshold on these singular values.
[u,s,v] = svd(M);
sigsvs = sum(diag(s)>tol);
N = u(:,sigsvs+1:size(u,2));

author: strawlab


text:
function t = uminus(t)
%UMINUS Unary minus (-) for tenmat.
%
%   See also TENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

t.data = -t.data;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
%{
    Copyright (C) Cambridge Electronic Design Limited 2014
    Author: James Thompson
    Web: www.ced.co.uk email: james@ced.co.uk, softhelp@ced.co.uk

    This file is part of CEDS64ML, a MATLAB interface to the SON64 library.

    CEDS64ML 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.

    CEDS64ML 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 CEDS64ML.  If not, see <http://www.gnu.org/licenses/>.
%}

function [ iOk, sTitleOut ] = CEDS64ChanTitle( fhand, iChan, sTitleIn )
%CEDS64CHANTITLE Gets and sets a channel title
%   [ iOk, sTitleOut ] = CEDS64ChanTitle( fhand, iChan {, sTitleIn} )
%   Input
%   fhand - An integer handle to an open file
%   iChan - A channel number
%   sTitleIn - (Optional) A string containing the new title
%   Output
%   iOk - 0 if the title was set correctly otherwise a negative error code
%   sTitleOut - A string containing the old title

if (nargin == 2 || nargin == 3) % always get the old title
    %step 1 find out how big the title is is
    dummystring = blanks(1);
    [iSize] = calllib('ceds64int', 'S64GetChanTitle', fhand, iChan, dummystring, -1);
    %step 2 create a string buffer of the correct size
    stringptr = blanks(iSize+1);
    [iOk, sTitleOut] = calllib('ceds64int', 'S64GetChanTitle', fhand, iChan, stringptr, 0);
else
    iOk = -22;
end

% has there been an error?
if iOk < 0
    return;
end

% if not set the new title if we're given one
if (nargin == 3)
iOk = calllib('ceds64int', 'S64SetChanTitle', fhand, iChan, sTitleIn);
end
end


text:
function x = log_1(x,rnd)
%LOG_         Rigorous calculation of  log(1+x)  for 0<=x<=1 according to rnd
%
%   y = log_1(x)
%
%Internal function
%

% written  12/30/98     S.M. Rump
% modified 08/31/98     S.M. Rump  improved accuracy
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  INTLAB_STDFCTS_LOG = getappdata(0,'INTLAB_STDFCTS_LOG');

  setround(0)
  xs = pow2( floor(x*2^13) , -13 );    % first 13 bits
  log1xs = log(1+xs);                  % 1+xs exactly representable in 14 bits

  setround(rnd)
  d = ( x - xs ) ./ (1+xs);            % 0 <= d < 2^-13

  % log(1+x) = log( (1+xs) * (1+d) ) ,   0 <= err <= d^5/5 < 4.45e-17*d
  if rnd==-1
    log1d = ((( (-d)/4 + 1/3 ).*d - 0.5 ).*d).*d + d;
    x = log1xs + ( log1d + (-INTLAB_STDFCTS_LOG.EPS)*abs(log1xs) );
  else
    log1d = (((( d/5 - .25 ).*d + 1/3 ).*d - 0.5 ).*d).*d + d;
    x = log1xs + ( log1d + INTLAB_STDFCTS_LOG.EPS*abs(log1xs) );
  end

  setround(0)

repo: OpenMAS


text:
function x = r8vec_scale ( s, n, x )

%*****************************************************************************80
%
%% R8VEC_SCALE multiplies an R8VEC by a scale factor.
%
%  Discussion:
%
%    An R8VEC is a vector of R8 values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer S, the scale factor.
%
%    Input, integer N, the length of the vector.
%
%    Input/output, real X(N), the vector to be scaled.
%
  x(1:n) = s * x(1:n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% ASTROTIK by Francesco Santilli
% Compute spherical coordinates of a given position.
%
% Usage: [r,alpha,delta] = spherical(X)
%
% where: X = position [L]
%        r = radius [L]
%        alpha = right ascension [rad]
%        delta = declination [rad]

function [r,alpha,delta] = spherical(X)

    if ~(nargin == 1)
        error('Wrong number of input arguments.')
    end

    D = check(X,1);

    if ~(D==3)
        error('Wrong size of an input argument.')
    end

    xx1 = X(1)^2;
    xx2 = X(2)^2;
    xx3 = X(3)^2;
    r = sqrt(xx1 + xx2 + xx3);
    xy = sqrt(xx1 + xx2);
    alpha = atan2(X(2),X(1));
    delta = atan(X(3)/xy);

end

author: Sable


text:
% Copyright (C) 1993-2013, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for MATLAB (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% We demonstrate computing Harris corner features on an image and an
% image sequence.

% We load an image
im = iread('eiffel2-1.jpg', 'grey');
% and display it
idisp(im)

% Now we compute the coordinates of the 200 strongest Harris corner features
h = icorner(im, 'nfeat', 200);
about h
% we see that the result is a vector of PointFeature objects, which have values
h(1)
% We overlay the features on the image
h.plot();
% We note that the features are located where there is very high image contrast
% such as the edges of trees and the tower

% We load an image sequence
im = iread('~/rvc/bridge-l/*.png', 'roi', [20 750; 20 480]);
% Note that We have used the 'roi' option to chop off the ragged edges
% of these image

% The size of the resulting image is
about im
% which is slightly unusual in that the pixels are uint16 types.  The image
% has three dimensions, the third dimension is the number in the sequence.

% Once again we compute the coordinates of the features
h = icorner(im, 'nfeat', 200);
about h
% and because the input was an image sequence the features are computed for
% every image in the sequence.  The result is a cell array, one cell per image
% in the sequence which contains a vector of corner feature objects
about(h{1})

% We can animate the image sequence and the features
ianimate(im, h, 'fps', 10)
% Note that frame to frame the Harris features tend to stick the same point
% in the world such as cars, trees, road markings.

% Finally we can compute Harris featurse on live video.  You will find the features
% tend to stick to corners of pictures, books, shelves and not so much to smoother
% objects likes hands and faces.

% Create an object which represents a video camera connected to your computer
% If you have a laptop you should see the video recording light come on
camera = VideoCamera(0, 'grey');
% Now create a loop to grab frames from the camera and display them
while true
    % grab the next frame from the attached camera
    im = camera.grab();
    % and display it
    idisp(im);
    % compute Harris corners
    h = icorner(im, 'nfeat', 100);
    % and overlay them
    h.plot();
    drawnow

    % test whether to quit, hit 'q'
    if get(gcf,'CurrentCharacter') == 'q' break; end;
end
clear camera

% See also demo/tracking.

author: petercorke

repo: machinevision-toolbox-matlab