text:
function cummomentplot(obj,tabgrouptag)
    %CUMMOMENTPLOT plot the Cumulative Moment into the specified tab group

    Tags.xs = 'CumMom xs contextmenu';
    Tags.bg = 'CumMom bg contextmenu';
    Tags.line = 'CumMom line contextmenu';
    Tags.ax = 'CumMom axes';

    myTab = findOrCreateTab(obj.fig, obj.fig, tabgrouptag, 'Moment');
    ax=findobj(myTab.Children,'flat','Tag',Tags.ax);
    if isempty(ax)
        ax              = axes(myTab);
        ax.TickDir      = 'out';
        ax.YMinorTick   = 'on';
        ax.Box          = 'on';
        ax.Tag          = Tags.ax;
        ax.UserData     = CumMomentAnalysisWindow(ax);
        grid(ax, 'on');
    end
    analy_win = ax.UserData;

    %% plot the main series
    lineProps.Color         = [0 0 0];
    lineProps.LineWidth     = 2.5;
    lineProps.UIContextMenu = line_plot_context_menu();
    lineProps.MinBigMag     = ZmapGlobal.Data.CatalogOpts.BigEvents.MinMag;
    lineProps.DisplayName   = obj.catalog.Name;

    analy_win.add_series(obj.catalog, 'catalog', lineProps);

    %% plot & synchronize cross sections

    plotted_xs = get(ax.Children,'Tag');
    plotted_xs = plotted_xs(startsWith(plotted_xs, 'Xsection'));

    if ~isempty(plotted_xs)
     existing_xs = obj.XSectionTitles;
        % delete cross sections that shouldn't exist
        todel = plotted_xs(~ismember(plotted_xs,existing_xs));
        analy_win.remove_series(todel);
    end
    %% if necessary, add context menu to figure
    cxs = findobj(obj.fig,'Tag',Tags.xs);
    if isempty(cxs)
        cxs = uicontextmenu(obj.fig,'tag',Tags.xs);
        uimenu(cxs,'Label','Open in new window', 'MenuSelectedFcn', @cb_xstimeplot);
    end

    obj.plot_xsections(@xsplotter, 'Xsection cummomplot');

    ax.YLabel.UIContextMenu = obj.sharedContextMenus.LogLinearYScale;

    %% add context menu to the axes
    cbg=findobj(obj.fig,'Tag',Tags.bg);

    if isempty(cbg)
        cbg=uicontextmenu(obj.fig,'Tag',Tags.bg);
        addLegendToggleContextMenuItem(cbg,'bottom','above');
        uimenu(cbg,'Label','Open in new window','Enable','off','MenuSelectedFcn',@obj.cb_timeplot);
    end

    if isempty(ax.UIContextMenu)
        ax.UIContextMenu=cbg;
    end

    function cln = line_plot_context_menu()
        % set up the context menu for these line plots
        cln=findobj(obj.fig,'Tag',Tags.line);
        if isempty(cln)
            cln=uicontextmenu(obj.fig,'tag',Tags.line);
            uimenu(cln, 'Label', 'start here',              'MenuSelectedFcn', @(~,~)obj.cb_starthere(ax));
            uimenu(cln, 'Label', 'end here',                'MenuSelectedFcn', @(~,~)obj.cb_endhere(ax));
            uimenu(cln, 'Label', 'trim to largest event',   'MenuSelectedFcn', @obj.cb_trim_to_largest);
            uimenu(cln, 'Label', 'Open in new window',      'MenuSelectedFcn', @obj.cb_timeplot);
        end
    end

    function h = xsplotter(xs, xscat)
        xsProps.LineWidth   = 1.5;
        xsProps.DisplayName = xs.Name;
        xsProps.Color       = xs.Color;
        mytag               = ['Xsection cummomplot ' xs.Name];
        h = analy_win.add_series(xscat, mytag, xsProps);

    end

    function cb_xstimeplot(~,~)
        myName = get(gco,'DisplayName');
        idx=strcmp(obj.XSectionTitles,myName);
        ZG=ZmapGlobal.Data;
        ZG.newt2=obj.catalog.subset(obj.CrossSections(idx).inside(obj.catalog));
        ZG.newt2.Name=sprintf('Events within %g km of %s',myName);
        ctp=CumTimePlot(ZG.newt2);
        ctp.plot();
    end
end

author: CelsoReyes


text:
function [E,L,G]=GenPCBasis(S,A)
% this function computes the conditional principal portfolios
% see A. Meucci - "Managing Diversification", Risk Magazine, June 2009
% available at www.ssrn.com

% Code by A. Meucci. This version March 2009.
% Last version available at MATLAB central as "Managing Diversification"

% inputs
% S : covariance matrix
% A : conditioning matrix

% outputs
% E : conditional principal portfolios composition
% L : conditional principal portfolios variances
% G : map weights -> conditional diversification distribution (square root of, not normalized)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(A)
    N=size(S,1);
    K=0;
    [E_,L_]=eig(S);
    E=E_;
    for n=1:N
        E(:,n)=E_(:,N-n+1);
        L(n)=L_(N-n+1,N-n+1);
    end

else

    [K,N]=size(A);
    E=[];
    B=A;
    for n=1:N-K
        if ~isempty(E)
            B=[A
                E'*S];
        end
        e=GenFirstEigVect(S,B);
        E=[E e];
    end

    for n=N-K+1:N
        B=E'*S;
        e=GenFirstEigVect(S,B);
        E=[E e];
    end

    % swap order
    E=[E(:,N-K+1:N) E(:,1:N-K)];
end

L=diag(E'*S*E);

G=diag(sqrt(L))*inv(E);
G=G(K+1:N,:);

author: Sable

repo: mcbench-benchmarks


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypemethod  @@sym {[@var{A}, @var{b}] =} equationsToMatrix (@var{eqns}, @var{vars})
%% @deftypemethodx @@sym {[@var{A}, @var{b}] =} equationsToMatrix (@var{eqns})
%% @deftypemethodx @@sym {[@var{A}, @var{b}] =} equationsToMatrix (@var{eq1}, @var{eq2}, @dots{})
%% @deftypemethodx @@sym {[@var{A}, @var{b}] =} equationsToMatrix (@var{eq1}, @dots{}, @var{v1}, @var{v2}, @dots{})
%% Convert set of linear equations to matrix form.
%%
%% In its simplest form, equations @var{eq1}, @var{eq2}, etc can be
%% passed as inputs:
%% @example
%% @group
%% syms x y z
%% [A, b] = equationsToMatrix (x + y == 1, x - y + 1 == 0)
%%   @result{} A = (sym 2×2 matrix)
%%
%%       ⎡1  1 ⎤
%%       ⎢     ⎥
%%       ⎣1  -1⎦
%%
%%   @result{} b = (sym 2×1 matrix)
%%
%%       ⎡1 ⎤
%%       ⎢  ⎥
%%       ⎣-1⎦
%% @end group
%% @end example
%% In this case, appropriate variables @emph{and their ordering} will be
%% determined automatically using @code{symvar} (@pxref{@@sym/symvar}).
%%
%% In some cases it is important to specify the variables as additional
%% inputs @var{v1}, @var{v2}, etc:
%% @example
%% @group
%% syms a
%% [A, b] = equationsToMatrix (a*x + y == 1, y - x == a)
%%   @print{} ??? ... nonlinear...
%%
%% [A, b] = equationsToMatrix (a*x + y == 1, y - x == a, x, y)
%%   @result{} A = (sym 2×2 matrix)
%%
%%       ⎡a   1⎤
%%       ⎢     ⎥
%%       ⎣-1  1⎦
%%
%%   @result{} b = (sym 2×1 matrix)
%%
%%       ⎡1⎤
%%       ⎢ ⎥
%%       ⎣a⎦
%% @end group
%% @end example
%%
%% The equations and variables can also be passed as vectors @var{eqns}
%% and @var{vars}:
%% @example
%% @group
%% eqns = [x + y - 2*z == 0, x + y + z == 1, 2*y - z + 5 == 0];
%% [A, B] = equationsToMatrix (eqns, [x y])
%%   @result{} A = (sym 3×2 matrix)
%%
%%       ⎡1  1⎤
%%       ⎢    ⎥
%%       ⎢1  1⎥
%%       ⎢    ⎥
%%       ⎣0  2⎦
%%
%%   B = (sym 3×1 matrix)
%%
%%       ⎡ 2⋅z ⎤
%%       ⎢     ⎥
%%       ⎢1 - z⎥
%%       ⎢     ⎥
%%       ⎣z - 5⎦
%% @end group
%% @end example
%% @seealso{@@sym/solve}
%% @end deftypemethod

function [A, b] = equationsToMatrix(varargin)

  % when Symbols are specified, this won't be used
  s = findsymbols (varargin);

  cmd = {'L, symvars = _ins'
         'if not isinstance(L[-1], MatrixBase):'
         '    if isinstance(L[-1], Symbol):'  % Symbol given, fill vars...
         '        vars = list()'
         '        for i in reversed(range(len(L))):'
         '            if isinstance(L[i], Symbol):'
         '                vars = [L.pop(i)] + vars'
         '            else:'  % ... until we find a non-Symbol
         '                break'
         '    else:'
         '        vars = symvars'
         'else:'
         '    if len(L) == 1:'  % we have only a list of equations
         '        vars = symvars'
         '    else:'
         '        vars = L.pop(-1)'
         'if len(L) == 1:'  % might be matrix of eqns, don't want [Matrix]
         '    L = L[0]'
         'vars = list(vars)'
         'A, B = linear_eq_to_matrix(L, vars)'
         'return True, A, B' };

  for i = 1:length(varargin)
    varargin{i} = sym (varargin{i});
  end

  [s, A, b] = pycall_sympy__ (cmd, varargin, s);

  if ~s
    error('Cannot convert to matrix; system may be nonlinear.');
  end

end

%!test
%! syms x y z
%! [A, B] = equationsToMatrix ([x + y - z == 1, 3*x - 2*y + z == 3, 4*x - 2*y + z + 9 == 0], [x, y, z]);
%! a = sym ([1 1 -1; 3 -2 1; 4 -2 1]);
%! b = sym ([1; 3; -9]);
%! assert (isequal (A, a))
%! assert (isequal (B, b))

%!test
%! syms x y z
%! A = equationsToMatrix ([3*x + -3*y - 5*z == 9, 4*x - 7*y + -3*z == -1, 4*x - 9*y - 3*z + 2 == 0], [x, y, z]);
%! a = sym ([3 -3 -5; 4 -7 -3; 4 -9 -3]);
%! assert (isequal (A, a))

%!test
%! syms x y
%! [A, B] = equationsToMatrix ([3*x + 9*y - 5 == 0, -8*x - 3*y == -2]);
%! a = sym ([3 9; -8 -3]);
%! b = sym ([5; -2]);
%! assert (isequal (A, a))
%! assert (isequal (B, b))

%!test
%! % override symvar order
%! syms x y
%! [A, B] = equationsToMatrix ([3*x + 9*y - 5 == 0, -8*x - 3*y == -2], [y x]);
%! a = sym ([9 3; -3 -8]);
%! b = sym ([5; -2]);
%! assert (isequal (A, a))
%! assert (isequal (B, b))

%!test
%! syms x y z
%! [A, B] = equationsToMatrix ([x - 9*y + z == -5, -9*y*z == -5], [y, x]);
%! a = sym ([[-9 1]; -9*z 0]);
%! b = sym ([-5 - z; -5]);
%! assert (isequal (A, a))
%! assert (isequal (B, b))

%!test
%! syms x y
%! [A, B] = equationsToMatrix (-6*x + 4*y == 5, 4*x - 4*y - 5, x, y);
%! a = sym ([-6 4; 4 -4]);
%! b = sym ([5; 5]);
%! assert (isequal (A, a))
%! assert (isequal (B, b))

%!test
%! % vertical list of equations
%! syms x y
%! [A, B] = equationsToMatrix ([-6*x + 4*y == 5; 4*x - 4*y - 5], [x y]);
%! a = sym ([-6 4; 4 -4]);
%! b = sym ([5; 5]);
%! assert (isequal (A, a))
%! assert (isequal (B, b))

%!test
%! syms x y
%! [A, B] = equationsToMatrix (5*x == 1, y, x - 6*y - 7, y);
%! a = sym ([0; 1; -6]);
%! b = sym ([1 - 5*x; 0; -x + 7]);
%! assert (isequal (A, a))
%! assert (isequal (B, b))

%!error <nonlinear>
%! syms x y
%! [A, B] = equationsToMatrix (x^2 + y^2 == 1, x - y + 1, x, y);

%!test
%! % single equation
%! syms x
%! [A, B] = equationsToMatrix (3*x == 2, x);
%! a = sym (3);
%! b = sym (2);
%! assert (isequal (A, a))
%! assert (isequal (B, b))

%!test
%! % single equation w/ symvar
%! syms x
%! [A, B] = equationsToMatrix (3*x == 2);
%! a = sym (3);
%! b = sym (2);
%! assert (isequal (A, a))
%! assert (isequal (B, b))

%!error <unique>
%! syms x
%! equationsToMatrix (3*x == 2, [x x])

repo: octsympy


text:
function DCM = spm_dcm_diagnose(DCM,varargin)
% Post hoc diagnosis of DCMs (under the Laplace approximation)
% FORMAT spm_dcm_diagnose(DCM,'field','field',...)
%
% DCM        - DCM stricture (inverted)
% field      - field name(s) of parameters to consider
%
%--------------------------------------------------------------------------
% This routine searches over all possible reduced models of a full model
% (DCM) and uses post hoc model selection to select the best. Reduced
% models mean all permutations of free parameter sets (parameters with non-
% zero prior covariance), where models are defined in terms of their prior
% covariance. The full model should be inverted prior to post hoc
% optimization. If there are more than 8 free-parameter sets, this routine
% will implement a greedy search: This entails searching over all
% permutations of the 8 sets whose removal (shrinking the prior
% variance to zero) produces the smallest reduction (greatest increase)
% in model evidence. This procedure is repeated until all sets
% are retained in the best model or there are no more parameters to
% consider.
%
% A parameter set is specified implicitly by the structure (DCM.Ep). Each
% set corresponds to a column of (the cell arrays or matrix) each field of
% DCM.Ep.
%
% if only one field is specified the log-evidence is computed as a function
% of the scaled prior variance. Redundant parameters have a log-evidence
% that keeps increasing as the prior variance shrinks.
%
% The outputs of this routine are graphics reporting the model reduction
% (optimisation). Red means weak evidence; blue strong evidence (> 3) and
% cyan very strong evidence (> 5)
%__________________________________________________________________________
% Copyright (C) 2008-2011 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_diagnose.m 7679 2019-10-24 15:54:07Z spm $

%-Get priors and posteriors
%==========================================================================
qE      = DCM.Ep;
qC      = DCM.Cp;
pE      = DCM.M.pE;
pC      = DCM.M.pC;

if isstruct(pC);
    pV  = pC;
    pC  = diag(spm_vec(pC));
end

% find the partition of parameters
%--------------------------------------------------------------------------
str   = {};
ind   = {};
if isstruct(qE)

    % get fields (parameters) to consider
    %----------------------------------------------------------------------
    if ~isempty(varargin)
        field = varargin;
    else
        field = fieldnames(qE);
    end
    for i = 1:length(field)
        q = getfield(qE,field{i});
        k = spm_fieldindices(qE,field{i});
        for j = 1:size(q,2)
            u            = 1:length(spm_vec(q(:,j)));
            str{end + 1} = [field{i} '(' num2str(j) ')'];
            ind{end + 1} = k(u); k(u) = [];
        end
    end

else
    field = {};
    k     = 1:length(spm_vec(qE));
    for j = 1:size(qE,2)
        u            = 1:length(spm_vec(qE(:,j)));
        str{end + 1} = ['qE(' num2str(j) ')'];
        ind{end + 1} = k(u); k(u) = [];
    end
end

% find free parameters
%--------------------------------------------------------------------------
ip    = {};
sp    = {};
for i = 1:length(str)
    j = find(diag(pC(ind{i},ind{i})));
    if ~isempty(j)
        ip{end + 1} = ind{i}(j);
        sp{end + 1} = str{i};
    end
end

% model search over shrinking priors
% -------------------------------------------------------------------------
n     = length(pC);
if length(field) == 1

    r     = (0:64)/32;
    for j = 1:length(r)
        for i = 1:length(ip)
            R      = speye(n,n) - sparse(ip{i},ip{i},1,n,n)*(1 - r(j));
            Z(i,j) = spm_log_evidence(qE,qC,pE,pC,pE,R*pC*R);

        end
    end

    subplot(2,1,1)
    plot(r,Z)
    title(['log-evidence for ' field{1}],'FontSize',16);
    xlabel('proportion of full prior')
    ylabel('negative free-energy')
    return

end

% model search over new prior without the i-th parameter set
% -------------------------------------------------------------------------
for i = 1:length(ip)
    R    = speye(n,n) - sparse(ip{i},ip{i},1,n,n);
    Z(i) = spm_log_evidence(qE,qC,pE,pC,pE,R*pC*R);
end

% find parameters with the least evidence
%--------------------------------------------------------------------------
[Z i] = sort(-Z);
ip    = ip(i);
sp    = sp(i);

% If there are too many subsets search use those with the least evidence
%--------------------------------------------------------------------------
if length(ip) > 8

    % flag a greedy search
    %----------------------------------------------------------------------
    ip     = ip(1:8);
    repeat = 1;

else
    repeat = 0;
end

% Create model space in terms of free parameter indices
%--------------------------------------------------------------------------
K     = spm_perm_mtx(length(ip));

% and get log-evidences
%==========================================================================
for i = 1:size(K,1)
    k    = spm_vec(ip(find(K(i,:))));
    R    = speye(n,n) - sparse(k,k,1,n,n);
    S(i) = spm_log_evidence(qE,qC,pE,pC,pE,R*pC*R);
end

% Find a good model with the most parameters
% =========================================================================
[i k]  = max(S);
k      = find(K(k,:));

nelim  = length(k);
repeat = repeat & nelim;

% Continue greedy search if any parameters have been eliminated
%--------------------------------------------------------------------------
if repeat

    fprintf('removing %i parameter set: ',nelim)
    fprintf('%s, ',sp{k});fprintf('\n')

    % Optimise selected model (prior covariance)
    %----------------------------------------------------------------------
    k         = spm_vec(ip(k));
    R         = speye(n,n) - sparse(k,k,1,n,n);
    rC        = R*pC*R;

    % Get posterior of selected model (rC)
    % ---------------------------------------------------------------------
    [F Ep Cp] = spm_log_evidence(qE,qC,pE,pC,pE,rC);

    DCM.M.pC  = rC;
    DCM.Ep    = Ep;
    DCM.Cp    = Cp;
    DCM.F     = F;
    DCM       = spm_dcm_diagnose(DCM,field{:});
    return

else

    % model posterior
    % ---------------------------------------------------------------------
    S     = S - S(end);
    p     = exp(S - max(S));
    p     = p/sum(p);

    fprintf('\n%i Irrelevant parameter sets: ',length(k))
    fprintf('%s, ',sp{k});fprintf('\n')

    % Organise field names and show results
    % ---------------------------------------------------------------------
    spm_figure('Getwin','Graphics');

    for i = 1:4
        j = (1:8) + (i - 1)*8;
        try
            tsr{i,1} = sprintf('%s,',sp{j});
        catch
            tsr{i,1} = sprintf('%s,',sp{j(1):end});
        end
    end

    subplot(2,1,1)
    bar(Z,'c'),hold on
    bar(Z(Z < 5),'b'),hold on
    bar(Z(Z < 3),'r'),hold off
    title('relative log evidence (with vs. without)','FontSize',16)
    xlabel(tsr,'FontSize',16)
    ylabel('negative free-energy','FontSize',12)

    subplot(2,2,3)
    if size(K,1) > 32, plot(S,'k'), else, bar(S,'c'), end
    title('log-posterior','FontSize',16)
    xlabel('model','FontSize',12)
    ylabel('log-probability','FontSize',12)
    axis square

    subplot(2,2,4)
    if size(K,1) > 32, plot(p,'k'), else, bar(p,'r'), end
    title('model posterior','FontSize',16)
    xlabel('model','FontSize',12)
    ylabel('probability','FontSize',12)
    axis square

end

author: spm


text:
function [tap] = hanning(n, str)

%HANNING   Hanning window.
%   HANNING(N) returns the N-point symmetric Hanning window in a column
%   vector.  Note that the first and last zero-weighted window samples
%   are not included.
%
%   HANNING(N,'symmetric') returns the same result as HANNING(N).
%
%   HANNING(N,'periodic') returns the N-point periodic Hanning window,
%   and includes the first zero-weighted window sample.
%
%   NOTE: Use the HANN function to get a Hanning window which has the
%          first and last zero-weighted samples.
%
%   See also BARTLETT, BLACKMAN, BOXCAR, CHEBWIN, HAMMING, HANN, KAISER
%   and TRIANG.
%
%   This is a drop-in replacement to bypass the signal processing toolbox

% Copyright (c) 2010, Jan-Mathijs Schoffelen, DCCN Nijmegen
%
% 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$

if nargin==1
  str = 'symmetric';
end

switch str
case 'periodic'
   % Includes the first zero sample
   tap = [0; hanningX(n-1)];
case 'symmetric'
   % Does not include the first and last zero sample
   tap = hanningX(n);
end

function tap = hanningX(n)

% compute taper
N   = n+1;
tap = 0.5*(1-cos((2*pi*(1:n))./N))';

% make symmetric
halfn = floor(n/2);
tap( (n+1-halfn):n ) = flipud(tap(1:halfn));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% this script plots the tire figure
alpha = -0.2:0.001:0.2;

figure; grid on; hold on; box on;
PacParam = [40,   1, 4500,  0];
PacParam2 = [25,   40/25, 4500,  0];
FyF_linear = PacParam(1)*PacParam(2)*PacParam(3).*alpha;
FyF_nonlinear1 = PacParam(3).*sin(PacParam(2).*atan(PacParam(1).*alpha - PacParam(4).*(PacParam(1).*alpha - atan(PacParam(1).*alpha))));
FyF_nonlinear2 = PacParam2(3).*sin(PacParam2(2).*atan(PacParam2(1).*alpha - PacParam2(4).*(PacParam2(1).*alpha - atan(PacParam2(1).*alpha))));
plot(alpha, FyF_linear);
plot(alpha, FyF_nonlinear1);
plot(alpha, FyF_nonlinear2);
xlabel('Side slip angle $\alpha$ in rad');
ylabel('Force $F_y$ in N');
ylim([-10000, 10000]);
legend('Linear', 'Pacejka 1', 'Pacejka 2');
matlab2tikz('TireModelComparison.tex', 'standalone', true);

author: TUMFTM


text:
%
% SCRIPT TO TEST THE KINEMATICS OF THE 3RRR robot
%
% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE 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 ARTE.  If not, see <http://www.gnu.org/licenses/>.

close all

%load the robot
%the robot is represented (interiorly) by three planar 2dof arms
robot=load_robot('example','3RRR');

adjust_view(robot)

%init T
T=eye(4);
%Represent a final point in space that will be achieved by the end effector
%P=[xA, yA]=[1.5 1.5]
T(1,4)=1.2;
T(2,4)=1.0;
%orientation
phi=pi/3;
T(1,1)=cos(phi)
T(2,1)=sin(phi)

%find 8 solutions for the INVERSE KINEMATIC in position
q=inversekinematic(robot, T)

%plot every solution
drawrobot3d(robot, q(:,1)), pause(2);
drawrobot3d(robot, q(:,2)), pause(2);
drawrobot3d(robot, q(:,3)), pause(2);
drawrobot3d(robot, q(:,4)), pause(2);
drawrobot3d(robot, q(:,5)), pause(2);
drawrobot3d(robot, q(:,6)), pause(2);
drawrobot3d(robot, q(:,7)), pause(2);
drawrobot3d(robot, q(:,8)), pause(2);

%Make a line in space with a variation in the angle phi
x=linspace(1,1.4,50);
y=linspace(1.1,1.15,50);
th=linspace(0,pi/3,50);
Q=[];
for i=1:length(x),
    T(1,4)=x(i);
    T(2,4)=y(i);
    %orientation

    T(1,1)=cos(th(i));
    T(2,1)=sin(th(i))

    q=inversekinematic(robot, T);

    drawrobot3d(robot, q(:,1), 1) %noclear activated
    drawrobot3d(robot, q(:,1)) %noclear deactivated
    pause(0.1)
end

%TEST NOW the directkinematic function
T(1,4)=1.0;
T(2,4)=1.0;
%orientation
phi=pi/3;
T(1,1)=cos(phi)
T(2,1)=sin(phi)

%find 8 solutions for the INVERSE KINEMATIC in position
q=inversekinematic(robot, T)

%now, for the first solution
T=directkinematic(robot, [q(1,1) q(3,1) q(5,1)])
%caution: there may exist less than 8 possible solutions
for i=1:8,
    drawrobot3d(robot, q(:,i))
    % 1 and 3 correspond to the joint variables of the robot
    T=directkinematic(robot, [q(1,i) q(3,i) q(5,i)])
    pause(2);
    close all
end

repo: ARTE


text:
function y = Fault_decision_wr_t2(u)
global st16

% Using model defined by the structure st16 to make decision about type2 wr fault

sker=st16.x2sup+(abs(u))'*abs(u)*ones(st16.Nlsup,1)-2*st16.xsup*abs(u);
y=(st16.w)'*exp(-sker./(2*(st16.sigma).^2))+st16.b;

if y>=-0.3
    y=1;
else
    y=0;
end
end


text:
function [im2, H] = repeatability_saltnpepper(im, n_pixels)
  seed = sum(im(:));
  s = RandStream('mt19937ar', 'Seed', seed);

  [h,w,~] = size(im);
  num_missing = n_pixels;
  coords = [];
  while num_missing > 0
    x = randi(s, [1 w], num_missing, 1);
    y = randi(s, [1 h], num_missing, 1);
    coords = unique([coords; [x y]], 'rows');
    num_missing = n_pixels - size(coords, 1);
  end

  im2 = im;
  H = eye(3);
  if isempty(coords)
    return;
  end

  im_gray = rgb2gray(im);
  coord_inds = sub2ind(size(im_gray), coords(:,2), coords(:,1));
  is_dark = im_gray(coord_inds) < 128;
  chan_stride = h * w;
  for ch = 1:3
    % dark goes bright
    im2(coord_inds(is_dark) + (ch - 1) * chan_stride) = 255;
    % bright goes dark
    im2(coord_inds(~is_dark) + (ch - 1) * chan_stride) = 0;
  end
end

author: hosang

repo: detection-proposals


text:
function [ i, j, v ] = find( x )

%Disciplined convex/geometric programming information for FIND:
%   When used in CVX models, FIND cannot deduce the numeric value of a
%   CVX variable. So it returns the indices of elements that are
%   "structurally" nonzero; i.e., that are not *identically* zero.
%   This is similar to the concept of structural nonzeros in sparse
%   matrix factorizations. To illustrate the distinction, consider:
%      variable x(3);
%      x(2) == 0;
%      y = [x(1);0;x(3)];
%   In this case, FIND(x) returns [1;2;3] even though x(2) has been set
%   to zero in an equality constraint. (After all, the overall model may
%   be infeasible, in which case x(2) is arguably not zero.) However,
%   FIND(y) will return [1;3], because y(2) is identically zero.
%
%   When X is a CVX variable, the first two outputs of [I,J,V]=FIND(X)
%   are constant, and the third is a CVX variable containing the
%   structural nonzeros of X.
%
%   FIND(X) places no convexity restrictions on its argument.

ndxs = find( any( x.basis_, 1 ) );
ndxs = ndxs( : );

if nargout > 1,
    i = ndxs - 1;
    j = floor( i / x.size_(1) ) + 1;
    i = rem( i, x.size_(1) ) + 1;
else
	i = ndxs;
end

if nargout > 2,
    v = reshape( cvx_subsref( x, ndxs ), length( ndxs ) );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function y = perform_thresholding(x, t, type)

% perform_thresholding - perform hard or soft thresholding
%
%   y = perform_thresholding(x, t, type);
%
%   type is either 'hard' or 'soft' or 'semisoft'
%   t is the threshold
%
%   works also for complex data, and for cell arrays.
%
%   if type is 'strict' then it keeps the t largest entry in each
%   column of x.
%
%   Copyright (c) 2006 Gabriel Peyre

if nargin<3
    type = 'hard';
end

if iscell(x)
    % for cell arrays
    for i=1:size(x,1)
        for j=1:size(x,2)
            y{i,j} = perform_thresholding(x{i,j},t, type);
        end
    end
    return;
end

switch lower(type)
    case {'hard', ''}
        y = perform_hard_thresholding(x,t);
    case 'soft'
        y = perform_soft_thresholding(x,t);
    case 'semisoft'
        y = perform_semisoft_thresholding(x,t);
    case 'strict'
        y = perform_strict_thresholding(x,t);
    otherwise
        error('Unkwnown thresholding type.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = perform_strict_thresholding(X,s)
%% keep only the s largest coefficients in each column of X

v = sort(abs(X)); v = v(end:-1:1,:);
v = v(round(s),:);
v = repmat(v, [size(X,1) 1]);
X = X .* (abs(X)>=v);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = perform_hard_thresholding(x,t)

t = t(1);
y = x .* (abs(x) > t);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = perform_soft_thresholding(x,t)

if not(isreal(x))
    % complex threshold
    d = abs(x);
    d(d<eps) = 1;
    x = x./d .* perform_soft_thresholding(d,t);
end

t = t(1);
s = abs(x) - t;
s = (s + abs(s))/2;
y = sign(x).*s;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = perform_semisoft_thresholding(x,t)

if length(t)==1
    t = [t 2*t];
end
t = sort(t);

y = x;
y(abs(x)<t(1)) = 0;
I = find(abs(x)>=t(1) & abs(x)<t(2));
y( I ) = sign(x(I)) .* t(2)/(t(2)-t(1)) .* (abs(x(I))-t(1));

author: gpeyre

repo: matlab-toolboxes


text:
%%***********************************************************************
%% fapproblem:
%%
%% (primal prob.) min  Tr C*X
%%                s.t. diag(X) = e
%%                    -Eij * X = 2/(k-1)  if (i,j) in  U
%%                    -Eij * X <= 2/(k-1) if (i,j) in GU
%%  where Eij = ei*ej' + ej*ei'
%%  b = e
%%  C = (1/2)*Diag(We) - (k-1)/(2k) *L(G)
%%  SignP(i,j) = 1 if (i,j) in GU, 2 if (i,j) in U; 0 otherwise
%%  M_ij = -1/(k-1) if (i,j) in  U or GU; 0 otherwise
%%-----------------------------------------------------------------------
%% (primal prob.) min  Tr C*X
%%                s.t. diag(X) = e
%%                    LL <= X <= UU  if (i,j) in  U
%%
%%  where L_ij = -1/(k-1) if (i,j) in  U or GU; -inf otherwise,
%%      U_ij = -1/(k-1) if (i,j) in  U; inf otherwise
%%
%%
%%
%% [blk,Avec,C,b,LL,UU,SignP,M,kparm] = fapread_lu(fname);
%%***********************************************************************
%% SDPNAL+
%% Copyright (c) 2014 by
%% Liuqin Yang, Defeng Sun, and Kim-Chuan Toh
%%***********************************************************************

 function [blk,Avec,C,b,LL,UU,SignP,M,kparm] = fapread_lu(fname)

%%
%% read fap data
%%
   if exist(fname)
      fid = fopen(fname,'r');
   elseif exist([fname,'.dat']);
      fid = fopen([fname,'.dat'],'r');
   else
      fprintf('** Problem not found. \n');
      blk = []; Avec = []; C = []; b = [];
      LL =[]; UU =[];SignP = []; M =[]; kparm = [];
      return;
   end
   [tmpr,count] = fscanf(fid,'%c');
   datavec = sscanf(tmpr,'%f'); clear tmpr;
   n = datavec(1);
   numedges = datavec(2);
   kparm = datavec(3);
   datavec = datavec(4:length(datavec));
   len = length(datavec);
   if (len ~= 3*numedges)
      error(' fapread: numedges and data do not match.');
   end
   I = datavec(1:3:len);
   J = datavec(2:3:len);
   w = datavec(3:3:len);

   idxU = find(w==1000);
   IU = I(idxU); JU = J(idxU); wU = w(idxU);
   U = spconvert([IU JU wU; n n 0]);
   U = U + U';

   idx2 = find(w~=1000);
   I2 = I(idx2); J2 = J(idx2); w2 = w(idx2);
   GU = spconvert([I2 J2 w2; n n 0]);
   GU = GU + GU';

   fclose(fid);
%%
%% blk, Avec, C, b
%%
    n = length(U);
    mU = nnz(triu(U,1));
    mGU = nnz(triu(GU,1));
    m = mGU + mU + n;

    %%b = [ones(n,1); (2/(kparm-1))*ones(mU+mGU,1)];
    b = ones(n,1);
    %b = [];
    LG = diag(GU*ones(n,1))-GU;
    C{1,1} = 0.5*(diag(GU*ones(n,1))) - ((kparm-1)/(2*kparm))*LG;
    %%C{2,1} = zeros(mGU,1);

    blk{1,1} = 's';  blk{1,2} = n;
    %%blk{2,1} = 'l';  blk{2,2} = mGU;
%%
%%
%%
    r2 = sqrt(2);
    I = zeros(m,1); J = zeros(m,1); w = zeros(m,1);
    cnt = 0;
    e = [1:n]';
    I(1:n) = e.*(e+1)/2;
    J(1:n) = e;
    w = ones(n,1);
    cnt = cnt+n;
    for i = 1:n
        idx = find(U(i,i+1:n));
        idx = idx+i;      %% adjust index.
        len = length(idx);
        I(cnt+[1:len]) = i + idx.*(idx-1)/2;
        J(cnt+[1:len]) = cnt+[1:len]';
        w(cnt+[1:len]) = -r2*ones(len,1);
        cnt = cnt + len;
    end
    for i = 1:n
        idx = find(GU(i,i+1:n));
        idx = idx+i;      %% adjust index.
        len = length(idx);
        I(cnt+[1:len]) = i + idx.*(idx-1)/2;
        J(cnt+[1:len]) = cnt+[1:len]';
        w(cnt+[1:len]) = -r2*ones(len,1);
        cnt = cnt + len;
    end
    Av = spconvert([I,J,w; n*(n+1)/2, m, 0]);
    Avec{1,1} = Av(:,1:n);

    sign0 = zeros(n,n);
    sign0(GU ~= 0) = 1;
    sign0(U ~= 0) = 2;
    SignP{1,1} = sign0;

    M0 = zeros(n,n);
    M0(sign0 ~= 0) = -1/(kparm-1);
    M{1,1} = M0;
    LL0 = -inf*ones(n,n);
    LL0(sign0 ~= 0) = -1/(kparm-1);
    LL{1,1} = LL0;
    UU0 = inf*ones(n,n);
    UU0(sign0 == 2) = -1/(kparm-1);
    UU{1,1} = UU0;
%%***********************************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function Test_fixedrank_3factors()
% function Test_LSR()
%
% Test file for the fixedrankLSRquotientfactory
% Problem: Low-rank matrix completion
% Paper link: http://arxiv.org/abs/1112.2318
%
% All intputs are optional.
%

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, July 11, 2013.
% Contributors: Nicolas Boumal
% Change log:

    clc; %close all;

    % Define the problem data
    m = 100;
    n = 100;
    r = 5;
    A = randn(m, r);
    S = randn(n, r);
    Xopt = A*S'; % Original low-rank matrix

    % Create the problem structure
    % quotient LSR geometry

%         M = fixedrankfactory_3factors(m, n, r);
    M = fixedrankfactory_3factors_preconditioned(m, n, r);

    problem.M = M;

    % Which entries do we observe ?
    df = problem.M.dim();
    p = 3*df/(m*n);
    mask = rand(m, n) <= p;
    % mask = ones(m, n);

    % Define the problem cost function
    problem.cost = @cost;
    function f = cost(X)
        f = .5 * norm(mask.*((X.L*X.S*X.R') - Xopt ), 'fro')^2;
    end

    problem.grad = @(X) problem.M.egrad2rgrad(X, egrad(X));
    function g = egrad(X)
        P = (mask.^2) .*(X.L*X.S*X.R' - Xopt);
        g.L= P*X.R*X.S';
        g.S = X.L'*P*X.R;
        g.R = P'*X.L*X.S;
    end

    problem.hess = @(X, L) problem.M.ehess2rhess(X, egrad(X), ehess(X, L), L);
    function Hess = ehess(X, eta)
        P = (mask.^2) .* ( X.L*X.S*X.R' - Xopt );
        Pdot  = (mask.^2).*(eta.L*X.S*X.R' + X.L*eta.S*X.R' + X.L*X.S*eta.R');

        Hess.L = P*(X.R*eta.S' + eta.R*X.S') + Pdot*X.R*X.S';
        Hess.R = P'*(X.L*eta.S + eta.L*X.S) + Pdot'*X.L*X.S;
        Hess.S = X.L'*Pdot*X.R + eta.L'*P*X.R + X.L'*P*eta.R;

    end

    % % Check numerically whether gradient and Hessian are correct
    checkgradient(problem);
    drawnow;
    pause;
    checkhessian(problem);
    drawnow;
    pause;

    % Initial guess
    X0 = [];

    % Options (not mandatory)
    options.maxiter = inf;
    options.maxinner = 30;
    options.maxtime = 120;
    options.tolgradnorm = 1e-5;
    options.Delta_bar = min(m, n)*r;
    options.Delta0 =  min(m, n)*r/64;

    % Pick an algorithm to solve the problem
    [Xsol costopt info] = trustregions(problem, X0, options);
    % [Xsol costopt info] = steepestdescent(problem, X0, options);
    % [Xsol costopt info] = conjugategradient(problem, X0, options);

    evs = hessianspectrum(problem, Xsol);
    evs = real(evs);
    max(evs)/min(evs)
    stairs(sort(evs));
    title(['Eigenvalues of the Hessian of the cost function ' ...
        'at the solution']);

end

author: NicolasBoumal

repo: manopt


text:
function fhat = convSO3(fhat1,fhat2)
%

% old sizes
s1 = size(fhat1);
s2 = size(fhat2);

% get bandwidth
L = min(dim2deg(s1(1)),dim2deg(s2(1)));

% new size
l=length(s2)-length(s1);
s = max([s1(2:end),ones(1,l);s2(2:end),ones(1,-l)]);

% compute Fourier coefficients of the convolution
fhat = zeros([deg2dim(L+1),s]);
if prod(s) == 1 %simple SO3Fun
  for l = 0:L
    ind = deg2dim(l)+1:deg2dim(l+1);
    fhat(ind) = reshape(fhat2(ind),2*l+1,2*l+1) * ...
      reshape(fhat1(ind),2*l+1,2*l+1) ./ sqrt(2*l+1);
  end
else % vector valued SO3Fun
  for l = 0:L
    ind = deg2dim(l)+1:deg2dim(l+1);
    fhat_l = pagemtimes( full(reshape(fhat2(ind,:),[2*l+1,2*l+1,s2(2:end)])) , ...
      full(reshape(fhat1(ind,:),[2*l+1,2*l+1,s1(2:end)])) ) ./ sqrt(2*l+1);
    fhat(ind,:) = reshape(fhat_l,[],prod(s));
  end
end

end

author: mtex-toolbox

repo: mtex

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


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

%*****************************************************************************80
%
%% F04_F0 returns the value of function 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N,1), Y(N,1), the evalution points.
%
%    Output, real F(N,1), the function values.
%
  f(1:n,1) = exp ( - 5.0625 * ( ( x(1:n,1) - 0.5 ).^2 ...
                              + ( y(1:n,1) - 0.5 ).^2 ) ) / 3.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function spline_test ( )

%*****************************************************************************80
%
%% SPLINE_TEST tests the SPLINE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 June 2013
%
%  Author
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPLINE_TEST\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Test the SPLINE library.\n' );

  spline_test001 ( );
  spline_test002 ( );
  spline_test003 ( );
  spline_test004 ( );
  spline_test005 ( );
  spline_test006 ( );

  spline_test01 ( );
  spline_test02 ( );
  spline_test03 ( );
  spline_test04 ( );
  spline_test05 ( );
  spline_test06 ( );
  spline_test07 ( );
  spline_test08 ( );
  spline_test09 ( );

  spline_test10 ( );
  spline_test11 ( );
  spline_test115 ( );
  spline_test116 ( );
  spline_test12 ( );
  spline_test125 ( );
  spline_test126 ( );
  spline_test127 ( );
  spline_test13 ( );
  spline_test14 ( );
  spline_test145 ( );
  spline_test15 ( );
  spline_test16 ( );
  spline_test17 ( );
  spline_test18 ( );
  spline_test19 ( );
  spline_test195 ( );

  spline_test20 ( );
  spline_test205 ( );
  spline_test21 ( );
  spline_test215 ( );
  spline_test22 ( );
  spline_test225 ( );
  spline_test23 ( );
  spline_test235 ( );
  spline_test24 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPLINE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% Test for @chebfun/comet.m.

function pass = test_comet(pref)

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

% Real scalar functions.
fsr1 = chebfun(@sin, [-1 0 1], pref);
fsr2 = chebfun(@cos, [-1 0.5 1], pref);

% Obviously, we can't check if the plots are correct without human intervention,
% so all these tests are meant to do is make sure that we don't crash.

hfig = figure('Visible', 'off');

% Check COMET and COMET3.
pass(1) = doesNotCrash(@() comet(fsr1));
pass(2) = doesNotCrash(@() comet(fsr1, fsr2));

close(hfig);

end

function pass = doesNotCrash(fn)

try
    fn();
    pass = true;
catch ME;
    pass = false;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function e = rmModelSearchFit_twoGaussiansDoGSigmasOnly(p,Y,XvYv,stim,t)
% rmModelSearchFit_twoGaussians - actual fit function of rmSearchFit
%
% error = rmModelSearchFit_twoGaussians(p,Y,trends,Gx,Gy,stim,rawrss);
%
% Basic barebones fit of a single time-series. Error is returned in
% percentage: 100% is RSS of unfitted time-series. This way we can quantify
% the improvement of the fit independend of the variation in the raw
% time-series.
%
% 2006/06 SOD: wrote it.
% 2006/12 SOD: modifications for fmincon, this is litterally called >>10000
% times so we cut every corner possible.

% make RF (taken from rfGaussian2d)
p = -2.*(p.^2);
RF = zeros(numel(XvYv),2);
RF(:,1) = exp( (XvYv) ./ p(1) );
RF(:,2) = exp( (XvYv) ./ p(2) );

% make prediction (taken from rfMakePrediction)
X = [stim * RF t];

% fit
%b = pinv(X)*Y;
[U,S,V] = svd(X,0);
s = diag(S);
tol = numel(X) * eps(max(s));
r = sum(s > tol);
if (r == 0)
    pinvX = zeros(size(X'));
else
    s = diag(ones(r,1)./s(1:r));
    pinvX = V(:,1:r)*s*U(:,1:r)';
end
b = pinvX*Y;

% force positive fit
b(1) = abs(b(1));

%force negative b2 fit
b(2) = -(abs(b(2)));

% The center of the pRF should be positive. Thus, b(1)+b(2)>=0, or
% b(2) should be larger than -b(1) (implemented).
b(2) = max(b(2),-b(1));

% force second Gaussian to be negative
b(2) = -abs(b(2));

% compute residual sum of squares (e)
e = norm(Y - X*b);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

% N- point IDFTs of a 4-point sequence

X=[10 -2+2j -2 -2-2j];

% 4-point IDFT
ifft(X)

% 4-point IDFT
idft(X)

% 6-point IDFT
ifft(X,6)

% 6-point IDFT
X(6)=0
ifft(X)


text:
function Data = findRxnMatch(cRxn, nMatch, scoreTotal)
% Finds the highest scoring matches for a given reaction.
% Called by `reactionCompareGUI`.
%
% USAGE:
%
%    Data = findRxnMatch(cRxn, nMatch, scoreTotal)
%
% INPUTS:
%    cRxn:          Number of the compared rxn.
%    nMatch:        Number of matches to return.
%    scoreTotal:    summed SCORE of reactions vs reactions.
%    CMODEL:        global input
%    TMODEL:        global input
%
% OUTPUTS:
%    Data:          Structure containing comparison information.
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

global CMODEL TMODEL % Declare variables.

% Structure comtaining comparison information
Data = ([]);

%% cRxn Reaction Info
Data.cRxnTable = cell(8,1) ;
Data.cRxnTable(:) = {''} ;
Data.cRxnTable{1} = num2str(cRxn) ;
Data.cRxnTable{2} = CMODEL.rxns{cRxn} ;
Data.cRxnTable{3} = CMODEL.rxnNames{cRxn} ;
Data.cRxnTable{4} = CMODEL.rxnEquations{cRxn} ;
Data.cRxnTable{5} = CMODEL.rxnECNumbers{cRxn} ;
Data.cRxnTable{6} = CMODEL.rxnKEGGID{cRxn} ;
Data.cRxnTable{7} = CMODEL.rxnSEEDID{cRxn} ;
Data.cRxnTable{8} = CMODEL.subSystems{cRxn} ;

%% cRxn Metabolite Info
Data.cMetTable = cell(14, 1);
Data.cMetTable(:) = {''};

% Reactants from model
Data.cMetTable{1} = num2str(CMODEL.metNums(cRxn, 3)); % # of reactants

metPos = find(CMODEL.S(:, cRxn) < 0) ; % find reactants
if ~isempty(metPos)
    metlist = CMODEL.mets{metPos(1)} ; % add in first metabolite Id
    names = CMODEL.metNames{metPos(1)} ;
    stoics = num2str(abs(CMODEL.S(metPos(1), cRxn))) ; % add in first stoic
    forms = CMODEL.metFormulas{metPos(1)} ; % first formula
    charges = num2str(CMODEL.metCharge(metPos(1))) ; % first charge
    keggs = CMODEL.metKEGGID{metPos(1)} ; % first charge

    % If there are more than one mets, add additional info.
    if length(metPos) > 1
        for j = 2:length(metPos)
            metlist = strcat(metlist, '; ', CMODEL.mets{metPos(j)});
            names = strcat(names, '; ', CMODEL.metNames{metPos(j)});
            stoics = strcat(stoics, '; ', ...
                            num2str(abs(CMODEL.S(metPos(j), cRxn))));
            forms = strcat(forms, '; ', CMODEL.metFormulas{metPos(j)});
            charges = strcat(charges, '; ', ...
                             num2str(CMODEL.metCharge(metPos(j))));
            keggs = strcat(keggs, '; ', CMODEL.metKEGGID{metPos(j)});
        end
    end

    % Assign to table.
    Data.cMetTable{2} = metlist ;
    Data.cMetTable{3} = names ;
    Data.cMetTable{4} = stoics ;
    Data.cMetTable{5} = forms ;
    Data.cMetTable{6} = charges ;
    Data.cMetTable{7} = keggs ;
end

% Products from crxn.
% Number of products.
Data.cMetTable{8} = num2str(CMODEL.metNums(cRxn, 5));

% Find product position.
metPos = find(CMODEL.S(:, cRxn) > 0) ;

if ~isempty(metPos)
    % First product.
    metlist = CMODEL.mets{metPos(1)};
    names = CMODEL.metNames{metPos(1)} ;
    stoics = num2str(abs(CMODEL.S(metPos(1), cRxn)));
    forms = CMODEL.metFormulas{metPos(1)};
    charges = num2str(CMODEL.metCharge(metPos(1)));
    keggs = CMODEL.metKEGGID{metPos(1)};

    % Additional products.
    if length(metPos) > 1
        for j = 2:length(metPos)
            metlist = strcat(metlist, '; ', CMODEL.mets{metPos(j)});
            names = strcat(names, '; ', CMODEL.metNames{metPos(j)});
            stoics = strcat(stoics, '; ', ...
                            num2str(abs(CMODEL.S(metPos(j), cRxn))));
            forms = strcat(forms, '; ', CMODEL.metFormulas{metPos(j)});
            charges = strcat(charges, '; ', ...
                             num2str(CMODEL.metCharge(metPos(j))));
            keggs = strcat(keggs, '; ', CMODEL.metKEGGID{metPos(j)});
        end
    end

    % Assign data.
    Data.cMetTable{9} = metlist ;
    Data.cMetTable{10} = names ;
    Data.cMetTable{11} = stoics ;
    Data.cMetTable{12} = forms ;
    Data.cMetTable{13} = charges ;
    Data.cMetTable{14} = keggs ;
    clear metlist; clear stoics; clear forms; clear charges; clear keggs
end

%% Match Reaction(s) Info
% Allocate match reaction table.
Data.tRxnTable = cell(8, nMatch);
Data.tRxnTable(:) = {''};

% Sort scores, I is the original index.
[y, I] = sort(scoreTotal(cRxn, :), 'descend');

for i = 1:nMatch
    Data.tRxnTable{1, i} = strcat(num2str(y(i)), '; ', num2str(I(i)));
    Data.tRxnTable{2, i} = TMODEL.rxns{I(i)} ;
    Data.tRxnTable{3, i} = TMODEL.rxnNames{I(i)} ;
    Data.tRxnTable{4, i} = TMODEL.rxnEquations{I(i)} ;
    Data.tRxnTable{5, i} = TMODEL.rxnECNumbers{I(i)} ;
    Data.tRxnTable{6, i} = TMODEL.rxnKEGGID{I(i)} ;
    Data.tRxnTable{7, i} = TMODEL.rxnSEEDID{I(i)} ;
    Data.tRxnTable{8, i} = TMODEL.subSystems{I(i)} ;
end

%% Match Metabolite Info
Data.tMetTable = cell(14, nMatch);
Data.tMetTable(:) = {''};

for i = 1:nMatch
    % Reactants from model.
    Data.tMetTable{1,i} = num2str(TMODEL.metNums(I(i),3));

    % First reactant.
    metPos = find(TMODEL.S(:,I(i)) < 0) ;

    if ~isempty(metPos)
        metlist = TMODEL.mets{metPos(1)};
        names = TMODEL.metNames{metPos(1)} ;
        stoics = num2str(abs(TMODEL.S(metPos(1),I(i))));
        forms = TMODEL.metFormulas{metPos(1)};
        charges = num2str(TMODEL.metCharge(metPos(1)));
        keggs = TMODEL.metKEGGID{metPos(1)};

        % Additioanl Reactants.
        if length(metPos) > 1
            for j = 2:length(metPos)
                metlist = strcat(metlist, '; ', TMODEL.mets{metPos(j)});
                names = strcat(names, '; ', TMODEL.metNames{metPos(j)});
                stoics = strcat(stoics, '; ', ...
                                num2str(abs(TMODEL.S(metPos(j), I(i)))));
                forms = strcat(forms, '; ', TMODEL.metFormulas{metPos(j)});
                charges = strcat(charges, '; ', ...
                                 num2str(TMODEL.metCharge(metPos(j))));
                keggs = strcat(keggs, '; ', TMODEL.metKEGGID{metPos(j)});
            end
        end

        Data.tMetTable{2, i} = metlist ;
        Data.tMetTable{3, i} = names ;
        Data.tMetTable{4, i} = stoics ;
        Data.tMetTable{5, i} = forms ;
        Data.tMetTable{6, i} = charges ;
        Data.tMetTable{7, i} = keggs ;
    end

    % Products from model.
    if TMODEL.metNums(I(i), 5) ~= 0
        Data.tMetTable{8, i} = num2str(TMODEL.metNums(I(i), 5));

        % Find product indexes.
        metPos = find(TMODEL.S(:, I(i)) > 0) ;

        % First product
        metlist = TMODEL.mets{metPos(1)};
        names = TMODEL.metNames{metPos(1)} ;
        stoics = num2str(abs(TMODEL.S(metPos(1), I(i))));
        forms = TMODEL.metFormulas{metPos(1)};
        charges = num2str(TMODEL.metCharge(metPos(1)));
        keggs = TMODEL.metKEGGID{metPos(1)};

        % Additional products.
        if length(metPos) > 1
            for j = 2:length(metPos)
                metlist = strcat(metlist, '; ', TMODEL.mets{metPos(j)});
                names = strcat(names, '; ', TMODEL.metNames{metPos(j)});
                stoics = strcat(stoics, '; ', ...
                                num2str(abs(TMODEL.S(metPos(j), I(i)))));
                forms = strcat(forms, '; ', TMODEL.metFormulas{metPos(j)});
                charges = strcat(charges, '; ', ...
                                 num2str(TMODEL.metCharge(metPos(j))));
                keggs = strcat(keggs, '; ', TMODEL.metKEGGID{metPos(j)});
            end
        end

        % Assign to table.
        Data.tMetTable{9, i} = metlist ;
        Data.tMetTable{10, i} = names ;
        Data.tMetTable{11, i} = stoics ;
        Data.tMetTable{12, i} = forms ;
        Data.tMetTable{13, i} = charges ;
        Data.tMetTable{14, i} = keggs ;
    end
end

author: opencobra

repo: cobratoolbox


text:
function format_string = guess_nitf_format(imagesubhdr)
% GUESS_NITF_FORMAT Given a NITF image subheader structure, this functions
% tries to determine which set of sensor specific adjustments might need to
% be applied
%
% Assumes a set of test functions (all of the format isnitf*.m) are in the
% same directory.  If format is unrecognized, an empty string is returned.
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

fullname=which(mfilename('fullpath'));
path=fileparts(fullname);
filelist=dir([path filesep 'isnitf*.m']);

format_string=[];
for i=1:length(filelist)
    str_ind=strfind(filelist(i).name,'.m');
    if ~isempty(str_ind)&&(str_ind(end)>1)
        test_format=filelist(i).name(1:(str_ind(end)-1)); % Strip off .m extension
        if(feval(test_format, imagesubhdr))
            format_string=test_format(7:end);
            return
        end
    end
end

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

author: ngageoint


text:
function tagS = general_study_module_tags_subset
%"general_study_module_tags_subset"
%   Return a subset of the tags specified in C.7.2.1 in PS3.3 of 2006 DICOM
%   specification required to distinguish different studies.
%
%   Tags are returned in a struct array with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array.
%
%AI 11/20/2020
%
%Usage:
%   tagS = general_study_module_tags_subset
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

tagC = {'0020000D','00080020','00080030','00200010','00080050'};
tagDecC = {2097165, 524320, 524336, 2097168, 524368};
tagTypeC = {'1','2','2','2','2'};
childC = cell(1,length(tagC));

tagS = struct('tag',tagC,'tagdec',tagDecC,'type',tagTypeC,'children',childC);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%%% [4dmatrix] = convert_video_to_4d(struct)
%
function [V] = convert_video_to_4d(video)
  for i = 1:video.nrFramesTotal
    frame = video.frames(i).cdata;
    V(:,:,:,i) = frame;
  end
end

author: andrewssobral

repo: lrslibrary


text:
function [ar,arp,aru,g]=v_lpcrf2ar(rf)
%V_LPCRF2AR Convert reflection coefs to autoregressive coefs [AR,ARP,ARU,G]=(RF)
%
% Input:  RF(:,p+1) gives reflection coefficients of one or more p-section lossless tubes
% Ouputs: G is the gain of the all-pole AR filter
%         AR/G is the transfer function from U_in to the glottal input wave, U_g.
%               AR(:,1)=1 always.
%         ARP*K is the transfer function from U_in to the pressure just after the glottis
%               where K = rho*c/Alips: rho = air density 1.23 kg/m^3, c=sound speed 340 m/s,
%               Alips = effective area of free space beyond the lips.
%         ARU is the transfer function from U_in to the total volume velocity through the glottis
%
%              where U_in=z^(p/2)*U_lips is the time-advanced volume velocity at the lips
%
%         Energy into the vcal tract is equal to K*filter(ARP,1,Ulips).*filter(ARU,1,Ulips)
%              reverse glottal flows divided by 1-r0 where r0 is the glottal reflection coefficient.
%              The scale factor is included to avoid a zero answer when the glottis is closed giving r0=1.
%
% The transfer functions have ar(:,1)=art(:,1)=1
% They should both be multiplied by z^(p/2)/prod(1+rf) to correct the absolute
% gain and to compensate for the delay of p/2 samples along the length of the tube.
%
% The energy into the vocal tract is given by ars(speech) * are(speech)
%
% Ref: D. M. Brookes and H. P. Loke. "Modelling energy flow in the vocal tract with
%           applications to glottal closure and opening detection." In Proc ICASSP'99, pages 213-216, Mar 1999.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: v_lpcrf2ar.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1]=size(rf);
p2=p1+1;
p=p1-1;
pm=p-1;
arf=[ones(nf,1) zeros(nf,p)];
arr=[zeros(nf,p) rf(:,p1)];
cr=zeros(nf,p);
for k=1:p-1
  rk=rf(:,(p1-k)*ones(1,k));
  cr(:,1:k)=arr(:,p2-k:p1);
  arr(:,p1-k:p)=arr(:,p1-k:p)+rk.*arf(:,1:k);
  arf(:,2:k+1)=arf(:,2:k+1)+rk.*cr(:,1:k);
end
r1=rf(:,1);
ar=arf+r1(:,ones(1,p1)).*arr;
if nargout>1
   kp=prod(1-rf(:,2:p1),2);
   arp=(arf-arr)./kp(:,ones(1,p1));
   if nargout>2
      g=prod(1+rf(:,2:p1),2);
      aru=(arf+arr)./g(:,ones(1,p1));
      if nargout>3
         g=g.*(1+rf(:,1));
      end
   end
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [rank,description,usedcards,unusedcards] = rankp(hands,handsize)
% A function that evaluates poker hands and returns a numerical rank and
% text description for each. The list of cards used to make the best
% possible hand is returned in 'usedcards'. If handsize is smaller than
% the length of hands given, then the cards not used in the best possible
% hand are returned in 'unusedcards'. Rank of 1 is best (Royal Flush),
% higher numbers are worse. Rank values are not contiguous across all
% hands, but are correctly ordered.
%
% The hand matrix expected is an mxn list of m hands with n cards.  The
% cards are numbered from 1 to 52.  Suit doesn't matter for ranking. Order
% of the cards in the input vector does not matter. Numbering starts with
% the Ace at position 1. Here is a suggested card assignment:
%
% A | K | Q | J | 10| 9 | 8 | 7 | 6 | 5 | 4 | 3 | 2 | Suit
% -------------------------------------------------------------
% 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10| 11| 12| 13| Spades
% 14| 15| 16| 17| 18| 19| 20| 21| 22| 23| 24| 25| 26| Hearts
% 27| 28| 29| 30| 31| 32| 33| 34| 35| 36| 37| 38| 39| Clubs
% 40| 41| 42| 43| 44| 45| 46| 47| 48| 49| 50| 51| 52| Diamonds
%
% Example call to the function:
% trick = [33,3,43,1,25,4,19;
%          6,36,20,10,35,21,38;
%          45,31,51,40,37,22,26;
%          28,49,46,32,17,16,39];
% [ranks,descriptions,usedcards,unusedcards] = rankp(trick,2);
%
% revision 5
% - fixed a conditional error in the numerical ranking section,
%   thanks to Tim C. for testing it
% revision 4
% - fixed a hole in the condition to wipe low aces if not used
% revision 3
% - made rankp "toolbox neutral" by using nchoosek instead of combnk,
%   thanks to James for pointing that out
% - fixed a hole in the conditional statements for unusedcards
% revision 2
% - added handsize input, usedcards and unusedcards outputs
% - improved speed of classification with help from Doug
% revision 1
% - written for Doug's LazyWeb by Rob

%% Setup variables

handlen = size(hands,2);
numhands = size(hands,1);

% quality check on handsize argument
if nargin > 1
    handsize = min([handsize handlen 5]);
else
    handsize = min([handlen 5]);
end

rank = zeros(numhands,1);
description = cell(numhands,1);
usedcards = zeros(numhands,handsize);
if handlen > handsize
    unusedcards = zeros(numhands,(handlen-handsize));
else
    unusedcards = [];
end

%% Main process loop

% test each hand, keep max rank
for i = 1:numhands
    handcombs = nchoosek(hands(i,:),handsize);
    a = zeros(size(handcombs,1),1);
    b = cell(size(a));
    for j = 1:size(handcombs,1)
        [a(j),b{j}] = rankhand(handcombs(j,:));
    end
    [a,ix] = sort(a);
    rank(i) = a(1);
    description{i} = b(ix(1));
    usedcards(i,:) = handcombs(ix(1),:);
    if handlen > handsize
        unusedcards(i,:) = setxor(usedcards(i,:),hands(i,:));
    end
end

%% Sub function

function [r,handclass] = rankhand(hand)

% Setup variables
if nargin < 1 || length(hand) > 5
    r = 0;
    handclass = 'Hand must have 1 to 5 cards!';
    return;
end

len = length(hand);
cards = zeros(13,4);
cards(hand) = 1;
cards(14,:) = cards(1,:);

colsum = sum(cards);
rowsum = sum(cards,2);
rownames = {'Ace';'King';'Queen';'Jack';'10';'9';'8';'7';'6';'5';'4'; ...
    '3';'2';'Ace'};

handclass = 'High';
done = false;

%% Classify the hand

% test for Straight and Straight-flush, wipe low Ace if not used
if len==5
    if any(sum(logical([rowsum(1:5) rowsum(2:6) rowsum(3:7) rowsum(4:8)...
            rowsum(5:9) rowsum(6:10) rowsum(7:11) rowsum(8:12)...
            rowsum(9:13) rowsum(10:14)]))==5)
        if sum(logical(colsum))==1
            handclass = 'Straigh Flush';
        else
            handclass = 'Straight';
        end
        if sum(logical(rowsum(10:14))) == 5
            cards(1,:) = 0;
        else
            cards(14,:) = 0;
        end
        colsum = sum(cards);
        rowsum = sum(cards,2);
        done = true;
    end
end
if ~done
    cards(14,:) = 0;
    colsum = sum(cards);
    rowsum = sum(cards,2);
end

% test for Four of a Kind
if len >= 4
    if ~done && sum(rowsum==4)==1
        handclass = 'Four';
        done = true;
    end
end

% test for full house
if len == 5
    if ~done && sum(rowsum==3)==1 && sum(rowsum==2)==1
        handclass = 'Full House';
        done = true;
    end
% test for Flush
    if ~done && sum(colsum==5)==1
        handclass = 'Flush';
        done = true;
    end
end

% test for 3 of a kind
if len >= 3
    if ~done && sum(rowsum==3)==1
        handclass = 'Three';
        done = true;
    end
end

% test for 2 pairs
if len >= 4
    if ~done && sum(rowsum==2)==2
        handclass = 'Two Pairs';
        done = true;
    end
end

% test for 1 pair
if len >= 2
    if ~done && sum(rowsum==2)==1
        handclass = 'Pair';
    end
end

% if nothing else, handclass stays 'High'

%% Numerically rank the hand within its class

switch handclass
    case 'Straigh Flush'
        r = find(rowsum,1,'first');
    case 'Four'
        r = 1e1*find(rowsum==4);
        if len == 5
            r = r + find(rowsum==1);
        end
    case 'Full House'
        r = 2e2*find(rowsum==3) + find(rowsum==2);
    case 'Flush'
        a = find(rowsum);
        r = 3e3*a(1)+1e2*a(2)+1e2/2*a(3)+1e1*a(4)+a(5);
    case 'Straight'
        r = 4e4*find(rowsum,1,'first');
    case 'Three'
        r = 5e5*find(rowsum==3);
        if len > 3
            a = find(rowsum==1);
            r = r+1e4*a(1);
            if len==5
                r = r+1e3*a(2);
            end
        end
    case 'Two Pairs'
        a = find(rowsum==2);
        r = 7e6*a(1)+1e5*a(2);
        if len == 5
            b = find(rowsum==1);
            r = r+1e4*b(1);
        end
    case 'Pair'
        a = find(rowsum==2);
        r = 2e8*a(1);
        if len > 2
            b = find(rowsum==1);
            r = r+1e6*b(1);
            if len > 3
                r = r+1e5*b(2);
                if len == 5
                    r = r+1e4*b(3);
                end
            end
        end
    case 'High'
        a = find(rowsum);
        r = 5e9*a(1);
        if len > 1
            r = r+1e7*a(2);
            if len > 2
                r = r+1e6*a(3);
                if len > 3
                    r = r+1e5*a(4);
                    if len == 5
                        r = r+a(5);
                    end
                end
            end
        end
    otherwise
        disp('There was an error classifying this hand!');
        return;
end

%% Generate text output to describe the hand

switch handclass
    case 'Straigh Flush'
        if r == 1
            handclass = 'Royal Flush!';
        else
            handclass = [handclass,' to the ',rownames{find(rowsum,1,'first')}];
        end
    case 'Four'
        handclass = [handclass,' ',rownames{find(rowsum==4)},'''s'];
        if len == 5
            handclass = [handclass,' and a ',rownames{find(rowsum==1)}];
        end
    case 'Full House'
        handclass = [handclass,', ',rownames{find(rowsum==3)},'''s and ',...
            rownames{find(rowsum==2)},'''s'];
    case 'Flush'
        handclass = [handclass,' with ',rownames{a(1)},', ',rownames{a(2)},...
            ', ',rownames{a(3)},', ',rownames{a(4)},', and ',rownames{a(5)}];
    case 'Straight'
        handclass = [handclass,' to the ',rownames{find(rowsum,1,'first')}];
    case 'Three'
        handclass = [handclass,' ',rownames{find(rowsum==3)},'s'];
        if len > 3
            handclass = [handclass,' with ',rownames{a(1)}];
            if len == 5
                handclass = [handclass,', and ',rownames{a(2)}];
            end
        end
    case 'Two Pairs'
        handclass = [handclass,', ',rownames{a(1)},'''s and ',rownames{a(2)},'''s'];
        if len == 5
            handclass = [handclass,' with a(n) ',rownames{b(1)}];
        end
    case 'Pair'
        handclass = [handclass,' of ',rownames{a(1)},'''s'];
        if len > 2
            handclass = [handclass,' with ',rownames{b(1)}];
            if len > 3
                handclass = [handclass,', ',rownames{b(2)}];
                if len == 5
                    handclass = [handclass,', and ',rownames{b(3)}];
                end
            end
        end
    case 'High'
        handclass = [rownames{a(1)},' ',handclass,' with ',rownames{a(2)}];
        if len > 2
            handclass = [handclass,', ',rownames{a(3)}];
            if len > 3
                handclass = [handclass,', ',rownames{a(4)}];
                if len == 5
                    handclass = [handclass,', and ',rownames{a(5)}];
                end
            end
        end
    otherwise
        disp('There was an error generating text for this hand.');
end

author: Sable

repo: mcbench-benchmarks


text:
% Fig. 8.19   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

numGs=1;
denGs=[1 0 0];  % s^2

numDs=.81*[1 .2];
denDs=[1 2];

numC=conv(numGs,numDs);
denC=conv(denGs,denDs);

[numCL,denCL]=feedback(numC,denC,1,1);
damp(denCL)

tf=30;
t=0:.2:tf;
y=step(numCL,denCL,t);

axis([0 30 0 1.5])
plot(t,y),grid
hold on

% discrete designs

T=1;
[numGz,denGz]=c2dm(numGs,denGs,T,'zoh');

numDz1=[.389 -.319];
denDz1=[1 -.135];

numDz2=.374*[1 -.85];
denDz2=[1 0];

numz1=conv(numGz,numDz1);
denz1=conv(denGz,denDz1);

numz2=conv(numGz,numDz2);
denz2=conv(denGz,denDz2);

[numCLz1,denCLz1]=feedback(numz1,denz1,1,1);
[numCLz2,denCLz2]=feedback(numz2,denz2,1,1);

ddamp(denCLz1,T)
ddamp(denCLz2,T)

N=tf/T+1;
td=0:1:tf;
yd1=dstep(numCLz1,denCLz1,N);
yd2=dstep(numCLz2,denCLz2,N);
plot(td,yd1,'-',td,yd1,'*')
plot(td,yd2,'-',td,yd2,'o')
text(10,.55,'----------- continuous design')
text(10,.75,'*----*----* discrete equivalent design')
text(10,.65,'o----o----o discrete design')
title('Fig. 8.19  Step responses of the various design methods')
xlabel('Time (sec)')
ylabel('Plant output')
hold off


text:
function hded1(currentfile)

[path,nam,ext] = fileparts(currentfile);
sig = signal(currentfile);
title = ['Axes Editor for ' nam];
lineNo = 1;
prompt = {};
defaults = {};
for i=1:ndim(sig)
	a = getaxis(sig, i);
	prompt{end+1} = ['Name of dimension ' num2str(i)];
	prompt{end+1} = ['Unit of dimension ' num2str(i)];
	prompt{end+1} = ['Start value and sampling rate for dimension ' num2str(i)];
	defaults{end+1} = name(a);
	defaults{end+1} = label(unit(a));
	defaults{end+1} = num2str([first(a) samplerate(a)]);
end

answer  = inputdialog(prompt,title,lineNo, defaults);
if ~isempty(answer)
	for i=1:ndim(sig)
		a = getaxis(sig, i);
		a = setunit(a, unit(answer{3*i-1}));
		if ~isempty(answer(3*i-2))
			a = setname(a, answer{3*i-2});
		end
		rs = str2num(answer{3*i});
		first = rs(1);
		rate =  rs(2);
		if rate > 0
			delta = 1/rate;
		else
			delta = 1;
		end
		a = setfirst(a, first);
		a = setdelta(a, delta);
		sig = setaxis(sig, i, a);
	end
	write(sig, currentfile);
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%% Deep Neural Network with Caffe models
% Load Caffe framework models.
%
% In this tutorial you will learn how to use DNN module for image
% classification by using GoogLeNet trained network from
% <http://caffe.berkeleyvision.org/model_zoo.html Caffe model zoo>.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.4.0/samples/dnn/caffe_googlenet.cpp>
% * <https://github.com/opencv/opencv/blob/3.4.0/samples/dnn/googlenet_python.py>
% * <https://docs.opencv.org/3.4.0/d5/de7/tutorial_dnn_googlenet.html>
%

%% BVLC GoogLeNet
% First, we download GoogLeNet model files:
%
% * |deploy.prototxt| and |bvlc_googlenet.caffemodel|
% * Also we need file with names of
%   <http://image-net.org/challenges/LSVRC/2012/browse-synsets ILSVRC2012>
%   classes: |synset_words.txt|.
%
dirDNN = fullfile(mexopencv.root(), 'test', 'dnn', 'GoogLeNet');
modelLabels = fullfile(dirDNN, 'synset_words.txt');
modelTxt = fullfile(dirDNN, 'deploy.prototxt');
modelBin = fullfile(dirDNN, 'bvlc_googlenet.caffemodel');  % 51 MB file
files = {modelLabels, modelTxt, modelBin};
urls = {
    'https://cdn.rawgit.com/opencv/opencv/3.4.0/samples/data/dnn/synset_words.txt';
    'https://cdn.rawgit.com/opencv/opencv/3.4.0/samples/data/dnn/bvlc_googlenet.prototxt';
    'http://dl.caffe.berkeleyvision.org/bvlc_googlenet.caffemodel';
};
if ~isdir(dirDNN), mkdir(dirDNN); end
for i=1:numel(files)
    if exist(files{i}, 'file') ~= 2
        disp('Downloading...')
        urlwrite(urls{i}, files{i});
    end
end

%% Load class labels
if ~mexopencv.isOctave()
    fid = fopen(modelLabels, 'rt');
    C = textscan(fid, '%*s %[^\n]');
    fclose(fid);
    labels = C{1};
else
    %HACK: textscan is buggy and unreliable in Octave!
    labels = textread(modelLabels, '%s', 'Delimiter','\n');
    labels = regexprep(labels, '^\w+\s*', '', 'once');
end
fprintf('%d classes\n', numel(labels));

%% Create and initialize network from Caffe model
net = cv.Net('Caffe', modelTxt, modelBin);
assert(~net.empty(), 'Cant load network');
if false
    net.setPreferableTarget('OpenCL');
end

%% Prepare blob from input image
% Read input image (BGR channel order)
img = cv.imread(fullfile(mexopencv.root(), 'test', 'space_shuttle.jpg'), ...
    'Color',true, 'FlipChannels',false);

%%
% we resize and convert the image to 4-dimensional blob (so-called batch)
% with 1x3x224x224 shape, because GoogLeNet accepts only 224x224 BGR-images.
% we also subtract the mean pixel value of the training dataset ILSVRC_2012
% (B: 104.0069879317889, G: 116.66876761696767, R: 122.6789143406786)
if true
    blob = cv.Net.blobFromImages(img, ...
        'Size',[224 224], 'Mean',[104 117 123], 'SwapRB',false);
else
    % NOTE: blobFromImages does crop/resize to preserve aspect ratio of image
    blob = cv.resize(img, [224 224]);                % Size
    blob = single(blob);                             % CV_32F
    blob = bsxfun(@minus, blob, cat(3,104,117,123)); % Mean (BGR)
    blob = permute(blob, [4 3 1 2]);                 % HWCN -> NCHW
end

%% Set the network input
% In |deploy.prototxt| the network input blob named as "data".
% Other blobs labeled as "name_of_layer.name_of_layer_output".
net.setInput(blob, 'data');

%% Make forward pass and compute output
% During the forward pass output of each network layer is computed,
% but in this example we need output from "prob" layer only.
tic
p = net.forward('prob');  % 1x1000 vector
toc

%% Gather output of "prob" layer
% We determine the best class by taking the output of "prob" layer, which
% contain probabilities for each of 1000 ILSVRC2012 image classes, and finding
% the index of element with maximal value in this one. This index correspond
% to the class of the image.
[~,ord] = sort(p, 'descend');  % ordered by maximum probability

%% Show predictions
if ~mexopencv.isOctave()
    %HACK: TABLE not implemented in Octave
    t = table(labels(:), p(:), 'VariableNames',{'Label','Probability'});
    t = sortrows(t, 'Probability', 'descend');
    disp('Top 5 predictions');
    disp(t(1:5,:));
end

subplot(3,1,1:2), imshow(flip(img,3))
title('space_shuttle.jpg', 'Interpreter','none')
subplot(3,1,3), barh(1:5, p(ord(1:5))*100)
set(gca, 'YTickLabel',labels(ord(1:5)), 'YDir','reverse')
axis([0 100 0 6]), grid on
title('Top-5 Predictions'), xlabel('Probability (%)')

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% NURBS Toolbox.
% Version 1.0
%
% demos        - NURBS demonstrations
%
% nrbmak       - Construct a NURBS from control points and knots.
% nrbtform     - Applying scaling, translation or rotation operators.
% nrbkntins    - Knot insertion/refinement.
% nrbdegelev   - Degree elevation.
% nrbderiv     - NURBS representation of the derivative.
% nrbdeval     - Evaluation of the NURBS derivative.
% nrbkntmult   - Find the multiplilicity of a knot vector.
% nrbreverse   - Reverse evaluation direction of the NURBS.
% nrbtransp    - Swap U and V for NURBS surface.
% nrbline      - Construct a straight line.
% nrbcirc      - Construct a circular arc.
% nrbrect      - Construct a rectangle.
% nrb4surf     - Surface defined by 4 corner points.
% nrbeval      - Evalution of NURBS curve or surface.
% nrbextrude   - Extrude a NURBS curve along a vector.
% nrbrevolve   - Construct surface by revolving a profile.
% nrbruled     - Ruled surface between twp NURBS curves.
% nrbcoons     - Construct Coons bilinearly blended surface patch.
% nrbplot      - Plot NURBS curve or surface.
%
% bspeval      - Evaluate a univariate B-Spline.
% bspderiv     - B-Spline representation of the derivative
% bspkntins    - Insert a knot or knots into a univariate B-Spline.
% bspdegelev   - Degree elevation of a univariate B-Spline.
%
% vecnorm      - Normalise the vectors.
% vecmag       - Magnitaude of the vectors.
% vecmag2      - Squared Magnitude of the vectors.
% vecangle     - Alternative to atan2 (0 <= angle <= 2*pi)
% vecdot       - Dot product of two vectors.
% veccross     - Cross product of two vectors.
% vecrotx      - Rotation matrix around the x-axis.
% vecroty      - Rotation matrix around the y-axis.
% vecrotz      - Rotation matrix around the z-axis.
% vecscale     - Scaling matrix.
% vectrans     - Translation matrix.
%
% deg2rad      - Convert degrees to radians.
% rad2deg      - Convert radians to degrees.


text:
function [x,y] = combinatorialDependentRows(A)

[m,n]=size(A);
% A is a binarized [F R] matrix (each entry of A is 0 or 1)
% if this program is feasible, then there exist 2 combinatorially dependent subsets of rows of A

%cvx_solver gurobi
cvx_solver mosek

cvx_begin sdp quiet

    variable x(m) binary;
    variable y(m) binary;

    % x,y are nonempty and disjoint subsets of rows of A
    x+y<=1;
    sum(x)>=1;
    sum(y)>=1;

    % trick to test for same support (idea: m > largest possible sum)
    m*x'*A>=y'*A;
    m*y'*A>=x'*A;

cvx_end

disp(find(x))
disp(find(y))

end


text:
function Q = costMat(tau, r)
%costMat constructs a cost matrix, Q for a single segment
% r: The order of the derivative subjected to optimization
% ex) r = 4: minimum snap, r = 2: minimum acceleration

% nth order poly
n = 2*r+1;

Q = zeros(n+1);
QQ = zeros(n+1,n+1,length(r));
% eg optimizing 4th derivative
% r = 4;
% row
for rr = 1:length(r)
    for i=0:size(Q,1)-1
        % column
        for j=0:size(Q,2)-1
            if i >= r(rr) && j >= r(rr)
                m = 0:r(rr)-1;
                QQ(i+1,j+1,rr) = 2*prod((i-m).*(j-m))*tau^(i+j-2*r(rr)+1)/(i+j-2*r(rr)+1);
            end
        end
    end
end

for i = 1:length(r)
    Q = Q + QQ(:,:,i);
end

author: yorgoon

repo: minimum-snap-geometric-control


text:
classdef LvdGeometry < matlab.mixin.SetGet
    %LvdGeometry Summary of this class goes here
    %   Detailed explanation goes here

    properties
        points GeometricPointSet = GeometricPointSet.empty(1,0)
        vectors GeometricVectorSet = GeometricVectorSet.empty(1,0)
        coordSyses GeometricCoordSysSet = GeometricCoordSysSet.empty(1,0)
        refFrames GeometricRefFrameSet = GeometricRefFrameSet.empty(1,0)
        angles GeometricAngleSet = GeometricAngleSet.empty(1,0);
        planes GeometricPlaneSet = GeometricPlaneSet.empty(1,0);

        lvdData LvdData
    end

    methods
        function obj = LvdGeometry(lvdData)
            obj.lvdData = lvdData;

            obj.points = GeometricPointSet(lvdData);
            obj.vectors = GeometricVectorSet(lvdData);
            obj.coordSyses = GeometricCoordSysSet(lvdData);
            obj.refFrames = GeometricRefFrameSet(lvdData);
            obj.angles = GeometricAngleSet(lvdData);
            obj.planes = GeometricPlaneSet(lvdData);
        end

        function tf = usesGroundObj(obj, groundObj)
            tf = false;
            tf = tf || obj.points.usesGroundObj(groundObj);
            tf = tf || obj.vectors.usesGroundObj(groundObj);
            tf = tf || obj.coordSyses.usesGroundObj(groundObj);
            tf = tf || obj.refFrames.usesGroundObj(groundObj);
            tf = tf || obj.angles.usesGroundObj(groundObj);
        end

        function tf = usesGeometricPoint(obj, point)
            tf = false;
            tf = tf || obj.points.usesGeometricPoint(point);
            tf = tf || obj.vectors.usesGeometricPoint(point);
            tf = tf || obj.coordSyses.usesGeometricPoint(point);
            tf = tf || obj.refFrames.usesGeometricPoint(point);
            tf = tf || obj.angles.usesGeometricPoint(point);
        end

        function tf = usesGeometricVector(obj, vector)
            tf = false;
            tf = tf || obj.points.usesGeometricVector(vector);
            tf = tf || obj.vectors.usesGeometricVector(vector);
            tf = tf || obj.coordSyses.usesGeometricVector(vector);
            tf = tf || obj.refFrames.usesGeometricVector(vector);
            tf = tf || obj.angles.usesGeometricVector(vector);
        end

        function tf = usesGeometricCoordSys(obj, coordSys)
            tf = false;
            tf = tf || obj.points.usesGeometricCoordSys(coordSys);
            tf = tf || obj.vectors.usesGeometricCoordSys(coordSys);
            tf = tf || obj.coordSyses.usesGeometricCoordSys(coordSys);
            tf = tf || obj.refFrames.usesGeometricCoordSys(coordSys);
            tf = tf || obj.angles.usesGeometricCoordSys(coordSys);
        end

        function tf = usesGeometricRefFrame(obj, refFrame)
            tf = false;
            tf = tf || obj.points.usesGeometricRefFrame(refFrame);
            tf = tf || obj.vectors.usesGeometricRefFrame(refFrame);
            tf = tf || obj.coordSyses.usesGeometricRefFrame(refFrame);
            tf = tf || obj.refFrames.usesGeometricRefFrame(refFrame);
            tf = tf || obj.angles.usesGeometricRefFrame(refFrame);
        end

        function tf = usesGeometricAngle(obj, angle)
            tf = false;
            tf = tf || obj.points.usesGeometricAngle(angle);
            tf = tf || obj.vectors.usesGeometricAngle(angle);
            tf = tf || obj.coordSyses.usesGeometricAngle(angle);
            tf = tf || obj.refFrames.usesGeometricAngle(angle);
            tf = tf || obj.angles.usesGeometricAngle(angle);
        end

        function tf = usesGeometricPlane(obj, plane)
            tf = false;
            tf = tf || obj.points.usesGeometricPlane(plane);
            tf = tf || obj.vectors.usesGeometricPlane(plane);
            tf = tf || obj.coordSyses.usesGeometricPlane(plane);
            tf = tf || obj.refFrames.usesGeometricPlane(plane);
            tf = tf || obj.angles.usesGeometricPlane(plane);
            tf = tf || obj.planes.usesGeometricPlane(plane);
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            if(isempty(obj.points))
                obj.points = GeometricPointSet(obj.lvdData);
            end

            if(isempty(obj.vectors))
                obj.vectors = GeometricVectorSet(obj.lvdData);
            end

            if(isempty(obj.coordSyses))
                obj.coordSyses = GeometricCoordSysSet(obj.lvdData);
            end

            if(isempty(obj.refFrames))
                obj.refFrames = GeometricRefFrameSet(obj.lvdData);
            end

            if(isempty(obj.angles))
                obj.angles = GeometricAngleSet(obj.lvdData);
            end

            if(isempty(obj.planes))
                obj.planes = GeometricPlaneSet(obj.lvdData);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ csd ] = bz_CSD (lfp, varargin)

% [ CSD ] = bz_CSD (lfp, varargin)
% Calculates the 1D approximation of current source density (CSD) from a
% linear array of LFPs

% INPUT
%    lfp            a buzcode structure with fields lfp.data,
%                                                   lfp.timestamps
%                                                   lfp.samplingRate
%                   -lfp can also be a [t x 1] timeseries signal. in which
%                   case you need to input 'samplingRate'
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%       channels    vector with channels to inlcude. If empty take all (default)
%       win         time interval to compute CSD. If empty take all (default)
%       spat_sm     degree of spatial smoothing. Default = 0.
%       temp_sm     degree of temporal smoothing. Default = 0.
%       plotCSD     true/false. Default true.
%       plotLFP     true/false. Default true.
%    =========================================================================

% OUTPUT:
%    CSD           a buzcode structure with fields csd.data,
%                                                   csd.timestamps
%                                                   csd.samplingRate
%                                                   csd.channels
%                                                   csd.params

% Antonio FR, 7/18

%% Parse inputs

p = inputParser;
addParameter(p,'channels',1:size(lfp.data,2),@isvector);
addParameter(p,'samplingRate',1250,@isnumeric);
addParameter(p,'win',[1 size(lfp.data,1)],@isnumeric);
addParameter(p,'spat_sm',11,@isnumeric);
addParameter(p,'temp_sm',11,@isnumeric);
addParameter(p,'doDetrend',false,@islogical);
addParameter(p,'plotCSD',true,@islogical);
addParameter(p,'plotLFP',true,@islogical);

parse(p,varargin{:});
channels = p.Results.channels;
samplingRate = p.Results.samplingRate;
spat_sm = p.Results.spat_sm;
temp_sm = p.Results.temp_sm;
doDetrend = p.Results.doDetrend;
plotCSD = p.Results.plotCSD;
plotLFP = p.Results.plotLFP;

%lfp input
if isstruct(lfp)
    data = lfp.data;
    timestamps = lfp.timestamps;
    samplingRate = lfp.samplingRate;
elseif iscell(lfp) %for multiple trials
    celllengths = cellfun(@length,lfp);
    data = vertcat(lfp{:});
elseif isnumeric(lfp)
    data = lfp;
    timestamps = [1:length(lfp)]'./samplingRate;
end

win = p.Results.win*samplingRate;

%% Compute CSD

lfp_frag = data(win(1):win(2),channels)*-1;

% detrend
if doDetrend
   lfp_frag = detrend(lfp_frag')';
end

% temporal smoothing
if temp_sm > 0
   for ch = 1:size(lfp_frag,2)
       lfp_frag(:,ch) = smooth(lfp_frag(:,ch),temp_sm,'sgolay');
   end
end

% spatial smoothing
if spat_sm > 0
   for t = 1:size(lfp_frag,1)
       lfp_frag(t,:) = smooth(lfp_frag(t,:),spat_sm,'lowess');
   end
end

% calculate CSD
CSD = diff(lfp_frag,2,2);

% generate output structure
csd.data = CSD;
csd.timestamps = timestamps(win(1):win(2));
csd.samplingRate = samplingRate;
csd.channels = channels;
csd.params.spat_sm = spat_sm;
csd.params.temp_sm = temp_sm;
csd.params.detrend = doDetrend;

%% Plot

if plotLFP

    cmax = max(max(CSD));

    figure;
    subplot(1,2,1);
    contourf(timestamps(win(1):win(2)),1:size(CSD,2),CSD',40,'LineColor','none');hold on;
    colormap jet; caxis([-cmax cmax]);
    set(gca,'YDir','reverse');xlabel('time (s)');ylabel('channel');title('CSD');

    subplot(1,2,2);
    for ch=1:size(lfp_frag,2)
        offset = 500*(ch-1);
        sh_tmp = 10e5*(lfp_frag(:,ch)) + offset;
        plot(timestamps(win(1):win(2)),sh_tmp,'k','LineWidth',1.5); hold on;
        clear sh_tmp
    end
    set(gca,'YDir','reverse','YTickLabel',[]);ylim([-500 offset+500]);xlim([timestamps(win(1)) timestamps(win(2))]);
    xlabel('time (s)');ylabel('channel');title('LFP');

elseif plotCSD

     cmax = max(max(CSD));

     figure;
     contourf(timestamps(win(1):win(2)),1:size(CSD,2),CSD',40,'LineColor','none');hold on;
     colormap jet; caxis([-cmax cmax]);
     set(gca,'YDir','reverse');xlabel('time (s)');ylabel('channel');title(CSD);

end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function net = train(net, tutor, varargin)

% TRAIN
%
% Train a max-win multi-class support vector classifier network using the
% specified tutor to train each component two-class network.
%
%    load data/iris x y;
%
%    C      = 100;
%    kernel = rbf(0.5);
%
%    net = train(maxwin, smosvctutor, x, y, C, kernel);

%
% File        : @maxwin/train.m
%
% Date        : Wednesday 13th September 2000
%
% Author      : Dr Gavin C. Cawley
%
% Description : Gateway function used to train a max-win multi-class support
%               vector classifier network using a given tutor.  Part of an
%               object-oriented implementation of Vapnik's Support Vector
%               Machine, as described in [1].
%
% References  : [1] V.N. Vapnik,
%                   "The Nature of Statistical Learning Theory",
%                   Springer-Verlag, New York, ISBN 0-387-94559-8,
%                   1995.
%
% History     : 13/09/2000 - v1.00
%
% Copyright   : (c) Dr Gavin C. Cawley, September 2000
%
%    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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

y = varargin{2};

if size(y, 2) == 1

   for i=1:max(y)

      varargin{2} = 2*(y == i) - 1;

      net.net(i) = train(tutor, varargin{:});

   end

else

   for i=1:size(y, 2)

      varargin{2} = y(:,i);

      net.net(i) = train(tutor, varargin{:});

   end

end

% bye bye...


text:
function prior = gaussian2PriorExpandParam(prior, params)

% GAUSSIAN2PRIOREXPANDPARAM Expand Gaussian prior structure from param vector.
% FORMAT
% DESC places the given parameters in a Gaussian prior structure.
% ARG prior : the structure to place the parameters in.
% ARG params : the parameters to place in the structure.
% RETURN prior : the structure with the parameters in place.
%
% SEEALSO priorExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2004
% MODIFICATIONS: Andreas C. Damianou, 2013

% SHEFFIELDML

prior.precision = params(1);
prior.mean = params(2:end);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function xy = perform_classical_mds(D,ndims)

N = size(D,1);
opt.disp = 0; opt.isreal = 1; opt.issym = 1;
M = -.5*(D.^2 - sum(D.^2)'*ones(1,N)/N - ones(N,1)*sum(D.^2)/N + sum(sum(D.^2))/(N^2));
[xy, val] = eigs(M, ndims, 'LR', opt);

for i=1:ndims
    xy(:,i) = xy(:,i)*sqrt(val(i,i));
end

author: gpeyre

repo: matlab-toolboxes


text:
classdef SigmoidNode < GraphNode

    methods
        function obj = SigmoidNode(dimOut)
            obj = obj@GraphNode('Sigmoid',dimOut);
        end

        function obj = forward(obj,prev_layers)
            obj = obj.preprocessingForward(prev_layers);

            input = prev_layers{1}.a;
            obj.a = sigmoid(input);
            obj = forward@GraphNode(obj, prev_layers);
        end

        function obj = backward(obj,prev_layers, future_layers)
            if obj.skipGrad || obj.skipBP
                return;
            end

            future_grad = obj.GetFutureGrad(future_layers);

            if obj.L1weight>0
                tmp = -obj.L1target./max(1e-3,obj.rho) + (1-obj.L1target)./max(1e-3,(1-obj.rho));
                future_grad = future_grad + repmat(obj.L1weight * tmp, 1, size(future_grad,2));
            end

            obj.grad{1} = future_grad .* obj.a .* (1-obj.a);

            obj = backward@GraphNode(obj, prev_layers, future_layers);
        end

    end

end

author: singaxiong

repo: SignalGraph


text:
function [S,t,f,Serr]=mtspecgramc(data,movingwin,params)
% Multi-taper time-frequency spectrum - continuous process
%
% Usage:
% [S,t,f,Serr]=mtspecgramc(data,movingwin,params)
% Input:
% Note units have to be consistent. Thus, if movingwin is in seconds, Fs
% has to be in Hz. see chronux.m for more information.
%       data        (in form samples x channels/trials) -- required
%       movingwin         (in the form [window winstep] i.e length of moving
%                                                 window and step size)
%                                                 Note that units here have
%                                                 to be consistent with
%                                                 units of Fs - required
%       params: structure with fields tapers, pad, Fs, fpass, err, trialave
%       - optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%                   Note that T has to be equal to movingwin(1).
%
%	        pad		    (padding factor for the FFT) - optional. Defaults to 0.
%			      	 e.g. For N = 500, if PAD = 0, we pad the FFT
%			      	 to 512 points; if PAD = 2, we pad the FFT
%			      	 to 2048 points, etc.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%           err  (error calculation [1 p] - Theoretical error bars; [2 p] - Jackknife error bars
%                                   [0 p] or 0 - no error bars) - optional. Default 0.
%           trialave (average over trials/channels when 1, don't average when 0) - optional. Default 0
% Output:
%       S       (spectrum in form time x frequency x channels/trials if trialave=0; in the form time x frequency if trialave=1)
%       t       (times)
%       f       (frequencies)
%       Serr    (error bars) only for err(1)>=1

if nargin < 2; error('Need data and window parameters'); end;
if nargin < 3; params=[]; end;

if length(params.tapers)==3 & movingwin(1)~=params.tapers(2);
    error('Duration of data in params.tapers is inconsistent with movingwin(1), modify params.tapers(2) to proceed')
end

[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
if nargout > 3 && err(1)==0;
%   Cannot compute error bars with err(1)=0. change params and run again.
    error('When Serr is desired, err(1) has to be non-zero.');
end;

N=size(data,1);
Nwin=round(Fs*movingwin(1)); % number of samples in window
Nstep=round(movingwin(2)*Fs); % number of samples to step through
nfft=2^(nextpow2(Nwin)+pad);

[f,findx]=getfgrid(Fs,nfft,fpass);
tapers=dpsschk(tapers,Nwin,Fs); % check tapers
[NC C]=size(data); % size of data
[NK K]=size(tapers); % size of tapers

tapers=tapers(:,:,ones(1,C)); % add channel indices to tapers

if NK~=Nwin ; error('length of tapers is incompatible with length of data'); end;

winstart=1:Nstep:N-Nwin+1;
nw=length(winstart);

S = zeros(nw,length(f),C);

for n=1:nw;
    J=fft((permute(data(winstart(n):winstart(n)+Nwin-1,:,ones(1,K)),[1 3 2]) .* tapers)  ,nfft)/Fs;   % fft of projected data
    J=J(findx,:,:);
    S(n,:,:)=squeeze(mean(conj(J).*J,2));
end

if nargout==4;Serr=squeeze(Serr);end;
winmid=winstart+round(Nwin/2);
t=winmid/Fs;


text:
% Find x that minimizes ||Ax|| under different constraints
%
% The problem can be solved with different constraints:
%  Reference: Alg 5.6, HZ2, p. 595
%   Find x that minimizes ||Ax|| subject to ||x||=1 and x=G*xHat, where G
%   has rank r
%  Reference: Alg 5.7, HZ2, p. 596
%   Find x that minimizes ||Ax|| subject to ||Cx||=1
%
% USAGE
%  x = solveLeastSqAx(A)
%  [x,xHat] = solveLeastSqAx(A,G,method)
%
% INPUTS 1 - Find x that minimizes ||Ax|| subject to ||x||=1
%  A       - constraint matrix, ||Ax|| to be minimized with ||x||=1
%
% INPUTS 2
%
% INPUTS 3 - Find x that minimizes ||Ax|| subject to ||x||=1 and x=G*xHat,
%            where G has rank r
%  A       - constraint matrix
%  G       - condition matrix
%  method  - method=2
%
% INPUTS 4 - Find x that minimizes ||Ax|| subject to ||Cx||=1
%  A       - constraint matrix
%  C       - condition matrix
%  method  - method=3
%
% OUTPUTS 1,2,4
%  x      - solution
%
% OUTPUTS 3
%  x      - solution
%  xHat   - vector such that x = G*xHat
%
% EXAMPLE
%
% See also
%
% Vincent's Structure From Motion Toolbox      Version 1.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

function [ x xHat ] = solveLeastSqAx(A,G,method)

if nargin==1; [U,D,V] = svd(A,0); x=V(:,end); return; end
if nargin==2; error('method argument required'); end

switch method
  case 1
    % Find x that minimizes ||Ax|| subject to ||x||=1
    % Reference: HZ2, Algorithm 5.4, p593
    if issparse(A)
      [ disc disc V ]=svds(A);
    else
      [ disc disc V ]=svd(A);
    end
    x = V(:,end);
  case 2
    % Find x that minimizes ||Ax|| subject to ||x||=1 and x=G*xHat, where G
    % has rank r
    % Reference: HZ2, Algorithm 5.6, p595
    % (i)
    [U,D,V] = svd(G,0);
    % (ii)
    r = rank(G); Up = U(:,1:r);
    % (iii)
    xp = solveLeastSqAx(A*U2);
    % (iv)
    x = Up*xp;
    % (v)
    if nargout==2; Vp=V(:,1:r); xHat=Vp*diag(1./diag(D(1:r,1:r)))*x2; end
  case 3
    %   Find x that minimizes ||Ax|| subject to ||Cx||=1
    %  Reference: Alg 5.7, HZ2, p. 596
    % (i)
    [U,D,V]=svd(C); Ap=A*V;
    % (ii)
    r=rank(D); A1p=Ap(:,1:r); A2p=Ap(:,r:end);
    % (iii)
    D1=D(1:r,1:r);
    % (iv)
    D1inv=diag(1./diag(D1)); A2pPinv=pinv(A2p);
    App=(A2p*A2pPinv-eye(size(A2p,1)))*A1p*D1inv;
    % (v)
    xpp=solveLeastSqAx(App);
    % (vi)
    x1p=D1inv*xpp; x2p=-A2pPinv*A1p*x1p; xp=[x1p;x2p];
    % (vii)
    x=V*xp;
end

author: vrabaud


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR:  How to use FAIRplots
%
% - loads (or generates data) based on E9_Hands_MLIR_SSD_mbElas
% - calls FAIRplotswith all options
%
% FAIRplots('clear');                  - clear the function:
% FAIRplots('reset','mode','PIR');     - reset and set mode to parametric image registration
% FAIRplots('init',para);              - initialize (T,R,omega,m,)
% FAIRplots('stop',para);              - display the stopping configuration
% FAIRplots('start',para);             - display the start configuration
% FAIRplots(k,para);                   - display the k'th configuration
%==============================================================================

clear; help(mfilename);

tempFile = fullfile(FAIRpath,'temp',[mfilename,'.mat'])
if ~exist(tempFile,'file'),

  % generate some dummy data to be visualized using FAIRplots
  % load data, initialize image viewer, interpolator, transformation
  E9_Hands_MLIR_SSD_mbElas
  viewImage('reset','viewImage','viewImage2D','colormap','gray(256)');
  imgModel('reset','imgModel','linearInter');
  trafo('reset','trafo','affine2D');

  % shortcurs for transformed image
  X = getCellCenteredGrid(omega,m);
  Rc = imgModel(dataR,omega,X);
  Tc = @(Y)  imgModel(dataT,omega,Y);
  Dc = @(Tc) SSD(Tc,Rc,omega,m);

  % compute the stopping, starting and final values
  builtin('clear','center')
  wStop = trafo('w0');
  Ystop = trafo(wStop,X);
  Tstop = Tc(Ystop);
  Y0    = trafo(wc,X);
  yOpt  = center(yc,m);
  T0    = Tc(Y0);
  Topt  = Tc(yOpt);
  save(tempFile,'dataT','dataR','Tc','Rc','omega','m','Tstop',...
      'Ystop','T0','Y0','Topt','yOpt')
  close all
end;
load(tempFile)

Dc = @(Tc) SSD(Tc,Rc,omega,m);
viewImage('reset','viewImage','viewImage2D','colormap',bone(256),'axis','off');
imgModel('reset','imgModel','linearInter');
trafo('reset','trafo','affine2D');

% reset FAIRplots to have a fresh start
FAIRplots('clear');

% setup the plotting functionality, no plots yet
FAIRplots('reset','mode','PIR')%mfilename);

% initialize the plots, show T and R
FAIRplots('init',struct('Tc',dataT,'Rc',dataR,'omega',omega,'m',m));
FAIRpause(2);

% show the stopping values
para  = struct('Tc',Tstop,'Rc',Rc,'omega',omega,'m',m,'yc',Ystop,'Jc',Dc(Tstop));
FAIRplots('stop',para);
FAIRpause(2);

% show the starting values
para.Tc = T0; para.yc = Y0; para.Jstop = para.Jc; para.Jc = Dc(T0);
FAIRplots('start',para);
FAIRpause(2);

% show the final values
para.Tc = Topt; para.yc = yOpt; para.Jc = Dc(Topt);
para.normdY = norm(yOpt-Y0)/norm(Ystop);
FAIRplots(1,para);
%==============================================================================