text:
%  internal function

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)

function [hdr, filetype, fileprefix, machine] = load_nii_hdr(fileprefix)

   if ~exist('fileprefix','var'),
      error('Usage: [hdr, filetype, fileprefix, machine] = load_nii_hdr(filename)');
   end

   machine = 'ieee-le';
   new_ext = 0;

   if ~isempty(findstr('.nii',fileprefix)) && ...
                        strcmp(fileprefix(end-3:end), '.nii')
      new_ext = 1;
      fileprefix(end-3:end)='';
   end

   if ~isempty(findstr('.hdr',fileprefix)) && ...
                        strcmp(fileprefix(end-3:end), '.hdr')
      fileprefix(end-3:end)='';
   end

   if ~isempty(findstr('.img',fileprefix)) && ...
                        strcmp(fileprefix(end-3:end), '.img')
      fileprefix(end-3:end)='';
   end

   if new_ext
      fn = sprintf('%s.nii',fileprefix);

      if ~exist(fn)
         msg = sprintf('Cannot find file "%s.nii".', fileprefix);
         error(msg);
      end
   else
      fn = sprintf('%s.hdr',fileprefix);

      if ~exist(fn)
         msg = sprintf('Cannot find file "%s.hdr".', fileprefix);
         error(msg);
      end
   end

   fid = fopen(fn,'r',machine);

   if fid < 0,
      msg = sprintf('Cannot open file %s.',fn);
      error(msg);
   else
      fseek(fid,0,'bof');

      if fread(fid,1,'int32') == 348
         hdr = read_header(fid);
         fclose(fid);
      else
         fclose(fid);

         %  first try reading the opposite endian to 'machine'
         %
         switch machine,
         case 'ieee-le', machine = 'ieee-be';
         case 'ieee-be', machine = 'ieee-le';
         end

         fid = fopen(fn,'r',machine);

         if fid < 0,
            msg = sprintf('Cannot open file %s.',fn);
            error(msg);
         else
            fseek(fid,0,'bof');

            if fread(fid,1,'int32') ~= 348

               %  Now throw an error
               %
               msg = sprintf('File "%s" is corrupted.',fn);
               error(msg);
            end

            hdr = read_header(fid);
            fclose(fid);
         end
      end
   end

   if strcmp(hdr.hist.magic, 'n+1')
      filetype = 2;
   elseif strcmp(hdr.hist.magic, 'ni1')
      filetype = 1;
   else
      filetype = 0;
   end

   return					% load_nii_hdr

%---------------------------------------------------------------------
function [ dsr ] = read_header(fid)

        %  Original header structures
	%  struct dsr
	%       {
	%       struct header_key hk;            /*   0 +  40       */
	%       struct image_dimension dime;     /*  40 + 108       */
	%       struct data_history hist;        /* 148 + 200       */
	%       };                               /* total= 348 bytes*/

    dsr.hk   = header_key(fid);
    dsr.dime = image_dimension(fid);
    dsr.hist = data_history(fid);

    %  For Analyze data format
    %
    if ~strcmp(dsr.hist.magic, 'n+1') && ~strcmp(dsr.hist.magic, 'ni1')
        dsr.hist.qform_code = 0;
        dsr.hist.sform_code = 0;
    end

    return					% read_header

%---------------------------------------------------------------------
function [ hk ] = header_key(fid)

    fseek(fid,0,'bof');

	%  Original header structures
	%  struct header_key                     /* header key      */
	%       {                                /* off + size      */
	%       int sizeof_hdr                   /*  0 +  4         */
	%       char data_type[10];              /*  4 + 10         */
	%       char db_name[18];                /* 14 + 18         */
	%       int extents;                     /* 32 +  4         */
	%       short int session_error;         /* 36 +  2         */
	%       char regular;                    /* 38 +  1         */
	%       char dim_info;   % char hkey_un0;        /* 39 +  1 */
	%       };                               /* total=40 bytes  */
	%
	% int sizeof_header   Should be 348.
	% char regular        Must be 'r' to indicate that all images and
	%                     volumes are the same size.

    v6 = version;
    if str2num(v6(1))<6
       directchar = '*char';
    else
       directchar = 'uchar=>char';
    end

    hk.sizeof_hdr    = fread(fid, 1,'int32')';	% should be 348!
    hk.data_type     = deblank(fread(fid,10,directchar)');
    hk.db_name       = deblank(fread(fid,18,directchar)');
    hk.extents       = fread(fid, 1,'int32')';
    hk.session_error = fread(fid, 1,'int16')';
    hk.regular       = fread(fid, 1,directchar)';
    hk.dim_info      = fread(fid, 1,'uchar')';

    return					% header_key

%---------------------------------------------------------------------
function [ dime ] = image_dimension(fid)

	%  Original header structures
	%  struct image_dimension
	%       {                                /* off + size      */
	%       short int dim[8];                /* 0 + 16          */
        %       /*
        %           dim[0]      Number of dimensions in database; usually 4.
        %           dim[1]      Image X dimension;  number of *pixels* in an image row.
        %           dim[2]      Image Y dimension;  number of *pixel rows* in slice.
        %           dim[3]      Volume Z dimension; number of *slices* in a volume.
        %           dim[4]      Time points; number of volumes in database
        %       */
	%       float intent_p1;   % char vox_units[4];   /* 16 + 4       */
	%       float intent_p2;   % char cal_units[8];   /* 20 + 4       */
	%       float intent_p3;   % char cal_units[8];   /* 24 + 4       */
	%       short int intent_code;   % short int unused1;   /* 28 + 2 */
	%       short int datatype;              /* 30 + 2          */
	%       short int bitpix;                /* 32 + 2          */
	%       short int slice_start;   % short int dim_un0;   /* 34 + 2 */
	%       float pixdim[8];                 /* 36 + 32         */
	%	/*
	%		pixdim[] specifies the voxel dimensions:
	%		pixdim[1] - voxel width, mm
	%		pixdim[2] - voxel height, mm
	%		pixdim[3] - slice thickness, mm
	%		pixdim[4] - volume timing, in msec
	%					..etc
	%	*/
	%       float vox_offset;                /* 68 + 4          */
	%       float scl_slope;   % float roi_scale;     /* 72 + 4 */
	%       float scl_inter;   % float funused1;      /* 76 + 4 */
	%       short slice_end;   % float funused2;      /* 80 + 2 */
	%       char slice_code;   % float funused2;      /* 82 + 1 */
	%       char xyzt_units;   % float funused2;      /* 83 + 1 */
	%       float cal_max;                   /* 84 + 4          */
	%       float cal_min;                   /* 88 + 4          */
	%       float slice_duration;   % int compressed; /* 92 + 4 */
	%       float toffset;   % int verified;          /* 96 + 4 */
	%       int glmax;                       /* 100 + 4         */
	%       int glmin;                       /* 104 + 4         */
	%       };                               /* total=108 bytes */

    dime.dim        = fread(fid,8,'int16')';
    dime.intent_p1  = fread(fid,1,'float32')';
    dime.intent_p2  = fread(fid,1,'float32')';
    dime.intent_p3  = fread(fid,1,'float32')';
    dime.intent_code = fread(fid,1,'int16')';
    dime.datatype   = fread(fid,1,'int16')';
    dime.bitpix     = fread(fid,1,'int16')';
    dime.slice_start = fread(fid,1,'int16')';
    dime.pixdim     = fread(fid,8,'float32')';
    dime.vox_offset = fread(fid,1,'float32')';
    dime.scl_slope  = fread(fid,1,'float32')';
    dime.scl_inter  = fread(fid,1,'float32')';
    dime.slice_end  = fread(fid,1,'int16')';
    dime.slice_code = fread(fid,1,'uchar')';
    dime.xyzt_units = fread(fid,1,'uchar')';
    dime.cal_max    = fread(fid,1,'float32')';
    dime.cal_min    = fread(fid,1,'float32')';
    dime.slice_duration = fread(fid,1,'float32')';
    dime.toffset    = fread(fid,1,'float32')';
    dime.glmax      = fread(fid,1,'int32')';
    dime.glmin      = fread(fid,1,'int32')';

    return					% image_dimension

%---------------------------------------------------------------------
function [ hist ] = data_history(fid)

	%  Original header structures
	%  struct data_history
	%       {                                /* off + size      */
	%       char descrip[80];                /* 0 + 80          */
	%       char aux_file[24];               /* 80 + 24         */
	%       short int qform_code;            /* 104 + 2         */
	%       short int sform_code;            /* 106 + 2         */
	%       float quatern_b;                 /* 108 + 4         */
	%       float quatern_c;                 /* 112 + 4         */
	%       float quatern_d;                 /* 116 + 4         */
	%       float qoffset_x;                 /* 120 + 4         */
	%       float qoffset_y;                 /* 124 + 4         */
	%       float qoffset_z;                 /* 128 + 4         */
	%       float srow_x[4];                 /* 132 + 16        */
	%       float srow_y[4];                 /* 148 + 16        */
	%       float srow_z[4];                 /* 164 + 16        */
	%       char intent_name[16];            /* 180 + 16        */
	%       char magic[4];   % int smin;     /* 196 + 4         */
	%       };                               /* total=200 bytes */

    v6 = version;
    if str2num(v6(1))<6
       directchar = '*char';
    else
       directchar = 'uchar=>char';
    end

    hist.descrip     = deblank(fread(fid,80,directchar)');
    hist.aux_file    = deblank(fread(fid,24,directchar)');
    hist.qform_code  = fread(fid,1,'int16')';
    hist.sform_code  = fread(fid,1,'int16')';
    hist.quatern_b   = fread(fid,1,'float32')';
    hist.quatern_c   = fread(fid,1,'float32')';
    hist.quatern_d   = fread(fid,1,'float32')';
    hist.qoffset_x   = fread(fid,1,'float32')';
    hist.qoffset_y   = fread(fid,1,'float32')';
    hist.qoffset_z   = fread(fid,1,'float32')';
    hist.srow_x      = fread(fid,4,'float32')';
    hist.srow_y      = fread(fid,4,'float32')';
    hist.srow_z      = fread(fid,4,'float32')';
    hist.intent_name = deblank(fread(fid,16,directchar)');
    hist.magic       = deblank(fread(fid,4,directchar)');

    fseek(fid,253,'bof');
    hist.originator  = fread(fid, 5,'int16')';

    return					% data_history

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
 function Asums = osem_sum(G, yi, ci, ri, block_order)
%	precompute some sums needed for OSEM
%	NOT DONE!!!!!!!!!!!
%	block_order is a vector of length nblock
%	that specifies the starting angle for each block
%	simplest case is just block_order = 1:nblock,
%	but that choice is suboptimal
%	if block_order is an empty matrix, then 1 block is used
%	if block_order is a scalar (nblock) power of 2 != 1,
%	then the "bit-reversal ordering" is used
%	(nblock can be 1 to na)
%
%	Copyright Apr 1999, Jeff Fessler

if nargin < 3, help osem_sum, error args, end

[nb, na] = size(yi);

if (nargin < 4 || isempty(ci))
	ci = ones(size(yi));
end
if (nargin < 5 || isempty(ri))
	ri = zeros(size(yi));
end
if (nargin < 6 || isempty(block_order))
	block_order = 1;
end
nblock = length(block_order);
if length(block_order) == 1
	nblock = block_order;
	block_order = 1:nblock;
end
if (nargin < 7 || isempty(Asums))
	% compute Asum

	not done!!!!!!!1
	for ii=1:nblock
		ib = block_order(ii);
	end
	Asum = zeros([size(x) nblock]);
end
if (nargin < 7 || isempty(nblock))
	nblock = 1;
end

Gt = G';	% easier to work with transpose

for ib=1:nblock
	ia = ib:nblock:na;
	Asum = zeros(size(x));
	eterm = zeros(size(x));

	%	Asum = G' * ci
	for ia=ib:nblock:na
		Asum = Asum + Gt(:,[1:nb]+(ia-1)*na) * ci(ia,:)';
	end
	return
	yp = ci .* (G * x) + ri;	% predicted measurements
	eterm = G' * (ci .* (yi ./ yp));
	x = x .* eterm ./ Asum;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [Lnew,Q] = bst_remove_silent(L)
% bst_remove_silent: removes silent component of single sphere head model.
%
% USAGE:  [Lnew,Q]=bst_remove_silent(L)
%
% DESCRIPTION:
%     Removes silent component of single sphere head model.
%
% INPUTS:
%     - L : Forward field matrix for all the channels
%
% OUTPUTS:
%     - Lnew : New forward field matrix for all the channels containg 2/3
%     the number of dipole components. The radial dipole components have
%     been removed, and only the gain vectors produced by the two
%     tangential dipole components are included.
%     - Q: A matrix containg the orientations of the two tangential dipole
%     components per source point.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Copyright (C) 2010 - Rey Rene Ramirez
% Authors:  Rey Rene Ramirez, Ph.D.   e-mail: rrramirez at mcw.edu

szL = size(L);
L = reshape(L,[szL(1) 3 szL(2)/3]);
Lnew = zeros(szL(1),szL(2)*(2/3));
sp = 0;
for spoint = 1:2:szL(2)*(2/3)
    sp = sp+1;
    [uL,sL,vL] = svd(L(:,:,sp),'econ');
    Lnew(:,spoint:spoint+1) = uL(:,1:2)*sL(1:2,1:2);
    Q(:,spoint:spoint+1) = vL(:,1:2);
end

author: brainstorm-tools


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

for i=1:length(X)
    distance = inf;
    for j=1:K
        kDist = norm(X(i, :) - centroids(j, :));
        if (kDist < distance)
            distance = kDist;
            idx(i) = j;
        end
    end
end

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

end

author: merwan

repo: ml-class

save_path: github-repos/MATLAB/merwan-ml-class


text:
%% function for get files full name by replacing relative locations like '.', '~'
function nam = get_fullname(nam)
tmp_dir = cd();
[dir_nm, file_nm, ext] = fileparts(nam);
if isempty(dir_nm)
    dir_nm = tmp_dir;
end
cd(dir_nm);
dir_nm = cd();
nam = [dir_nm,filesep, file_nm, ext];
cd(tmp_dir);
end

author: zhoupc


text:
% DEMSWISSROLL1 Model the face swiss roll with a 2-D GPLVM.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'brendan';
experimentNo = 1;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('fitc');
%options.optimiser = 'conjgrad';

latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, dataSetName, experimentNo);
end

% Load the results and display dynamically.
lvmResultsDynamic(model.type, dataSetName, experimentNo, 'image', [20 28], 1, 0, 1)

% Display results
fgplvmScatterPlotColor(model, model.y(:, 2));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [index, window] = FineMeyerWindow(dyadic_points,deg);
% FineMeyerWindow: Evaluates a highpass Meyer window
%  Usage:
%    [index, window] = FineMeyerWindow(L,deg);
%  Inputs:
%    dyadic-points   Pair of the form 2^j, 2^(j+1)
%    deg    Degree of the polynomial
%  Outputs:
%    index  location of the window on the time axis; contains [2^j 2^(j+1)]
%    window highpass Meyer window

pio2 = pi/2;

eps    = floor(dyadic_points(1)/3);
epsp   = dyadic_points(1) - eps - 1;
farlftind = [ 1 : dyadic_points(1)-eps];
lftind     = [ dyadic_points(1)-eps+1 : dyadic_points(1)];
lmidind    = [ dyadic_points(1)+1 : dyadic_points(1)+eps+1];
rmidind    = [ dyadic_points(2)-epsp+1 : dyadic_points(2) ];

farlft = zeros(1,length(farlftind));
lft  = sin(pio2*WindowMeyer(3*((lftind-1)/dyadic_points(2))-1,deg));
lmid = sin(pio2*WindowMeyer(3*((lmidind-1)/dyadic_points(2))-1,deg));
rmid = ones(1,length(rmidind));

index = [farlftind lftind lmidind rmidind];
window = [farlft lft lmid rmid];

author: thomaskuestner


text:
function geometry_test20325 ( )

%*****************************************************************************80
%
%% TEST20325 tests TETRAHEDRON_SIZE_3D, TETRAHEDRON_SHAPE_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 July 2007
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST20325\n' );
  fprintf ( 1, '  For the tetrahedron,\n' );
  fprintf ( 1, '  TETRAHEDRON_SIZE_3D returns dimension information;\n' );
  fprintf ( 1, '  TETRAHEDRON_SHAPE_3D returns face and order info.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We will use this information to compute the\n' );
  fprintf ( 1, '  areas and centers of each face.\n' );

  [ point_num, edge_num, face_num, face_order_max ] = ...
    tetrahedron_size_3d ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of points =   %d\n', point_num );
  fprintf ( 1, '  Number of edges =    %d\n', edge_num );
  fprintf ( 1, '  Number of faces =    %d\n', face_num );
  fprintf ( 1, '  Maximum face order = %d\n', face_order_max );

  [ point_coord, face_order, face_point ] = ...
    tetrahedron_shape_3d ( point_num, face_num, face_order_max );

  shape_print_3d ( point_num, face_num, face_order_max, ...
    point_coord, face_order, face_point );
%
%  Compute the area of each face.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Face  Order  Area\n' );
  fprintf ( 1, '\n' );

  for face = 1 : face_num

    for j = 1 : face_order(face)
      point = face_point(j,face);
      v(1:dim_num,j) = point_coord(1:dim_num,point);
    end

    [ area, normal ] = polygon_area_3d ( face_order(face), v );

    fprintf ( 1, '  %6d  %5d  %8f\n', face, face_order(face), area );

  end
%
%  Find the center of each face.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Face  Center\n' );
  fprintf ( 1, '\n' );

  for i = 1 : face_num

    vave(1:dim_num) = 0.0;

    for j = 1 : face_order(i)
      k = face_point(j,i);
      vave(1:dim_num) = vave(1:dim_num) + point_coord(1:dim_num,k)';
    end

    vave(1:dim_num) = vave(1:dim_num) / face_order(i);

    fprintf ( 1, '  %6d  %10f  %10f  %10f\n', i, vave(1:dim_num) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Chemical master equation for 2H + O <-> H2O

L = 7;  % QTT dimension in state variables

Trange = 0:0.5:5;
d0ts = 12*ones(1,numel(Trange)-1);

tol = 1e-7;
eps = 1e-12;

% initial distribution
l = [13,15,0]; % H, O, H2O

% One shift z=-1 (for O)
S1 = tt_shf(L);
S2 = round(S1*S1, eps); %  double shift z=-2 (for H)
% z=1 for H2O will be S1'

I = tt_eye(2,3*L);
e1 = tt_tensor(num2cell([1;0]*ones(1,L), 1));
e2 = zeros(2^L,1); e2(2)=1; e2 = tt_tensor(reshape(e2, 2*ones(1,L)),eps);

x = tt_x(2,L);
o = tt_ones(2,L);

w1 = 2*mtkron(o-e1-e2, o-e1, o);
w2 = 1*mtkron(o, o, o-e1);
% w = mtkron(x, x, o);

% CME matrix
A = -(mtkron(S2,S1,S2')*diag(w1) - diag(w1)); % forward
A = A - (mtkron(S2',S1',S2)*diag(w2) - diag(w2)); % backward
A = round(A, eps);

% initial state
e1 = zeros(2^L,1); e1(l(1)+1)=1;
u0 = tt_tensor(e1);
e1 = zeros(2^L,1); e1(l(2)+1)=1;
u0 = tkron(u0, tt_tensor(e1));
e1 = zeros(2^L,1); e1(l(3)+1)=1;
u0 = tkron(u0, tt_tensor(e1));
% make QTT
u0 = tt_reshape(u0, 2*ones(1,3*L), eps);

Nt = numel(Trange)-1;
u = u0;
ll = zeros(Nt, 3);
tic;
for t=1:Nt
    d0t = d0ts(t);
    tau = (Trange(t+1)-Trange(t))/(2^d0t);

    Grad_t = IpaS(d0t,-1)/tau;
%     CN_term = IpaS(d0t,1)*0.5;   % Crank-Nicolson
    CN_term = tt_eye(2,d0t);   % Backw Euler
    e1t = tt_tensor(num2cell([1;0]*ones(1,d0t), 1));
    et = tt_ones(2, d0t);

    M = tkron(I, Grad_t)+tkron(A, CN_term);

    rhs = u/tau; %-0.5*A*u;
    rhs = tkron(rhs, e1t);
    U = tkron(u, et);

    U = amen_solve2(M, rhs, tol, 'x0', U);

    % Extract the final solution
    ext = tt_unit(2,d0t,2*ones(d0t,1));
    u = dot(ext, U, 3*L+1, U.d);

    mass = dot(u, tt_ones(u.n))
    u = u/mass;

    u2 = tt_reshape(u, 2^L*ones(1,3));
    % It should be the delta-function
    % Find the final distribution
    [v,l]=tt_stat(u2, 'lm');
    l=l-1
    ll(t,:)=l;
end
toc;

sp = zeros(2^L, 4);
for z=1:2^L
    [v,l]=tt_stat(u2(:,:,z), 'lm');
    if (v>tol)
        sp(z, 1)=l(1)-1;
        sp(z, 2)=l(2)-1;
        sp(z, 3)=z-1;
        sp(z, 4)=log2(v);
    else
        sp(z,:)=NaN;
    end
end
z = z-1;
scatter3(sp(1:z,1), sp(1:z,2), sp(1:z,3), 100, sp(1:z,4), 'filled');
xlabel('H');
ylabel('O');
zlabel('H2O');
title('log2(P)');
colorbar;

author: oseledets

repo: TT-Toolbox


text:
function [cl1,cl2,dat1,dat2] = iimg_cluster_intersect(dat1,dat2,volInfo)
% Prunes two image data vectors (dat1 and dat2) assumed to contain suprathreshold
% contiguous clusters ("blobs") by saving only those blobs that have one or
% more significant (nonzero) elements in both images
%
% :Usage:
% ::
%
%     [cl1,cl2,dat1,dat2] = iimg_cluster_intersect(dat1,dat2,volInfo)
%
% :Inputs:
%
%   **volInfo.xyzlist:**
%        must contain xyz coordinates corresponding to dat1 and dat2
%
%   **dat1 and dat2:**
%        can be either image-length or in-mask length vectorized images
%
% Try iimg_threshold to create dat vectors from Analyze images (e.g.,
% statistic images)
%
% The outputs cl1 and cl2 are the overlap (intersection) clusters, in the
% space of dat1 (cl1) and dat2 (cl2).  The outputs dat1 and dat2 are the
% 'pruned' intersection data vectors.
%
% ..
%    tor wager, july 06
% ..

% prune dat1 with sig values in dat2
[dat1] = iimg_cluster_prune(dat1,dat2,volInfo);

% do the reverse
[dat2] = iimg_cluster_prune(dat2,dat1,volInfo);

% make clusters
cl1 = iimg_indx2clusters(dat1,volInfo);
cl2 = iimg_indx2clusters(dat2,volInfo);

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% RATE OF CONVERGENCE OF MIXED FINITE ELEMENT METHOD IN 2D
%
% This example is to show the rate of convergence of mixed finite element
% (RT0-P0) approximation of the Poisson equation on the unit square:
%
% $$- \Delta u = f \; \hbox{in } (0,1)^2$$
%
% for the following boundary condition:
%
% # Pure Dirichlet boundary condition. $\Gamma _D = \partial \Omega$.
% # Pure Neumann boundary condition. $\Gamma _N = \partial \Omega$.
% # Mix Dirichlet and Neumann boundary condition. $u=g_D \hbox{ on }\Gamma_D,
% \quad \nabla u\cdot n=g_N \hbox{ on }\Gamma_N. \Gamma _N = \{(x,y): x=0,
% y\in [0,1]\}, \; \Gamma _D = \partial \Omega \backslash \Gamma _N$.
%
% Written by Ming Wang.

%%
clear all; close all;
[node,elem] = squaremesh([0,1,0,1],0.25);
pde = mixBCdata;
option.L0 = 2;
option.maxIt = 4;
option.printlevel = 1;
option.elemType = 'RT0-P0';
% option.solver = 'uzawapcg';
option.solver = 'dmg';

%% Pure Dirichlet boundary condition.
% option.plotflag = 0;
bdFlag = setboundary(node,elem,'Dirichlet');
mfemPoisson3(node,elem,pde,bdFlag,option);

%% Pure Neumann boundary condition.
% option.plotflag = 0;
bdFlag = setboundary(node,elem,'Neumann');
mfemPoisson(node,elem,pde,bdFlag,option);

%% Mix Dirichlet and Neumann boundary condition.
option.plotflag = 1;
option.solver = 'uzawapcg';
bdFlag = setboundary(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
mfemPoisson(node,elem,pde,bdFlag,option);

%% Conclusion
%
% The optimal rates of convergence for u and sigma are observed, namely,
% 1st order for L2 norm of u, L2 norm of sigma and H(div) norm of sigma.
% The 2nd order convergent rates between two discrete functions ||uI-uh||
% and ||sigmaI-sigmah|| are known as superconvergence.
%
% dmg and uzawapcg converges uniformly in all cases. Distributive MG is two
% times faster than uzawapcg.

repo: Aero-matlab


text:
function  bz_ConcatenatedTimes(basepath,sortFiles)
%   bz_ConcatenatedTimes(basepath)
%   This function generates the events.mat file that saves all the merge
%   information. Useful if you dind't do this when you merged the dats.
%
%  USAGE
%
%    bz_ConcatenateDats(basepath,sortFiles)
%
%  INPUTS
%
%    basepath          computer path to session folder.  Defaults to
%                      current folder if no input given
%    sortFiles               - boolean denoting whether to sort files according
%                              to time of recording (1) or
%                              not (0) and thus sort them alphabetically
%                              Default = 1.
%
%  OUTPUT
%     Operates on files in specified folder.  No output variable
%
%   Antonio FR 10/2018

%% Handling inputs
if ~exist('basepath','var')
    basepath = cd;
end
basename = bz_BasenameFromBasepath(basepath);

if ~exist('sortFiles','var')
    sortFiles = 1;
end

%% If the dats are not merged quit
if ~exist(fullfile(basepath,[basename,'.dat']),'file')
    disp('no merged dat in directory, quitting')
    return
end

if exist(fullfile(basepath,[basename,'.MergePoints.events.mat']),'file')
    disp('MergedPoints already in directory')
    return
end

%% Find all .dat paths in subfolders

otherdattypes = {'analogin';'digitalin';'auxiliary';'time';'supply'};
bad_otherdattypes = [];
for odidx = 1:length(otherdattypes)
    %eval(['new' otherdattypes{odidx} 'path = fullfile(basepath,''' otherdattypes{odidx} '.dat'');'])
    newpaths.(otherdattypes{odidx}) = fullfile(basepath,[otherdattypes{odidx}, '.dat']);
end

d = dir(basepath);
datpaths = {};
datsizes.amplifier = [];
recordingnames = {};
rcount = 0; %Count of good subfolders
for a = 1:length(d)
    %look in each subfolder
    if d(a).isdir
        %Check for amplifier.dat or subfolderbaseName.dat
        if exist(fullfile(basepath,d(a).name,[d(a).name,'.dat']),'file')
            ampfile = fullfile(basepath,d(a).name,[d(a).name,'.dat']);
        else
            ampfile = fullfile(basepath,d(a).name,'amplifier.dat');
        end

        if exist(ampfile,'file')
            rcount = rcount+1;
            datpaths.amplifier{rcount} = ampfile;
            t = dir(ampfile);
            datsizes.amplifier(rcount) = t.bytes;
            recordingnames{rcount} = d(a).name;

            for odidx = 1:length(otherdattypes)%loop through other .dat types found here
               % eval([otherdattypes{odidx} 'datpaths{rcount} = fullfile(basepath,recordingnames{rcount},''' otherdattypes{odidx} '.dat'');'])
                datpaths.(otherdattypes{odidx}){rcount} = fullfile(basepath,recordingnames{rcount},[otherdattypes{odidx} '.dat']);
                %eval(['d2 = dir(' otherdattypes{odidx} 'datpaths.amplifier{rcount});'])
                d2 = dir(datpaths.(otherdattypes{odidx}){rcount});
                if isempty(d2)
                    bad_otherdattypes(odidx) = 1;
                else
                    %eval([otherdattypes{odidx} 'datsizes(rcount) = d2(1).bytes;'])
                    datsizes.(otherdattypes{odidx})(rcount) = d2(1).bytes;
                end
            end
        end
    end
end
otherdattypes(find(bad_otherdattypes)) = [];%if there weren't analogin or digitalin in some recording
if isempty(datpaths.amplifier)
    disp('No .dats found in subfolders.  Exiting bz_ConcatenateDats.')
    return
end

%% Get the XML
try
    %Look for xml/sessionInfo in topfolder
    sessionInfo = bz_getSessionInfo(basepath,'noPrompts',true);
catch
    %If none exists, look for xml in any of the subpaths
    disp('No .xml or .sessionInfo in top folder, trying subfolders')
    for ff = 1:length(recordingnames)
        try
            sessionInfo = LoadParameters(fullfile(basepath,recordingnames{ff}));
            xmlfilename = fullfile(sessionInfo.session.path,[sessionInfo.session.name,'.xml']);
            [SUCCESS,MESSAGE,MESSAGEID] = copyfile(...
                xmlfilename,...
                fullfile(basepath,[basename,'.xml']),'f');
            display(['Copied xml from ',recordingnames{ff}])
            break
        catch
        end
    end
end

%% Sort files according to time of recording

if sortFiles

    try
        names2sort = cellfun(@(X) str2num(X(end-5:end)),recordingnames,'UniformOutput',false);
        names2sort = cell2mat(names2sort);
        %disp('Assuming the last 6 digits reflect recording time.')
        %disp('Don''t like it? Write in some new options for sorting.')
    catch
        disp('Last 6 digits not numeric... sorting alphabetically')
    end

    [~,I] = sort(names2sort);
    recordingnames = recordingnames(I);
    datpaths.amplifier = datpaths.amplifier(I);
    datsizes.amplifier = datsizes.amplifier(I);
    for odidx = 1:length(otherdattypes)
        datpaths.(otherdattypes{odidx}) = datpaths.(otherdattypes{odidx})(I);
        datsizes.(otherdattypes{odidx}) = datsizes.(otherdattypes{odidx})(I);
    end
else
    disp('not implemented for non temporally sorted concatenation')
end

%% Check that size of resultant .dat is equal to the sum of the components
newdatpath = fullfile(basepath,[basename,'.dat']);
t = dir(newdatpath);
if t.bytes ~= sum(datsizes.amplifier)
    error('New .dat size not right.  Exiting')
    return
else
    sizecheck.amplifier = true;
    disp('Primary .dats concatenated and size checked')
end

%% Get time points from the time.dat
%Use the timestamps from time.dat to get the sort order
%Number of samples in time.dat. First timepoint, last timepoint
%Convert from number of samples to recording time of start/ends
for ff = 1:length(datpaths.time)

	f = fopen(datpaths.time{ff},'r');
    % Determine total number of samples in file
    fileStart = ftell(f);

    %Read the first time point
    firsttimepoint = fread(f,1,'int32');
    status = fseek(f,-4,'eof'); %int32 = 4 bytes
    lasttimepoint = fread(f,1,'int32');
    fileStop = ftell(f);

    firstlasttimepoints(ff,:) = [firsttimepoint lasttimepoint];
    numsamples(ff) = fileStop./4;
    if ff==1
        transitiontimes_samp = firstlasttimepoints(ff,:);
    else
        transitiontimes_samp(ff,:) = firstlasttimepoints(ff,:)+transitiontimes_samp(ff-1,2)+1;
    end
end

disp(['Calculating merge times based on wideband samplingRate of ',num2str(sessionInfo.rates.wideband),'Hz.'])
transitiontimes_sec = transitiontimes_samp./sessionInfo.rates.wideband; %convert to seconds

%% Make the events.mat file that saves all the merge information
% filesT.int == MergePoint.timestamps

eventsfilename = fullfile(basepath,[basename,'.MergePoints.events.mat']);

MergePoints.timestamps = transitiontimes_sec;
MergePoints.timestamps_samples = transitiontimes_samp;
MergePoints.firstlasttimpoints_samples = firstlasttimepoints;
MergePoints.foldernames = recordingnames;
MergePoints.filesmerged = datpaths;
MergePoints.filesizes = datsizes;
MergePoints.sizecheck = sizecheck;
MergePoints.detectorinfo.detectorname = 'bz_ConcatenateDats';
MergePoints.detectorinfo.detectiondate = datestr(now,'yyyy-mm-dd');

save(eventsfilename,'MergePoints');

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function value = ch_is_control ( c )

%*****************************************************************************80
%
%% CH_IS_CONTROL is TRUE if C is a control character.
%
%  Discussion:
%
%    A "control character" has ASCII code <= 31 or 127 <= ASCII code.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the character to be tested.
%
%    Output, logical VALUE, TRUE if C is a control character, and
%    FALSE otherwise.
%
  if ( c <= 31 || 127 <= c )
    value = 1;
  else
    value = 0;
  end

  return
end


text:
function options = dnetOptions(mappingType, latentPoints, mappingOptions)

% DNETOPTIONS Options for a density network.
% FORMAT
% DESC returns the default options for a density network.
% RETURN options : default options structure for density network.
%
% FORMAT
% returns the default options for a density network given a
% number of hidden units.
% ARG mappingType : type of mapping to use in density network.
% ARG latentPoints : number of latent points to use. If give as a vector,
% layout is assumed to be a grid with rows and columns given by the vector.
% ARG mappingOptions : the options for the mapping model.
% RETURN options : default options structure for density network with the
% specified number of hidden units.
%
% SEEALSO : dnetCreate, modelOptions
%
% COPYRIGHT : Neil D. Lawrence, 2008

% MLTOOLS

if nargin < 1
  % This is the GTM settings
  mappingType = 'rbf';
end
if nargin < 2
  latentPoints = [16 16];
end
if nargin < 3
  options = modelOptions(mappingType);
end

options.mappingOptions = options;
options.mappingType = mappingType;
options.latentSamples = latentPoints;

options.alpha = 1;

options.activeFunc = 'linear';
options.optimiser = optimiDefaultOptimiser;

options.initX = 'ppca'; % Initialise X by ppca by default.
if length(latentPoints)>1
  options.M = prod(latentPoints);
  options.grid = latentPoints;
  options.basisStored = true;
else
  options.grid =[];
  options.M = latentPoints;
  options.basisStored = false;
end


text:
function value = c8mat_norm_fro ( m, n, a )

%*****************************************************************************80
%
%% C8MAT_NORM_FRO returns the Frobenius norm of a C8MAT.
%
%  Discussion:
%
%    The Frobenius norm is defined as
%
%      C8MAT_NORM_FRO = sqrt (
%        sum ( 1 <= I <= M ) sum ( 1 <= j <= N ) A(I,J) * A(I,J) )
%
%    The matrix Frobenius norm is not derived from a vector norm, but
%    is compatible with the vector L2 norm, so that:
%
%      c8vec_norm_l2 ( A * x ) <= c8mat_norm_fro ( A ) * c8vec_norm_l2 ( x ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, complex A(M,N), the matrix whose norm is desired.
%
%    Output, real VALUE, the norm of A.
%
  value = ...
    sqrt ...
    ( ...
      sum ...
      ( ...
        sum ...
        ( ...
          ( ...
            abs ...
            ( ...
              a(1:m,1:n) ...
            ) ...
          ) .^ 2 ...
        ) ...
      ) ...
   );

  return
end


text:
%kvmindis 'Minimum Distance Classifier (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vmindis.pane file
%
% Parameters:
% InputFile: i1 'Input Image ', required: 'input image to be classified'
% InputFile: i2 'Input Center/Class Image', required: 'input center/class image'
% Integer: b 'Border Width ', default: 0: ' specifies the border width in pixels (default = 0)'
% OutputFile: o 'Output Classified Image ', required: 'output specifying which vector belongs to which cluster'
%
% Example: o = kvmindis({i1, i2}, {'i1','';'i2','';'b',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vmindis - Minimum Distance Classifier  (K1)
%
%  DESCRIPTION
% .I vmindis
% is a simple minimum distance classifier. The distance metric used
% to the Euclidean distance. This routine accepts two images as the
% input. The image that corresponds to the -i1 argument is the image
% that needs to be classified. This image can be have as many data
% bands as necessary. The number of data bands depicts the dimensionality
% of the vectors. Each pixel in the image is a vector with dimensionality
% of the number of data bands.
%
% The other input image which corresponds to the -i2 argument is
% the prototype image. The last data band of this image is the
% class data band. This image must contain the same number of data
% bands as the other input image plus an extra data band that represents
% the class mapping. This
% image would most likely have been created by vkmeans or some other
% routine that will give cluster centers. This image contains vectors that
% correspond to the prototype of each class.
%
% As stated above the center image's last data band is a class data band. The
% class data band simply maps each vector in the center image to the final class.
%
% At this point most class images must be created manually.
% A class image can be created by using pseudo color in editimage
% on the resulting clustered image from the
% vkmeans routine. Vcustom can then be used to create the class image. Finally
% kinset can be used to add the class image data to the center image. Normally
% one would over cluster using vkmeans then reduce the space down. So vkmeans
% might produce 100 clusters but, really only 3 clusters are desired. Pseudo
% color in editimage and kasc2val can be used to reduce the space down and
% create a class image. The class data might map cluster center vectors 1-50
% to class 1, vectors 51-60 to class 2 and vectors 61-100 to class 3.
% When the vmindis routine is run, the result would be a classified image of
% three classes.
%
% The border option (-b) allows the user to specify a border width,
% in pixels, encompassing the image. The border region is skipped
% by vmindis when classification is performed.
% This useful if neighborhood operators have been used
% previously and have not updated edge pixels.
%
% All input images must be of data storage type FLOAT.
%
%
%
%  EXAMPLES
% vmindis -i1 aerial_image.xv -i2 vkmeans_centers -b 5 -o test
%
% aerial_image.xv is a 5 band image; thus, each vector has dimensionality of 5,
% vkmeans_centers image is a 6 band image, 5 bands are the cluster center values
% produced by vkmeans and 1 class band.
% The border that vmindis will ignore is 5 pixels wide. the output image will be
% a single band image called test.
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% All input images must be of data storage type FLOAT.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef Dataset < handle
%DATASET Parent class for datasets (e.g. datasets.CIFAR10)
%   The datasets.Dataset class makes it easier to partition a dataset into
%   training/validation sets, as well as in mini-batches.
%
%   It can be used either as a base class (e.g. datasets.CIFAR10), or on
%   its own for custom datasets, without creating a subclass. See the
%   'autonn/examples/cnn/' directory for different examples.
%
%   D = datasets.Dataset() creates an empty, default dataset.
%
%   D.train() returns a cell array with mini-batches, from the shuffled
%   training set. Each mini-batch consists of a set of indexes drawn from
%   the vector D.trainSet.
%
%   D.val() returns a cell array with mini-batches, from the validation set
%   (without shuffling). Each mini-batch consists of a set of indexes drawn
%   from the vector D.valSet.
%
%   D.partition(IDX) partitions a given set (vector) IDX into mini-batches.
%
%   datasets.Dataset('option', value, ...) sets the following properties:
%
%   `batchSize`:: 128
%     The batch size.
%
%   `partialBatches`:: false
%     Whether partial batches are returned (which can happen for the last
%     batch in a set, if the batch size does not divide the set size).
%
%   `trainSet`:: []
%     The training set. Note that subclasses (e.g. datasets.CIFAR10)
%     usually override this property with the default training set.
%
%   `valSet`:: []
%     The validation set. Note that subclasses (e.g. datasets.CIFAR10)
%     usually override this property with the default validation set.

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

  properties
    batchSize = 128  % batch size
    trainSet = []  % indexes of training samples
    valSet = []  % indexes of validation samples
    partialBatches = false  % whether to return a partial batch (if the batch size does not divide the dataset size)
  end

  methods
    function dataset = Dataset(varargin)
      % parse generic Dataset arguments on non-subclassed construction.
      % allows Dataset to be used as a stand-alone class.
      varargin = dataset.parseGenericArgs(varargin) ;
      assert(isempty(varargin), 'Unknown arguments.') ;
    end

    function args = parseGenericArgs(o, args)
      % called by subclasses to parse generic Dataset arguments
      args = vl_parseprop(o, args, {'batchSize', 'trainSet', 'valSet', 'partialBatches'}) ;
    end

    function batches = train(o)
      % shuffle training set, and partition it into batches
      assert(~isempty(o.trainSet), ['To use the default dataset.train() ' ...
        'method, the training set must be specified (dataset.trainSet).']) ;

      batches = o.partition(o.trainSet(randperm(end))) ;
    end

    function batches = val(o)
      % partition validation set into batches
      assert(~isempty(o.valSet), ['To use the default dataset.val() ' ...
        'method, the validation set must be specified (dataset.valSet).']) ;

      batches = o.partition(o.valSet) ;
    end

    function batches = partition(o, idx, batchSz)
      % partition indexes into batches (stored in a cell array).
      % if IDX is a matrix, each column is a distinct sample.
      if nargin < 3  % allow overriding batch size
        batchSz = o.batchSize ;
      end
      if isvector(idx)
        idx = idx(:)' ;  % ensure row-vector
      end
      batchSz = min(batchSz, size(idx,2));  % guard against big batch size
      batches = cell(1, ceil(size(idx,2) / batchSz)) ;
      b = 1 ;
      for start = 1 : batchSz : size(idx,2)
        batches{b} = idx(:, start : min(start + batchSz - 1, end)) ;
        b = b + 1 ;
      end

      % delete last partial batch if needed
      if ~o.partialBatches && size(batches{end}, 2) < batchSz
        batches(end) = [] ;
      end
    end
  end

end

repo: MFT


text:
classdef libsvm_CL

%  libsvm_CL is a Support Vector Machine (SVM) classifier object (CL)
%   that uses the LIBSVM package (i.e., this function is Neural Decoding Toolbox
%   wrapper for the LIBSVM package).  A support vector machine is a classifier
%   that learns a function f that minimize the hinge loss between the training
%   data and labels, while also applying a penalty to more complex f (the penalty
%   is based on the norm of f in a RKHS).  The SVM has a parameter C, that controls
%   the trade off between a smaller empirical loss (i.e., smaller prediction error
%   on the training set), a more complex function f.  Support vector machine can use
%   different kernels to create nonlinear decision boundaries (several are supported here).
%
%  SVMs are designed to work on binary classification problems.  In order to
%   support multi-class classification, we use two different methods.  The first method is
%   called 'all-pairs' and works by training separate classifiers for all pairs of labels
%   (i.e., if there are 100 different classes then nchoosek(100, 2) = 4950 different classifiers
%   are trained.  Testing the classifier in all-pairs involves having all classifiers
%   classify the test point, and then the class label is given to the class the was chosen most
%   often by the binary classifiers (in the case of a tie in the number of classes that won a contest
%   the class label is randomly chosen).  The decision values for all-pairs are the number
%   of contests won by each class (for each test point).
%
%  The second multi-class method is called 'one-vs-all' classification
%   and works by training one classifier for each class (thus if there are 100 classes
%   there will be 100 classifiers) using data from one class as the positive labels,
%   and data from all the other classes as the negative labels.  A test point is
%   then run through these different classifiers and the class that has the largest
%   SVM prediction value is returned as the predicted label
%   (i.e., SVMs create a function f(x) = y, and the class label is
%   usually given as sign(y), however here we are comparing the actual y values to determine
%   the label).  The decision values here are the f(x) = y values returned by the SVM. Our limited
%   tests have found all pairs is faster and gives slightly more accuracy results so it is the
%   default (although the decision values might be considered more crude).
%
%
%  Function options and default values:
%
%  The following values can be set to change this classifiers behavior (here we
%  are assuming that svm = libsvm_CL)
%
%  - svm.C: scalar  (default svm.C = 1). This is the 1/regularization constant that determines
%       the trade off between the fit to the training data and the amount of regularization/simplicity
%       of the learned function.  A larger value of C means more emphasis on a better fit to the training data.
%
%  - svm.kernel:  'linear',  'polynomial'/'poly', or 'gaussian'/'rbf' (default 'linear).  The type of kernel used
%       which controls the class of functions that classifier is built from.
%
%       -  If svm.kernel = 'polynomial', then the following options must be set:
%           - svm.poly_degree: scalar (there is no default value so this must be set).  The degree of the of the polynomial.
%           - svm.poly_offset: scalar (default svm.poly_offset = 0);  a constant offset if a polynomntal kernel is used
%
%       -  If svm.kernel = 'gaussian', then the following options must be set:
%           - svm.gaussian_gamma (there is no default value so this must be set).  This controls the fall off
%             of the radial basis function kernel (larger values mean a slower fall off).
%
%  - svm.additional_libsvm_options (default svm.additional_libsvm_options = '').  This allows one to add additional
%     control of the classifiers behavior using a string that is in LIBSVM format.  For more details see:
%     http://www.csie.ntu.edu.tw/~cjlin/libsvm/
%
%  - svm.multiclass_classificaion_scheme (default value is 'all_pairs').  This field can
%    be set to 'all_pairs' or 'one_vs_all' and determines if a one-vs-all or an all-pairs
%    multi-class classification scheme is used, as described above (for binary problems
%    both all-pairs and one-vs-all will return the same predicted labels, although the
%    decision values will differ).
%
% Like all CL objects, there are two main methods, which are:
%
%  1.  cl = train(cl, XTr, YTr)
%       This method takes the training data (XTr, YTr) to train the support vector machine.
%        XTr and XTe are in the form [num_features x num_examples]
%        YTr is in the form [num_examples x 1]
%
%  2.  [predicted_labels decision_values] = test(cl, XTe)
%        This method takes the test data and produces a [num_test_points x 1] vector of
%         predicted labels based on the  function learned from the training set.  It also
%         returns a matrix of dimension [num_test_points x num_classes] of decision values.
%
%
%  Notes:
%   1.  If there is a tie among the decision values, then one of tied the classes
%          is chosen randomly as the predicted label.
%

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

%     This code is part of the Neural Decoding Toolbox.
%     Copyright (C) 2011 by Ethan Meyers (emeyers@mit.edu)
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

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

  properties
        C = 1;  % inverse of the regularization constant - higher values cause more emphasis to be place on the
                %  emperical loss (and cause the classifier use more complex functions, which could lead to overfitting).
        kernel = 'linear';  % the type of kernel used in the SVM.  Choices are 'linear', 'polynomial'/'poly', or 'gaussian'/'rbf'.
        poly_degree = [];  % if a polynomial classifier is used, the degree of the polynomial must be set
        poly_offset = 1;  % a constant offset if a polynomntal kernel is used (default is to use an offset of 1)
        gaussian_gamma = [];  % this controlls the falloff of the Gaussian kernel (larger values mean slower falloff = wider Gaussian functions)
        additional_libsvm_options = '';  % can set this to use additional libsvm options that are not explicitly supported by this wrapper
        multiclass_classificaion_scheme = 'all_pairs';   % options are 'one_vs_all' or 'all_pairs'.
  end

  properties (GetAccess = 'public', SetAccess = 'private')
        labels = [];  % all the unique labels for each class
        model = [];  % the model learned from the training data
  end

    methods

        % constructor
        function cl = libsvm_CL
        end

        function cl = train(cl, XTr, YTr)

            % added sanity check
            if size(YTr, 2) ~= size(XTr, 2)  &&  size(YTr, 1) ~= size(XTr, 2)
               error('Number of columns in YTr, and XTr must be the same (i.e., there must be one and exactly one label for each data point)')
            end

            cl.labels = unique(YTr);

            % create the libsvm arguments for the kernel
            cl.kernel = lower(cl.kernel);
            if strcmp(cl.kernel, 'poly'), cl.kernel = 'polynomial'; end
            if strcmp(cl.kernel, 'rbf'), cl.kernel = 'gaussian'; end

            switch cl.kernel
              case 'linear'
                kernelstring = '-t 0 ';
              case 'polynomial'
                kernelstring = ['-t 1 -g 1 -r ' cl.poly_offset ' -d ' cl.poly_degree ' '];
              case 'gaussian'
                kernelstring = ['-t 2 -g ' cl.gaussian_gamma ' '];
            end

            cstring = ['-c ' num2str(cl.C) ' '];  %  create the libsvm arguments for the (inverse of the) regularization constant
            basicstring = ['-s 0 '];   % use C-SVC  (which is the default anyway)
            probstring = '';  %['-b 1 '];

            paramstring = [basicstring cstring kernelstring probstring cl.additional_libsvm_options];

            if strcmp(cl.multiclass_classificaion_scheme, 'all_pairs')

                cl.model = svmtrain2(YTr, double(XTr'), paramstring);  % renamed the libsvm function svmtrain2 because the Matlab bioinformatics toolbox now has a funciton named svmtrain
                %cl.model = svmtrain(YTr, double(XTr'), paramstring);

            elseif strcmp(cl.multiclass_classificaion_scheme, 'one_vs_all')

                % do one vs. all training
               for iClass = 1:length(cl.labels)

                   reformatted_test_data = [XTr(:, (YTr == cl.labels(iClass))) XTr(:, (YTr ~= cl.labels(iClass)))];
                   curr_labels = zeros(size(YTr));
                   curr_labels(1:length(find(YTr == cl.labels(iClass)))) = 1;
                   the_weight = length(find(YTr ~= cl.labels(iClass)))./length(find(YTr == cl.labels(iClass)));
                   cl.model{iClass} = svmtrain2(curr_labels, double(reformatted_test_data'), paramstring);
                   %cl.model{iClass} = svmtrain(curr_labels, double(reformatted_test_data'), paramstring);

               end
            end

        end

        function [predicted_labels decision_values] = test(cl, XTe)

            junk_fake_xte_labels = ones(size(XTe,2),1);

            if strcmp(cl.multiclass_classificaion_scheme, 'all_pairs')

                [predic_label, junk_fake_accuracy, libsvm_format_all_pairs_decision_vals] = svmpredict(junk_fake_xte_labels, double(XTe'), cl.model);   % edited by Ethan to work with libsvm 2.86 (compared to 2.8)

                predicted_labels = predic_label(:,1);
                if (size(predic_label,2)>1)
                    libsvm_format_all_pairs_decision_vals = predic_label(:,2);
                end

                % Creating the decision values that are the number of 'all pairs' classifiers that are in favor of a given class
                %  (rather than the f(x) values that classifier returns)

                % For each test point, all of the 'all pairs' decision values are in each row of the libsvm_format_all_pairs_decision_vals matrix
                % the order goes:  [(1 vs 2), (1 vs 3), (1 vs 4), ... (1 vs num_classes), (2 vs 3), (2 vs 4)... (2 vs num_classes), ... etc... (num_classes -1 vs num_classes)]
                % Below I reparse this into a matrix of the form:  allpairs_decision_vals(iClass1, iClass2, iTestPoint) = libsvm_decision_value for iTestPoint destinguishing between iClass1 and iClass2
                % more negative values indicate preference for iClass1 and positive values indicate iClass2

                num_labels = length(cl.model.Label);
                end_inds = cumsum((cl.model.nr_class -1):-1:1);
                start_inds = end_inds + 1;
                start_inds = [1 start_inds(1:end-1)];

                start_inds(end + 1) = 2;
                end_inds(end + 1) = 1;

                allpairs_decision_vals(1, :, :) = libsvm_format_all_pairs_decision_vals(:, start_inds(1):end_inds(1))';
                for iLabel = 2:num_labels

                    for iPrevLabel = 1:(iLabel -1)
                        curr_decision_vals(iPrevLabel, :) = squeeze(allpairs_decision_vals(iPrevLabel, iLabel-1, :) .* -1);
                    end
                    allpairs_decision_vals(iLabel, :, :) = [curr_decision_vals; libsvm_format_all_pairs_decision_vals(:, start_inds(iLabel):end_inds(iLabel))'];

                end

                % the final decision values are the number of all-pairs wins that a given class has
                for iPoints = 1:size(allpairs_decision_vals, 3)
                    decision_values(iPoints, :) = (sum(sign((allpairs_decision_vals(:, :, iPoints))), 2)');
                end

                % LibSVM's final result is based on the class that has the most 'all pairs wins' (and for ties it takes the lowest class number)
                % since this can introduce a bias (particularly in the confusion matrix), I will randomly select between classes that are tied by
                % adding a small amount of noise to the number of all pairs wins (I am doing this instead of using the randmax function b/c I want
                % this noise in the decision_values when the rank results and other functions are calculated from the decision values)

                decision_values = decision_values + eps .* rand(size(decision_values));

                [vals ind] = max(decision_values');
                predicted_labels = ind';

                predicted_labels = cl.labels(predicted_labels);

            elseif strcmp(cl.multiclass_classificaion_scheme, 'one_vs_all')

               % do one vs. all testing
               for iClass = 1:length(cl.labels)

                   [predic_label, junk_fake_accuracy, libsvm_format_all_pairs_decision_vals] = svmpredict(junk_fake_xte_labels, double(XTe'), cl.model{iClass});
                   decision_values(:, iClass) = libsvm_format_all_pairs_decision_vals;

                   % % linear f(x)  % the decision values are the same as these when a linear kernel is used...
                   %w = cl.model{iClass}.SVs' * cl.model{iClass}.sv_coef;
                   %linear_pred_vals(:, iClass) = (XTe' * w) - cl.model{iClass}.rho;
                   %all_weights(:, iClass) = w;
                   %all_b(:, iClass) = -1 .* cl.model{iClass}.rho;
               end

               [val ind] = randmax(decision_values');
               predicted_labels = ind';

               predicted_labels = cl.labels(predicted_labels);

            end   % end one vs all

        end  % end test method

    end  % end public methods

end   % end classdef


text:
%% Example
% An example script, it is used to show how to use the arcs
% functions of the KUKA iiwa matlab toolbox
% First run the following script in Matlab
% Then start the client on the KUKA iiwa controller

% Mohammad SAFEEA, 24th of April 2018

close all,clear all;clc;
warning('off');
ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
global t_Kuka;
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  warning('Connection could not be establised, script aborted');
  return;
else

      %% Go to initial configuration
      relVel=0.25; % over ride relative joint velocities

      pos={0, -pi / 180 * 10, 0, -pi / 180 * 100, pi / 180 * 90,pi / 180 * 90, 0};   % initial cofiguration

      movePTPJointSpace( t_Kuka , pos, relVel); % go to home position
      %% Move in an arc, the orientation of EEF changes while performing the motion,
      % The function utilized (movePTPCirc1OrintationInter)
      % f2 is the final frame, at which the arc motion ends
      % f1 is an intermidiate frame, through wich the robot passes while
      % performing the motion.

      f1=getEEFPos( t_Kuka );
      f2=f1;
      r=75;
      f1{2}=f1{2}+r;
      f1{3}=f1{3}-r;
      f1{6}=f1{6}+pi/8;

      f2{3}=f2{3}-2*r;
      f2{6}=f2{6}+pi/2;

      vel=150; % linear velocity of end-effector mm/sec
      movePTPCirc1OrintationInter( t_Kuka , f1,f2, vel)

          %% Move robot in joint space to some initial configuration
pinit={0,pi*20/180,0,-pi*70/180,0,pi*90/180,0}; % joint angles of initial confuguration
relVel=0.15; % the relative velocity
movePTPJointSpace( t_Kuka , pinit, relVel); % point to point motion in joint space

%% Move EEF -100 mm in Z direction
deltaX=0.0;deltaY=0;deltaZ=-100.; % relative displacemnets of end-effector
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;
movePTPLineEefRelBase( t_Kuka , Pos, vel);

%% Store the current position in the memory
Cen=getEEFPos(t_Kuka); % Concider the current position as the center of the arcs

%% Move EEF 50mm in X direction
deltaX=100;deltaY=0;deltaZ=0.;
Pos{1}=deltaX;
Pos{2}=deltaY;
Pos{3}=deltaZ;
movePTPLineEefRelBase( t_Kuka , Pos, vel);
%% Store the current position in the memory
circle_Starting_Point=getEEFPos( t_Kuka ); % Consider the current point as circle starting point

%% Move in an arc, the arc is drawn on an incliend plane
% using the function ((movePTPArc_AC))
theta=pi/2; % the angle subtended by the arc at the center ((c))
k=[1;1;1]; % normal vector of the plane, on which the circle is drawn
c=[Cen{1};Cen{2};Cen{3}]; % the center of the arc
vel=100; % the motion velocity mm/sec
movePTPArc_AC(t_Kuka,theta,c,k,vel)

%% Go back to ((circle_Starting_Point)) coordinates
vel=150;
movePTPLineEEF( t_Kuka , circle_Starting_Point, vel);
%% Move in an arc, the arc is drawn in XY plane
% using the function ((movePTPArcXY_AC))
theta=1.98*pi; % the angle subtended by the arc at the center ((c))
c=[Cen{1};Cen{2}]; % the XY coordinate of the center of the arc
vel=150; % the motion velocity mm/sec
movePTPArcXY_AC(t_Kuka,theta,c,vel)

      %% turn off the server
       net_turnOffServer( t_Kuka );

       fclose(t_Kuka);
end

warning('on');

repo: KST-Kuka-Sunrise-Toolbox


text:
%% Harmonic Representation of Rotational Functions
%

%%
% Similarly as periodic functions may be represented as weighted sums of
% sines and cosines a rotational function $f\colon \mathcal{SO}(3)\to\mathbb C$
% can be written as a series of the form
%
% $$ f({\bf R}) = \sum_{n=0}^N \sum_{k,l = -n}^n \hat f_n^{k,l} \, \mathrm{D}_n^{k,l}({\bf R}) $$
%
% with respect to Fourier coefficients $\hat f_n^{k,l}$ and the so called
% <WignerFunctions.html Wigner-D functions> $D_n^{k,l}$.
%
% There exists various normalizations for the <WignerFunctions.html Wigner-D functions>.
% In MTEX they are $L_2$ normalized, which means
%
% $$\| D_n^{k,l} \|_2 = 1$$
%
% for all $n,k,l$. For more information take a look on
% <WignerFunctions.html Wigner-D functions> and
% <SO3FunOperations.html#6 Integration of SO3Fun's>.
%
%%
%
% We construct an arbitrary ODF which generally is an SO3Fun:
mtexdata dubna
odf = calcODF(pf,'resolution',5*degree,'zero_Range')
%%
% Now we may transform an arbitrary SO3Fun into its Fourier representation
% using the command <SO3FunHarmonic.SO3FunHarmonic.html SO3FunHarmonic>

f = SO3FunHarmonic(odf,'bandwidth',32)

%% Fourier Coefficients
%
% Within the class |@SO3FunHarmonic| rotational functions are represented by
% their complex valued Fourier coefficients which are stored in the field
% |fun.fhat|.
% They are stored in a linear order, which means |f.fhat(1)| is the
% zero order Fourier coefficient, |f.fhat(2:10)| are the first order
% Fourier coefficients that form a 3x3 matrix and so on.
% Accordingly, we can extract the second order Fourier coefficients by

reshape(f.fhat(11:35),5,5)

%%
% As an additional example lets define a harmonic function by its Fourier
% coefficients $\hat f_0^{0,0} = 0.5$ and
% $\hat f_1 = \begin{array}{rrr}
% 1 & 4 & 7 \\
% 2 & 5 & 8 \\
% 3 & 6 & 9 \\
% \end{array}$

f2 = SO3FunHarmonic([0.5,1:9]')

plot(f2)
%%
% The Fourier coefficients $\hat f_n^{k,l}$ allow us a complete
% characterization of the rotational function. They are of particular
% importance for the calculation of mean macroscopic properties e.g.
% the second order Fourier coefficients characterize thermal expansion,
% optical refraction index, and electrical conductivity whereas the
% fourth order Fourier coefficients characterize the elastic properties
% of the specimen.
%
% Moreover, the decay of the Fourier coefficients is directly related to
% the smoothness of the SO3Fun. The decay of the Fourier coefficients might
% also hint for the presents of a ghost effect. See
% <PoleFigure2ODFGhostCorrection.html Ghost Correction>.

%%
% The decay of the Fourier coefficients is shown in the plot
close all;
plotSpektra(f)

%% ODFs given by Fourier coefficients
%
% In order to define an ODF by it *Fourier coefficients* ${\bf \hat{f}}$,
% they has to be given as a literally ordered, complex valued
% vector of the form
%
% $$ {\bf \hat{f}} = [\hat{f}_0^{0,0},\hat{f}_1^{-1,-1},\ldots,\hat{f}_1^{1,1},\hat{f}_2^{-2,-2},\ldots,\hat{f}_N^{N,N}] $$
%
% where $n=0,\ldots,N$ denotes the order of the Fourier coefficients.

cs   = crystalSymmetry('1');    % crystal symmetry
fhat = [1;reshape(eye(3),[],1);reshape(eye(5),[],1)]; % Fourier coefficients
odf = SO3FunHarmonic(fhat,cs)

plot(odf,'sections',6,'silent','sigma')

%%

plotPDF(odf,[Miller(1,0,0,cs),Miller(1,1,0,cs)],'antipodal')

%% TODO: Add some non ODF example for an SO3Fun
%
%

author: mtex-toolbox

repo: mtex

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


text:
function a = complex_i ( )

%*****************************************************************************80
%
%% COMPLEX_I returns the COMPLEX_I matrix.
%
%  Formula:
%
%    0 1
%   -1 0
%
%  Properties:
%
%    A is integral: int ( A ) = A.
%
%    A is anti-involutional: A * A = - I
%
%    A * A * A * A = I
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A(2,2), the matrix.
%
  a(1,1) =  0.0;
  a(1,2) =  1.0;
  a(2,1) = -1.0;
  a(2,2) =  0.0;

  return
end


text:
function [mergePsi, logQ] = getMergeConfig_SeqAllocation( Psi, data_struct, hyperparams, model, anchorObjIDs, featIDs, TargetPsi )
% Merge Move launch
%   MERGE all blocks assigned to either featIDs into single feature kmerge
if ~exist( 'TargetPsi', 'var' )
    TargetPsi = [];
end

if isfield( Psi, 'theta' );
    Psi = rmfield( Psi, 'theta');
    Psi = rmfield( Psi, 'TS');
end

logQ.F   = 0;
logQ.z   = 0;
obsModel = model.obsModel;

F = Psi.F;
nObj = size(F,1);
ii = anchorObjIDs(1);
jj = anchorObjIDs(2);

allObjIDs = find( sum( F(:,featIDs), 2 ) > 0 )';
otherObjIDs = setdiff( allObjIDs, anchorObjIDs );
otherFeatIDs = setdiff( 1:size(F,2), featIDs );

propF = F;
propF(:, featIDs) = 0;
kmerge = size(propF,2)+1;

propFeatIDs = [kmerge];
propF(allObjIDs, kmerge) = 1;

propStateSeq = Psi.stateSeq;
for aa = anchorObjIDs
    aIDs =  propStateSeq(aa).z == featIDs(1);
    bIDs =  propStateSeq(aa).z == featIDs(2);
    propStateSeq(aa).z( aIDs | bIDs ) = kmerge;
end
propTS = getEta_PriorMean( propF, [], hyperparams, model, allObjIDs, propFeatIDs );

% Ustats  = getStateSeqSuffStats_C1( propStateSeq, propF, data_struct, model, 1:nObj, [otherFeatIDs propFeatIDs] );
% myTheta.logp = zeros( size(Ustats.Nkv)  );
% myTheta = getTheta_PosteriorMean( myTheta, Ustats, obsModel, [otherFeatIDs propFeatIDs] );

[Xstats] = getXSuffStats( propF, propStateSeq, data_struct, model, 1:nObj, [otherFeatIDs propFeatIDs] );
myTheta = setThetaToPosteriorMean( [], Xstats, obsModel, [otherFeatIDs propFeatIDs] );

% Keep *only* the params for the feature kM being merged
%   this is all we update
AXstats  = getXSuffStats( propF, propStateSeq, data_struct, model, [ii jj], [propFeatIDs] );

if exist( 'TargetPsi', 'var' ) && ~isempty( TargetPsi )
    TargetPsi.externalFeatIDs = kmerge;
end
for aa = [ otherObjIDs( randperm( length(otherObjIDs) ) ) anchorObjIDs]
    if aa == ii || aa == jj
        Xdelta = getXSuffStats( propF, propStateSeq, data_struct, model, aa, propFeatIDs );
        AXstats = decXStats( AXstats, Xdelta, model.obsModel );
    end

    % Construct soft evidence
    activeFeatIDs = find( propF(aa,:)==1 );
    LL = calcLogCondPrObsGivenTheta( data_struct(aa), myTheta, obsModel, activeFeatIDs, size(propF,2) );
    LL = LL( activeFeatIDs, :);
    % Propose new State Seq. for current item
    [propStateSeq(aa).z, logPrQ_z] = sampleStateSeq_PrecompSoftEv( aa, propTS.obj(aa), LL, TargetPsi);
    %[propStateSeq, logQaa] = sampleZ_RGS(data_struct, propTS, propF, myTheta, model.obsModel, propStateSeq, aa, TargetPsi);
    logQ.z = logQ.z + logPrQ_z;

    % Update Emission Params for new proposed features
    %   using the newly allocated data from the current sequence
    Xdelta = getXSuffStats( propF, propStateSeq, data_struct, model, aa, propFeatIDs );
    AXstats = incXStats( AXstats, Xdelta, model.obsModel );
    myTheta = setThetaToPosteriorMean( myTheta, AXstats, obsModel, propFeatIDs );
end

mergePsi.F = propF;
mergePsi.stateSeq = propStateSeq;
mergePsi.activeFeatIDs = propFeatIDs;
mergePsi.anchorObjIDs  = anchorObjIDs;

author: michaelchughes

repo: NPBayesHMM


text:
%% Machine Learning Online Class
%  Exercise 1: Linear regression with multiple variables
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear regression exercise.
%
%  You will need to complete the following functions in this
%  exericse:
%
%     warmUpExercise.m
%     plotData.m
%     gradientDescent.m
%     computeCost.m
%     gradientDescentMulti.m
%     computeCostMulti.m
%     featureNormalize.m
%     normalEqn.m
%
%  For this part of the exercise, you will need to change some
%  parts of the code below for various experiments (e.g., changing
%  learning rates).
%

%% Initialization

%% ================ Part 1: Feature Normalization ================

%% Clear and Close Figures
clear ; close all; clc

fprintf('Loading data ...\n');

%% Load Data
data = load('ex1data2.txt');
X = data(:, 1:2);
y = data(:, 3);
m = length(y);

% Print out some data points
fprintf('First 10 examples from the dataset: \n');
fprintf(' x = [%.0f %.0f], y = %.0f \n', [X(1:10,:) y(1:10,:)]');

fprintf('Program paused. Press enter to continue.\n');
pause;

% Scale features and set them to zero mean
fprintf('Normalizing Features ...\n');

[X mu sigma] = featureNormalize(X);

% Add intercept term to X
X = [ones(m, 1) X];

%% ================ Part 2: Gradient Descent ================

% ====================== YOUR CODE HERE ======================
% Instructions: We have provided you with the following starter
%               code that runs gradient descent with a particular
%               learning rate (alpha).
%
%               Your task is to first make sure that your functions -
%               computeCost and gradientDescent already work with
%               this starter code and support multiple variables.
%
%               After that, try running gradient descent with
%               different values of alpha and see which one gives
%               you the best result.
%
%               Finally, you should complete the code at the end
%               to predict the price of a 1650 sq-ft, 3 br house.
%
% Hint: By using the 'hold on' command, you can plot multiple
%       graphs on the same figure.
%
% Hint: At prediction, make sure you do the same feature normalization.
%

fprintf('Running gradient descent ...\n');

% Choose some alpha value
alpha = 0.01;
num_iters = 400;

% Init Theta and Run Gradient Descent
theta = zeros(3, 1);
[theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters);

% Plot the convergence graph
figure;
plot(1:numel(J_history), J_history, '-b', 'LineWidth', 2);
xlabel('Number of iterations');
ylabel('Cost J');

% Display gradient descent's result
fprintf('Theta computed from gradient descent: \n');
fprintf(' %f \n', theta);
fprintf('\n');

% Estimate the price of a 1650 sq-ft, 3 br house
% ====================== YOUR CODE HERE ======================
% Recall that the first column of X is all-ones. Thus, it does
% not need to be normalized.
price = 0; % You should change this

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

fprintf(['Predicted price of a 1650 sq-ft, 3 br house ' ...
         '(using gradient descent):\n $%f\n'], price);

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ================ Part 3: Normal Equations ================

fprintf('Solving with normal equations...\n');

% ====================== YOUR CODE HERE ======================
% Instructions: The following code computes the closed form
%               solution for linear regression using the normal
%               equations. You should complete the code in
%               normalEqn.m
%
%               After doing so, you should complete this code
%               to predict the price of a 1650 sq-ft, 3 br house.
%

%% Load Data
data = csvread('ex1data2.txt');
X = data(:, 1:2);
y = data(:, 3);
m = length(y);

% Add intercept term to X
X = [ones(m, 1) X];

% Calculate the parameters from the normal equation
theta = normalEqn(X, y);

% Display normal equation's result
fprintf('Theta computed from the normal equations: \n');
fprintf(' %f \n', theta);
fprintf('\n');

% Estimate the price of a 1650 sq-ft, 3 br house
% ====================== YOUR CODE HERE ======================
price = 0; % You should change this

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

fprintf(['Predicted price of a 1650 sq-ft, 3 br house ' ...
         '(using normal equations):\n $%f\n'], price);

author: Ayatans

repo: Machine-Learning-homework


text:
%DEMO_HURDLE Demonstration of Logit Negative-binomial hurdle model
%            using Gaussian process prior
%
%  Description
%    Hurdle models can be used to model excess number of zeros
%    compared to usual Poisson and negative binomial count models.
%    Hurdle models assume a two-stage process, where the first
%    process determines whether the count is larger than zero, and
%    the second process determines the non-zero count. Both
%    processes are given a zero mean Gaussian process prior. The
%    two stage model formulation makes it possible to make
%    inference for the two latent processes separately.
%
%    In this demo we construct logit negative binomial hurdle model
%    by using logit and zero truncated negative binomial models.
%    The posterior inference is made with parallel-EP
%    approximation.
%
%    See also  DEMO_SPATIAL2, DEMO_CLASSIFIC1

% Copyright (c) 2008-2010 Jarno Vanhatalo
% Copyright (c) 2010-2011 Aki Vehtari
% Copyright (c) 2011 Jaakko Riihimäki

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

% load the data
S = which('demo_spatial1');
data = load(strrep(S,'demo_spatial1.m','demodata/spatial1.txt'));

x = data(:,1:2);
ye = data(:,3);
y = data(:,4);

x0=x;
x=bsxfun(@rdivide,bsxfun(@minus,x0,mean(x0)),std(x0));

% for zero-process
yz=double(y>0)*2-1;
% for count-process
ci=find(y>0);
yc=y(ci);
xc=x(ci,:);
yec=ye(ci,:);

% Create the covariance functions
pl = prior_t('s2',10);
pm = prior_sqrtunif();
cf = gpcf_matern32('lengthScale', 5, 'magnSigma2', 0.05, ...
                   'lengthScale_prior', pl, 'magnSigma2_prior', pm);

% Create the zero part
likz=lik_probit();
gpz=gp_set('lik',likz,'cf',cf,'jitterSigma2',1e-6,'latent_method','EP','latent_opt',struct('parallel','on'));
% Create the count part
likc=lik_negbinztr();
gpc=gp_set('lik',likc,'cf',cf,'jitterSigma2',1e-6,'latent_method','EP','latent_opt',struct('parallel','on'));
% Note that although both parts use 'cf', they have separate hyperparameters

% Set the options for the quasi-Newton optimization
opt=optimset('TolFun',1e-2,'TolX',1e-2,'Display','iter');
gpz=gp_optim(gpz,x,yz,'opt',opt,'optimf',@fminlbfgs);
gpc=gp_optim(gpc,xc,yc,'z',yec,'opt',opt,'optimf',@fminlbfgs);

% make prediction to the data points
[Efz, Varfz] = gp_pred(gpz, x, yz, x);
[Efc, Varfc] = gp_pred(gpc, xc, yc, xc, 'z', yec);

% Define help parameters for plotting
xii=sub2ind([60 35],x0(:,2),x0(:,1));
[X1,X2]=meshgrid(1:35,1:60);

% Plot the figures
figure
subplot(1,2,1)
G=repmat(NaN,size(X1));
G(xii)=Efz(:);
pcolor(X1,X2,G),shading flat
colorbar
axis equal
axis([0 35 0 60])
title('Posterior mean of latent zero process')

subplot(1,2,2)
G=repmat(NaN,size(X1));
G(xii(ci))=Efc(:);
pcolor(X1,X2,G),shading flat
colorbar
axis equal
axis([0 35 0 60])
title('Posterior mean of latent count process')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Population = EnvironmentalSelection(Population,N,kappa)
% The environmental selection of IBEA

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    Next = 1 : length(Population);
    [Fitness,I,C] = CalFitness(Population.objs,kappa);
    while length(Next) > N
        [~,x]   = min(Fitness(Next));
        Fitness = Fitness + exp(-I(Next(x),:)/C(Next(x))/kappa);
        Next(x) = [];
    end
    Population = Population(Next);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%--- help for contains ---
%
% CONTAINS True if pattern is found in text.
%    TF = CONTAINS(STR,PATTERN) returns 1 (true) if STR contains PATTERN,
%    and returns 0 (false) otherwise.
%
%    STR can be a string array, a character vector, or a cell array of
%    character vectors. So can PATTERN. PATTERN and STR need not be the same
%    size. If PATTERN is a string array or cell array, then CONTAINS returns
%    true if it finds any element of PATTERN in STR. If STR is a string
%    array or cell array, then TF is a logical array that is the same size.
%
%    TF = CONTAINS(STR,PATTERN,'IgnoreCase',IGNORE) ignores case when searching
%    for PATTERN in STR if IGNORE is true. The default value of IGNORE is false.
%
%    Examples
%        STR = "data.tar.gz";
%        P = "tar";
%        contains(STR,P)                   returns  1
%
%        STR = ["abstracts.docx","data.tar.gz"];
%        P = 'tar';
%        contains(STR,P)                   returns  [0 1]
%
%        STR = 'data.tar.gz';
%        P = {'docx','tar'};
%        contains(STR,P)                   returns  1
%
%        STR ={'DATA.TAR.GZ','SUMMARY.PPT'};
%        P = "tar";
%        contains(STR,P,'IgnoreCase',true) returns  [1 0]
%
%    See also endsWith, startsWith.
%
%    Reference page in Doc Center
%       doc contains
%
%

author: jmaih


text:
function assign(X,value,ls)
%Obsolete, see WARMSTART

if nargin<3
    ls = 0;
end

if ~isa(X,'sdpvar')
    error('First argument should be an SDPVAR object.');
end

if isa(value,'logical')
    value = double(value);
end

if ~isnumeric(value)
    error('Second argument should be NUMERIC.');
end

if prod(size(value)) == 1
    value = repmat(value,size(X));
end

if isempty(value)
    return
end

if min(size(X))==1 && min(size(value))==1 && (length(size(X))==2) && (length(size(value))==2)
    X = reshape(X,[],1);
    value = reshape(value,[],1);
end

if ~isequal(size(X),size(value))
    error('Both arguments must have same size')
end
if ~isa(X,'sdpvar')
    error('First arguments must be an sdpvar object')
end

if is(X,'complex')
    assign(real(X),real(value));
    assign(imag(X),imag(value));
    return
end

x_lmi_variables = X.lmi_variables;
b = value(:)-X.basis(:,1);
A = X.basis(:,2:end);
feas_var = A\b;
% Improve
e = A*feas_var-b;
de = A\e;
feas_var = feas_var-de;

if ~ls
    if norm(A*feas_var-b)>sqrt(eps)
        error('Inconsistent assignment')
    end
end

sol = yalmip('getsolution');
keep_these = find(~ismember(sol.variables,x_lmi_variables));
sol.optvar = [sol.optvar(keep_these);feas_var(:)];
sol.variables = [sol.variables(keep_these);x_lmi_variables(:)];
yalmip('setallsolution',sol);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%

sigmo = sigmoid(X * theta);

J = 1 / m * (-y' * log(sigmo) - (1 - y') * log(1 - sigmo));
grad = 1 / m * ((sigmo - y)' * X)';

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

end

author: benoitvallon

repo: coursera-machine-learning


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
% patch_color.m
%
% Goal: surface rendering with color
%       need to open figure or subplot before use
%
% Li Shen
% 02/21/2006 - create

function patch_color(vertices,faces,cdata)

patches = patch('faces', faces, 'vertices', vertices, 'FaceVertexCData', cdata, ...
		'FaceColor', 'interp', 'EdgeColor', 'k', 'FaceAlpha', 1); % transparancy can be adjusted by FaceAlpha

c = [1 1 1]*0.5;
light('Position',[ 1  0 0],'Style','infinite', 'Color', c);
light('Position',[-1  0 0],'Style','infinite', 'Color', c);
light('Position',[ 0  1 0],'Style','infinite', 'Color', c);
light('Position',[ 0 -1 0],'Style','infinite', 'Color', c);
lighting gouraud

axis image; box on; view(3);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
 function ob = Gdown(Nd, varargin)
%function ob = Gdown(Nd, varargin)
%|
%| Construct Gdown object, which performs downsampling.
%| This is useful perhaps for iterative zooming methods.
%| See Gdown_test.m for example usage.
%|
%| in
%|	Nd	[1 D]		input signal dimensions
%|
%| options
%|	mask	logical			default: true(Nd)
%|	down	1|2|3|...		down sampling factor (default 2)
%|	type	func|Gsparse|...	not done.  default: 'func'
%| out
%|	ob			fatrix2 object
%|
%| Copyright 2006-8-25, Jeff Fessler, University of Michigan

if nargin == 1 && streq(Nd, 'test'), Gdown_test, return, end
if nargin < 1, ir_usage, end

% defaults
arg.mask = [];
arg.idim = Nd;
arg.down = 2;
arg.type = 'func';
arg.class = 'fatrix2';

% options
arg = vararg_pair(arg, varargin);

arg.odim = arg.idim ./ arg.down;
tmp = round(arg.odim) * arg.down;
if any(tmp ~= arg.idim)
	error 'only integer multiple image size supported'
end

if length(arg.idim) > 2
	error 'only 2d implemented due to downsample2() limitation'
end

if isempty(arg.mask)
	arg.mask = true(arg.idim);
end

arg.up_scale = 1 / arg.down^numel(arg.idim);

switch arg.class
case 'Fatrix'
	arg.np = sum(arg.mask(:));
	arg.dim = [prod(arg.odim) arg.np]; % nd x np
	ob = Fatrix(arg.dim, arg, ...
		'abs', @(ob) ob, ...
		'forw', @Gdown_forw_Fatrix, 'back', @Gdown_back_Fatrix);
case 'fatrix2'
	forw = @(arg, x) downsample2(x, arg.down);
	back = @(arg, y) fatrix2_maskit(arg.mask, ...
			upsample_rep(arg.up_scale * y, arg.down));
	ob = fatrix2('arg', arg, 'forw', forw, 'back', back, ...
		'abs', @(ob) ob);

otherwise
	fail 'not done'
end

% Gdown_forw_Fatrix(): y = A * x
% in
%	x	[np L] or [(Nd) L]
% out
%	y	[M L]
%
function y = Gdown_forw_Fatrix(arg, x)

[x ei] = embed_in(x, arg.mask, arg.np);

LL = size(x, 1+length(arg.idim)); % *L
if LL == 1
	y = downsample2(x, arg.down);
else
	y = [];
	for ll=1:LL
		y(:,:,ll) = downsample2(x(:,:,ll), arg.down); % fix: generalize
	end
end

y = ei.shape(y);

% Gdown_back_Fatrix(): x = A' * y
% in
%	y	[M L]
% out
%	x	[np L]
%
function x = Gdown_back_Fatrix(arg, y)

[y eo] = embed_out(y, arg.odim);

LL = size(y, 1+length(arg.odim)); % *L
if LL == 1
	x = upsample_rep(y, arg.down);
else
	for ll=1:LL
		x(:,:,ll) = upsample_rep(y(:,:,ll), arg.down); % fix: generalize
	end
end
x = x * arg.up_scale;

x = eo.shape(x, arg.mask, arg.np);

%
% Gdown_test
%
function Gdown_test

dim = [4 6];
down = 2;
A = Gdown(dim, 'down', down);

x = reshape(1:prod(dim), dim);
y1 = downsample2(x, down);
y2 = A * x;
jf_equal(y1, y2)

y = y1;
x1 = upsample_rep(y, down) * A.arg.up_scale;
x2 = A' * y;
jf_equal(x1, x2)

mask = true(dim);
if isa(A, 'Fatrix')
	Fatrix_test_basic(A, mask)
else
	fatrix2_tests(A)
end

test_adjoint(A);

author: JeffFessler

repo: mirt