text:
%% Structured grid
% This is a sample code showing how to use ConstructProjector3D().
% This example is going to setup a transformation matrix based on forth
% order polynomials in 3D.
% Pay attention that no data is needed to construct the interpolator. Only
% the coordinates of the points on source and destination grid.
clear;clc;
%% Initializing Part I
disp('- Initializing')
xMin=0;
xMax=2*pi;
yMin=0;
yMax=2*pi;
zMin=0;
zMax=2*pi;

nx=20;
ny=20;
nz=20;

dx=(xMax-xMin)/(nx-1);
dy=(yMax-yMin)/(ny-1);
dz=(zMax-zMin)/(nz-1);

nPoly=2; % We gonna fit a fourth order polynomial,
         % i.e. sum_{n,m,t=0}^{n,m,t=nPoly} x^n*y^m*z^t
nInterp=27:37;

%% Turning off the warnings
% Depending on the parameter that you choose the matrices might be nearly
% singular. In this example, despite being nearly singular, they are not
% affecting the final interpolation. Therefore, I set the warning to off.
% But generally do not turn it off and if you got this warning check if it
% is going to affect your results or not.
% disp('- Turning the warnings off')
% warning('OFF','MATLAB:nearlySingularMatrix');

%% Generating the Source Grid
disp('- Generating the source grid.')
[xn,yn,zn]=meshgrid(linspace(xMin,xMax,nx),linspace(yMin,yMax,ny),linspace(zMin,zMax,nz));

%% Finding the cell centers
disp('- Generating the destination grid')
xc=xn(1:ny-1,1:nx-1,1:nz-1)+dx/2;
yc=yn(1:ny-1,1:nx-1,1:nz-1)+dy/2;
zc=zn(1:ny-1,1:nx-1,1:nz-1)+dz/2;

% figure
% plot3(xn(:),yn(:),zn(:),'k.','MarkerSize',20);
% hold on
% plot3(xc(:),yc(:),zc(:),'b.','MarkerSize',20);
% axis tight;
% axis square;
% title('Source Grid, black dots, & the destination grid, blue dots.', ...
%       'FontName','Arial','FontSize',12,'FontWeight','Bold');

%% Generarting some Data
disp('- Generating some data and interpolating')
F1=@(x,y,z) (sqrt(x.^2+y.^2+z.^2));
Vn=F1(xn,yn,zn);
Vc_Analytic=F1(xc,yc,zc);

%% Generating the interpolant
disp('- Generating the interpolant and testing nInterp')
RMSE=zeros(size(nInterp));
for i=1:numel(nInterp)
  disp(['-- nInterp= ' num2str(nInterp(i))])
  P=ConstructProjector3D(xn(:),yn(:),zn(:),xc(:),yc(:),zc(:),nPoly,nInterp(i));
  % Interpolating
  Vc_interp=P*Vn(:);

  % Calculating the RMSE
  RMSE(i)=sqrt(mean( (Vc_interp-Vc_Analytic(:)).^2 ));
  disp(['   RMSE= ' num2str(RMSE(i))])
end
%% plotting

figure
plot(nInterp,RMSE,'b.-')
xlabel('nInterp');
ylabel('RMSE');
axis tight

disp('All done!')


% Valid ending arguments:
% - SequenceNumber (0 - 65535)
% - StateInfo0 (1 byte)
% - StateInfo1 (1 byte)
% - SerialNumber (14-16 characters)
% - Longitude (double)
% - Latitude (double)
% - Height (-32768 to 32767)
% - Altitude (-32768 to 32767)
% - VelocityNorth (-32768 to 32767)
% - VelocityEast (-32768 to 32767)
% - VelocityUp (-32768 to 32767)
% - Yaw (-32768 to 32767)
% - PhoneAppGPSTime (milliseconds since Epoch)
% - PhoneAppLatitude (double)
% - PhoneAppLongitude (double)
% - HomeLatitude (double)
% - HomeLongitude (double)
% - ProductType (1 byte)
% - UUID (19 bytes)

function [frame] = create_frame_bytes(varargin)
%     MIN_SAMPLE_RATE = 15.36e6;
%
%     assert(sample_rate >= MIN_SAMPLE_RATE, "Sample rate must be at least 15.36 MSPS");
%     assert(round(sample_rate / MIN_SAMPLE_RATE) * MIN_SAMPLE_RATE == sample_rate, ...
%         "Sample rate must be multiple of %d", MIN_SAMPLE_RATE);

    assert(isempty(varargin) || mod(length(varargin), 2) == 0, ...
        "Variable argument list must have even number of elements");

    sequence_number = 0;
    state_info = [0, 0];
    serial_number = '0123456789abcd';
    longitude = 0.0;
    latitude = 0.0;
    height = 0;
    altitude = 0;
    velocity_north = 0;
    velocity_east = 0;
    velocity_up = 0;
    yaw = 0;
    phone_app_gps_time = 0;
    phone_app_latitude = 0.0;
    phone_app_longitude = 0.0;
    home_latitude = 0.0;
    home_longitude = 0.0;
    product_type = 0;
    uuid = zeros(1, 19);

    for idx=1:2:length(varargin)
        key = varargin{idx};
        val = varargin{idx + 1};

        switch(key)
            case 'SequenceNumber'
                assert(isnumeric(val), 'Sequence number must be an integer');
                assert(val >= -32768 && val <= 32767, 'Invalid sequence number');

                sequence_number = uint16(val);
            case 'StateInfo0'
                assert(isinteger(val), 'StateInfo0 must be an integer');
                assert(val >= 0 && val <= 255, 'Invalid StateInfo0 value');

                state_info(1) = uint8(val);
            case 'StateInfo1'
                assert(isinteger(val), 'StateInfo1 must be an integer');
                assert(val >= 0 && val <= 255, 'Invalid StateInfo1 value');

                state_info(2) = uint8(val);
            case 'SerialNumber'
                assert(isstring(val) || ischar(val), 'SerialNumber must be a string')
                assert(length(val) >= 14 && length(val) <= 16, ...
                    'SerialNumber length %d is not valid.  Must be between 14 and 16 chars', length(val))

                serial_number = val;
            case 'Longitude'
                assert(isnumeric(val), 'Longitude must be a number')
                assert(abs(val) <= 180, 'Longitude must be between -180 and 180')

                longitude = double(val);
            case 'Latitude'
                assert(isnumeric(val), 'Latitude must be a number')
                assert(abs(val) <= 90, 'Latitude must be between -90 and 90')

                latitude = double(val);
            case 'Height'
                assert(isinteger(val), 'Height must be an integer')
                assert(val >= -32768 && val <= 32767, 'Height must be between -32768 and 32767')

                height = int16(val);
            case 'Altitude'
                assert(isinteger(val), 'Altitude must be an integer')
                assert(val >= -32768 && val <= 32767, 'Altitude must be between -32768 and 32767')

                altitude = int16(val);
            case 'VelocityNorth'
                assert(isinteger(val), 'VelocityNorth must be an integer')
                assert(val >= -32768 && val <= 32767, 'VelocityNorth must be between -32768 and 32767')

                velocity_north = int16(val);
            case 'VelocityEast'
                assert(isinteger(val), 'VelocityEast must be an integer')
                assert(val >= -32768 && val <= 32767, 'VelocityEast must be between -32768 and 32767')

                velocity_east = int16(val);
            case 'VelocityUp'
                assert(isinteger(val), 'VelocityUp must be an integer')
                assert(val >= -32768 && val <= 32767, 'VelocityUp must be between -32768 and 32767')

                velocity_up = int16(val);
            case 'Yaw'
                assert(isinteger(val), 'Yaw must be an integer')
                assert(val >= -32768 && val <= 32767, 'Yaw must be between -32768 and 32767')

                yaw = int16(val);
            case 'PhoneAppGPSTime'
                assert(isinteger(val), 'PhoneAppGPSTime must be an integer')
                assert(val >= 0 && val <= (2^64)-1, 'PhoneAppGPSTime must be between 0 and (2^64)-1')

                phone_app_gps_time = uint64(val);
            case 'PhoneAppLatitude'
                assert(isnumeric(val), 'PhoneAppLatitude must be a number')
                assert(abs(val) <= 90, 'PhoneAppLatitude must be between -90 and 90')

                phone_app_latitude = double(val);
            case 'PhoneAppLongitude'
                assert(isnumeric(val), 'PhoneAppLongitude must be a number')
                assert(abs(val) <= 180, 'PhoneAppLongitude must be between -180 and 180')

                phone_app_longitude = double(val);
            case 'HomeLatitude'
                assert(isnumeric(val), 'HomeLatitude must be a number')
                assert(abs(val) <= 90, 'HomeLatitude must be between -90 and 90')

                home_latitude = double(val);
            case 'HomeLongitude'
                assert(isnumeric(val), 'HomeLongitude must be a number')
                assert(abs(val) <= 180, 'HomeLongitude must be between -180 and 180')

                home_longitude = double(val);
            case 'ProductType'
                assert(isinteger(val), 'ProductType must be an integer')
                assert(val >= 0 && val <= 255, 'ProductType must be between 0 and 255')

                product_type = uint8(val);
            case 'UUID'
                assert(isstring(val) || ischar(val), 'UUID must be a string')
                assert(length(val) == 19, 'UUID must be 19 characters long (THIS IS TEMPORARY!!)')

                uuid = val;
            otherwise
                error('Invalid key "%s"', key);
        end
    end

    coord_adj = 10000000 / 57.2957795785523;

    % Pack the fields into bytes
    sequence_number_bytes = to_bytes(sequence_number, 2);
    state_info_bytes = state_info;
    serial_number_bytes = [uint8(serial_number) zeros(1, 16 - length(serial_number))];
    longitude_bytes = to_bytes(int32(round(longitude * coord_adj)), 4);
    latitude_bytes = to_bytes(int32(round(latitude * coord_adj)), 4);
    height_bytes = to_bytes(height, 2);
    altitude_bytes = to_bytes(altitude, 2);
    velocity_north_bytes = to_bytes(velocity_north, 2);
    velocity_east_bytes = to_bytes(velocity_east, 2);
    velocity_up_bytes = to_bytes(velocity_up, 2);
    yaw_bytes = to_bytes(yaw, 2);
    phone_app_gps_time_bytes = to_bytes(phone_app_gps_time, 8);
    phone_app_latitude_bytes = to_bytes(int32(round(phone_app_latitude * coord_adj)), 4);
    phone_app_longitude_bytes = to_bytes(int32(round(phone_app_longitude * coord_adj)), 4);
    home_latitude_bytes = to_bytes(int32(round(home_latitude * coord_adj)), 4);
    home_longitude_bytes = to_bytes(int32(round(home_longitude * coord_adj)), 4);
    product_type_bytes = product_type;
    uuid_bytes = uint8(uuid);

    version = uint8(2);
    message_type = uint8(16);
    uuid_len = uint8(length(uuid_bytes));

    frame = [ ...
        message_type, version, sequence_number_bytes, state_info_bytes, serial_number_bytes, longitude_bytes, ...
        latitude_bytes, height_bytes, altitude_bytes, velocity_north_bytes, velocity_east_bytes, velocity_up_bytes, ...
        yaw_bytes, phone_app_gps_time_bytes, phone_app_latitude_bytes, phone_app_longitude_bytes, home_longitude_bytes, ...
        home_latitude_bytes, product_type_bytes, uuid_len, uuid_bytes, uint8(0) ...
    ];

    frame = [uint8(length(frame)) frame];

    crc_bytes = to_bytes(calculate_crc(frame), 2);
    frame = [frame, crc_bytes];
end


text:
function z = sparse( i, j, x, m, n )

%   Disciplined convex/geometric programming information for SPARSE:
%       For a CVX variable X, SPARSE(X) just returns X. In the three-
%       and five-argument versions SPARSE(I,J,V,M,N), the index
%       arguments I and J and size arguments M and N must be constant.
%       If any of the index pairs (I(k),J(k)) are repeated, then the
%       corresponding elements V(k) will be added together. Therefore,
%       those elements must be compatible with addition; i.e., they
%       must have the same curvature.

switch nargin,
	case 1,
		z = i;
		return
	case {2,4}
		cvx_throw( 'Not enough arguments.' );
end

%
% Check sizes and indices
%

if ( ~isnumeric( i ) && ~ischar( i ) ) || ( ~isnumeric( j ) && ~ischar( j ) ),
    cvx_throw( 'The first two arguments must be numeric.' );
end
nn = [ numel( i ), numel( j ), prod( x.size_ ) ];
nz = max( nn );
if any( nn(nn~=nz) ~= 1 ),
    cvx_throw( 'Vectors must be the same lengths.' );
end
i = i( : ); j = j( : );
if nargin == 3,
    m = max( i );
    n = max( j );
elseif ( ~isnumeric( m ) && ~ischar( m ) ) || ( ~isnumeric( n ) && ~ischar( n ) ) || n < 0 || m < 0 || n ~= floor( n ) || m ~= floor( m ),
    cvx_throw( 'Sparse matrix sizes must be positive integers.' );
elseif any( i > m ) || any( j > n ),
    cvx_throw( 'Index exceeds matrix dimensions.' );
end

%
% Reconstruct basis matrices and indices
%

[ ix, jx, vx ] = find( x.basis_ );
ij = i + m * ( j - 1 );
if length(ij) > 1,
    ij = ij(jx);
end
xb = sparse( ix, ij, vx, max(ix), m * n );
z = cvx( [ m, n ], xb );

%
% Verify vexity is preserved
%

if ~cvx_isvalid(z),
    cvx_throw( 'Disciplined convex programming error:\n    Sparse matrix construction produced invalid sums of convex and concave terms.' );
end

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

author: yu-jiang

repo: radpbook

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


text:
report_this_filefun(mfilename('fullpath'));

clf
[X,Y,Z] = meshgrid(gx,gy,gz);
v = [];
main = [ -174.83 51.412 33];
i = 0.05
ii = 0.05
for i3 = z0:1*dz:z1
    l = loc(3,:) == i3;
    re3 = reshape(as(l),length(gy),length(gx));
    r = reshape(loc(4,l),length(gy),length(gx));
    l = r > tresh;
    re3(l) = zeros(1,length(find(l)))*nan;
    rect = [i ii 0.15 0.15];
    i = i+0.20;
    if i > 0.9 ; i = 0.05 ; ii = ii + 0.20; end
    axes('position',rect)
    pcolor(gx,gy,re3)
    caxis([min(as) max(as)])
    colormap(jet)
    title([ 'Depth ' num2str(i3) ' km'], 'FontSize',8)
    set(gca,'FontSize',8')

    hold on
    l = a.Depth > i3-dz/2 & a.Depth < i3+dz/2;
    plot3(a(l,1),a(l,2),a(l,7)*0,'k.')
    shading interp
    % axis('off')
    pl =plot3(main(:,1),main(:,2),main(:,3)*0,'xk');
    set(pl,'LineWidth',2)
    hold on

end

author: CelsoReyes


text:
function []=hddisp(n,endo,Y,decimaldates1,hd_estimates,stringdates1,T,pref,IRFt,signreslabels,Fstartlocation)

% function []=hddisp(n,endo,decimaldates1,hd_estimates,stringdates1,T,datapath)
% plots the results for the historical decomposition
% inputs:  - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - cell 'endo': list of endogenous variables of the model
%          - vector 'decimaldates1': dates converted into decimal values, for the sample period
%          - cell 'hd_estimates': lower bound, point estimates, and upper bound for the historical decomposition
%          - cell 'stringdates1': date strings for the sample period
%          - integer 'T': number of sample time periods (defined p 7 of technical guide)
%          - string 'datapath': user-supplied path to excel data spreadsheet
% outputs: none

% transpose the cell of records (required for the plot function in order to obtain correctly ordered plots)
hd_estimates=hd_estimates';

% loop over all shocks for all endogenous variables
for ii=1:n*(n+1)
% positive and negative distributions are calculated seperately for
% graphical Matlab issues
contributions(:,ii)=[hd_estimates{ii}(2,:)'];
contribpos(:,ii)=contributions(:,ii);
contribpos(contribpos<0)=0;
contribneg(:,ii)=contributions(:,ii);
contribneg(contribneg>0)=0;
end

% calculate the sum of all contributions to proxy the actual development of
% the endogenous variable

Total = zeros(length(decimaldates1),n);
for i=1:n
    Total(:,i) = sum(contributions(:,(n+1)*(i-1)+1:(n+1)*i-1)');
end

% Colors
myC= [0 0 1
1 1 0
1 0.4 0
0 0.8 1
0 1 0
1 0 1
0 1 1
1 0.2 0.5
0.7 1 0.7
0.7 0.7 1
0.8 0.3 0.6
0.3 0.8 0.6
0.4 0.7 0.8
];

if pref.plot
ncolumns=ceil(n^0.5);
nrows=ceil(n/ncolumns);
hd1=figure('Tag','BEARresults');
for i=1:n
set(hd1,'name','Historical decomposition');
subplot(nrows,ncolumns,i)
hdpos=bar(decimaldates1, contribpos(:,n*(i-1)+i:(n+1)*i-1), 0.8, 'stacked');
hold on
hdneg=bar(decimaldates1, contribneg(:,n*(i-1)+i:(n+1)*i-1), 0.8, 'stacked');
hold on
plot(decimaldates1,Total(:,i),'k','LineWidth',2.8);
hold on
for k=1:n
    set(hdpos(k),'facecolor', myC(k,:), 'Edgecolor', 'none');
    set(hdneg(k),'facecolor', myC(k,:), 'Edgecolor', 'none');
end
if nargin==11
   yy=get(gca,'YLim');
   Xpatch=[ones(1,2)*decimaldates1(Fstartlocation) ones(1,2)*decimaldates1(T)];
   Ypatch=[yy fliplr(yy)];
   FHDpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
   set(FHDpatch,'facealpha',0.2);
   set(FHDpatch,'edgecolor','k','linewidth',0.1);
end
axis tight
hold off
% label the endogenous variables
title(endo{i,1})
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
%box off
end
if IRFt==4
   legend(signreslabels,'Location','Northoutside','Orientation','Vertical')
else
   legend(endo);
end
legend boxoff
end % pref.plot

% finally, record results in excel
% retranspose the cell of records
hd_estimates=hd_estimates';
bear.data.excelrecord7fcn(n, T ,IRFt, signreslabels, endo, stringdates1, hd_estimates, pref)

author: european-central-bank

repo: BEAR-toolbox


text:
function B = filter_base( )
%%%%%%%%%%%%%%%%%%%%DCT base
 config;
 fS = nnconfig.FilterSize ;
 fN = nnconfig.FilterNumber;
 fS_sqrt = fS^2;

 DCT = dctmtx(fS);
 DCT = kron(DCT, DCT);
 B = zeros(fS_sqrt, fN);
 for i = 2 : fS_sqrt
 B(:, i-1) = DCT(i, :);
 end

end

repo: Deep-ADMM-Net


text:
function tform = geotrans(varargin)
%GEOTRANS Make affine and projective geometric transformations.
%   TFORM = GEOTRANS(TYPE,P1,P2,___) makes a geometric transformation
%   with the specified type and parameters.
%
%   Valid values for TYPE and P1, P2, ..., are:
%
%   'scale'     Scale. If only P1 is provided, it is the scale
%               factor in both directions. If P1 and P2 are provided,
%               then P1 is the horizontal scale factor and P2 is the
%               vertical scale factor.
%
%   'rotate'    Rotation. P1 is the rotation angle, measured in
%               degrees counterclockwise from the positive horizontal
%               axis.
%
%   'translate' Translation. P1 is the horizontal translation and P2
%               is the vertical translation.
%
%   'h-reflect' Horizontal reflection.
%
%   'v-reflect' Vertical reflection.
%
%   'h-shear'   Horizontal shear. P1 is the sheer angle in degrees,
%               measured counterclockwise from the downward-pointing
%               y-axis.
%
%   'v-shear'   Vertical shear. P1 is the sheer angle in degrees,
%               measured clockwise from the horizontal axis.
%
%
%   'compose'   Composition of transformations. P1, P2, ..., are
%               affine2d and projective2d transformations.
%
%   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

if ischar(varargin{1}) || isstring(varargin{1})
   operation = varargin{1};
   switch operation
      case 'scale'
         narginchk(2,3);
         tform = scale(varargin{2:end});
      case 'rotate'
         narginchk(2,2);
         tform = rotate(varargin{2});
      case 'translate'
         narginchk(3,3);
         tform = translate(varargin{2:end});
      case 'h-reflect'
         narginchk(1,1);
         tform = reflectHorizontal;
      case 'v-reflect'
         narginchk(1,1);
         tform = reflectVertical;
      case 'h-shear'
         narginchk(2,2);
         tform = shearHorizontal(varargin{2});
      case 'v-shear'
         narginchk(2,2);
         tform = shearVertical(varargin{2});
      case 'compose'
         narginchk(2,Inf);
         tform = compose(varargin{2:end});
      otherwise
         error('Unrecognized operation: %s',operation);
   end
else
   error('Input argument TYPE must be a string or char vector.');
end

%----------------------------------------------------------------------%
function tform = scale(sx,sy)
if nargin < 2
   sy = sx;
end
tform = affine2d([sx 0 0; 0 sy 0; 0 0 1]);

%----------------------------------------------------------------------%
function tform = rotate(angle)
theta = (pi/180) * angle;
T = [ ...
   cos(theta)  sin(theta)  0
   -sin(theta) cos(theta)  0
   0           0           1 ];
tform = affine2d(T);

%----------------------------------------------------------------------%
function tform = translate(dx,dy)
T = [ ...
   1  0  0
   0  1  0
   dx dy 1 ];
tform = affine2d(T);

%----------------------------------------------------------------------%
function tform = reflectHorizontal
T = [ ...
   -1  0  0
   0  1  0
   0  0  1 ];
tform = affine2d(T);

%----------------------------------------------------------------------%
function tform = reflectVertical
T = [ ...
   1  0  0
   0  -1 0
   0  0  1];
tform = affine2d(T);

%----------------------------------------------------------------------%
function tform = shearVertical(angle_degrees)
alpha = tand(angle_degrees);
T = [ ...
   1 alpha 0
   0 1     0
   0 0     1 ];
tform = affine2d(T);

%----------------------------------------------------------------------%
function tform = shearHorizontal(angle_degrees)
alpha = tand(angle_degrees);
T = [ ...
   1     0  0
   alpha 1  0
   0     0  1];
tform = affine2d(T);

%----------------------------------------------------------------------%
function tform = compose(varargin)
all_affine = true;
for k = 1:nargin
   if ~(isa(varargin{k},'affine2d') || ...
         isa(varargin{k},'projective2d'))
      error('Input transforms must be affine2d or projective2d.');
   end

   if ~isa(varargin{k},'affine2d')
       all_affine = false;
   end

end
T = eye(3);
for k = 1:nargin
   tform_k = varargin{k};
   T = T * tform_k.T;
end
if all_affine
    tform = affine2d(T);
else
    tform = projective2d(T);
end

author: dipum

repo: dipum-toolbox

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


text:
function [f, finv, pol, polinv] = conformal(C, varargin)
%% CONFORMAL  Conformal map to unit disk
%   [F, FINV] = CONFORMAL(C, ctr) computes a conformal map F of the region
%   bounded by the complex periodic chebfun C to the unit disk and its inverse
%   FINV, with F(ctr) = 0 and F'(ctr) > 0.  Both maps are represented by
%   function handles evaluating rational functions.  If ctr is omitted it is
%   set to 0.
%
%   CONFORMAL(..., 'tol', tol) uses tolerance tol instead of the default 1e-5.
%
%   CONFORMAL(..., 'plots') produces plots of the map and its inverse.
%
%   CONFORMAL(..., 'numbers') prints various quantities.
%
%   CONFORMAL(..., 'poly') uses a less robust algorithm based on polynomials
%                          instead of the Kerzman-Stein integral equation.
%
%   [F, FINV, POL, POLINV] = CONFORMAL(...) returns the poles of F and FINV.
%
%   This experimental code is good for smooth simple regions, but easy to break.
%
%   Examples:
%
%   C = chebfun('exp(pi*1i*t)','trig')*(1+.15*randnfun(.2,'trig'));
%   [f, finv] = conformal(C, 'plots');                               % random
%
%   C = chebfun('2*cos(t)+1i*(sin(t)+2*cos(t).^3)',[0 2*pi],'trig'); % Ellacott
%   conformal(C, 'poly', 'plots');                                   % blade
%
%   C = chebfun('exp(pi*1i*t)*(1+.3*cos(6*pi*t))','trig');
%   conformal(C, 'plots', 'tol', 1e-8, 'numbers');                   % snowflake
%
%   s = chebfun('s'); C = join(s-.5i,1+.5i*s,.5i-s,-1-.5i*s);
%   conformal(C, 'poly', 'plots');                                   % rectangle
%
% References:
%
% A. Gopal and L. N. Trefethen, Representation of conformal maps by
% rational functions, Numer. Math. 142 (2019), 359--382.
%
% L. N. Trefethen, Numerical conformal mapping with rational
% functions, Comp. Meth. Funct. Th. 20 (2020), 369-387.
%
% See also CONFORMAL2.

%%
%   Default algorithm:
%     (1) Solve discretized Kerzman-Stein integral equation
%     (2) Use AAA to approximate f and its inverse by rational functions
%
%   Alternative algorithm if 'poly' is specified:
%     (1) Use least-squares to find harmonic u s.t. u(z) = -log(abs(z-ctr)) on C
%     (2) Set w = u+iv (v = harmonic conjugate of u)
%     (3) Set f = z*exp(w(z-ctr));
%     (4) Use AAA to approximate f and its inverse by rational functions
%
%   Although in principle these algorithms should be embedded in the
%   Chebfun constructor, for simplicity in this numerically challenging
%   area we have not done that.
%
% This code was written by L. N. Trefethen in September 2019.  The
% Kerzman-Stein part originates with Anne Greenbaum and Trevor Caldwell.

t1 = tic;
[ctr, tol, plots, numbers, poly] = parseinputs(C, varargin{:});
err = Inf;
scl = norm(C-ctr, inf);

if poly == 0               % DEFAULT ALGORITHM: KERZMAN-STEIN INTEGRAL EQUATION

    M = 300;
    while err > tol
        M = M + 300;
        [g, Z, W] = kerzstein((C-ctr)/scl, M, 0);
        Z = Z*scl + ctr;
        gc = trigcoeffs(g);
        err = norm(gc([1:10 end-9:end]));           % a crude error measure
        if (err > tol) && (M >= 1200)
            warning('CONFORMAL did not converge')
            break
        end
    end

else                       % ALTERNATIVE ALGORITHM: POLYNOMIAL EXPANSION

    w1 = warning('off', 'MATLAB:rankDeficientMatrix');
    logn = 4;
    dom = domain(C);
    dom = dom([1 end]);
    while err > tol
        n = round(2^logn);                    % degree of polynomial
        M = 8*n;                              % number of sample points
        logn = logn + 1/2;
        Z = C(dom(1) + (1:M)'*diff(dom)/M);   % sample points
        Zscl = (Z-ctr)/scl;                   % rescale to improve conditioning
        G = -log(abs(Zscl));                  % RHS for Dirichlet problem
        n1 = 0:n;                             % exponents for real part
        n2 = 1:n;                             % of exponents for imag part
        Q = ones(size(Zscl));                 % Arnoldi.  Q has orthog cols of norm sqrt(M).
        H = zeros(n+1,n);
        for k = 1:n
            v = Zscl.*Q(:,k);
            v = v - Q*(Q'*v)/M;               % or execute twice for better orthogonality!
            H(k+1,k) = norm(v)/sqrt(M);       % At end, Zscl.*Q(:,1:n)/M = Q*H
            Q = [Q v/H(k+1,k)];               % Q has orthog cols of norm sqrt(M)
        end
        A = [real(Q) imag(Q(:,2:end))];
        c = A\G;                              % soln of least-squares problem
        err = norm(A*c-G, inf);               % max error
        cc = c(n1+1)-1i*[0; c(n+1+n2)];       % coeffs for harmonic -> analytic
        W = Zscl.*exp(Q*cc);                  % images on boundary
        if (err > tol) && (logn >= 9.5)
            warning('CONFORMAL did not converge')
            break
        end
    end
    warning(w1.state, 'MATLAB:rankDeficientMatrix')

end

w2 = warning('off', 'CHEBFUN:aaa:Froissart');
[f0, pol] = aaa(W, Z, 'tol', tol);            % forward map
zz  = 1e-4*scl*[1 1i -1 -1i];                 % finite diff for simplicity
dwdz = sum(f0(ctr+zz)./zz);                   % derivative at ctr
rot = exp(-1i*angle(dwdz));                   % rotation for f'(ctr) > 0
f = @(z) rot*f0(z);                           % rotate mapping function
W = rot*W;                                    % rotate points on circle
[finv, polinv] = aaa(Z, W, 'tol', tol);       % inverse map
warning(w2.state, 'CHEBFUN:aaa:Froissart')

inC = inpolygon(real(pol), imag(pol), ...     % check for poles of f in region
	real(Z), imag(Z));
if max(inC) > 0
warning('CONFORMAL: pole in region')
end
if min(abs(polinv)) < 1                       % check for poles of finv in disk
    warning('CONFORMAL: pole in disk')
end
tcomp = toc(t1);

% Plot solution if requested
if plots
    t2 = tic;
    clf
    LW = 'linewidth'; PO = 'position'; MS = 'markersize';
    FW = 'fontweight'; NO = 'normal'; XT = 'xtick'; YT = 'ytick';
    circ = exp(2i*pi*(0:300)/300);

    h1 = axes(PO, [.04 .38 .45 .53]);         % plot C and poles of f
    plot(C, 'b', LW, 1)
    axis([real(ctr)+1.4*scl*[-1 1] imag(ctr)+1.4*scl*[-1 1]])
    axis square, hold on, plot(pol, '.r', MS, 8)
    title([num2str(length(pol)) ' poles'],FW,NO)

    h2 = axes(PO, [.52 .38 .45 .53]);         % plot disk and poles of finv
    plot(circ, 'b', LW, 1), axis(1.6*[-1 1 -1 1])
    axis square, hold on, set(gca,XT,-1:1,YT,-1:1)
    plot(polinv, '.r', MS, 8)
    title([num2str(length(polinv)) ' poles'],FW,NO)

    ncirc = 8;                  % plot concentric circles and their images
    for r = (1:ncirc-1)/ncirc
        axes(h1), plot(finv(r*circ), '-k', LW, .5)
        axes(h2), plot(r*circ, '-k', LW, .5)
    end

    nrad = 16;                  % plot radii and their images
    ray = chebpts(301);
    ray = ray(ray>=0);
    for k = 1:nrad
        axes(h1), plot(finv(ray*exp(2i*pi*k/nrad)), '-k', LW, .5)
        axes(h2), plot(ray*exp(2i*pi*k/nrad), '-k', LW, .5)
    end

    axes(h1), hold off
    axes(h2), hold off
    tplot = toc(t2);

end

% Print various quantities if requested
if numbers
    disp(' ')
    fprintf('                         computation time in seconds:%6.2f\n', tcomp)
    if plots
        fprintf('                            plotting time in seconds:%6.2f\n', tplot)
    end
    fprintf('            number of sample points Z on boundary, M:  %d\n', M)
    fprintf('                      numbers of poles of f and finv:  %d, %d\n', ...
        length(pol), length(polinv))
    fprintf('back-and-forth boundary error norm(Z-finv(f(Z)),inf):  %6.1e\n',...
        norm(Z-finv(f(Z)), inf))
    fprintf(' inverse back-and-forth error norm(W-f(finv(W)),inf):  %6.1e\n',...
        norm(W-f(finv(W)), inf))
    fprintf(' interior inverse error norm(.9*W-f(finv(.9*W)),inf):  %6.1e\n',...
        norm(.9*W-f(finv(.9*W)), inf))
    if poly
        fprintf('             max error of least-squares problem, err:  %6.1e\n', err)
        fprintf('                                polynomial degree, n:  %d\n', n)
        fprintf('                number of real degrees of freedom, N:  %d\n', 2*n+1)
        fprintf('          condition number of least-squares matrix A:  %6.1e\n', cond(A))
    else
        fprintf('                            rough error measure, err:  %6.1e\n', err)
    end
    disp(' ')
end

end   % of conformal

function [ctr, tol, plots, numbers, poly] = parseinputs(C, varargin)
ctr = 0;                        % defaults
tol = 1e-5;
plots = 0;
numbers = 0;
poly = 0;
j = 1;
while j < nargin
    j = j+1;
    v = varargin{j-1};
    if isnumeric(v)
        ctr = v;
    elseif strcmp(v, 'tol')
        j = j+1;
        tol = varargin{j-1};
    elseif strcmp(v, 'plots')
        plots = 1;
    elseif strcmp(v, 'numbers')
        numbers = 1;
    elseif strcmp(v, 'poly')
        poly = 1;
    end
end
winding_number = sum(diff(C)/(C-ctr))/(2i*pi);
if abs(winding_number - 1) > tol
    error('CONFORMAL:parseinputs','C must wind once counterclockwise around ctr')
end
end   % end of parseinputs

function [g, Z, W] = kerzstein(C, M, ctr)
%
% Given a smooth periodic chebfun C defining the boundary of a region Omega
% in the counterclockwise direction, this function computes the boundary
% correspondence function g for a conformal mapping from Omega to the unit
% disk, mapping ctr to 0.  It first finds M points (Z) on C that are equally
% spaced in arclength.  It then solves an integral equation to determine the
% images of these points (W) on the unit circle.  Finally, a periodic
% chebfun g is defined to map the arclengths associated with Z to their
% images on the unit circle.  The original version of this code cones
% from Anne Greenbaum and Trevor Caldwell.
% Reference:  Kerzman and Trummer, "Numerical conformal mapping via the Szego
% kernel," Journal of Computational and Applied Mathematics 14 (1986), 111-123.

% Express C as a function of arclength and compute its derivative

S = arcLength(C);
s = cumsum(abs(diff(C)));
dC_arg = angle(diff(C));
dC_arg = unwrap(dC_arg - dC_arg(0));
u = inv(s);
C = newDomain(C, minandmax(u));
D = C(u);
Dprime = diff(D);

%  Set up and solve the integral equation

ds = S/M;
svec = [0:M-1]'*ds;
Dvec = D(svec);
Z = Dvec;
Dprimevec = Dprime(svec);
gamdotvec = Dprimevec ./ abs(Dprimevec);   % unit tangents
d = 1/(2i*pi);
gvec = d*conj(gamdotvec./(ctr - Dvec));
A = eye(M);
for j = 1:M
    w = Dvec(j);
    for i = 1:M
        z = Dvec(i);
        if i ~= j
            A(i,j) = A(i,j) - d*(conj(gamdotvec(i)/(z-w)) + gamdotvec(j)/(w-z))*ds;
        end
    end
end
fvec = A\gvec;

%  Compute boundary correspondence function g

Rprimevec = fvec.^2;
Rvec = -1i*gamdotvec.*(Rprimevec./abs(Rprimevec));
W = Rvec;
g = chebfun(Rvec, [0,S], 'trig');

end   % of kerzstein

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [recall, rankloss, allRecalls, opts]= testFromFn(db, dbFeatFn, qFeatFn, opts, varargin)

    if nargin<4 || isempty(opts)
        % a bit hacky but fine..
        opts= struct(...
            'nTestRankSample', 0, ...
            'nTestSample', inf, ...
            'recallNs', [1:5, 10:5:100], ...
            'margin', 0.1, ...
            'nNegChoice', 1000, ...
            'cropToDim', 0 ...
            );
    end
    opts= vl_argparse(opts, varargin);

    relja_display('testFromFn:\n%s\n%s', dbFeatFn, qFeatFn);

    qFeat= fread( fopen(qFeatFn, 'rb'), inf, 'float32=>single');
    qFeat= reshape(qFeat, [], db.numQueries);
    nDims= size(qFeat, 1);
    dbFeat= fread( fopen(dbFeatFn, 'rb'), [nDims, db.numImages], 'float32=>single');
    assert(size(dbFeat,2)==db.numImages);

    if isfield(opts, 'cropToDim') && opts.cropToDim>0
        if opts.cropToDim > nDims
            warning('cropToDim (%d) larger than the dimensionality (%d) -- ignoring', opts.cropToDim, nDims);
        else
            qFeat= relja_l2normalize_col( qFeat(1:opts.cropToDim, :) );
            dbFeat= relja_l2normalize_col( dbFeat(1:opts.cropToDim, :) );
        end
    end

    if opts.nTestRankSample>0
        rankloss= testCoreRank(db, qFeat, dbFeat, opts.margin, opts.nNegChoice, 'nTestSample', opts.nTestRankSample);
    else
        rankloss= [];
    end
    [recall, allRecalls]= testCore(db, qFeat, dbFeat, 'nTestSample', opts.nTestSample, 'recallNs', opts.recallNs);
end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
function [ Ix, Iy, Iz, It ] = imageDerivatives3D( image1, image2 )
%This fucntion computes 3D derivatives between two 3D images.
%
%   Description :
%
%   There are four derivatives here; three along X, Y, Z axes and one along
%   timeline axis.
%
%   -image1, image2 :   two subsequent images or frames
%   -dx, dy, dz : vectors along X, Y and Z axes respectively
%   -dt : vectors along timeline axis
%   -Ix, Iy, Iz : derivatives along X, Y and Z axes respectively
%   -It : derivatives along timeline axis
%
%   Author : Mohammad Mustafa
%   By courtesy of The University of Nottingham and
%   Mirada Medical Limited, Oxford, UK
%
%   Published under a Creative Commons Attribution-Non-Commercial-Share Alike
%   3.0 Unported Licence http://creativecommons.org/licenses/by-nc-sa/3.0/
%
%   June 2012

dx=zeros(2,2,2);
dx(:,:,1)=[-1 1; -1 1 ]; dx(:,:,2)=[-1 1; -1 1 ];
dx=0.25*dx;

dy=zeros(2,2,2);
dy(:,:,1)=[-1 -1; 1 1 ]; dy(:,:,2)=[-1 -1; 1 1 ];
dy=0.25*dy;

dz=zeros(2,2,2);
dz(:,:,1)=[-1 -1; -1 -1 ]; dz(:,:,2)=[1 1; 1 1 ];
dz=0.25*dz;

dt=ones(2,2,2);
dt=0.25*dt;

% Computing derivatives
Ix = 0.5 * (convn(image1,dx) + convn(image2,dx) );
Iy = 0.5 * (convn(image1,dy) + convn(image2,dy) );
Iz = 0.5 * (convn(image1,dz) + convn(image2,dz) );
It = 0.5 * (convn(image1,dt) - convn(image2,dt) );

% Adjusting sizes
Ix=Ix(1:size(Ix,1)-1, 1:size(Ix,2)-1, 1:size(Ix,3)-1);
Iy=Iy(1:size(Iy,1)-1, 1:size(Iy,2)-1, 1:size(Iy,3)-1);
Iz=Iz(1:size(Iz,1)-1, 1:size(Iz,2)-1, 1:size(Iz,3)-1);
It=It(1:size(It,1)-1, 1:size(It,2)-1, 1:size(It,3)-1);

end

author: Sable

repo: mcbench-benchmarks


text:
function logProb = vargplvmProbabilityCompute(model, y, display, iters)

% VARGPLVMPROBABILITYCOMPUTE description

% VARGPLVM

% Takes an input a trained vargplvm and a test data point (with possibly missing values)
% Computes the probability density in the test data point

% Indices of missing dimension
indexMissing = isnan(y);
indexPresent = ~indexMissing;
nTestSamples = size(y,1);
logProb = zeros(nTestSamples,1);

% if ~matlabpool('size')
%     if display, fprintf('vargplvmProbabilityCompute.m: line 17 -- Starting up MatLab pool for parallel computing.'); end
%     matlabpool
% end
try
    nWorkers = matlabpool('size');
catch e
    nWorkers = 1;
end
if nWorkers==0 % If matlabpool has not started, don't try to use it.
    nWorkers = 1;
end

% compute the variational lower without the new data point
Fold = vargplvmLogLikelihood(model);

if all(all(repmat(indexMissing(1,:), [nTestSamples, 1]) == indexMissing))
    % Faster option:
    if display,
        fprintf('Optimisation will be done all at once as the missing variables are always the same (or if there is no missing variable).\n');
        fprintf('Initialising the latent point using the nearest neighbour from the training data.\n');
    end
    dst = dist2(y(:,indexPresent(1,:)), model.y(:,indexPresent(1,:)));
    [~, mini] = min(dst,[],2);

    if display, fprintf('Creating the variational distribtion for the test latent point.\n'); end
    vardistx = vardistCreate(model.vardist.means(mini,:), model.q, 'gaussian');

    if display, fprintf('Optimising over the latent points.\n'); end
    model.vardistx = vardistx;
    [X, varX] = vargplvmOptimisePoint(model, vardistx, y, display, iters);

    if display, fprintf('Computing the variational lower with the new data points included.\n'); end
    vardistx.means = X;
    vardistx.covars = varX;

    if display, fprintf('Generating results for each sample separately (so optimisation is done only once for all test samples).\n'); end
    parfor t = 1:nTestSamples,
        vardistxTmp = vardistx;
        vardistxTmp.numData = 1;
        % if display, fprintf('Testing sample %d/%d\n', t, nTestSamples);
        vardistxTmp.means = vardistx.means(t,:);
        vardistxTmp.covars = vardistx.covars(t,:) ;
        Fnew = vargplvmPointLogLikelihood(model, vardistxTmp, y(t,:));

        % compute the probability
        logProb(t) = Fnew - Fold;
    end

else % This version is slower, but also uses multithreading:
    c = clock;
    tmpDir = sprintf('/tmp/%s/%d_%d_%d_%d_%d', getenv('USER'), c(1), c(2), c(3), c(4), c(5));
    if display,
        fprintf('There are missing variables but they are not consistent. Optimising for each test point separately.\n');
        fprintf('In order to reduce the memory usage by each worker, I will use temporary files in %s\n', tmpDir);
    end
    if ~exist(tmpDir,'dir')
        mkdir(tmpDir);
    end

    logProbCells = cell(nWorkers,1);
    for w=1:nWorkers
        tmpFile = sprintf('%s/Ytest_cpu%02d.mat', tmpDir, w);
        if ~exist(tmpFile,'file')
            currSamplesBegin = 1 + (w-1)*floor(nTestSamples/nWorkers);
            currSamplesEnd = w*floor(nTestSamples/nWorkers);
            if w==nWorkers
                currSamplesEnd = nTestSamples;
            end
            sampleIndices = [currSamplesBegin currSamplesEnd];
            if display, fprintf('Writing samples for worker %d/%d: from %d to %d...', w, nWorkers, currSamplesBegin, currSamplesEnd); end
            currY = y(currSamplesBegin:currSamplesEnd,:);
            currIndexPresent = indexPresent(currSamplesBegin:currSamplesEnd,:);
            save(tmpFile, 'currY', 'currIndexPresent', 'sampleIndices');
            if display, fprintf(' done!\n'); end
        end
    end
    clear y;
    clear indexPresent;

    parfor w=1:nWorkers,
      tmpFile = sprintf('%s/Ytest_cpu%02d.mat', tmpDir, w);
      if display, fprintf('Loading %s...', tmpFile); end
      testData = load(tmpFile);
      if display, fprintf(' done!\n'); end

      currSamples = (testData.sampleIndices(1):testData.sampleIndices(end))';
      currLogProb = zeros(currSamples(end)-currSamples(1)+1,1);

      for tIdx=1:length(currSamples),
        modelCopy = model;
        t = currSamples(tIdx);
        if display, fprintf('-----\nvargplvmProbabilityCompute(): evaluating sample %d of %d\n', t, nTestSamples); end;

        dst = dist2(testData.currY(tIdx,testData.currIndexPresent(tIdx,:)), modelCopy.y(:,testData.currIndexPresent(tIdx,:)));
        [~, mini] = min(dst);

        % create the variational distribtion for the test latent point
        vardistx = vardistCreate(modelCopy.vardist.means(mini,:), modelCopy.q, 'gaussian');

        % optimize over the latent point
        modelCopy.vardistx = vardistx;
        [X, varX] = vargplvmOptimisePoint(modelCopy, vardistx, testData.currY(tIdx,:), display, iters);

        % compute the variational lower with the new data point included
        vardistx.means = X;
        vardistx.covars = varX;

        % Old way:
        Fnew = vargplvmPointLogLikelihood(modelCopy, vardistx, testData.currY(tIdx,:));

        % compute the probability
        currLogProb(tIdx) = Fnew - Fold;
      end
      logProbCells{w} = currLogProb;
    end

    % Concatenating the results of all workers:
    for w=1:nWorkers
      currSamples = ((1 + (w-1)*floor(nTestSamples/nWorkers)) : w*floor(nTestSamples/nWorkers))';
      if w==nWorkers && (currSamples(end) < nTestSamples),
          currSamples = [currSamples; ((currSamples(end)+1):nTestSamples)'];
      end
      logProb(currSamples) = logProbCells{w};
    end
end
end

author: SheffieldML

repo: vargplvm


text:
function tests = test_spm_get_lm
% Unit Tests for spm_get_lm
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_get_lm.m 6534 2015-08-24 16:02:56Z guillaume $

tests = functiontests(localfunctions);

function test_spm_get_lm_2D(testCase)
dim     = [5 5];
vol     = rand(dim(1),dim(2));
subvol  = vol(2:4,2:4);
vol(3,3) = max(subvol(:)) + 0.1;
[X,Y] = ndgrid(1:dim(1),1:dim(2));
gm      = sub2ind(size(vol),3,3);
idx     = spm_get_lm(vol,[X(:)';Y(:)']);
testCase.verifyTrue(any(idx == gm));

function test_spm_get_lm_3D(testCase)
dim     = [5 5 5];
vol     = rand(dim(1),dim(2),dim(3));
subvol  = vol(2:4,2:4,2:4);
vol(3,3,3) = max(subvol(:)) + 0.1;
[X,Y,Z] = ndgrid(1:dim(1),1:dim(2),1:dim(3));
gm      = sub2ind(size(vol),3,3,3);
idx     = spm_get_lm(vol,[X(:)';Y(:)';Z(:)']);
testCase.verifyTrue(any(idx == gm));

author: spm


text:
function binv = beta_inv (x, a, b)
%  BETA_INV Inverse of the beta cumulative distribution function (cdf).
%     X = BETA_INV(P,A,B) returns the inverse of the beta cdf with
%     parameters A and B at the values in P.

% Copyright (c) 1995-1997, 2005-2007 Kurt Hornik

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if (nargin ~= 3)
  error('must provide 3 parameters')
end

if (~isscalar (a) || ~isscalar(b))
  if ((size(a,1) ~= size(b,1)) || (size(b,1) ~= size(x,1)))
    error ('betainv: x, a and b must be of common size or scalars');
  end
end

[n,nin] = size(x);
binv = zeros(n,nin);

k = find ((x<0) || (x>1) || (a<0) || (b<0) || isnan(x));
if (any(k))
  binv(k) = NaN;
end

k = find ((x==1) && (a>0) && (b>0));
if (any (k))
  binv(k) = 1;
end

k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0));
if (any (k))
  if (~isscalar(a) || ~isscalar(b))
    a = a(k);
    b = b(k);
    y = a./(a+b);
  else
    y = a/(a + b)*ones(size(k));
  end
  x = x(k);
  l = find (y < eps);
  if (any (l))
    y(l) = sqrt (eps) * ones (length (l), 1);
  end
  l = find (y > 1 - eps);
  if (any (l))
    y(l) = 1 - sqrt (eps) * ones (length (l), 1);
  end

  y_old = y;
  for i = 1 : 10000
    h     = (beta_cdf (y_old, a, b) - x) ./ beta_pdf (y_old, a, b);
    y_new = y_old - h;
    ind   = find (y_new <= eps);
    if (any (ind))
      y_new (ind) = y_old (ind) / 10;
    end
    ind = find (y_new >= 1 - eps);
    if (any (ind))
      y_new (ind) = 1 - (1 - y_old (ind)) / 10;
    end
    h = y_old - y_new;
    if (max (abs (h)) < sqrt (eps))
      break;
    end
    y_old = y_new;
  end

  binv(k) = y_new;
end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function bti_test06 ( )

%*****************************************************************************80
%
%% BTI_TEST06 tests BURGERS_TIME_INVISCID.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BTI_TEST06\n' );
  fprintf ( 1, '  Test BURGERS_TIME_INVISCID\n' );

  method = 6;
  nx = 81;
  nt = 200;
  t_max = 2.0;
  bc = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Method: 6, Godonov.\n' );
  fprintf ( 1, '  Initial condition: expansion\n' );
  fprintf ( 1, '  Number of space nodes = %d\n', nx );
  fprintf ( 1, '  Number of time steps = %d\n', nt );
  fprintf ( 1, '  Final time T_MAX = %g\n', t_max );
  fprintf ( 1, '  Boundary condition = %d\n', bc );

  U = burgers_time_inviscid ( method, @ic_expansion, nx, nt, t_max, bc );

  x = linspace ( -1.0, +1.0, nx );

  figure ( 6 )

  plot ( x, U(1:50:(nt+1),:), 'Linewidth', 3 )
  grid on
  xlabel ( '<-- X -->' )
  ylabel ( '<-- U(X,T) -->' )
  title ( 'Burgers, inviscid, initial expansion, Godonov' )

  filename = 'bti_test06.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plot as "%s"\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function value = t_polynomial_value ( n, x )

%*****************************************************************************80
%
%% T_POLYNOMIAL_VALUE: returns the single value T(n,x).
%
%  Discussion:
%
%    In cases where calling T_POLYNOMIAL is inconvenient, because it returns
%    a vector of values for multiple arguments X, this simpler interface
%    may be appropriate.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the polynomial.
%
%    Input, real X, the argument of the polynomial.
%
%    Output, real VALUE, the value of T(n,x).
%
  m = 1;

  v_vec = t_polynomial ( m, n, x );

  value = v_vec(n+1);

  return
end


text:
data_path = '/home/ftp/Coldog/Dataset/toflow/vimeo_septuplet/sequences';
noise_path = '/home/ftp/Coldog/Dataset/toflow/vimeo_septuplet/sequences_blur';
if ~exist(noise_path, 'dir')
  mkdir(noise_path);
end
videos = dir(data_path);
total_num = 91701;
count = 0;
disp('Start bluring...')
for i=1:length(videos)
    if strcmp(videos(i).name, '.') || strcmp(videos(i).name, '..')
        continue
    end
    seps = dir(fullfile(data_path, videos(i).name));
    for j=1:length(seps)
        if strcmp(seps(j).name, '.') || strcmp(seps(j).name, '..')
            continue
        end
        blur_input(fullfile(data_path, videos(i).name, seps(j).name), fullfile(noise_path, videos(i).name, seps(j).name))
        count = count + 1;
        if mod(count,500)==0
            disp(['  Processed ',num2str(count/total_num*100), '%'])
        end
    end
end
disp(['  Processed ',num2str(count/total_num*100), '%'])

repo: pytoflow


text:
function [K, grad] = covSEard(hyp, x, z, i)

% Squared Exponential covariance function with Automatic Relevance Detemination
% (ARD) distance measure. The covariance function is parameterized as:
%
% k(x^p,x^q) = sf2 * exp(-(x^p - x^q)'*inv(P)*(x^p - x^q)/2)
%
% where the P matrix is diagonal with ARD parameters ell_1^2,...,ell_D^2, where
% D is the dimension of the input space and sf2 is the signal variance. The
% hyperparameters are:
%
% hyp = [ log(ell_1)
%         log(ell_2)
%          .
%         log(ell_D)
%         log(sqrt(sf2)) ]
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2010-09-10.
%
% See also COVFUNCTIONS.M.

if nargin<2, K = '(D+1)'; return; end              % report number of parameters
if nargin<3, z = []; end                                   % make sure, z exists
xeqz = numel(z)==0; dg = strcmp(z,'diag') && numel(z)>0;        % determine mode

[n,D] = size(x);
ell = exp(hyp(1:D));                               % characteristic length scale
sf2 = exp(2*hyp(D+1));                                         % signal variance

% precompute squared distances
if dg                                                               % vector kxx
  K = zeros(size(x,1),1);
else
  if xeqz                                                 % symmetric matrix Kxx
    K = sq_dist(diag(1./ell)*x');
  else                                                   % cross covariances Kxz
    K = sq_dist(diag(1./ell)*x',diag(1./ell)*z');
  end
end

K = sf2*exp(-K/2);                                                  % covariance

if nargout >= 2
	grad = bsxfun(@times, diag(1./(ell.^2))*bsxfun(@minus, x', z'), K');
end

if nargin>3                                                        % derivatives
  if i<=D                                              % length scale parameters
    if dg
      K = K*0;
    else
      if xeqz
        K = K.*sq_dist(x(:,i)'/ell(i));
      else
        K = K.*sq_dist(x(:,i)'/ell(i),z(:,i)'/ell(i));
      end
    end
  elseif i==D+1                                            % magnitude parameter
    K = 2*K;
  else
    error('Unknown hyperparameter')
  end
end

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo


text:
filename='Cantilever_tetrahedra_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
%   - data                 Hand, Omega=(0,20)x(0,25),
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             SSD
%   - regularizer          mfCurvature
%   - optimizer            Gauss-Newton
% ===============================================================================

clear; close all; help(mfilename)

setup2DhandData
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',0.01);
distance('reset','distance','SSD');
trafo('reset','trafo','affine2D');
regularizer('reset','regularizer','mfCurvature','alpha',1000);

yc = MLIR(ML,'parametric',1);
showResults(ML,yc);

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


text:
function X = fast_pca(in_X, K)

in_X = in_X - repmat(mean(in_X),size(in_X,1),1);
[U, S, ~] = rsvd(in_X, K);
K = min(size(S,2),K);
X = U(:,1:K)*diag(sqrt(diag(S(1:K,1:K))));
X = X./repmat(sqrt(sum(X.*X,2)),1,K);
end

function [U,S,V] = rsvd(A,K)
%-------------------------------------------------------------------------------------
% random SVD
% Extremely fast computation of the truncated Singular Value Decomposition, using
% randomized algorithms as described in Halko et al. 'finding structure with randomness
%
% usage :
%
%  input:
%  * A : matrix whose SVD we want
%  * K : number of components to keep
%
%  output:
%  * U,S,V : classical output as the builtin svd matlab function
%-------------------------------------------------------------------------------------
% Antoine Liutkus  (c) Inria 2014

[M,N] = size(A);
P = min(2*K,N);
X = randn(N,P);
Y = A*X;
W1 = orth(Y);
B = W1'*A;
[W2,S,V] = svd(B,'econ');
U = W1*W2;
K=min(K,size(U,2));
U = U(:,1:K);
S = S(1:K,1:K);
V=V(:,1:K);
end

author: BatzoglouLabSU

repo: SIMLR


text:
%% housekeeping
close all
clear
clc
%% load the data

load tut01_data

%% set up VAR model

warning('This system is nonstationary because of LGDP !!!')

endog={'rtwi','LGDP','Aldp','r','LRER'};

nlags=4;

const=true;

exog={'du93Q1','du95Q4','du92Q3'};

v=rfvar(endog,exog,nlags,const); % formerly redfvar

%%
clc
ve=estimate(v,db,{db.LGDP.start,db.LGDP.finish});

%% should the Feds fund rate react to domestic variables?
linres={};

for ilag=1:nlags

    for iv=2:numel(endog)

        y=endog{iv};

        linres=[linres;{sprintf('b%0.0f(1,%s)=0',ilag,y)}];

    end

end

%%

ve_lr=estimate(v,db,{db.LGDP.start,db.LGDP.finish},[],linres);

%% save estimated models for later use

models=struct('ve',ve,'ve_lr',ve_lr);

save('tut02_estimation','models')

author: jmaih


text:
function t = bartlett_unit_sample ( m, df )

%*****************************************************************************80
%
%% BARTLETT_UNIT_SAMPLE samples the unit Bartlett distribution.
%
%  Discussion:
%
%    If the matrix T is sampled from the unit Bartlett distribution, then
%    the matrix W = T' * T is a sample from the unit Wishart distribution.
%
%    This function requires functions from the PDFLIB and RNGLIB libraries.
%
%    The "initialize()" function from RNGLIB must be called before using
%    this function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Patrick Odell, Alan Feiveson,
%    A numerical procedure to generate a sample covariance matrix,
%    Journal of the American Statistical Association,
%    Volume 61, Number 313, March 1966, pages 199-203.
%
%    Stanley Sawyer,
%    Wishart Distributions and Inverse-Wishart Sampling,
%    Washington University,
%    30 April 2007, 12 pages.
%
%  Parameters:
%
%    Input, integer M, the order of the matrix.
%
%    Input, integer DF, the number of degrees of freedom.
%    M <= DF.
%
%    Output, real T(M,M), the sample matrix from the unit Bartlett distribution.
%
  if ( df < m )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BARTLETT_UNIT_SAMPLE - Fatal error!\n' );
    fprintf ( 1, '  DF = %d < M = %d.\n', df, m );
    error ( 'BARTLETT_UNIT_SAMPLE - Fatal error!\n' );
  end

  t = zeros ( m, m );

  for i = 1 : m - 1
    df_chi = df - i + 1;
    t(i,i) = sqrt ( r8_chi_sample ( df_chi ) );
    for j = i + 1 : m
      t(i,j) = r8_normal_01_sample ( );
    end
  end

  return
end


text:
function [HotellingT2] = HotellingT2(X,alpha)
%Hotelling T-Squared testing procedures for multivariate samples.
%
%   Syntax: function [HotellingT2] = HotellingT2(X,alpha)
%
%     Inputs:
%          X - multivariate data matrix.
%      alpha - significance level (default = 0.05).
%
%     Outputs:
%          It depends of the Hotelling's T-Squared multivariate test of interest,
%          being able to be:
%
%            |-One-sample
%            |                          |-Homoskedasticity (to test)
%            |            |-Independent |
%            |            |             |-Heteroskedasticity (to test)
%            |-Two-sample |
%                         |
%                         |-Dependent
%
%          Each case calls to a corresponding function that contains a complete
%          explanation.
%
%  Created by A. Trujillo-Ortiz and R. Hernandez-Walls
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.mx
%             And the special collaboration of the post-graduate students of the 2002:2
%             Multivariate Statistics Course: Karel Castro-Morales, Alejandro Espinoza-Tenorio,
%             Andrea Guia-Ramirez.
%
%  Copyright (C) December 2002.
%

if nargin < 2,
    alpha = 0.05; %(default)
end;

if nargin < 1,
   error('Requires at least one input argument.');
end;

sam = input('Do you have one multivariate sample (1) or two multivariate samples (2)?: ');
if sam == 1;
   T2Hot1(X,alpha)
else
   id = input('They are independent (1) or dependent (2)?: ');
   if id == 1;
      disp('The covariance matrix homogeneity will be testing.:');
      MBoxtest(X,alpha);
      disp(' ')
      dc = input('Are they significant? (y/n): ','s');
      if dc == 'y'
         T2Hot2ihe(X,alpha);
      else
         T2Hot2iho(X,alpha);
      end;
   else
      T2Hot2d(X,alpha);
   end;
end;

return;

author: xiuyechen

repo: FishExplorer


text:
%SUBTRACT  Calculates the per-element difference between two arrays or array and a scalar
%
%     dst = cv.subtract(src1, src2)
%     dst = cv.subtract(src1, src2, 'OptionName',optionValue, ...)
%
% ## Input
% * __src1__ first input array or a scalar.
% * __src2__ second input array or a scalar.
%
% ## Output
% * __dst__ output array of the same size and number of channels as the input
%   array(s). The depth is defined by `DType` or that of `src1`/`src2`.
%
% ## Options
% * __Mask__ optional operation mask; this is an 8-bit single channel array
%   that specifies elements of the output array to be changed. Not set by
%   default.
% * __Dest__ Used to initialize the output `dst` when a mask is used. Not set
%   by default.
% * __DType__ optional depth of the output array: `uint8`, `int16`, `double`,
%   etc. (see the discussion below). Must be specified if input arrays are of
%   different types. default -1
%
% The function cv.subtract calculates:
%
% * Difference between two arrays, when both input arrays have the same size
%   and the same number of channels:
%
%       dst(I) = saturate(src1(I) - src2(I)) if mask(I) != 0
%
% * Difference between an array and a scalar, when `src2` is constructed from
%   Scalar or has the same number of elements as `size(src1,3)`:
%
%       dst(I) = saturate(src1(I) - src2) if mask(I) != 0
%
% * Difference between a scalar and an array, when `src1` is constructed from
%   Scalar or has the same number of elements as `size(src2,3)`:
%
%       dst(I) = saturate(src1 - src2(I)) if mask(I) != 0
%
% * The reverse difference between a scalar and an array in the case of
%   `SubRS`:
%
%       dst(I) = saturate(src2 - src1(I)) if mask(I) != 0
%
% where `I` is a multi-dimensional index of array elements. In case of
% multi-channel arrays, each channel is processed independently.
%
% The first function in the list above can be replaced with matrix expressions:
%
%     dst = src1 - src2;
%
% The input arrays and the output array can all have the same or different
% depths. For example, you can subtract to 8-bit unsigned arrays and store the
% difference in a 16-bit signed array. Depth of the output array is determined
% by `DType` parameter. In the second and third cases above, as well as in the
% first case, when `class(src1) == class(src2)`, `DType` can be set to the
% default -1. In this case the output array will have the same depth as the
% input array, be it `src1`, `src2` or both.
%
% Note: Saturation is not applied when the output array has the depth `int32`.
% You may even get result of an incorrect sign in the case of overflow.
%
% See also: cv.add, cv.addWeighted
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   demo script for surface repairing using surf2vol and remeshsurf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% preparation
% user must add the path of iso2mesh to matlab path list
% addpath('../');

% user need to add the full path to .../iso2mesh/bin directory
% to windows/Linux/Unix PATH environment variable

%% load the sample data
load rat_head.mat

% volimage is a volumetric image such as an X-ray or MRI image
% A,b are registration matrix and vector, respectively
%% perform mesh generation

[node,face]=v2s(volimage,0.5,2,'cgalmesh');

node=node(:,1:3);
face=face(:,1:3);

plotmesh(node,face);
axis equal

[newno,newfc]=remeshsurf(node,face,1);

newno=sms(newno,newfc(:,1:3),3,0.5);

figure;
plotmesh(newno,newfc(:,1:3));
axis equal

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text: %MAKEFIGS_TRIVARIATE has been moved to JLAB_MAKEFIGS.

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function name = p32_name ( )

%*****************************************************************************80
%
%% P32_NAME returns the name of problem 32.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string NAME, the name of the problem.
%
  name = 'Genz #6, Discontinuous';

  return
end


text:
function bvec_mul_test ( )

%*****************************************************************************80
%
%% BVEC_MUL_TEST tests BVEC_MUL;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2015
%
%  Author:
%
%    John Burkardt
%
  n = 15;
  seed = 123456789;
  test_num = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BVEC_MUL_TEST\n' );
  fprintf ( 1, '  BVEC_MUL multiplies binary vectors \n' );
  fprintf ( 1, '  representing integers;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        I        J        I * J   BVEC_MUL\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    [ i, seed ] = i4_uniform_ab ( -100, 100, seed );
    [ j, seed ] = i4_uniform_ab ( -100, 100, seed );

    k = i * j;

    bvec1 = i4_to_bvec ( i, n );
    bvec2 = i4_to_bvec ( j, n );
    bvec3 = bvec_mul ( n, bvec1, bvec2 );
    l = bvec_to_i4 ( n, bvec3 );

    fprintf ( 1, '  %8d  %8d  %8d  %8d\n', i, j, k, l );

  end

  return
end


text:
function [ a, b ] = p01_ab ( m )

%*****************************************************************************80
%
%% P01_AB returns bounds for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Output, real A(M,1), B(M,1), lower and upper bounds.
%
  a(1:m,1) = 0.0;
  b(1:m,1) = 1.0;

  return
end


text:
function [v,varargout] = load(fname,varargin)
% import directions
%
% Description
%
% vector3d.load is a high level method for importing vector data from
% external files. It autodetects the format of the file. As parameters the
% method requires a filename and the column positions of either the x, y, z
% coordinates or the polar angles
%
% Syntax
%   v = vector3d.load(fname,'ColumnNames',{'x','y','z'})
%   v = vector3d.load(fname,'ColumnNames',{'latitude','longitude'})
%   v = vector3d.load(fname,'ColumnNames',{'polar angle','azimuth'},'columns',[2,3])
%
% Input
%  fname - file name (text files only)
%
% Options
%  columnNames       - names of the colums to be imported
%  columns           - postions of the columns to be imported
%  radians           - treat input in radiand
%  delimiter         - delimiter between numbers
%  interface  - specific interface to be used
%
% Output
%  v - @vector3d
%
% See also
%

if iscell(fname), fname = fname{1};end

%  determine interface
if check_option(varargin,'interface')
  interface = get_option(varargin,'interface');
  options = delete_option(varargin,'interface',1);
  if isempty(interface), return; end
elseif check_option(varargin,'columnames')
  interface = 'generic';
  options = varargin;
else
  [interface,options] = check_interfaces(fname,'Vector3d',varargin{:});
end

% load tensor
[v,varargout{1:nargout-1}] = feval(['loadVector3d_',char(interface)],fname,options{:});

author: mtex-toolbox

repo: mtex

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


text:
function D = COPAR_updateD(Y, Y_range, D, X, opts) % and DCp1
% function D = COPAR_updateD(Y, Y_range, D, X, opts)
% update D in COPAR, including both PARTICULAR dictionaries and the
% COMMON dictionary.
% The algorithm used here is the efficient algorithm presented in LRSDL paper
% -----------------------------------------------
% Author: Tiep Vu, thv102@psu.edu, 5/11/2016
%         (http://www.personal.psu.edu/thv102/)
% -----------------------------------------------
    if nargin == 0
        clc;
        addpath(fullfile('..', 'utils'));
        addpath(fullfile('..', 'DLSI'));
        C = 10;    N = 7;    d = 300;
        k = 7;
        opts.k0          = 10;
        opts.lambda      = 0.01;
        opts.eta         = 0.1;
        opts.D_range     = k* (0:C);
        opts.D_range_ext = [opts.D_range opts.D_range(end)+opts.k0];
        opts.max_iter    = 10;
        opts.verbose      = true;

%         Y       = normc(rand(d, C*N));
%         Y_range = N * (0:C);
%         D       = normc(rand(d, opts.D_range_ext(end)));
%         X       = 0.01*rand(size(D,2), size(Y,2));
%         save('tmp2.mat', 'Y', 'Y_range', 'D', 'X');

        load('tmp2.mat', 'Y', 'Y_range', 'D', 'X');
    end

    %%
    C              = numel(Y_range) - 1;
    D_range_ext    = opts.D_range_ext;
    DCp1           = get_block_col(D, C+1, D_range_ext);
    optsD          = opts;
    optsD.verbose   = false;
    optsD.max_iter = 100;
    Yhat           = zeros(size(Y));
    %% ========= update Dc ==============================
    for c = 1: C
        %Dc = arg\min_Dc \|Ychat - Dc*Xcc\|_F^2 + \|Ycbar - Dc*Xcc\| + 2*eta\|A*Dc\|_F^2
        % = \arg\min_Dc \| [Ychat Ycbar] - Dc*[Xcc Xcc]\|_F^2 + 2*eta\|A*Dc\|_F^2
        % and solved using DLSI_updateD
        Dc_range = D_range_ext(c)+1: D_range_ext(c+1);
        Yc_range = Y_range(c)+1: Y_range(c+1);
        Yc       = Y(:, Yc_range);
        Dc       = D(:, Dc_range);
        Xc       = X(:, Yc_range);
        Xcc      = get_block_row(Xc, c, D_range_ext);
        XCp1c    = get_block_row(Xc, C+1, D_range_ext);
        Ychat    = Yc - D*Xc + Dc*Xcc;
        Ycbar    = Yc - DCp1*XCp1c;
        E        = (Ychat + Ycbar)*Xcc';
        F        = 2*Xcc*Xcc';
        A        = D;
        A(:,Dc_range)     = [];
        D(:, Dc_range)    = DLSI_updateD(Dc, E, F, A', opts.eta, optsD);
        Yhat(:, Yc_range) = Yc - D(:, Dc_range)*Xcc;
    end
    %% ========= DCp1 ==============================
    XCp1       = X(D_range_ext(C+1) + 1 : D_range_ext(C+2), :);
    Ybar       = Y - D(:, 1: D_range_ext(end-1))*X(1: D_range_ext(end-1), :);
    E          = (Ybar + Yhat)*XCp1';
    F          = 2*XCp1*XCp1';
    A          = D(:, 1: D_range_ext(C+1));
    DCp1_range = D_range_ext(C+1) + 1: D_range_ext(C+2);
    D(:, DCp1_range) = DLSI_updateD(D(:, DCp1_range), E, F, A', opts.eta,  optsD);
    %%
    if nargin == 0
        D = [];
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function C = xor(A,B)
%XOR Logical XOR for sptensors.
%
%   See also SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

%% Observations for sparse matrix case.
% The result of xor(a,5) is dense!
% The result of xor(a,0) is dense!
% The result of xor(a,full(a)) is dense!
% The result of xor(a,b) is sparse.

%% Case 1: One argument is a scalar
if isscalar(B) || isa(B,'tensor')
    C = xor(full(A),B);
    return;
end
if isscalar(A)
    C = xor(A,full(B));
    return;
end

%% Case 2: Both x and y are tensors of some sort
if ~isequal(size(A),size(B))
    error('Must be tensors of the same size');
end

if isa(A,'sptensor') && isa(B,'sptensor')
    C = sptensor([A.subs; B.subs], 1, size(A), @(x) length(x) == 1);
    return;
end

if isa(B,'tensor')
    Bsubs = find(B ~= 0);
    C = sptensor([A.subs; Bsubs], 1, size(A), @(x) length(x) == 1);
    return;
end

%% Otherwise
error('The arguments must be two sptensors or an sptensor and a scalar.');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% MEDIAN Median value.
%    For vectors, MEDIAN(x) is the median value of the elements in x.
%    For matrices, MEDIAN(X) is a row vector containing the median value
%    of each column.  For N-D arrays, MEDIAN(X) is the median value of the
%    elements along the first non-singleton dimension of X.
%
%    MEDIAN(X,DIM) takes the median along the dimension DIM of X.
%
%    MEDIAN(...,NANFLAG) specifies how NaN (Not-A-Number) values
%    are treated. The default is 'includenan':
%
%    'includenan' - the median of a vector containing NaN values is also NaN.
%    'omitnan'    - the median of a vector containing NaN values is the
%                   median of all its non-NaN elements. If all elements
%                   are NaN, the result is NaN.
%
%    Example:
%        X = [1 2 4 4; 3 4 6 6; 5 6 8 8; 5 6 8 8]
%        median(X,1)
%        median(X,2)
%
%    Class support for input X:
%       float: double, single
%       integer: uint8, int8, uint16, int16, uint32, int32, uint64, int64
%
%    See also MEAN, STD, MIN, MAX, VAR, COV, MODE.
%
%    Reference page in Doc Center
%       doc median
%
%    Other functions named median
%
%       categorical/median    duration/median    timeseries/median
%       datetime/median       tall/median        ts/median
%       distributed/median
%


text:
function [wpos,ypos]=rlasymp(G)
G=zpk(G);
zer=G.z{1}; pol=G.p{1}; gain=G.k;
ii=find(abs(zer)<1e10);
if length(ii)>0, zer=zer(ii); end
nExcess=length(pol)-length(zer);
if nExcess>0
   pp=(sum(real(pol))-sum(real(zer)))/nExcess;
   deltaP=pi/nExcess;
   xx=get(gca,'Xlim');
   yy=get(gca,'YLim');
   xx1=(xx(1)-pp)*tan(deltaP);
end
wpos=[pp*ones(1,nExcess); zeros(1,nExcess)];
ypos=zeros(2,nExcess);
for i=1:nExcess
   PAngle=(2*i-1)*deltaP;
   Kslp=tan(PAngle);
   if (pi/2>=PAngle & PAngle>=0)
      xP=xx(2); yP=yy(2);
   elseif (pi>=PAngle&PAngle>pi/2)
      xP=xx(1); yP=yy(2);
   elseif (3*pi/2>=PAngle&PAngle>pi)
      xP=xx(1); yP=yy(1);
   else
      xP=xx(2); yP=yy(1);
   end
   xx1=xP; yy1=Kslp*(xx1-pp);
   if yy1>yy(2) | yy1<yy(1),
      yy1=yP; xx1=yy1/Kslp+pp;
   end
   wpos(2,i)=xx1; ypos(2,i)=yy1;
end


text:
function [f,k,s0,gof] = relaxFitMtDist(data,delta,t1,s0,tr,flipAngle,brainMask,fitMethod,outMontage)
%
% [f,k,s0,gof] = relaxFitMtDist(data,delta,t1,s0,tr,flipAngle,brainMask,fitMethod,outMontage)
%
% Computes a nonlinear fit of the bound-pool (f) map.
%
% Returns:
%
% SEE ALSO:
%
% relaxFitT1.m to fit the t1 and pd maps used by this function.
%
% HISTORY:
% 2008.02.26 RFD: wrote it.

if(~exist('fitMethod','var')||isempty(fitMethod))
    fitMethod = 'fmin';
end
if(lower(fitMethod(1))=='p')
    fitMethod = fitMethod(2:end);
    useParfor = true;
else
    useParfor = false;
end
fitMethod = lower(fitMethod(1));

if(~exist('outMontage','var')||isempty(outMontage))
    outMontage = 'f.png';
end

flipAngle = flipAngle*pi/180;
tr = tr/1000;
sz = size(data);

% disp('smoothing s0 map...');
% deltaT = 1/44;
% kappa = 200;
% s0 = dtiSmoothAnisoPM(s0, 4, deltaT, kappa, 1, [1 1 1]);
%
%disp('smoothing MT measurments...');
%kappa = 50;
%for(ii=1:sz(4))
%    data(:,:,:,ii) = dtiSmoothAnisoPM(data(:,:,:,ii), 2, deltaT, kappa, 1, [1 1 1]);
%end

brainInds = find(brainMask);
numVoxelsPerUpdate = 16000;
nVoxAll = length(brainInds);
for(ii=1:size(data,4))
  tmpVol = data(:,:,:,ii);
  tmpMT(ii,:) = tmpVol(brainInds);
end
clear tmpVol;
data = tmpMT;
clear tmpMT;
r1 = 1./t1(brainInds);
s0 = s0(brainInds);

f = zeros(1,nVoxAll);
k = zeros(1,nVoxAll);
gof = zeros(1,nVoxAll);
totalSecs = 0;

tmpImg = zeros(sz(1:3));
tmpName = tempname
nSteps = ceil(nVoxAll/numVoxelsPerUpdate);
totalTime = 0;
tic;
for(ii=1:nSteps)
    if(useParfor), matlabpool; end
    curInd = (ii-1)*numVoxelsPerUpdate+1;
    endInd = min(curInd+numVoxelsPerUpdate,nVoxAll);
    [tf,tk,tgof] = relaxFitMt(data(:,curInd:endInd),delta,r1(curInd:endInd),s0(curInd:endInd),tr,flipAngle,fitMethod,useParfor);
    prevSecs = toc;
    totalTime = totalTime+prevSecs;
    f(curInd:endInd) = tf;
    k(curInd:endInd) = tk;
    gof(curInd:endInd) = tgof;
    secsPerVox = prevSecs/(endInd-curInd+1);
    estTime = secsPerVox*(nVoxAll-endInd);
    if(estTime>5400) estTime=estTime./3600; estTimeUnits='hours';
    elseif(estTime>90) estTime=estTime./60; estTimeUnits='minutes';
    else estTimeUnits='seconds';
    end
    fprintf('Processed %d of %d voxels- %0.1f %s remaining (%0.3f secs per vox)...\n',endInd,nVoxAll,estTime,estTimeUnits,secsPerVox);
    tmpImg(brainInds) = f;
    m = makeMontage(tmpImg);
    if(max(f)>0), m = m./max(f); end
    m = uint8(round(m.*255));
    imwrite(m,outMontage);
    save(tmpName,'k','f','gof','brainMask');
    tic;
    if(useParfor), matlabpool close; end
end
fprintf('Processed %d voxels in %0.2f hours.\n',totalTime/3600);

% f = vertcat(tf);
% k = vertcat(tk);
% gof = vertcat(tgof);

im=zeros(size(brainMask)); im(brainInds) = f; f = im;
im=zeros(size(brainMask)); im(brainInds) = k; k = im;
im=zeros(size(brainMask)); im(brainInds) = gof; gof = im;
im=zeros(size(brainMask)); im(brainInds) = s0; s0 = im;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function M1 = spm_eeg_inv_headcoordinates(nas, lpa, rpa)
% Returns the homogenous coordinate transformation matrix
% that converts the specified fiducials in any coordinate system (e.g. MRI)
% into the rotated and translated headccordinate system.
%
% M1 = headcoordinates(nas, lpa, rpa)
%
% The headcoordinate system in CTF is defined as follows:
% the origin is exactly between lpa and rpa
% the X-axis goes towards nas
% the Y-axis goes approximately towards lpa, orthogonal to X and in the plane spanned by the fiducials
% the Z-axis goes approximately towards the vertex, orthogonal to X and Y
%_______________________________________________________________________
% Copyright (C) 2003 Robert Oostenveld

% Robert Oostenveld
% $Id: spm_eeg_inv_headcoordinates.m 3589 2009-11-20 17:17:41Z guillaume $

% ensure that they are row vectors
lpa = lpa(:)';
rpa = rpa(:)';
nas = nas(:)';

% compute the origin and direction of the coordinate axes in MRI coordinates

% follow CTF convention
origin = [lpa+rpa]/2;
dirx = nas-origin;
dirx = dirx/norm(dirx);
dirz = cross(dirx,lpa-rpa);
dirz = dirz/norm(dirz);
diry = cross(dirz,dirx);

% compute the rotation matrix
rot = eye(4);
rot(1:3,1:3) = inv(eye(3) / [dirx; diry; dirz]);
% compute the translation matrix
tra = eye(4);
tra(1:4,4)   = [-origin(:); 1];
% compute the full homogenous transformation matrix from these two
M1 = rot * tra;


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%                  1         , -2<=t<=2
% 	Graph of f(t)= 0         , 2<t<5
%                t*sin(4pi*t), 5<=t<=8

t1=-2:.1:2;
t2=2.1:.1:4.9;
t3=5:.1:8;

f1=ones(size(t1));
f2=zeros(size(t2));
f3=t3.*sin(4*pi*t3);

t=[t1 t2 t3];
f=[f1 f2 f3];

plot(t,f)
title('Multi-part function f(t)')

author: Sable


text:
% TESTFITLINE - demonstrates RANSAC line fitting
%
% Usage: testfit2Dline(outliers, sigma, t, feedback)
%
% Arguments:
%               outliers - Fraction specifying how many points are to be
%                          outliers.
%               sigma    - Standard deviation of inlying points from the
%                          true line.
%               t        - Distance threshold to be used by the RANSAC
%                          algorithm for deciding whether a point is an
%                          inlier.
%               feedback - Optional flag 0 or 1 to turn on RANSAC feedback
%                          information.
%
%  Try using:  testfit2Dline(0.3, 0.05, 0.05)
%
% See also: RANSACFITPLANE, FITPLANE

% Copyright (c) 2003-2006 Peter Kovesi and Felix Duvallet (CMU)
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% 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, 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.

% August 2006  testfitline created from testfitplane
%              author: Felix Duvallet

function testfit3Dline(outliers, sigma, t, feedback)

    close all;
    if nargin == 0
        outliers = 0.3;
        sigma = 0.05;
        t = 0.05;
        feedback = 1;
    end
    if nargin == 3
        feedback = 0;
    end

    % Hard wire some constants - vary these as you wish

    npts = 100;  % Number of 3D data points

    % Define a line:
    %    Y = m*X
    %    Z = n*X + Y + b
    % This definition needs fixing, but it works for now

    m = 6;
    n = -3;
    b = -4;

    outsigma = 30*sigma;  % outlying points have a distribution that is
                          % 30 times as spread as the inlying points

    vpts = round((1-outliers)*npts);  % No of valid points
    opts = npts - vpts;               % No of outlying points

    % Generate npts points in the line
    X = rand(1,npts);

    Y = m*X + b;

    XY =  [X; Y];

    % Add uniform noise of +/-sigma
    XY = XY + (2*rand(size(XY))-1)*sigma;

    % Generate opts random outliers

    n = length(XY);
    ind = randperm(n);  % get a random set of point indices
    ind = ind(1:opts);  % ... of length opts

    % Add uniform noise of outsigma to the points chosen to be outliers.
    XY(:,ind) = XY(:,ind)  +   sign(rand(2,opts)-.5).*(rand(2,opts)+1)*outsigma;

    % Perform RANSAC fitting of the line
    [V, P, inliers] = ransacfit2Dline(XY, t, feedback);

    if(feedback)
        disp(['Number of Inliers: ' num2str(length(inliers)) ]);
    end

    % We want to plot the inlier points blue, with the outlier points in
    % red.  In order to do that, we must find the outliers.
    % Use setxor on all the points, and the inliers to find outliers
    %  (plotting all the points in red and then plotting over them in blue
    %  does not work well)
    oulier_points = setxor(transpose(XY), transpose(XY(:, inliers)), 'rows');
    oulier_points = oulier_points';

    % Display the cloud of outlier points
    figure(1); clf
    hold on;
    plot(oulier_points(1,:),oulier_points(2,:) , 'r*');

    % Plot the inliers as blue points
    plot(XY(1,inliers), XY(2, inliers), 'b*');

    % Display the line formed by the 2 points that gave the
    % line of maximum consensus as a green line
    line(P(1,:), P(2,:), 'Color', 'green', 'LineWidth', 4);

    %Display the line formed by the covariance fitting in magenta
    line(V(1,:), V(2, :), 'Color', 'magenta', 'LineWidth', 5);
    grid('on');
    % box('on');

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
% *************************************************************************
% Video Super-Resolution with Convolutional Neural Networks
%
% main VSRnet script, includes the interface to the caffe framework
% do not run this script directly. call it from VSRnet_demo instead
%
% Version 1.0
%
% Created by:   Armin Kappeler
% Date:         02/19/2016
%
% *************************************************************************

%% more parameters (do not change, unless you know what you are doing)

BORDERSIZE = 8; %4+2+2; %border size = sum of all zeropaddings in model def file

% parameters for LOW_MEMORY_MODE=2
Param.patchSize = 36;
Param.stride = 20;
Param.outputPatchSize = 20;
Param.zeropadding = 1;

%% Init VSRnet configuration files

if (MOTIONCOMPENSATION==false && UPSCALE_FACTOR ~= 3 && ADAPTIVEMOTIONCOMPENSATION==false)
    error('Error: MOTIONCOMPENSATION=false is only available for upscale factor 3. Please either set UPSCALEFACTOR=3 or MOTIONCOMPENSATION=true')
end
if (MOTIONCOMPENSATION==false && ADAPTIVEMOTIONCOMPENSATION==true)
    display('The MOTIONCOMPENSATION flag is ignored because ADAPTIVEMOTIONCOMPENSATION=true')
end

if MOTIONCOMPENSATION
    model_def_file_template = ['models/u' num2str(UPSCALE_FACTOR) '/superres_deploy.prototxt'];
    model_file =    ['models/u' num2str(UPSCALE_FACTOR) '/SRnet_train_iter_200000.caffemodel'];
else
    model_def_file_template = ['models/u' num2str(UPSCALE_FACTOR) '_noMotionCompensation/superres_deploy.prototxt'];
    model_file =    ['models/u' num2str(UPSCALE_FACTOR) '_noMotionCompensation/SRnet_train_iter_200000.caffemodel'];
end
model_def_file = 'models/tmp/tmp_superres_deploy.prototxt';

%% preprocess input data

TESTVIDEO_PATH
if PREPROCESSED_INPUT
    load(TESTVIDEO_PATH);

    if TESTONLY1FRAME
        input_data{1} = input_data{1}(:,:,:,TESTONLY1FRAME);
        im_gt = im_gt(:,:,:,TESTONLY1FRAME);
    end
else
    data = load(TESTVIDEO_PATH);
    if TESTONLY1FRAME
        data.frames = data.frames(:,:,TESTONLY1FRAME:TESTONLY1FRAME+4);
    end
    [input_data,idx_gt] = preprocess_frames(data.frames,UPSCALE_FACTOR,MOTIONCOMPENSATION,ADAPTIVEMOTIONCOMPENSATION);

end
im_bic = permute(input_data{1}(:,:,3,:),[2,1,3,4])*255;

%% caffe setup

batchsize = size(input_data{1},4);
imageSize = [ size(input_data{1},2), size(input_data{1},1)];
Param.inputImgSize = imageSize;

caffe.reset_all();

if USE_GPU
  caffe.set_mode_gpu();
  caffe.set_device(GPU_ID);
else
  caffe.set_mode_cpu();
end

if LOW_MEMORY_MODE==1
    imagecut = imageSize(1)/2;
    imageSize(1) = imageSize(1)/2 + BORDERSIZE;
elseif LOW_MEMORY_MODE==2
    imageSize = [Param.patchSize Param.patchSize];
end

change_caffe_image_input_size(model_def_file_template,model_def_file,imageSize,1)
net = caffe.Net(model_def_file, model_file, 'test');

%% do forward pass

tic;
output = zeros(size(im_bic,2),size(im_bic,1),size(im_bic,3),size(im_bic,4));
for imgIdx = 1:size(im_bic,4)
    display(['Superresolve frame ' num2str(imgIdx)]);
    if LOW_MEMORY_MODE==1 % divide image in two parts to reduce memory usage
        tmp1 = net.forward({input_data{1}(:,1:imagecut + BORDERSIZE,:,imgIdx)});
        tmp2 = net.forward({input_data{1}(:,imagecut+1-BORDERSIZE:end,:,imgIdx)});
        output(:,:,:,imgIdx) = cat(2,tmp1{1}(:,1:imagecut,:,:),tmp2{1}(:,BORDERSIZE+1:end,:,:));
    elseif LOW_MEMORY_MODE==2 % process the image patchwise
        input_data_new = patchify(input_data{1}(:,:,:,imgIdx), Param);
        if Param.zeropadding==1
            output_raw = zeros(size(input_data_new,1),size(input_data_new,2),1,size(input_data_new,4));
        else
            output_raw = zeros(Param.outputPatchSize,Param.outputPatchSize,1,size(input_data_new,4));
        end
        for pIdx = 1:size(input_data_new,4)
            output_raw_tmp = net.forward({input_data_new(:,:,:,pIdx)});
            output_raw(:,:,:,pIdx) = output_raw_tmp{1};
            if mod(pIdx,1000)==0
                fprintf('.');
            end
        end
        fprintf('\n');
        output(:,:,:,imgIdx) = depatchify(output_raw,Param);
    else
        tmp = net.forward({input_data{1}(:,:,:,imgIdx)});
        output(:,:,:,imgIdx) = tmp{1};
    end
end
toc;

im_SR = output;
im_SR = permute(im_SR,[2 1 3 4]);
im_SR = im_SR*255;

save(RESULTFILE_PATH,'im_SR','model_file','TESTVIDEO_PATH');

author: thomas-koehler

repo: SupER


text:
function map = readCbMap(fileName)
% readCbMap reads in a map from a textfile and returns a map struct
%
% map = readCbMap(fileName)
%
%INPUT
% fileName is the text file exported from BiGG database
%
%OUTPUT
% map                   Map structure
%   molPosition             the positions of the metabolite nodes
%   molIndex                index for molecule node(same as connection index)
%   molName                 full name of metabiolites
%   molAbbreviation         The abbreviation formetabolites(labels)
%   molLabelPos             The positions to display the metabolite's labels
%   molPrime                Y if metabolite is primary and N if not
%   rxnPosition             the positions of the reaction nodes
%   rxnLabelPosition        the positions of the reaction nodes' labels if the
%                           node is a midpoint
%   rxnIndex                the index of reaction nodes used in connection
%   connection              the node connection matrix
%   connectionAbb           the reaction abbreviations assigned to each
%                           segment
%   connectionName          the name of reactions
%   connectionReversibile   1 if reaction is reversible and 0 otherwise
%

% The general format of the import file in 9 cols
if nargin < 1
    [fileName,filePath] = uigetfile({'*.txt'});
else
    filePath = '';
    if isempty(regexp(fileName,'.txt$', 'once'))
        fileName = strcat(fileName,'.txt');
    end
end

format = '%s %n %s %f %f %f %f %n %s';
fid = fopen(strcat(filePath,fileName));
molecules = textscan(fid, format,'delimiter', '\t','HeaderLines',1);
fclose(fid);

molAbb = molecules{1,1};
molComp = molecules{1,2};
molPrim = molecules{1,3};
molLabelx = molecules{1,4};
molLabely = molecules{1,5};
molPosx = molecules{1,6};
molPosy = molecules{1,7};
molName = molecules{1,9};
molId = molecules{1,8};
% s = size(molLabelx);            %lenght of a general column
% i = 1;

% while i<= s(1)
%     if strcmp(molAbb(i),'Reactions Nodes')
%         rxnN = i;               % number of molecules
%     elseif strcmp(molAbb(i),'Reactions')
%         rxnEnd = i;
%     elseif strcmp(molAbb(i),'Texts')
%         textBeg = i;
%     end
%     i= i+1;
% end

rxnN = (strmatch('Reactions Nodes',molAbb));
rxnEnd = strmatch('Reactions',molAbb,'exact');
textBeg = strmatch('Texts',molAbb);

indexMols = (1:(rxnN-1));
indexRxnNodes = ((rxnN+1):(rxnEnd-1));
indexRxnNames = ((rxnEnd+1):(textBeg-1));
indexTexts = ((textBeg+1):length(molAbb));

rxnNum = (rxnEnd-rxnN-1);       % number of rxn nodes
connectionNum = textBeg-rxnEnd;    % number of connections: connectionNum

% set up the position matrix for molecules
molPos(:,1) = molPosx(indexMols);
molPos(:,2) = molPosy(indexMols);

% set up the position matrix for molecules' labels
molLabelPos = [molLabelx(indexMols) molLabely(indexMols)];

% set up the position matrix for reactions
%rxnPos = [molLabelx(indexRxnNodes) molLabely(indexRxnNodes)];
%Updated BiGG to output files in a less disorganized fashion
rxnPos = [molPosx(indexRxnNodes) molPosy(indexRxnNodes)];

% set up the position matrix for reactions' labels
%rxnLabelPos = [molPosx(indexRxnNodes) molPosy(indexRxnNodes)];
%Updated BiGG to output files in a less disorganized fashion
rxnLabelPos = [molLabelx(indexRxnNodes) molLabely(indexRxnNodes)];

% reaction ids
rxnId = molId(indexRxnNodes);

% set up the connection matrix
node = [molLabelx(indexRxnNames) molLabely(indexRxnNames)];

% set up the reversibility for each connection
r = zeros(connectionNum,1);
r(strmatch('Reversible',molPrim(indexRxnNames),'exact'))=1;

rxnAbb = molAbb(indexRxnNames);
rxnName = molName(indexRxnNames);

%% Handling other shapes and texts
shapeIDs = [strmatch('Circle',molAbb); strmatch('Rect',molAbb); strmatch('Line',molAbb)];
indexShapes = logical(sparse(length(molAbb),1));
indexShapes(shapeIDs) = 1;
map.shapeType = molAbb(indexShapes);
map.shapePos = [molLabelx(indexShapes) molLabely(indexShapes)];
map.shapeSize = [molPosx(indexShapes) molPosy(indexShapes)];
for i=1:length(shapeIDs)
    str = molName(shapeIDs(i));
    [colorStr, str] = strtok(str,':');
    [c1,c2] = strtok(colorStr,'/');
    [c2,c3] = strtok(c2,'/');
    c3 = regexprep(c3,'\/','');
    if isnan(str2double(c3))
        c3 = 0;
    end
    [thickness, style] = strtok(str,'@');
    thickness = regexprep(thickness,'(\w*):','');
    style = regexprep(style,'\@','');
    map.shapeStyle(i,1) = style;
    map.shapeThickness(i,1) = str2num(thickness{1,1});
    map.shapeColor(i,1:3) = [str2double(c1{1,1}) str2double(c2{1,1}) str2double(c3{1,1})];
end
indexShapes(1:textBeg) = 1;
map.text = molName(~indexShapes);
map.textFont = molAbb(~indexShapes);
map.textPos = [molLabelx(~indexShapes) molLabely(~indexShapes)];
map.textSize = molPosx(~indexShapes);

map.molPosition = molPos';
map.molIndex = molId(1:(rxnN-1));
map.molName = molName(1:(rxnN-1));
map.molAbbreviation = molAbb(1:(rxnN-1));
map.molLabelPos = molLabelPos;
map.molPrime = molPrim(1:(rxnN-1));
map.molCompartment = molComp(1:(rxnN-1));
map.rxnPosition = rxnPos';
map.rxnLabelPosition = rxnLabelPos';
map.rxnIndex = rxnId;
map.connection = node;
map.connectionAbb = rxnAbb;
map.connectionName = rxnName;
map.connectionReversible = r;
end

author: opencobra

repo: cobratoolbox


text:
function Y = DGradient(X, Dim, Spacing, Method)  %#ok<STOUT,INUSD>
% Gradient along a dimension
% Y = DGradient(X, Dim, Spacing, Method)
% INPUT:
%   X:   Real DOUBLE array.
%   Spacing: Scalar or vector of the length SIZE(X, Dim).
%        A scalar value is the distance between all points, while a vector
%        contains all time points, such that DIFF(Spacing) are the distances.
%        For equally spaced input a scalar Spacing is much faster.
%        Optional, default: 1.0
%   Dim: Dimension to operate on.
%        Optional, default: [] (1st non-singelton dimension).
%   Method: String, order of the applied method for unevenly spaced X:
%        '1stOrder', faster centered differences as in Matlab's GRADIENT.
%        '2ndOrder', 2nd order accurate centered differences.
%        On the edges forward and backward difference are used.
%        Optional, default: '1stOrder'.
%
% OUTPUT:
%   Y:   Gradient of X, same size as X.
%
% EXAMPLES:
%   t = cumsum(rand(1, 100)) + 0.01;  t = 2*pi * t ./ max(t);
%   x = sin(t);
%   dx1 = DGradient(x, t, 2, '1stOrder');
%   dx2 = DGradient(x, t, 2, '2ndOrder');
%   dx  = cos(t);          % Analytic solution
%   h = plot(t, dx, t, dx1, 'or', t, dx2, 'og');  axis('tight');
%   title('cos(x) and DGradient(sin(x))');
%   legend(h, {'analytic', '1st order', '2nd order'}, 'location', 'best');
%
% NOTES:
% - There are a lot of other derivation tools in the FEX. This function is
%   faster, e.g. 25% faster than dqdt and 10 to 16 times faster than Matlab's
%   GRADIENT. In addition it works with multi-dim arrays, on a specific
%   dimension only, and can use a 2nd order method for unevenly spaced data.
% - This function does not use temporary memory for evenly spaced data and if
%   a single vector is processed. Otherwise the 1st-order method needs one and
%   the 2nd-order method 3 temporary vectors of the length of the processed
%   dimension.
% - Matlab's GRADIENT processes all dimensions ever, while DGradient operates on
%   the specified dimension only.
% - 1st order centered difference:
%     y(i) = (x(i+1) - x(i-1) / (s(i+1) - s(i-1))
% - 2nd order centered difference:
%     y(i) = ((x(i+1) * (s(i)-s(i-1)) / (s(i+1)-s(i))) -
%             (x(i-1) * (s(i+1)-s(i)) / (s(i)-s(i-1)))) / (s(i+1)-s(i-1))
%            + x(i) * (1.0 / (s(i)-s(i-1)) - 1.0 / (s(i+1)-s(i)))
%   For evenly spaced X, both methods reply equal values.
%
% COMPILE:
%   mex -O DGradient.c
% Consider C99 comments on Linux:
%   mex -O CFLAGS="\$CFLAGS -std=c99" DGradient.c
% Pre-compiled Mex: http://www.n-simon.de/mex
% Run the unit test uTest_DGradient after compiling.
%
% Tested: Matlab 6.5, 7.7, 7.8, WinXP, 32bit
%         Compiler: LCC2.4/3.8, BCC5.5, OWC1.8, MSVC2008
% Assumed Compatibility: higher Matlab versions, Mac, Linux, 64bit
% Author: Jan Simon, Heidelberg, (C) 2011 matlab.THISYEAR(a)nMINUSsimon.de
%
% See also GRADIENT, DIFF.
% FEX: central_diff (#12 Robert A. Canfield)
%      derivative (#28920, Scott McKinney)
%      movingslope (#16997, John D'Errico)
%      diffxy (#29312, Darren Rowland)
%      dqdt (#11965, Geoff Wawrzyniak)

% $JRev: R0c V:003 Sum:RkNcmbCGEJGT Date:02-Jan-2010 02:41:46 $
% $License: BSD $
% $File: Tools\Mex\Source\DGradient.c $
% History:
% 001: 30-Dec-2010 22:42, First version published under BSD license.

% This is a dummy file to support Matlab's HELP.

error(['JSimon:', mfilename, ':NoMex'], ...
   'Cannot find compiled Mex.\nPlease compile at first:  mex -O DGradient.c');

repo: mcbench-benchmarks


text:
function center = polylineCentroid(varargin)
%POLYLINECENTROID Computes the centroid of a curve defined by a series of points.
%
%   PT = polylineCentroid(POINTS);
%   Computes center of mass of a polyline defined by POINTS. POINTS is a
%   N-by-D array of double, representing a set of N points in a
%   D-dimensional space.
%
%   PT = polylineCentroid(PTX, PTY);
%   PT = polylineCentroid(PTX, PTY, PTZ);
%   Specifies points as separate column vectors
%
%   PT = polylineCentroid(..., TYPE);
%   Specifies if the last point is connected to the first one. TYPE can be
%   either 'closed' or 'open'.
%
%   Example
%   poly = [0 0;10 0;10 10;20 10];
%   polylineCentroid(poly)
%   ans =
%       [10 5]
%
%   See also
%     polygons2d, centroid, polygonCentroid, polylineLength
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2006-05-22
% Copyright 2006-2022 INRA - TPV URPOI - BIA IMASTE

%% process input arguments

% check whether the curve is closed
closed = false;
var = varargin{end};
if ischar(var)
    if strcmpi(var, 'closed')
        closed = true;
    end
    % remove last argument
    varargin(end) = [];
end

% extract point coordinates
if length(varargin)==1
    points = varargin{1};
elseif length(varargin)==2
    points = [varargin{1} varargin{2}];
end

%% Main computation

% compute centers and lengths composing the curve
if closed
    centers = (points + points([2:end 1],:))/2;
    lengths = sqrt(sum(diff(points([1:end 1],:)).^2, 2));
else
    centers = (points(1:end-1,:) + points(2:end,:))/2;
    lengths = sqrt(sum(diff(points).^2, 2));
end

% centroid of edge centers weighted by edge length
%weigths = repmat(lengths/sum(lengths), [1 size(points, 2)]);
center = sum(centers .* repmat(lengths, [1 size(points, 2)]), 1) / sum(lengths);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ D ] = D06( f,h1,h2 )

    Df = 0.000389*f*h1*h2;
    Dh = 4.1*(sqrt(h1)+sqrt(h2));
    D = Df*Dh/(Df+Dh);
end


text:
function data_num = p01_data_num ( )

%*****************************************************************************80
%
%% P01_DATA_NUM returns the dimension of the data vector for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DATA_NUM, the dimension of the data vector.
%
  data_num = 10;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef (Abstract) AbstractTrajectoryViewType < matlab.mixin.SetGet & matlab.mixin.Heterogeneous
    %AbstractTrajectoryViewType Summary of this class goes here
    %   Detailed explanation goes here

    methods
        [hCBodySurf, childrenHGs] = plotStateLog(obj, orbitNumToPlot, lvdData, handles)
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function peakIdxV = findFirstPeak(X)
% AI  10/17/16
% ========================================================================
% X  : Input data matrix (nVox x nTimePts)
% ========================================================================

nVox = size(X,1);

% Find pts of local maxima 't' where x(t)>x(t-1) & x(t)>x(t+1)
test1_M = [zeros(nVox,1) diff(X,1,2)];
test2_M = [-diff(X,1,2) zeros(nVox,1)];
test12_M = test1_M>=0 & test2_M>=0;
% Retain local maxima that are at least 80% of the max. signal intensity
test3_M = bsxfun(@gt,X,.8*max(X,[],2));
allPeaksM = test12_M & test3_M;
% First peak
[~,peakIdxV] = max(allPeaksM,[],2);   %Max returns the index corresponding to the
                                      %first occurrence of maximum (here 1)
                                      %If no peaks are found, returns first point)

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR

path: github-repos/MATLAB/cerr-CERR/CERR-d320754abad9dcb78508ab69f33ae9f644202114/CERR_core/PlanAnalysis/DCE-MR analysis/TTHP analysis for DCE-MR/findFirstPeak.m


text:
function varargout = prtUtilSubplotTitle(varargin)
%   prtUtilSubplotTitle add titles to subplots
%
%   prtUtilSubplotTitle(STRING) makes string STRING the title*(see note) across all
%   available sub-plots in the current figure.
%
%   prtUtilSubplotTitle(AXESLIST,STRING) place the title STRING over only the
%   axes specified in AXESLIST (a list of axes handles as generated by
%   subplot or equivalent).
%
%   For example:
%
%       clear all; close all;
%       a1 = subplot(3,3,1); plot(1:10);
%       a1_2 = subplot(3,3,2); plot(sin(linspace(0,2*pi)*5));
%       a1_3 = subplot(3,3,3); plot(cos(linspace(0,2*pi)*5));
%       a2 = subplot(3,3,[4,7]); h = surf(peaks); set(h,'EdgeAlpha',0);
%       a3 = subplot(3,3,[5,8]); h = surf(-peaks); set(h,'EdgeAlpha',0);
%       a4 = subplot(3,3,[6,9]); mesh(randn(10));
%
%       th1 = subplot_title(a1,'Plot of 1-10');
%       th2 = subplot_title([a1_2,a1_3],'Plot of Sinusoids');
%       th3 = subplot_title([a2 a3],'SURF Pics');
%       th4 = subplot_title([a4],'Mesh');
%
%   *Note that subplot_title creates invisible axes to make these new
%   titles, the handles to which are lost after subplot_title exits.  As a
%   result, multiple calls to subplot_title in a row over the same set of
%   axes will result in multiple overlaying titles.  This is undesirable.
%   To mediate these effects, use title_handle to erase current titles
%   prior to calling subplot_title multiple times.

% Originally written as subplot_title, part of FROI:
%   This function is part of froi, available from http://froi.sourceforge.net,
%   and is governed under the terms of the Artistic License.

%Mod History:
ax = gca;
cfig = gcf;

if nargin == 1
    title_string = varargin{1};
    relevantAxes = findobj(gcf, 'Type', 'axes', 'Box', 'on');
elseif nargin == 2
    relevantAxes = varargin{1};
    title_string = varargin{2};
elseif nargin > 2
    if ischar(varargin{1}) || iscell(varargin{1})
        title_string = varargin{1};
        relevantAxes = findobj(gcf, 'Type', 'axes', 'Box', 'on');
    else
        relevantAxes = varargin{1};
        title_string = varargin{2};
    end
end

if isempty(relevantAxes);
    error('No current axes available for titling');
end

for i = 1:length(relevantAxes);
    P(i,:) = get(relevantAxes(i),'Position');
end

minX = min(P(:,1));
maxX = max(P(:,3) + P(:,1));
minY = min(P(:,2));
maxY = max(P(:,4) + P(:,2));

axes('position',[minX minY maxX-minX maxY-minY],'Box','off','Visible','off');

if nargin > 2
    if ~any(ishandle(varargin{1}))
        title_handle = title(title_string,varargin{2:end});
    else
        title_handle = title(title_string,varargin{3:end});
    end
else
    title_handle = title(title_string);
end

set(title_handle,'visible','on');

children = get(gcf,'Children');
set(gcf,'Children',[children(2:end); children(1)])

figure(cfig);

if nargout>0
    varargout{1} = title_handle;
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT