text:
function varargout = cylinderMesh(cyl, varargin)
%CYLINDERMESH  Create a 3D mesh representing a cylinder
%
%   [V F] = cylinderMesh(CYL)
%   Computes vertex coordinates and face vertex indices of a mesh
%   representing a 3D cylinder given as [X1 Y1 Z1 X2 Y2 Z2 R].
%
%   Example
%     % Draw a rotated cylinder
%     cyl = [0 0 0 10 20 30 5];
%     [v f] = cylinderMesh(cyl);
%     figure;drawMesh(v, f, 'FaceColor', 'r');
%     view(3); axis equal;
%
%     % Draw three mutually intersecting cylinders
%       p0 = [30 30 30];
%       p1 = [90 30 30];
%       p2 = [30 90 30];
%       p3 = [30 30 90];
%       [v1 f1] = cylinderMesh([p0 p1 25]);
%       [v2 f2] = cylinderMesh([p0 p2 25]);
%       [v3 f3] = cylinderMesh([p0 p3 25]);
%       figure; hold on;
%       drawMesh(v1, f1, 'FaceColor', 'r');
%       drawMesh(v2, f2, 'FaceColor', 'g');
%       drawMesh(v3, f3, 'FaceColor', 'b');
%       view(3); axis equal
%       set(gcf, 'renderer', 'opengl')
%
%   See also
%     drawCylinder, torusMesh, sphereMesh

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

% extract cylinder data
p1 = cyl(:, 1:3);
p2 = cyl(:, 4:6);
r  = cyl(:, 7);

% compute length and orientation
[theta, phi, rho] = cart2sph2d(p2 - p1);

% parametrisation on x
t = linspace(0, 2*pi, 20);
lx = r * cos(t);
ly = r * sin(t);

% parametrisation on z
lz = linspace(0, rho, 10);

% generate surface grids
x = repmat(lx, [length(lz) 1]);
y = repmat(ly, [length(lz) 1]);
z = repmat(lz', [1 length(t)]);

% transform points
trans   = localToGlobal3d(p1, theta, phi, 0);
[x, y, z] = transformPoint3d(x, y, z, trans);

% convert to FV mesh
[vertices, faces] = surfToMesh(x, y, z, 'xPeriodic', true);

% format output
varargout = formatMeshOutput(nargout, vertices, faces);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function sonorobot(r,ns,ds,mapoff)
if nargin<4
mapoff=false;
if nargin<3
ds=.1;
if nargin<2
ns=24;
if nargin<1;
r=.37;
end
end
end
end
clf reset
shg
set(gcf,'menubar','none',...%[0    0.2850    0.3985]
    'numbertitle','off','name','Sono Robot');
t_n=uicontrol('style','toggle','position',[35,65,25,25],...
    'string','N');
t_s=uicontrol('style','toggle','position',[35,5,25,25],...
    'string','S');
t_w=uicontrol('style','toggle','position',[5,35,25,25],...
    'string','W');
t_e=uicontrol('style','toggle','position',[65,35,25,25],...
    'string','E');
t_ne=uicontrol('style','toggle','position',[65,65,25,25],...
    'string','NE');
t_se=uicontrol('style','toggle','position',[65,5,25,25],...
    'string','SE');
t_sw=uicontrol('style','toggle','position',[5,5,25,25],...
    'string','SW');
t_nw=uicontrol('style','toggle','position',[5,65,25,25],...
    'string','NW');
t_start=uicontrol('style','toggle','position',[505,35,50,25],...
    'string','Auto');
t_quit=uicontrol('style','toggle','position',[505,5,50,25],...
    'string','Quit');
msg1=uicontrol('style','text','position',[120,5,345,55]);
msg2=uicontrol('style','text','position',[120,65,345,25],...
    'string','Measurement results...','fontsize',11.5,...
    'horizontal','left');
% orderdisp=['1st';'2nd';'3th';'4th';'5th';'6th';'7th';'8th';'9th';'0th'];
drawnow
hold on
axisx=12.5;
axisy=10;
axis([0 axisx 0 axisy]);
axis off

%Note:
%distance from one point to line:
%x=(x0*(x1-x2)^2+(x2*(-y0+y1)+x1*(y0-y2))*(y1-y2))/(x1^2-2*x1*x2+x2^2+(y1-y2)^2);
%y=(x2^2*y1+y0*y1^2+x0*(x1-x2)*(y1-y2)+x1^2*y2-2*y0*y1*y2+y0*y2^2-x1*x2*(y1+y2))/(x1^2-2*x1*x2+x2^2+(y1-y2)^2);

%Input environment & robot position
%Environment
nmax=50;
stpmax=100;
segset=zeros(nmax,4);%line segment set:[Ax,Ay,Bx,By]
areax=1;
areay=.5;
dftarea=[axisx-areax,axisy-areay;axisx-areax,axisy;axisx,axisy;axisx,axisy-areay];
okarea=[axisx-areax,axisy-2*areay;axisx-areax,axisy-areay;axisx,axisy-areay;axisx,axisy-2*areay];
fill(okarea(:,1),okarea(:,2),'k')
fill(dftarea(:,1),dftarea(:,2),'w')
stt=true;%true: on, ready to insert; false: off, waiting for the 2nd point
polygon=false;%%
n=0;%number of segments
%debug%
mdldisp=['A_1';'A_2';'A_3';'B_1';'B_2';'B_3';'B_4';'B_5';'B_6';'C  ';...
    'D_1';'D_2';'D_3';'D_4';'E  ';'nil';'DK '];
%--%
while n<=nmax
    %display current operation
    if n==0
        set(msg1,'string',['Insert the first edge. Or click the white '...
            'area and load default data.'],'fontsize',11.5,'horizontal','left');
    else
        if stt
            set(msg1,'string',['Insert a new line: The ',orderdisp(n+1),...
                ' edge. First end point. Or click the black area to finish.'],...
                'fontsize',11.5,'horizontal','left');
        else
            set(msg1,'string',['Insert a new line: The ',orderdisp(n),...
                ' edge. Second end point.'],'fontsize',11.5,...
                'horizontal','left');
        end
    end
    %ginput
    [x,y]=ginput(1);
    if (x-(axisx-areax))*(x-axisx)<0&&(y-(axisy-2*areay))*(y-(axisy-areay))<0||((x-(axisx-areax))*(x-axisx)<0&&(y-(axisy-areay))*(y-axisy)<0&&n==0)
        break
    end
    plot(x,y,'ro')
    if stt
        n=n+1;
        segset(n,1:2)=[x,y];
    else
        segset(n,3:4)=[x,y];
        plot(segset(n,[1,3]),segset(n,[2,4]),'b');
    end
    stt=~stt;
end
if n~=nmax
    segset=segset(1:n,:);
end

%default map
dftmap=[2,2,2,7;2,7,3,8;3,8,8,8;8,8,8,2;8,2,4,2;4,2,4,6;4,6,6,6;6,6,6,4];
if n==0
    segset=dftmap;
    n=8;
end

%Robot
% display current operation
set(msg1,'string','Set original position of the robot. Or click the white area and load default data.',...
    'fontsize',11.5,'horizontal','left');
[rx,ry]=ginput(1);
if (rx-(axisx-areax))*(rx-axisx)<0&&(ry-(axisy-areay))*(ry-axisy)<0
    rx=5;
    ry=5;
end
plot(rx,ry,'ro')
a=2*pi/ns;

%Probing & Moving
set(msg1,'string','Looking for exit...','fontsize',11.5,'horizontal','left');
rtrack=zeros(stpmax,2);rtrack(1,:)=[rx,ry];
start=false;
stp=0;
istouch=false;
while ~get(t_quit,'value')
    if ~get(t_start,'value')
        v_n=get(t_n,'value');
        v_s=get(t_s,'value');
        v_w=get(t_w,'value');
        v_e=get(t_e,'value');
        v_ne=get(t_ne,'value');
        v_se=get(t_se,'value');
        v_sw=get(t_sw,'value');
        v_nw=get(t_nw,'value');
        if any([v_n,v_s,v_w,v_e,v_ne,v_se,v_sw,v_nw])||~start
            start=true;
            dir=find([v_n,v_s,v_w,v_e,v_ne,v_se,v_sw,v_nw]);
            if dir
            %movement detected
            rxo=rx;
            ryo=ry;
            switch dir(1)
                case 1
                    ry=ry+ds;
                    set(t_n,'value',0)
                case 2
                    ry=ry-ds;
                    set(t_s,'value',0)
                case 3
                    rx=rx-ds;
                    set(t_w,'value',0)
                case 4
                    rx=rx+ds;
                    set(t_e,'value',0)
                case 5
                    rx=rx+ds;
                    ry=ry+ds;
                    set(t_ne,'value',0)
                case 6
                    rx=rx+ds;
                    ry=ry-ds;
                    set(t_se,'value',0)
                case 7
                    rx=rx-ds;
                    ry=ry-ds;
                    set(t_sw,'value',0)
                case 8
                    rx=rx-ds;
                    ry=ry+ds;
                    set(t_nw,'value',0)
            end
            if approachwall(segset,[rx,ry],2*r)
                set(msg1,'string','Approaching wall...',...
                    'fontsize',11.5,'horizontal','left');
            else
                set(msg1,'string','Searching exit...',...
                    'fontsize',11.5,'horizontal','left');
            end
            [istouch,~]=approachwall(segset,[rx,ry],r);
            if istouch
                set(msg1,'string','Robot has touched wall... Move back.',...
                    'fontsize',11.5,'horizontal','left');
                rx=rxo;
                ry=ryo;
            end
            end
            if ~istouch
                stp=stp+1;
                if stp>stpmax
                    rtrack=[rtrack;zeros(stpmax,2)];
                    stpmax=2*stpmax;
                end
%                 %debug%
%                 disp(stp)
%                 %--%
                rtrack(stp,:)=[rx,ry];
                %reprobe
                rs=sonomeasure(segset,[rx,ry],polygon,ns);
                [modelpt,mptadd,finished,mdl]=sonomodel(rs,r);%!!
                %redraw
                cla
                % environment
                if ~mapoff
                    for k=1:n
                        plot(segset(k,[1,3]),segset(k,[2,4]),'--r','LineWidth',2)
                    end
                end
                % probing data
                for k=1:ns
                    plot(rx+rs(k)*cos(((k-1):.1:k)*a),ry+rs(k)*sin(((k-1):.1:k)*a),'b:')
                end
                plot(rx,ry,'ro');%text(rx-.15,ry-.15,'Robot');
                plot(rx+r*cos((0:.1:2)*pi),ry+r*sin((0:.1:2)*pi),'b','LineWidth',1.45);
                % model report
                mdlrpt=[0 0];
                for k=1:ns
                    if finished(k)
                        plot(rx+modelpt(k,[1,3]),ry+modelpt(k,[2,4]),'k','LineWidth',1.5)
                        switch mdl(k)
                            case {1,2,3}
                                plot(rx+(rs(k)-.2)*cos((k-1:.1:k)*a),ry+(rs(k)-.2)*sin((k-1:.1:k)*a),'c','LineWidth',2.5);
                                mdlrpt(1)=1;
                            case {4,5,6,7,8,9}
                                plot(rx+(rs(k)-.2)*cos((k-1:.1:k)*a),ry+(rs(k)-.2)*sin((k-1:.1:k)*a),'g','LineWidth',2.5);
                                mdlrpt(2)=1;
                        end
                    end
                end
                for k=1:size(mptadd,1)
                    plot(rx+mptadd(k,[1,3]),ry+mptadd(k,[2,4]),'k','LineWidth',1.5)
                end
                if mdlrpt(1)&&mdlrpt(1)
                    set(msg2,'string','Faultages found...',...
                        'fontsize',11.5,'horizontal','left');
                elseif mdlrpt(1)
                    set(msg2,'string','Faultages on two sides...',...
                        'fontsize',11.5,'horizontal','left');
                elseif mdlrpt(2)
                    set(msg2,'string','Faultages on one side...',...
                        'fontsize',11.5,'horizontal','left');
                else
                    set(msg2,'string','Measurement results...',...
                        'fontsize',11.5,'horizontal','left');
                end
                % robot track
                plot(rtrack(1:stp,1),rtrack(1:stp,2),'b:')
                %debug%
                % report model
                for k=1:ns
                    text(rx+rs(k)*cos((k-.5)*a),ry+rs(k)*sin((k-.5)*a),mdldisp(mdl(k),:))
                end
            end
        end
    else
%         %test
%         pause(0.5)
%         beep
        %--%
        %Auto Move

    end
    pause(0.01);
end
close
end

author: Sable

repo: mcbench-benchmarks


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

%%
clear; close all; clc;

%% Syntax
% |FUNCTION_NAME|

%% Description
% FUNCTION_NAME

%% Examples

%%
% Plot settings

%% Example:

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%SVMR SVM regression
%
%      W = SVMR(X,NU,KTYPE,KPAR,EP)
%      W = X*SVMR([],NU,KTYPE,KPAR,EP)
%      W = X*SVMR(NU,KTYPE,KPAR,EP)
%
% INPUT
%   X      Regression dataset
%   NU     Fraction of objects outside the 'data tube'
%   KTYPE  Kernel type (default KTYPE='p', for polynomial)
%   KPAR   Extra parameter for the kernel
%   EP     Epsilon, with of the 'data tube'
%
% OUTPUT
%   W      Support vector regression
%
% DESCRIPTION
% Train an nu-Support Vector Regression on dataset X with parameter NU.
% The kernel is defined by kernel type KTYPE and kernel parameter KPAR.
% For the definitions of these kernels, have a look at proxm.m.
%
% SEE ALSO (PRTools Guide)
%  LINEARR, PROXM, GPR, SVC

% Copyright: D.M.J. Tax, D.M.J.Tax@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function y = svmr(varargin)

	argin = shiftargin(varargin,{'scalar','char'},1);
	argin = shiftargin(argin,'char',2);
  argin = setdefaults(argin,[],0.01,'p',1,0.2);

  if mapping_task(argin,'definition')
    y = define_mapping(argin,'untrained');
    y = setname(y,'SVM regression');

  elseif mapping_task(argin,'training')			% Train a mapping.

    [x,nu,ktype,kpar,ep] = deal(argin{:});
    [n,d] = size(x);
    y = gettargets(x);
    % kernel mapping:
    wk = proxm(+x,ktype,kpar);
    K = +(x*wk);
    % setup optimization:
    tol = 1e-6;
    C = 1/(n*nu);
    H = [K -K; -K K];
    f = repmat(ep,2*n,1) - [y; -y];
    A = [];
    b = [];
    Aeq = [ones(1,n) -ones(1,n)];
    beq = 0;
    lb = zeros(2*n,1);
    ub = repmat(C,2*n,1);
    %do the optimization:
    if exist('qld')
      alf = qld(H,f,Aeq,beq,lb,ub,[],1);
    else
      alf = quadprog(H,f,A,b,Aeq,beq,lb,ub);
    end
    % find SV's
    w = alf(1:n)-alf((n+1):end);  % the real weights
    Isv = find(abs(w)>tol);
    Ibnd = find(abs(w(Isv))<C-tol);
    % for the 'real' sv's, we want the classifier output to compute the
    % offset:
    I = Isv(Ibnd);
    out1 = sum(repmat(w',length(Ibnd),1).*K(I,:),2);
    sw = sign(w);
    out = y(I) - out1 - sw(I).*ep;

    % store the useful parameters:
    W.b = mean(out);
    W.wk = proxm(+x(Isv,:),ktype,kpar);
    W.w = w(Isv);

    y = prmapping(mfilename,'trained',W,1,d,1);
    y = setname(y,'Linear regression');
  else
    % evaluation
    [x,v] = deal(argin{1:2});
    W = getdata(v);
    [n,d] = size(x);
    Kz = +(x*W.wk);
    out = sum(repmat(W.w',n,1).*Kz,2) + repmat(W.b,n,1);
    y = setdat(x,out);

  end

return

author: marianux

repo: ecg-kit

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


text:
function elliptic_km_values_test ( )

%*****************************************************************************80
%
%% ELLIPTIC_KM_VALUES_TEST demonstrates the use of ELLIPTIC_KM_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELLIPTIC_KM_VALUES_TEST:\n' );
  fprintf ( 1, '  ELLIPTIC_KM_VALUES stores values of\n' );
  fprintf ( 1, '  the complete elliptic integral of the first\n' );
  fprintf ( 1, '  kind, with parameter modulus M.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      M            KM(M)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = elliptic_km_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function slines = readln(fid,nl)
%
% slines = readln(fid,<nlines>)
%
%
% $Id: readln.m,v 1.1 2004/03/11 19:41:26 sayres Exp $

F = fread(fid);
lF = length(F);

if( isempty(F) )
  warning('File appears to be empty');
  slines = [];
  return;
end

indNL = find(F == 10); % indicies of new lines
nNL = length(indNL);

if(max(indNL) == lF) % last char is NL
  nLines = nNL;
else
  nLines = nNL + 1;
end

if(nLines == 1)
  slines = setstr(F');
  return;
end

if(nargin == 2) nlmax = nl;
else            nlmax = nLines;
end

if(nlmax > nLines) nlmax = nLines; end

sF = setstr(F');
clear F;
slines = [];
b = 1;
for n = 1:nlmax

  % need this in case file does not end in a newline
  if(n > nNL) e = lF;
  else        e = indNL(n);
  end

  if(e > b + 1)  % skip lines with only NLs
    s = sF(b:e);
    slines = strvcat(slines,s);
  end

  if(n <= nNL) b = indNL(n) + 1; end

end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function F = exclude(X,Y)
%EXCLUDE Excludes a binary solution
%
%    F = exclude(X,value)
%
%EXCLUDE is used to avoid a particular binary solution. This can be used
% to repeatedly solve MILP problems while exluding all past solutions
%
% A = randn(30,15);
% b = 25*rand(30,1);
% c = randn(15,1);
% x = binvar(15,1);
% Model = A*x <= b;
% sol = solvesdp(Model,c'*x);
% while sol.problem == 0
%    Model = [Model, exclude(x,double(x))];
%    sol = solvesdp(Model,c'*x);
% end

if isa(X,'sdpvar') & is(X,'binary') &  isnumeric(Y) &  ismember(Y,[0 1])

    if isequal(size(X),size(Y))
    else
        error('Dimension mismatch in EXCLUDE')
    end

    zv = find((Y == 0));
    ov = find((Y == 1));
    lhs = 0;
    if ~isempty(zv)
        lhs = lhs + sum(extsubsref(X,zv));
    end
    if ~isempty(ov)
        lhs = lhs + sum(1-extsubsref(X,ov));
    end
    F = [lhs >=1];

else
    error('EXCLUDE only applicable to binary variables and data');
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function w=pde(u)
% PDE  W = PDE(U) is the residual for the convection-diffusion
%      problem in Chapter 3
%
%
global rhsf;
%
% C=20
%
v=20*u.*(dxmf(u)+dymf(u));
%
% uncomment this line for the unpreconditioned problem
%
%w=-lapmf(u)+20*u.*(dxmf(u)+dymf(u))-rhsf;
%
% uncomment this line for the preconditioned problem
%
w=u+fish2d(v)-rhsf;

repo: Aero-matlab


text:
function [] = extract_ros_data(file_name, load_folder, save_folder, cam_name)
%EXTRACT_ROS_DATA reads a ROS bag, and saves the events, IMU messages and
%camera info in a mat file.
%
% Syntax: extract_ros_data(file_name, load_folder, save_folder)
%
% You will need the matlab_rosbag package:
% https://github.com/bcharrow/matlab_rosbag
%
% Inputs:
%    file_name   - name of the bag file (no extension).
%    load_folder - folder of the bag. No trailing slash. Set to . by
%                  default.
%    save_folder - folder the mat is saved in. Set to load_folder by
%                  default.
%    cam_name    - OPTIONAL name of the event camera topic if not davis or
%                  dvs.
%
% See also GENERATE_UNDISTORT_MAP
%
% Author: Alex Zihao Zhu, University of Pennsylvania
% Email: alexzhu(at)seas.upenn.edu
% Copyright 2018 University of Pennsylvania
% Alex Zihao Zhu, Nikolay Atanasov, Kostas Daniilidis
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS, CONTRIBUTORS, AND THE
% TRUSTEES OF THE UNIVERSITY OF PENNSYLVANIA "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, CONTRIBUTORS OR THE TRUSTEES OF
% THE UNIVERSITY OF PENNSYLVANIA 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.

if nargin==1
    load_folder = '.';
end

if nargin < 3
    save_folder = load_folder;
end

if nargin == 4
    fprintf(['Warning: You have set the cam_name argument to ''%s''. ', ...
        'If this is intentional, please ignore this message.\n'], cam_name)
end

bag_name = strcat(load_folder,'/', file_name, '.bag');

fprintf('Reading ROS bag %s. This may take some time.\n', bag_name)

bag = ros.Bag.load(bag_name);

event_topics = bag.topics('.*events');

if isempty(event_topics)
    fprintf('Error, no topic with events found in bag %s\n', bag_name);
    return
end

if contains(event_topics{1}, 'davis')
    cam_name = 'davis';
elseif contains(event_topics{1}, 'dvs')
    cam_name = 'dvs';
elseif nargin < 4
    fprintf(['Could not find the default camera names ''davis'' or ''dvs'' in the bag.', ...
        'Please enter the camera name as a fourth argument to this function.\n']);
    return
end

if length(event_topics) == 1
    event_topic = ['/' cam_name '/events'];
    cinfo_topic = ['/' cam_name '/camera_info'];
    imu_topic = ['/' cam_name '/imu'];
elseif length(event_topics) == 2
    if contains(event_topics{1}, 'left') || contains(event_topics{2}, 'left')
        event_topic = ['/' cam_name '/left/events'];
        cinfo_topic = ['/' cam_name '/left/camera_info'];
        imu_topic = ['/' cam_name '/left/imu'];
    else
        fprintf('No left event topic found despite two event topics found.\n')
    end
else
    fprintf('Can''t find monocular or stereo events.\n')
end

[event_msgs, ~] = bag.readAll(event_topic);
[imu_msgs, ~] = bag.readAll(imu_topic);
[cinfo_msgs, ~] = bag.readAll(cinfo_topic);

num_event = 0;

for i=1:length(event_msgs)
    num_event = num_event + size(event_msgs{i}.events, 2);
end

cinfo = cinfo_msgs{1};

events = zeros(4, num_event);

iter = 1;

for i=1:length(event_msgs)
    events(:, iter:iter+size(event_msgs{i}.events,2)-1) = event_msgs{i}.events;
    iter = iter+size(event_msgs{i}.events,2);
end

num_imu = length(imu_msgs);

lin_acc = zeros(3, num_imu);
ang_vel = zeros(3, num_imu);
imu_time = zeros(num_imu, 1);

for i=1:num_imu
    lin_acc(:,i) = imu_msgs{i}.linear_acceleration;
    ang_vel(:,i) = imu_msgs{i}.angular_velocity;
    imu_time(i) = imu_msgs{i}.header.stamp.time;
end

output_name = strcat(save_folder,'/',file_name);

fprintf('Saving data to %s. This may take some time.\n', output_name)

save(output_name,'file_name','events',...
    'lin_acc','ang_vel','imu_time', 'cinfo', '-v7.3');
end

author: daniilidis-group


text:
function [ list_files ] = read_list_of_files( list_file_path )

list_files = {};
fid = fopen(list_file_path);
tline = fgets(fid);
c = 0;
while ischar(tline)
    c = c + 1;
    list_files{c} = tline;
%     fprintf('%s\n',list_files{c})
    tline = fgets(fid);
end
fclose(fid);
end

author: gidariss

repo: mrcnn-object-detection


text:
%
%  Context-Aware Correlation Filters
%
%  Written by Matthias Mueller, 2016
%
%  This function takes care of setting up parameters,
%  and interfacing with the online tracking benchmark.

function results=run_DCF_CA(seq, res_path, bSaveImage)

%default settings
kernel.type =  'linear';

padding = 2;  %extra area surrounding the target
lambda1 = 1e-4;     %regularization
lambda2 = 25;
interp_factor = 0.015;  %linear interpolation factor for adaptation
output_sigma_factor = 0.1; %spatial bandwidth (proportional to target)

features.gray = false;
features.hog = true;
features.hog_orientations = 9;

cell_size = 4;

target_sz = seq.init_rect(1,[4,3]);
pos = seq.init_rect(1,[2,1]) + floor(target_sz/2);
img_files = seq.s_frames;
video_path = [];

%call tracker function with all the relevant parameters
[positions, time] = tracker(video_path, img_files, pos, target_sz, ...
    padding, kernel, lambda1, lambda2, output_sigma_factor, interp_factor, ...
    cell_size, features, false);

%return results to benchmark, in a workspace variable
rects = [positions(:,2) - target_sz(2)/2, positions(:,1) - target_sz(1)/2];
rects(:,3) = target_sz(2);
rects(:,4) = target_sz(1);
fps = numel(seq.s_frames) / time;
disp(['fps: ' num2str(fps)])
results.type = 'rect';
results.res = rects;
results.fps = fps;
%assignin('base', 'res', res);

end

repo: Context-Aware-CF-Tracking


text:
% SP_DRCHLT_L2_PROJ: assign the degrees of freedom of Dirichlet boundaries through an L2 projection.
%
%   [u, dofs] = sp_drchlt_l2_proj (sp, msh, h, sides)
%
% INPUT:
%
%  sp:    object defining the space of discrete functions (see sp_scalar)
%  msh:   object defining the domain partition and the quadrature rule (see msh_cartesian)
%  h:     function handle to compute the Dirichlet condition
%  sides: boundary sides on which a Dirichlet condition is imposed
%
% OUTPUT:
%
%  u:    assigned value to the degrees of freedom
%  dofs: global numbering of the corresponding basis functions
%
% Copyright (C) 2010 Carlo de Falco
% Copyright (C) 2010, 2011, 2015 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function [u, dofs] = sp_drchlt_l2_proj (sp, msh, h, sides)

  rhs  = zeros (sp.ndof, 1);

  % In the 1D case, with an open knot vector,  it is not necessary to compute a projection.
  % For now it only works for scalars
  if (msh.ndim == 1)
    dofs = []; u = zeros (numel(sides), 1);
    for ii = 1:numel(sides)
      iside = sides(ii);
      dofs = union (dofs, sp.boundary(iside).dofs);
      if (iside == 1)
        geo_map = msh.map(msh.breaks{1}(1));
      else
        geo_map = msh.map(msh.breaks{1}(end));
      end
      for jj = 1:msh.rdim; xx{jj} = geo_map(jj); end
      u(ii) = h(xx{:}, iside);
    end
    u = u(:);
    return
  end

  dofs = [];
  nent = 0;
  for iside = sides
    nent = nent + msh.boundary(iside).nel * sp.boundary(iside).nsh_max^2;
    dofs = union (dofs, sp.boundary(iside).dofs);
  end

  rows = zeros (nent, 1);
  cols = zeros (nent, 1);
  vals = zeros (nent, 1);

  ncounter = 0;
  for iside = sides
% Restrict the function handle to the specified side, in any dimension, hside = @(x,y) h(x,y,iside)
    hside = @(varargin) h(varargin{:},iside);
    [rs, cs, vs] = op_u_v_tp (sp.boundary(iside), sp.boundary(iside), msh.boundary(iside));

    bnd_dofs = sp.boundary(iside).dofs;

    rows(ncounter+(1:numel(rs))) = bnd_dofs(rs);
    cols(ncounter+(1:numel(rs))) = bnd_dofs(cs);
    vals(ncounter+(1:numel(rs))) = vs;
    ncounter = ncounter + numel (rs);

    rhs(bnd_dofs) = rhs(bnd_dofs) + op_f_v_tp (sp.boundary(iside),msh.boundary(iside), hside);
  end

  M = sparse (rows(1:ncounter), cols(1:ncounter), vals(1:ncounter));
  u = M(dofs, dofs) \ rhs(dofs, 1);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [x,y,xy,xyp,ee] = q2p1grid(x,y,xy,mv,bound);
%q2p1grid   Q2-P1 element grid generator
%   [x,y,xy,xyp,ee] = q2p1grid(x,y,xy,mv,bound);
%   input
%          x          x coordinate vector
%          y          y coordinate vector
%          xy         nodal coordinate vector
%          mv         Q2 macroelement mapping matrix
%          bound      boundary vertex vector
%   output
%          xyp        centroid coordinate vector
%          ee         element edge connection matrix
%
%   IFISS function: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
xx=xy(:,1); yy=xy(:,2); nvtx=length(xx);
nel=length(mv(:,1));
%% recompute mid-side points in the case of stretched grids
% y-direction
yv=yy; ny=length(y);
for k=2:2:ny;
yold=y(k); ynew=0.5*(y(k+1)+y(k-1));
l=find(yy==yold); yv(l)=ynew; y(k)=ynew;
end
% x-direction
xv=xx; nx=length(x);
for k=2:2:nx;
xold=x(k); xnew=0.5*(x(k+1)+x(k-1));
l=find(xx==xold); xv(l)=xnew; x(k)=xnew;
end
xy=[xv,yv];
% centroid coordinates
for ielem=1:nel
xc(ielem)=mean(xx(mv(ielem,1:4))); yc(ielem)=mean(yy(mv(ielem,1:4)));
end
xyp=[xc',yc'];
%
%% compute edge to edge connection array ee
      np=nel;
% initialise global matrices
      adj = sparse(nvtx,nvtx);
      ee = zeros(nel,4);
%
% evaluate element number on each edge in turn
% and assemble into adjacency matrix
%% nx= 0, ny=-1
		 adj=adj + sparse(mv(:,1),mv(:,2),1:np,nvtx,nvtx);
%% nx= 1, ny= 0
		 adj=adj + sparse(mv(:,2),mv(:,3),1:np,nvtx,nvtx);
%% nx= 0, ny= 1
		 adj=adj + sparse(mv(:,3),mv(:,4),1:np,nvtx,nvtx);
%% nx=-1, ny= 0
		 adj=adj + sparse(mv(:,4),mv(:,1),1:np,nvtx,nvtx);
%
       for el=1:nel
		   [ii,jj]=find(adj==el);
           ee(el,:)=diag(adj(jj,ii))';
	   end
       ee=ee(:,[2,4,3,1]);

%
% plotting of the grid
%if nel <=64,
adj=sparse(nvtx,nvtx);
mel=length(mv(:,1));
for i=1:nel
adj(mv(i,1),mv(i,2))=1;
adj(mv(i,2),mv(i,3))=1;
adj(mv(i,3),mv(i,4))=1;
adj(mv(i,4),mv(i,1))=1;
end
figure(30)
gplot(adj,xy,'b')
axis('square')
hold on
plot(xy(:,1),xy(:,2),'ro')
xybd=xy(bound,:);
plot(xybd(:,1),xybd(:,2),'ko')
plot(xyp(:,1),xyp(:,2),'k*',xyp(:,1),xyp(:,2),'ro')
hold off
title('Q2-P1 finite element subdivision')
drawnow
%end
return


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

    properties
        %Tolerances
        tolCon(1,1) double = 1E-6;

        %Maximums
        maxFuncEvals(1,1) double = NaN;

        %Parallel
        useParallel(1,1) SurrogateOptUseParallelEnum = SurrogateOptUseParallelEnum.DoNotUseParallel;
        numWorkers(1,1) double = feature('numCores');

        %Other Options
        batchUpdateInterval(1,1) double = NaN;
        minSampleDistance(1,1) double = NaN
        minSurrogatePoints(1,1) double = NaN;
    end

    methods
        function obj = SurrogateOptimizerOptions()

        end

        function options = getOptionsForOptimizer(obj, x0)
            options = optimoptions(@surrogateopt, 'Display','iter');

            x0 = x0(:)'; %needs to be a row
            options = optimoptions(options, 'InitialPoints', x0);

            %%%

            if(not(isnan(obj.tolCon)))
                options = optimoptions(options, 'ConstraintTolerance', obj.tolCon);
            end

            if(not(isnan(obj.maxFuncEvals)))
                options = optimoptions(options, 'MaxFunctionEvaluations', obj.maxFuncEvals);
            end

            %%%

            if(not(isnan(obj.batchUpdateInterval)))
                options = optimoptions(options, 'BatchUpdateInterval', obj.batchUpdateInterval);
            end

            if(not(isnan(obj.minSampleDistance)))
                options = optimoptions(options, 'MinSampleDistance', obj.minSampleDistance);
            end

            if(not(isnan(obj.minSurrogatePoints)))
                options = optimoptions(options, 'MinSurrogatePoints', obj.minSurrogatePoints);
            end

            options = optimoptions(options, ...
                                   'UseVectorized', false, ...
                                   'PlotFcn', {}, ...
                                   'UseParallel', obj.useParallel.optionVal);
        end

        function tf = usesParallel(obj)
            tf = obj.useParallel;
        end

        function numWorkers = getNumParaWorkers(obj)
            numWorkers = obj.numWorkers;
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            if(obj.numWorkers < 1 || obj.numWorkers > feature('numCores'))
                obj.numWorkers = feature('numCores');
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function hDataTip = createDataTipROE(hFig,target)
% Create new datatip or expand existing one.
%
% AI 10/15/21

ud = guidata(hFig);
cMode = datacursormode(hFig);
cursorInfoS = getCursorInfo(cMode);
%target = protS(pNum).(selTypeC{k})(selNum);

%Check if datatip exists
if ~isfield(ud.handle,'datatips')
    ud.handle.datatips = [];
end

evt.Target = target;
evt.Position = [evt.Target.XData(1),evt.Target.YData(1)];

existing = 0;
for c = 1:length(cursorInfoS)
    if isequal(get(cursorInfoS(c).Target,'userdata'),...
            get(evt.Target,'userdata'))
        hDataTip = ud.handle.datatips(c);
        guidata(hFig,ud);
        expandDataTipROE(hDataTip,evt,hFig);
        existing = 1;
        break
    end
end
if ~existing
    hDataTip = cMode.createDatatip(target);
    set(hDataTip,'Marker','^','OrientationMode','Manual',...
        'Tag','criteria','userdata',get(target,'userdata'),'UpdateFcn',...
        @(hObj,hEvt)expandDataTipROE(hObj,hEvt,hFig)); %'MarkerSize,',7,
    ud.handle.datatips = [ud.handle.datatips,hDataTip];
    guidata(hFig,ud);
    expandDataTipROE(hDataTip,evt,hFig);
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%THIS FUMCTION ILLUSTRATES HOW TO USE THE USER-DEFINED ALGORITHM
function [AH,XH] = user_alg3(Y,r,X)
%
% The example of implementing the user-defined algorithm (this is the regularized Fixed-Point algorithm based)
%
% INPUTS:
% Y - mixed signals (matrix of size [m by T])
% r - number of estimated signals
% X - true source signals
%
% OUTPUTS
% AH - estimated mixing matrix (matrix of size [m by r])
% XH - estimated source signals (matrix of size [r by T])
%
% #########################################################################
% Initialization
[m,T]=size(Y);
Y(Y <=0) = eps; % this enforces the positive value in the data
AH=rand(m,r);
XH=rand(r,T);

IterNo = 1000; % number of alternating steps

% Iterations
for k = 1:IterNo

    alpha_reg = 20*exp(-k/10); % regularization parameter
    XH = max(1E6*eps,pinv(AH'*AH +  alpha_reg)*AH'*Y);
    AH = max(1E6*eps, Y*XH'*pinv(XH*XH' + alpha_reg));
    AH = AH*diag(1./(sum(AH,1) + eps));

end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
%   - data                 3d-brains, Omega=(0,20)x(0,10)x(0,20), level=3:6, m=[128,64,128]
%   - viewer               imgmontage
%   - interpolation        linearInterMex
%   - distance             SSD
%   - pre-registration     affine3D
%   - regularizer          mfElastic
%   - optimizer            Gauss-Newton
% ===============================================================================

setup3DbrainData;

% extract data of level 4
level = 4; omega = ML{level}.omega; m = ML{level}.m;
imgModel('reset','imgModel','linearInterMex');
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);
xc = getCellCenteredGrid(omega,m);
Rc = imgModel(R,omega,xc);

% initialize distance measure
distance('reset','distance','SSD');

% initialize the transformation and a starting guess
% trafo('reset','trafo','affine3Dsparse');
trafo('reset','trafo','splineTransformation3Dsparse',...
  'omega',omega,'p',[ 3 4 5],'m',m);
w0 = trafo('w0');

% setup plots and initialize
FAIRplots('reset','mode','PIR-GN','omega',omega,'m',m,'fig',1,'plots',1);
FAIRplots('init',struct('Tc',T,'Rc',R,'omega',omega,'m',m));

%------------------------------------------------------------------------------
% build objective function
% note: T  is data for template image
%       Rc is sampled reference image
%       optional Tikhonov-regularization is disabled by setting m = [], wRef = []
%       beta = 0, M = [], wRef = []:
%       disables additional regularization of Hessian approximation
beta = 0; M = []; wRef = [];
fctn = @(wc) PIRobjFctn(T,Rc,omega,m,beta,M,wRef,xc,wc);
fctn([]);   % report status

%------------------------------------------------------------------------------
%% -- solve the optimization problem on one level
[wc,his] = GaussNewton(fctn,w0,'solver','CG','Plots',@FAIRplots); return;

%------------------------------------------------------------------------------
%% finally: run the MultiLevel Non-Parametric Image Registration
[wc,his] = MLPIR(ML);

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


text:
function d = det2x2(x)

% DET2X2 computes determinant of matrix x, using explicit analytic definition
% if size(x,1) < 4, otherwise use MATLAB det-function

% Copyright (C) 2012, Donders Centre for Cognitive Neuroimaging, Nijmegen, NL
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

siz = size(x);
if all(siz(1:2)==2),
  d = x(1,1,:,:).*x(2,2,:,:) - x(1,2,:,:).*x(2,1,:,:);
elseif all(siz(1:2)==3),
  d = x(1,1,:,:).*x(2,2,:,:).*x(3,3,:,:) - ...
      x(1,1,:,:).*x(2,3,:,:).*x(3,2,:,:) - ...
      x(1,2,:,:).*x(2,1,:,:).*x(3,3,:,:) + ...
      x(1,2,:,:).*x(2,3,:,:).*x(3,1,:,:) + ...
      x(1,3,:,:).*x(2,1,:,:).*x(3,2,:,:) - ...
      x(1,3,:,:).*x(2,2,:,:).*x(3,1,:,:);
elseif numel(siz)==2,
  d = det(x);
else
  ft_error('not implemented');
  % write for loop for the higher dimensions, using normal inv
end

author: spm


text:
function inside = triangle_contains_point_1 ( t, p )

%*****************************************************************************80
%
%% TRIANGLE_CONTAINS_POINT_1 finds if a point is inside a triangle.
%
%  Discussion:
%
%    It is conventional to list the triangle vertices in counter clockwise
%    order.  However, this routine does not require a particular order
%    for the vertices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, logical INSIDE, is TRUE if the point is inside
%    the triangle or on its boundary.
%
  xsi = triangle_barycentric ( t, p );

  if ( any ( xsi(1:3,1) < 0.0 ) )
    inside = 0;
  else
    inside = 1;
  end

  return
end


text:
function [G, Y, optinf] = bpdndl(D0, S, lambda, opt)

% bpdndl -- BPDN Dictionary Learning
%
%         argmin_{D,X} (1/2)||D X - S||_2^2 + lambda ||X||_1
%
%         Dictionary learning consists of interleaved alternation of
%         the ADMM (see boyd-2010-distributed) steps for the BPDN (see
%         chen-1998-atomic) and MOD (see engan-1999-method) problems.
%
% Usage:
%       [D, X, optinf] = bpdndl(D0, S, lambda, opt)
%
% Input:
%       D0          Initial dictionary
%       S           Input image
%       lambda      Regularization parameter
%       opt         Options/algorithm parameters structure (see below)
%
% Output:
%       D           Dictionary
%       X           Coefficients
%       optinf      Details of optimisation
%
%
% Options structure fields:
%   Verbose          Flag determining whether iteration status is displayed.
%                    Fields are iteration number, functional value,
%                    data fidelity term, l1 regularisation term, and
%                    primal and dual residuals (see Sec. 3.3 of
%                    boyd-2010-distributed). The values of rho and sigma
%                    are also displayed if options request that they are
%                    automatically adjusted.
%   MaxMainIter      Maximum main iterations
%   AbsStopTol       Absolute convergence tolerance (see Sec. 3.3.1 of
%                    boyd-2010-distributed)
%   RelStopTol       Relative convergence tolerance (see Sec. 3.3.1 of
%                    boyd-2010-distributed)
%   L1Weight         Weight matrix for L1 norm
%   Y0               Initial value for Y
%   U0               Initial value for U
%   G0               Initial value for G (overrides D0 if specified)
%   H0               Initial value for H
%   rho              Augmented Lagrangian penalty parameter
%   AutoRho          Flag determining whether rho is automatically updated
%                    (see Sec. 3.4.1 of boyd-2010-distributed)
%   AutoRhoPeriod    Iteration period on which rho is updated
%   RhoRsdlRatio     Primal/dual residual ratio in rho update test
%   RhoScaling       Multiplier applied to rho when updated
%   AutoRhoScaling   Flag determining whether RhoScaling value is
%                    adaptively determined (see wohlberg-2015-adaptive). If
%                    enabled, RhoScaling specifies a maximum allowed
%                    multiplier instead of a fixed multiplier
%   sigma            Augmented Lagrangian penalty parameter
%   AutoSigma        Flag determining whether sigma is automatically
%                    updated (see Sec. 3.4.1 of boyd-2010-distributed)
%   AutoSigmaPeriod  Iteration period on which sigma is updated
%   SigmaRsdlRatio   Primal/dual residual ratio in sigma update test
%   SigmaScaling     Multiplier applied to sigma when updated
%   AutoSigmaScaling Flag determining whether SigmaScaling value is
%                    adaptively determined (see wohlberg-2015-adaptive). If
%                    enabled, SigmaScaling specifies a maximum allowed
%                    multiplier instead of a fixed multiplier.
%   StdResiduals     Flag determining whether standard residual definitions
%                    (see Sec 3.3 of boyd-2010-distributed) are used instead
%                    of normalised residuals (see wohlberg-2015-adaptive)
%   XRelaxParam      Relaxation parameter (see Sec. 3.4.3 of
%                    boyd-2010-distributed) for X update
%   DRelaxParam      Relaxation parameter (see Sec. 3.4.3 of
%                    boyd-2010-distributed) for D update
%   NonNegCoef       Flag indicating whether solution should be forced to
%                    be non-negative
%   AuxVarObj        Flag determining whether objective function is computed
%                    using the auxiliary (split) variable
%   ZeroMean         Force learned dictionary entries to be zero-mean
%
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2015-07-30
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'License' file distributed with
% the library.

if nargin < 4,
  opt = [];
end
checkopt(opt, defaultopts([]));
opt = defaultopts(opt);
Nx = size(D0,2)*size(S,2);
Nd = prod(size(D0));

% Set up status display for verbose operation
hstr = ['Itn   Fnc       DFid      l1        Cnstr     '...
        'r(X)      s(X)      r(D)      s(D) '];
sfms = '%4d %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e';
nsep = 84;
if opt.AutoRho,
  hstr = [hstr '     rho  '];
  sfms = [sfms ' %9.2e'];
  nsep = nsep + 10;
end
if opt.AutoSigma,
  hstr = [hstr '     sigma  '];
  sfms = [sfms ' %9.2e'];
  nsep = nsep + 10;
end
if opt.Verbose && opt.MaxMainIter > 0,
  disp(hstr);
  disp(char('-' * ones(1,nsep)));
end

% Mean removal and normalisation projections
Pzmn = @(x) bsxfun(@minus, x, mean(x,1));
Pnrm = @(x) normalise(x);

% Projection of dictionary filters onto constraint set
if opt.ZeroMean,
  Pcn = @(x) Pnrm(Pzmn(x));
else
  Pcn = @(x) Pnrm(x);
end

% Start timer
tstart = tic;

% Project initial dictionary onto constraint set
D = Pnrm(D0);

% Set up algorithm parameters and initialise variables
rho = opt.rho;
if isempty(rho), rho = 50*lambda+1; end;
sigma = opt.sigma;
if isempty(sigma), sigma = size(S,2)/200; end;
optinf = struct('itstat', [], 'opt', opt);
rx = Inf;
sx = Inf;
rd = Inf;
sd = Inf;
eprix = 0;
eduax = 0;
eprid = 0;
eduad = 0;

% Initialise main working variables
X = [];
if isempty(opt.Y0),
  Y = zeros(size(D,2), size(S,2));
else
  Y = opt.Y0;
end
Yprv = Y;
if isempty(opt.U0),
  if isempty(opt.Y0),
    U = zeros(size(D,2), size(S,2), class(S));
  else
    U = (lambda/rho)*sign(Y);
  end
else
  U = opt.U0;
end
if isempty(opt.G0),
  G = D;
else
  G = opt.G0;
end
Gprv = G;
if isempty(opt.H0),
  if isempty(opt.G0),
    H = zeros(size(G), class(S));
  else
    H = G;
  end
else
  H = opt.H0;
end
GS = G'*S;

% Main loop
k = 1;
while k <= opt.MaxMainIter && (rx > eprix|sx > eduax|rd > eprid|sd >eduad),

  % Solve X subproblem, using G as the dictionary for improved stability
  [luLx, luUx] = factorise(G, rho);
  X = linsolveX(G, rho, luLx, luUx, GS + rho*(Y - U));

  % See pg. 21 of boyd-2010-distributed
  if opt.XRelaxParam == 1,
    Xr = X;
  else
    Xr = opt.XRelaxParam*X + (1-opt.XRelaxParam)*Y;
  end

  % Solve Y subproblem
  Y = shrink(Xr + U, (lambda/rho)*opt.L1Weight);
  if opt.NonNegCoef,
    Y(Y < 0) = 0;
  end
  SY = S*Y';

  % Update dual variable corresponding to X, Y
  U = U + Xr - Y;

  % Compute primal and dual residuals and stopping thresholds for X update
  nX = norm(X(:)); nY = norm(Y(:)); nU = norm(U(:));
  if opt.StdResiduals,
    % See pp. 19-20 of boyd-2010-distributed
    rx = norm(vec(X - Y));
    sx = norm(vec(rho*(Yprv - Y)));
    eprix = sqrt(Nx)*opt.AbsStopTol+max(nX,nY)*opt.RelStopTol;
    eduax = sqrt(Nx)*opt.AbsStopTol+rho*nU*opt.RelStopTol;
  else
    % See wohlberg-2015-adaptive
    rx = norm(vec(X - Y))/max(nX,nY);
    sx = norm(vec(Yprv - Y))/nU;
    eprix = sqrt(Nx)*opt.AbsStopTol/max(nX,nY)+opt.RelStopTol;
    eduax = sqrt(Nx)*opt.AbsStopTol/(rho*nU)+opt.RelStopTol;
  end

  % Solve D subproblem, using Y as the coefficients for improved stability
  [luLd, luUd] = factorise(Y, sigma);
  D = linsolveD(Y, sigma, luLd, luUd, SY + sigma*(G - H));

  % See pg. 21 of boyd-2010-distributed
  if opt.DRelaxParam == 1,
    Dr = D;
  else
    Dr = opt.DRelaxParam*D + (1-opt.DRelaxParam)*G;
  end

  % Solve G subproblem
  G = Pcn(Dr + H);
  GS = G'*S;

  % Update dual variable corresponding to D, G
  H = H + Dr - G;

  % Compute primal and dual residuals and stopping thresholds for D update
  nD = norm(D(:)); nG = norm(G(:)); nH = norm(H(:));
  if opt.StdResiduals,
    % See pp. 19-20 of boyd-2010-distributed
    rd = norm(vec(D - G));
    sd = norm(vec(sigma*(Gprv - G)));
    eprid = sqrt(Nd)*opt.AbsStopTol+max(nD,nG)*opt.RelStopTol;
    eduad = sqrt(Nd)*opt.AbsStopTol+sigma*nH*opt.RelStopTol;
  else
    % See wohlberg-2015-adaptive
    rd = norm(vec(D - G))/max(nD,nG);
    sd = norm(vec(Gprv - G))/nH;
    eprid = sqrt(Nd)*opt.AbsStopTol/max(nD,nG)+opt.RelStopTol;
    eduad = sqrt(Nd)*opt.AbsStopTol/(sigma*nH)+opt.RelStopTol;
  end

  % Objective function
  if opt.AuxVarObj,
    Jdf = sum(vec(abs(G*Y - S).^2))/2;
    Jl1 = sum(abs(vec(opt.L1Weight .* Y)));
  else
    Jdf = sum(vec(abs(D*X - S).^2))/2;
    Jl1 = sum(abs(vec(opt.L1Weight .* X)));
  end
  Jfn = Jdf + lambda*Jl1;
  Jcn = norm(vec(Pcn(D) - D));

  % Record and display iteration details
  tk = toc(tstart);
  optinf.itstat = [optinf.itstat; ...
        [k Jfn Jdf Jl1 rx sx rd sd eprix eduax eprid eduad rho sigma tk]];
  if opt.Verbose,
    dvc = [k Jfn Jdf Jl1 Jcn rx sx rd sd];
    if opt.AutoRho,
      dvc = [dvc rho];
    end
    if opt.AutoSigma,
      dvc = [dvc sigma];
    end
    disp(sprintf(sfms, dvc));
  end

  % See wohlberg-2015-adaptive and pp. 20-21 of boyd-2010-distributed
  if opt.AutoRho,
    if k ~= 1 && mod(k, opt.AutoRhoPeriod) == 0,
      if opt.AutoRhoScaling,
        rhomlt = sqrt(rx/sx);
        if rhomlt < 1, rhomlt = 1/rhomlt; end
        if rhomlt > opt.RhoScaling, rhomlt = opt.RhoScaling; end
      else
        rhomlt = opt.RhoScaling;
      end
      rsf = 1;
      if rx > opt.RhoRsdlRatio*sx, rsf = rhomlt; end
      if sx > opt.RhoRsdlRatio*rx, rsf = 1/rhomlt; end
      rho = rsf*rho;
      U = U/rsf;
    end
  end
  if opt.AutoSigma,
    if k ~= 1 && mod(k, opt.AutoSigmaPeriod) == 0,
      if opt.AutoSigmaScaling,
        sigmlt = sqrt(rd/sd);
        if sigmlt < 1, sigmlt = 1/sigmlt; end
        if sigmlt > opt.SigmaScaling, sigmlt = opt.SigmaScaling; end
      else
        sigmlt = opt.SigmaScaling;
      end
      ssf = 1;
      if rd > opt.SigmaRsdlRatio*sd, ssf = sigmlt; end
      if sd > opt.SigmaRsdlRatio*rd, ssf = 1/sigmlt; end
      sigma = ssf*sigma;
      H = H/ssf;
    end
  end

  Yprv = Y;
  Gprv = G;
  k = k + 1;

end

% Record run time and working variables
optinf.runtime = toc(tstart);
optinf.X = X;
optinf.Y = Y;
optinf.U = U;
optinf.D = D;
optinf.G = G;
optinf.H = H;
optinf.lambda = lambda;
optinf.rho = rho;
optinf.sigma = sigma;

if opt.Verbose && opt.MaxMainIter > 0,
  disp(char('-' * ones(1,nsep)));
end

return

function u = vec(v)

  u = v(:);

return

function u = shrink(v, lambda)

  u = sign(v).*max(0, abs(v) - lambda);

return

function u = normalise(v)

  vn = sqrt(sum(v.^2, 1));
  vn(vn == 0) = 1;
  u = bsxfun(@rdivide, v, vn);

return

function [L,U] = factorise(A, c)

  [N,M] = size(A);
  % If N < M it is cheaper to factorise A*A' + cI and then use the
  % matrix inversion lemma to compute the inverse of A'*A + cI
  if N >= M,
    [L,U] = lu(A'*A + c*eye(M,M));
  else
    [L,U] = lu(A*A' + c*eye(N,N));
  end

return

function x = linsolveX(A, c, L, U, b)

  [N,M] = size(A);
  if N >= M,
    x = U \ (L \ b);
  else
    x = (b - A'*(U \ (L \ (A*b))))/c;
  end

return

function x = linsolveD(A, c, L, U, b)

  [N,M] = size(A);
  if N >= M,
    x = (b - (((b*A) / U) / L)*A')/c;
  else
    x = (b / U) / L;
  end

return

function opt = defaultopts(opt)

  if ~isfield(opt,'Verbose'),
    opt.Verbose = 0;
  end
  if ~isfield(opt,'MaxMainIter'),
    opt.MaxMainIter = 1000;
  end
  if ~isfield(opt,'AbsStopTol'),
    opt.AbsStopTol = 1e-6;
  end
  if ~isfield(opt,'RelStopTol'),
    opt.RelStopTol = 1e-4;
  end
  if ~isfield(opt,'L1Weight'),
    opt.L1Weight = 1;
  end
  if ~isfield(opt,'Y0'),
    opt.Y0 = [];
  end
  if ~isfield(opt,'U0'),
    opt.U0 = [];
  end
  if ~isfield(opt,'G0'),
    opt.G0 = [];
  end
  if ~isfield(opt,'H0'),
    opt.H0 = [];
  end
  if ~isfield(opt,'rho'),
    opt.rho = [];
  end
  if ~isfield(opt,'AutoRho'),
    opt.AutoRho = 0;
  end
  if ~isfield(opt,'AutoRhoPeriod'),
    opt.AutoRhoPeriod = 10;
  end
  if ~isfield(opt,'RhoRsdlRatio'),
    opt.RhoRsdlRatio = 10;
  end
  if ~isfield(opt,'RhoScaling'),
    opt.RhoScaling = 2;
  end
  if ~isfield(opt,'AutoRhoScaling'),
    opt.AutoRhoScaling = 0;
  end
  if ~isfield(opt,'sigma'),
    opt.sigma = [];
  end
  if ~isfield(opt,'AutoSigma'),
    opt.AutoSigma = 0;
  end
  if ~isfield(opt,'AutoSigmaPeriod'),
    opt.AutoSigmaPeriod = 10;
  end
  if ~isfield(opt,'SigmaRsdlRatio'),
    opt.SigmaRsdlRatio = 10;
  end
  if ~isfield(opt,'SigmaScaling'),
    opt.SigmaScaling = 2;
  end
  if ~isfield(opt,'AutoSigmaScaling'),
    opt.AutoSigmaScaling = 0;
  end
  if ~isfield(opt,'StdResiduals'),
    opt.StdResiduals = 0;
  end
  if ~isfield(opt,'XRelaxParam'),
    opt.XRelaxParam = 1;
  end
  if ~isfield(opt,'DRelaxParam'),
    opt.DRelaxParam = 1;
  end
  if ~isfield(opt,'NonNegCoef'),
    opt.NonNegCoef = 0;
  end
  if ~isfield(opt,'AuxVarObj'),
    opt.AuxVarObj = 1;
  end
  if ~isfield(opt,'ZeroMean'),
    opt.ZeroMean = 0;
  end

return

author: thfylsty


text:
function [path,name,ext] = myfileparts(s)

path = '';
name = '';
if isempty(s),
  return;
end

if s(end) == '/' || ispc && s(end) == '\',
  s = s(1:end-1);
end

[path,name,ext] = fileparts_platind(s);
if nargout < 3,
  name = [name,ext];
end

author: kristinbranson

repo: JAABA


text:
function pass = test_solharm( pref )
% Test the equality for P^m_l, 0 <= l <= n and -m <= l <= m

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e6*pref.techPrefs.chebfuneps;

n = 10;
Max_difference = 0;
Max_norm = 0;
Max_laplacian = 0;

for l = 0:n
    for m = -l:l
        Yml = ballfun.solharm(l,m);
        Y = Yml(1,:,:,'spherical')/sqrt(2*l+3);
        Z = spherefun.sphharm(l,m);
        % Check that the Legendre polynomials is the same as in Spherefun
        Max_difference = max(norm(Y-Z),Max_difference);
        % Check that the 2-norm is 1
        Max_norm = max(abs(norm(Yml)-1),Max_norm);
        % Check that they are solutions to the Laplace equation
        Max_laplacian = max(norm(laplacian(Yml)),Max_laplacian);
    end
end
pass(1) = Max_difference < tol;
pass(2) = Max_norm < tol;
pass(3) = Max_laplacian < tol;

% Example 1 : Y^0_0
m = 0; l = 0;
Normalization = sqrt(2*l+3);
Exact = ballfun(@(r,lam,th)0.5*sqrt(1/pi)*Normalization, 'spherical');
f = ballfun.solharm(l,m,'complex');
pass(4) = norm(Exact - f) < tol;

% Example 2 : Y^-1_1
m = -1; l = 1;
Normalization = sqrt(2*l+3);
Exact = ballfun(@(r,lam,th)0.5*sqrt(3/(2*pi))*r.^l.*exp(1i*m*lam).*sin(th)*Normalization, 'spherical');
f = ballfun.solharm(l,m,'complex');
pass(5) = norm(Exact - f) < tol;

% Example 4 : Y^0_1
m = 0; l = 1;
Normalization = sqrt(2*l+3);
Exact = ballfun(@(r,lam,th)0.5*sqrt(3/pi)*r.^l.*exp(1i*m*lam).*cos(th)*Normalization, 'spherical');
f = ballfun.solharm(l,m,'complex');
pass(6) = norm(Exact - f) < tol;

% Example 5 : Y^1_1
m = 1; l = 1;
Normalization = sqrt(2*l+3);
Exact = ballfun(@(r,lam,th)-0.5*sqrt(3/(2*pi))*r.^l.*exp(1i*m*lam).*sin(th)*Normalization, 'spherical');
f = ballfun.solharm(l,m,'complex');
pass(7) = norm(Exact - f) < tol;

% Example 6 : Y^-2_2
m = -2; l = 2;
Normalization = sqrt(2*l+3);
Exact = ballfun(@(r,lam,th)0.25*sqrt(15/(2*pi))*r.^l.*exp(1i*m*lam).*sin(th).^2*Normalization, 'spherical');
f = ballfun.solharm(l,m,'complex');
pass(8) = norm(Exact - f) < tol;

% Example 7 : Y^-1_2
m = -1; l = 2;
Normalization = sqrt(2*l+3);
Exact = ballfun(@(r,lam,th)0.5*sqrt(15/(2*pi))*r.^l.*exp(1i*m*lam).*sin(th).*cos(th)*Normalization, 'spherical');
f = ballfun.solharm(l,m,'complex');
pass(9) = norm(Exact - f) < tol;

% Example 8 : Y^0_2
m = 0; l = 2;
Normalization = sqrt(2*l+3);
Exact = ballfun(@(r,lam,th)0.25*sqrt(5/pi)*r.^l.*exp(1i*m*lam).*(3*cos(th).^2-1)*Normalization, 'spherical');
f = ballfun.solharm(l,m,'complex');
pass(10) = norm(Exact - f) < tol;

% Example 9 : Y^1_2
m = 1; l = 2;
Normalization = sqrt(2*l+3);
Exact = ballfun(@(r,lam,th)-0.5*sqrt(15/(2*pi))*r.^l.*exp(1i*m*lam).*sin(th).*cos(th)*Normalization, 'spherical');
f = ballfun.solharm(l,m,'complex');
pass(11) = norm(Exact - f) < tol;

% Example 10 : Y^2_2
m = 2; l = 2;
Normalization = sqrt(2*l+3);
Exact = ballfun(@(r,lam,th)0.25*sqrt(15/(2*pi))*r.^l.*exp(1i*m*lam).*sin(th).^2*Normalization, 'spherical');
f = ballfun.solharm(l,m,'complex');
pass(12) = norm(Exact - f) < tol;

% Example 11: real solid harmonics
f = ballfun.solharm(3,2);
pass(13) = abs(norm(f)-1) < tol;

% Example 12: Y^0_5
f = ballfun.solharm(5,0);
pass(14) = abs(norm(f)-1) < tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function psnr= PSNR(X,Y)
%Calculates the Peak-to-peak Signal to Noise Ratio of two images X and Y
[M,N]=size(X);
m=double(0);
X=cast(X,'double');
Y=cast(Y,'double');
for i=1:M
    for j=1:N
        m=m+((X(i,j)-Y(i,j))^2);
    end
end
m=m/(M*N);
psnr=10*log10(255*255/m);


text:
function out=prox_Huber(x,mu,alpha)
%PROX_HUBER computes the proximal operator of the function alpha*H_(mu)
%                      where H_(mu)=huber function with parameter mu
%
%  Usage:
%  out = PROX_HUBER(x,mu,alpha)
%  ===========================================
%  INPUT:
%  x - point to be projected (vector/matrix)
%  mu - positive scalar
%  alpha - positive scalar
%  ===========================================
%  Assumptions:
%  mu is positive
%  ===========================================
%  Output:
%  out - proximal operator at x

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

if (nargin < 3)
    error ('usage: prox_huber(x,mu,alpha)') ;
end

if (alpha < 0)
    error('usage: prox_huber(x,mu,alpha) - alpha should be positive')
end

if (mu < 0)
    error('usage: prox_huber(x,mu,alpha) - mu should be positive')
end

out = x*( 1  - alpha/max(norm(x,'fro'),mu+alpha))   ;

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Dh=hammingDist(B1, B2)
%
% Written by Rob Fergus
% Compute hamming distance between two sets of samples (B1, B2)
%
% Dh=hammingDist(B1, B2);
%
% Input
%    B1, B2: compact bit vectors. Each datapoint is one row.
%    size(B1) = [ndatapoints1, nwords]
%    size(B2) = [ndatapoints2, nwords]
%    It is faster if ndatapoints1 < ndatapoints2
%
% Output
%    Dh = hamming distance.
%    size(Dh) = [ndatapoints1, ndatapoints2]
%
% example query
% Dhamm = hammingDist(B2, B1);
% this will give the same result than:
%    Dhamm = distMat(U2>0, U1>0).^2;
% the size of the distance matrix is:
%    size(Dhamm) = [Ntest x Ntraining]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% loop-up table:
bit_in_char = uint16([...
    0 1 1 2 1 2 2 3 1 2 2 3 2 3 3 4 1 2 2 3 2 3 ...
    3 4 2 3 3 4 3 4 4 5 1 2 2 3 2 3 3 4 2 3 3 4 ...
    3 4 4 5 2 3 3 4 3 4 4 5 3 4 4 5 4 5 5 6 1 2 ...
    2 3 2 3 3 4 2 3 3 4 3 4 4 5 2 3 3 4 3 4 4 5 ...
    3 4 4 5 4 5 5 6 2 3 3 4 3 4 4 5 3 4 4 5 4 5 ...
    5 6 3 4 4 5 4 5 5 6 4 5 5 6 5 6 6 7 1 2 2 3 ...
    2 3 3 4 2 3 3 4 3 4 4 5 2 3 3 4 3 4 4 5 3 4 ...
    4 5 4 5 5 6 2 3 3 4 3 4 4 5 3 4 4 5 4 5 5 6 ...
    3 4 4 5 4 5 5 6 4 5 5 6 5 6 6 7 2 3 3 4 3 4 ...
    4 5 3 4 4 5 4 5 5 6 3 4 4 5 4 5 5 6 4 5 5 6 ...
    5 6 6 7 3 4 4 5 4 5 5 6 4 5 5 6 5 6 6 7 4 5 ...
    5 6 5 6 6 7 5 6 6 7 6 7 7 8]);

n1 = size(B1,1);
[n2, nwords] = size(B2);

Dh = zeros([n1 n2], 'uint16');
for j = 1:n1
    for n=1:nwords
        y = bitxor(B1(j,n),B2(:,n));
        Dh(j,:) = Dh(j,:) + bit_in_char(y+1);
    end
end

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
classdef ODE45Integrator < AbstractFirstOrderIntegrator
    %ODE45Integrator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        options(1,1) BuiltInIntegratorOptions = BuiltInIntegratorOptions();
        integratorEnum = IntegratorEnum.ODE45;
    end

    methods
        function obj = ODE45Integrator(options)
            if(nargin > 0)
                obj.options = options;
            else
                obj.options = BuiltInIntegratorOptions();
            end
        end

        function [t,y,te,ye,ie] = integrate(obj, odefun, tspan, y0, evtsFunc, odeOutputFun)
            odeSetOptions = obj.options.getIntegratorOptions();
            optionsToUse = odeset(odeSetOptions, 'Events',evtsFunc, 'OutputFcn',odeOutputFun);

            [t,y,te,ye,ie] = ode45(odefun, tspan, y0, optionsToUse);
        end

        function options = getOptions(obj)
            options = obj.options;
        end
    end
end


text:
function [y,dzdx2] = vl_nnscosinesim(x1, x2, varargin)
% VL_NNCOSINESIM Compute the cosine similariy between features
%   Y = VL_NNCOSINESIM(X1, X2) computes the cosine similarity between X1
%   and X2 where X1 has shape H x W x C x N, X2 has shape
%   H x W x C x N and Y is has shape 1 x 1 x 1 x N (producing a scalar
%   per batch element).
%
%   [DZDX1, DZDX2] = VL_NNCOSINESIM(X1, X2, DZDY) computes the
%   derivatives of the block projected onto DZDY. DZDX1 and DZDX2
%   have the same dimensions as X1, X2 and B respectively.
%
%   NOTES: The cosine similarity is defined between vectors A and B of
%   dimension N x 1 to be the scalar given by:
%   cos_sim = a'b / (norm(a,2) * norm(b,2))
%   To operate with tensors, the first three channels (H, W and C) are
%   treated as the vector elements, while the fourth dimension denotes
%   batch indicies.
%
%   VL_NNCOSINELOSS(.., 'option', value, ...) accepts the following options:
%
%   `eps`:: 1e-6
%    A small value to avoid possible division by zero.
%
% Copyright (C) 2018 Samuel Albanie.
% Licensed under The MIT License [see LICENSE.md for details]

  opts.eps = 1e-6 ;
  [~, dzdy] = vl_argparsepos(struct(), varargin) ;

  sz1 = size(x1) ; sz2 = size(x2) ;
  assert(all(sz1 == sz2), 'tensor sizes do not match') ;
  x1 = reshape(x1, [], sz1(4)) ;
  x2 = reshape(x2, [], sz1(4)) ;
  dots = dot(x1, x2, 1) ;
  x1_norms = sqrt(sum(x1 .* x1, 1)) ;
  x2_norms = sqrt(sum(x2 .* x2, 1)) ;

  if isempty(dzdy)
    y = dots ./ max(x1_norms .* x2_norms, opts.eps) ;
    y = reshape(y, 1, 1, 1, sz1(4)) ;
  else
    dzdy = dzdy{1} ; dsize = size(dzdy) ;
    assert(numel(dsize) == 4 & all(dsize(1:3) == ones(1,3)) ...
           & dsize(4) == sz1(4), 'DZDY has an unexpected size') ;
    t1 = bsxfun(@rdivide, x2, max(x1_norms .* x2_norms, opts.eps)) ;
    t2 = bsxfun(@rdivide, bsxfun(@times, x1, dots), ...
                              max(x1_norms .^3 .* x2_norms, opts.eps)) ;
    dzdx1 = t1 - t2 ;

    t1 = bsxfun(@rdivide, x1, max(x2_norms .* x1_norms, opts.eps)) ;
    t2 = bsxfun(@rdivide, bsxfun(@times, x2, dots), ...
                              max(x2_norms .^3 .* x1_norms, opts.eps)) ;
    dzdx2 = t1 - t2 ;

    % reshape to match input and compute projected derivatives
    dzdx1 = reshape(dzdx1, sz1) ;
    dzdx2 = reshape(dzdx2, sz2) ;
    dzdx1 = bsxfun(@times, dzdx1, dzdy) ;
    dzdx2 = bsxfun(@times, dzdx2, dzdy) ;
    y = dzdx1 ;
  end

repo: MFT


text:
function npart = setpart_enum ( m )

%*****************************************************************************80
%
%% SETPART_ENUM enumerates the partitions of a set of M elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer M, the number of elements in the set.
%    M must be positive.  However, for the enumeration routine only,
%    it is legal to call with any value of M.
%
%    Output, integer NPART, the number of partitions of the set.
%
  if ( m < 0 )

    npart = 0;

  elseif ( m == 0 )

    npart = 1;

  else

    offset = 1;
    b(0+offset) = 1;
    for j = 1 : m
      b(j+offset) = 0;
      for i = 0 : j - 1
        b(j+offset) = b(j+offset) + i4_choose ( j - 1, i ) * b(i+offset);
      end
    end

    npart = b(m+offset);

  end

  return
end


text:
%krp_randptlist 'randomize point list'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros rp_randptlist.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% OutputFile: o 'Output', required: 'Resulting output data object'
% Integer: seed 'seed', default: 12345: 'seed for random number generator'
%
% Example: o = krp_randptlist(i, {'i','';'o','';'seed',12345})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% rp_randptlist - randomize point list
%
%  DESCRIPTION
% This program randomizes the order of the points listed in the
% "kPointList" attribute.
%
%
%
%  EXAMPLES
% See the html tutorial in $REGISTER/examples/html/README_FIRST.html
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996 - 1997, University of New Mexico.  All rights reserved.
%

function varargout = krp_randptlist(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = krp_randptlist(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output';'seed', 12345};
maxval={0,0,2};
minval={0,0,2};
istoggle=[0,0,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Integer'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'rp_randptlist"  '],Inputs,narglist);


text:
function spm_mvb_cvk_display(MVB)
% model display for MVB with cross-validation
% FORMAT spm_mvb_cvk_display(MVB)
% MVB  - multivariate Bayes structure, select one if not provided
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Christophe Phillips
% $Id: spm_mvb_cvk_display.m 3806 2010-04-06 14:42:32Z ged $

if nargin<1
    load(spm_select(1,'^MVB.*\.mat','Select MVB to display'))
end
if ~isfield(MVB,'cvk')
    error(['No crossvalidation data available. ' ...
        'Select another file or perform the crossvalidation']);
end

%-Get figure handles and set title
%--------------------------------------------------------------------------
Fmvb = spm_figure('GetWin','MVB');
spm_clf(Fmvb);

% get stuff in place for display
%--------------------------------------------------------------------------
K     = MVB.K;
X     = K*MVB.X;
X0    = orth(K*MVB.X0);
R     = speye(length(X)) - X0*X0';
R     = orth(R);
pX     = R*R'*X;

% plot validation
%--------------------------------------------------------------------------
subplot(2,2,1)
s      = 1:length(pX);
plot(s,pX,s,MVB.cvk.qX,'-.')
xlabel('sample')
ylabel('response (adjusted)')
title('cross-validation')
axis square

subplot(2,2,2)
plot(pX,MVB.cvk.qX,'.')
xlabel('true')
ylabel('predicted')
title(sprintf('p-value (parametric) = %.5f',MVB.p_value))
axis square
abc = axis;
hold on
plot([max(abc([1 3])) min(abc([2 4]))],[max(abc([1 3])) min(abc([2 4]))],'k')

% plot feature weights
%--------------------------------------------------------------------------
subplot(2,2,3)
imagesc(corrcoef(MVB.cvk.qE))
colorbar
caxis([0 1])
xlabel('bipartition (k)')
title({'correlations among';'k-fold feature weights'})
axis square

subplot(2,2,4)
if ~isempty(MVB.XYZ) && ~isempty(MVB.VOX)
    if isfield(MVB.cvk, 'P')
        spm_mip(prod(MVB.cvk.P,2),MVB.XYZ(1:3,:),MVB.VOX)
        title({[MVB.name ' (' MVB.contrast ')'];'prod( P(|weights| > 0) )'})
    else % reproduce plot from original spm_mvb_cvk2
        qe = mean(MVB.cvk.qE,2);
        qe = qe.*(qe > 0);
        spm_mip(qe,MVB.XYZ(1:3,:),MVB.VOX)
        title({[MVB.name ' (' MVB.contrast ')'];'mean (positive) weights'})
    end
else
    % (Allows MVB for non-spatial data, excluding smooth or compact priors)
    title('No spatial info in MVB')
end
axis square

fprintf('\np-value = %.4f; classification: %.1f%%; R-squared %.1f%%\n', ...
    MVB.p_value,MVB.percent,MVB.R2)


text:
function reduction_example ( n )

%*****************************************************************************80
%
%% REDUCTION_EXAMPLE shows an example of a parallel reduction operation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the size of the array.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'REDUCTION_EXAMPLE:\n' );
  fprintf ( 1, '  The MATLAB PARFOR command can carry out certain\n' );
  fprintf ( 1, '  "reduction" operations in parallel.\n' );

  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    n = input ( '  Enter N, the size of the vector.\n' );
  end
%
%  Set up a random data vector.
%
  x = rand ( n );
%
%  Sum its entries.
%

  total1 = 0.0;
  for i = 1 : n
    total1 = total1 + x(i);
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sum all entries of X, in sequential mode:\n' );
  fprintf ( 1, '  TOTAL1 = %f\n', total1 );
%
%  Sum its entries in parallel.
%
  total2 = 0.0;
  parfor i = 1 : n
    total2 = total2 + x(i);
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sum all entries of X, in parallel mode:\n' );
  fprintf ( 1, '  TOTAL2 = %f\n', total2 );
%
%  Sum some of its entries.
%
  total3 = 0.0;
  for i = 1 : n
    if ( 0.5 < x(i) )
      total3 = total3 + x(i);
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sum entries of X greater than 1/2, in sequential mode:\n' );
  fprintf ( 1, '  TOTAL3 = %f\n', total3 );
%
%  Sum some of its entries in parallel.
%
  total4 = 0.0;
  parfor i = 1 : n
    if ( 0.5 < x(i) )
      total4 = total4 + x(i);
    end
  end

  fprintf ( 1, '  Sum entries of X greater than 1/2, in parallel mode:\n' );
  fprintf ( 1, '  TOTAL4 = %f\n', total4 );

  return
end


text:
function H = lbfgsadd(H,k,dscale)
%lbfgsadd  adds k variables to the L-BFGS matrix H.
%
%   H = lbfgsadd(H,k,dscale) expands H to accomodate k new variables.
%   The new diagonal part of H (H0) gets set to dscale.
%
%   See also lbfgsdel, lbfgshprod, lbfgsupdate, lbfgsinit.

% $Id: lbfgsadd.m 385 2008-09-02 22:22:41Z mpf $

jMax = H.jMax;

H.S  = [ H.S ;    zeros(k,jMax)   ];
H.Y  = [ H.Y ;    zeros(k,jMax)   ];
H.H0 = [ H.H0; repmat(dscale,k,1) ];

author: mpf


text:
function dx = expmall(J,f,t,EP)

EP(EP==1) = J;
EP(EP==2) = f;

dx = expm(EP*t);


text:
function val = mrViewGet(varargin);
%
% mrViewGet: Get the value of the specified property for
% a mrViewer UI.
%
% Can be used a number of ways:
%
% mrViewGet without any arguments returns the ui struct
% corresponding to the most current mrViewer UI. The most
% current UI is the one the selected figure belongs to;
% if the selected figure doesn't belong to a mrViewer UI,
% then it looks for the most recently-created UI. If no
% UIs are found, displays a warning and returns an empty
% matrix.
%
% mrViewGet(ui) lists the properties and set values for
% the specified ui, similar to GET.
%
% mrViewGet('PropertyName') returns the value of the
% specified property for the most current UI.
%
% mrViewGet(ui,'PropertyName') returns the value of the
% specified property for the specified UI.
%
% Properties include (non-case-sensitive):
%   'ImageCoords': Coordinates of the images being displayed
%   by the UI (in the ui.display struct, see mrViewUnderlay), in terms
%   of the UI's currently-selected space. Specifying
%       mrViewGet(ui,'imagecoords',N)
%   will return only the coordinates pertaining to the Nth image being
%   displayed. Otherwise, a cell array of image coordinates will be
%   returned. In all cases, the coords for each image are a 3xnPixels
%   matrix, with the rows being the row, column, and slice indices into
%   the base MR data.
%
%   'DataCoords': Coordinates of the images being displayed
%   by the UI (in the ui.display struct, see mrViewUnderlay), in terms
%   of the base MR data. Specifying  mrViewGet(ui,'datacoords',N) will
%   return only the coordinates pertaining to the Nth image being
%   displayed. Otherwise, a cell array of image coordinates will be
%   returned. In all cases, the coords for each image are a 3xnPixels
%   matrix, with the rows being the row, column, and slice indices into
%   the base MR data.
%
%   'ROI', [r]: Returns ROI struct for the selected ROI. If the index r is
%   omitted, returns the selected ROI.
%
%   'MapVals', [m], [roi]: Get values of the data map m in the selected
%   ROI. m is an index into the loaded data maps. ROI can be the name
%   of a loaded ROI, index into the ui.rois field, or 3xN list of
%   coordinates in the base mr space. If omitted, returns the map value
%   at the cursor position. The returned values are of size 1 x nCoords,
%   where nCoords is the number of columns in roi.coords. For coordinates
%   in the ROI for which data are not defined, returns a NaN.
%
%   'OverlayVals', [roi]: Get values of the data map currently being
%   displayed as an overlay. If more than one overlays are displayed, the
%   returned value is a matrix of size nOverlays x nCoords. If omitted,
%   uses the cursor position as the ROI.
%
%	'OverlayMask', [whichOverlays]: Returns a 3-D logical matrix, the size of
%	ui.mr.data, with a 1 where a given voxel passes the overlay thresholds
%	and a 0 otherwise. If the optional 'whichOverlays' argument is passed,
%	will only threshold according to the selected overlays; otherwise,
%	thresholds according to all overlays.
%
%	'CursorPosition': position of the mrViewer cursor, in the "base" space
%	(units of (rows, cols, slices) in mr.ui.data).
%
% ras, 07/05.
val = [];

if nargin==0
    % find most current UI:
    % first check current figure to see if it has a UI
    tag = get(gcf,'Tag');
    if strncmp(tag,'mrViewer',8)==1
        val = get(gcf,'UserData');
    else
        ud = get(gcf,'UserData');
        if isstr(ud) & strncmp(ud,'mrViewer',8)==1
            % The user data is pointing to the main UI figure
            fig = findobj('Tag',ud);
            val= get(fig,'UserData');
        else
            % last thing to try: search for figures w/
            % tag 'mrViewer[#]', and get highest #
            tag = 'mrViewer1';
            if isempty(findobj('Tag',tag))
                warning('No mrViewer UIs found!')
                val = [];
            else
                i = 1;
                while ~isempty(findobj('Tag',sprintf('mrViewer%i',i+1)))
                    i = i + 1;
                end
                fig = findobj('Tag',sprintf('mrViewer%i',i));
                val = get(fig,'UserData');
            end
        end
    end

    return
end

if length(varargin)==1
    if isstruct(varargin{1})
        % request to print out UI info
        disp('mrViewer UI Info: ')
        ui = varargin{1};
        fprintf('Base MR file: %s\n', ui.mr.name);
        fprintf('Base MR path: %s\n\n', ui.mr.path);
        for i = 1:length(ui.maps)
            fprintf('Map %i file: %s\n', i, ui.maps(i).name);
            fprintf('Map %i path: %s\n\n', i, ui.maps(i).path);
        end
        disp('UI Settings:')
        disp(ui.settings)
        val = ui.settings;
    else
        ui = mrViewGet;
        val = mrViewGet(ui,varargin{1});
    end
    return
end

% if we've gotten here, we have both a UI struct and an argument to get
ui = varargin{1};
prop = varargin{2};
if length(varargin)>2
    args = varargin(3:end);
else
    args = {};
end

% allow an empty UI struct to be passed, meaning get the current one
if isempty(ui), ui = mrViewGet; end
if ishandle(ui), ui = get(ui, 'UserData'); end

switch lower(prop)
case {'imgcoords','imagecoords','imagecoordinates'},
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % coordinates of displayed images in selected data space %
    % (args can specify an image #, otherwise return cell of %
    % coords for all displayed images)                       %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isempty(args)       % get coords for specified image
        selectedImages = args{1};
    else                    % get coords for all images
        selectedImages = 1:length(ui.display.images);
    end

    b = ui.settings.bounds;
    val = {};
    for i = selectedImages
        % get the coords corresponding to this slice
        rng = {[b(1,1):b(1,2)] [b(2,1):b(2,2)] [b(3,1):b(3,2)]};
        rng{ui.display.oris(i)} = ui.display.slices(i);
        [X Y Z] = meshgrid(rng{2},rng{1},rng{3});
        val{end+1} = [Y(:) X(:) Z(:)]'; clear X Y Z
    end

    % if an image was specified, return the coords directly
    if ~isempty(args),  val = val{1};  end

case {'datacoords'},
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % coordinates of displayed images in base data space     %
    % (args can specify an image #, otherwise return cell of %
    % coords for all displayed images)                       %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~isempty(args)       % get coords for specified image
        selectedImages = args{1};
    else                    % get coords for all images
        selectedImages = 1:length(ui.display.images);
    end

    val = {};

    for i = selectedImages
        coords = mrViewGet(ui,'ImageCoords',i);
        if ui.settings.space ~= 1
            %%%%%%%%%%%%%%%%%%%%%
            % Transformed Space %
            %%%%%%%%%%%%%%%%%%%%%
            % xform from image space -> data space
            % (this is the inverse of the xform we
            % store in spaces)
            s = ui.settings.space;
            if ~isempty(ui.spaces(s).coords)
                % do the direct-coord method, or
                % use an external function
            else
                xform = inv(ui.spaces(s).xform);
                coords = xform * [coords; ones(1,size(coords,2))];
                coords = coords(1:3,:);
            end
        end

        % return coords
        val{end+1} = coords;
    end

    % if an image was specified, return the coords directly
    if ~isempty(args),  val = val{1};  end

case {'cursorcoords' 'cursorpos' 'cursorposition' 'cursordatacoords'}
    % position of cursor in base space
    val = ui.settings.cursorLoc;

    if ui.settings.space ~= 1
        %%%%%%%%%%%%%%%%%%%%%
        % Transformed Space %
        %%%%%%%%%%%%%%%%%%%%%
        % xform from image space -> data space
        % (this is the inverse of the xform we
        % store in spaces)
        s = ui.settings.space;
        if ~isempty(ui.spaces(s).coords)
            % do the direct-coord method, or
            % use an external function
        else
            xform = inv(ui.spaces(s).xform);
            val = xform * [val(:); 1];
            val = val(1:3,:)';
        end
	end

case {'displayformat'},		val = ui.settings.displayFormat;
case {'curslice', 'slice', 'selectedslice', 'firstslice'}
    val = ui.settings.slice;

case {'numslices' 'nslices'}
    val = ui.mr.dims(3);

case {'base' 'mr' 'anat' 'underlay'},   val = ui. mr;
case {'maps'}, val = ui.maps;
case {'overlay' 'curmap' 'selectedmap'},
	if ~isempty(args), o = args{1}; else, o = 1; end
    val = ui.maps(ui.overlays(o).mapNum);

case {'overlaysettings' 'overlaystruct'}
	if ~isempty(args), o = args{1}; else, o = 1:length(ui.overlays); end
	val = ui.overlays(o);

case {'roiviewmode' 'roidrawmode'}, val = ui.settings.roiViewMode;

case {'curroinum' 'roinum' 'selroinum' 'selectedroinum'},
	val = ui.settings.roi;
case {'roi' 'selroi' 'curroi' 'selectedroi'}
    if isempty(args), r = ui.settings.roi; else, r = args{1}; end
    val = ui.rois(r);

case {'mesh' 'curmesh'}
    s = ui.settings.segmentation;
    if s==0 | isempty(ui.segmentation), val = []; return; end

    if isempty(ui.segmentation(s).mesh)     % load a mesh
        ui = mrViewLoad(ui, '', 'mesh');
    end

    val = segGet(ui.segmentation(s), 'mesh');

case {'zoom'}, val = ui.settings.zoom;
case {'imagezoom'}, % zoom in format [xmin xmax ymin ymax] for displayed iamges
    val = ui.settings.zoom;
    if ui.settings.displayFormat==2 % 3-axis view: zoom different for each
        warning('Zoom is 3-D in 3-axis view -- returning 3 x 2 matrix');
        return
    end
    otherDims = setdiff(1:3, ui.settings.ori);
    val = [val(otherDims(2),:) val(otherDims(1),:)];

case {'displayimages' 'images'}, val = ui.display.images;
case {'croppedimages', 'zoomedimages'},  % images that are zoomed in on
     images = ui.display.images;
     zoom = round( mrViewGet(ui, 'imagezoom') );
     if ui.settings.displayFormat==2        % different zoom for each image
         for i = 1:3
             otherDims = setdiff(1:3, i); % orientation
             rows = zoom(otherDims(2),:);
             cols = zoom(otherDims(1),:);
             val{i} = images{i}(rows,cols,:);
         end
     else                                   % same zoom for each image
         for i = 1:length(images)
             val{i} = images{i}(zoom(3):zoom(4),zoom(1):zoom(2),:);
         end
     end

case {'colorbar' 'cbar'}
    if ~isempty(args), o = args{1}; else, o = 1; end
    val = ui.overlays(o).cbar;

case {'segmentationnum' 'segnum' 'cursegnum' 'cursegmentationnum'}
    val = ui.settings.segmentation;

case {'segmentation' 'seg' 'segstruct' 'cursegmentation' 'selectedsegmentation'}
    val = ui.segmentation( ui.settings.segmentation );

case {'curmeshnum' 'curmeshn' 'meshnum' 'selectedmeshnum'}
    s = ui.settings.segmentation;
    val = ui.segmentation(s).settings.mesh;

case {'curmesh' 'currentmesh' 'selectedmesh' 'mesh'}
    s = ui.settings.segmentation;
    m = ui.segmentation(s).settings.mesh;
    val = ui.segmentation(s).mesh{m};

case {'allmeshnames' 'meshlist'}
	val = {};
	for s = 1:length(ui.segmentation)
		for m = 1:length(ui.segmentation(s).mesh)
			val{end+1} = sprintf('%i %s', ...
								ui.segmentation(s).mesh{m}.id, ...
								ui.segmentation(s).mesh{m}.name);
		end
	end

case {'curmeshdir' 'selectedmeshdir' 'meshdir'}
    s = ui.settings.segmentation;
    val = ui.segmentation(s).params.meshDir;

case {'mapvals' 'mapdata'}
    if length(args) < 1, error('Need to specify a map index.'); end
    m = args{1};
    if length(args) >= 2,       roi = roiParse(args{2}, ui);
    else,                       roi = roiParse(ui.settings.cursorLoc');
	end
	% get xform from viewer coordinates -> map coordinates
	viewer2MapXform = inv(ui.maps(m).baseXform);

	% xform ROI coords to map coords
	roi = roiXformCoords(roi, viewer2MapXform, ui.maps(m).voxelSize(1:3));

	% sample the data at the map coordinates
	val = mrGet(ui.maps(m), 'data', roi.coords);

case {'overlayvals' 'overlaydata'}
    if isempty(args), roi = roiParse(ui.settings.cursorLoc');
    else,             roi = roiParse(args{1}, ui);
    end

    for o = 1:length(ui.overlays)
        val(o,:) = mrViewGet(ui, 'mapvals', ui.overlays(o).mapNum, roi);
	end

case {'overlaymask' 'datamask'}
    %% 3-D binary matrix with 1 where voxels pass the overlay thresholds,
    %% and 0 otherwise.
    if isempty(args), whichOverlays = find( ~([ui.overlays.hide]) );
    else,             whichOverlays = args{1};
	end

	% initialize empty matrix
	val = logical( zeros(ui.mr.dims(1:3)) );

	% input check
	if isempty(whichOverlays)
		warning('[mrViewGet(''overlaymask'')]: No overlays. Returning emtpy matrix.')
		return
	end

	% get the coords for all slices in the base data space
	[X Y Z] = meshgrid( 1:size(val,2), 1:size(val,1), 1:size(val,3) );
	coords = [Y(:) X(:) Z(:)]';

	% restrict coords to those which pass the threshold
	% of this overlay
	[coords, imageMask] = mrViewRestrict(ui, coords, whichOverlays);

	% set those coords which pass the threshold to 1 (true)
	I = sub2ind( size(val), coords(1,:), coords(2,:), coords(3,:) );
	val(I) = true;

otherwise, warning('Unknown mrViewer property.')

end

return


text:
function sparse_test04 ( )

%*****************************************************************************80
%
%% SPARSE_TEST04 demonstrates incrementing a sparse matrix.
%
%  Discussion:
%
%    For this example, we show that you don't need to have the entire
%    sparse matrix data available in order to start defining the matrix.
%    (Of course, it is more efficient to make fewer calls, with more data,
%    but it is not the only way to work.)
%
%    For simplicity, we assume the matrix is 10 x 10, and we
%    generate 50 random coordinate pairs (I,J), and increment them by 1.
%    Initially, the sparse matrix is empty, but each time we call,
%    we are either creating a new entry, or incrementing an existing one.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Lindfield, John Penny,
%    Numerical Methods Using MATLAB,
%    Prentice Hall, 1999
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_TEST04:\n' );
  fprintf ( 1, '  Demonstrate the use of MATLAB''s SPARSE facility\n' );
  fprintf ( 1, '  to set or increment a matrix one entry at a time.\n' );
%
%  This defines an "empty" sparse matrix.
%
  i = [];
  j = [];
  v = [];
  m = 10;
  n = 10;

  a = sparse ( i, j, v, m, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I     J  New A(I,J)\n' );
  fprintf ( 1, '\n' );

  for test = 1 : 50

    i = round ( 10.0 * rand ( ) + 0.5 );
    j = round ( 10.0 * rand ( ) + 0.5 );
    a(i,j) = a(i,j) + 1;
%
%  FPRINTF won't print a sparse matrix reference,
%  so you have to convert it to FULL first.
%
    fprintf ( 1, '  %4d  %4d  %g\n', i, j, full ( a(i,j) ) );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nonzero entries is %d\n', nnz ( a ) );
%
%  FPRINT won't print the SUM SUM of A,
%  so you have to convert it to FULL first.
%
  fprintf ( 1, '  Sum of entries is %f\n', full ( sum ( sum ( a ) ) ) );
%
%  Typing "a" will list the sparse triplet form of the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix printed in sparse triplet form:\n' );
  fprintf ( 1, '\n' );

  a
%
%  To see the usual layout, you must convert it to FULL.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix printed in usual full matrix form:\n' );
  fprintf ( 1, '\n' );

  full ( a )

  return
end


text:
function anal = er_wtaClassifier(amps1,amps2,subtractDC,plotFlag,names);
%
% anal = er_wtaClassifier(amps1,amps2,[plotFlag],[names]);
%
% Apply a winner-take-all classifier to event-related
% data.
%
% [more description forthcoming]
%
% Input args: amps1 and amps2 should be of shape voxels x conditions:
% each column represents the pattern of response across voxels to a
% given condition. (See mv_reliability for an example of how
% to compute these.)
%
% plotFlag: 1 x 3 vector to flag each of the following
% possible output plots:
%   1) image the amplitudes for each subset
%   2) summarize the reliability analysis
%   3) 'omnibus' contrast: regress mean response, across conditions,
%   from each subset
%
% subtractDC: option to subtract the baseline response
% across conditions for
%
% ras 05/05. Based on exp5_haxby, by ras 03/05.
if ieNotDefined('subtractDc')
    subtractDC = 0;
else
    subtractDC
end

if ieNotDefined('plotFlag')
    plotFlag = [0 0 0];
else
    plotFlag
end

if ieNotDefined('names')
    for i = 1:size(amps1,2)
        names{i} = num2str(i);
    end
end

if size(amps1) ~= size(amps2)
    error('amps1 and amps2 must be same size.')
end

nVoxels = size(amps1,1);
nConds = size(amps1,2);
font = 'Arial';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% omnibus test                    %
% (always do this before zeroing, %
% or you remove the effect)       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu1 = mean(amps1,2);
mu2 = mean(amps2,2);
% ignore zeroed voxels
notzeroed = find(mu1~=0 & mu2~=0);
mu1 = mu1(notzeroed);
mu2 = mu2(notzeroed);
[R P] = corrcoef(mu1(:),mu2(:));
anal.omnibusR = R(1,2);
anal.omnibusP = P(1,2);
anal.mu1 = mu1;
anal.mu2 = mu2;

if subtractDC==1
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Normalize responses, in a manner similar   %
	% to Haxby et al 2001 -- subtract mean       %
	% response across conditions from each voxel %
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	bsl1 = repmat(mean(amps1,2),1,size(amps1,2));
	usedAmps1 = amps1 - bsl1;
	bsl2 = repmat(mean(amps2,2),1,size(amps2,2));
	usedAmps2 = amps2 - bsl2;
else
    % use the non-corrected amps
    usedAmps1 = amps1;
    usedAmps2 = amps2;
end

% assign amplitude fields to anal struct
anal.amps1 = amps1;
anal.amps2 = amps2;
anal.usedAmps1 = usedAmps1;
anal.usedAmps2 = usedAmps2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform pairwise correlations between maps %
% in the first and second data sets          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:nConds
    for j = 1:nConds
        training = usedAmps1(:,i);
        test = usedAmps2(:,j);

        % remove NaNs
        ok = ~isnan(training) & ~isnan(test);
        if sum(ok)==0
            error('Too many NaNs in amplitudes.');
        end
        training = training(ok); test = test(ok);

        [R p] = corrcoef(training,test);
        anal.corrRvals(i,j) = R(1,2);
        anal.corrPvals(i,j) = p(1,2);
    end
end

% % 02/22/05: kalanit suggested thresholding
% % the R vals, based on the p -- let's try it:
% anal.corrRvals(anal.corrPvals>0.05) = 0;

% as an ideal observer of activity patterns,
% 'guess' the image shown in the 2nd data set
% by looking at what provoked the most similar
% (by correlation) response from the 1st data set
for i = 1:nConds
    anal.guess(i,:) = zeros(1,nConds);
    guess = find(anal.corrRvals(i,:)==max(anal.corrRvals(i,:)));
    anal.guess(i,guess) = 1;

    % guesses along the diagonal are correct -- the
    % highest-correlated pattern comes from the same image.
    % other guesses are mistakes.
    if isempty(guess), guess = 0; end
    anal.correct(i) = (guess(1)==i);
end
anal.pctCorrect = 100 * sum(anal.correct)/length(anal.correct);

% let's also compute the 'voxel reliability',the
% correlation value for each voxel of the set of nConds
% response amplitudes between the two subsets:
for v = 1:nVoxels
    [R p] = corrcoef(anal.amps1(v,:),anal.amps2(v,:));
    voxR(v) = R(1,2);
end
anal.voxR = voxR;

% calculate mean correlations for same image,
% diff image but same cat, diff cat:
nImgs = nConds; %length(1:nConds);
group = 3*ones(nImgs,nImgs); % diff cat
for i = 1:nImgs
    for j = 1:nImgs
        if ceil(j/4)==ceil(i/4)
            group(i,j) = 2; % same cat
        end
    end
    group(i,i) = 1;  % same img
end

anal.group = group;
anal.meanR = [];
anal.semR = [];
for j = 1:3
    anal.meanR(j) = mean(anal.corrRvals(group==j));
    nObs = sum(group(:)==j);
    anal.semR(j) = std(anal.corrRvals(group==j))/sqrt(nObs-1);
    try
        [H p] = ttest2(anal.corrRvals(group==j),anal.corrRvals(group==3));
    catch
        H = 0; p = 1;
    end
    anal.sigDiffR(j) = H;
    anal.pvalR(j) = p;
end

%%%%%%%%%%%%%%%%%%%%%
% plot the results  %
%%%%%%%%%%%%%%%%%%%%%
% Amplitude plots
if plotFlag(1)==1
	anal.h1 = figure('Name','Amplitudes From Each Half of Data',...
                    'Units','Normalized','Position',[.7 .8 .3 .18],...
                    'MenuBar','none',...
                    'Color','w'); % [0 .4 .7 .4]
	colormap jet

    minVal = min([min(anal.amps1(:)) min(anal.amps2(:))]);
	maxVal = max([max(anal.amps1(:)) max(anal.amps2(:))]);

    subplot('Position',[.1 .15 .35 .75])
	imagesc(anal.amps1,[minVal maxVal]);
	set(gca,'XTick',1:nConds,'XTickLabel',names);
	ylabel('Voxels');
    title('Subset 1 (Training)')
    set(gca,'FontName',font,'FontSize',10);

    subplot('Position',[.5 .15 .35 .75])
	imagesc(anal.amps2,[minVal maxVal]);
	set(gca,'XTick',1:nConds,'XTickLabel',names,'YTick',[]);
	ylabel('Voxels');
    title('Subset 2 (Test)');
    set(gca,'FontName',font,'FontSize',10);

    hcb = subplot('Position',[.9 .2 .03 .6]);
    colorbar(hcb);
    ylabel('Response Amplitude, % Signal');
    set(gca,'FontName',font,'FontSize',10);
end

% reliability summary
if plotFlag(2)==1
	anal.h2 = figure('Name','Reliability Amps Results',...
                'Color','w','MenuBar','none','Units','Normalized',...
                'Position',[.7 .6 .3 .18]); % [.1 .2 .7 .7]
	subplot(2,2,1);
	imagesc(anal.corrRvals); % colorbar;
	set(gca,'XTick',1:nConds,'XTickLabel',names,...
            'YTick',1:nConds,'YTickLabel',names);
	set(gca,'FontName',font,'FontSize',10);
	xlabel('First Half');
	ylabel('Second Half');
	title('Correlation Coefficient R');

    subplot(2,2,2);
	imagesc(anal.guess); % colorbar;
	set(gca,'XTick',1:nConds,'XTickLabel',names,...
            'YTick',1:nConds,'YTickLabel',names);
	set(gca,'FontName',font,'FontSize',10);
	xlabel('First Half');
	ylabel('Second Half');
	title('Ideal Observer Best Guess');

    subplot(2,2,3);
	starbar(anal.meanR,anal.semR,anal.sigDiffR);
	grps = {'Same Image' 'Same Category' 'Different Category'};
	set(gca,'XTickLabel',grps);
	set(gca,'FontName',font,'FontSize',10);
	ylabel('Mean R Value');

    subplot(2,2,4);
	table = {'Percent Correct:'; num2str(anal.pctCorrect)};
	plotTable(table);
end

% 'omnibus' test
if plotFlag(3)==1
	anal.h3 = figure('Color','w','MenuBar','none','Units','Normalized',...
                     'Position',[.7 .4 .2 .18]);
	regressPlot(mu1,mu2);
	set(gca,'FontName',font,'FontSize',10);
	axis equal
	axis square
	xlabel('Amplitude, Odd Runs','FontName','Helvetica','FontSize',nConds);
	ylabel('Amplitude, Even Runs','FontName','Helvetica','FontSize',nConds);
	ttltxt = sprintf('Amplitudes Across All Images');
	title(ttltxt,'FontName','Helvetica','FontSize',14);
end

return


text:
function R1f = computeR1(Param, R1obs)

kf  = Param.kf;
F  = Param.F;
R1r = Param.R1r;

R1f = R1obs - kf*(R1r - R1obs) / (R1r - R1obs + kf/F);

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab