text:
% This function belongs to Piotr Dollar's Toolbox
% http://vision.ucsd.edu/~pdollar/toolbox/doc/index.html
% Please refer to the above web page for definitions and clarifications
%
% Calculates the distance between sets of vectors.
%
% Let X be an m-by-p matrix representing m points in p-dimensional space
% and Y be an n-by-p matrix representing another set of points in the same
% space. This function computes the m-by-n distance matrix D where D(i,j)
% is the distance between X(i,:) and Y(j,:).  This function has been
% optimized where possible, with most of the distance computations
% requiring few or no loops.
%
% The metric can be one of the following:
%
% 'euclidean' / 'sqeuclidean':
%   Euclidean / SQUARED Euclidean distance.  Note that 'sqeuclidean'
%   is significantly faster.
%
% 'chisq'
%   The chi-squared distance between two vectors is defined as:
%    d(x,y) = sum( (xi-yi)^2 / (xi+yi) ) / 2;
%   The chi-squared distance is useful when comparing histograms.
%
% 'cosine'
%   Distance is defined as the cosine of the angle between two vectors.
%
% 'emd'
%   Earth Mover's Distance (EMD) between positive vectors (histograms).
%   Note for 1D, with all histograms having equal weight, there is a simple
%   closed form for the calculation of the EMD.  The EMD between histograms
%   x and y is given by the sum(abs(cdf(x)-cdf(y))), where cdf is the
%   cumulative distribution function (computed simply by cumsum).
%
% 'L1'
%   The L1 distance between two vectors is defined as:  sum(abs(x-y));
%
%
% USAGE
%  D = pdist2( X, Y, [metric] )
%
% INPUTS
%  X        - [m x p] matrix of m p-dimensional vectors
%  Y        - [n x p] matrix of n p-dimensional vectors
%  metric   - ['sqeuclidean'], 'chisq', 'cosine', 'emd', 'euclidean', 'L1'
%
% OUTPUTS
%  D        - [m x n] distance matrix
%
% EXAMPLE
%  [X,IDX] = demoGenData(100,0,5,4,10,2,0);
%  D = pdist2( X, X, 'sqeuclidean' );
%  distMatrixShow( D, IDX );
%
% See also PDIST, DISTMATRIXSHOW

% Piotr's Image&Video Toolbox      Version 2.0
% Copyright (C) 2007 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

function D = pdist2(X, Y, metric)

if( nargin<3 || isempty(metric) ); metric=0; end;

switch metric
  case 'sqeuclidean'
    D = distEucSq( X, Y );
  case {0, 'euclidean'}
    D = sqrt(distEucSq( X, Y ));
  case 'L1'
    D = distL1( X, Y );
  case 'cosine'
    D = distCosine( X, Y );
  case 'emd'
    D = distEmd( X, Y );
  case 'chisq'
    D = distChiSq( X, Y );
  otherwise
    error(['pdist2 - unknown metric: ' metric]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = distL1( X, Y )

m = size(X,1);  n = size(Y,1);
mOnes = ones(1,m); D = zeros(m,n);
for i=1:n
  yi = Y(i,:);  yi = yi( mOnes, : );
  D(:,i) = sum( abs( X-yi),2 );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = distCosine( X, Y )

if( ~isa(X,'double') || ~isa(Y,'double'))
  error( 'Inputs must be of type double'); end;

p=size(X,2);
XX = sqrt(sum(X.*X,2)); X = X ./ XX(:,ones(1,p));
YY = sqrt(sum(Y.*Y,2)); Y = Y ./ YY(:,ones(1,p));
D = 1 - X*Y';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = distEmd( X, Y )

Xcdf = cumsum(X,2);
Ycdf = cumsum(Y,2);

m = size(X,1);  n = size(Y,1);
mOnes = ones(1,m); D = zeros(m,n);
for i=1:n
  ycdf = Ycdf(i,:);
  ycdfRep = ycdf( mOnes, : );
  D(:,i) = sum(abs(Xcdf - ycdfRep),2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = distChiSq( X, Y )

%%% supposedly it's possible to implement this without a loop!
m = size(X,1);  n = size(Y,1);
mOnes = ones(1,m); D = zeros(m,n);
for i=1:n
  yi = Y(i,:);  yiRep = yi( mOnes, : );
  s = yiRep + X;    d = yiRep - X;
  D(:,i) = sum( d.^2 ./ (s+eps), 2 );
end
D = D/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = distEucSq( X, Y )

%if( ~isa(X,'double') || ~isa(Y,'double'))
 % error( 'Inputs must be of type double'); end;
m = size(X,1); n = size(Y,1);
%Yt = Y';
XX = sum(X.*X,2);
YY = sum(Y'.*Y',1);
D = XX(:,ones(1,n)) + YY(ones(1,m),:) - 2*X*Y';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function D = distEucSq( X, Y )
%%%% code from Charles Elkan with variables renamed
% m = size(X,1); n = size(Y,1);
% D = sum(X.^2, 2) * ones(1,n) + ones(m,1) * sum(Y.^2, 2)' - 2.*X*Y';

%%% LOOP METHOD - SLOW
% [m p] = size(X);
% [n p] = size(Y);
%
% D = zeros(m,n);
% onesM = ones(m,1);
% for i=1:n
%   y = Y(i,:);
%   d = X - y(onesM,:);
%   D(:,i) = sum( d.*d, 2 );
% end

%%% PARALLEL METHOD THAT IS SUPER SLOW (slower then loop)!
% % From "MATLAB array manipulation tips and tricks" by Peter J. Acklam
% Xb = permute(X, [1 3 2]);
% Yb = permute(Y, [3 1 2]);
% D = sum( (Xb(:,ones(1,n),:) - Yb(ones(1,m),:,:)).^2, 3);

%%% USELESS FOR EVEN VERY LARGE ARRAYS X=16000x1000!! and Y=100x1000
% call recursively to save memory
% if( (m+n)*p > 10^5 && (m>1 || n>1))
%   if( m>n )
%     X1 = X(1:floor(end/2),:);
%     X2 = X((floor(end/2)+1):end,:);
%     D1 = distEucSq( X1, Y );
%     D2 = distEucSq( X2, Y );
%     D = cat( 1, D1, D2 );
%   else
%     Y1 = Y(1:floor(end/2),:);
%     Y2 = Y((floor(end/2)+1):end,:);
%     D1 = distEucSq( X, Y1 );
%     D2 = distEucSq( X, Y2 );
%     D = cat( 2, D1, D2 );
%   end
%   return;
% end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function s = yjf_to_s_hebrew ( y, j, f )

%*****************************************************************************80
%
%% YJF_TO_S_HEBREW writes a Hebrew YJF date into a string.
%
%  Format:
%
%    AM YYYY/JJJ.FF
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, J, real F, the YJF date.
%
%    Output, string S, the representation of the date.
%
  s = sprintf ( 'AM %d/%f', y, j + f );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function tname = treepart(item, dflag)

% function tname = treepart(item, dflag)
% tree part to search - for cfg_repeat/cfg_choice this is val for filled
% cfg_items and values for default cfg_items.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: treepart.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

if dflag
    tname = 'values';
else
    tname = 'val';
end;

author: spm


text:
function termimpl(term_typ,term_opr,term_loc)
%
% Utility Function: TERMIMPL
%
% The purpose of this function is to implement the add, edit, or iterate
% operations

% Author: Craig Borghesani
% Date: 8/8/94
% Revised:
% Copyright (c) 1999, Prentice-Hall

% obtain handle information
f = gcf;
ui_data = get(f,'userdata');
ui_han = ui_data{1};
ui_obj = ui_data{2};
%popup      = ui_han(1);
%txt_top    = ui_han(2);
%txt_btm    = ui_han(3);
%edt_top    = ui_han(4);
%edt_btm    = ui_han(5);
%sld_top    = ui_han(6);
%sld_btm    = ui_han(7);
%lft_btn    = ui_han(8);
%rht_btn    = ui_han(9);
popup = ui_obj.SysUI(5);
txt_top = ui_obj.SysUI(6);
txt_btm = ui_obj.SysUI(9);
edt_top = ui_obj.SysUI(7);
edt_btm = ui_obj.SysUI(10);
sld_top = ui_obj.SysUI(8);
sld_btm = ui_obj.SysUI(11);
lft_btn = ui_obj.SysUI(12);
rht_btn = ui_obj.SysUI(13);
plant      = ui_han(30);
for_cont   = ui_han(31);
bac_cont   = ui_han(32);
set_to_one1= ui_han(33);
set_to_one2= ui_han(34);

% file modification
set(ui_han(47),'userdata',1);

cur_sys = get(ui_han(30),'userdata');
cur_page = get(ui_han(31),'userdata');
sys_state = get(ui_han(33),'userdata');
page_setup = get(ui_han(34),'userdata');

temp1_mat = [1,0,NaN,1;0,NaN,NaN,2];

% variable term_mat alternates as either the current controller matrix or
% current plant matrix
if strcmp(get(plant,'checked'),'on'),
 term_mat = get(ui_han(3+cur_sys),'userdata');
elseif strcmp(get(for_cont,'checked'),'on'),
 term_mat = get(ui_han(6+cur_sys),'userdata');
else
 term_mat = get(ui_han(9+cur_sys),'userdata');
end

% get data from ADD/UPDATE and ITERATE uicontrols
if any(term_opr == [1,2]), % Add/Update operation
 val1 = str2num(get(edt_top,'string'));
 val2 = str2num(get(edt_btm,'string'));

 if ~length(val1) & any(term_typ==[1,2,4,5,9,10]),
  katcherr(2,term_typ); return;
 elseif (~length(val1) | (~length(val2))) & any(term_typ == [6,7,8,9,11]),
  katcherr(2,term_typ); return;
 elseif any(term_typ==[1:5,11]),
  if (length(val1) & val1==0 & any(term_typ==[1:3])) | (length(val2) & val2==0 & any(term_typ==[6,7,11])),
   katcherr(4,term_typ); return;
  end
 elseif term_typ == 8,
  if abs(val1) > 88,
   katcherr(5,8); return;
  elseif val2 == 0,
   katcherr(6,8); return;
  end
 end

 if ~length(val2), val2 = NaN; end

% update top and bottom sliders
 if term_opr == 2,
  if term_typ ~= 9,
   set(sld_top,'min',0.5*val1,'max',1.5*val1,'value',val1);
   if ~isnan(val2),
    set(sld_btm,'min',0.5*val2,'max',1.5*val2,'value',val2);
   end
  end
 end

elseif term_opr == 3, % Iterate operation

 val1 = get(sld_top,'value');
 set(edt_top,'string',num2str(val1));
 min1 = get(sld_top,'min'); max1 = get(sld_top,'max');
 if abs(val1) <= abs(min1) | abs(val1) >= abs(max1),
  set(sld_top,'min',0.5*val1,'max',1.5*val1,'value',val1);
 end

 val2 = get(sld_btm,'value');
 if strcmp(get(sld_btm,'enable'),'off'),
  val2 = NaN;
 else
  set(edt_btm,'string',num2str(val2));
  min2 = get(sld_btm,'min'); max2 = get(sld_btm,'max');
  if abs(val2) <= abs(min2) | abs(val2) >= abs(max2),
   set(sld_btm,'min',0.5*val2,'max',1.5*val2,'value',val2);
  end
 end

elseif term_opr == 4, % Delete operation
 val1 = term_mat(term_loc,1);
 val2 = term_mat(term_loc,2);

end

if any(term_typ == [1,4:8]), % standard elements
 if term_typ == 1,
  new_term = [val1,term_mat(1,2),NaN,1];
 else
  if imag(val1)~=0,
   zeta = -cos(atan2(imag(-val1),real(-val1)));
   wn = abs(val1);
   term_typ = term_typ + 2;
   val1 = zeta;
   val2 = wn;
  end
  new_term = [val1,val2,NaN,term_typ];
 end

elseif any(term_typ==[2,3]), % integrators/differentiators
 new_term = [val1*((term_typ==2)-(term_typ==3)),val2,NaN,2];

elseif term_typ == 9, % num/den
 new_term = termpars(val1,val2);

elseif term_typ == 10, % delay time
 new_term = [term_mat(1,1),val1,NaN,1];

elseif term_typ == 15, % setting = 1

 set_to_one = get(f,'currentmenu');
 if set_to_one == set_to_one1,
  loc1 = cur_sys*2 - 1;
  loc2 = 4;
 else
  loc1 = cur_sys*2;
  loc2 = 5;
 end

 set_label = get(set_to_one,'label');
 if length(find(set_label=='P')), sel = 3;
 elseif length(find(set_label=='G')), sel = 6;
 else sel = 9; end

 if strcmp(get(ui_han(29+sel/3),'enable'),'off'), % retrieving
  new_term = get(ui_han(12+loc1),'userdata');
  set(ui_han(12+loc1),'userdata',temp1_mat);
  set(ui_han(sel+cur_sys),'userdata',new_term);
  set(set_to_one,'checked','off');
  term_opr = 1;
  set(ui_han(29+sel/3),'enable','on');
  sys_state(loc2,cur_sys) = 0;
 else % storing
  new_term = get(ui_han(sel+cur_sys),'userdata');
  set(ui_han(12+loc1),'userdata',new_term);
  set(ui_han(sel+cur_sys),'userdata',temp1_mat);
  set(set_to_one,'checked','on');
  term_opr = 4;
  set(ui_han(29+sel/3),'enable','off');
  sys_state(loc2,cur_sys) = 1;
 end
end
set(ui_han(33),'userdata',sys_state);

% reset root locus storage
if (term_typ ~= 1) | (term_typ == 1 & sign(val1) ~= sign(term_mat(1,1))),
 set(ui_han(37),'userdata',[]);
end

if term_opr == 1, % Add

 if any(term_typ == [1,4:8]), % standard elements
  term_mat = [term_mat;new_term];
  term_loc = length(term_mat(:,1));

 elseif any(term_typ == [2,3]), % integrator/differentiator
  term_mat(2,:) = [term_mat(2,1)+new_term(1),new_term(2:4)];
  term_loc=2;

 elseif term_typ == 9, % num/den
  term_mat(1,1) = term_mat(1,1)*new_term(1,1);
  term_mat(2,1) = term_mat(2,1)+new_term(2,1);
  term_mat = [term_mat;new_term(3:length(new_term(:,1)),:)];

 elseif term_typ == 10, % delay time
  term_mat(1,2) = new_term(1,2);

 end

elseif any(term_opr == [2,3]), % Update/Iterate

 old_term = term_mat(term_loc,:);
 term_mat(term_loc,:) = new_term;

else % Delete

 if term_typ == 2,
  term_mat(2,1)=0;
 elseif any(term_typ == [1,4:8]),
  term_mat(term_loc,:)=[];
 end

end

if strcmp(get(plant,'checked'),'on'),
 set(ui_han(3+cur_sys),'userdata',term_mat);
elseif strcmp(get(for_cont,'checked'),'on'),
 set(ui_han(6+cur_sys),'userdata',term_mat);
else
 set(ui_han(9+cur_sys),'userdata',term_mat);
end

term_res = get(ui_han(cur_sys),'userdata');
w = term_res(1,:);

new_res = termcplx(new_term,w);

if term_opr == 1, % Add
 old_res = ones(1,length(w));
% if term_typ < 9,
%  set(lft_btn,'string','Update',...
%      'callback',['termimpl(',int2str(term_typ),',2,',int2str(term_loc),')']);
% end

elseif any(term_opr == [2,3]), % Edit/Iterate
 old_res = termcplx(old_term,w);

else % delete
 old_res = new_res;
 new_res = ones(1,length(w));

end

term_res(2,:) = term_res(2,:).*(new_res./old_res);
set(ui_han(cur_sys),'userdata',term_res);

% display laplacian transfer function
if term_loc >= 0,
 dispfrac
end

pageplot(0,[]);

if term_typ < 10,
 if term_opr == 1,
%  termmang(term_typ,2,term_loc);
  set(rht_btn,'callback',['termimpl(',int2str(term_typ),',4,',int2str(term_loc),')']);
 elseif term_opr == 4,
  termmang(term_typ,1,0);
 end
end

author: Sable

repo: mcbench-benchmarks


text:
function inside = sphere_imp_contains_point_3d ( r, center, p )

%*****************************************************************************80
%
%% SPHERE_IMP_CONTAINS_POINT_3D: point in implicit sphere in 3D?
%
%  Discussion:
%
%    An implicit sphere in 3D satisfies the equation:
%
%      sum ( ( P(1:DIM_NUM) - CENTER(1:DIM_NUM) )**2 ) = R**2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real CENTER(3), the center of the sphere.
%
%    Input, real P(3), the point to be checked.
%
%    Output, logical INSIDE, is TRUE if the point is inside the sphere.
%
  dim_num = 3;

  if ( sum ( ( p(1:dim_num) - center(1:dim_num) ).^2 ) <= r * r )
    inside = 1;
  else
    inside = 0;
  end

  return
end


text:
% VOCLABELCOLORMAP Creates a label color map such that adjacent indices have different
% colors.  Useful for reading and writing index images which contain large indices,
% by encoding them as RGB images.
%
% CMAP = VOCLABELCOLORMAP(N) creates a label color map with N entries.
function cmap = labelcolormap(N)

if nargin==0
    N=256
end
cmap = zeros(N,3);
for i=1:N
    id = i-1; r=0;g=0;b=0;
    for j=0:7
        r = bitor(r, bitshift(bitget(id,1),7 - j));
        g = bitor(g, bitshift(bitget(id,2),7 - j));
        b = bitor(b, bitshift(bitget(id,3),7 - j));
        id = bitshift(id,-3);
    end
    cmap(i,1)=r; cmap(i,2)=g; cmap(i,3)=b;
end
cmap = cmap / 255;

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
function [ prediction, value ] = LRClassify( weights, xTest )
    [dataSize, ~] = size(xTest);

    % Add the column of ones for the bias feature
    xTest = [ones(dataSize, 1) xTest];

    % Find the probability of the image being in each class
    probability = sigmoid(xTest * weights');

    % Find the best class label for the image
    % The best class is one with the highest probability
    [value, index] =  max(probability, [], 2);
    prediction = index - 1;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [kHz] = radps2kHz(radps)
% Convert frequency from radians per second to kilohertz.
% Chad A. Greene 2012
kHz = radps/(2*pi)/1000;


text:
% Digital Video Stabilization and Rolling Shutter Correction using Gyroscopes
% Copyright (C) 2011 Alexandre Karpenko
%
% 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
% 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 [mov num_frames vid_height vid_width frame_rate duration] = import_video(file)

% read corresponding movie
xyloObj = mmreader(['data/' file '.mov']);
display(xyloObj);

num_frames = xyloObj.NumberOfFrames;
vid_height = xyloObj.Height;
vid_width = xyloObj.Width;
frame_rate = xyloObj.FrameRate;
duration = xyloObj.Duration;

% Preallocate movie structure.
mov(1:num_frames) = ...
    struct('cdata', zeros(vid_height, vid_width, 3, 'uint8'),...
           'colormap', []);

% Read one frame at a time.
for k = 1 : num_frames
    mov(k).cdata = read(xyloObj, k);
end
% Play back the movie once at the video's frame rate.
% movie(mov, 1, xyloObj.FrameRate);

%save(['mat/' file], 'gyro', 'frame_time', 'mov', 'num_frames', 'vid_height', 'vid_width', 'frame_rate', 'duration');

author: alex-golts

repo: Video-Stabilization


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

% DTFT of a discrete time signal x[n]

%x[n]=0.8^n , 0<=n<=20
syms w
n= 0:20;
x=0.8.^n;
X=sum(x.*exp(-j*w*n));
ezplot(abs(X),[-pi pi])
title(' Magnitude of DTFT')
ylim([0 5.4])

figure
w1=-pi:.01:pi;
XX=subs(X,w,w1)
plot(w1,angle(XX));
xlim([-pi pi])
title('Phase of DTFT')

figure
ezplot(abs(X),[-5*pi 5*pi]);
title('Magnitude of DTFT in 5 periods')
ylim([0 5.8])

figure
w1=-5*pi:.01:5*pi;
XX=subs(X,w,w1);
plot(w1,angle(XX));
xlim([-5*pi 5*pi])
title(' Phase of DTFT in 5 periods')

%x[n]=0.8^n u[n]
figure
syms n w
x=0.6^n
X=symsum(x*exp(-j*w*n),n,0,inf)
w1=-pi:.01:pi;
X_=subs(X,w,w1);
plot(w1,abs(X_));
legend('|X(\omega)|')
xlim([-pi pi])

figure
plot(w1,angle(X_));
xlim([-pi pi])
legend('\angle X(\omega)')


text:
function [ i, j ] = i4mat_min_index ( m, n, a )

%*****************************************************************************80
%
%% I4MAT_MIN_INDEX returns the location of the minimum of an I4AMT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, integer A(M,N), the M by N matrix.
%
%    Output, integer I, J, the indices of the minimum entry of A.
%
  i = -1;
  j = -1;

  for jj = 1 : n
    for ii = 1 : m
      if ( ii == 1 && jj == 1 )
        i = ii;
        j = jj;
      elseif ( a(ii,jj) < a(i,j) )
        i = ii;
        j = jj;
      end
    end
  end

  return
end


text:
function AConstraint = assignschur(AConstraint,thecompiler,varargin)
AConstraint = flatten(AConstraint);
AConstraint.clauses{1}.schurfun  = thecompiler;
AConstraint.clauses{1}.schurdata = varargin;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%This function sorts the input array in ascending order using the Comb Sort algorithm
%For details, refer https://en.wikipedia.org/wiki/Comb_sort

function y = comb_sort(array)

len = length(array);
k = len;
isSwapped = true;
% value of shrink should be greater than 1
shrink = 1.4;
while ((k > 1) || (isSwapped == true))
    k = max(floor(k / shrink),1);
    % Bubble sort with given value of k
    i = 1;
    isSwapped = false;
    while ((i + k) <= len)
        if (array(i) > array(i + k))
            array = swap(array,i,i + k);
            isSwapped = true;
        end
        i = i + 1;
    end
end
y = array;
end

function array = swap(array,i,j)
value = array(i);
array(i) = array(j);
array(j) = value;
% Note: In practice, array should be passed by reference
end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function G = approxfcn(F,range)
%APPROXFCN Approximation function.
%   G = APPROXFCN(F,RANGE) returns a function handle, G, that
%   approximates the function handle F by using a lookup table. RANGE is
%   an M-by-2 matrix specifying the input range for each of the M inputs
%   to F.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

num_inputs = size(range,1);
max_table_elements = 10000;
max_table_dim = 100;
table_dim = min(floor(max_table_elements^(1/num_inputs)), ...
   max_table_dim);

% Compute the input grid values.
inputs = cell(1,num_inputs);
grid = cell(1,num_inputs);
for k = 1:num_inputs
   grid{k} = linspace(range(k,1),range(k,2),table_dim);
end

if num_inputs > 1
   [inputs{:}] = ndgrid(grid{:});
else
   inputs = grid;
end

% Initialize the lookup table.
table = zeros(size(inputs{1}));

% Initialize the waitbar.
bar = waitbar(0,'Working...');

% Initialize the cell array used to pass inputs to F.
Zk = cell(1,num_inputs);
L = numel(inputs{1});
% Update the progress bar at 2% intervals.
waitbar_update_interval = ceil(0.02 * L);

for p = 1:L
   for k = 1:num_inputs
      Zk{k} = inputs{k}(p);
   end
   table(p) = F(Zk{:});
   if (rem(p,waitbar_update_interval) == 0)
      % Update the progress bar.
      waitbar(p/L);
   end
end
close(bar)

G = @tableLookupFcn;

   %-------------------------------------------------------------------%
   function out = tableLookupFcn(varargin)
      if num_inputs > 1
         out = interpn(grid{:},table,varargin{:});
      else
         out = interp1(grid{1},table,varargin{1});
      end
   end

end

author: dipum

repo: dipum-toolbox

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


text:
sigmoid1 = @(x) 1./(1+exp(-x));

kl_tmp=zeros(1,maxiter);
beta_tmp=zeros(1,maxiter);
dm_tmp=zeros(1,maxiter);
iter=0;
dm=1;
eta=eta0;
m=m+1e-13*randn(1,n);
m=max(m,1e-10);
m=min(m,1-1e-10);
m0 = m;
if (beta>=beta_max), beta=beta_max-1; end;
while ((dm>dmmin)&&(iter<maxiter)&&(beta<beta_max))
	iter=iter+1;
	z=(1-m)./m.*diag(C)';
	C1=C+diag(z);
    v=C1\b;
	w1=v./m';
	betaold=beta;
	beta=1/(sigmay-b'*v);

    mold=m;
    m=(1-eta)*m+eta*sigmoid1(gamma+0.5*p*beta*(w1').^2.*diag(C)');
    m=max(m,1e-10);
    m=min(m,1-1e-10);
    dm=max(abs(m-mold));
    kl_tmp(iter)=-p/2*log(beta)+beta*p/2*(v'*C*v+sum((1-m')./m'.*v.^2.*diag(C))-2*b'*v)+beta*p/2*sigmay- gamma*sum(m)+sum(m.*log(m)+(1-m).*log(1-m));
    if iter>1
        dkl =(kl_tmp(iter)-kl_tmp(iter-1));
        if dkl>1e-10,
            eta=eta/2;
            fprintf('\t[%d] dkl=%e eta=%.3e, dm=%.3e\n',iter,kl_tmp(iter)-kl_tmp(iter-1),eta,dm);
            m=mold;
        end
    end
end;
if iter==maxiter, gamma, iter,
	figure(3)
	plot(1:iter,kl_tmp(1:iter));
end;
kl1=-p/2*log(beta)+beta*p/2*(v'*C*v+sum((1-m')./m'.*v.^2.*diag(C))-2*b'*v)+beta*p/2*sigmay- gamma*sum(m)+sum(m.*log(m)+(1-m).*log(1-m));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [qq]=jacobiana(robot,Xpunto)
    i=1;
    q = [-0.1440 -0.9810 -1.2830 0.0070 -0.0545 1.0050]';
    pause on;

    while i<99
       J = compute_jacobian(robot,q);
       qd = pinv(J)*Xpunto(:,i);
       q = q + qd;

       qq(i,:)=q; %almacenamos en una matriz todas las 'q'

       i = i + 1;
    end
end

repo: ARTE


text:
function [z,e]=exactPairMult(xIn,yIn)
%%EXACTPAIRMULT Calculate the product of the floating point double values
%               x and y such that the result z+e (if added with infinite
%               precision) is exactly the product x*y, where z is the
%               floating point product and e is the part of the result that
%               cannot be expressed exactly. This is implemented for x and
%               y being floating point doubles.
%
%INPUTS xIn, yIn Two real arrays of sizes that can be added.
%
%OUTPUTS: z The real product x.*y as limited by floating point precision.
%         e The real error such that if done with inifite precision
%           (z+e)==x*y. This does not necessarily hold when numbers become
%           denormalized.
%
%This function implements Dekker's algorithm, which is given as Mul12 in
%[1]. Note that additional scaling is performed for the instance where
%Dekker's algorithm might result in an overflow. That is when x or y is
%above 2^(1023 (maximum exponent) -53 (mantissa length))=2^970.
%
%The output pair (z,e) can be considered to be a "doublelength" floating
%point number, as defined in [1]. This means that
%abs(e)<=abs(z+e)*2^(-t)/(1+2^(-t)) (considering exact arithmetic), where t
%is the number of bits in the mantissa of the floating point number. Here,
%that is 53.
%
%Note that this assumes that the processor rounding mode has been set to
%round to "nearest," which is the default in Matlab, but which can be
%changed using the function setProcRoundingMode.
%
%EXAMPLE:
%Consider the product
% [z,e]=exactPairMult(1+2^50,1+2^(-15))
%One will get z=1.125934266580993e+15 and e=3.051757812500000e-05.
%The exact value of the product can be found to be
%1125934266580993.000030517578125
%If one adds z and e exactly, one gets precisely that result.
%
%REFERENCES:
%[1] T. J. Dekker, "A Floating Point Technique for Extending the Available
%    Precision," Numerische Mathematik, vol. 18, no. 3, Jun. 1971, pp.
%    224-242.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(~isa(xIn,'double')||~isa(yIn,'double'))
    error('This function is only implemented for x and y being double floating point values.')
end

if(isempty(xIn)||isempty(yIn))
    %Returning empty matrices if one input is an empty matrix is consistent
    %with what Matlab does when trying to multiply an empty matrix with
    %something.
    z=[];
    e=[];
    return
end

%This is 2^(1023 (maximum exponent) -53)=2^970
maxUnscaled=9.9792015476735990582818635651842e291;
two53=9007199254740992;%2^53;
twom53=1.1102230246251565404236316680908e-16;%2^(-53)

isScalX=isscalar(xIn);
isScalY=isscalar(yIn);

if(xor(isScalX,isScalY))
    %If one is scalar and the other is not.
    if(isScalX)
        dims=size(yIn);
        xIn=repmat(xIn,dims);
    else
        dims=size(xIn);
        yIn=repmat(yIn,dims);
    end
else
    %If neither is scalar or both are scalar, assume that they are the same
    %size.
    dims=size(xIn);
end

z=zeros(dims);
e=zeros(dims);
numEl=numel(z);

for curEl=1:numEl
    x=xIn(curEl);
    y=yIn(curEl);
    if(x>maxUnscaled)
        x=x*twom53;
        invScalFactX=two53;
    else
        invScalFactX=1;
    end

    if(y>maxUnscaled)
        y=y*twom53;
        invScalFactY=two53;
    else
        invScalFactY=1;
    end
    invScalFact=invScalFactX*invScalFactY;

    [zCur,eCur]=exactPairMultPrescaled(x,y);

    z(curEl)=zCur*invScalFact;
    e(curEl)=eCur*invScalFact;
end
end

function [z,e]=exactPairMultPrescaled(x,y)
    %The double mantissa is 53 bits (including the implicit one).
    const=134217729;%This is 2^ceil(53/2)+1=2^27+1

    p=x*const;
    hx=(x-p)+p;
    tx=x-hx;
    p=y*const;
    hy=(y-p)+p;
    ty=y-hy;
    p=hx*hy;
    q=(hx*ty)+(tx*hy);
    z=p+q;
    e=p-z+q+tx*ty;

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function hRD = plugin_rd_rxtx
% Reference design definition

%   Copyright 2014-2015 The MathWorks, Inc.

% Call the common reference design definition function
hRD = AnalogDevices.adrv9009.common.plugin_rd('ZCU102', 'Rx & Tx');
AnalogDevices.adrv9009.zcu102.add_io(hRD, 'Rx & Tx');

author: analogdevicesinc


text:
function [voxels,keep] = carve( voxels, camera )
%CARVE: remove voxels that are not in the silhouette
%
%   VOXELS = CARVE(VOXELS,CAMERA) carves away voxels that are not inside
%   the silhouette contained in CAMERA. The resulting voxel array is
%   returned.
%
%   [VOXELS,KEEP] = CARVE(VOXELS,CAMERA) also returns the indices of the
%   voxels that were retained.
%
%   Example:
%   >> camera = loadcameradata(1);
%   >> camera.Silhouette = getsilhouette( camera.Image );
%   >> voxels = carve( makevoxels(50), camera );
%   >> showscene( camera, voxels );
%
%   See also: LOADCAMERADATA
%             MAKEVOXELS
%             CARVEALL

%   Copyright 2005-2009 The MathWorks, Inc.
%   $Revision: 1.0 $    $Date: 2006/06/30 00:00:00 $

% Project into image
[x,y] = spacecarving.project( camera, voxels.XData, voxels.YData, voxels.ZData );

% Clear any that are out of the image
[h,w,d] = size(camera.Image); %#ok<NASGU>
keep = find( (x>=1) & (x<=w) & (y>=1) & (y<=h) );
x = x(keep);
y = y(keep);

% Now clear any that are not inside the silhouette
ind = sub2ind( [h,w], round(y), round(x) );
keep = keep(camera.Silhouette(ind) >= 1);

voxels.XData = voxels.XData(keep);
voxels.YData = voxels.YData(keep);
voxels.ZData = voxels.ZData(keep);
voxels.Value = voxels.Value(keep);


text:
%  INTERNAL FUNCTION: computes C=A*(BxBx...xB)
%
%

author: jmaih


text:
%INITCAMERAMATRIX2D  Finds an initial camera matrix from 3D-2D point correspondences
%
%     cameraMatrix = cv.initCameraMatrix2D(objectPoints, imagePoints, imageSize)
%     [...] = cv.initCameraMatrix2D(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __objectPoints__ Vector of vectors of the calibration pattern points in
%   the calibration pattern coordinate space. Cell array of cell array of
%   3-element vectors are accepted `{{[x,y,z],...}, ...}`.
% * __imagePoints__ Vector of vectors of the projections of the calibration
%   pattern points. Cell array of cell array of 2-element vectors are accepted
%   `{{[x,y],...}, ...}`.
% * __imageSize__ Image size in pixels used to initialize the principal point
%   `[w,h]`.
%
% ## Output
% * __cameraMatrix__ Camera matrix 3x3, `A = [fx 0 cx; 0 fy cy; 0 0 1]`
%
% ## Options
% * __AspectRatio__ If it is zero or negative, both `fx` and `fy` are
%   estimated independently. Otherwise, `fx = fy * AspectRatio`. default 1.0
%
% The function estimates and returns an initial camera matrix for the
% camera calibration process. Currently, the function only supports planar
% calibration patterns, which are patterns where each object point has
% z-coordinate=0.
%
% See also: cv.calibrateCamera
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function x = ellipse_rs ( r, s, i )

%*****************************************************************************80
%
%% ELLIPSE_RS maps the boundary of the unit square to an ellipse.
%
%  Discussion:
%
%    The ellipse is ( 3 * cos ( THETA ), 2 * sin ( THETA ) ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, S, the coordinates of a point that lies on the
%    boundary of the square.
%
%    Input, integer I, the component of the data which is to be returned.
%
%    Output, real X, the I-th component of the data vector X, evaluated
%    at the point (R,S), which is on a corner, or edge, of the unit square.
%
  if ( r == 0.0 )
    theta = 0.25 * pi * ( 5.0 * ( 1.0 - s ) + 3.0 * s  );
  elseif ( r == 1.0 ) then
    theta = 0.25 * pi * ( - 1.0 * ( 1.0 - s ) + 1.0 * s );
  elseif ( s == 0.0 ) then
    theta = 0.25 * pi * ( 5.0 * ( 1.0 - r ) + 7.0 * r );
  elseif ( s == 1.0 ) then
    theta = 0.25 * pi * ( 3.0 * ( 1.0 - r ) + 1.0 * r );
  end

  if ( i == 1 )

    x = 3.0 * cos ( theta );

  elseif ( i == 2 ) then

    x = 2.0 * sin ( theta );

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'ELLIPSE_RS - Fatal error!\n' );
    fprintf ( 1, '  Illegal component index I = %d\n', i );
    error ( 'ELLIPSE_RS - Fatal error!' );

  end

  return
end


text:
%  Figure 10.45      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% fig10_45.m is a script to generate Fig. 10.45 the step response for the SRL
% design for the aircraft with state feedback augmenting the
% inner loop stabilization
clf;
%system matrices
ftq =[-0.0064    0.0263         0  -32.2000      0;
   -0.0941   -0.6240  761.1400 -196.2000         0;
   -0.0002   -0.0015   -4.4120  -12.4800         0;
         0         0    1.0000         0         0;
         0   -1.0000         0  830.0000         0];
g=[0; -32.7; -2.08; 0; 0];
h=[0 0 0 0 1];
j=0;
k=[-.0011 .0016 -.0833 -1.613 -.0010];
fc=ftq-g*k;
% dc gain
dcgain=-h*inv(fc)*g;
r=100/dcgain;
t=0:.3:30;
% step response
y=step(fc,g*r,h,j,1,t);
plot(t,y);
v=[0, 30, 0, 100];
axis(v);
xlabel('Time (sec)');
ylabel('Altitude, h (ft)');
title( 'Fig. 10.45 Step response of the altitude autopilot')
%grid
nicegrid


text:
function x = equidistant1 ( n )

%*****************************************************************************80
%
%% EQUIDISTANT1 returns the Type 1 Equidistant points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2018
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), the points.
%
  x(1:n) = ( 1 - n : 2 : n - 1 ) / ( n + 1 );

  return
end


text:
function d = size(a,varargin)
% Method 'size' for file_array objects
%__________________________________________________________________________
% Copyright (C) 2005-2017-2012 Wellcome Trust Centre for Neuroimaging

%
% $Id: size.m 7440 2018-10-10 17:28:26Z john $

sa  = struct(a);
nd  = 0;
for i=1:numel(sa)
    nd = max(nd,numel(sa(i).dim));
    nd = max(nd,max(find(sa(i).pos==1)));
end
nd = nd+1;

dim = ones(length(sa),nd);
pos = ones(length(sa),nd);

for i=1:length(sa)
    sz = sa(i).dim;
    dim(i,1:length(sz)) = sz;
    ps = sa(i).pos;
    pos(i,1:length(ps)) = ps;
end

tmp = pos==1;
d   = zeros(1,nd);
for i=1:nd
    ind  = all(tmp(:,[1:(i-1) (i+1):nd]),2);
    d(i) = sum(dim(ind,i));
end
lim = max(max(find(d~=1)),2);
d   = d(1:lim);

if nargin > 1
    d_tmp  = d;
    d      = ones(size(varargin{1}));
    msk    = varargin{1}<=length(d_tmp);
    d(msk) = d_tmp(varargin{1}(msk));
end


text:
function gamma_inc_test ( )

%*****************************************************************************80
%
%% GAMMA_INC_TEST tests R4_GAMIC and R8_GAMIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAMMA_INC_TEST:\n' );
  fprintf ( 1, '  Test GAMMA_INC_VALUES, R4_GAMIC, R8_GAMIC.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             A               X     GAMIC(A,X)\n' );
  fprintf ( 1, '                                R4_GAMIC(A,X)         Diff\n' );
  fprintf ( 1, '                                R8_GAMIC(A,X)         Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, x, fx1 ] = gamma_inc_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_gamic ( single ( a ), single ( x ) );
    fx3 = r8_gamic ( a, x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.4f  %14.6g\n', a, x, fx1 );
    fprintf ( 1, '                                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end


text:
function stapx = stnddv2(XYdata,numpnts,reps,dis)
%This function approximates the standard deviation by using neighboring
%values

%This provides a limit on the number of total replicates to consider on
% neighboring points.
maxval= 15;
x= reshape([XYdata.x],length(XYdata(1).x),[])';
n = size(x,1);

if ~exist('dis','var') || isempty(dis)
    dis = pdist2(x,x);
end
%Ensure there are enough points available
if numpnts>n, numpnts=n;end
stapx=zeros(n,1);
parfor i=1:n %parfor
    %Finds the closest n points
    [~,pntsloc] = sort(dis(i,:));
    pntsloc=pntsloc(1:numpnts);
    y=[];
    leny = length(y); %this has to be done to remove a warning
    for j=1:numpnts
        %Grab y values
        y=[y XYdata(pntsloc(j)).y];
        leny = length(y);
        %If the sampling has many replicates only use close samples
        if leny>maxval
            %Only give the limit of points
            if j~=1, leny=maxval;y=y(1:maxval);end
            break;
        end
    end
    stsamps = zeros(reps,1);
    for j=1:reps
        %Calculate the standard deviation - std func is not used because
        % the mean has already been removed
        samps = randi(leny,leny,1);
        stsamps(j)=std(y(samps));
    end
    %Find the mean standard deviation
    stapx(i)=mean(stsamps);
end

end


text:
function [X, Sigma2] = Truncate(Z, tau)
% [Y, n, Sigma2] = Pro2TraceNorm(Z, tau);
% X = Z-Y;

%% new
[m, n] = size(Z);
if 2*m < n
    AAT = Z*Z';
    [S, Sigma2, D] = svd(AAT);
    Sigma2 = diag(Sigma2);
    V = sqrt(Sigma2);
    tol = max(size(Z)) * eps(max(V));
    n = sum(V > max(tol, tau));
    mid = max(V(1:n)-tau, 0) ./ V(1:n) ;
    X = (eye(m)-S(:, 1:n) * diag(mid) * S(:, 1:n)') * Z;
    return;
end
if m > 2*n
    [X, Sigma2] = Truncate(Z', tau);
    X = X';
    return;
end
[S,V,D] = svd(Z, 0);
Sigma2 = diag(V).^2;
n = sum(diag(V) > tau);
X = Z - S(:, 1:n) * diag(diag(V(1:n,1:n))-tau) * D(:, 1:n)';

author: andrewssobral


text:
%% <../index.html MATLAB xUnit Test Framework>: How to Put Multiple Test Cases in One M-file
% The Quick Start example showed how you can write a simple M-file
% to be a single test case.  This example shows you how to put multiple
% test cases in one M-file.
%
% Name your M-file beginning or ending with "test", like
% "testMyFunc".  Start by putting the following two lines at the
% beginning of the file.  It's important that the output variable
% name on line 1 be |test_suite|.
%
%    function test_suite = testMyFunc
%    initTestSuite;
%
% Next, add subfunctions to the file.  Each subfunction beginning
% or ending with "test" becomes an individual test case.
%
% The directory example_subfunction_tests contains a test M-file
% containing subfunction test cases for the |fliplr| function.

cd example_subfunction_tests

type testFliplr

%%
% As usual, run the test cases using |runtests|:

runtests

%%
% <../index.html Back to MATLAB xUnit Test Framework>

%%
% Copyright 2008-2010 The MathWorks, Inc.

author: ricedsp


text:
%% This file will be used to plot the performance comparison of SINDy-PI and PDE-Find
% Coded By: K
% Last Updated: 2019/06/23
%%
clc;clear all; close all;

%% Add path
addpath('./SINDy_PI_Results')
addpath('./PDE_Find_Results')
%% Load data
Files=dir('./SINDy_PI_Results/*.mat');

Error_DL=zeros(length(Files),1);
for k=1:length(Files)
   FileNames=strcat('SINDy_PI_Results/',Files(k).name);
   load(FileNames)
   Error_DL(k,1)=minVal_Final;
   digits(4)
   Expression
end

%%
Files=dir('./PDE_Find_Results/*.mat');

Error_PDE_Find=zeros(length(Files),1);
for k=1:length(Files)
   FileNames=strcat('PDE_Find_Results/',Files(k).name);
   load(FileNames)
   Error_PDE_Find(k,1)=minVal_Final;
   digits(4)
   Expression
end

%% Plot
close all
% Note that the first point should have zero noise, we set it up as 1e-5
% for the convenience of plotting
g0=[1e-5;0.0001;0.001;0.01;0.1;1];
figure(1)
hold on
set(gca,'FontSize',18);
set(gcf,'Position',[100 100 800 150]);
set(gcf,'PaperPositionMode','auto');
grid on

plot(g0,Error_PDE_Find,'linewidth',3,'linestyle','-','color','red')
scatter(g0,Error_PDE_Find,100,'filled','MarkerEdgeColor',[0 0 0],'MarkerFaceColor',1*[1 0 0])

plot(g0,Error_DL,'linewidth',3,'linestyle','-','color','blue')
scatter(g0,Error_DL,100,'filled','MarkerEdgeColor',[0 0 0],'MarkerFaceColor',1*[0 0 1])

xticks([1e-5 1e-4 1e-3 1e-2 1e-1 1])
xticklabels({'0','10^{-4}','10^{-3}','10^{-2}','10^{-1}','1'})

box('on')
set(gca, 'XScale', 'log', 'YScale', 'log');

author: dynamicslab

repo: SINDy-PI


text:
function SetPerFrameData(obj,fn,x,varargin)

nidx = numel(varargin);
if nidx < 1,
  error('Index not specified');
end
if nidx > 2,
  error('At most two indexing parameters allowed');
end
if ischar(varargin{1}),
  n = find(strcmp(obj.expdirs,varargin{1}),1);
  if isempty(n),
    error('Unknown experiment %s',varargin{1});
  end
elseif isnumeric(varargin{1}),
  n = varargin{1};
  if numel(n) > 1 || n < 1 || round(n) ~= n || n > obj.nexpdirs,
    error('Illegal experiment index %s',mat2str(n));
  end
else
  error('Illegal experiment index');
end
if nidx == 1,
  flies = 1:obj.nfliespermovie(n);
else
  flies = varargin{2};
  if ~isnumeric(flies) || any(flies < 1) || any(flies > obj.nfliespermovie(n)) || ...
      any(flies ~= round(flies)),
    error('Illegal fly indices');
  end
end

if iscell(x),
  ndataadd = sum(cellfun(@numel,x(1:numel(flies))));
else
  ndataadd = numel(x);
end
obj.FreeDataCache(ndataadd);

obj.nfnscached(n) = obj.nfnscached(n) + 1;
j = obj.nfnscached(n);
%fprintf('Incremented nfnscached(%d) to %d for %s\n',n,obj.nfnscached(n),fn);
obj.fnscached{n}{j} = fn;
obj.datacached{n}{j} = cell(1,numel(flies));

for flyidx = 1:numel(flies),
  fly = flies(flyidx);
  % delete data from cache if necessary
  if iscell(x),
    xcurr = x{flyidx};
  else
    xcurr = x;
  end
%   ndataadd = numel(xcurr);
%   obj.FreeDataCache(ndataadd);

  % add to cache
  %fprintf('Adding %s to data cache for fly %d at idx %d\n',fn,fly,j);
  obj.datacached{n}{j}{fly} = xcurr;

end

% update cache size
obj.ndatacached = obj.ndatacached + ndataadd;
obj.ndatacachedperexp(n) = obj.ndatacachedperexp(n) + ndataadd;

% update access time for this property
j = find(strcmp(obj.perframehistory(:,1),fn),1);
if isempty(j),
  j = size(obj.perframehistory,1)+1;
end
obj.perframehistory{j,1} = fn;
obj.perframehistory{j,2} = now;

author: kristinbranson

repo: JAABA


text:
function J = ObjectiveFCN(u,x,Ts,N,Nu,xref,u0,p,Q,R,Ru)
%% Cost function of nonlinear MPC for Lotka-Volterra system
%
% Inputs:
%   u:      optimization variable, from time k to time k+N-1
%   x:      current state at time k
%   Ts:     controller sample time
%   N:      prediction horizon
%   Nu:     control horizon [not implemented]
%   xref:   state references, varying from time k+1 to k+N
%   u0:     previous controller output at time k-1
%    p:      Parameters for model
%    Q:      State weights
%    R:      Penalization/weights on rate of change in u, uk - uk-1
%    Ru:     Control input u penalization/weights
%
% Output:
%   J:      objective function cost
%

%% Nonlinear MPC design parameters

%% Cost Calculation
% Set initial plant states, controller output and cost
xk = x;
uk = u(1);
J = 0;

% Loop through each prediction step
for ct=1:N

    % Obtain plant state at next prediction step
    xk1 = rk4u(@F8Sys,xk,uk,Ts,1,[],p);

    % Accumulate state tracking cost from x(k+1) to x(k+N)
    J = J + (xk1-xref(:,ct))'*Q*(xk1-xref(:,ct));

    % Accumulate rate of change cost from u(k) to u(k+N-1)
    if ct==1
        J = J + (uk-u0)'*R*(uk-u0) + uk'*Ru*uk;
    else
        J = J + (uk-u(ct-1))'*R*(uk-u(ct-1)) + uk'*Ru*uk;
    end

    % Update xk and uk for the next prediction step
    xk = xk1;
    if ct<N
        uk = u(ct+1);
    end
end

author: eurika-kaiser

repo: SINDY-MPC


text:
function [h_x,h_u] = sys_Dh(neq,t,x,u)

global sys_params A B C D

% h_x must be an n by n matrix.
% h_u must be an n by m matrix.

h_u = A+B*(C).*t;

h_x = B;

author: Sable


text:
function F = spm_Tcdf(x,v)
% Cumulative Distribution Function (CDF) of Students t distribution
% FORMAT p = spm_Tcdf(x,v)
%
% x - T-variate (Student's t has range (-Inf,Inf)
% v - degrees of freedom (v>0, non-integer d.f. accepted)
% F - CDF of Student's t-distribution with v degrees of freedom at points x
%__________________________________________________________________________
%
% spm_Tcdf implements the Cumulative Distribution of the Students t-distribution.
%
% Definition:
%--------------------------------------------------------------------------
% The CDF F(x) of the Student's t-distribution with v degrees of
% freedom is the probability that a realisation of a t random variable
% X has value less than x; F(x)=Pr{X<x} for X~G(h,c). Student's
% t-distribution is defined for real x and positive integer v (See
% Evans et al., Ch37).
%
% This implementation is not restricted to whole (positive integer) df
% v, rather it will compute for any df v>0.
%
% Variate relationships: (Evans et al., Ch37 & 7)
%--------------------------------------------------------------------------
% The Student's t distribution with 1 degree of freedom is the Standard
% Cauchy distribution, which has a simple closed form CDF.
%
% Algorithm:
%--------------------------------------------------------------------------
% The CDF of the Student's t-distribution with v degrees of freedom
% is related to the incomplete beta function by:
%       Pr(|X|<x) = betainc(v/(v+x^2),v/2,1/2)
% so
%              {     betainc(v/(v+x^2),v/2,1/2) / 2      for x<0
%       F(x) = |   0.5                                   for x=0
%              { 1 - betainc(v/(v+x^2),v/2,1/2) / 2      for x>0
%
% See Abramowitz & Stegun, 26.5.27 & 26.7.1; Press et al., Sec6.4 for
% definitions of the incomplete beta function. The relationship is
% easily verified by substituting for v/(v+x^2) in the integral of the
% incomplete beta function.
%
% MATLAB's implementation of the incomplete beta function is used.
%
%
% References:
%--------------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%
%__________________________________________________________________________
% Copyright (C) 1992-2011 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_Tcdf.m 4182 2011-02-01 12:29:09Z guillaume $

%-Format arguments, note & check sizes
%--------------------------------------------------------------------------
if nargin<2, error('Insufficient arguments'), end

ad = [ndims(x);ndims(v)];
rd = max(ad);
as = [[size(x),ones(1,rd-ad(1))];...
      [size(v),ones(1,rd-ad(2))]];
rs = max(as);
xa = prod(as,2)>1;
if all(xa) && any(diff(as(xa,:)))
    error('non-scalar args must match in size');
end

%-Computation
%--------------------------------------------------------------------------
%-Initialise result to zeros
F = zeros(rs);

%-Only defined for strictly positive v. Return NaN if undefined.
md = ( ones(size(x))  &  v>0 );
if any(~md(:))
    F(~md) = NaN;
    warning('Returning NaN for out of range arguments');
end

%-Special case: f is 0.5 when x=0 (where betainc involves log of zero)
F( md  &  x==0 ) = 0.5;

%-Special case: Standard Cauchy distribution when v=1
ml = ( md  &  v==1 ); if xa(1), mlx=ml; else mlx=1; end
F(ml) = 0.5 + atan(x(mlx))/pi;

%-Compute where defined & not special cases
Q  = find( md  &  x~=0  &  v~=1 );
if isempty(Q), return, end
if xa(1), Qx=Q; else Qx=1; end
if xa(2), Qv=Q; else Qv=1; end

%-Compute
xQxPos = x(Qx)>0;
F(Q) = xQxPos -(xQxPos*2-1).*0.5.*betainc(v(Qv)./(v(Qv)+x(Qx).^2),v(Qv)/2,1/2);

author: spm


text:
function r = isnan(u)
%ISNAN        Array of 1's for NaN components
%
%   r = isnan(u)
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  r = reshape( any(isnan(u.r),2) | any(isnan(u.s),2) , u.size );

repo: OpenMAS


text:
%ISHOMOG Test if SE(3) homogeneous transformation matrix
%
% ISHOMOG(T) is true (1) if the argument T is of dimension 4x4 or 4x4xN, else
% false (0).
%
% ISHOMOG(T, 'check') as above, but also checks the validity of the rotation
% sub-matrix.
%
% Notes::
% - A valid rotation sub-matrix has determinant of 1.
% - The first form is a fast, but incomplete, test for a transform is SE(3).
%
% See also ISROT, ISHOMOG2, ISVEC.

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

function h = ishomog(T, rtest)
    h = false;
    d = size(T);

    if ndims(T) >= 2
        if ~(all(d(1:2) == [4 4]))
            return %false
        end

        if nargin > 1
            for i = 1:size(T,3)
                % check rotational part
                R = T(1:3,1:3,i);
                e = R'*R - eye(3,3);
                if norm(e) > 10*eps
                    return %false
                end
                e = abs(det(R) - 1);
                if norm(e) > 10*eps
                    return %false
                end
                % check bottom row
                if ~all(T(4,:,i) == [0 0 0 1])
                    return %false
                end
            end
        end
    end

    h = true;
end

author: petercorke

repo: spatialmath-matlab


text:
function [kPa] = Pa2kPa(Pa)
% Convert units of pressure from pascals to kilopascals.
% Chad A Greene 2012
kPa = Pa/1000;


text:
function [N,E,Zone,lcm]=ell2utm(lat,lon,a,e2,lcm)
% ELL2UTM  Converts ellipsoidal coordinates to UTM coordinates.
%   UTM northing and easting coordinates must be in a 6 degree
%   system. Zones begin with zone 1 at longitude 180E to 186E
%   and increase eastward.  Formulae from E.J. Krakiwsky,
%   "Conformal Map Projections in Geodesy", Dept. Surveying
%   Engineering Lecture Notes No. 37, University of New Brunswick,
%   Fredericton, N.B, 1973. Krakwisky meridian arc length (S)
%   replaced with Hermert (1880) expansion. Vectorized.
% Version: 2011-11-13
% Useage:  [N,E,Zone,lcm]=ell2utm(lat,lon,a,e2,lcm)
%          [N,E,Zone,lcm]=ell2utm(lat,lon,a,e2)
%          [N,E,Zone,lcm]=ell2utm(lat,lon,lcm)
%          [N,E,Zone,lcm]=ell2utm(lat,lon)
% Input:   lat - vector of latitudes (rad)
%          lon - vector of longitudes (rad)
%          a   - optional ref. ellipsoid major semi-axis (m);
%                default GRS80
%          e2  - optional ref. ellipsoid eccentricity squared;
%                default GRS80
%          lcm - optional vector of non-standard central meridians
%                (rad); scalar if same for all; default = UTM def'n
% Output:  N   - vector of UTM northings (m)
%          E   - vector of UTM eastings (m)
%          Zone- vector of UTM zones (zeros if lcm specified)
%          lcm - central meridian(s) used in conversion (rad)

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

if nargin ~= 2 & nargin ~= 3 & nargin ~= 4 & nargin ~= 5
  warning('Incorrect number of input arguments');
  return
end

if nargin == 3
  lcm=a;
end
if nargin == 2 | nargin == 3
  [a,b,e2,finv]=refell('grs80');
  f=1/finv;
else
  f=1-sqrt(1-e2);
end
if nargin == 3 | nargin==5
  Zone=zeros(size(lat));
else
  Zone=floor((rad2deg(lon)-180)/6)+1;
  Zone=Zone+(Zone<0)*60-(Zone>60)*60;
  lcm=deg2rad(Zone*6-183);
end

if abs(lat)>deg2rad(80)
  warning('Latitude outside 80N/S limit for UTM');
end

ko=0.9996;           % Scale factor
No=zeros(size(lat)); % False northing (north)
No(lat<0)=1e7;       % False northing (south)
Eo=500000;           % False easting

lam=lon-lcm;
lam=lam-(lam>=pi)*(2*pi);

RN=a./(1-e2*sin(lat).^2).^0.5;
RM=a*(1-e2)./(1-e2*sin(lat).^2).^1.5;
h2=e2*cos(lat).^2/(1-e2);
t=tan(lat);
n=f/(2-f);

%----- Hinks (1927) Bessel series expansion used in DMA definition of UTM
%A0=1 - n + n.^2*5/4 - n.^3*5/4 + n.^4*81/64 - n.^5*81/64;
%A2=3/2*( n - n.^2 + n.^3*7/8 - n.^4*7/8 + n.^5*55/64 );
%A4=15/16*( n.^2 - n.^3 + n.^4*3/4 - n.^5*3/4 );
%A6=35/48*( n.^3 - n.^4 + n.^5*11/16 );
%A8=315/512*( n.^4 - n.^5 );
%S=a*(A0*lat-A2*sin(2*lat)+A4*sin(4*lat)-A6*sin(6*lat)+A8*sin(8*lat));

%----- Helmert (1880) expansion & simplification of Bessel series (faster)
A0=1+n^2/4+n^4/64;
A2=3/2*(n-n^3/8);
A4=15/16*(n^2-n^4/4);
A6=35/48*n^3;
A8=315/512*n^4;
S=a/(1+n)*(A0*lat-A2*sin(2*lat)+A4*sin(4*lat)-A6*sin(6*lat)+A8*sin(8*lat));

%----- Krakiwsky (1973) expansion
%A0=1-(e2/4)-(e2^2*3/64)-(e2^3*5/256)-(e2^4*175/16384);
%A2=(3/8)*( e2+(e2^2/4)+(e2^3*15/128)-(e2^4*455/4096) );
%A4=(15/256)*( e2^2+(e2^3*3/4)-(e2^4*77/128) );
%A6=(35/3072)*( e2^3-(e2^4*41/32) );
%A8=-(315/131072)*e2^4;
%S=a*(A0*lat-A2*sin(2*lat)+A4*sin(4*lat)-A6*sin(6*lat)+A8*sin(8*lat));

E1=lam.*cos(lat);
E2=lam.^3.*cos(lat).^3/6*(1-t.^2+h2);
E3=lam.^5.*cos(lat).^5/120.*(5-18*t.^2+t.^4+14*h2-58*t.^2.*h2+ ...
   13*h2.^2+4*h2.^3-64*t.^2.*h2.^2-24*t.^2.*h2.^3);
E4=lam.^7.*cos(lat).^7/5040.*(61-479*t.^2+179*t.^4-t.^6);
E=Eo+ko*RN.*(E1+E2+E3+E4);

N1=S./RN;
N2=lam.^2/2.*sin(lat).*cos(lat);
N3=lam.^4/24.*sin(lat).*cos(lat).^3.*(5-t.^2+9*h2+4*h2.^2);
N4=lam.^6/720.*sin(lat).*cos(lat).^5.*(61-58*t.^2+t.^4+ ...
   270*h2-330*t.^2.*h2+445*h2.^2+324*h2.^3-680*t.^2.*h2.^2+ ...
   88*h2.^4-600*t.^2.*h2.^3-192*t.^2.*h2.^4);
N5=lam.^8/40320.*sin(lat).*cos(lat).^7.*(1385-311*t.^2+543*t.^4-t.^6);
N=No+ko*RN.*(N1+N2+N3+N4+N5);


text:
function test_suite=test_correlation_measure
% tests for cosmo_correlation_measure
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_correlation_measure_basis()
    ds3=cosmo_synthetic_dataset('nchunks',3,'ntargets',4);
    ds=cosmo_slice(ds3,ds3.sa.chunks<=2);

    ds.sa.chunks=ds.sa.chunks+10;
    ds.sa.targets=ds.sa.targets+20;
    x=ds.samples(ds.sa.chunks==11,:);
    y=ds.samples(ds.sa.chunks==12,:);

    cxy=atanh(corr(x',y'));

    diag_msk=eye(4)>0;
    c_diag=mean(cxy(diag_msk));
    c_off_diag=mean(cxy(~diag_msk));

    delta=c_diag-c_off_diag;

    c1=cosmo_correlation_measure(ds);
    assertElementsAlmostEqual(delta,c1.samples,'relative',1e-5);
    assertEqual(c1.sa.labels,{'corr'});

    % reset state and do not show warnings
    orig_warning_state=cosmo_warning();
    warning_cleaner=onCleanup(@()cosmo_warning(orig_warning_state));
    cosmo_warning('reset');
    cosmo_warning('off');

    c2=cosmo_correlation_measure(ds,'output','correlation');
    assertElementsAlmostEqual(reshape(cxy,[],1),c2.samples);
    assertEqual(kron((1:4)',ones(4,1)),c2.sa.half2);
    assertEqual(repmat((1:4)',4,1),c2.sa.half1);

    i=7;
    assertElementsAlmostEqual(cxy(c2.sa.half1(i),c2.sa.half2(i)),...
                                        c2.samples(i));

    assertEqual({'half1','half2'},c2.a.sdim.labels);
    assertEqual({20+(1:4)',20+(1:4)'},c2.a.sdim.values);

    c4=cosmo_correlation_measure(ds3,'output','mean_by_fold');
    %
    for j=1:3
        train_idxs=(3-j)*4+(1:4);
        test_idxs=setdiff(1:12,train_idxs);

        ds_sel=ds3;
        ds_sel.sa.chunks(train_idxs)=2;
        ds_sel.sa.chunks(test_idxs)=1;

        c5=cosmo_correlation_measure(ds_sel,'output','mean');
        assertElementsAlmostEqual(c5.samples, c4.samples(j));
    end

    % test permutations
    ds4=cosmo_synthetic_dataset('nchunks',2,'ntargets',10);
    rp=randperm(20);

    ds4_perm=cosmo_slice(ds4,rp);
    assertEqual(cosmo_correlation_measure(ds4),...
                    cosmo_correlation_measure(ds4_perm));
    opt=struct();
    opt.output='correlation';
    assertEqual(cosmo_correlation_measure(ds4,opt),...
                    cosmo_correlation_measure(ds4_perm,opt));

function test_correlation_measure_single_target
    for ntargets=2:6
        ds=cosmo_synthetic_dataset('nchunks',2,'ntargets',ntargets);
        ds.samples=randn(size(ds.samples));
        ds.sa.targets(:)=1;

        idxs=cosmo_index_unique(mod(ds.sa.chunks,2));
        assert(numel(idxs)==2);

        x=mean(ds.samples(idxs{1},:));
        y=mean(ds.samples(idxs{2},:));
        r_xy=atanh(corr(x',y'));

        r_ds=cosmo_correlation_measure(ds,'template',1);

        assertElementsAlmostEqual(r_xy, r_ds.samples);
    end

function test_correlation_measure_regression()
    helper_test_correlation_measure_regression(false);

function test_correlation_measure_regression_spearman()
    if cosmo_skip_test_if_no_external('@stats')
        return;
    end
    helper_test_correlation_measure_regression(true);

function helper_test_correlation_measure_regression(test_spearman)
    % reset state and do not show warnings
    orig_warning_state=cosmo_warning();
    warning_cleaner=onCleanup(@()cosmo_warning(orig_warning_state));
    cosmo_warning('reset');
    cosmo_warning('off');

    ds=cosmo_synthetic_dataset('ntargets',3,'nchunks',5,'sigma',.5);
    params=get_regression_test_params(test_spearman);

    n_params=numel(params);
    for k=1:n_params
        param=params{k};

        args=param{1};
        samples=param{2};
        sa=param{3};
        sdim=param{4};
        res=cosmo_correlation_measure(ds,args{:});

        % test samples
        assertElementsAlmostEqual(res.samples,samples','absolute',5e-3);

        % test sa
        keys=fieldnames(res.sa);
        assertEqual(sort(keys(:)), sort(sa(1:2:end))');
        for j=1:2:numel(sa)
            key=sa{j};
            value=sa{j+1};
            assertEqual(res.sa.(key),value(:));
        end

        % test sa
        if isempty(sdim)
            assertFalse(isfield(res,'a'))
        else
            keys=fieldnames(res.a.sdim);
            assertEqual(sort(keys(:)), sort(sdim(1:2:end))');
            for j=1:2:numel(sa)
                key=sdim{j};
                value=sdim{j+1};
                sdim_value=res.a.sdim.(key);
                assertEqual(sdim_value(:),value(:));
            end
        end
    end

function params=get_regression_test_params(test_spearman)
    % contents
    % 1) input arguments
    % 2) samples
    % 3) sample attributes
    % 4) sdim
    if test_spearman
        params={{{ 'corr_type' 'Spearman' },...
                -0.228,...
                { 'labels' { 'corr' } },...
                []}};
    else
        params={{{ },...
                    -0.24,...
                    { 'labels' { 'corr' } },...
                    []},...
                {{ 'template' [ -2 2 3; -1 1 2; 2 -4 -3  ] },...
                    2.48,...
                    { 'labels' { 'corr' } },...
                    []},...
                {{ 'merge_func' @(x)sum(abs(x),1) },...
                    0.567,...
                    { 'labels' { 'corr' } },...
                    []},...
                {{ 'post_corr_func' @(x)x+1 },...
                    -0.204,...
                    { 'labels' { 'corr' } },...
                    []},...
                {{ 'output' 'mean_by_fold' },...
                    [ -0.289 -0.274 -0.532 -0.112 -0.269 ...
                                -0.0535 -0.203 -0.3 -0.198 -0.173 ],...
                    { 'partition' [ 1 2 3 4 5 6 7 8 9 10 ] },...
                    []},...
                {{ 'output' 'correlation' },...
                    [ -0.649 0.0675 0.345 0.126 0.643 ...
                                 0.413 0.266 0.399 0.0933 ],...
                    { 'half1', [ 1 2 3 1 2 3 1 2 3 ],...
                                'half2' [ 1 1 1 2 2 2 3 3 3 ] },...
                    { 'labels' { 'half1' 'half2' },...
                                'values' { [ 1 2 3 ]' [ 1 2 3 ]' } }
                }};

    end

function test_correlation_measure_exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                cosmo_correlation_measure(varargin{:}),'');

    % reset state and do not show warnings
    orig_warning_state=cosmo_warning();
    warning_cleaner=onCleanup(@()cosmo_warning(orig_warning_state));
    cosmo_warning('reset');
    cosmo_warning('off');

    ds=cosmo_synthetic_dataset('nchunks',2);
    aet(ds,'template',eye(4));
    aet(ds,'output','foo');
    aet(ds,'output','one_minus_correlation');

    % single target throws exception
    ds.sa.targets(:)=1;
    aet(ds);
    aet(ds,'template',2);
    aet(ds,'template',eye(2));

    ds.sa.targets(1)=2;
    aet(ds);

function x=identity(x)

function test_correlation_measure_warning_shown_if_no_defaults()
    orig_warning_state=cosmo_warning();
    cleaner=onCleanup(@()cosmo_warning(orig_warning_state));

    % reset state and do not show warnings
    cosmo_warning('reset');
    cosmo_warning('off');

    funcs={[],@atanh};
    outputs={[],'mean','raw','correlation'};

    for k=1:numel(funcs)
        func=funcs{k};
        for j=1:numel(outputs)
            output=outputs{j};

            is_default_func=k<=2;
            is_default_output=j<=2;
            expect_warning=~(is_default_func && is_default_output);

            opt=struct();
            if ~isempty(func)
                opt.post_corr_func=func;
            end

            if ~isempty(output)
                opt.output=output;
            end

            cosmo_warning('reset');
            cosmo_warning('off');

            ds=cosmo_synthetic_dataset('nchunks',2);
            cosmo_correlation_measure(ds,opt);

            s=cosmo_warning();

            showed_warning=numel(s.shown_warnings)>0;

            assertEqual(expect_warning,showed_warning);
        end
    end

function test_correlation_measure_wrong_template_size()
    ds=cosmo_synthetic_dataset('nchunks',2,'ntargets',2);

    measure_args=struct();
    measure_args.template = [ 1 -1  0  0 ; ...
                             -1  1  0  0 ; ...
                              0  0  1 -1 ;
                              0  0 -1  1 ];
    measure=@cosmo_correlation_measure;
    assertExceptionThrown(@()measure(ds,measure_args));

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function planC_for_XNAT(dicomPath,cerrPath,xhost,xproj,xsubj,xexp, rebuildRS)
% function planC_for_XNAT(dicomPath,cerrPath,xhost,xproj,xsubj,xexp,rebuildRS)
%
% This function imports DICOM files from dicomPath and adds header info
% with XNAT addressing metadata

if ~exist('rebuildRS','var')
    rebuildRS = 0;
end

if strcmpi(rebuildRS,'Y')
    rebuildRS = 1;
end

disp(['importing DICOM from ' dicomPath]);
initFlag = init_ML_DICOM;
importDICOM(dicomPath,cerrPath);

C = dir(fullfile(cerrPath, '*.mat'));
cerrFile = fullfile(C.folder, C.name);

planC = loadPlanC(cerrFile);

planC = annotatePlanCForXNAT(planC, xhost,xexp,xproj,xsubj);

strE = strsplit(xexp,'E');
numExt = str2num(strE{end});

if rebuildRS
    planC = reviveRS(planC,cerrPath,numExt);
end

save_planC(planC,[],'passed',cerrFile);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function out = iszero(f)
%ISZERO   Check if a CHEBFUN3 is identically zero on its domain.
%   OUT = ISZERO(F) returns logical 1 if the CHEBFUN3 object F is exactly
%   the zero function, and logical 0 otherwise.
%
% See also CHEBFUN2/ISZERO.

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

% Get data:
[fCore, fCols, fRows, fTubes] = tucker(f);

% Trivial check: If all the pivots are zero, then the SEPARABLEAPPROX is zero:
if ( norm(fCore(:), inf) == 0 )
    out = 1;
    return
end

% Quick check: Evaluate on a grid. If the tensor is nonzero then the
% CHEBFUN3 is nonzero.
dom = f.domain;
x = linspace(dom(1), dom(2), 10);
y = linspace(dom(3), dom(4), 10);
z = linspace(dom(5), dom(6), 10);
vals = fevalt(f, x, y, z);

if ( norm(vals(:), inf) > 0 )
   out = 0;
   return
end

% Slower check: The core may be nonzero, but the columns, rows, or tubes
% may be zero:
[rk1, rk2, rk3] = rank(f);
out_cols = zeros(rk1, 1);
out_rows = zeros(rk2, 1);
out_tubes = zeros(rk3, 1);
for j = 1:rk1
    out_cols(j) = iszero(fCols(:, j));
end
for j = 1:rk2
    out_rows(j) = iszero(fRows(:, j));
end
for j = 1:rk3
    out_tubes(j) = iszero(fTubes(:, j));
end
out = ( all(out_cols) || all(out_rows) || all(out_tubes) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [data, Ht] = rarch_simulate(T,C,parameters,p,q,type)
% Simulation of RARCH(p,q) multivariate volatility model of Noureldin, Shephard and Sheppard
%
% USAGE:
%  [DATA,HT] = rarch_simulate(T,C,PARAMETERS,P,Q,TYPE)
%
% INPUTS:
%   T          - Either a scalar containing the length of the series to simulate, or a T by K matrix
%                  of simulated random variables.  The default is to use standard normal random
%                  variables.  Providing a T by K matrix allows other distributions to be used.
%   C          - Unconditional covariance of the data
%   PARAMETERS - Vector of parameters governing the dynamics.  The form of the parameters depends on the TYPE.
%                  'Scalar':
%                  [a(1) ... a(p) b(1) ... b(q)]'  (all scalars)
%                  'CP' :
%                  [diag(A(:,:,1))' ... diag(A(:,:,p))' theta]' (theta is the scalar persistence)
%                  'Diagonal'
%                  [diag(A(:,:,1))' ... diag(A(:,:,p))' diag(B(:,:,1))' ... diag(B(:,:,p))']'
%   P          - Positive, scalar integer representing the number of symmetric innovations
%   Q          - Non-negative, scalar integer representing the number of conditional covariance
%                  lags.  When using 'CP' model, 0<=q<=1
%   TYPE       - String, one of :
%                  'Scalar' (Default)
%                  'CP' (Common Persistence)
%                  'Diagonal'
%
% OUTPUTS:
%   DATA   - A T by K matrix of simulated data
%   HT     - A [K K T] dimension matrix of conditional covariances
%
% COMMENTS:
%   The dynamics of a RARCH model are identical to that of a BEKK, except
%   that the model evolves in the rotated space.
%
%   G(:,:,t) = (eye(K) - sum(A.^2,3) - sum(B.^2,3)) +
%       A(:,:,1)*OP(:,:,t-1)*A(:,:,1) + ... A(:,:,p)*OP(:,:,t-1)*A(:,:,p) +
%       B(:,:,1)*G(:,:,t-1)*B(:,:,1) + ... B(:,:,p)*OP(:,:,t-1)*B(:,:,p)
%
%   where in the scalar model A(:,:,i) = a(i)*eye(K), B(:,:,j)=b(j)*eye(K)
%   and in the CP model, B(:,:,j) = theta - sum(A.^2,3).  OP is the outer product of the
%   unconditionally standardized data.
%
% EXAMPLES:
%   % Scalar with A.^2=.05 and B.^2=.93
%   [data,Ht] = rarch_simulate(1000,eye(2)+1,sqrt([.05,.93]),1,1,'Scalar')
%   % Diagonal
%   [data,Ht] = rarch_simulate(1000,eye(2)+1,sqrt([.05 .07 .93 .88]),1,1,'Diagonal')
%   % Common Persistence, note uses theta (sqrt(A.^2+B.^2) not B)
%   [data,Ht] = rarch_simulate(1000,eye(2)+1,sqrt([.05 .07 .99]),1,1,'CP')
%
% See also RARCH

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 3/27/2012

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k = size(C,1);
if isscalar(T)
    e = randn(2*T,k);
else
    e = T;
    if size(e,2)~=k
        error('T must have K columns when providing simulated random numbers.')
    end
    T = size(e,1);
    e = [e(ceil(rand(T,1)*T),:);e];
end

if strcmpi(type,'Scalar')
    type = 1;
elseif strcmpi(type,'CP')
    type = 2;
elseif strcmpi(type,'Diagonal')
    type = 3;
else
    error('TYPE must be ''Scalar'', ''CP'' or  ''Diagonal''.')
end

if type==2 && q>1
    error('Q must be either 0 or 1 for the ''CP'' model.')
end

switch type
    case 1
        count = p+q;
    case 2
        count = p*k+q;
    case 3
        count = (p+q)*k;
end
if length(parameters)~=count
    error('PARAMETERS does not have the expected number of elements.')
end
[C,A,B] = rarch_parameter_transform(parameters,p,q,k,C,type,false);
if max(diag(sum(A.^2,3)+sum(B.^2,3)))>=1
    warning('MFE:nonstationary','The parameters do not correspond to the stationary region.')
end
if type==2 && q==1 && min(B(:))<0
    error('When using ''CP'', the common persistence parameter Theta must satisfy Theta^2>=sum(A.^2,3)')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Gt = repmat(eye(k),[1 1 2*T]);
intercept = eye(k) - sum(A.^2,3) - sum(B.^2,3);
backCast = eye(k);
for i=1:(2*T)
    Gt(:,:,i) = intercept;
    for j=1:p
        if (i-j)<=0
            Gt(:,:,i) = Gt(:,:,i) + A(:,:,j)*backCast*A(:,:,j);
        else
            Gt(:,:,i) = Gt(:,:,i) + A(:,:,j)*(e(i-j,:)'*e(i-j,:))*A(:,:,j);
        end
    end
    for j=1:q
        if (i-j)<=0
            Gt(:,:,i) = Gt(:,:,i) + B(:,:,j)*backCast*B(:,:,j);
        else
            Gt(:,:,i) = Gt(:,:,i) + B(:,:,j)*Gt(:,:,i-j)*B(:,:,j);
        end
    end
    Gt12 = Gt(:,:,i)^(0.5);
    e(i,:) = e(i,:)*Gt12;
end

data = zeros(size(e));
Ht = zeros(k,k,T);
C12 = C^(0.5);
for i=1:length(e)
    data(i,:) = e(i,:)*C12;
    Ht(:,:,i) = C12*Gt(:,:,i)*C12;
end

data = data(T+1:2*T,:);
Ht = Ht(:,:,T+1:2*T);

author: bashtage

repo: mfe-toolbox


text:
function showScale(hAxis, i)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS
global planC

% check for comparemode if selected
compareMode = getappdata(hAxis,'compareMode');

pos = get(hAxis, 'position');
wid = pos(3); hgt = pos(4);

[scanSet, doseSet, viewType, coord] = getAxisInfo(uint8(i), 'scanSets', 'doseSets', 'view', 'coord');

scanText = ['S: ' sprintf('%d',scanSet)];
doseText = ['D: ' sprintf('%d',doseSet)];

if isempty(coord)
    set(stateS.handle.CERRAxisLabel1(i), 'string', '', 'visible', 'on', 'color', [0 0 0], 'hittest', 'off');
    set(stateS.handle.CERRAxisLabel2(i), 'string', '', 'visible', 'on', 'color', [0 0 0], 'hittest', 'off');
    return;
end

if ~isempty(scanSet)
    indexS = planC{end};
    [xV, yV, zV] = getScanXYZVals(planC{indexS.scan}(scanSet(1)));
else
    %warning('No Scan Selected')
    return
end

if stateS.currentAxis == i;
    color = [0.5 1 0.5];
else
    color = [0.9 0.9 0.5];
end
if stateS.contourAxis == i
    color = [1 0 0];
end
%try
    transM = getTransM('scan', scanSet(1), planC);
%catch
    %transM = [];
%end
if ~isempty(transM) || isequal(transM,eye(4))
    [nCoordX nCoordY nCoordZ] = applyTransM(inv(transM),coord,coord,coord);

else
    nCoordX = coord;
    nCoordY = coord;
    nCoordZ = coord;
end

switch viewType
    case 'transverse'
        viewTxt        = 'Tra: '; dim1 = 'z: '; dim2 = '\Deltax:'; dim3 = '\Deltay:';
        if isempty(zV)
            numSlices = []; zVal = [];
        else
            sliceNum    = findnearest(zV, nCoordZ);
            numSlices   = sprintf('%d',size(getScanArray(planC{indexS.scan}(scanSet(1))), 3));
            %xLim        = get(hAxis, 'xLim');
            %deltaX      = num2str(xLim(2) - xLim(1), '%0.4g');
            %yLim        = get(hAxis, 'yLim');
            %deltaY      = num2str(yLim(2) - yLim(1), '%0.4g');
            zVal        = sprintf('%.2f',coord);
            stateS.transverse.ZCoord = str2double(zVal);
%             %Show 5cm bar to display zoom-level
%             len = 5; %cm
%             dx = xLim(2)-xLim(1);
%             xStart = xLim(1) + dx * 0.05;
%             xEnd = xStart + len;
%             dy = yLim(2)-yLim(1);
%             yStart = yLim(1) + dy * 0.05;
%             yEnd = yStart + len;
%             %Delete previous handles
% %             hScale = findobj(hAxis,'tag','scale');
% %             delete(hScale)
% %             line([xStart xEnd], [yStart yStart], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %             line([xStart xStart], [yStart yEnd], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
%             xAll = linspace(xStart,xEnd,6);
%             yAll = linspace(yStart,yEnd,6);
% %             if wid/dx < 6 || hgt/dy < 6
% %                 line([xAll(6) xAll(6)], [yStart-dy*0.005 yStart+dy*0.005], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                 line([xStart-dx*0.005 xStart+dx*0.005], [yAll(6) yAll(6)], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                 text('parent', hAxis, 'string', '5', 'position', [xAll(6) yStart-dy*0.02 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %                 text('parent', hAxis, 'string', '5', 'position', [xStart-dx*0.02 yAll(6) 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %             else
% %                 for iS = 2:length(xAll)
% %                     line([xAll(iS) xAll(iS)], [yStart-dy*0.005 yStart+dy*0.005], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                     line([xStart-dx*0.005 xStart+dx*0.005], [yAll(iS) yAll(iS)], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                     text('parent', hAxis, 'string', num2str(iS-1), 'position', [xAll(iS) yStart-dy*0.02 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %                     text('parent', hAxis, 'string', num2str(iS-1), 'position', [xStart-dx*0.02 yAll(iS) 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %                 end
% %             end
%             %text('parent', hAxis, 'string', 'cm', 'position', [xStart-dx*0.02 yStart-dy*0.02 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale','rotation',0);
        end

    case 'sagittal'
        viewTxt        = 'Sag: ';  dim1 = 'x: '; dim2 = '\Deltay:'; dim3 = '\Deltaz:';
        if isempty(xV)
            numSlices = []; zVal = [];
        else
            sliceNum    = findnearest(xV, nCoordX);
            numSlices   = sprintf('%d',size(getScanArray(planC{indexS.scan}(scanSet(1))), 2));
            %xLim        = get(hAxis, 'xLim');
            %deltaX      = num2str(xLim(2) - xLim(1), '%0.4g');
            %yLim        = get(hAxis, 'yLim');
            %deltaY      = num2str(yLim(2) - yLim(1), '%0.4g');
            zVal        = sprintf('%.2f',coord);
            stateS.sagittal.ZCoord = str2double(zVal);
%             %Show 5cm bar to display zoom-level
%             len = 5; %cm
%             dx = xLim(2)-xLim(1);
%             xStart = xLim(2) - dx * 0.05;
%             xEnd = xStart - len;
%             dy = yLim(2)-yLim(1);
%             yStart = yLim(2) - dy * 0.05;
%             yEnd = yStart - len;
%             %Delete previous handles
% %             hScale = findobj(hAxis,'tag','scale');
% %             delete(hScale)
% %             line([xStart xEnd], [yStart yStart], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %             line([xStart xStart], [yStart yEnd], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
%             xAll = linspace(xStart,xEnd,6);
%             yAll = linspace(yStart,yEnd,6);
% %             if wid/dx < 6 || hgt/dy < 6
% %                 line([xAll(6) xAll(6)], [yStart-dy*0.005 yStart+dy*0.005], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                 line([xStart-dx*0.005 xStart+dx*0.005], [yAll(6) yAll(6)], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                 text('parent', hAxis, 'string', '5', 'position', [xAll(6) yStart+dy*0.02 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %                 text('parent', hAxis, 'string', '5', 'position', [xStart+dx*0.02 yAll(6) 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %             else
% %                 for iS = 2:length(xAll)
% %                     line([xAll(iS) xAll(iS)], [yStart-dy*0.005 yStart+dy*0.005], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                     line([xStart-dx*0.005 xStart+dx*0.005], [yAll(iS) yAll(iS)], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                     text('parent', hAxis, 'string', num2str(iS-1), 'position', [xAll(iS) yStart+dy*0.02 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %                     text('parent', hAxis, 'string', num2str(iS-1), 'position', [xStart+dx*0.02 yAll(iS) 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %                 end
% %             end
%             %text('parent', hAxis, 'string', 'cm', 'position', [xStart+dx*0.02 yStart+dy*0.02 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale','rotation',0);
        end

    case 'coronal'
        viewTxt        = 'Cor: ';  dim1 = 'y: '; dim2 = '\Deltax:'; dim3 = '\Deltaz:';
        if isempty(yV)
            numSlices = []; zVal = [];
        else
            sliceNum    = findnearest(yV, nCoordY);
            numSlices   = sprintf('%d',size(getScanArray(planC{indexS.scan}(scanSet(1))), 1));
            %xLim        = get(hAxis, 'xLim');
            %deltaX      = num2str(xLim(2) - xLim(1), '%0.4g');
            %yLim        = get(hAxis, 'yLim');
            %deltaY      = num2str(yLim(2) - yLim(1), '%0.4g');
            zVal        = sprintf('%.2f',coord);
            stateS.coronal.ZCoord = str2double(zVal);
%             %Show 5cm bar to display zoom-level
%             len = 5; %cm
%             dx = xLim(2)-xLim(1);
%             xStart = xLim(1) + dx * 0.05;
%             xEnd = xStart + len;
%             dy = yLim(2)-yLim(1);
%             yStart = yLim(2) - dy * 0.05;
%             yEnd = yStart - len;
%             %Delete previous handles
% %             hScale = findobj(hAxis,'tag','scale');
% %             delete(hScale)
% %             line([xStart xEnd], [yStart yStart], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %             line([xStart xStart], [yStart yEnd], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
%             xAll = linspace(xStart,xEnd,6);
%             yAll = linspace(yStart,yEnd,6);
% %             if wid/dx < 6 || hgt/dy < 6
% %                 line([xAll(6) xAll(6)], [yStart-dy*0.005 yStart+dy*0.005], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                 line([xStart-dx*0.005 xStart+dx*0.005], [yAll(6) yAll(6)], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                 text('parent', hAxis, 'string', '5', 'position', [xAll(6) yStart+dy*0.02 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %                 text('parent', hAxis, 'string', '5', 'position', [xStart-dx*0.02 yAll(6) 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %             else
% %                 for iS = 2:length(xAll)
% %                     line([xAll(iS) xAll(iS)], [yStart-dy*0.005 yStart+dy*0.005], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                     line([xStart-dx*0.005 xStart+dx*0.005], [yAll(iS) yAll(iS)], [2 2], 'parent', hAxis, 'color', 'y', 'tag', 'scale', 'hittest', 'off');
% %                     text('parent', hAxis, 'string', num2str(iS-1), 'position', [xAll(iS) yStart+dy*0.02 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %                     text('parent', hAxis, 'string', num2str(iS-1), 'position', [xStart-dx*0.02 yAll(iS) 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale');
% %                 end
% %             end
%             %text('parent', hAxis, 'string', 'cm', 'position', [xStart-dx*0.02 yStart+dy*0.02 0], 'color', 'y', 'units', 'data', 'visible', 'on','fontSize',8, 'tag', 'scale','rotation',0);
        end

    case 'legend'
        set(stateS.handle.CERRAxisLabel1(i), 'string', 'Legend', 'visible', 'on', 'color', color, 'hittest', 'off','position',[0.02 0.98 0]);
        return;
    otherwise
        return;
end

% APA commented for RIVIEW - begin

% set(stateS.handle.CERRAxisScale1(i),'xData',[xStart xEnd], 'yData', [yStart yStart],'visible','on')
% set(stateS.handle.CERRAxisScale2(i),'xData',[xStart xStart], 'yData', [yStart yEnd],'visible','on')

% set(stateS.handle.CERRAxisLabel3(i),'position',[xEnd yStart+dy*0.02 0],'visible','on');
% set(stateS.handle.CERRAxisLabel4(i),'position',[xStart-dx*0.02 yEnd 0],'visible','on');
% for j = 1:size(stateS.handle.CERRAxisTicks1,2)
%     set(stateS.handle.CERRAxisTicks1(i,j),'xData',[xAll(j) xAll(j)], 'yData', [yStart-dy*0.0025 yStart+dy*0.0025],'visible','on')
%     set(stateS.handle.CERRAxisTicks2(i,j),'xData',[xStart-dx*0.0025 xStart+dx*0.0025], 'yData', [yAll(j) yAll(j)],'visible','on')
% end

if wid < 100 || hgt < 100
    set(stateS.handle.CERRAxisLabel1(i), 'string', [viewTxt '        ' compareMode], 'visible', 'on', 'color', color, 'hittest', 'off','position',[0.02 0.98 0]);
    set(stateS.handle.CERRAxisLabel2(i),'string','');
else
    viewTxt = [viewTxt sprintf('%d',sliceNum) '/' numSlices '          ' compareMode];
    %         set(stateS.handle.CERRAxisLabel(i), 'string', {[viewTxt num2str(sliceNum) '/' numSlices], [dim1 zVal 'cm'], [dim2 deltaX 'cm'], [dim3 deltaY 'cm']}, 'visible', 'on', 'erasemode', 'none', 'color', color, 'hittest', 'off');
    set(stateS.handle.CERRAxisLabel1(i), 'string', {viewTxt, [dim1 zVal 'cm']}, 'visible', 'on', 'color', color, 'hittest', 'off','position',[0.02 0.98 0]);

    set(stateS.handle.CERRAxisLabel2(i), 'string', {scanText, doseText}, 'visible', 'on', 'color', color, 'hittest', 'off');
end

% kids = get(hAxis, 'children');
% index1 = find(kids == stateS.handle.CERRAxisLabel1(i));
% kids(index1) = [];
% index2 = find(kids == stateS.handle.CERRAxisLabel2(i));
% kids(index2) = [];
% index3 = ismember(kids, findobj(hAxis,'tag', 'scale'));
% kids(index3) = [];
% set(hAxis, 'children', [findobj(hAxis,'tag', 'scale');stateS.handle.CERRAxisLabel1(i);stateS.handle.CERRAxisLabel2(i);kids]);

% %%%%%% APA commented to test whether this is needed?
% kids = get(hAxis, 'children');
% index1 = kids == stateS.handle.CERRAxisLabel1(i);
% %kids(index1) = [];
% index2 = kids == stateS.handle.CERRAxisLabel2(i);
% %kids(index2) = [];
% %scaleKids = findobj(hAxis,'tag', 'scale');
% %index3 = ismember(kids, scaleKids);
% %kids(index3) = [];
% index3 = kids == stateS.handle.CERRAxisLabel3(i);
% index4 = kids == stateS.handle.CERRAxisLabel4(i);
% index5 = kids == stateS.handle.CERRAxisScale1(i);
% index6 = kids == stateS.handle.CERRAxisScale2(i);
% % index7 = ismember(kids, stateS.handle.CERRAxisTicks1(i,:));
% % index8 = ismember(kids, stateS.handle.CERRAxisTicks2(i,:));
% index7 = zeros(length(kids),1);
% index8 = index7;
% for j = 1:size(stateS.handle.CERRAxisTicks1,2)
%     index7 = index7 | kids == stateS.handle.CERRAxisTicks1(i,j);
%     index8 = index8 | kids == stateS.handle.CERRAxisTicks2(i,j);
% end
%
% reOrderedKids = kids;
% index = index1 | index2 | index3 | index4 | index5 | index6 | index7 | index8;
% reOrderedKids(index) = [stateS.handle.CERRAxisLabel1(i) stateS.handle.CERRAxisLabel2(i) stateS.handle.CERRAxisLabel3(i) stateS.handle.CERRAxisLabel4(i) stateS.handle.CERRAxisScale1(i) stateS.handle.CERRAxisScale2(i) stateS.handle.CERRAxisTicks1(i,:) stateS.handle.CERRAxisTicks2(i,:)];
% numFilled = sum(index);
% reOrderedKids(numFilled+1:end) = kids(~index);
%
% %%%%%% APA commented to test whether this is needed? ends

%set(hAxis, 'children', reOrderedKids);

% APA commented for RIVIEW - end


text:
function A = randsym(n,cnd)
%RANDSYM      Random symmetric matrix
%
%   res = randsym(n,cnd)
%
% cnd, if specified, is the approximate condition number of the generated matrix
%

% written  11/05/12     S.M. Rump
% modified 12/06/12     S.M. Rump  lower case
%

  if nargin==1
    A = randn(n);
  else
    s = exp( [ 0 rand(1,n-2) 1 ] * log(cnd) ) / cnd;
    A = randorth(n);
    A = A' * diag(s.*sign(randn(1,n))) * A;
  end

  % make sure A is symmetric
  A = tril(A)+tril(A,-1)';


text:
%% Example
% used to test the timing characteristic of the connection between IIWA and PC

% First start the server on the KUKA iiwa controller
% Then run this script using Matlab

% The script plots the timing data

% Copyright: Mohammad SAFEEA, 15th of June 2017

close all,clear all;clc;

ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
global t_Kuka;
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  disp('Connection could not be establised, script aborted');
  return;
else
    disp('Connection established')
end
pause(1);
disp('This will take up a minute or so to finish, please wait')

[time_stamps,time_delay]=net_updateDelay(t_Kuka);

net_turnOffServer( t_Kuka );

repo: KST-Kuka-Sunrise-Toolbox


text:
function imdb = getTIDDatabase_train(Dir)

fileID = fopen(fullfile(Dir,'\mos_with_names.txt'));
mos_name = textscan(fileID,'%f %s');
mos = mos_name{1,1};
name = mos_name{1,2};

path = cell(1,3000);
for i = 1:3000
    namet = name{i,1};
%     path{1,i} = fullfile(Dir,'distorted_images',namet);
    path{1,i} = fullfile('distorted_images',namet);
end

imdb.mos = mos';
imdb.label = imdb.mos;
imdb.name = name;
imdb.imgpath = path;
imdb.filenum = 3000;
imdb.refnum = 25;
imdb.dataset = 'TID2013';

refpath = fullfile(Dir,'reference_images');
refpath = strcat(refpath,'\');
dir_rf = dir([refpath '*.BMP']);
imdb.refname = cell(1,25);
for i = 1:25
    file_name = dir_rf(i).name;
    %imdb.refname{i} = string(file_name);
    imdb.refname{i} = file_name;
end

imdb.refnames_all = cell(1,3000);
for i = 1:25
    for j = 1:120
        imdb.refnames_all{1,(i-1)*120+j} = imdb.refname{i};
    end
end

sel = randperm(25);
train_path = [];
train_mos = [];
for i = 1:20
    train_sel = strcmpi(imdb.refname(sel(i)),imdb.refnames_all );
    train_sel = find(train_sel == 1);
    train_path = [train_path, imdb.imgpath(train_sel)];
    train_mos = [train_mos,imdb.label(train_sel)];
end

test_path = [];
test_mos = [];
for i = 21:25
    test_sel = strcmpi(imdb.refname(sel(i)),imdb.refnames_all );
    test_sel = find(test_sel == 1);
    test_path = [test_path, imdb.imgpath(test_sel)];
    test_mos = [test_mos,imdb.label(test_sel)];
end

imdb.images.id = 1:3000 ;
imdb.images.set = [ones(1,size(train_path,2)),2*ones(1,size(test_path,2))];
imdb.images.label = [train_mos,test_mos];
imdb.images.label = imdb.images.label * 100/9 ;
imdb.classes.description = {'TID2013'};
imdb.images.name = [train_path,test_path] ;
imdb.imageDir = Dir ;
fclose(fileID);

repo: DBCNN


text:
%JSINGU Show the linearly dependent joints in a Jacobian matrix
%
% JSINGU(J) displays the linear dependency of joints in a Jacobian matrix.
% This dependency indicates joint axes that are aligned and causes singularity.
%
% See also SerialLink.jacobn.

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

function jsingu(J)

    % convert to row-echelon form
    [R, jb] = rref(J);
    R(abs(R) < 100*eps) = 0;

    depcols = setdiff( 1:numcols(J), jb);

    fprintf('%d linearly dependent joints:\n', length(depcols));
    for d=depcols
        fprintf('  q%d depends on: ', d)
        for k=find(R(:,d))
            fprintf('q%d ', k);
        end
        fprintf('\n');
    end

repo: robotics-toolbox-matlab


text:
% The signless Laplacian matrix of a graph
% Def: the sum of the diagonal degree matrix and the adjacency matrix
%
% INPUTS: adjacency matrix, nxn
% OUTPUTs: signless Laplacian matrix, nxn
%
% GB: last updated, Dec 6 2015

function L=signlessLaplacian(adj)

L=diag(sum(adj))+adj;

author: aeolianine

repo: octave-networks-toolbox


text:
function transferModel(modelFile, srcVocabFile_new, tgtVocabFile_new, srcCharPrefix_new, tgtCharPrefix_new, outModelFile, varargin)
  addpath(genpath(sprintf('%s/../..', pwd)));

  %% Argument Parser
  p = inputParser;
  % required
  addRequired(p,'modelFile',@ischar);
  addRequired(p,'srcVocabFile_new',@ischar);
  addRequired(p,'tgtVocabFile_new',@ischar);
  addRequired(p,'srcCharPrefix_new',@ischar);
  addRequired(p,'tgtCharPrefix_new',@ischar);
  addRequired(p,'outModelFile',@ischar);

  % optional
  addOptional(p,'gpuDevice', 0, @isnumeric); % choose the gpuDevice to use: 0 -- no GPU

  p.KeepUnmatched = true;
  parse(p, modelFile, srcVocabFile_new, tgtVocabFile_new, srcCharPrefix_new, tgtCharPrefix_new, outModelFile, varargin{:})
  transferParams = p.Results;

  % GPU settings
  transferParams.isGPU = 0;
  if transferParams.gpuDevice
    n = gpuDeviceCount;
    if n>0 % GPU exists
      fprintf(2, '# %d GPUs exist. So, we will use GPUs.\n', n);
      transferParams.isGPU = 1;
      gpuDevice(transferParams.gpuDevice)
      transferParams.dataType = 'single';
    else
      transferParams.dataType = 'double';
    end
  else
    transferParams.dataType = 'double';
  end
  printParams(2, transferParams);

  modelFile = transferParams.modelFile;
  srcVocabFile_new = transferParams.srcVocabFile_new;
  tgtVocabFile_new = transferParams.tgtVocabFile_new;
  srcCharPrefix_new = transferParams.srcCharPrefix_new;
  tgtCharPrefix_new = transferParams.tgtCharPrefix_new;
  outModelFile = transferParams.outModelFile;

  %% Transfer
  % load current data
  fprintf(2, '# Loading %s\n', modelFile);
  savedData = load(modelFile);
  model = savedData.model;
  params = savedData.params;

  % src word
  fprintf(2, '# Transfering src word\n');
  [model.W_emb_src, params.srcVocab, ~, params.srcVocabFile] = transferMatrix(model.W_emb_src, params.srcVocab, params.srcCharShortList, ...
    srcVocabFile_new, params, 1, 0);

  % tgt word
  fprintf(2, '# Transfering tgt word\n');
  [model.W_emb_tgt, ~, ~, ~] = transferMatrix(model.W_emb_tgt, params.tgtVocab, params.tgtCharShortList, ...
    tgtVocabFile_new, params, 1, 0);

  % tgt word soft
  fprintf(2, '# Transfering tgt word soft\n');
  [model.W_soft, params.tgtVocab, ~, params.tgtVocabFile] = transferMatrix(model.W_soft, params.tgtVocab, params.tgtCharShortList, ...
    tgtVocabFile_new, params, 0, 0);
  % IMPORTANT: only update tgtVocab in the last call

  % src char
  if isfield(params, 'charOpt') && params.charOpt == 1
    fprintf(2, '# Transfering src char\n');
    [model.W_emb_src_char, params.srcCharVocab, params.srcCharVocabSize, params.srcCharVocabFile] = transferMatrix(model.W_emb_src_char, params.srcCharVocab, ...
      params.srcCharVocabSize, [srcCharPrefix_new '.char.vocab'], params, 1, 1);
    params.srcCharMapFile = [srcCharPrefix_new '.char.map'];
    params.srcCharMap = loadWord2CharMap(params.srcCharMapFile, params.charMaxLen);
  end

  if isfield(params, 'charOpt') && (params.charOpt == 2 || params.charOpt == 3)
    % tgt char
    fprintf(2, '# Transfering tgt char\n');
    [model.W_emb_tgt_char, ~, ~, ~] = transferMatrix(model.W_emb_tgt_char, params.tgtCharVocab, ...
      params.tgtCharVocabSize, [tgtCharPrefix_new '.char.vocab'], params, 1, 1);
    params.tgtCharMapFile = [tgtCharPrefix_new '.char.map'];
    params.tgtCharMap = loadWord2CharMap(params.tgtCharMapFile, params.charMaxLen);

    % tgt char soft
    fprintf(2, '# Transfering tgt char soft\n');
    [model.W_soft_char, params.tgtCharVocab, params.tgtCharVocabSize, params.tgtCharVocabFile] = transferMatrix(model.W_soft_char, params.tgtCharVocab, ...
      params.tgtCharVocabSize, [tgtCharPrefix_new '.char.vocab'], params, 0, 1);
    % IMPORTANT: only update tgtCharVocab in the last call
  end

  % save model
  save(outModelFile, 'model', 'params');
  model
  fprintf(2, '# Done saving to %s\n', outModelFile);
end

function [W_new, vocab_new, vocabSize_new, vocabFile_new] = transferMatrix(W, vocab, shortList, vocabFile_new, params, isCol, isChar)
  fprintf(2, '  W [%s], isCol=%d, isChar=%d, original vocab size %d, short list %d\n', num2str(size(W)), isCol, isChar, length(vocab), shortList);
  [vocab_new] = loadVocab(vocabFile_new);
  if isChar % for chars, we learn for the entire char vocab, so shortList is NOT used
    vocab_new{end+1} = '<c_s>'; % not learn
    vocab_new{end+1} = '</c_s>';

    if isCol
      W_new = initMatrixRange(params.initRange, [params.lstmSize, length(vocab_new)], params.isGPU, params.dataType);
    else
      W_new = initMatrixRange(params.initRange, [length(vocab_new), params.lstmSize], params.isGPU, params.dataType);
    end
    learnedVocab_new = vocab_new;
    vocabMap = cell2map(vocab);
  else % for words, we only learn embeddings for a small shortlist
    assert(shortList < length(vocab));
    if isCol
      W_new = initMatrixRange(params.initRange, [params.lstmSize, shortList], params.isGPU, params.dataType);
    else
      W_new = initMatrixRange(params.initRange, [shortList, params.lstmSize], params.isGPU, params.dataType);
    end
    learnedVocab_new = vocab_new(1:shortList);
    vocabMap = cell2map(vocab(1:shortList));
  end
  vocabSize_new = length(vocab_new);

  % mapping
  flags = isKey(vocabMap, learnedVocab_new);
  indices = values(vocabMap, learnedVocab_new(flags));
  indices = [indices{:}];
  if isCol
    W_new(:, flags) = W(:, indices);
  else
    W_new(flags, :) = W(indices, :);
  end

  fprintf(2, '  vocabSize_new %d, num overlap %d\n  new words:', vocabSize_new, sum(flags));
  fprintf(2, ' %s', learnedVocab_new{~flags});
  fprintf(2, '\n');
end

author: lmthang