text:
function foreground = extractForeground(video,frameSkip,dataRange,minSigmaLevel,shadowLevel,significanceThreshold,useGraph,graphAlpha,morphStructureElement,returnLargestComponent)

% extractForeground:  Extract the foreground from a video.
%
% This file shows sample usage for extracting the moving foreground
% from a static background.  It is intended as an example only;
% More sophisitcated background modeling techniques, particularly
% adaptive ones, may be applied instead.  However, this should
% suffice for many simple videos.
%
% A number of parameters are available to control the behavior
% of the algorithm.  If these are unsupplied or empty, reasonable
% default values will be chosen.
%
% video:  This is the only required parameter.  It may consist
%   of a string containing the name of an avi file to load.
%   Alternately, it may contain the video frames as a cell
%   array of rgb images, stored as uint8 data.
%
% frameSkip:  Long videos may use enormous amounts of memory
%   during the intermediate processing stages.  If this becomes
%   a problem, frameSkip allows the background model to be
%   computed using every nth frame.
%
% dataRange:  A cheap way of approximating robust statistics.
%   Mean/deviation are computed only for data falling within
%   the specified percentile ranges.  For example, with the
%   default value of [.25 .75], outlier pixels in up to 25% of
%   the video frames will be ignored.
%
% minSigmaLevel:  Pixels with very low (or zero) deviation cause
%   trouble.  This parameter specifies a percentile level,
%   which gives a floor deviation value for all pixels.
%
% shadowLevel:  Amount of shadowing to account for.  Positive
%   values allow some darkening of the pixel intensity without
%   treating it as evidence of a foreground object.
%
% significanceThreshold:  Number of standard deviations required
%   for a pixel to be considered foreground, in isolation.
%   This is perhaps the single most important parameter.  If pieces
%   of the foreground are being left out, try decreasing its value.
%   If too much background is included, try increasing its value.
%   (Lax dataRange bounds can also cause pieces of foreground to be
%   left out, if an object obscures the background for too long.
%   Increasing the shadowLevel can also help prevent background
%   inclusion, particularly around subjects' legs.)
%
% useGraph:  Boolean control parameter.  If true, use graph cut
%   to find foreground.  If false, use morphology.
%
% graphAlpha:  Parameter controlling "clumpiness" of graph cut
%   operation.
%
% morphStructureElement:  Parameter controlling the morphological
%   operations if these are used.
%
% returnLargestComponent:  Boolean control parameter.  If true,
%   then return single largest connected foreground component.
%   If false, return raw labels.
%
% Usage:
%   function foreground = extractForeground(video,frameSkip,dataRange,minSigmaLevel,shadowLevel,significanceThreshold,useGraph,graphAlpha,morphStructureElement,returnLargestComponent)

% Apply default values if required
if (nargin < 2)|isempty(frameSkip), frameSkip = 1; end;
if (nargin < 3)|isempty(dataRange), dataRange = [.25 .75]; end;
if (nargin < 4)|isempty(minSigmaLevel), minSigmaLevel = 0.1; end;
if (nargin < 5)|isempty(shadowLevel), shadowLevel = 0.05; end;
if (nargin < 6)|isempty(significanceThreshold), significanceThreshold = 6; end;
if (nargin < 7)|isempty(useGraph), useGraph = 1; end;
if (nargin < 8)|isempty(graphAlpha), graphAlpha = 1; end;
if (nargin < 9)|isempty(morphStructureElement), morphStructureElement = strel('disk',1); end;
if (nargin < 10)|isempty(returnLargestComponent), returnLargestComponent = 1; end;

if ischar(video)
    % Load the video from an avi file.
    avi = aviread(video);
    pixels = double(cat(4,avi(1:frameSkip:end).cdata))/255;
    clear avi
else
    % Compile the pixel data into a single array
    pixels = double(cat(4,video{1:frameSkip:end}))/255;
    clear video
end;
% NOTE:  The result of the cat operation above may be very large.
% If it will not fit into memory, try increasing frameSkip.
% Another alternative is to decrease the video resolution.

% Convert to HSV color space.
nFrames = size(pixels,4);
for f = 1:nFrames
    pixels(:,:,:,f) = rgb2hsv(pixels(:,:,:,f));
end;

% Generate Gaussian background model in hsv space for each pixel.
[backgroundMean, backgroundDeviation] = hsvGaussModel(pixels,dataRange);
backgroundDeviation(:,:,1) = max(backgroundDeviation(:,:,1),pctile(nonzeros(backgroundDeviation(:,:,1)),minSigmaLevel));
backgroundDeviation(:,:,2) = max(backgroundDeviation(:,:,2),pctile(nonzeros(backgroundDeviation(:,:,2)),minSigmaLevel));
backgroundDeviation(:,:,3) = max(backgroundDeviation(:,:,3),pctile(nonzeros(backgroundDeviation(:,:,3)),minSigmaLevel));

% Do frame-by-frame differencing
connections = pixCon(size(backgroundMean(:,:,1)));
foreground = cell(1,nFrames);
for f = 1:nFrames
    % Find scaled deviation of this frame from background
    deviation = sum(hsvCompare(pixels(:,:,:,f),backgroundMean,shadowLevel)./backgroundDeviation,3);

    % Compare with threshold to generate labeling
    if useGraph
        % Use graph cuts
        alpha = graphAlpha/significanceThreshold;
        label = reshape(graphLabel(alpha*[max(2*significanceThreshold-deviation(:)',0);deviation(:)'],connections),size(deviation));
    else
        % Use morphological operations
        label = double(imopen(imclose(deviation > significanceThreshold,morphStructureElement),morphStructureElement));
    end;

    % clean up:  find largest connected component, etc.
    if returnLargestComponent
        label = fgRanked(label,1);
    end;
    foreground{f} = logical(label);  % Could use bwpack if desired
end;
% end of extractForeground

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fits a gaussian to the specified distribution, in a robust manner.
%
% Returns the mean and standard deviation of the input arguments
% using the specified portion of the data.  Uses HSV color space
% for better shadow control.

function [mu,sigma] = hsvGaussModel(pixels,dataRange)

% statistics of s & v dimensions are simple.
[mu(:,:,2:3),sigma(:,:,2:3)] = fitGauss(pixels(:,:,2:3,:),dataRange,4);
hmu = polarMean(pixels(:,:,1,:)*2*pi,4);
h = zeros(size(pixels(:,:,1,:)));
nFrames = size(pixels,4);
for i = 1:nFrames
    h(:,:,1,i) = mod(pixels(:,:,1,i)*2*pi-hmu+pi,2*pi);
end;
[mu(:,:,1),sigma(:,:,1)] = fitGauss(h,dataRange,4);
mu(:,:,1) = mod(mu(:,:,1)-pi+hmu,2*pi)/(2*pi);
sigma(:,:,1) = sigma(:,:,1)/(2*pi);
% end of hsvGaussModel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fits a gaussian to the specified distribution, in a robust manner.
%
% Returns the mean and standard deviation of the input arguments
% using the specified portion of the data

function [mu,sigma] = fitGauss(m,dataRange,dim)

x = -3:0.01:3;
y = exp(-x.^2);
y = y./sum(y);
yy = cumsum(y);
fdev = sqrt(sum(y.*x.^2));
xbd = min(find(yy > dataRange(1))):min(find(yy > dataRange(2)));
pdev = sqrt(sum(y(xbd).*x(xbd).^2)/sum(y(xbd)));

dimlen = size(m,dim);
if all(dataRange == [0 1])|(dim > length(size(m)))
    % skip sorting step if we'll be using everything anyway
    ms = m;
else
    ms = sort(m,dim);
end;
dimind = round(dataRange*(dimlen-1)+1);
[refstr{1:length(size(m))}] = deal(':');
if (dim <= length(size(m)))
    refstr{dim} = dimind;
end;
sref = struct('type','()','subs',{refstr});
bds = subsref(ms,sref);
if (dim <= length(size(m)))
    sref.subs{dim} = dimind(1):dimind(2);
end;
data = subsref(ms,sref);
mu = mean(data,dim);
sigma = std(data,0,dim)*fdev/pdev;
% end of fitGauss

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% function d = hsvCompare(a,b,shadow)
%
% Returns the component differences between images in HSV color space.

function d = hsvCompare(a,b,shadow)

aclr = length(a(:))/3;
bclr = length(b(:))/3;
sa = a(aclr+1:2*aclr);
sb = b(bclr+1:2*bclr);
va = a(2*aclr+1:3*aclr);
vb = b(2*bclr+1:3*bclr);
d = zeros(size(a));
d(1:aclr) = (angdiff(2*pi*a(1:aclr),2*pi*b(1:bclr)).*min(sa,sb))';
d(aclr+1:2*aclr) = abs(sa-sb);
d(2*aclr+1:3*aclr) = va-vb;
neg = find(d(2*aclr+1:3*aclr)<0);
d(2*aclr+neg) = max(0,vb(neg)-va(neg)-shadow);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%function d = angdiff(th1, th2)
%
% Returns the difference between angles expressed in radians.
% This will always be <= pi.

function d = angdiff(th1, th2)

d = abs(mod(th1-th2+pi,2*pi)-pi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Returns the "mean" of a set of angular variables

function pm = polarMean(th,dim)

th = mod(th,2*pi);
x = cos(th);
y = sin(th);
if (nargin < 2)
    pm = mod(atan2(mean(y),mean(x)),2*pi);
else
    pm = mod(atan2(mean(y,dim),mean(x,dim)),2*pi);
end;
% end of polarMean

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PCTILE Percentile
%
% For vectors, PCTILE(X,Y) returns the element at the Yth percentile of X.
% For matrices, PCTILE(X,Y) returns a row vector containing the
% Yth percentile of each column.  Y must be between 0 and 1.
%
% PCTILE(X,Y,DIM) takes the percentile along dimension DIM of X.
%
% See also MIN, MAX, MEDIAN

function pct = pctile(x,pct,dim)

if (nargin < 3)
    [x,n] = shiftdim(x);
    dim = 1;
end;

sx = sort(x,dim);
subs = cell(1,length(size(x)));
[subs{:}] = deal(':');
subs{dim} = round(pct*(size(x,dim)-1))+1;
ref = struct('type','()','subs',{subs});
pct = subsref(sx,ref);

if (nargin < 3)
    pct = shiftdim(pct,-n);
end;
% end of pctile

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Returns the foreground connected component in the binary image
% supplied that have the specified ranked size(s).

function fgr = fgRanked(bin,rank)

fg = bwlabel(bin,4);
maxfg = max(fg(:));
h = hist(fg(find(bin)),[1:maxfg]);
[sh,sr] = sort(-h);
if (rank < 1)|(rank > max(find(sh > 0)))
    fgr = zeros(size(img))
else
    fgr = (fg==sr(rank));
end;
% end of fgRanked

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function bernoulli_number_test ( )

%*****************************************************************************80
%
%% BERNOULLI_NUMBER_TEST tests BERNOULLI_NUMBER.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BERNOULLI_NUMBER_TEST\n' );
  fprintf ( 1, '  BERNOULLI_NUMBER computes Bernoulli numbers;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I      Exact        Bernoulli\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, c0 ] = bernoulli_number_values ( n_data );

    if ( n_data == 0 )
      break
    end

    c1(1:n+1) = bernoulli_number ( n );

    fprintf ( 1, '  %2d  %14e  %14e\n', n, c0, c1(n+1) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function alpha = anglePoints3d(varargin)
%ANGLEPOINTS3D Compute angle between three 3D points
%
%   ALPHA = anglePoints3d(P1, P2)
%   Computes angle (P1, O, P2), in radians, between 0 and PI.
%
%   ALPHA = anglePoints3d(P1, P2, P3)
%   Computes angle (P1, P2, P3), in radians, between 0 and PI.
%
%   ALPHA = anglePoints3d(PTS)
%   PTS is a 3x3 or 2x3 array containing coordinate of points.
%
%   See also
%   points3d, angles3d
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005.
%

%   HISTORY
%   20/09/2005: add case of single argument for all points
%   04/01/2007: check typo
%   27/05/2014: adjust known vector sizes n1, n0, n2 once corrected for

p2 = [0 0 0];
if length(varargin) == 1
    pts = varargin{1};
    if size(pts, 1)==2
        p1 = pts(1,:);
        p0 = [0 0 0];
        p2 = pts(2,:);
    else
        p1 = pts(1,:);
        p0 = pts(2,:);
        p2 = pts(3,:);
    end

elseif length(varargin) == 2
    p1 = varargin{1};
    p0 = [0 0 0];
    p2 = varargin{2};

elseif length(varargin) == 3
    p1 = varargin{1};
    p0 = varargin{2};
    p2 = varargin{3};
end

% ensure all data have same size
n1 = size(p1, 1);
n2 = size(p2, 1);
n0 = size(p0, 1);

if n1 ~= n0
    if n1 == 1
        p1 = repmat(p1, [n0 1]);
        n1 = n0;
    elseif n0==1
        p0 = repmat(p0, [n1 1]);
    else
        error('Arguments P1 and P0 must have the same size');
    end
end

if n1 ~= n2
    if n1 == 1
        p1 = repmat(p1, [n2 1]);
    elseif n2 == 1
        p2 = repmat(p2, [n1 1]);
    else
        error('Arguments P1 and P2 must have the same size');
    end
end

% normalized vectors
p1 = normalizeVector3d(p1 - p0);
p2 = normalizeVector3d(p2 - p0);

% compute angle
alpha = acos(dot(p1, p2, 2));

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ h, hmin, hmax ] = p19_stepsize ( option )

%*****************************************************************************80
%
%% P19_STEPSIZE returns step sizes for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, real H, HMIN, HMAX, suggested values for the
%    initial step, the minimum step, and the maximum step.
%
  h =    0.125000;
  hmin = 0.015625;
  hmax = 4.000000;

  return
end


text:
function [z_local] = interpolate2local(X_regular,Y_regular,Z_regular,xy_local,int_method)
% [z_local] = interpolate2local(X_regular,Y_regular,Z_regular,xy_local,int_method)
% Interpolates thed ata to a refular grid and fill the gaps outside the convexhull.
% The data is outputed by the new regular grid in X, in Y. The data is a regular grid Z
% with in its thrid dimentsion the the different datasets.
% input:
% X_regular         Regular X-grid in km
% Y_regular         Regular Y_grid in km
% Z_regular         Regular Z-grid (multi-dimensional when having more than 1 dataset)
% int_method  		The interpoaltion method to be used. This is an optional input
%                   argument. By default a triangular interpolation is performed.
% xy_local  		A local xy grid (preferably rotated to reduce interpolation
%                   effects for points outside the convex hull). Needs to be specified
%                   as a 2 column matrix in km.
% output:
% z_local           The data observations that have been interpolated, with each
%                   dataset repressented by a column.
%
%     Copyright (C) 2015  Bekaert David - University of Leeds
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%

% setting the default interpolation method when needed
flag_control_fig = 0;       % when 1 plot the control figures
if nargin < 5
	int_method = [];
end
if isempty(int_method)==1
	int_method = 'linear';
end

% Convert the position grid back to a column vector
x_regular = reshape(X_regular,[],1);
clear X_regular
y_regular = reshape(Y_regular,[],1);
clear Y_regular

% Converting the datagrid back to a column vector or a matrix
n_datasets = size(Z_regular,3);
n_gridpoints = length(y_regular);
z_regular = reshape(Z_regular,n_gridpoints,n_datasets,1);

z_local = NaN([size(xy_local,1) n_datasets]);
for k=1:n_datasets
%     % alternative replace by tri interpolation
%     F = TriScatteredInterp(x_regular,y_regular,z_regular(:,k));
%     z_local_tri(:,k) = F(xy_local);
	z_local(:,k) = griddata(x_regular,y_regular,z_regular(:,k),xy_local(:,1),xy_local(:,2),int_method);
end

if flag_control_fig==1
    figure
    scatter3(xy_local(:,1),xy_local(:,2),z_local(:,1),3,z_local(:,1),'filled');
    hold on
    scatter3(x_regular(:,1),y_regular(:,1),z_regular(:,1),50,z_regular(:,1),'filled');
    view(0,90)
    axis xy
    axis equal
    axis tight
    colorbar
    clear x_regular y_regular z_regular
end

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function line_data = sphere_cubed_lines ( n, line_num )

%*****************************************************************************80
%
%% SPHERE_CUBED_LINES computes the lines on a cubed sphere grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of sections into which each face of
%    the cube is to be divided.
%
%    Input, integer LINE_NUM, the number of lines.
%
%    Output, real LINE_DATA(3,2,LINE_NUM), distinct points on the unit sphere
%    generated by a cubed sphere grid.
%
  line_data = zeros ( 3, 2, line_num );

  l = 0;
%
%  If N = 1, the corners form 12 lines.
%
  if ( n == 1 )
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, n, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, n, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, 0 );

    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, n, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, n, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, n );

    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, n, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, n, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, n );
    return
%
%  If 1 < N, each of 8 corners connects to three neighboring edges.
%
  else
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 1, 0, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, 1, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, 1 );

    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n-1, 0, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, 1, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, 1 );

    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, n, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n-1, n, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, n, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, n-1, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, n, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, n, 1 );

    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 1, n, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, n-1, 0 );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, 0 );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, 1 );

    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 1, 0, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, 1, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, n-1 );

    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n-1, 0, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, 1, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, n-1 );

    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, n, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n-1, n, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, n, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, n-1, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, n, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, n, n-1 );

    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 1, n, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, n-1, n );
    l = l + 1;
    line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, n );
    line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, n-1 );
  end
%
%  If 2 < N, then each of the 12 edges includes lines.
%
  if ( 2 < n )

    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, i,   0, 0 );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, i+1, 0, 0 );
    end
    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n,   i, 0 );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, i+1, 0 );
    end
    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n-i,   n, 0 );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n-i-1, n, 0 );
    end
    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n-i,   0 );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, n-i-1, 0 );
    end

    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, i,   0, n );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, i+1, 0, n );
    end
    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n,   i, n );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, i+1, n );
    end
    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n-i,   n, n );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n-i-1, n, n );
    end
    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n-i,   n );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, n-i-1, n );
    end

    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, i   );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, 0, i+1 );
    end
    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, i   );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, 0, i+1 );
    end
    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, n, i   );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, n, i+1 );
    end
    for i = 1 : n - 2
      l = l + 1;
      line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, i   );
      line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, n, i+1 );
    end

  end
%
%  Lines that belong to one of the six faces.
%
  if ( 1 < n )
%% 000 : nn0
    for i = 1 : n - 1
      for j = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, i, j,   0 );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, i, j+1, 0 );
      end
    end
    for j = 1 : n - 1
      for i = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, i,   j, 0 );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, i+1, j, 0 );
      end
    end
%% 00n : nnn
    for i = 1 : n - 1
      for j = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, i, j,   n );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, i, j+1, n );
      end
    end
    for j = 1 : n - 1
      for i = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, i,   j, n );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, i+1, j, n );
      end
    end
%% 000:n0n
    for i = 1 : n - 1
      for j = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, i, 0, j   );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, i, 0, j+1 );
      end
    end
    for j = 1 : n - 1
      for i = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, i,   0, j );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, i+1, 0, j );
      end
    end
%% 0n0:nnn
    for i = 1 : n - 1
      for j = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, i, n, j   );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, i, n, j+1 );
      end
    end
    for j = 1 : n - 1
      for i = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, i,   n, j );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, i+1, n, j );
      end
    end
%% 000:0nn
    for i = 1 : n - 1
      for j = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, i, j   );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, i, j+1 );
      end
    end
    for j = 1 : n - 1
      for i = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, 0, i,   j );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, 0, i+1, j );
      end
    end
%% n00:nnn
    for i = 1 : n - 1
      for j = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, i, j   );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, i, j+1 );
      end
    end
    for j = 1 : n - 1
      for i = 0 : n - 1
        l = l + 1;
        line_data(1:3,1,l) = sphere_cubed_ijk_to_xyz ( n, n, i,   j );
        line_data(1:3,2,l) = sphere_cubed_ijk_to_xyz ( n, n, i+1, j );
      end
    end

  end

  if ( l ~= line_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPHERE_CUBED_LINES - Fatal error!\n' );
    fprintf ( 1, '  LINE_NUM = %d\n', line_num );
    fprintf ( 1, '  L = %d\n', l );
    error ( 'SPHERE_CUBED_LINES - Fatal error!' );
  end

  return
end


text:
%% Selections of the dataset, model, and category
% Choices of the dataset. You may select 'ilsvrcanimalpart', 'vocpart', or 'cub200'.
dataset='ilsvrcanimalpart'; %'vocpart'; 'cub200';
% If you select dataset='ilsvrcanimalpart', then you may choose of the following categories.
% 'n01443537','n01503061','n01639765','n01662784','n01674464','n01882714','n01982650','n02084071','n02118333','n02121808','n02129165','n02129604','n02131653','n02324045','n02342885','n02355227','n02374451','n02391049','n02395003','n02398521','n02402425','n02411705','n02419796','n02437136','n02444819','n02454379','n02484322','n02503517','n02509815','n02510455'
% If you select dataset='vocpart', then you may choose of the following categories.
% 'bird','cat','cow','dog','horse','sheep'
% If you select dataset='vocpart', then you need to choose categoryName='cub200'.
categoryName='n02118333';
% Choices of the networks. You may select 'vgg-vd-16', 'alexnet', 'vgg-m', or 'vgg-s'.
model='vgg-vd-16'; % 'alexnet'; 'vgg-m'; 'vgg-s';
dropoutRate=0.8; %0.5; 0.6; 0.7; 0.8; 0.9; %when using a small number of training samples, avoid over-fitting.
% Learn a CNN for multi-class classification or single-class classification
isMultiClassClassification=false; %true;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Network configurations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

setup_matconvnet(); % setup matconvnet

switch(dataset) % setup the dataset
    case 'ilsvrcanimalpart'
        setup_ilsvrcanimalpart();
    case 'vocpart'
        setup_vocpart();
    case 'cub200'
        setup_cub200();
    otherwise
        errors('Cannot find the target dataset.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Learning an interpretable CNN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(isMultiClassClassification)
    if(strcmp(dataset,'vocpart'))
        categoryName={'bird','cat','cow','dog','horse','sheep'};
        lossType='ourloss_logistic'; % 'ourloss_softmaxlog';
        learn_icnn_multiclass(model,categoryName,lossType,dropoutRate);
    end
else
    learn_icnn(model,categoryName,dropoutRate);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Showing results and evaluation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(isMultiClassClassification)
    if(strcmp(dataset,'vocpart'))
        categoryName={'bird','cat','cow','dog','horse','sheep'};
    end
end
showResults(categoryName,isMultiClassClassification,model,dataset);

repo: interpretableCNN


text:
% This code is in association with the following paper
% "Ma J, Zhou Z, Wang B, et al. Infrared and visible image fusion based on visual saliency map and weighted least square optimization[J].
% Infrared Physics & Technology, 2017, 82:8-17."
% Authors: Jinlei Ma, Zhiqiang Zhou, Bo Wang, Hua Zong
% Code edited by Jinlei Ma, email: majinlei121@163.com

clear all
close all

label='';
for k = 1:5
if k==1
    label='_gau_001';
end
if k==2
    label='_gau_0005';
end
if k==3
    label='_sp_01';
end
if k==4
    label='_sp_02';
end
if k==5
    label='_poi';
end
disp(label);

for i=1:10
index = i;

% path1 = ['./MF_images/image',num2str(index),'_left.png'];
% path2 = ['./MF_images/image',num2str(index),'_right.png'];
% fused_path = ['./fused_mf/fused',num2str(index),'_wls.png'];

path1 = ['./mf_noise_images/image',num2str(i),label,'_left.png'];
path2 = ['./mf_noise_images/image',num2str(i),label,'_right.png'];
fused_path = ['./fused_mf_noise/fused',num2str(i),label,'_wls.png'];

% I1 is a visible image, and I2 is an infrared image.
I1 = imread(path1);
I2 = imread(path2);

I1 = im2double(I1);
I2 = im2double(I2);

% figure;imshow(I1);
% figure;imshow(I2);
tic
fused = WLS_Fusion(I1,I2);
toc

% figure;imshow(fused);
imwrite(fused,fused_path,'png');

end
end

author: thfylsty


text:
function [C] = collide(map, points)
% COLLIDE Test whether points collide with an obstacle in an environment.
%   C = collide(map, points).  points is an M-by-3 matrix where each
%   row is an (x, y, z) point.  C in an M-by-1 logical vector;
%   C(i) = 1 if M(i, :) touches an obstacle and is 0 otherwise.
map3d = map.map3d_collision;
sz = size(map3d); nx = sz(1); ny = sz(2); nz = sz(3);
% vectorized
points = points_to_idx(map, points);
idx = (points(:,3)-1)*nx*ny + (points(:,2)-1)*nx + points(:,1);
C = map3d(idx) ~= 255 | map3d(nx*ny*nz+idx) ~= 255 | map3d(nx*ny*nz*2+idx) ~= 255;
end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
function net = gp(nin, covar_fn, prior)
%GP	Create a Gaussian Process.
%
%	Description
%
%	NET = GP(NIN, COVARFN) takes the number of inputs NIN  for a Gaussian
%	Process model with a single output, together with a string COVARFN
%	which specifies the type of the covariance function, and returns a
%	data structure NET. The parameters are set to zero.
%
%	The fields in NET are
%	  type = 'gp'
%	  nin = number of inputs
%	  nout = number of outputs: always 1
%	  nwts = total number of weights and covariance function parameters
%	  bias = logarithm of constant offset in covariance function
%	  noise = logarithm of output noise variance
%	  inweights = logarithm of inverse length scale for each input
%	  covarfn = string describing the covariance function:
%	      'sqexp'
%	      'ratquad'
%	  fpar = covariance function specific parameters (1 for squared exponential,
%	   2 for rational quadratic)
%	  trin = training input data (initially empty)
%	  trtargets = training target data (initially empty)
%
%	NET = GP(NIN, COVARFN, PRIOR) sets a Gaussian prior on the parameters
%	of the model. PRIOR must contain the fields PR_MEAN and PR_VARIANCE.
%	If PR_MEAN is a scalar, then the Gaussian is assumed to be isotropic
%	and the additional fields NET.PR_MEAN and PR_VARIANCE are set.
%	Otherwise,  the Gaussian prior has a mean defined by a column vector
%	of parameters PRIOR.PR_MEAN and covariance defined by a column vector
%	of parameters PRIOR.PR_VARIANCE. Each element of PRMEAN corresponds
%	to a separate group of parameters, which need not be mutually
%	exclusive. The membership of the groups is defined by the matrix
%	PRIOR.INDEX in which the columns correspond to the elements of
%	PRMEAN. Each column has one element for each weight in the matrix, in
%	the order defined by the function GPPAK, and each element is 1 or 0
%	according to whether the parameter is a member of the corresponding
%	group or not.  The additional field NET.INDEX is set in this case.
%
%	See also
%	GPPAK, GPUNPAK, GPFWD, GPERR, GPCOVAR, GPGRAD
%

%	Copyright (c) Ian T Nabney (1996-2001)

net.type = 'gp';
net.nin = nin;
net.nout = 1;  % Only do single output GP

% Store log parameters
net.bias = 0;
net.min_noise = sqrt(eps);  % Prevent output noise collapsing completely
net.noise = 0;
net.inweights = zeros(1,nin);  % Weights on inputs in covariance function

covarfns = {'sqexp', 'ratquad'};

if sum(strcmp(covar_fn, covarfns)) == 0
  error('Undefined activation function. Exiting.');
else
  net.covar_fn = covar_fn;
end

switch covar_fn

  case 'sqexp'		% Squared exponential
    net.fpar = zeros(1,1);  % One function specific parameter

  case 'ratquad' 	% Rational quadratic
    net.fpar = zeros(1, 2); % Two function specific parameters

  otherwise
    error(['Unknown covariance function ', covar_fn]);
end

net.nwts = 2 + nin + length(net.fpar);

if nargin >= 3
  if size(prior.pr_mean) == [1 1]
    net.pr_mean = prior.pr_mean;
    net.pr_var = prior.pr_var;
  else
    net.pr_mean = prior.pr_mean;
    net.pr_var = prior.pr_var;
    net.index = prior.index;
  end
end

% Store training data as needed for gpfwd
net.tr_in = [];
net.tr_targets = [];

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function c = tapas_logrt_linear_whatworld_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the linear log-reaction time response model according to as
% developed with Louise Marshall and Sven Bestmann
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The Gaussian noise observation model assumes that responses have a Gaussian distribution around
% the inferred mean of the relevant state. The only parameter of the model is the noise variance
% (NOT standard deviation) zeta.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2014 Christoph Mathys, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'Linear log-reaction time for WhatWorld models';

% Sufficient statistics of Gaussian parameter priors
%
% Beta_0
c.be0mu = log(500);
c.be0sa = 4;

% Beta_1
c.be1mu = 0;
c.be1sa = 4;

% Beta_2
c.be2mu = 0;
c.be2sa = 4;

% Beta_3
c.be3mu = 0;
c.be3sa = 4;

% Zeta
c.logzemu = log(log(20));
c.logzesa = log(2);

% Gather prior settings in vectors
c.priormus = [
    c.be0mu,...
    c.be1mu,...
    c.be2mu,...
    c.be3mu,...
    c.logzemu,...
         ];

c.priorsas = [
    c.be0sa,...
    c.be1sa,...
    c.be2sa,...
    c.be3sa,...
    c.logzesa,...
         ];

% Model filehandle
c.obs_fun = @tapas_logrt_linear_whatworld;

% Handle to function that transforms observation parameters to their native space
% from the space they are estimated in
c.transp_obs_fun = @tapas_logrt_linear_whatworld_transp;

return;

author: translationalneuromodeling

repo: tapas


text:
maxepoch=5;
fprintf(1,'\nFine-tuning deep autoencoder by minimizing cross entropy error. \n');

load batchdata
load mnistvh
load mnisthp
load mnisthp2
load mnistpo

[numcases numdims numbatches]=size(batchdata);
w1=[vishid; hidrecbiases];
w2=[hidpen; penrecbiases];
w3=[hidpen2; penrecbiases2];
w4=[hidtop; toprecbiases];
w5=[hidtop'; topgenbiases];
w6=[hidpen2'; hidgenbiases2];
w7=[hidpen'; hidgenbiases];
w8=[vishid'; visbiases];
l1=size(w1,1)-1;
l2=size(w2,1)-1;
l3=size(w3,1)-1;
l4=size(w4,1)-1;
l5=size(w5,1)-1;
l6=size(w6,1)-1;
l7=size(w7,1)-1;
l8=size(w8,1)-1;
l9=l1;

for epoch = 1:maxepoch
    tt=0;
    for batch = 1:numbatches/10
        fprintf(1,'epoch %d batch %d\r',epoch,batch);
        tt=tt+1;
        data=[];
        for kk=1:10
            data=[data
                batchdata(:,:,(tt-1)*10+kk)];
        end
        max_iter=3;
        VV = [w1(:)' w2(:)' w3(:)' w4(:)' w5(:)' w6(:)' w7(:)' w8(:)']';
        Dim = [l1; l2; l3; l4; l5; l6; l7; l8; l9];

        [X, fX] = minimize(VV,'CG_MNIST',max_iter,Dim,data);

        w1 = reshape(X(1:(l1+1)*l2),l1+1,l2);
        xxx = (l1+1)*l2;
        w2 = reshape(X(xxx+1:xxx+(l2+1)*l3),l2+1,l3);
        xxx = xxx+(l2+1)*l3;
        w3 = reshape(X(xxx+1:xxx+(l3+1)*l4),l3+1,l4);
        xxx = xxx+(l3+1)*l4;
        w4 = reshape(X(xxx+1:xxx+(l4+1)*l5),l4+1,l5);
        xxx = xxx+(l4+1)*l5;
        w5 = reshape(X(xxx+1:xxx+(l5+1)*l6),l5+1,l6);
        xxx = xxx+(l5+1)*l6;
        w6 = reshape(X(xxx+1:xxx+(l6+1)*l7),l6+1,l7);
        xxx = xxx+(l6+1)*l7;
        w7 = reshape(X(xxx+1:xxx+(l7+1)*l8),l7+1,l8);
        xxx = xxx+(l7+1)*l8;
        w8 = reshape(X(xxx+1:xxx+(l8+1)*l9),l8+1,l9);
    end
    savepath=['.\data\mnist_weights' num2str(j) num2str(i)];
    save(savepath, 'w1', 'w2', 'w3', 'w4', 'w5', 'w6', 'w7', 'w8', 'vishid', 'hidrecbiases', 'visbiases', 'hidpen', 'penrecbiases', 'hidgenbiases', 'hidpen2', 'penrecbiases2', 'hidgenbiases2', 'hidtop', 'toprecbiases', 'topgenbiases');
end
%

repo: DeepFi


text:
classdef Perimeter3DCompliance < handle

    properties (Access = public)

    end

    properties (Access = private)
        fileName
        topOptSet
        topOptProblem
    end

    methods (Access = public)

        function obj = Perimeter3DCompliance()
            obj.init();
            obj.createTopOptSettings();
            obj.solveProblem();
        end

    end

    methods (Access = private)

        function init(obj)
          %  obj.fileName = 'CantileverTetra';
            obj.fileName = 'CantileverFlat';
        end

        function createTopOptSettings(obj)
          t = SettingsTopOptProblem(obj.fileName);
          obj.topOptSet = t;
        end

        function solveProblem(obj)
            obj.topOptProblem = TopOpt_Problem(obj.topOptSet);
            obj.topOptProblem.computeVariables();
            obj.topOptProblem.postProcess();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%  Author: Relja Arandjelovic (relja@relja.info)

% This file contains a few examples of how to train and test CNNs for place recognition, refer to README.md for setup instructions, and to our project page for all relevant information (e.g. our paper): http://www.di.ens.fr/willow/research/netvlad/

%  The code samples use the GPU by default, if you want to use the CPU instead (very slow especially for training!), add `'useGPU', false` to the affected function calls (trainWeakly, addPCA, serialAllFeats, computeRepresentation)

% For a tiny example of running the training on a small dataset, which takes only a few minutes to run, refer to the end of this file.

% Set the MATLAB paths
setup;

error('Don''t run this script, it is only meant as a collection of useful commands');

% ---------- Use/test our networks

% Load our network
netID= 'vd16_tokyoTM_conv5_3_vlad_preL2_intra_white';
paths= localPaths();
load( sprintf('%s%s.mat', paths.ourCNNs, netID), 'net' );
net= relja_simplenn_tidy(net);

%  Compute the image representation by simply running the forward pass using
% the network `net` on the appropriately normalized image
% (see `computeRepresentation.m`).

im= vl_imreadjpeg({which('football.jpg')}); im= im{1}; % slightly convoluted because we need the full image path for `vl_imreadjpeg`, while `imread` is not appropriate - see `help computeRepresentation`
feats= computeRepresentation(net, im); % add `'useGPU', false` if you want to use the CPU

% To compute representations for many images, use the `serialAllFeats` function
% which is much faster as it uses batches and it moves the network to
% the GPU only once:
%
%          serialAllFeats(net, imPath, imageFns, outputFn);
%
%  `imageFns` is a cell array containing image file names relative to the `imPath` (i.e. `[imPath, imageFns{i}]` is a valid JPEG image), the representations are saved in binary format (single 4-byte floats). Batch size used for computing the forward pass can be changed by adding the `batchSize` parameter, e.g. `'batchSize', 10`. Note that if your input images are not all of same size (they are in place recognition datasets), you should set `batchSize` to 1.

%  To test the network on a place recognition dataset, set up the test dataset
dbTest= dbTokyo247();

% Set the output filenames for the database/query image representations
paths= localPaths();
dbFeatFn= sprintf('%s%s_%s_db.bin', paths.outPrefix, netID, dbTest.name);
qFeatFn = sprintf('%s%s_%s_q.bin', paths.outPrefix, netID, dbTest.name);

% Compute db/query image representations
serialAllFeats(net, dbTest.dbPath, dbTest.dbImageFns, dbFeatFn, 'batchSize', 10); % adjust batchSize depending on your GPU / network size
serialAllFeats(net, dbTest.qPath, dbTest.qImageFns, qFeatFn, 'batchSize', 1); % Tokyo 24/7 query images have different resolutions so batchSize is constrained to 1

% Measure recall@N
[recall, ~, ~, opts]= testFromFn(dbTest, dbFeatFn, qFeatFn);
plot(opts.recallNs, recall, 'ro-'); grid on; xlabel('N'); ylabel('Recall@N');

% ---------- Full train and test example: Tokyo
% Train: Tokyo Time Machine, Test: Tokyo 24/7

% Set up the train/val datasets
dbTrain= dbTokyoTimeMachine('train');
dbVal= dbTokyoTimeMachine('val');
lr= 0.0001;

% --- Train the VGG-16 network + NetVLAD, tuning down to conv5_1
sessionID= trainWeakly(dbTrain, dbVal, ...
    'netID', 'vd16', 'layerName', 'conv5_3', 'backPropToLayer', 'conv5_1', ...
    'method', 'vlad_preL2_intra', ...
    'learningRate', lr, ...
    'doDraw', true);

% Get the best network
% This can be done even if training is not finished, it will find the best network so far
[~, bestNet]= pickBestNet(sessionID);

% Either use the above network as the image representation extractor (do: finalNet= bestNet), or do whitening (recommended):
finalNet= addPCA(bestNet, dbTrain, 'doWhite', true, 'pcaDim', 4096);

% --- Test

% Set up the test dataset
dbTest= dbTokyo247();

% Set the output filenames for the database/query image representations
paths= localPaths();
dbFeatFn= sprintf('%s%s_ep%06d_%s_db.bin', paths.outPrefix, finalNet.meta.sessionID, finalNet.meta.epoch, dbTest.name);
qFeatFn = sprintf('%s%s_ep%06d_%s_q.bin', paths.outPrefix, finalNet.meta.sessionID, finalNet.meta.epoch, dbTest.name);

% Compute db/query image representations
serialAllFeats(finalNet, dbTest.dbPath, dbTest.dbImageFns, dbFeatFn, 'batchSize', 10); % adjust batchSize depending on your GPU / network size
serialAllFeats(finalNet, dbTest.qPath, dbTest.qImageFns, qFeatFn, 'batchSize', 1); % Tokyo 24/7 query images have different resolutions so batchSize is constrained to 1

% Measure recall@N
[recall, ~, ~, opts]= testFromFn(dbTest, dbFeatFn, qFeatFn);
plot(opts.recallNs, recall, 'ro-'); grid on; xlabel('N'); ylabel('Recall@N');

% --- Test smaller dimensionalities:

% All that needs to be done (only valid for NetVLAD+whitening networks!)
% to reduce the dimensionality of the NetVLAD representation below 4096 to D
% is to keep the first D dimensions and L2-normalize.
% This is done automatically in `testFromFn` using the `cropToDim` option:

cropToDims= [64, 128, 256, 512, 1024, 2048, 4096];
recalls= [];
plotN= 5;
figure;

for iCropToDim= 1:length(cropToDims)
    cropToDim= cropToDims(iCropToDim);
    relja_display('D= %d', cropToDim);
    [recall, ~, ~, opts]= testFromFn(dbTest, dbFeatFn, qFeatFn, [], 'cropToDim', cropToDim);

    whichRecall= find(opts.recallNs==plotN);
    recalls= [recalls, recall(whichRecall)];
    hold off;
    semilogx( cropToDims(1:iCropToDim), recalls, 'bo-');
    set(gca, 'XTick', cropToDims(1:iCropToDim));
    xlabel('Number of dimensions'); ylabel(sprintf('Recall@%d', plotN)); grid on;
    drawnow;
end

% ---------- Pittsburgh examples

% Pitts250k or Pitts30k?
doPitts250k= false;

if doPitts250k
    % Pittsburgh 250k
    lr= 0.0001;
else
    % Pittsburgh 30k
    lr= 0.001;
end

dbTrain= dbPitts(doPitts250k, 'train');
dbVal= dbPitts(doPitts250k, 'val');
dbTest= dbPitts(doPitts250k, 'test');

% Now just run the same code as above for Tokyo

% ---------- Miscellaneous

% --- Other potentially useful commands for training

% caffe vlad: backprop down to conv2

sessionID= trainWeakly(dbTrain, dbVal, ...
    'netID', 'caffe', 'layerName', 'conv5', 'backPropToLayer', 'conv2', ...
    'method', 'vlad_preL2_intra', ...
    'learningRate', lr, ...
    'doDraw', true);

% vgg16 max: backprop down to conv5

sessionID= trainWeakly(dbTrain, dbVal, ...
    'netID', 'vd16', 'layerName', 'conv5_3', 'backPropToLayer', 'conv5_1', ...
    'method', 'max', ...
    'learningRate', lr, ...
    'doDraw', true);

% --- Constructing an off-the-shelf network:

dbTrain= dbPitts(false, 'train'); % or dbTokyoTimeMachine('train');
opts.netID= 'vd16'; % or 'caffe'
opts.layerName= 'conv5_3'; % or whatever layer you want, for caffe: 'conv5'
opts.method= 'vlad_preL2_intra'; % or 'max', 'avg'
net= loadNet(opts.netID, opts.layerName); % load the network and crop at desired layer
net= addLayers(net, opts, dbTrain); % add the NetVLAD/Max/Avg layer
net= addPCA(net, dbTrain, 'doWhite', true, 'pcaDim', 4096, 'batchSize', 10); % add whitening (for non-vlad methods pick a smaller pcaDim value)

% ---------- Tiny dummy training example

% This is just to see if you got all the dependencies and configurations set up correctly.
% Get a tiny version of the Tokyo Time Machine dataset from our research page ( www.di.ens.fr/willow/research/netvlad/ ) as well as the dataset specification. Point paths.dsetRootTokyoTM in localPaths.m to its location. Run the code below to train max pooling on top of AlexNet for this tiny dataset:

dbTrain= dbTiny('train'); dbVal= dbTiny('val');

trainWeakly(dbTrain, dbVal, ...
    'netID', 'caffe', 'layerName', 'conv5', ...
    'method', 'max', 'backPropToLayer', 'conv5', ...
    'margin', 0.1, ...
    'batchSize', 4, 'learningRate', 0.01, 'lrDownFreq', 3, 'momentum', 0.9, 'weightDecay', 0.1, 'compFeatsFrequency', 10, ...
    'nNegChoice', 30, 'nNegCap', 10, 'nNegCache', 10, ...
    'nEpoch', 10, ...
    'epochTestFrequency', 1, 'test0', true, ...
    'nTestSample', inf, 'nTestRankSample', 40, ...
    'saveFrequency', 15, 'doDraw', true, ...
    'useGPU', true, 'numThreads', 12, ...
    'info', 'tiny test');

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
function varargout = pdeset(varargin)
%PDESET    Set options for PDE15S
%   PDESET('NAME1', VALUE1, 'NAME2', VALUE2,...) creates options for the
%   CHEBFUN/PDE15S() routine. It acts as a gateway to ODESET() for the usual ODE
%   options for use in advancing through time, in addition to some new options.
%
%   OPTIONS = PDESET(OLDOPTS, 'NAME1', VALUE1,...) alters an existing options
%   structure OLDOPTS.
%
%   PDESET PROPERTIES (In addition to ODESET properties)
%
%       Eps - Tolerance to use in solving the PDE [ positive scalar {1e-6} ].
%
%       N - Turn off spacial adaptivity. [{NaN} | positive integer  ]
%           Use a fixed spacial grid of size N. If N is NaN, then the automatic
%           procedure is used.
%
%       Plot - Plot the solution at the end of every time chunk. [ on | {off} ]
%              Turning this off can improve speed considerably. Note that if the
%              PDE methods are called with no output arguments, the default
%              behaviour is to plot the solution at the end of every time chunk.
%
%       HoldPlot - Hold the plots after each chunk. [ on | {off} ]
%
%       YLim - Fix the limits of the Y axis if plotting. [ 2x1 vector | {NaN} ]
%           If Ylim is NaN then the imits are determined automatically.
%
%       PlotStyle - Change the plotting options. [ string | ''-'' ].
%
%       PDEflag - Specify which entries correspond to time derivatives.
%           [ vector of logicals {true} ].
%
%       AdjustBCs - Adjust boundary conditions to be satisfied by initial
%                   conditions.
%           [ logical {true} ].
%
%       ODESolver - Function handle determining the ODE solver to use for time
%       stepping the method of lines solution. @ode15s is used by default, but
%       this struggles for non-diffusive problems. In this case, @ode23t is
%       recommended instead.
%           [ {@ode15s} | @ode23t ].
%
% See also ODESET.

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

%  PDE specific options
names = ['Eps      '
         'N        '
         'Plot     '
         'HoldPlot '
         'YLim     '
         'PlotStyle'
         'PDEflag  '
         'AdjustBCs'
         'ODESolver'];

m = size(names, 1);
shortNames = cell(m, 1);
for k = 1:m
    shortNames{k} = strtrim(names(k, :));
end

% Initialise:
opts = {};
pdeOpts = {};

if ( nargin == 0 )
    if ( nargout == 0 )
        odeset;
        fprintf('             Eps: [ positive scalar {1e-6} ]\n')
        fprintf('               N: [ {NaN} | positive integer  ]\n')
        fprintf('            Plot: [ on | {off} ]\n')
        fprintf('        HoldPlot: [ on | {off} ]\n')
        fprintf('            YLim: [ 2x1 vector | {NaN} ]\n')
        fprintf('       PlotStyle: [ string | ''-'']\n')
        fprintf('         PDEflag: [ vector of logicals {true} ]\n')
        fprintf('       AdjustBCs: [ logical {true} ]\n')
        fprintf('       ODESolver: [ {@ode15s} | @ode23t ]\n')
    else
        % Get the ODE opts:
        opts = odeset();
        % Add empty PDE opts:
        for j = 1:m
            opts.(shortNames{j}) = [];
        end
        varargout{1} = opts;
    end

    return
end

% Is an ODESET() / PDESET() structure being passed?
if ( isstruct(varargin{1}) )
    opts = varargin{1};
    varargin(1) = [];
end

% Remember the old PDEOPT values:
for k = 1:m
    namek = shortNames{k};
    if ( isfield(opts, namek) )
        pdeOpts = [ pdeOpts {namek, opts.(namek)}];
    end
end

% Parse the remaining input and update pdeopts entries
k = 1;
while ( k < length(varargin) )
    if ( ~any(strcmpi(fieldnames(odeset), varargin{k})) )
        if ( strcmpi(varargin{k}, 'Plot') || strcmpi(varargin{k}, 'HoldPlot') )
            varargin{k + 1} = onoff(varargin{k + 1});
        end
        pdeOpts = [pdeOpts varargin(k : k + 1)];
        varargin(k : k + 1) = [];
    else
        k = k + 2;
    end
end

% Get the ODE opts:
opts = odeset(opts, varargin{:});

% Add empty PDE opts:
for j = 1:m
    opts.(shortNames{j}) = [];
end

% Attach the PDE opts:
for k = 1:2:length(pdeOpts)
    for j = 1:m
        if ( strcmpi(pdeOpts{k}, shortNames{j}) )
            opts.(shortNames{j}) = pdeOpts{k + 1};
            break
        end
        if ( j == m )
            error('CHEBFUN:pdeset:unknownOption', ...
                ['Unrecognized property name ', pdeOpts{k}, '.'])
        end
    end
end

% Output the opts:
varargout{1} = opts;

end

function out = onoff(in)
% Convert logical values to 'on' or 'off'.
    if ( ~ischar(in) )
        if ( logical(in) )
            out = 'on';
        else
            out = 'off';
        end
    else
        out = in;
    end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [I,V,N]=integrateStochCubAdaptive(func,mu,SR,algorithm,epsVal2,NMax)
%%INTEGRATESTOCHCUBADAPTIVE Perform approximate Monte Carlo numerical
%                integration of a function times a Gaussian PDF using a
%                method that can be compared to importance sampling in that
%                random sets of approximate cubature points are generated
%                rather than more traditional Monte Carlo appraoches.
%
%INPUTS: func The handle to the function that is multiplies by a Gaussian
%             PDF before integration is performed. The function can take a
%             multidimensional input, but must produce a real univariate
%             output.
%          mu The numDimX1 mean of the Gaussian PDF.
%          SR The numDimXnumDim lower-triangular square root of the
%             covariance matrix of the Gaussian PDF.
%   algorithm An optional parameter specifying the algorithm to use.
%             Possible values are:
%             0 Use the first-order method of [1].
%             1 Use the third-order method of [1].
%             2 (The default if omitted or an empty matrix is passed) Use
%               the fifth-order method of [1].
%     epsVal2 An optional convergence bound based on the approximate
%             variance of the estimate. The default if this parameter is
%             omitted or an empty matrix is passed is 1e-2.
%        NMax The maximum number of iterations to perform. The default
%             value if this parameter is omitted or an empty matrix is
%             passed is 10e3.
%
%OUTPUTS: I The approximate scalar value of the integral.
%         V The approximate variance of the estimate.
%         N The number of iterations performed.
%
%All three algorithms arise from similar derivations in [1].
%
%EXAMPLE:
%Here, we choose a function whose true integral with a normal weighting
%function can be easily found. In this instance, we just choose a bivariate
%polynomial.
% f=@(x)(x(2)^4*(x(1)+1)^3-(x(1)+2)^2-(x(2)-4)^3);
% mu=[0;1/2];
% R=[2,-1;
%   -1, 2];
% SR=chol(R,'lower');
% algorithm=2;
% %The exact solution to the problem is known to be
% exactSol=1917/16
% [I,V,N]=integrateStochCubAdaptive(f,mu,SR,algorithm)
%One will find the result is relatively close, but not equal to the exact
%solution.
%
%REFERENCES:
%[1] A. Genz and J. Monahan, "Stochastic integration rules for infinite
%    regions," SIAM Journal on Scientific Computing, vol. 19, no. 2, pp.
%    426-439, Mar. 1998.
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(NMax))
   NMax=10e3;
end

if(nargin<5||isempty(epsVal2))
   epsVal2=1e-2;
end

if(nargin<4||isempty(algorithm))
    algorithm=2;
end

numDim=length(mu);

%Deal with the distribution not being a standard normal distribution.
f=@(x)func(SR*x+mu);

switch(algorithm)
    case 0 %Degree 1 spherical-radial integration rule
        %Step 2
        N=0;
        I=0;
        V=0;
        while(1)
           %Step 3a
           N=N+1;

           %Step3b
           x=randn(numDim,1);

           %Step 3c
           SR=(f(-x)+f(x))/2;
           D=(SR-I)/N;
           I=I+D;
           V=(N-2)*V/N+D^2;

           if(V<epsVal2||N==NMax)
               break;
           end
        end

        %V is approximate sample variance....
    case 1 %Degree 3 spherical-radial integration rule
        %Step 2
        N=0;
        I=0;
        V=0;
        F0=f(zeros(numDim,1));
        e=zeros(numDim,1);
        while(1)
            %Step 3a
            N=N+1;
            SR=0;
            %Step 3b
            Q=randOrthoMat(numDim);
            %Step 3c

            rho=ChiD.rand(1,numDim+2);

            %Step 3d
            for j=1:numDim
                e(j)=1;
                x=rho*Q*e;

                SR=SR+f(-x)+f(x);
                e(j)=0;
            end

            %Step 3e
            SR=F0*(1-numDim/rho^2)+SR/(2*rho^2);
            D=(SR-I)/N;
            I=I+D;
            V=(N-2)*V/N+D^2;

            if(V<epsVal2||N==NMax)
                break;
            end
        end
    case 2%Degree 5 spherical-radial integration rule
        m=numDim;
        %Step 2
        N=0;
        I=0;
        V=0;
        F0=f(zeros(m,1));

        v=regularNSimplexCoords(m);
        while(1)
            %Step 3a
            N=N+1;
            %Step 3b
            Q=randOrthoMat(m);
            vTilde=Q*v;
            %Step 3c
            r=ChiD.rand(1,2*m+7);
            q=BetaD.rand(1,m+2,3/2);
            rho=r*sin(asin(q)/2);
            delta=r*cos(asin(q)/2);
            Fv=0;
            Fy=0;

            %Step 3d
            for j=1:(m+1)
                x1=rho*vTilde(:,j);
                x2=delta*vTilde(:,j);

                Fv=Fv+(m+2-delta^2)*(f(-x1)+f(x1))/(rho^2*(rho^2-delta^2))...
                     +(m+2-rho^2)*(f(-x2)+f(x2))/(delta^2*(delta^2-rho^2));

                for i=1:(j-1)
                    y=(vTilde(:,j)+vTilde(:,i))/norm(vTilde(:,j)+vTilde(:,i));

                    x1=rho*y;
                    x2=delta*y;

                    Fy=Fy+(m+2-delta^2)*(f(-x1)+f(x1))/(rho^2*(rho^2-delta^2))...
                         +(m+2-rho^2)*(f(-x2)+f(x2))/(delta^2*(delta^2-rho^2));
                end
            end

            %Step 3e
            SR=F0*(1-m*(rho^2+delta^2-(m+2))/(rho^2*delta^2))+(Fv*(7-m)*m^2+4*Fy*(m-1)^2)/(2*(m+1)^2*(m+2));
            D=(SR-I)/N;
            I=I+D;
            V=(N-2)*V/N+D^2;

            if(V<epsVal2||N==NMax)
                break;
            end
        end
    otherwise
        error('Unknown algorithm specified.')
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function varargout = mesh(varargin)
%MESH   3-D mesh surface of a CHEBFUN2.
%   MESH is not supported for CHEBFUN2 objects, and so returns an error.
%
% See also SURF.

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

[varargout{1:nargout}] = mesh@separableApprox(varargin{:});

end


text:
function refPh = referencePhase(sessions,roiName,viewType,dataType,fileNameSuffix)

% function refPh = referencePhase(sessions,ROIname,[viewType],[dataType],[fileNameSuffix])
%
% Loops through sessions, averaging across scans to compute
% the average response phase. Uses the output of analyzeSessions so
% you have to run that first on the same sessions and ROI.
%
% datadr: top level data directory, e.g., datadr = /usr/local/mri/mrLoadRet2/'
%
% sessions: struct array with subdir and referenceScans fields, e.g.,
%     sessions(1).subdir = '121897';
%     sessions(1).referenceScans = [1:10];
% where the referenceScans field lists all the scans in that session that
% will be included in the calculation of the referencePhase. This
% might include only one scan per session or it might include all of the scans
% to get the grand mean.
%
% roiName: name of the ROI
%
% fileNameSuffix: specifies data file output from analyzeSessions.
% Default is:
%     fileNameSuffix = '';
% so that the default filename is: roiName.mat
%
% djh, 2/2/2001

mrGlobals

if ~exist('viewType','var')
    viewType = 'Inplane';
end
if ~exist('dataType','var')
    dataType = 'Original';
end
if ~exist('fileNameSuffix','var')
    fileNameSuffix = '';
end

% Loop through sessions
repeats = [];
for s = 1:length(sessions)
	HOMEDIR = sessions(s).path;
    loadSession(HOMEDIR); %changed from loadsession to loadSession (BZL)

    % Open hidden view of the chosen viewType (either 'Inplane' or 'Gray')
    switch viewType
    case 'Inplane'
        view = initHiddenInplane;
    case 'Gray'
        view = initHiddenGray;
    otherwise
        error('analyzeSessions: supports only Inplane and Gray views.');
    end

    % Select chosen data type
    view = selectDataType(view,existDataType(dataType));

    fileName = [roiName,fileNameSuffix];
    load(fullfile(dataDir(view),'Analysis',fileName));
    % Loop through scans
    for scan=sessions(s).referenceScans
        amp = analysis.amps(scan);
        ph = analysis.phs(scan);
        repeats = [repeats, amp*exp(j*ph)];
    end
end
meanRef = mean(repeats(isfinite(repeats)));
refPh = angle(meanRef);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef My_util_dagnn

properties

end

methods(Static)

% -------------------------------------------------------------------------
    function layers = simpleFindLayersOfType(net, type)
    % -------------------------------------------------------------------------
    layers = find(cellfun(@(x)strcmp(x.type, type), net.layers)) ;
    end

    % -------------------------------------------------------------------------
    function net = simpleRemoveLayersOfType(net, type)
    % -------------------------------------------------------------------------
    layers = simpleFindLayersOfType(net, type) ;
    net.layers(layers) = [] ;
    end

    % -------------------------------------------------------------------------
    function layers = dagFindLayersWithOutput(net, outVarName)
    % -------------------------------------------------------------------------
    layers = {} ;
    for l = 1:numel(net.layers)
      if any(strcmp(net.layers(l).outputs, outVarName))
        layers{1,end+1} = net.layers(l);
      end
    end

    end

    % -------------------------------------------------------------------------
    function layers = dagFindLayersOfType(net, type)
    % -------------------------------------------------------------------------
    layers = [] ;
    for l = 1:numel(net.layers)
      if isa(net.layers(l).block, type)
        layers{1,end+1} = net.layers(l).name ;
      end
    end

    end

    % -------------------------------------------------------------------------
    function dagRemoveLayersOfType(net, type)
    % -------------------------------------------------------------------------
    names = dagFindLayersOfType(net, type) ;
    for i = 1:numel(names)
      layer = net.layers(net.getLayerIndex(names{i})) ;
      net.removeLayer(names{i}) ;
      net.renameVar(layer.outputs{1}, layer.inputs{1}, 'quiet', true) ;
    end

    end

    % -------------------------------------------------------------------------
    function [layers, layer_idxes] = dagFindLayersWithInput(net, varName)
    % -------------------------------------------------------------------------
    layers = {} ;
    layer_idxes=[];
    for l = 1:numel(net.layers)
      if any(strcmp(net.layers(l).inputs, varName))
        layers{1,end+1} = net.layers(l) ;
        layer_idxes(1, end+1)=l;
      end
    end

    end

    function f = bilinear_u(k, numGroups, numClasses)
    %BILINEAR_U  Create bilinear interpolation filters
    %   BILINEAR_U(K, NUMGROUPS, NUMCLASSES) compute a square bilinear filter
    %   of size k for deconv layer of depth numClasses and number of groups
    %   numGroups

    factor = floor((k+1)/2) ;
    if rem(k,2)==1
      center = factor ;
    else
      center = factor + 0.5 ;
    end
    C = 1:k ;
    if numGroups ~= numClasses
      f = zeros(k,k,numGroups,numClasses) ;
    else
      f = zeros(k,k,1,numClasses) ;
    end

    for i =1:numClasses
      if numGroups ~= numClasses
        index = i ;
      else
        index = 1 ;
      end
      f(:,:,index,i) = (ones(1,k) - abs(C-center)./factor)'*(ones(1,k) - abs(C-center)./(factor));
    end
    end

end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function pass = test_iszero( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end

% Test with vector 1-1
f = ballfun(ones(19,20,22));
F = ballfunv(f,f,f);
G = F-F;
pass(1) = iszero(G);
end


text:
function sigma = PropError(f,varlist,vals,errs)
%SIGMA = PROPERROR(F,VARLIST,VALS,ERRS)
%
%Finds the propagated uncertainty in a function f with estimated variables
%"vals" with corresponding uncertainties "errs".
%
%varlist is a row vector of variable names. Enter in the estimated values
%in "vals" and their associated errors in "errs" at positions corresponding
%to the order you typed in the variables in varlist.
%
%Example using period of a simple harmonic pendulum:
%
%For this example, lets say the pendulum length is 10m with an uncertainty
%of 1mm, and no error in g.
%syms L g
%T = 2*pi*sqrt(L/g)
%type the function T = 2*pi*sqrt(L/g)
%
%PropError(T,[L g],[10 9.81],[0.001 0])
%ans =
%
%    [       6.3437]    '+/-'    [3.1719e-004]
%    'Percent Error'    '+/-'    [     0.0050]
%
%(c) Brad Ridder 2007. Feel free to use this under the BSD guidelines. If
%you wish to add to this program, just leave my name and add yours to it.
n = numel(varlist);
sig = vpa(ones(1,n));
for i = 1:n
    sig(i) = diff(f,varlist(i),1);
end
error1 =sqrt((sum((subs(sig,varlist,vals).^2).*(errs.^2))));
error = double(error1);
sigma = [{subs(f,varlist,vals)} {'+/-'} {error};
         {'Percent Error'} {'+/-'} {abs(100*(error)/subs(f,varlist,vals))}];

author: Sable

repo: mcbench-benchmarks


text:
%%% AttitudeControl_HIL
clear
path('./icon/',path);

%Constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
%throttle when UAV is hovering
THR_HOVER = 0.609;

%% control parameter
%Attitude PID parameters
Kp_YAW_AngleRate = 0.5;
Ki_YAW_AngleRate = 0.01;
Kd_YAW_AngleRate = 0.00;
%integral saturation
Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
%max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
%max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
%Throttle amplitude
MAX_MAN_THR=0.9;
MIN_MAN_THR=0.05;
%% run simulink model
AttitudeControl_HIL

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function [feven,fodd] = filterGabor1d( r, sig, omega, show )
% Creates an even/odd pair of 1D Gabor filters.
%
% USAGE
%  [feven,fodd] = filterGabor1d( r, sig, omega, [show] )
%
% INPUTS
%  r       - final filter will be 2r+1 (good choice for r is r=2*sig)
%  sig     - standard deviation of Gaussian mask
%  omega   - frequency of underlying sine/cosine in [1/(2r+1) r/(2r+1)]
%  show    - [0] figure to use for optional display
%
% OUTPUTS
%  feven   - even symmetric filter (-cosine masked with Gaussian)
%  fodd    - odd symmetric filter (-sine masked with Gaussian)
%
% EXAMPLE
%  sig = 15; f=filterGabor1d(2*sig,sig,1/sig,1);
%
% See also FILTERGABOR2D
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<4 || isempty(show) ); show=0; end;

r = ceil(r);  n=2*r+1;
if( omega<1/(n) || omega>r/(n) )
  error(['omega=' num2str(omega) ' out of range =[' num2str([1 r]/n) ']']);
end;

% create even and odd pair
x = -r:r;
feven = -cos(2*pi*x*omega) .* exp(-(x.^2)/sig^2);
fodd  = -sin(2*pi*x*omega) .* exp(-(x.^2)/sig^2); %=imag(hilbert(feven));

% normalize to mean==0, but only in locs that are nonzero
inds = abs(feven)>.00001;  feven(inds) = feven(inds) - mean(feven(inds));
inds = abs(fodd)>.00001;  fodd(inds) = fodd(inds) - mean(fodd(inds));

% set L1norm to 0
feven = feven/norm(feven(:),1);
fodd = fodd/norm(fodd(:),1);

% visualization
if( show )
  filterVisualize( feven, show );
  filterVisualize( fodd, show+1 );
end

author: zouchuhang

repo: LayoutNet


text:
function I = surfacearea(f , varargin )
%SURFACEAREA    Surface area of a SEPARABLEAPPROX.
%
%   SURFACEAREA(F) computes the surface area of the SEPARABLEAPPROX in the domain of F.

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

if ( nargin > 1 )
    if ( isa(varargin{1}, 'double') )
        if ( length(varargin{1}) == 4 )
            % Restricted surface area over rectangular region.
            dom = varargin{1};
            f = restrict( f, dom );
        else
            error('CHEBFUN:SEPARABLEAPPROX:surfacearea:domain', 'Bad domain.');
        end
    elseif ( isa(varargin{1}, 'chebfun') )
        f = restrict( f, varargin{1} );
        % Surface area is now just the arc length.
        I = sum( sqrt( 1 + diff( f ).^2 ) );
        return
    else
        error('CHEBFUN:SEPARABLEAPPROX:surfacearea:domain','Bad restricting domain.');
    end
end

% First order derivatives:
fx = diff(f, 1, 2);
fy = diff(f, 1, 1);

% Integrand:
G = 1 + fx.^2 + fy.^2;
S = sqrt( G );

% Surface area:
I = integral2( S );

end


text:
classdef scaleBar < handle
% Inserts a scale bar on the current ebsd or grain map.
%
% Syntax
%   hg = scaleBar(ebsd, scanunits)
%
% Input
%  ebsd      - an mtex ebsd or grain object
%  scanunits - units of the xy coordinates of the ebsd scan (e.g., 'um')
%
% Output
%  oval  - output value
%  ounit - output unit
%
% Options
%  BACKGROUNDCOLOR - background color (ColorSpec)
%  BACKGROUNDALPHA - background transparency (scalar 0<=a<=1)

%
% Example
%  Use a scale bar on the aachen mtexdata.
%
%   mtexdata aachen
%   plot(ebsd)
%   scaleBar(ebsd,'um','BackgroundColor','k','LineColor','w', ...
%                 'Border','off','BackgroundAlpha',0.6,'Location','nw')
%
% Bugs/Issues
%  [1] Using the hardware OpenGL renderer can sometimes cause the text not
%  to appear on the scale bar. This can usually be fixed by adding the line
%
%   opengl software
%
%  _Before_ the call to scaleBar().
%
%  [2] The font and the bounding box do not scale with the figure window.
%  Therefore, the desired figure window size needs to be set before the
%  call to scaleBar().
%
% Authors
%  Eric Payton, eric.payton@bam.de
%  Philippe Pinard, pinard@gfe.rwth-aachen.de
%
% Revision History
%  2012.07.17 EJP - First version submitted for mtex commit.
%  2012.07.27 EJP - Added option for specifying scale bar lengths.

properties (Hidden = true)
  hgt %
end

properties
  txt                      % handle of the text
  shadow                   % handle of the shadow
  ruler                    % handle of the ruler
end

properties (SetObservable)
  backgroundColor = 'k'    % background color (ColorSpec)
  backgroundAlpha = 0.6    % background transparency (scalar 0<=a<=1)
  scanUnit        = 'um'   % units of the xy coordinates of the ebsd scan (e.g., 'um')
  lineColor       = 'w'    % border color and text color (ColorSpec)
  length          = NaN    % desired scale bar length
  % border          = 'off ' % controls whether the box has a border ('on', 'off')
  % borderWidth     = 2      % width of border (scalar)
  % location        = 'sw'   % location of the scale bar ('nw,'ne','sw','se', or [x,y] coordinate vector (len(vector)==2)
end

properties (Dependent = true)
  visible
end

methods

  function sB = scaleBar(mP,scanUnit,varargin)

    sB.scanUnit = scanUnit;
    sB.hgt = hgtransform('parent',mP.ax);
    sB.shadow = patch('parent',sB.hgt,'Faces',1,'Vertices',[NaN NaN NaN]);
    sB.txt = text('parent',sB.hgt,'string','1mm','position',[NaN,NaN],...
      'Interpreter',getMTEXpref('textInterpreter'),'FontSize',getMTEXpref('FontSize'));
    sB.ruler = patch('parent',sB.hgt,'Faces',1,'Vertices',[NaN NaN NaN]);

    % set resize function
    hax = handle(mP.ax);
    try
      hListener(1) = handle.listener(hax, findprop(hax, 'Position'), ...
        'PropertyPostSet',@(a,b) sB.update);
      % save listener, otherwise  callback may die
      setappdata(hax, 'updatePos', hListener);
    catch
      if ~isappdata(hax, 'updatePos')
        hListener = addlistener(hax,'Position','PostSet',...
          @(obj,events) sB.update);
        setappdata(hax, 'updatePos', hListener);
      end
    end

    try
      addlistener(sB,'length','PostSet', @(x,y) sB.update);
      addlistener(sB,'scanUnit','PostSet', @(x,y) sB.update);
      addlistener(sB,'lineColor','PostSet', @(x,y) sB.update);
      addlistener(sB,'backgroundColor','PostSet', @(x,y) sB.update);
      addlistener(sB,'backgroundAlpha','PostSet', @(x,y) sB.update);
      %addlistener(sB,'location','PostSet', @(x,y) sB.update);
      %addlistener(sB,'borderWidth','PostSet', @(x,y) sB.update);
      %addlistener(sB,'border','PostSet', @(x,y) sB.update);
    end

  end

  function setOnTop(sB)
    uistack(sB.hgt,'top')
  end

  function set.visible(sB,value)
    set(sB.hgt,'visible',value);
  end

  function value = get.visible(sB)
    value = get(sB.hgt,'visible');
  end

  function update(sB)

    % get axes orientation
    ax = get(sB.hgt,'parent');
    [el,az] = view(ax);
    xDir = mod((-1)^(az<0) * round(el / 90),4); % E-S-W-N is 0 1 2 3
    yDir = mod(xDir - round(az / 90),4); % E-S-W-N is 0 1 2 3

    % get extent
    dx = xlim(ax); dy = ylim(ax);
    if any(xDir == [1,2]), dx= fliplr(dx); end
    if any(yDir == [1,2]), dy= fliplr(dy); end
    if mod(xDir,2), [dx,dy] = deal(dy,dx); end

    % Find the range in meters for later determination of magnitude
    % We do this so that we never display 10000 nm and always something like
    % 10 microns. Also, the correct choice of units will avoid decimals.
    [sBLength, sBUnit, factor] = switchUnit(0.1*abs(diff(dx)), sB.scanUnit);
    if strcmpi(sBUnit,'um'), sBUnit = '$\mu$m';end

    % we would like to have SBlength beeing a nice number
    if isnan(sB.length)
      goodValues = [1 2 5 10 15 20 25 50 75 100 125 150 200 500 750]; % Possible values for scale bar length
      [~,ind] = min(abs(sBLength-goodValues));
      sB.length = goodValues(ind);
    end
    rulerLength = sB.length * factor * sign(diff(dx));

    % A gap around the bar of 1% of bar length looks nice
    set(sB.txt,'position',[dx(1),dy(1)])
    textHeight = get(sB.txt, 'Extent');
    textHeight = min(textHeight(3:4)) * sign(diff(dy));
    gapY = textHeight/3;
    gapX = abs(gapY) * sign(diff(dx));

    % Box position
    boxWidth = rulerLength + 2.0 * gapX;
    boxx = dx(1) + gapX;
    boxy = dy(1) + gapY;

    % Make bounding box. The z-coordinate is used to put the box under the
    % line.
    verts = [boxx, boxy;
      boxx, boxy +  3*gapY + textHeight;
      boxx + boxWidth, boxy + 3*gapY + textHeight;
      boxx + boxWidth, boxy];
    set(sB.shadow,'Vertices', cP(verts), ...
      'Faces', [1 2 3 4], ...
      'FaceColor', sB.backgroundColor , 'EdgeColor', 'none', ...
      'LineWidth', 1, 'FaceAlpha', sB.backgroundAlpha);

    % update text
    set(sB.txt,'string',['\rm{\textbf{' num2str(sB.length) ' ' sBUnit '}}'],...
      'HorizontalAlignment', 'Center',...
      'VerticalAlignment', 'baseline','color',sB.lineColor,...
      'Position', cP([boxx+boxWidth/2,boxy+3*gapY]));

    % Create line as a patch. The z-coordinate is used to layer the patch over
    % top of the bounding box.
    set(sB.ruler,'Vertices',cP([boxx+gapX, boxy+gapY; ...
      boxx + gapX, boxy+2*gapY; ...
      boxx + gapX + rulerLength, boxy + 2*gapY; ...
      boxx + gapX + rulerLength, boxy + gapY]), ...
      'Faces',[1 2 3 4], 'FaceColor',sB.lineColor, 'FaceAlpha',1);

    sB.setOnTop;

    function pos = cP(pos)
      % interchange x and y if needed
      if mod(xDir,2), pos(:,[1,2]) = pos(:,[2,1]); end
    end

  end

end

end

author: mtex-toolbox

repo: mtex

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


text:
function element_node = gmsh_mesh2d_element_data_example ( element_num, ...
  element_order )

%*****************************************************************************80
%
%% GMSH_MESH2D_ELEMENT_DATA_EXAMPLE returns element information for the example.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Output, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM),
%    the indices of the nodes that make up each element.
%
  element_node = [ ...
    1,  2,  6; ...
    7,  6,  2; ...
    2,  3,  7; ...
    8,  7,  3; ...
    3,  4,  8; ...
    9,  8,  4; ...
    4,  5,  9; ...
   10,  9,  5; ...
    6,  7, 11; ...
   12, 11,  7; ...
    7,  8, 12; ...
   13, 12,  8; ...
    8,  9, 13; ...
   14, 13,  9; ...
    9, 10, 14; ...
   15, 14, 10; ...
   11, 12, 16; ...
   17, 16, 12; ...
   12, 13, 17; ...
   18, 17, 13; ...
   16, 17, 19; ...
   20, 19, 17; ...
   17, 18, 20; ...
   21, 20, 18 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   PRRP example arm with 4 DOF.
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   03/01/2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%Kinematic parameters
robot.DH.theta= '[0 q(2) q(3) 0]';
robot.DH.d='[q(1) 0  0    q(4)]';
robot.DH.a='[0  1  1   0]';
robot.DH.alpha= '[pi/2  0 pi  0]';

%Jacobian matrix. Variation of (X, Y, Z) as a function of (w1, w2, w3)
robot.J='[];';
robot.name='PRRP';

robot.inversekinematic_fn = 'inversekinematic_PRRP(robot, T)';

%number of degrees of freedom
robot.DOF = 4;

%rotational: R, translational: T
robot.kind=['T' 'R' 'R' 'T'];

%minimum and maximum rotation angle in rad
robot.maxangle =[-1 1; %Axis 1, minimum, maximum
                -pi pi;
                -pi pi; %Axis 3, translational
                0 0.5]; %Axis 4

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(200); %Axis 1, rad/s
                deg2rad(200); %Axis 2, m/s
                2; %Axis 3, m/s
                deg2rad(360)]; %Axis 4, rad/s

% end effectors maximum velocity
robot.linear_velmax = 1.0; %m/s

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time

%base reference system
robot.T0 = eye(4);

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

% GRAPHICS
robot.graphical.has_graphics=0;
robot.graphical.color = [25 20 40];
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot
robot.axis=[-3 3 -3 3 -1 1];
%read graphics files
robot = read_graphics(robot);

%DYNAMICS
robot.has_dynamics=0;

repo: ARTE


text:
function ll = fgplvmLogLikelihood(model)

% FGPLVMLOGLIKELIHOOD Log-likelihood for a GP-LVM.
% FORMAT
% DESC returns the log likelihood for a given GP-LVM model.
% ARG model : the model for which the log likelihood is to be
% computed. The model contains the data for which the likelihood is
% being computed in the 'y' component of the structure.
% RETURN ll : the log likelihood of the data given the model.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006, 2009
%
% SEEALSO : gpLogLikelihood, fgplvmCreate
%
% MODIFICATIONS : Carl Henrik Ek, 2008, 2009
%
% FGPLVM

ll = gpLogLikelihood(model);

if isfield(model, 'dynamics') && ~isempty(model.dynamics)
  % A dynamics model is being used.
  ll = ll + modelLogLikelihood(model.dynamics);
elseif isfield(model, 'prior') &&  ~isempty(model.prior)
  for i = 1:model.N
    ll = ll + priorLogProb(model.prior, model.X(i, :));
  end
end

switch model.approx
  case {'dtc', 'dtcvar', 'fitc', 'pitc'}
   if isfield(model, 'inducingPrior') && ~isempty(model.inducingPrior)
     for i = 1:model.k
       ll = ll + priorLogProb(model.inducingPrior, model.X_u(i, :));
     end
   end
 otherwise
  % do nothing
end

if(isfield(model,'constraints')&&~isempty(model.constraints))
  for(i = 1:1:model.constraints.numConstraints)
    ll = ll + constraintLogLikelihood(model.constraints.comp{i},model.X);
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ element_area, mesh_area ] = area_q4_mesh ( node_num, element_num, ...
  node_xy, element_node )

%*****************************************************************************80
%
%% AREA_Q4_MESH computes areas of elements in a Q4 mesh.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, real NODE_XY(2,NODE_NUM), the node coordinates.
%
%    Input, integer ELEMENT_NODE(4,ELEMENT_NUM), lists the
%    nodes that make up each element, in counterclockwise order.
%
%    Output, real ELEMENT_AREA(ELEMENT_NUM), the element areas.
%
%    Output, real MESH_AREA, the mesh area.
%
  for element = 1 : element_num
    for node = 1 : 4
      q4(1:2,node) = node_xy(1:2,element_node(node,element));
    end
    element_area(element) = area_quad ( q4 );
  end

  mesh_area = sum ( element_area(1:element_num) );

  return
end


text:
function L=computeLap(X)
% Compute graph Laplacian

optL=[];
optL.NeighborMode='KNN';
optL.k=5;
optL.WeightMode='Cosine';   % 'Cosine' or 'HeatKernel'
optL.t=1;
L=constructW(X,optL);

author: zzstefan

repo: BrainNetClass


text:
function h = plot_tensor_field(H, M, options)

% plot_tensor_field - display a tensor field
%
%   h = plot_tensor_field(H, M, options);
%
%   options.sub controls sub-sampling
%   options.color controls color
%
%   Copyright (c) 2006 Gabriel Peyre

if nargin<3
    options.null = 0;
end
if not( isstruct(options) )
    sub = options;
    clear options;
    options.sub = sub;
end

% sub = getoptions(options, 'sub', 1);
sub = getoptions(options, 'sub', round(size(H,1)/30) );
color = getoptions(options, 'color', 'r');

if nargin<2
    M = [];
end

if not(isempty(M)) && size(M,3)==1
    M = repmat(M, [1 1 3]); % ensure B&W image
end

if size(H,3)==3 && size(H,4)==1
    H = cat(3, H(:,:,1), H(:,:,3), H(:,:,3), H(:,:,2) );
    H = reshape(H, size(H,1), size(H,2), 2, 2);
    if 0
    % flip the main eigen-axes
    [e1,e2,l1,l2] = perform_tensor_decomp(H);
    H = perform_tensor_recomp(e2,e1,l1,l2);
    end
    h = plot_tensor_field(H, M, sub);
    return;
end

% swap X and Y axis
%%% TODO
a = H(:,:,2,2);
H(:,:,2,2) = H(:,:,1,1);
H(:,:,1,1) = a;

hold on;
if ~isempty(M)
    imagesc(rescale(M)); drawnow;
end
h = fn_tensordisplay(H(:,:,1,1),H(:,:,1,2), H(:,:,2,2), 'sub', sub, 'color', color);
axis image; axis off;
colormap jet(256);
% hold off;

function h = fn_tensordisplay(varargin)

% function h = fn_tensordisplay([X,Y,]Txx,Txy,Tyy[,'sigma',sigma][,'sub',sub][,color][,patch options...]])
% function h = fn_tensordisplay([X,Y,]e[,'sigma',sigma][,'sub',sub][,color][,patch options...]])

% X,Y,Txx,Txy,Tyy
if isstruct(varargin{1}) || isstruct(varargin{3})
    if isstruct(varargin{1}), nextarg=1; else nextarg=3; end
    e = varargin{nextarg};
    Txx = e.ytyt;
    Txy = -e.ytyx;
    Tyy = e.yxyx;
    if nextarg==1
        [nj ni] = size(Txx);
        [X Y] = meshgrid(1:ni,1:nj);
    else
        [X Y] = deal(varargin{1:2});
    end
    nextarg = nextarg+1;
else
    [nj ni] = size(varargin{3});
    if nargin<5 || ischar(varargin{4}) || ischar(varargin{5}) || any(size(varargin{5})~=[nj ni])
        [X Y] = meshgrid(1:ni,1:nj);
        nextarg = 1;
    else
        [X Y] = deal(varargin{1:2});
        nextarg = 3;
    end
    [Txx Txy Tyy] = deal(varargin{nextarg:nextarg+2});
    nextarg = nextarg+3;
end
if any(size(X)==1), [X Y] = meshgrid(X,Y); end
[nj,ni] = size(X);
if any(size(Y)~=[nj ni]) || ...
        any(size(Txx)~=[nj ni]) || any(size(Txy)~=[nj ni]) || any(size(Tyy)~=[nj ni])
    error('Matrices must be same size')
end
% sigma, sub, color
color = 'r';
while nextarg<=nargin
    flag = varargin{nextarg};
    nextarg=nextarg+1;
    if ~ischar(flag), color = flag; continue, end
    switch lower(flag)
        case 'sigma'
            sigma = varargin{nextarg};
            nextarg = nextarg+1;
            switch length(sigma)
                case 1
                    sigmax = sigma;
                    sigmay = sigma;
                case 2
                    sigmax = sigma(1);
                    sigmay = sigma(2);
                otherwise
                    error('sigma definition should entail two values');
            end
            h = fspecial('gaussian',[ceil(2*sigmay) 1],sigmay)*fspecial('gaussian',[1 ceil(2*sigmax)],sigmax);
            Txx = imfilter(Txx,h,'replicate');
            Txy = imfilter(Txy,h,'replicate');
            Tyy = imfilter(Tyy,h,'replicate');
        case 'sub'
            sub = varargin{nextarg};
            nextarg = nextarg+1;
            switch length(sub)
                case 1
                    [x y] = meshgrid(1:sub:ni,1:sub:nj);
                    sub = y+nj*(x-1);
                case 2
                    [x y] = meshgrid(1:sub(1):ni,1:sub(2):nj);
                    sub = y+nj*(x-1);
            end
            X = X(sub); Y = Y(sub);
            Txx = Txx(sub); Txy = Txy(sub); Tyy = Tyy(sub);
            [nj ni] = size(sub);
        case 'color'
            color = varargin{nextarg};
            nextarg = nextarg+1;
        otherwise
            break
    end
end
% options
options = {varargin{nextarg:end}};

npoints = 50;
theta = (0:npoints-1)*(2*pi/npoints);
circle = [cos(theta) ; sin(theta)];
Tensor = cat(3,Txx,Txy,Txy,Tyy);        % jdisplay x idisplay x tensor
Tensor = reshape(Tensor,2*nj*ni,2);     % (display x 1tensor) x 2tensor
Ellipse = Tensor * circle;              % (display x uv) x npoints
Ellipse = reshape(Ellipse,nj*ni,2,npoints);         % display x uv x npoints
XX = repmat(X(:),1,npoints);                        % display x npoints
YY = repmat(Y(:),1,npoints);                        % display x npoints
U  = reshape(Ellipse(:,1,:),nj*ni,npoints);         % display x npoints
V  = reshape(Ellipse(:,2,:),nj*ni,npoints);         % display x npoints
umax = max(U')'; vmax = max(V')';
umax(umax==0)=1; vmax(vmax==0)=1;
if ni==1, dx=1; else dx = X(1,2)-X(1,1); end
if nj==1, dy=1; else dy = Y(2,1)-Y(1,1); end
fact = min(dx./umax,dy./vmax)*.35;
fact = repmat(fact,1,npoints);
U = XX + fact.*U;
V = YY + fact.*V;

%-----------
MM = mmax(Txx+Tyy);
mm = mmin(Txx+Tyy);
[S1, S2] = size(U);
Colormap = zeros(S2, S1, 3);
Map = colormap(jet(256));
for k =1:npoints
    Colormap(k,:,1) = Map(floor(255*(Txx(:)+Tyy(:)-mm)/(MM-mm)) + 1, 1);
    Colormap(k,:,2) = Map(floor(255*(Txx(:)+Tyy(:)-mm)/(MM-mm)) + 1, 2);
    Colormap(k,:,3) = Map(floor(255*(Txx(:)+Tyy(:)-mm)/(MM-mm)) + 1, 3);
end
%-----------
%h = fill(U',V',color,'EdgeColor',color,options{:});
h = fill(U',V',Colormap,'EdgeColor', 'interp');

axis ij;

if nargout==0, clear h, end

function a=mmax(a)
a = max(a(:));
function a=mmin(a)
a = min(a(:));

author: gpeyre

repo: matlab-toolboxes


text:
% GAUSSIAN_DLOGPDF_COV - The gradient of the log pdf of a multivariate
%                        normal distribution with respect to the covariance
%                        matrix.
%
% Computes the gradient of the logarithm of the multivariate normal
% probability density function N(Y|MU,COV), where Y is the random variable,
% MU is the mean and COV is the covariance. The gradient is evaluated with
% respect to the covariance matrix COV and returned in matrix form.
%
%   d(LOG(N(Y|MU,COV))) / dCOV
%
% Usage:
%
%   DX = GAUSSIAN_DLOGPDF_COV(INVCOV_Y, INVCOV_MU, INVCOV)
%
% where
%
% INVCOV_Y   : INV(COV)*Y
% INVCOV_MU  : INV(COV)*MU
% INVCOV     : INV(COV)
%
% Letting the user to compute the terms allows greater efficiency and
% flexibility.

% Last modified 2010-11-19
% Copyright (c) Jaakko Luttinen

function dx = gaussian_dlogpdf_cov(invCov_y, invCov_mu, invCov)

% $$$ if nargin == 3
% $$$   y = varargin{1};
% $$$   mu = varargin{2};
% $$$   L = varargin{3};
% $$$   invCov_y = linsolve_chol(L, y, 'lower');
% $$$   invCov_mu = linsolve_chol(L, mu, 'lower');
% $$$   invCov = linsolve_chol(L, eye(size(L)), 'lower');
% $$$ end

dx = - 0.5 * invCov ...
     + 0.5 * (invCov_y * invCov_y'  - ...
              invCov_y * invCov_mu' - ...
              invCov_mu * invCov_y' + ...
              invCov_mu * invCov_mu');

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%  Figure 7.53      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% script to generate Fig. 7.53
%
clf;
dp=[1 1 0];
np=[1];
nc=conv([1 1.001],[8.32 0.8]);
dc=conv([1 4.08],[1 0.0196]);
num=conv(np,nc);
den=conv(dp,dc);
dcl=[0 0 num]+den; % closed-loop denominator
t=0:.1:5;
y=step(num,dcl,t);
plot(t,y);
xlabel('Time (sec)');
ylabel('y(t)');
nicegrid;
title('Fig.7.53: Step response for lag compensation design')


text:
function [n, ts_left, ts_right] = nex_int(filename, varname)
% nex_int(filename, varname): Read interval variable from a .nex file
%
% [n, ts_left, ts_right] = nex_int(filename, varname)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   varname - variable name
% OUTPUT:
%   n - number of intervals
%   ts_left - array of left ends of the intervals (in seconds)
%   ts_right - array of right ends of the intervals (in seconds)

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

n = 0;
ts_left = 0;
ts_right = 0;

if(nargin ~= 2)
   disp('2 input arguments are required')
   return
end

if(ischar(filename) == 0)
   disp('input arguments should be character arrays')
   return
end

if(ischar(varname) == 0)
   disp('input arguments should be character arrays')
   return
end

if(length(filename) == 0)
   [fname, pathname] = uigetfile('*.nex', 'Select a Nex file');
    filename = strcat(pathname, fname);
end

fid = fopen(filename, 'r', 'ieee-le');
if(fid == -1)
    disp('cannot open file');
   return
end

disp(strcat('file = ', filename));
magic = fread(fid, 1, 'int32');
version = fread(fid, 1, 'int32');
comment = fread(fid, 256, 'char');
freq = fread(fid, 1, 'double');
tbeg = fread(fid, 1, 'int32');
tend = fread(fid, 1, 'int32');
nvar = fread(fid, 1, 'int32');
fseek(fid, 260, 'cof');
name = zeros(1, 64);
found = 0;
for i=1:nvar
    type = fread(fid, 1, 'int32');
    var_version = fread(fid, 1, 'int32');
    name = fread(fid, [1 64], 'char');
    offset = fread(fid, 1, 'int32');
    n = fread(fid, 1, 'int32');
    name = char(name);
    name = deblank(name);
    k = strcmp(name, deblank(varname));
    if(k == 1)
        if type ~= 2
            disp(sprintf('%s is not an interval variable', deblank(varname)));
            return;
        end
        found = 1;
        fseek(fid, offset, 'bof');
        ts_left = fread(fid, [1 n], 'int32');
        ts_right = fread(fid, [1 n], 'int32');
        break
    end
    dummy = fread(fid, 128, 'char');
end

fclose(fid);

if found == 0
    disp('did not find variable in the file');
else
    ts_left = ts_left/freq;
    ts_right = ts_right/freq;
    disp(strcat('number of intervals = ', num2str(n)));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [event] = read_brainstorm_event(filename)

% READ_BRAINSTORM_EVENT reads the event information from .EEG files
% that have been generated by the Nihon Kohden system. The function
% constitutes a wrapper around BrainStorm3 functionalities
%
% Use as
%   [event] = read_brainstorm_event(filename)
%
% See also READ_NK1200_HEADER, READ_NK1200_DATA

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

% call BrainStorm's in_fopen_nk
[sFile] = in_fopen_nk(filename);

% create the event file by pulling values from the output
event = [];
for e = 1:numel(sFile.events)
  event(e).type      = sFile.events(e).label;
  event(e).sample    = sFile.events(e).samples(1);
  event(e).value     = sFile.events(e).label;
  event(e).offset    = 0;
  if numel(sFile.events(e).samples)>1
    event(e).duration  = sFile.events(e).samples(2)-sFile.events(e).samples(1)+1;
  else
    event(e).duration  = 0;
  end
end


text:
%% MD_RUN is a script to run MD with prescribed input.
%
  filename = 'md_output.txt';

  delete ( filename );

  diary ( filename );

  nd = 3;
  np = 500;
  step_num = 500;
  dt = 0.1;

  md ( nd, np, step_num, dt );

  diary ( 'off' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Output written to "%s"\n', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Normal end of execution.\n' );


text:
function rgbim = binary2rgb(binaryim,rgbtriplet)
%BINARY2RGB Converts high values in a binary image to one RGB color.
%   RGBIM = BINARY2RGB(BINARYIM,RGBTRIPLET) converts the high values in
%   bivalued (binary) image BINARYIM to a single color specified in
%   triplet RGBTRIPLET, which is of the form [r,g,b], where r, g, and b
%   are 0 or 1, and define the color of the output. For example, [1,0,0]
%   specifies red and [1 1 0] specifies yellow.
%
%   This is a DIPUM3E utility function.
%
%   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 the input is not already logical, convert to logical, treating
% anything greater than the minimum value as high (foreground).
if ~islogical(binaryim)
   binaryim = binaryim > min(binaryim(:));
end

% Initialize three component images.
red = zeros(size(binaryim));
green = red;
blue = red;

% Insert the corresponding RGB triplet value into the component images
% based on the foreground pixel locations in binaryim.
red(binaryim) = rgbtriplet(1);
green(binaryim) = rgbtriplet(2);
blue(binaryim) = rgbtriplet(3);

% Form RGB output image from the three component images.
rgbim = cat(3,red,green,blue);

author: dipum

repo: dipum-toolbox

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


text:
function const_pts = rectQAM_const(M)
%

% M-PAM bit/symbol ordering
[bo so] = PAM_Gray_Code;

k = log2(M);
I = 2^ceil(k/2);
Q = 2^floor(k/2);

% Produce PAM constellation points for in-phase and quadrature
I_t = -(I-1):2:I-1;
Q_t = Q-1:-2:-(Q-1);

% Re-order based on the PAM Gray mapping
for i=1:length(I_t)
    I_pts(so{log2(I)}(i)+1) = I_t(i);
end
for i=1:length(Q_t)
    Q_pts(so{log2(Q)}(i)+1) = Q_t(i);
end

for i=0:M-1
    MSBs = floor(i/I);
    LSBs = i-MSBs*I;
    QAM_I(i+1) = I_pts(LSBs+1);
    QAM_Q(i+1) = Q_pts(MSBs+1);
end

const_pts=complex(QAM_I,QAM_Q);