text:
function MV = ffmInteractionsSparse(X,Xbox,Y,Ybox,V,Ibox,green,k,edg,tol)
%+========================================================================+
%|                                                                        |
%|         OPENFFM - LIBRARY FOR FAST AND FREE MEMORY CONVOLUTION         |
%|           openFfm is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2019.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : ffmInteractionsSparse.m                       |
%|    #    |   VERSION    : 0.61                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : Sparse product for low-frequency compressible |
%|  `---'  |                leaves                                        |
%+========================================================================+

% Initialisation du produit Matrice-Vecteur
MV = zeros(size(X,1),1,class(V));

% Quadrature des interpolations lagrangiennes
[Xq,ii,jj,kk,xq] = ffmQuadrature(green,k,edg,tol);

% Unicite des vecteurs de translation
XY        = Ybox.ctr(Ibox(:,2),:) - Xbox.ctr(Ibox(:,1),:);
[~,Il,It] = unique(floor(XY*1e6),'rows','stable');
Nt        = length(Il);

% Fonctions de transfert
Tx = cell(Nt,1);
Ty = cell(Nt,1);
for i = 1:Nt
    [Tx{i},Ty{i}] = ffmTransfert(Xq,XY(Il(i),:),green,k,edg,tol);
end

% Interpolations en Y
Vy = cell(size(Ybox.ind));
for i = unique(Ibox(:,2)')
    % Boite centree en Y dans le cube unitaire d'interpolation
    iy = Ybox.ind{i};
    ny = length(iy);
    Ym = (2/edg) .* (Y(iy,:) - ones(ny,1)*Ybox.ctr(i,:));

    % Interpolation de Lagrange (Ym->Xq)
    Vy{i} = ffmInterp(ii,jj,kk,xq,Ym,V(iy),-1);
end

% Translations
TVy = cell(size(Xbox.ind));
for i = 1:length(TVy)
    TVy{i} = 0;
end
for i = 1:size(Ibox,1)
    TVy{Ibox(i,1)} = TVy{Ibox(i,1)} + Tx{It(i)} * (Ty{It(i)} * Vy{Ibox(i,2)});
end

% Interpolations en X
for i = unique(Ibox(:,1)')
    % Boite centree en X dans le cube unitaire d'interpolation
    ix = Xbox.ind{i};
    nx = length(ix);
    Xm = (2/edg) .* (X(ix,:) - ones(nx,1)*Xbox.ctr(i,:));

    % Interpolation de Lagrange (Xq->Xm)
    MV(ix) = ffmInterp(ii,jj,kk,xq,Xm,TVy{i},+1);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Xq,ii,jj,kk,xq] = ffmQuadrature(green,k,edg,tol)
% Nuages de reference emetteurs X
x       = (-0.5*edg:edg/(10-1):0.5*edg)';
[x,y,z] = ndgrid(x,x,x);
Xref    = [x(:) y(:) z(:)];
Nref    = size(Xref,1);

% Nuage de reference transmetteur avec translation minimale
r0   = [2*edg 0 0];
Yref = [Xref(:,1)+r0(1) , Xref(:,2:3)];

% Xref dans le cube unitaire d'interpolation
a    = [-0.5 -0.5 -0.5]*edg;
b    = [0.5 0.5 0.5]*edg;
Xuni = ones(Nref,1)*(2./(b-a)) .* (Xref-ones(Nref,1)*(b+a)./2);

% Solution de reference
V     = ones(Nref,1) + 1i;
[I,J] = ndgrid(1:Nref,1:1:Nref);
ref   = reshape(ffmGreenKernel(Xref(I,:),Yref(J,:),green,k),Nref,Nref) * V;

% Initialisation
sol = 1e6;
nq  = 1;

% Boucle sur l'ordre de quadrature
while norm(ref-sol)/norm(ref) > tol
    % Incrementation
    nq = nq + 1;

    % Indices de construction de l'interpolateur de Lagrange
    [ii,jj,kk] = ndgrid(1:nq,1:nq,1:nq);

    % Points de quadratures Tchebitchev (1D et 3D)
    xq = cos( (2*(nq:-1:1)-1)*pi/(2*nq))';
    Xq = [xq(ii(:)) xq(jj(:)) xq(kk(:))];

    % Vecteur de la translation
    [TA,TB] = ffmTransfert(Xq,r0,green,k,edg,tol);

    % Interpolation de Lagrange (Ym->Yq)
    Vy = ffmInterp(ii,jj,kk,xq,Xuni,V,-1);

    % Transfert (Yq->Xq)
    TVy = TA * (TB * Vy);

    % Interpolation de Lagrange (Xq->Xm)
    sol = ffmInterp(ii,jj,kk,xq,Xuni,TVy,+1);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [TA,TB] = ffmTransfert(Xq,xy,green,k,edg,tol)
% Dimensions
Nq = size(Xq,1);

% Interpolants en X
a   = [0 0 0];
b   = [1 1 1]*edg;
Txq = (ones(Nq,1)*(b-a)/2).*Xq + ones(Nq,1)*(b+a)/2;

% Interpolants en Y
a   = xy + a;
b   = xy + b;
Tyq = (ones(Nq,1)*(b-a)/2).*Xq + ones(Nq,1)*(b+a)/2;

% Noyaux de green avec compression aux interpolants
[TA,TB,flag] = hmxACA(Txq,Tyq,green,k,tol/10);

% Calcul direct si necessaire
if ~flag
    [I,J] = ndgrid(1:Nq,1:1:Nq);
    TA    = reshape(ffmGreenKernel(Txq(I,:),Tyq(J,:),green,k),Nq,Nq);
    TB    = 1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function MV = ffmInterp(ii,jj,kk,xq,X,V,iflag)
% Dimensions
Nx = size(X,1);
nq = length(xq);

% Interpolateur de Lagrange par dimension d'espace
Aq = cell(1,nq);
for l = 1:nq^2
    if isempty(Aq{ii(l)})
        Aq{ii(l)} = ones(Nx,3,class(V));
    end
    if ii(l) ~= jj(l)
        Aq{ii(l)} = Aq{ii(l)} .* (X-xq(jj(l)))./(xq(ii(l))-xq(jj(l)));
    end
end

% Interpolation (X->Xq)
if (iflag == -1)
    MV = zeros(nq^3,1,class(V));
    for l = 1:nq^3
        MV(l) = (Aq{ii(l)}(:,1) .* Aq{jj(l)}(:,2) .* Aq{kk(l)}(:,3)).' * V;
    end

% Interpolation (Xq->X)
elseif (iflag == +1)
    MV = zeros(Nx,1,class(V));
    for l = 1:nq^3
        MV = MV + Aq{ii(l)}(:,1) .* Aq{jj(l)}(:,2) .* Aq{kk(l)}(:,3) .* V(l);
    end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function R = AngleAxis2RotationMatrix(angle_axis)

theta2 = dot(angle_axis,angle_axis);
if (theta2 > 0.0)
    % We want to be careful to only evaluate the square root if the
    % norm of the angle_axis vector is greater than zero. Otherwise
    % we get a division by zero.

    theta = sqrt(theta2);
    wx = angle_axis(1) / theta;
    wy = angle_axis(2) / theta;
    wz = angle_axis(3) / theta;

    costheta = cos(theta);
    sintheta = sin(theta);

    R(1+0, 1+0) =     costheta   + wx*wx*(1 -    costheta);
    R(1+1, 1+0) =  wz*sintheta   + wx*wy*(1 -    costheta);
    R(1+2, 1+0) = -wy*sintheta   + wx*wz*(1 -    costheta);
    R(1+0, 1+1) =  wx*wy*(1 - costheta)     - wz*sintheta;
    R(1+1, 1+1) =     costheta   + wy*wy*(1 -    costheta);
    R(1+2, 1+1) =  wx*sintheta   + wy*wz*(1 -    costheta);
    R(1+0, 1+2) =  wy*sintheta   + wx*wz*(1 -    costheta);
    R(1+1, 1+2) = -wx*sintheta   + wy*wz*(1 -    costheta);
    R(1+2, 1+2) =     costheta   + wz*wz*(1 -    costheta);
else
    % At zero, we switch to using the first order Taylor expansion.
    R(1+0, 1+0) =  1;
    R(1+1, 1+0) = -angle_axis(3);
    R(1+2, 1+0) =  angle_axis(2);
    R(1+0, 1+1) =  angle_axis(3);
    R(1+1, 1+1) =  1;
    R(1+2, 1+1) = -angle_axis(1);
    R(1+0, 1+2) = -angle_axis(2);
    R(1+1, 1+2) =  angle_axis(1);
    R(1+2, 1+2) = 1;
end

author: jianxiongxiao

repo: ProfXkit


text:
function linplus_test426 ( )

%*****************************************************************************80
%
%% TEST426 tests R8CC_PRINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  m = 5;
  n = 5;
  nz_num = 12;
  colptr = [ 1, 4, 6, 8, 10, 13 ];
  rowind = [ 1, 2, 4, 1, 2, 3, 5, 4, 5, 1, 2, 5 ];
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST426\n' );
  fprintf ( 1, '  R8CC_PRINT prints a R8CC matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =          %d\n', m );
  fprintf ( 1, '  Matrix columns N =       %d\n', n );
  fprintf ( 1, '  Matrix nonzeros NZ_NUM = %d\n', nz_num );
%
%  Set the matrix.
%
  [ a, seed ] = r8cc_random ( m, n, nz_num, colptr, rowind, seed );
%
%  Print the matrix.
%
  r8cc_print ( m, n, nz_num, colptr, rowind, a, '  The R8CC matrix:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [P,T] = cubic_flat_eval(C,tol)
  % CUBIC_FLAT_EVAL Recursively subdivide a cubic Bezier curve until each
  % segment spans a region of the curve that is locally flat up to a given
  % tolerance (i.e., this computes an adaptive refinement).
  %
  % [P,T] = cubic_flat_eval(C,tol)
  %
  % Inputs:
  %   C  4 by dim list of control points
  % Outputs:
  %   P  #P by dim list of evaluated points
  %   T  #T list of corresponding parameter values
  %
  % See also: cubic_eval, cubic_is_flat, cubic_split
  %
  if cubic_is_flat(C,tol)
    P = C([1 4],:);
    T = [0 1];
  else
    [C1,C2] = cubic_split(C,0.5);
    [P1,T1] = cubic_flat_eval(C1,tol);
    [P2,T2] = cubic_flat_eval(C2,tol);
    P = [P1;P2(2:end,:)];
    T = [T1*0.5;0.5+T2(2:end,:)];
  end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
 function ob = subsasgn(ob, sub, arg)
%function ob = subsasgn(ob, sub, arg)
% method for "ob.sub = arg" etc.

% trick: for struct.sub = object
% e.g., ob.arg.field = ob;
if isstruct(ob)
%	ob = subsasgn(ob, sub, arg); % no, causes recursion
%	ob = builtin('subsasgn', ob, sub, arg);	% samit basu dislikes this

	% recurse a.b.c.d... to avoid builtin:
	if streq(sub(1).type, '.')
		if length(sub) > 1
%			printf('descending to %s', sub(1).subs)
			ob.(sub(1).subs) = subsasgn(ob.(sub(1).subs), sub(2:end), arg);
		else
			ob.(sub(1).subs) = arg;
		end
	else
		error 'not done'
	end

% trick: just convert it to a structure and then back!
else
	cl = class(ob);
	ob = struct(ob);
	name = sub.subs;
	if isfield(ob.data, name)
%		ob.data.(name) = arg;
		% 2011-03-12 changed to handle ig.mask(1)=0
		ob.data = subsasgn(ob.data, sub, arg);
	elseif isfield(ob.meth, name)
		ob.meth.(name) = arg;
	else
		error(['"' name '"is not in data or meth'])
	end
%	ob = subsasgn(ob, sub, arg);
	ob = class(ob, cl);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%SAVE2PDF Saves a figure as a properly cropped pdf
%
%   save2pdf(pdfFileName,handle,dpi)
%
%   - pdfFileName: Destination to write the pdf to.
%   - handle:  (optional) Handle of the figure to write to a pdf.  If
%              omitted, the current figure is used.  Note that handles
%              are typically the figure number.
%   - dpi: (optional) Integer value of dots per inch (DPI).  Sets
%          resolution of output pdf.  Note that 150 dpi is the Matlab
%          default and this function's default, but 600 dpi is typical for
%          production-quality.
%
%   Saves figure as a pdf with margins cropped to match the figure size.

%   (c) Gabe Hoffmann, gabe.hoffmann@gmail.com
%   Written 8/30/2007
%   Revised 9/22/2007
%   Revised 1/14/2007

function save2pdf(pdfFileName,handle,dpi)

% Verify correct number of arguments
error(nargchk(0,3,nargin));

% If no handle is provided, use the current figure as default
if nargin<1
    [fileName,pathName] = uiputfile('*.pdf','Save to PDF file:');
    if fileName == 0; return; end
    pdfFileName = [pathName,fileName];
end
if nargin<2
    handle = gcf;
end
if nargin<3
    dpi = 150;
end

% Backup previous settings
prePaperType = get(handle,'PaperType');
prePaperUnits = get(handle,'PaperUnits');
preUnits = get(handle,'Units');
prePaperPosition = get(handle,'PaperPosition');
prePaperSize = get(handle,'PaperSize');

% Make changing paper type possible
set(handle,'PaperType','<custom>');

% Set units to all be the same
set(handle,'PaperUnits','inches');
set(handle,'Units','inches');

% Set the page size and position to match the figure's dimensions
paperPosition = get(handle,'PaperPosition');
position = get(handle,'Position');
set(handle,'PaperPosition',[0,0,position(3:4)]);
set(handle,'PaperSize',position(3:4));

% Save the pdf (this is the same method used by "saveas")
print(handle,'-dpdf',pdfFileName,sprintf('-r%d',dpi))

% Restore the previous settings
set(handle,'PaperType',prePaperType);
set(handle,'PaperUnits',prePaperUnits);
set(handle,'Units',preUnits);
set(handle,'PaperPosition',prePaperPosition);
set(handle,'PaperSize',prePaperSize);

author: Sable

repo: mcbench-benchmarks


text:
function [ball, labels] = imInscribedBall(lbl, varargin)
% Maximal ball inscribed in a 3D region.
%
%   BALL = imInscribedBall(IMG)
%   Computes the maximal ball inscribed in a given 3D particle, or
%   around each labeled particle in the input image.
%
%   BALL = imInscribedBall(IMG, LABELS)
%   Specify the labels for which the inscribed ball needs to be computed.
%   The result is a N-by-3 array with as many rows as the number of labels.
%
%   Examples
%   % Test with a discretized ball
%     img = discreteBall(1:100, 1:100, 1:100, [40 50 60 35]);
%     ball = imInscribedBall(img)
%     ball =
%         40    50    60    35
%
%   % Check with a an octant of ball
%     img = discreteBall(1:100, 1:100, 1:100, [90 90 90 80]);
%     img(91:end, :,:) = 0;
%     img(:, 91:end, :) = 0;
%     img(:, :, 91:end) = 0;
%     ball = imInscribedBall(img)
%     ball =
%         61    61    61    30
%
%   See also
%     drawSphere, imInscribedCircle, imInertiaEllipsoid
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2013-07-05,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2013 INRA - Cepia Software Platform.

% check if labels are specified
labels = [];
if ~isempty(varargin) && size(varargin{1}, 2) == 1
    labels = varargin{1};
end

% extract the set of labels, without the background
if isempty(labels)
    labels = imFindLabels(img);
end
nLabels = length(labels);

% allocate memory for result (3 coords + 1 radius)
ball = zeros(nLabels, 4);

for i = 1:nLabels
    % compute distance map from background
    distMap = bwdist(lbl ~= labels(i));

    % find value and position of the maximum
    [maxi, inds] = max(distMap(:));
    [yb, xb, zb] = ind2sub(size(distMap), inds);

    ball(i,:) = [xb yb zb maxi];
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function triangulation_test33 ( )

%*****************************************************************************80
%
%% TEST33 tests VORONOI_POLYGON_VERTICES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2006
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  neighbor_num = 4;
  node_num = 5;

  center = 5;
  neighbor_index = [ 1, 2, 3, 4 ];
  node_xy = [ ...
    0.0, 0.0; ...
    1.0, 0.0; ...
    1.0, 1.0; ...
    0.0, 1.0; ...
    0.5, 0.5 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST33\n' );
  fprintf ( 1, '  VORONOI_POLYGON_VERTICES computes the vertices of\n' );
  fprintf ( 1, '  a finite Voronoi polygon.\n' );

  v = voronoi_polygon_vertices ( center, neighbor_num, neighbor_index, ...
    node_num, node_xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Voronoi Polygon Vertex coordinates:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : neighbor_num
    fprintf ( 1, '  %4d  %14f  %14f\n', i, v(1:2,i) );
  end

  return
end


text:
function cluoverl(var1)
    %cluoverl option to display equivalent events or biggest events or clear cluster events
    %                           A.Allmann

    %

    global bgevent plot1_h plot2_h
    global equi %[IN]
    global cluscat backequi newclcat a
    global dplo1_h dplo2_h dplo3_h dep1 dep2 dep3
    global  file1 clu h5
    global ty stri2 strib
    global after_h fore_h main_h after_button fore_button
    global foresh aftersh mainsh calll66
    global mainfault main faults clus_button coastline
    global SizMenu TypMenu new
    ZG=ZmapGlobal.Data;

    if var1==1                            %hide biggest events
        set(plot1_h,'Visible','off')
    elseif var1==2                         %plot biggest evens
        if isempty(plot1_h)                 %first time
            plot1_h=plot(bgevent(:,1),bgevent(:,2),'xm');
            set(plot1_h,'MarkerSize',5)
            set(plot1_h,'LineWidth',2)
        else
            set(plot1_h,'Visible','on')     %show plot that already exists(biggest events)
        end
    elseif var1==3                    %plot equivalent events
        if isempty(plot2_h)
            plot2_h=plot(equi(:,1),equi(:,2),'xg');
            set(plot2_h,'MarkerSize',5)
            set(plot2_h,'LineWidth',2)
        else
            set(plot2_h,'Visible','on')
        end
    elseif var1==4                  %hide equivalent events
        set(plot2_h,'Visible','off')
    elseif var1==5                  %hide clustered events
        set(dplo1_h,'Visible','off')
        set(dplo2_h,'Visible','off')
        set(dplo3_h,'Visible','off')
    elseif var1==6                  %show clustered events
        set(dplo1_h,'Visible','on')
        set(dplo2_h,'Visible','on')
        set(dplo3_h,'Visible','on')

    elseif var1==7                  %plot clusters and faults for the first time
        set(clus_button,'Value',1)
        if isempty(newclcat)
            replaceMainCatalog(cluscat);
        else
            replaceMainCatalog(newclcat);
        end
        cla
        set(gca,'Visible','off')
        set(gca,'NextPlot','replace')
        minde = 0.;
        maxde = max(ZG.primeCatalog.Depth);

        dep1 = round(0.333*maxde);
        dep2 = round(0.666*maxde);
        dep3 = maxde;

        stri1 = [file1];

        % find min and Maximum axes points
        s1_east = max(ZG.primeCatalog.Longitude);
        s2_west = min(ZG.primeCatalog.Longitude);
        s3_north = max(ZG.primeCatalog.Latitude);
        s4_south = min(ZG.primeCatalog.Latitude);
        ni = ZG.ni;
        orient landscape
        rect = [0.15,  0.12, 0.75, 0.75];
        axes('position',rect)
        %
        % find start and end time of catalogue "primeCatalog"
        %
        [t0b, teb] = bounds(ZG.primeCatalog.Date) ;
        n = ZG.primeCatalog.Count;
        tdiff = round(teb - t0b)/days(ZG.bin_dur);

        n = ZG.primeCatalog.Count;

        % plot earthquakes (differnt colors for varous depth layers) as
        % defined in "startzmap"
        %
        set(gca,'NextPlot','add')
        dplo1_h =plot(...
            ZG.primeCatalog.Longitude(ZG.primeCatalog.Depth<=dep1),...
            ZG.primeCatalog.Latitude(ZG.primeCatalog.Depth<=dep1),'.b');
        set(dplo1_h,'MarkerSize',ZG.ms6,'Marker',ty)
        dplo2_h =plot(...
            ZG.primeCatalog.Longitude(ZG.primeCatalog.Depth<=dep2&ZG.primeCatalog.Depth>dep1),...
            ZG.primeCatalog.Latitude(ZG.primeCatalog.Depth<=dep2&ZG.primeCatalog.Depth>dep1),'.y');
        set(dplo2_h,'MarkerSize',ZG.ms6,'Marker',ty);
        dplo3_h =plot(...
            ZG.primeCatalog.Longitude(ZG.primeCatalog.Depth<=dep3&ZG.primeCatalog.Depth>dep2),...
            ZG.primeCatalog.Latitude(ZG.primeCatalog.Depth<=dep3&ZG.primeCatalog.Depth>dep2),'.r');
        set(dplo3_h,'MarkerSize',ZG.ms6,'Marker',ty)

        axis([ s2_west s1_east s4_south s3_north])
        xlabel('Longitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        ylabel('Latitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        if isempty(backequi)
            strib = [  ' Map of   '  file1 ];
            %ti2 =  title(strib,'FontWeight','bold',...
            %             'FontSize',ZmapGlobal.Data.fontsz.l,'Color','r')
        else
            delete ti2;
        end
        %make depth legend
        %
        s = sprintf('Depth < %3.1f km',dep1);
        text('Color','b','units','normalized', 'Position',[0.05 0.15 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m,'String',s);
        s = sprintf('Depth < %3.1f km',dep2);
        text('Color','y','units','normalized', 'Position',[0.05 0.10 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m,'String',s);
        s = sprintf('Depth < %3.1f km',dep3);
        text('Color','r','units','normalized', 'Position',[0.05 0.05 0 ],...
            'FontSize',ZmapGlobal.Data.fontsz.m,'String',s);

        %  h5 is the graphic handle to the main figure in window 1
        %
        h5 = gca;
        set(gca,'box','on',...
            'SortMethod','childorder','TickDir','out','FontWeight',...
            'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

        %
        %  Plots epicenters  and faults

        overlay

        % Make the figure visible
        %
    elseif var1==8
        replaceMainCatalog(ZG.ttcat);
        cla
        set(gca,'visible','off')
        set(gca,'NextPlot','replace')
        s1_east = max(ZG.primeCatalog.Longitude);       %limits for area plot
        s2_west = min(ZG.primeCatalog.Longitude);
        s3_north = max(ZG.primeCatalog.Latitude);
        s4_south = min(ZG.primeCatalog.Latitude);
        if s1_east==s2_west                %to avoid error when all earthquakes have
            s1_east=s1_east+0.05;              %same longitude and/or latitude
            s2_west=s2_west-0.05;
        end
        if s3_north==s4_south
            s3_north=s3_north+0.05;
            s4_south=s4_south-0.05;
        end
        orient landscape
        rect = [0.15,  0.12, 0.75, 0.75];
        axes('position',rect)
        %
        % find start and end time of catalogue "primeCatalog"
        %
        [t0b, teb] = bounds(ZG.primeCatalog.Date) ;
        n = ZG.primeCatalog.Count;
        tdiff = round(teb - t0b)/days(ZG.bin_dur);
        %define fore and aftershocks
        %
        tmp = find(ZG.primeCatalog.Magnitude==max(ZG.primeCatalog.Magnitude));     %index in a of first mainshock
        if length(tmp)>1
            tmp=tmp(1,1);
        end
        overlay;
        %plot fore and aftershocks in different colors
        %
        set(gca,'NextPlot','add')
        if tmp-1>=1
            fore_h=plot(ZG.primeCatalogLongitude(1:tmp-1),ZG.primeCatalog.Latitude(1:tmp-1),'.b');
            if isempty(aftersh)                           %only at first call
                foresh=aZG.primeCatalog.subset(1:tmp-1);
            end
        else
            if exist('fore_h')
                fore_h=[];
            end
        end
        main_h=plot(ZG.primeCatalog.Longitude(tmp),ZG.primeCatalog.Latitude(tmp),'xm');
        mainsh=ZG.primeCatalog.subset(tmp);
        set(main_h,'MarkerSize',10);
        set(main_h,'LineWidth',2);
        if tmp+1<=n
            after_h=plot(ZG.primeCatalog.Longitude(tmp+1:n),ZG.primeCatalog.Latitude(tmp+1:n),'.r');
            if isempty(aftersh)
                aftersh=ZG.primeCatalog.subset(tmp+1:n);
            end
        else
            if exist('after_h')
                after_h=[];
            end
        end
        set(after_button,'value',1);
        set(fore_button,'value',1);
        axis([ s2_west s1_east s4_south s3_north])
        xlabel('Longitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        ylabel('Latitude [deg]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
        strib = [  ' Map of   '  file1 ' #' num2str(new(10))];
        title(strib,'FontWeight','bold',...
            'FontSize',ZmapGlobal.Data.fontsz.l,'Color','r')
        set(gca,'box','on',...
            'SortMethod','childorder','TickDir','out','FontWeight',...
            'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)

        set(gca,'visible','on')
    end

author: CelsoReyes


text:
function calpak_test010 ( )

%*****************************************************************************80
%
%% CALPAK_TEST010 tests JED_TO_YJF_HEBREW and YJF_TO_JED_HEBREW.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST010\n' );
  fprintf ( 1, '  For the Hebrew calendar:\n' );
  fprintf ( 1, '  JED_TO_YJF_HEBREW: JED -> YJF.\n' );
  fprintf ( 1, '  YJF_TO_JED_HEBREW: YJF -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  JED (in)    YJF                JED (out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_hebrew ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ y2, j2, f2 ] = jed_to_yjf_hebrew ( jed1 );

      s2 = yjf_to_s_hebrew ( y2, j2, f2 );

      jed3 = yjf_to_jed_hebrew ( y2, j2, f2 );

      fprintf ( 1, '  %11.2f  %20s  %11.2f\n', jed1, s2, jed3 );

    end

  end

  return
end


text:
function [version, released] = manopt_version()
% Returns the version of the Manopt package you are running, as a vector.
%
% function [version, released] = manopt_version()
%
% version(1) is the primary version number.
% released is the date this version was released, in the same format as the
% date() function in Matlab.

    version = [2, 0, 0];
    released = '06-Jul-2015';

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function obj = updateBoxplot(obj, boxIndex)

% y: ...[DONE]
% x0: ...[DONE]
% x: ...[DONE]
% name: ...[DONE]
% boxmean: ...[NOT SUPPORTED IN MATLAB]
% boxpoints: ...[NOT SUPPORTED IN MATLAB]
% jitter: ...[NOT SUPPORTED IN MATLAB]
% pointpos: ...[NOT SUPPORTED IN MATLAB]
% whiskerwidth: ........................[TODO]
% fillcolor: ...[DONE]
% opacity: ---[TODO]
% xaxis: ...[DONE]
% yaxis: ...[DONE]
% showlegend: ...[DONE]
% stream: ...[HANDLED BY PLOTLY STREAM]
% visible: ...[DONE]
% type: ...[DONE]

% MARKER
% color: ...[NA]
% width: ...[NA]
% dash: ...[NA]
% opacity: ...[NA]
% shape: ...[NA]
% smoothing: ...[NA]
% outliercolor: ...[NOT SUPPORTED IN MATLAB]
% outlierwidth: ...[NOT SUPPORTED IN MATLAB]

% LINE
% color: ...[DONE]
% width: ...[DONE]
% dash: ...[DONE]
% opacity: ---[TODO]
% shape: ...[DONE]
% smoothing: ...[NOT SUPPORTED IN MATLAB]

%-------------------------------------------------------------------------%

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(boxIndex).AssociatedAxis);

%-BOX DATA STRUCTURE-%
box_data = get(obj.State.Plot(boxIndex).Handle);

%-BOX CHILDREN-%
box_child = box_data.Children;

%-------------------------------------------------------------------------%

%-CONFIRM PROPER BOXPLOT STRUCTURE-%

% check for compact boxplot
isCompact = ~isempty(findobj(obj.State.Plot(boxIndex).Handle,'Tag','Whisker'));

% number of boxplots
if isCompact
    bpcompnum = 6;
    bpnum = length(box_child)/bpcompnum;
    % check for assumed box structure
    if mod(length(box_child), bpcompnum)~=0
        return
    end
else
    bpcompnum = 8;
    bpnum = length(box_child)/bpcompnum;
    % check for assumed box structure
    if mod(length(box_child),bpcompnum)~=0
        return
    end
end

% initialize ydata
ydata = [];

%-------------------------------------------------------------------------%

%-box groupgap-%
obj.layout.bargroupgap = 1/bpnum;

%-------------------------------------------------------------------------%

%-box name-%
obj.data{boxIndex}.name = box_data.DisplayName;

%-------------------------------------------------------------------------%

% iterate through box plot children in reverse order
for bp = bpnum:-1:1

    %-CHECK FOR MULTIPLE AXES-%
    [xsource, ysource] = findSourceAxis(obj,axIndex);

    %-AXIS DATA-%
    eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
    eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

    %---------------------------------------------------------------------%

    %-box xaxis-%
    obj.data{boxIndex}.xaxis = ['x' num2str(xsource)];

    %---------------------------------------------------------------------%

    %-box yaxis-%
    obj.data{boxIndex}.yaxis = ['y' num2str(ysource)];

    %---------------------------------------------------------------------%

    %-box type-%
    obj.data{boxIndex}.type = 'box';

    %---------------------------------------------------------------------%

    %-box visible-%
    obj.data{boxIndex}.visible = strcmp(box_data.Visible,'on');

    %---------------------------------------------------------------------%

    %-box fillcolor-%
    obj.data{boxIndex}.fillcolor = 'rgba(0, 0, 0, 0)';

    %---------------------------------------------------------------------%

    %-box showlegend-%
    leg = get(box_data.Annotation);
    legInfo = get(leg.LegendInformation);

    switch legInfo.IconDisplayStyle
        case 'on'
            showleg = true;
        case 'off'
            showleg = false;
    end

    obj.data{boxIndex}.showlegend = showleg;

    %-boxplot components-%
    Q1 = [];
    Q3 = [];
    median = [];
    outliers = [];
    uwhisker = [];
    lwhisker = [];

    % iterate through boxplot components
    for bpc = 1:bpcompnum

        %get box child data
        box_child_data = get(box_child(bp+bpnum*(bpc-1)));

        %box name
        if strcmp(box_child_data.Type,'text')
            if iscell(box_child_data.String)
                boxname =  box_child_data.String{1};
            else
                boxname =  box_child_data.String;
            end
        end

        % parse boxplot tags
        switch box_child_data.Tag

            %-median-%
            case 'Median'

                median =  box_child_data.YData(1);

                %-upper whisker-%
            case 'Upper Whisker'

                uwhisker = box_child_data.YData(2);

                %-boxplot whisker width-%
                obj.data{boxIndex}.whiskerwidth = 1;

                %-lower whisker-%
            case 'Lower Whisker'

                lwhisker = box_child_data.YData(1);

            case 'Box'

                %-Q1-%
                Q1 = min(box_child_data.YData);

                %-Q3-%
                Q3 = max(box_child_data.YData);

                %-boxplot line style-%

                if isCompact
                    col = 255*box_child_data.Color;
                    obj.data{boxIndex}.fillcolor = ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')'];
                else
                    obj.data{boxIndex}.line = extractLineLine(box_child_data);
                end

                %-outliers-%
            case 'Outliers'

                if ~isnan(box_child_data.YData)
                    %-outlier marker data-%
                    outliers = box_child_data.YData;

                    %-outlier marker style-%
                    obj.data{boxIndex}.marker = extractLineMarker(box_child_data);
                end

                %-compact whiskers-%
            case 'Whisker'

                %-boxplot line style-%
                obj.data{boxIndex}.line = extractLineLine(box_child_data);

                %-boxplot whisker width-%
                obj.data{boxIndex}.whiskerwidth = 0;

                %-whisker data-%
                uwhisker = box_child_data.YData(2);
                lwhisker = box_child_data.YData(1);

                %-compact median-%
            case 'MedianInner'

                median = box_child_data.YData(1);

        end
    end

    %-generate boxplot data-%
    gendata = generateBoxData(outliers, lwhisker, Q1, median, Q3, uwhisker);

    %-boxplot y-data-%
    obj.data{boxIndex}.y(length(ydata)+1:length(ydata)+length(gendata)) = ...
        generateBoxData(outliers, lwhisker, Q1, median, Q3, uwhisker);

    %-boxplot x-data-%
    if (bpnum > 1)
        for n = (length(ydata)+1):(length(ydata)+length(gendata))
            obj.data{boxIndex}.x{n} = boxname;
        end
    end

    %-update ydata-%
    ydata = obj.data{boxIndex}.y;

end

%----------------------------!AXIS UPDATE!--------------------------------%

% take first text object as prototype for axis tick style/layout
text_child = findobj(obj.State.Plot(boxIndex).Handle,'Type','text');

%-STANDARDIZE UNITS-%
fontunits = get(text_child(1),'FontUnits');
set(text_child(1),'FontUnits','points');

%-text data -%
text_data = get(text_child(1));

%-------------------------------------------------------------------------%

%-xaxis tick font size-%
xaxis.tickfont.size = text_data.FontSize;

%-------------------------------------------------------------------------%

%-xaxis tick font family-%
xaxis.tickfont.family = matlab2plotlyfont(text_data.FontName);

%-------------------------------------------------------------------------%

%-xaxis tick font color-%
xaxis.tickfont.color = text_data.Color;

%-------------------------------------------------------------------------%

%-axis type-%
xaxis.type = 'category';

%-------------------------------------------------------------------------%

%-show tick labels-%
xaxis.showticklabels = true;

%-------------------------------------------------------------------------%

%-autorange-%
xaxis.autorange = true;

%-------------------------------------------------------------------------%

%-set the layout axis field-%
obj.layout = setfield(obj.layout,['xaxis' num2str(xsource)],xaxis);

%-------------------------------------------------------------------------%

%-REVERT UNITS-%
set(text_child(1),'FontUnits',fontunits);

end

author: goGPS-Project


text:
function a = vand1_inverse ( n, x )

%*****************************************************************************80
%
%% VAND1_INVERSE returns the inverse of the VAND1 matrix.
%
%  Formula:
%
%    A(I,J) = coefficient of X^(J-1) in I-th Lagrange basis polynomial.
%
%  Example:
%
%    N = 5,
%    X = ( 2, 3, 4, 5, 6 )
%
%     15.00  -14.25    4.96  -0.75   0.04
%    -40.00   44.67  -17.33   2.83  -0.17
%     45.00  -54.00   22.75  -4.00   0.25
%    -24.00   30.00  -13.33   2.50  -0.17
%      5.00   -6.42    2.96  -0.58   0.04
%
%  Properties:
%
%    The sum of the entries of A is
%
%      1 - product ( 1 <= I <= N ) ( 1 - 1 / X(I) ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real X(N), the values that define A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n
      if ( j == 1 )
        a(i,j) = 1.0;
      else
        a(i,j) = 0.0;
      end
    end
  end

  for i = 1 : n

    index = 0;

    for k = 1 : n

      if ( k ~= i )

        index = index + 1;

        for j = index + 1: - 1 : 1

          a(i,j) = - x(k) * a(i,j) / ( x(i) - x(k) );

          if ( 1 < j )
            a(i,j) = a(i,j) + a(i,j-1) / ( x(i) - x(k) );
          end

        end

      end

    end

  end

  return
end


text:
function [xsecx,xsecy] = lc_xsec2(eqlat,eqlon,depth,width,length,...
        lat1,lon1,lat2,lon2)

    %LC_XSEC2 make a cross section of data points on a map
    % Based off of LC_XSECTION
    %
    %	[xsecx, xsecy] = LC_xsec2(eqlat,eqlon,depth,width,length,...
    %                                        lat1,lon1,lat2,lon2)        (1)
    %
    %	[xsecx, xsecy] = LC_xsec2(eqlat,eqlon,depth,width,length,...
    %                                        lat0,lon0,azimuth)          (2)
    %
    %	[xsecx, xsecy] = LC_xsec2(eqlat,eqlon,depth,width)         (3)
    %
    %	Function to make a cross section of data points on a map
    %	created by LC_MAP (Lambert Conformal).
    %	The WIDTH of the zone from which the data points is given
    %	in "km" and represent the total width (1/2 on one side, 1/2 on
    %	the other side).
    %	The LENGTH of the xsection in "km" is only used in method (2),
    %	but this argument is still neccessary in argument list of
    %	method (1) in order to keep the argument list in the right order.
    %	The data points location are given by EQLAT, EQLON and DEPTH.
    %	The cross section location can be given in any one of three ways:
    %
    %	  (1) given latitudes and longitudes of two points on the map,
    %             using the arguments as described above.
    %
    %	  (2) given latitude and longitude of a center point and an azimuth.
    %
    %	  (3) using the cursor to select two points on the map by clicking
    %	      a mouse button above the desired points.
    %
    %	If the output argument is used, the distance-depth data is kept
    %	in variables for other use; otherwise the xsection will be plotted
    %	on a new figure window.
    %
    %	It is possible to set the symbol type, size and line width by
    %	setting the following global variables:
    %	"symb_type", "symb_size" and "symb_width" respectively. Otherwise
    %	it will use the defaults.
    %
    %	It is also possible to set the minimum and maximum depth of the
    %	cross-section by setting the following global variables:
    %	"mindepth" and "maxdepth".  If either or both are not set, it
    %	will use 0 km as the minimum depth and/or the depth of the deepest
    %	data point as the maximum depth.
    %
    %	NOTE:
    %	It is assumed that LC_MAP was used before using this function!
    %	This is neccessary to set global variables used by this function.

    global sine_phi0 lambda0
    global maxlatg minlatg maxlong minlong
    global symb_type symb_size symb_width
    global label1 label2
    global mindepth maxdepth

    report_this_filefun();

    ZG = ZmapGlobal.Data;
    todeg = 180 / pi;
    torad = ZG.torad;

    if nargin < 9

        if nargin == 8	% method 2: given lat & lon of center point and angle

            lat0 = lat1;
            lon0 = lon1;
            [x0, y0] = lc_tocart(lat0,lon0);
            azimuth = lat2;

            if azimuth >= 180, azimuth = azimuth - 180; end
            theta0 = ((lon0*torad - lambda0) * sine_phi0) * todeg;
            alpha = azimuth - theta0;
            beta = (90 - azimuth) + theta0;

            x2 = ((length / 2) * cos(beta*torad));
            y2 = ((length / 2) * sin(beta*torad));
            x1 = x0 - x2;
            y1 = y0 - y2;
            x2 = x0 + x2;
            y2 = y0 + y2;

            [lat1, lon1] = lc_froca(x1,y1);
            [lat2, lon2] = lc_froca(x2,y2);

        elseif nargin == 4	% method 3: selection of the end points by mouse

            limits = ginput(2);
            x1 = limits(1,1);
            y1 = limits(1,2);
            x2 = limits(2,1);
            y2 = limits(2,2);

            if x1 > x2
                xtemp = x1; ytemp = y1;
                x1 = x2; y1 = y2;
                x2 = xtemp; y2 = ytemp;
            end

            [lat1, lon1] = lc_froca(x1,y1);
            [lat2, lon2] = lc_froca(x2,y2);

            x0 = (x1 + x2) / 2;
            y0 = (y1 + y2) / 2;
            [lat0, lon0] = lc_froca(x0,y0);
            dx = x2 - x1;
            dy = y2 - y1;

            alpha = 90 - (atan(dy/dx)*todeg);
            length = sqrt(dx^2 + dy^2);

        else
            disp('ERROR: incompatible number of arguments')
            help lc_xsec2
            return
        end

    elseif nargin == 9	% method 1: given lat & lon of the two end points

        [x1, y1] = lc_tocart(lat1,lon1);
        [x2, y2] = lc_tocart(lat2,lon2);

        if x1 > x2
            xtemp = x1; ytemp = y1;
            x1 = x2; y1 = y2;
            x2 = xtemp; y2 = ytemp;
        end

        x0 = (x1 + x2) / 2;
        y0 = (y1 + y2) / 2;
        [lat0, lon0] = lc_froca(x0,y0);
        dx = x2 - x1;
        dy = y2 - y1;

        alpha = 90 - (atan(dy/dx)*todeg);
        length = sqrt(dx^2 + dy^2);

    else

        disp('ERROR: incompatible number of arguments')
        help lc_xsec2
        return
    end

    % correction factor to correct for longitude away from the center meridian
    theta0 = ((lon0*torad - lambda0) * sine_phi0) * todeg;

    % correct the XY azimuth of the Xsection line with the above factor to obtain
    % the true azimuth
    azimuth = alpha + theta0;
    if azimuth < 0, azimuth = azimuth + 180; end

    % convert XY coordinate azimuth to a normal angle like we used to deal with
    sigma = 90 - alpha;

    % transformation matrix to rotate the data coordinate w.r.t the Xsection line
    transf = [cos(sigma*torad) sin(sigma*torad)
        -sin(sigma*torad) cos(sigma*torad)];

    % inverse transformation matrix to rotate the data coordinate back
    invtransf = [cos(-sigma*torad) sin(-sigma*torad)
        -sin(-sigma*torad) cos(-sigma*torad)];

    % convert the map coordinate of the events to cartesian coordinates
    idx_map = find(minlatg < eqlat & eqlat < maxlatg & ...
        minlong < eqlon & eqlon < maxlong);
    [eq(1,:) eq(2,:)] = lc_tocart(eqlat,eqlon);

    % create new coordinate system at center of Xsection line
    eq0(1,:) = eq(1,:) - x0;
    eq0(2,:) = eq(2,:) - y0;

    % rotate this last coordinate system so that X-axis correspond to Xsection line
    eq0p = transf * eq0;

    % project the event data to the Xsection line
    eq1(1,:) = eq0p(1,:);
    eq1(2,:) = eq0p(2,:) * 0;

    % convert back to the original coordinate system
    eq1p = invtransf * eq1;
    eq2(1,:) = eq1p(1,:) + x0;
    eq2(2,:) = eq1p(2,:) + y0;

    % find index of all events which are within the given box width
    idx_box = find(abs(eq0p(2,:)) <= width/2 & abs(eq0p(1,:)) <= length/2);

    % Plot events on cross section figure
    xdist = eq1(1,idx_box) + (length / 2);

    % check if need to plot data (nargout = 0), or need to transfer data back
    % (nargout > 0)

    if nargout > 0
        xsecx = xdist;
        xsecy = depth(idx_box);

    else
        % plot the Xsection line on the map
        plot([x1 x2],[y1 y2],'--','LineWidth',1.5)

        % label the Xsection end points
        xlim = get(gca,'XLim');
        ylim = get(gca,'YLim');
        label_dist = (2 / 100) * sqrt((2*xlim(2))^2 + (2*ylim(2))^2);
        label_pt(1,1) = -(length/2 + label_dist);
        label_pt(2,1) = 0;
        label_pt(1,2) = length/2 + label_dist;
        label_pt(2,2) = 0;
        rlabel_pt = invtransf * label_pt;
        label_pt(1,:) = rlabel_pt(1,:) + x0;
        label_pt(2,:) = rlabel_pt(2,:) + y0;
        if isempty(label1), label1 = 'A'; end
        if isempty(label2), label2 = 'A'''; end
        lbl1_h = text(label_pt(1,1),label_pt(2,1),label1,'FontSize',14,...
            'Vertical','middle','Horizontal','center','FontWeight','bold');
        lbl2_h = text(label_pt(1,2),label_pt(2,2),label2,'FontSize',14,...
            'Vertical','middle','Horizontal','center','FontWeight','bold');

        % create a box of width "width" around the Xsection line and plot it
        box(1,1) = -length/2; box(2,1) = width/2;
        box(1,2) = length/2; box(2,2) = width/2;
        box(1,3) = length/2; box(2,3) = -width/2;
        box(1,4) = -length/2; box(2,4) = -width/2;
        xbox = [box(1,:) box(1,1)];
        ybox = [box(2,:) box(2,1)];
        rbox = invtransf * [xbox ; ybox];
        rbox(1,:) = rbox(1,:) + x0;
        rbox(2,:) = rbox(2,:) + y0;
        plot(rbox(1,:),rbox(2,:),'-.','LineWidth',0.8)

        % check if symbol parameters global variables are set, if not --> defaults
        if isempty(symb_type), symb_type = '+'; end
        if isempty(symb_size), symb_size = 6; end
        if isempty(symb_width), symb_width = [0.5]; end

        % plot the events on the map
        plot(eq(1,idx_map),eq(2,idx_map),symb_type,'MarkerSize',symb_size,...
            'LineWidth',symb_width)

        % Open another graphic window for the cross section
        figure(xsec_fig());
        set(xsec_fig,'PaperPosition',[1 .5 9 6.9545])

        %global Xwbz Ywbz
        %Xwbz = xdist;
        %Ywbz = depth(idx_box);

        plot(xdist,-depth(idx_box),symb_type,'MarkerSize',symb_size,...
            'LineWidth',symb_width)

        if isempty(maxdepth)
            maxZ = max(depth(idx_box));
        else
            maxZ = maxdepth;
        end

        if isempty(mindepth)
            minZ = 0;
        else
            minZ = mindepth;
        end

        if length > (maxZ - minZ)*11/8.5
            position = [.1 .1 .7 ((maxZ-minZ)/length)*0.7*11/8.5];
        else
            position = [.1 .1 (length/(maxZ-minZ))*0.7*8.5/11 .7];
        end
        set(gca,'Position',position,'XLim',[0 length],'Ylim',[-maxZ -minZ],...
            'LineWidth',2)

        % Plot labels
        Xstring = ['Distance from ' label1 ' (km)'];
        set(gca,'XLabel',text(0,0,Xstring),'YLabel',text(0,0,'Depth (km)'))
        label_base1 = 1 + .04;
        label_base2 = 1 + .06;
        label_base3 = 1 + .08;
        lbl3_h = text(0,label_base2,label1,'FontSize',14,'Horizontal','center',...
            'FontWeight','bold','Vertical','middle','Units','norm');
        lat1_dm = sprintf('    %2.2i N %4.2f''',fix(lat1),(frac(lat1)*60));
        lon1_dm = sprintf('   %3.3i W %4.2f''',abs(fix(lon1)),abs(frac(lon1)*60));
        lbl5_h = text(0,label_base1,lat1_dm,'FontSize',12,'Horizontal','left',...
            'Vertical','bottom','Units','norm');
        lbl6_h = text(0,label_base3,lon1_dm,'FontSize',12,'Horizontal','left',...
            'Vertical','bottom','Units','norm');

        lbl4_h = text(1,label_base2,label2,'FontSize',14,'Horizontal','center',...
            'FontWeight','bold','Vertical','middle','Units','norm');
        lat2_dm = sprintf('%2.2i N %4.2f''    ',fix(lat2),(frac(lat2)*60));
        lon2_dm = sprintf('%3.3i W %4.2f''    ',abs(fix(lon2)),abs(frac(lon2)*60));
        lbl7_h = text(1,label_base1,lat2_dm,'FontSize',12,...
            'Horizontal','right','Vertical','bottom','Units','norm');
        lbl8_h = text(1,label_base3,lon2_dm,'FontSize',12,...
            'Horizontal','right','Vertical','bottom','Units','norm');

        % Go back to map figure
        %figure(map_fig);

    end
end


text:
% 迭代法求解方程
clear;
format long;
tol = 1e-5;
N = 100;
x0 = 10.5;

phi = @(x) exp(-x);
for k = 1 : N
    x1 = phi(x0);
    if abs(x1 - x0) < tol
        fprintf('方程的正根: %10.8f\n', x1);
        fprintf('迭代次数: %d\n', k);
        break;
    end
    x0 = x1;
end
if k == N
    fprintf('迭代方法失败\n');
end

repo: NumericalAnalysis


text:
function s = hms2point(s)
% convert Herman – Mauguin (International) Symbol to point group
%
% Input
%  s - Herman – Mauguin (International) space group symbol
%
% Output
%  s - international point group symbol
%
% See also
%

% ignore first letter
assert(any(strcmp(s(1),{'A','B','C','F','I','P','R'})));
s = s(2:end);

% convert all planes to mirror planes
rep = {'a','b','c','d','n'};
for i = 1:length(rep)
  s = strrep(s,rep{i},'m');
end

% remove two consecutive numbers if the second is a divisor of the first
% eg. 63 62 42 32
if length(regexp(s,'[\dm]'))>=4
  ind = regexp(s,'\d\d');
  if ~isempty(ind)
    ind = ind(1);
    n1 = str2num(s(ind)); %#ok<ST2NM>
    n2 = str2num(s(ind+1)); %#ok<ST2NM>
    if n1 > n2, s(ind+1) = [];end
  end
end

% ignore 1
if ~isempty(regexp(s,'[m2-6]')), s = strrep(s,'1','');end

% clear all spaces
s = strrep(s,' ','');

% remove all upcase letter
ind = s==upper(s) & s~=lower(s);
s(ind) = [];

s= regexprep(s,'[:.]','');  % occurs sometimes

if isempty(s), s = '1';end

author: mtex-toolbox

repo: mtex

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


text:
function a = i4vec_sort_bubble_a ( n, a )

%*****************************************************************************80
%
%% I4VEC_SORT_BUBBLE_A ascending sorts an I4VEC using bubble sort.
%
%  Discussion:
%
%    Bubble sort is simple to program, but inefficient.  It should not
%    be used for large arrays.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, integer A(N), an unsorted array.
%
%    Output, integer A(N), the array has been sorted.
%
  for i = 1 : n-1
    for j = i+1 : n
      if ( a(j) < a(i) )
        t    = a(i);
        a(i) = a(j);
        a(j) = t;
      end
    end
  end

  return
end


text:
function [c, ceq] = rendezvousNonlconFunc(x, iniOrbit, finOrbit, gmuXfr, minPe, maxAp, onlyOptBurn1)
%rendezvousNonlconFunc Summary of this function goes here
%   Detailed explanation goes here
    [~, ~, ~, ~, ~, ~, xfrOrbit, ~, ~] = rendezvousObjFunc(x, iniOrbit, finOrbit, gmuXfr, [1,1], onlyOptBurn1);

	[rAp, rPe] = computeApogeePerigee(xfrOrbit(1), xfrOrbit(2));
    c(1) = rAp - maxAp;
    c(2) = minPe - rPe;
    ceq = [];
end


text:
function labels=dtiGetMoriLabels(short)

%ER 11/2009 wrote it

tdir = fullfile(fileparts(which('mrDiffusion.m')), 'templates');
labels = readTab(fullfile(tdir,'MNI_JHU_tracts_prob.txt'),',',false);
labels = labels(1:20,2);
if exist('short', 'var') && short
labels={'ATR_L', 'ATR_R', 'CST_L', 'CST_R', 'CCG_L', 'CCG_R', 'CHC_L','CHC_R', 'Forceps Major', 'Forceps Minor', 'IFO_L', 'IFO_R', 'ILF_L', 'ILF_R', 'SLF_L', 'SLF_R', 'UF_L', 'UF_R', 'SLF_t_L', 'SLF_t_R' }';
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function diffOddEven = compute_diff_odd_even(this, varargin)
% computes difference image between first and last volume of time series
% NOTE: short-cut for compute_stat_image('diff_last_first')
%
%   Y = MrImage()
%   diffOddEven = Y.compute_diff_odd_even('PropertyName', PropertyValue)
%
% This is a method of class MrImage.
%
% IN
%   'PropertyName'
%               'selectedVolumes'       [1,nVols] vector of selected
%                                       volumes for statistical calculation
% OUT
%   diffOddEven               MrImage holding difference image between
%                               odd and even
%                               time series volume, characterizing "image
%                               noise" as in FBIRN paper (Friedman and
%                               Glover, JMRI 2006)
% EXAMPLE
%   Y = MrImage()
%   diffOddEven = Y.compute_diff_odd_even('selectedVolumes', [6:100])
%
%   See also MrImage compute_stat_image

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2020-09-23
% Copyright (C) 2020 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

diffOddEven = this.compute_stat_image('diff_last_first', varargin{:});

author: translationalneuromodeling

repo: tapas


text:
function T = full(T)
%FULL Convert formatted data set to an array.
%   T = fmt(T) converts the formatted data set T into a MATLAB array. If T
%   is incomplete, unknown entries are represented by the value NaN.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if isnumeric(T), return; end
if ~isstruct(T)
    error('full:T','T must be a data set formatted by fmt.');
end

val = T.val;
size_tens = T.size;
if 8*prod(size_tens) > 8e9
    error('full:size','T is too large to be stored as an array.');
end
if ~isfield(T,'ind')
    ind = sub2ind(size_tens,T.sub{:});
else
    ind = T.ind;
end

if isfield(T,'sparse') && T.sparse
    T = zeros(size_tens);
elseif isfield(T,'incomplete') && T.incomplete
    T = nan(size_tens);
else
    error('full')
end
T(ind) = val;

author: andrewssobral

repo: lrslibrary


text:
function [filename, pathname] = getFilename(varargin)
% DESCRIPTION:
%       Function to automate collecting a single filename input using the
%       built in uigetfile gui.
%
% USAGE:
%       [filename pathname] = getFilename(start_dir)
%       [filename pathname] = getFilename()
%
% INPUTS:
%       start_dir   - the directory to start in
%
% OUTPUTS:
%       filename / pathname of selected file
%
% ABOUT:
%       author          - Bradley Treeby
%       date            - 5th November 2008
%       last update     - 30th June 2009
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% get the filename
if nargin == 1
    start_dir = varargin{1};
    [filename, pathname, filterindex] = uigetfile({'*.jpg; *.gif; *.bmp; *.png', 'Image Files (*.jpg, *.gif, *.bmp, *.png)'; '*.*',...
        'All Files (*.*)'}, 'Select data file', 'MultiSelect', 'off', start_dir);
else
    [filename, pathname, filterindex] = uigetfile({'*.jpg; *.gif; *.bmp; *.png', 'Image Files (*.jpg, *.gif, *.bmp, *.png)'; '*.*',...
        'All Files (*.*)'}, 'Select data file', 'MultiSelect', 'off');
end

% check and end with an error if no file has been selected
if filterindex == 0
    error('Please select a data file');
end

repo: Aero-matlab


text:
function [M, inliers] = RANSAC2(x, fittingfn, distfn, degenfn, s, t, index, validMask)

    maxTrials = 1000;
    maxDataTrials = 100;

    [rows, npts] = size(x);

    % Desired probability of choosing at least one sample free from outliers (probably should be a parameter)
    p = 0.99;

    bestM = NaN;
    trialcount = 0;
    bestscore =  0;
    N = 1;

    while N > trialcount

        % Select at random s datapoints to form a trial model, M.
        degenerate = 1;
        count = 1;
        while degenerate
            % Generate s random indicies in the range 1..npts
            ind = zeros(s, 1);
%             ind = randsample(npts, s);
            for i = 1:s
                pick = randsample(npts, 1);
                while validMask(index(pick)) == 0 || any(ind == pick) || any(index(ind(ind > 0)) == index(pick))
                    pick = randsample(npts, 1);
                end
                ind(i) = pick;
            end

            % Test that these points are not a degenerate configuration.
            degenerate = feval(degenfn, x(:,ind));

            if ~degenerate
                M = feval(fittingfn, x(:,ind));
                if isempty(M)
                    degenerate = 1;
                end
            end

            % Safeguard against being stuck in this loop forever
            count = count + 1;
            if count > maxDataTrials
                disp('Unable to select a nondegenerate data set');
                break
            end
        end

%         if scoreX1 ~= 0
%             l = min([scoreX1(ind) scoreX2(ind)]);
%             t = (l + 0.6) * t;
%         end

        [inliers, score, M] = feval(distfn, M, x, t, index, validMask);

        % Find the number of inliers to this model.
        ninliers = length(inliers);

        if score >= bestscore
            bestscore = score;
            bestinliers = inliers;
            bestM = M;

            % Update estimate of N
            fracinliers =  ninliers/npts;
            pNoOutliers = 1 -  fracinliers^s;
            pNoOutliers = max(eps, pNoOutliers);
            pNoOutliers = min(1-eps, pNoOutliers);
            N = log(1-p)/log(pNoOutliers);
        end

        trialcount = trialcount+1;

        % Safeguard against being stuck in this loop forever
        if trialcount > maxTrials
            fprintf('ransac reached the maximum number of %d trials\n', maxTrials);
            break
        end

        %fprintf('Iter: %d \n', trialcount);
    end

    if ~isnan(bestM)
        M = bestM;
        inliers = bestinliers;
    else
        M = [];
        inliers = [];
        disp('ransac was unable to find a useful solution');
    end
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function [cfg, newCfg ] = loaddefault(cfg, defaultCfg, appendStruct)
%loaddefault will be removed, use ezInput instead.
%
% keepStruct: recursive retain the struct
%LOAD DEFAULT CONFIG
% defaultCfg = struct('a',1,'b',2,'c',3);
% cfg = struct('b',0);
% [cfg, newCfg ] = loaddefault(cfg, defaultCfg)
%
% cfg = struct('b',0, 'd', 1, 'e', 2); % unknown paramvalue
% [cfg, newCfg ] = loaddefault(cfg, defaultCfg)

%     existingfields = isfield(cfg,fieldnames(defaultCfg));
    %if cfg
    if nargin < 3
        appendStruct = false;
    end

     if iscell(cfg)
        cfg = struct(cfg{:});
     end
    if iscell(defaultCfg)
        defaultCfg = struct(defaultCfg{:});
    end

    newCfg = cfg;

    % both are struct

    % defaultCfg --> cfg
    fields = fieldnames(defaultCfg);
    nfields = numel(fields);
    for n = 1:nfields
        if ~isfield(cfg,fields(n))
            cfg.(fields{n}) = defaultCfg.(fields{n});
        else
            if appendStruct && isstruct(defaultCfg.(fields{n}))
                cfg.(fields{n}) = loaddefault(cfg.(fields{n}), ...
                                              defaultCfg.(fields{n}) );
            end
            newCfg = rmfield(newCfg, fields(n));
        end
    end

    % newCfg --> cell
    if nargout > 1
        newCfg = convertToCell(newCfg);
    end
end

function c = convertToCell(s)
    twoRow = [fieldnames(s), struct2cell(s)]';
    c = twoRow(:)';
end

% scatter3_args = rmfield(param, default(1:2:end));
%

author: baidut

repo: BIMEF

save_path: github-repos/MATLAB/baidut-BIMEF


text:
function test_bug2986

% MEM 3gb
% WALLTIME 00:20:00
% DEPENDENCY ft_volumerealign ft_volumereslice

load standard_mri

% convert to 'cm' to show what happens
mri = ft_convert_units(mri, 'cm');

% reslice to standard (improves solutions)
mri =  ft_volumereslice([], mri);

% load headshape
load(fullfile(dccnpath('/home/common/matlab/fieldtrip/data/test'),'bug2986.mat'));

% for quick check feed in fiducial positions
fiducial = [];
fiducial.nas = [131 215 85];
fiducial.lpa =  [41 108 84];
fiducial.rpa = [213 107 80];
fiducial.zpoint = [131 125 195];

% coarse realignment based on fiducials first
cfg          = [];
cfg.method   = 'fiducial';
cfg.fiducial = fiducial;
cfg.coordsys = 'neuromag';
mri          = ft_volumerealign(cfg,mri);

% make automatic coregistration based on headshape
% the coregistration looks fine (depending on how well you defined the
% fiducials of course) - don't change anything in the interactive menu to
% avoid that that's the problem
cfg           = [];
cfg.method    = 'headshape';
cfg.headshape.headshape = shape;
cfg.headshape.interactive    = 'no';
cfg.headshape.icp = 'yes';
cfg.coordsys  = 'neuromag';
mri_aligned   = ft_volumerealign(cfg, mri);

cfg.headshape.shape = ft_convert_units(shape, 'm');
mri_aligned2  = ft_volumerealign(cfg, mri);
mri_aligned3  = ft_volumerealign(cfg, ft_convert_units(mri, 'mm'));

%% check based on fiducial + headshape alignment
% segment the scalp to check the alignment
cfg        = [];
cfg.output = 'scalp';
cfg.smooth = 2;
seg_align  = ft_volumesegment(cfg, mri_aligned);
seg_align2  = ft_volumesegment(cfg, mri_aligned2);
seg_align3  = ft_volumesegment(cfg, mri_aligned3);

% build headmodel
cfg             = [];
cfg.method      = 'singleshell';
cfg.numvertices = 2000;
hdm    = ft_prepare_headmodel(cfg, seg_align);
hdm2   = ft_prepare_headmodel(cfg, seg_align2);
hdm3   = ft_prepare_headmodel(cfg, seg_align3);

% now plot
figure; hold on
ft_plot_headmodel(ft_convert_units(hdm, 'cm'),'edgecolor','none','facecolor','w');
ft_plot_headshape(shape);

figure; hold on
ft_plot_headmodel(ft_convert_units(hdm2, 'cm'),'edgecolor','none','facecolor','w');
ft_plot_headshape(shape);

figure; hold on
ft_plot_headmodel(ft_convert_units(hdm3, 'cm'),'edgecolor','none','facecolor','w');
ft_plot_headshape(shape);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = polynomialCurveSetFit(seg, varargin)
%POLYNOMIALCURVESETFIT Fit a set of polynomial curves to a segmented image
%
%   COEFS = polynomialCurveSetFit(IMG);
%   COEFS = polynomialCurveSetFit(IMG, DEG);
%   Result is a cell array of matrices. Each matrix is DEG+1-by-2, and
%   contains coefficients of polynomial curve for each coordinate.
%   IMG is first binarised, then skeletonized. Each cure
%
%   [COEFS LBL] = polynomialCurveSetFit(...);
%   also returns an image of labels for the segmented curves. The max label
%   is the number of curves, and the length of COEFS.
%
%   Requires the toolboxes:
%   - Optimization
%   - Image Processing
%
%   Example
%     % Fit a set of curves to a binary skeleton
%     img = imread('circles.png');
%     % compute skeleton, and ensure one-pixel thickness
%     skel = bwmorph(img, 'skel', 'Inf');
%     skel = bwmorph(skel, 'shrink');
%     figure; imshow(skel==0)
%     coeffs = polynomialCurveSetFit(skel, 2);
%     % Display segmented image with curves
%     figure; imshow(~img); hold on;
%     for i = 1:length(coeffs)
%         hc = drawPolynomialCurve([0 1], coeffs{i});
%         set(hc, 'linewidth', 2, 'color', 'g');
%     end
%
%   See also
%   polynomialCurves2d, polynomialCurveFit
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2007-03-21
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%% Initialisations

% default degree for curves
deg = 2;
if ~isempty(varargin)
    deg = varargin{1};
end

% ensure image is binary
seg = seg > 0;

%% Extract branching points and terminating points

% compute image of end points
imgEndPoints = imfilter(double(seg), ones([3 3])) .* seg == 2;

% compute centroids of end points
lblEndPoints    = bwlabel(imgEndPoints, 4);
regEndPoints    = bwconncomp(imgEndPoints, 4);
struct   = regionprops(regEndPoints, 'Centroid');

endPoints = cat(1, struct.Centroid);

% compute image of multiple points (intersections between curves)
imgBranching    = imfilter(double(seg), ones([3 3])) .* seg > 3;

% compute coordinate of nodes, as centroids of the multiple points
lblBranching = bwlabel(imgBranching, 4);
regBranching = bwconncomp(imgBranching, 4);
struct   = regionprops(regBranching, 'Centroid');

branchPoints = cat(1, struct.Centroid);

% list of nodes (all categories)
nodes = [branchPoints; endPoints];

% image of node labels
lblNodes = lblBranching;
lblNodes(lblEndPoints > 0) = lblEndPoints(lblEndPoints > 0) + size(branchPoints, 1);

% isolate branches
imgBranches = seg & ~imgBranching & ~imgEndPoints;
lblBranches = bwlabel(imgBranches, 8);

% number of curves
nBranches = max(lblBranches(:));

% allocate memory
coefs = cell(nBranches, 1);

% For each curve, find interpolated polynomial curve
for i = 1:nBranches
    %disp(i);

    % extract points corresponding to current curve
    imgBranch = lblBranches == i;
    points = chainPixels(imgBranch);

    % if number of points is not sufficient, simply create a line segment
    if size(points, 1) < max(deg+1-2, 2)
        % find labels of nodes
        inds = unique(lblNodes(imdilate(imgBranch, ones(3,3))));
        inds = inds(inds~=0);

        if length(inds)<2
            disp(['Could not find extremities of branch number ' num2str(i)]);
            coefs{i} = [0 0;0 0];
            continue;
        end

        % consider extremity nodes
        node0 = nodes(inds(1),:);
        node1 = nodes(inds(2),:);

        % use only a linear approximation
        xc = zeros(1, deg+1);
        yc = zeros(1, deg+1);
        xc(1) = node0(1);
        yc(1) = node0(2);
        xc(2) = node1(1)-node0(1);
        yc(2) = node1(2)-node0(2);

        % assigne au tableau de courbes
        coefs{i} = [xc;yc];

        % next branch
        continue;
    end

    % find nodes closest to first and last points of the current curve
    [dist, ind0] = minDistancePoints(points(1, :), nodes); %#ok<*ASGLU>
    [dist, ind1] = minDistancePoints(points(end, :), nodes);

    % add nodes to the curve.
    points = [nodes(ind0,:); points; nodes(ind1,:)]; %#ok<AGROW>

    % parametrization of the polyline
    t = parametrize(points);
    t = t/max(t);

    % fit a polynomial curve to the set of points
    [xc, yc] = polynomialCurveFit(...
        t, points, deg, ...
        0, {points(1,1), points(1,2)},...
        1, {points(end,1), points(end,2)});

    % stores result
    coefs{i} = [xc ; yc];
end

%% Post-processing

% manage outputs
if nargout == 1
    varargout = {coefs};
elseif nargout == 2
    varargout = {coefs, lblBranches};
end

function points = chainPixels(img, varargin)
%CHAINPIXELS return the list of points which constitute a curve on image
%   output = chainPixels(input)
%
%   Example
%   chainPixels
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2007-03-21
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

conn = 8;
if ~isempty(varargin)
    conn = varargin{1};
end

% matrice de voisinage
if conn==4
    f = [0 1 0;1 1 1;0 1 0];
elseif conn==8
    f = ones([3 3]);
end

% find extremity points
nb = imfilter(double(img), f).*img;
imgEnding = nb==2 | nb==1;
[yi, xi] = find(imgEnding);

% extract coordinates of points
[y, x] = find(img);

% index of first point
if isempty(xi)
    % take arbitrary point
    ind = 1;
else
    ind = find(x==xi(1) & y==yi(1));
end

% allocate memory
points  = zeros(length(x), 2);

if conn==8
    for i=1:size(points, 1)
        % avoid multiple neighbors (can happen in loops)
        ind = ind(1);

        % add current point to chained curve
        points(i,:) = [x(ind) y(ind)];

        % remove processed coordinate
        x(ind) = [];    y(ind) = [];

        % find next candidate
        ind = find(abs(x-points(i,1))<=1 & abs(y-points(i,2))<=1);
    end
else
    for i=1:size(points, 1)
        % avoid multiple neighbors (can happen in loops)
        ind = ind(1);

        % add current point to chained curve
        points(i,:) = [x(ind) y(ind)];

        % remove processed coordinate
        x(ind) = [];    y(ind) = [];

        % find next candidate
        ind = find(abs(x-points(i,1)) + abs(y-points(i,2)) <=1 );
    end
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
% RES = RCONV2(MTX1, MTX2, CTR)
%
% Convolution of two matrices, with boundaries handled via reflection
% about the edge pixels.  Result will be of size of LARGER matrix.
%
% The origin of the smaller matrix is assumed to be its center.
% For even dimensions, the origin is determined by the CTR (optional)
% argument:
%      CTR   origin
%       0     DIM/2      (default)
%       1     (DIM/2)+1

% Eero Simoncelli, 6/96.

function c = rconv2(a,b,ctr)

if (exist('ctr') ~= 1)
  ctr = 0;
end

if (( size(a,1) >= size(b,1) ) & ( size(a,2) >= size(b,2) ))
    large = a; small = b;
elseif  (( size(a,1) <= size(b,1) ) & ( size(a,2) <= size(b,2) ))
    large = b; small = a;
else
  error('one arg must be larger than the other in both dimensions!');
end

ly = size(large,1);
lx = size(large,2);
sy = size(small,1);
sx = size(small,2);

%% These values are one less than the index of the small mtx that falls on
%% the border pixel of the large matrix when computing the first
%% convolution response sample:
sy2 = floor((sy+ctr-1)/2);
sx2 = floor((sx+ctr-1)/2);

% pad with reflected copies
clarge = [
    large(sy-sy2:-1:2,sx-sx2:-1:2), large(sy-sy2:-1:2,:), ...
	large(sy-sy2:-1:2,lx-1:-1:lx-sx2); ...
    large(:,sx-sx2:-1:2),    large,   large(:,lx-1:-1:lx-sx2); ...
    large(ly-1:-1:ly-sy2,sx-sx2:-1:2), ...
      large(ly-1:-1:ly-sy2,:), ...
      large(ly-1:-1:ly-sy2,lx-1:-1:lx-sx2) ];

c = conv2(clarge,small,'valid');

repo: SelfExSR


text:
function test_graphWin()
%Validate opening and closing of graph windows
%
%   test_graphWin()
%
% Tests: newGraphWin(), closeGraphWin()
%
% INPUTS
%  No inputs
%
% RETURNS
%  No returns
%
% Example:  test_graphWin()
%
% See also MRVTEST
%
% Copyright Stanford team, mrVista, 2015

% Because the figure handles changed from numbers (Matlab2014a and prior)
% to a FIGURE class (Matlab 2014b and later), it is important to ensure
% that opening and closing of figures works properly.

% open a figure using the newGraphWin command
newGraphWin();

% close the figure using closeGraphWin
closeGraphWin;

return


text:
function prediction = AlexNet(varargin)
%AlexNet Returns an AlexNet model for ImageNet
%   M = models.AlexNet() returns the model proposed in:
%
%     Krizhevsky et al., "ImageNet classification with deep convolutional
%     neural networks", NIPS 2012.
%
%   models.AlexNet(..., 'option', value, ...) accepts the following
%   options:
%
%   `pretrained`:: false
%     If true, returns a model pre-trained on ImageNet (using the
%     MatConvNet example code).
%
%   `input`:: default input
%     Specifies an input (images) layer for the network. If unspecified, a
%     new one is created.
%
%   `numClasses`:: 1000
%     Number of output classes.
%
%   `batchNorm`:: true
%     Whether to use batch normalization.
%
%   `normalization`:: [5 1 0.0001/5 0.75]
%     Parameters for vl_nnnormalize layer (only used without batch-norm).
%
%   Any other options will be passed to models.ConvBlock(), and can be used
%   to change the activation function, weight initialization, etc.
%
%   Suggested SGD training options are also returned in the struct M.meta.

% Copyright (C) 2018 Joao F. Henriques, Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  % parse options. unknown arguments will be passed to ConvBlock (e.g.
  % activation).
  opts.pretrained = false ;  % whether to fetch a pre-trained model
  opts.input = Input('name', 'images', 'gpu', true) ;  % default input layer
  opts.numClasses = 1000 ;  % number of predicted classes
  opts.batchNorm = true ;  % whether to use batch normalization
  opts.normalization = [5 1 0.0001/5 0.75] ;  % for LRN layer (vl_nnnormalize)
  [opts, convBlockArgs] = vl_argparse(opts, varargin, 'nonrecursive') ;

  % default training options for this network (returned as output.meta)
  meta.batchSize = 256 ;
  meta.imageSize = [227, 227, 3] ;
  meta.augmentation.crop = 227 / 256;
  meta.augmentation.location = true ;
  meta.augmentation.flip = true ;
  meta.augmentation.brightness = 0.1 ;
  meta.augmentation.aspect = [2/3, 3/2] ;
  meta.weightDecay = 0.0005 ;

  % the default learning rate schedule
  if ~opts.pretrained
    if ~opts.batchNorm
      meta.learningRate = logspace(-2, -4, 60) ;
    else
      meta.learningRate = logspace(-1, -4, 20) ;
    end
    meta.numEpochs = numel(meta.learningRate) ;
  else  % fine-tuning has lower LR
    meta.learningRate = 1e-5 ;
    meta.numEpochs = 20 ;
  end

  % return a pre-trained model
  if opts.pretrained
    if opts.batchNorm
      warning('The pre-trained model does not include batch-norm (set batchNorm to false).') ;
    end
    if opts.numClasses ~= 1000
      warning('Model options are ignored when loading a pre-trained model.') ;
    end
    prediction = models.pretrained('imagenet-matconvnet-alex') ;

    % return prediction layer (not softmax)
    assert(isequal(prediction{1}.func, @vl_nnsoftmax)) ;
    prediction = prediction{1}.inputs{1} ;

    % replace input layer with the given one
    input = prediction.find('Input', 1) ;
    prediction.replace(input, opts.input) ;

    prediction.meta = meta ;
    return
  end

  % get conv block generator with the given options. default activation is
  % ReLU, with pre-activation batch normalization (can be overriden).
  conv = models.ConvBlock('batchNorm', opts.batchNorm, ...
    'preActivationBatchNorm', true, 'weightScale', 0.01, convBlockArgs{:}) ;

  % build network
  images = opts.input ;

  % first conv block
  x = conv(images, 'size', [11, 11, 3, 96], 'stride', 4) ;
  if ~opts.batchNorm
    x = vl_nnnormalize(x, opts.normalization) ;
  end
  x = vl_nnpool(x, 3, 'stride', 2) ;

  % second conv block
  x = conv(x, 'size', [5, 5, 48, 256], 'pad', 2) ;
  if ~opts.batchNorm
    x = vl_nnnormalize(x, opts.normalization) ;
  end
  x = vl_nnpool(x, 3, 'stride', 2) ;

  % conv blocks 3-5
  x = conv(x, 'size', [3, 3, 256, 384], 'pad', 1) ;
  x = conv(x, 'size', [3, 3, 192, 384], 'pad', 1) ;
  x = conv(x, 'size', [3, 3, 192, 256], 'pad', 1) ;
  x = vl_nnpool(x, 3, 'stride', 2) ;

  % first fully-connected block
  x = conv(x, 'size', [6, 6, 256, 4096]) ;
  if ~opts.batchNorm
    x = vl_nndropout(x) ;
  end

  % second fully-connected block
  x = conv(x, 'size', [1, 1, 4096, 4096]) ;
  if ~opts.batchNorm
    x = vl_nndropout(x) ;
  end

  % prediction layer
  prediction = conv(x, 'size', [1, 1, 4096, opts.numClasses], ...
    'batchNorm', false, 'activation', 'none') ;

  prediction.meta = meta ;

end

repo: MFT


text:
% Local Regression and Likelihood, Figure 6.1.
%
% Derivative (local slope) estimation, for the Old Faithful Geyser Data.
% The 'deriv' argument specifies derivative estimation,
%   'deriv',1        First-order derivative.
%   'deriv',[1 1]    Second-order derivative.
%   'deriv',2        For bivariate fits, partial deriv. wrt second variable.
%   'deriv',[1 2]    Mixed second-order derivative.
%
% Density estimation is done on the log-scale. That is, the estimate
% is of g(x) = log(f(x)), where f(x) is the density.
%
% The relation between derivatives is therefore
%     f'(x) = f(x)g'(x) = g'(x)exp(g(x)).
% To estimate f'(x), we must estimate g(x) and g'(x) (fit1 and fit2 below),
% evaluate on a grid of points (p1 and p2), and apply the back-transformation.
%
% Disclaimer: I don't consider derivative estimation from noisy data
% to be a well-defined problem. Use at your own risk.
%
% Author: Catherine Loader
%
% NEED: m argument passed to lfmarg().

load geyser;
fit1 = locfit(geyser,'alpha',[0.1 0.6],'ll',1,'ur',6);
fit2 = locfit(geyser,'alpha',[0.1 0.6],'ll',1,'ur',6,'deriv',1);
z = lfmarg(fit1);
p1 = predict(fit1,z);
p2 = predict(fit2,z);
figure('Name','fig6_1: slope estimation: Old faithful data' );
plot(z{1},p2.*exp(p1));
xlabel('Eruption Duration (Minutes)');
ylabel('Density Derivative');

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% OtSkimCubeDelete Delete a matrix from a skim cube.
%
% Usage:
%    OtSkimCubeDelete(hscube, pmturi)
% where:
%    hscube skim cube handle as created by OtMatrixCubeNew()
%    pmturi is a PMTURI struct as created by OtPmturi()
%

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 20006-04-24 Created (Johan Meijdam, VORtech)

function iret = OtSkimCubeDelete(hscube, pmturi)
	iret = Ot('skimcube.delete', hscube, pmturi);
end


text:
function R = quantize_ucm_or(ucm, nori, angSpan)
if nargin<2, nori = 8; end
if nargin<3, angSpan = 1; end % in [1,...,8]. angSpan = 1: no overlap. angSpan = 3: overlap of nori/pi. angSpan = 8: todos los canales tienen full ucm.

strength = ucm.strength(3:2:end,3:2:end);
[tx, ty] = size(strength);
R = zeros(tx, ty, nori);

for o = 0 : nori-1,

    angMin = (2*o-angSpan)/nori/2*pi;
    angMax = (2*o+angSpan)/nori/2*pi;

    if (angMin >= 0) && (angMax <= pi),
        bw = (ucm.orient > angMin) & (ucm.orient <= angMax);
    elseif (angMin < 0) && (angMax <= pi)
        bw = (ucm.orient > (pi+angMin)) | (ucm.orient <= angMax);
    elseif (angMin >= 0) && (angMax > pi)
        bw = (ucm.orient > angMin) | (ucm.orient <= (angMax-pi) );
    else
        bw = (ucm.orient > (pi+angMin)) | (ucm.orient <= (angMax-pi) );
    end

    R(:, :, o+1) = strength .* ( bw & (ucm.orient~=0));
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function msh = meshApplySettings(msh, settings);
% Apply stored settings to a mrMesh window.
%
% msh = meshApplySettings(<msh>, <settings=get from UI>);
%
% See meshSettings for info on the settings struct. This contains info
% about the lights, camera, main mesh rotation, cursor, and ROI locations.
%
% If the first 'msh' argument is omitted, gets the selected mesh from the
% selected gray view. If settings is omitted, assumes this function was
% called from a GUI, and searches for a listbox uicontrol with the tag
% 'MeshSettingsList'.
%
% ras, 03/14/06.
if notDefined('msh'), msh = viewGet(getSelectedGray, 'selectedMesh');  end

if notDefined('settings')
    % get from UI -- first need to implement the UI :)
    return
end

verbose = prefsVerboseCheck;

if verbose
	hmsg = msgbox('Applying Stored Mesh Settings...');
end

% set background color
msh = mrmSet(msh, 'background', settings.backgroundColor);

% set camera settings
p = settings.camera;
p.actor = 0;
mrMesh(msh.host, msh.id, 'set', p);

% % get an actor list for the lights, mesh below
% actorList = mrmGet(msh, 'ListOfActors');
% actors = actorList.objectList;
%  (temp disabled, very slow, faster if we assume lights are 33 and 34)

% check there are enough actors for the settings (maybe fewer lights, etc)
% N.Y.I.

% set lights
for i = 1:2 % min(length(settings.lights), length(actors)-1)
    if ~isempty(settings.lights(i))
        p = settings.lights(i);
        p.actor = i+32; % actors(i+1);
		try
			% get lighting values from the settings -- sometimes this
			% hangs:
			p.ambient = settings.lights(i).ambient;
			p.diffuse = settings.lights(i).diffuse;
		catch
			% use default values
			p.ambient = [.4 .4 .3];
			p.diffuse = [.5 .5 .6];
		end
        mrMesh(msh.host, msh.id, 'set', p);
    end
end

% set cursor position
if settings.cursor > 0
    try
        mrmSet(msh, 'CursorVertex', settings.cursor);
    catch
        disp('Could not set cursor')
    end
end

% set main mesh rotation, origin
p = mrmGet(msh, 'ActorData', 32);
p.rotation = settings.mesh.rotation;
p.origin = settings.mesh.origin;
p.actor = 32;
mrMesh(msh.host, msh.id, 'set', p);

% % set ROI (n.y.i.)
% if ~isfield(settings.roi, 'error')
%     try
%
%     catch
%         disp('Could not set ROI')
%     end
% end

if verbose, close(hmsg); end

return


text:
function y = rms(f,varargin)
%RMS RMS value of signal
%   Usage: y = rms(f);
%          y = rms(f,...);
%
%   `RMS(f)` computes the RMS (Root Mean Square) value of a finite sampled
%   signal sampled at a uniform sampling rate. This is a vector norm
%   equal to the $l^2$ averaged by the length of the signal.
%
%   If the input is a matrix or ND-array, the RMS is computed along the
%   first (non-singleton) dimension, and a vector of values is returned.
%
%   The RMS value of a signal *x* of length *N* is computed by
%
%   ..                       N
%      rms(f) = 1/sqrt(N) ( sum |f(n)|^2 )^(1/2)
%                           n=1
%
%   .. math:: rms(f) = \frac{1}{\sqrt N} \left( \sum_{n=1}^N |f(n)|^2
%      \right)^{\frac{1}{2}}
%
%   `RMS` takes the following flags at the end of the line of input
%   parameters:
%
%     'ac'       Consider only the AC component of the signal (i.e. the mean is
%                removed).
%
%     'dim',d    Work along specified dimension. The default value of `[]`
%                means to work along the first non-singleton one.
%

%   AUTHOR : Peter L. Søndergaard

%% ------ Checking of input parameters ---------

if ~isnumeric(f)
  error('%s: Input must be numerical.',upper(mfilename));
end;

if nargin<1
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.keyvals.dim=[];
definput.flags.mean={'noac','ac'};
[flags,kv]=ltfatarghelper({'dim'},definput,varargin);

%% ------ Computation --------------------------

% It is better to use 'norm' instead of explicitly summing the squares, as
% norm (hopefully) attempts to avoid numerical overflow.

[f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,[],kv.dim, ...
                                                  upper(mfilename));
permutedsize(1)=1;
y=zeros(permutedsize);
if flags.do_ac

  for ii=1:W
    y(1,ii) = norm(f(:,ii)-mean(f(:,ii)))/sqrt(L);
   end;

else

  for ii=1:W
    y(1,ii)=norm(f(:,ii))/sqrt(L);
  end;

end;

y=assert_sigreshape_post(y,kv.dim,permutedsize,order);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [prediction rfParams] = rmPlotGUI_scaleFactor(M, pred, rfParams, recompFit);
%
%
%
%
% ras, 11/2008.
model     = M.model{M.modelNum};
modelName = rmGet(model, 'desc');
voxel     = get(M.ui.voxel.sliderHandle, 'Value');
coords	  = M.coords(:,voxel);
if isequal(M.viewType, 'Gray')  % convert coords into an index
    coords = M.coords(voxel);
end

%% get/make trends
[trends, ntrends, dcid] = rmMakeTrends(M.params, 0);
if isfield(M.params.analysis,'allnuisance')
    trends = [trends M.params.analysis.allnuisance];
end

%% Compute prediction; compute beta scale factors according to the model
% we also add this to the rfParams, to report later
switch modelName,
    case {'2D pRF fit (x,y,sigma, positive only)',...
            '2D RF (x,y,sigma) fit (positive only)'};
        if recompFit==0,
            beta = rmCoordsGet(M.viewType, model, 'b', coords);
            beta = beta([1 dcid+1]);
        else
            beta = pinv([pred trends(:,dcid)])*M.tSeries(:,voxel);
            beta(1) = max(beta(1),0);

            % also recompute variance explained
            prediction = [pred trends(:,dcid)] * beta;
            prediction = prediction(:,1);
        end

        prediction = [pred trends(:,dcid)] * beta;
        RFs        = RFs .* (beta(1) .* M.params.analysis.HrfMaxResponse);

        rfParams(4) = beta(1);

    case {'2D pRF fit (x,y,sigma_major,sigma_minor)' ...
			'oval 2D pRF fit (x,y,sigma_major,sigma_minor,theta)'};
        if recompFit==0,
            beta = rmCoordsGet(M.viewType, model, 'b', coords);
            beta = beta([1 dcid+1]);
        else
            beta = pinv([pred trends(:,dcid)])*M.tSeries(:,voxel);
            beta(1) = max(beta(1),0);

            % also recomput variance explained
            prediction = [pred trends(:,dcid)] * beta;
            prediction = prediction(:,1);
            rss = sum((M.tSeries(:,voxel)-prediction).^2);
            rawrss = sum(M.tSeries(:,voxel).^2);
            varexp = 1 - rss./rawrss;
        end

        prediction = [pred trends(:,dcid)] * beta;
        RFs        = RFs .* (beta(1) .* M.params.analysis.HrfMaxResponse);

        rfParams(4) = beta(1);

    case 'unsigned 2D pRF fit (x,y,sigma)';
        if recompFit==0,
            beta = rmCoordsGet(M.viewType, model, 'b', coords);
            beta = beta([1 dcid+1]);
        else
            beta = pinv([pred trends(:,dcid)])*M.tSeries(:,voxel);

            % also recomput variance explained
            prediction = [pred trends(:,dcid)] * beta;
            prediction = prediction(:,1);
        end

        prediction = [pred trends(:,dcid)] * beta;
        RFs        = RFs .* (beta(1) .* M.params.analysis.HrfMaxResponse);

        rfParams(4) = beta(1);

   case {'Double 2D pRF fit (x,y,sigma,sigma2, center=positive)'},
        if recompFit==0,
            beta = rmCoordsGet(M.viewType, model, 'b', coords);
            beta = beta([1 2 dcid+2]);
        else
            beta = pinv([pred trends(:,dcid)])*M.tSeries(:,voxel);
            beta(1) = max(beta(1),0);
            beta(2) = max(beta(2),-abs(beta(1)));

            % also recomput variance explained
            prediction = [pred trends(:,dcid)] * beta;
            prediction = prediction(:,1);
        end

        prediction = [pred trends(:,dcid)] * beta;
        RFs        = RFs * (beta(1:2).*M.params.analysis.HrfMaxResponse);

        rfParams(:,4) = beta(1);

    case {'Two independent 2D pRF fit (2*(x,y,sigma, positive only))'},
        if recompFit==0,
            beta = rmCoordsGet(M.viewType, model, 'b', coords);
            beta = beta([1 2 dcid+2]);
        else
            beta = pinv([pred trends(:,dcid)])*M.tSeries(:,voxel);
            beta(1:2) = max(beta(1:2),0);

            prediction = [pred trends(:,dcid)] * beta;
            prediction = prediction(:,1);
        end

        prediction = [pred trends(:,dcid)] * beta;
        RFs        = RFs * (beta(1:2) .* M.params.analysis.HrfMaxResponse);

        rfParams(:,4) = beta(1:2);
        rfParams = rfParams(1,:);

  case {'Sequential 2D pRF fit (2*(x,y,sigma, positive only))'},
        if recompFit==0,
            beta = rmCoordsGet(M.viewType, model, 'b', coords);
            beta = beta([1 2 dcid+2]);
        else
            beta = pinv([pred trends(:,dcid)])*M.tSeries(:,voxel);
            beta(1:2) = max(beta(1:2),0);
            % also recomput variance explained
            prediction = [pred trends(:,dcid)] * beta;
            prediction = prediction(:,1);
        end

        prediction = [pred trends(:,dcid)] * beta;
        RFs        = RFs * (beta(1:2) .* M.params.analysis.HrfMaxResponse);

        rfParams(:,4) = beta(1:2);
        rfParams = rfParams(1,:);

   otherwise,
        error('Unknown modelName: %s', modelName);
end;

return


text:
function Y = sptenrand(sz,nz)
%SPTENRAND Sparse uniformly distributed random tensor.
%
%   R = SPTENRAND(sz,density) creates a random sparse tensor of the
%   specified sz with approximately density*prod(sz) nonzero
%   entries.
%
%   R = SPTENRAND(sz,nz) creates a random sparse tensor of the
%   specified sz with approximately nz nonzero entries.
%
%   Example: R = sptenrand([5 4 2],12);
%
%   See also SPTENSOR, TENRAND, RAND.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

% Error check on siz
if ndims(sz) ~= 2 || size(sz,1) ~= 1
    error('Size must be a row vector');
end

% Error check on nz
if ~exist('nz','var') || (nz < 0)
    error('2nd argument must be positive');
end

% Is nz an number or a fraction? Ultimately want a number.
if (nz < 1)
    nz = prod(sz) * nz;
end

% Make sure nz is an integer
nz = ceil(nz);

% Keep iterating until we find enough unique nonzeros or we
% give up
subs = [];
cnt = 0;
while (size(subs,1) < nz) && (cnt < 10)
    subs = ceil( rand(nz, size(sz,2)) * diag(sz) );
    subs = unique(subs, 'rows');
    cnt = cnt + 1;
end

% Extract nnz subscipts and create a corresponding list of
% values
nz = min(nz, size(subs,1));
subs = subs(1:nz,:);
vals = rand(nz,1);

Y = sptensor(subs,vals,sz);
return;

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
% eeg_rejsuperpose() - superpose rejections of a EEG dataset.
%
% Usage:
%   >> EEGOUT = eeg_rejsuperpose( EEGIN, typerej, Rmanual, Rthres, ...
%                    Rconst, Rent, Rkurt, Rfreq, Rothertype);
%
% Inputs:
%   EEGIN      - input dataset
%   typerej    - type of rejection (1=raw data; 0=ica).
%   Rmanual    - include manual rejection (0|1).
%   Rthres     - include threshold rejection (0|1).
%   Rconst     - include rejection of constant activity (0|1).
%   Rent       - include entropy rejection (0|1).
%   Rkurt      - include kurtosis rejection (0|1).
%   Rfreq      - include frequcy based rejection (0|1).
%   Rothertype - include manual rejection (0|1).
%
% Outputs:
%   EEGOUT     - with rejglobal and rejglobalE fields updated
%
% See also: eeglab()

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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

function [EEG, com] = eeg_rejsuperpose( EEG, typerej, Rmanual, Rthres, Rconst, ...
                              Rent, Rkurt, Rfreq, Rothertype);

if nargin < 9
	help eeg_rejsuperpose;
	return;
end;

typerej = ~typerej;
rejglobal  = zeros( 1, EEG.trials);
if typerej == 0
    rejglobalE = zeros( EEG.nbchan, EEG.trials);
else
    rejglobalE = zeros( size(EEG.icaweights,1), EEG.trials);
end

if typerej == 0 | Rothertype
	if Rmanual
		rejglobal  = rejarray( rejglobal,  EEG.reject.rejmanual); % see bottom for the
		rejglobalE = rejarray( rejglobalE, EEG.reject.rejmanualE); % function rejarray
	end;
	if Rthres
		rejglobal  = rejarray( rejglobal,  EEG.reject.rejthresh);
		rejglobalE = rejarray( rejglobalE, EEG.reject.rejthreshE);
	end;
	if Rfreq
		rejglobal  = rejarray( rejglobal,  EEG.reject.rejfreq);
		rejglobalE = rejarray( rejglobalE, EEG.reject.rejfreqE);
	end;
	if Rconst
		rejglobal  = rejarray( rejglobal,  EEG.reject.rejconst);
		rejglobalE = rejarray( rejglobalE, EEG.reject.rejconstE);
	end;
	if Rent
		rejglobal  = rejarray( rejglobal,  EEG.reject.rejjp);
		rejglobalE = rejarray( rejglobalE, EEG.reject.rejjpE);
	end;
	if Rkurt
		rejglobal  = rejarray( rejglobal,  EEG.reject.rejkurt);
		rejglobalE = rejarray( rejglobalE, EEG.reject.rejkurtE);
	end;
end;

% ---------------

if typerej == 1 | Rothertype
	if Rmanual
		rejglobal  = rejarray( rejglobal,  EEG.reject.icarejmanual);
		rejglobalE = rejarray( rejglobalE, EEG.reject.icarejmanualE);
	end;
	if Rthres
		rejglobal  = rejarray( rejglobal,  EEG.reject.icarejthresh);
		rejglobalE = rejarray( rejglobalE, EEG.reject.icarejthreshE);
	end;
	if Rfreq
		rejglobal  = rejarray( rejglobal,  EEG.reject.icarejfreq);
		rejglobalE = rejarray( rejglobalE, EEG.reject.icarejfreqE);
	end;
	if Rconst
		rejglobal  = rejarray( rejglobal,  EEG.reject.icarejconst);
		rejglobalE = rejarray( rejglobalE, EEG.reject.icarejconstE);
	end;
	if Rent
		rejglobal  = rejarray( rejglobal,  EEG.reject.icarejjp);
		rejglobalE = rejarray( rejglobalE, EEG.reject.icarejjpE);
	end;
	if Rkurt
		rejglobal  = rejarray( rejglobal,  EEG.reject.icarejkurt);
		rejglobalE = rejarray( rejglobalE, EEG.reject.icarejkurtE);
	end;
end;

EEG.reject.rejglobal = rejglobal;
EEG.reject.rejglobalE = rejglobalE;

com =sprintf('%s = eeg_rejsuperpose( %s, %d, %d, %d, %d, %d, %d, %d, %d);', ...
	  inputname(1), inputname(1), ~typerej, Rmanual, Rthres, Rconst, Rent, Rkurt, Rfreq, Rothertype);

return;

% subfunction rejecting an array ------
function dest = rejarray( dest, ori)
    if isempty(dest)
        dest = ori;
    elseif ~isempty(ori)
		dest = dest | ori;
	end;
return;


text:
function[res] = varargToStr(separator, varargin)
% [res] = varargToStr([separator], varargin)
%
% Convert a cell of sth. to a single str.
% Output is similar to evalc('disp(proposalsVars)'), but without spaces and
% fancy characters.
%
% Copyright by Holger Caesar, 2015

res = '';
if ~exist('separator', 'var') || isempty(separator),
    separator = '';
end;

for i = 1 : numel(varargin),
    cur = varargin{i};

    if ischar(cur),
        % Do nothing
    elseif isnumeric(cur),
        cur = num2str(cur);
    elseif iscell(cur),
        cur = varargToStr(separator, cur{:});
    elseif isfunc(cur),
        cur = func2str(cur);
    else
        error('Error: Invalid type!');
    end;

    % Concatenate result
    if mod(i, 2) == 1 && i > 1,
        res = [res, separator, cur];
    else
        res = [res, cur];
    end;
end;

author: nightrome

repo: matconvnet-calvin


text:
function value = p34_autonomous ( )

%*****************************************************************************80
%
%% P34_AUTONOMOUS reports whether problem P34 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%*******************************************************************
%% HSDNTrhsfun: compute the right-hand side vector of the
%%              Schur complement equation for the NT direction.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

function [rhs,EinvRc,hRd] = HSDNTrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ)
global spdensity

m = par.m;
if (nargin > 8)
    corrector = 1;
else
    corrector = 0;
    hRd = zeros(m+2,1);
end
hEinvRc = zeros(m+2,1);
EinvRc  = cell(size(blk,1),1);
if length(sigmu)==1; sigmu = sigmu*ones(1,size(blk,1)); end
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    n = sum(pblk{2});  numblk = length(pblk{2});
    if strcmp(pblk{1},'l')
        if (corrector)
            Rq = dX{p}.*dZ{p};
        else
            Rq = sparse(n,1);
            tmp  = par.dd{p}.*Rd{p};
            tmp2 = mexMatvec(At{p},tmp,1);
            hRd = hRd + tmp2;
        end
        EinvRc{p} = sigmu(p)./Z{p}-X{p} -Rq./Z{p};
        tmp2 = mexMatvec(At{p},EinvRc{p},1);
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'q')
        w = sqrt(par.gamz{p}./par.gamx{p});
        if (corrector)
            hdx = qops(pblk,w,par.ff{p},5,dX{p});
            hdz = qops(pblk,w,par.ff{p},6,dZ{p});
            hdxdz = Arrow(pblk,hdx,hdz);
            vv = qops(pblk,w,par.ff{p},5,X{p});
            Vihdxdz = Arrow(pblk,vv,hdxdz,1);
            Rq = qops(pblk,w,par.ff{p},6,Vihdxdz);
        else
            Rq = sparse(n,1);
            tmp  = par.dd{p}.*Rd{p} + qops(pblk,qops(pblk,Rd{p},par.ee{p},1),par.ee{p},3);
            tmp2 = mexMatvec(At{p},tmp,1);
            hRd = hRd + tmp2;
        end
        EinvRc{p} = qops(pblk,-sigmu(p)./(par.gamz{p}.*par.gamz{p}),Z{p},4)-X{p}-Rq;
        tmp2 = mexMatvec(At{p},EinvRc{p},1);
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'s')
        n2 = pblk{2}.*(pblk{2}+1)/2;
        if (corrector)
            hdZ = Prod3(pblk,par.G{p},dZ{p},par.G{p}',1);
            hdX = spdiags(-par.sv{p},0,n,n)-hdZ;
            tmp = Prod2(pblk,hdX,hdZ,0);
            tmp = 0.5*(tmp+tmp');
            if (numblk == 1)
                d = par.sv{p};
                e = ones(pblk{2},1);
                Rq = 2*tmp./(d*e'+e*d');
                if (nnz(Rq) <= spdensity*n2); Rq = sparse(Rq); end
            else
                Rq = sparse(n,n);
                s = [0, cumsum(pblk{2})];
                for i = 1:numblk
                    pos = s(i)+1 : s(i+1);
                    d = par.sv{p}(pos); e = ones(length(pos),1);
                    Rq(pos,pos) = 2*tmp(pos,pos)./(d*e' + e*d'); %#ok
                end
            end
        else
            Rq = sparse(n,n);
            EinvRc{p} = Prod3(pblk,par.W{p},Rd{p},par.W{p},1,par.nzlistAy{p});
            tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),EinvRc(p));
            hRd = hRd + tmp2;
        end
        tmp = spdiags(sigmu(p)./par.sv{p} -par.sv{p},0,n,n);
        EinvRc{p} = Prod3(pblk,par.G{p}',tmp-Rq,par.G{p},1);
        tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),EinvRc(p));
        hEinvRc = hEinvRc + tmp2;
    end
end
%%
rhs = rp + hRd - hEinvRc;
rhs(m+1) = rhs(m+1) + (par.mu/par.tau - par.kap);
if (corrector)
    rhs(m+1) = rhs(m+1) - par.dtau*par.dkap/par.tau;
end
%%*******************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP