text:
function [] = aps_spectrometer_PWV_meris(batchfile)
% aps_spectrometer_PWV_comparison(datelist)
% Scipt to load meris data, mask out clouds. The meris data is assumed to
% be structured in date folders. The batchfile contains the full path to the
% meris files in these folders. Note that the first line of the batchfile should read "files".
%
%
% INPUTS:
% batchfile             A txt file containing the full path and file names of the
%                       meris data that needs to be processed. The first
%                       line of this file should read "files". The data
%                       should be structured in date folders.
% xlims                 Limits in the x-direction, either in degrees
% ylims                 Limits in the y-direction, either in degrees
% smpres                The output resolution, either in degrees
%                       Units needs to be consistend with xlims and ylims.
%
%
% By David Bekaert - University of Leeds
% August 2014

fig_test = 1;                   % when 1 show the dem as debug figure

if nargin<1
    fprintf('aps_spectrometer_PWV_meris(batchfile) \n')
    error('myApp:argChk', ['Not enough input arguments...\n'])
end

smpres = getparm_aps('region_res',1); % in degrees
xlims = getparm_aps('region_lon_range',1);
ylims = getparm_aps('region_lat_range',1);
stamps_processed = getparm_aps('stamps_processed',1);

%% the actual scripting
%bounds for all ifgs in degrees or in meters
xmin = xlims(1);
xmax = xlims(2);
ymin = ylims(1);
ymax = ylims(2);

% getting the number of files to be processed
files = char(textread(batchfile,'%s','headerlines',1));
ndates = size(files,1);

% loading the date information
if strcmp(stamps_processed,'y')
   ps = load(getparm_aps('ll_matfile',1));
   ifgs_dates = ps.day;
   fprintf('Stamps processed structure \n')
else
    ifgday_matfile = getparm_aps('ifgday_matfile',1);
    ifgs_dates = load(ifgday_matfile);
    ifgs_dates = ifgs_dates.ifgday;
    ifgs_dates = reshape(ifgs_dates,[],1);
    ifgs_dates = unique(ifgs_dates);
end

fprintf(['Lon range: ' num2str(xmin) ' -- '  num2str(xmax) ' degrees\n'])
fprintf(['Lat range: ' num2str(ymin) ' -- '  num2str(ymax) ' degrees\n'])
fprintf(['Output resolution is assumed to be ' num2str(smpres) ' degrees \n'])

% extracting the dates from the filenames
for k=1:ndates
    [path,filename_temp,ext_temp] = fileparts(files(k,:));
    clear filename_temp ext_temp
    [path_temp,date,ext_temp] = fileparts(path);
    clear path_temp ext_temp

    % save the paths as structures to allow for variable path lengths
    pathlist{k} = path;
    clear path
    % saving the date information
    datelist(k,:) =date;
    clear date
end

%start loop here to calculate atmos correction for each date
fprintf('Starting the masking and writing of the PWV for each SAR date \n')
for n = 1:ndates

    file = [files(n,:)];
    outfile_watervapor = [pathlist{n} filesep datelist(n,:) '_ZPWV_nointerp.xyz'];

    % get info on meris file from tif
    tifinfo = geotiffinfo(file);
    mer_xmin = tifinfo.CornerCoords.X(1);
    mer_xmax = tifinfo.CornerCoords.X(2);
    mer_ymin = tifinfo.CornerCoords.Y(3);
    mer_ymax = tifinfo.CornerCoords.Y(1);
    mer_cols = tifinfo.Width;
    mer_rows = tifinfo.Height;
    mer_x = (mer_xmax - mer_xmin)/mer_cols;
    mer_y = (mer_ymax - mer_ymin)/mer_rows;
    % read in tif
    meris=imread(file,'tif');

    %% Calculate wet delay if desired

    %do masking, interpolating etc here.
    %convert flag matrix to long list of binary numbers e.g. 1000010000010000
    flagbin = dec2bin(meris(:,:,33));

    %% Generating cloud mask
    %take certain bits of binary numbers (e.g. second bit, 1=cloud), and reshape
    %into matrix mask
    cloud = bin2dec(flagbin(:,2));
    cloudmat=reshape(cloud,mer_rows,mer_cols);
    cloudmat(cloudmat==1)=NaN;

    % turn ones (i.e. mask pixels) into NaNs so when we multiply masks together,
    % NaNs will penetrate through the stack
    pconf = bin2dec(flagbin(:,23));
    pconfmat=reshape(pconf,mer_rows,mer_cols);
    pconfmat(pconfmat==1)=NaN;
    lowp = bin2dec(flagbin(:,24));
    lowpmat=reshape(lowp,mer_rows,mer_cols);
    lowpmat(lowpmat==1)=NaN;
    maskmat= cloudmat .*pconfmat .*lowpmat;
    maskmat(isnan(maskmat))=1;

    %extend mask edges
    se = strel('square', 3);
    maskmat = imdilate(maskmat,se);
    maskmat(maskmat==1)=NaN;
    maskmat(maskmat==0)=1;

    %% loading and masking of the data
    %load water vapour
    watervap=meris(:,:,14);   % g/cm^2  water vapour
    %mask out dodgy pixels
    corwatervap = watervap .*maskmat;

    %% mask water additional
    %%% aditonal masking
    % water mask
    water = bin2dec(flagbin(:,3));
    water=reshape(water,mer_rows,mer_cols);

    %decrease water mask edges, to make sure we do not remove land points
    se = strel('square', 3);
    water_new = imerode(water,se);
    water_new(water_new==1)=NaN;
    water_new(water_new==0)=1;
    clear water
    corwatervap = corwatervap.*water_new;

    %% saving the data
    fid = fopen('out.bin','w');
    corwatervap = corwatervap';
    fwrite(fid,corwatervap,'real*4');
    fclose(fid);

    % do gmt resample and cut
    xyz2grd_cmd = ['xyz2grd -R',num2str(mer_xmin),'/',num2str(mer_xmax),'/',num2str(mer_ymin),'/',num2str(mer_ymax),' -I',num2str(mer_cols),'+/',num2str(mer_rows),'+ out.bin -Gtmp.grd  -F -ZTLf'];
    [a,b] = system(xyz2grd_cmd);

    % downsample and output without interpolation files
    grdsmp_cmd = ['grdsample -R',num2str(xmin),'/',num2str(xmax),'/',num2str(ymin),'/',num2str(ymax),' -I',num2str(smpres),' -F tmp.grd -Gtmp_smp.grd'];
    [a,b]=system(grdsmp_cmd);
    grd2xyz_cmd = ['grd2xyz -R',num2str(xmin),'/',num2str(xmax),'/',num2str(ymin),'/',num2str(ymax),' tmp_smp.grd -bo >' outfile_watervapor];
    [a,b]=system(grd2xyz_cmd);

    % load gaussian and nointerp back in, combine, and output as outfile_gauss
    % opening the data file (not-interpolated)
    nointfid = fopen(outfile_watervapor,'r');
    data_vector = fread(nointfid,'double');
    fclose(nointfid);
    % reshaping into the right n column matrix
    data = reshape(data_vector,3,[])';
    noint = data(:,3);
    xy = data(:,[1:2]);
    clear data data_vector

%     figure; scatter3(xy(:,1),xy(:,2),noint,15,noint,'filled'); view(0,90); axis equal; axis tight

    % writing out the date again as a binary table
    data_write = [xy noint]';
    clear noint xy
    %output
    fid = fopen(outfile_watervapor,'w');
    fwrite(fid,data_write,'double');
    fclose(fid);
    clear data_write

    fprintf([num2str(n) ' completed out of ' num2str(ndates) '\n'])

end

[a b] = system('!rm tmp.grd tmp.xyz out.bin tmp2.grd tmp_fil.grd tmp_smp.grd tmp_smp.xyz tmp_smp2.grd tmp_smp_fil.grd');

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function [ R ] = R_rpy( roll, pitch, yaw )
%R_RPY この関数の概要をここに記述
%   Camera rotates (1) yaw (y-axis) (2)pitch (x'-axis) (3)roll (z''-axis)

rrad = -roll * pi / 180.0;
prad = -pitch * pi / 180.0;
yrad = -yaw * pi / 180.0;

Rr = [cos(rrad), -sin(rrad), 0;...
      sin(rrad), cos(rrad), 0;...
      0, 0, 1];

Rp = [1, 0, 0;...
      0, cos(prad), -sin(prad);...
      0, sin(prad), cos(prad)];

Ry = [cos(yrad), 0, sin(yrad);...
      0, 1, 0;...
      -sin(yrad), 0, cos(yrad)];

R = Rr * Rp * Ry;

end

author: HajimeTaira


text:
function node_xyz = sphere01_triangle_project ( a_xyz, b_xyz, c_xyz, f1, ...
  f2, f3 )

%*****************************************************************************80
%
%% SPHERE01_TRIANGLE_PROJECT projects from plane to spherical triangle.
%
%  Discussion:
%
%    We assume that points A, B and C lie on the unit sphere, and they
%    thus define a spherical triangle.
%
%    They also, of course, define a planar triangle.
%
%    Let (F1,F2,F3) be the barycentric coordinates of a point in this
%    planar triangle.
%
%    This function determines the coordinates of the point in the planar
%    triangle identified by the barycentric coordinates, and returns the
%    coordinates of the projection of that point onto the unit sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A_XYZ(3), B_XYZ(3), C_XYZ(3), the coordinates
%    of the points A, B, and C.
%
%    Input, integer F1, F2, F3, the barycentric coordinates
%    of a point in the triangle ABC.  Normally, these coordinates would
%    be real numbers, and would sum to 1.  For convenience, we allow these
%    to be integers which must be divided by F1+F2+F3.
%
%    Output, real NODE_XYZ(3), the coordinates of the
%    point on the unit sphere which is the projection of the point on the plane
%    whose barycentric coordinates with respect to A, B, and C is
%    (F1,F2,F3)/(F1+F2+F3).
%

%
%  Destroy all row vectors!
%
  a_xyz = a_xyz(:);
  b_xyz = b_xyz(:);
  c_xyz = c_xyz(:);

  node_xyz(1:3,1) = ...
    ( ( f1           ) * a_xyz(1:3)   ...
    + (      f2      ) * b_xyz(1:3)   ...
    + (           f3 ) * c_xyz(1:3) ) ...
    / ( f1 + f2 + f3 );

  node_norm = r8vec_norm ( 3, node_xyz(1:3,1) );

  node_xyz(1:3,1) = node_xyz(1:3,1) / node_norm;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [germs, germPaths] = centroidalVoronoi2d(germs, poly, varargin)
%CENTROIDALVORONOI2D Centroidal Voronoi tesselation within a polygon
%
%   PTS = centroidalVoronoi2d(NPTS, POLY)
%   Generate points in a polygon based on centroidal voronoi tesselation.
%   Centroidal germs can be computed by using the Llyod's algorithm:
%   1) initial germs are chosen at random within polygon
%   2) voronoi polygon of the germs is computed
%   3) the centroids of each domain are computed, and used as germs of the
%   next iteration
%
%   [PTS, PATHLIST] = centroidalVoronoi2d(NPTS, POLY)
%   Also returns the path of each germs at each iteration. The result
%   PATHLIST is a cell array with as many cells as the number of germs,
%   containing in each cell the successive positions of the germ.
%
%   PTS = centroidalVoronoi2d(.., PARAM, VALUE)
%   Specify one or several optional arguments. PARAM can be one of:
%   * 'nIter'   specifies the number of iterations of the algorithm
%       (default is 50)
%   * 'verbose' display iteration number. Default is false.
%
%   Example
%     poly = ellipseToPolygon([50 50 40 30 20], 200);
%     nGerms = 100;
%     germs = centroidalVoronoi2d(nGerms, poly);
%     figure; hold on;
%     drawPolygon(poly, 'k');
%     drawPoint(germs, 'bo');
%     axis equal; axis([0 100 10 90]);
%     % extract regions of the CVD
%     box = polygonBounds(poly);
%     [n, e] = boundedVoronoi2d(box, germs);
%     [n2, e2] = clipGraphPolygon(n, e, poly);
%     drawGraphEdges(n2, e2, 'b');
%
%   See also
%   graphs, boundedVoronoi2d, centroidalVoronoi2d_MC
%
%   Rewritten from programs found in
%   http://people.scs.fsu.edu/~burkardt/m_src/cvt/cvt.html
%
%   Reference:
%    Qiang Du, Vance Faber, and Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review, Volume 41, 1999, pages 637-676.
%

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

%% Parse input arguments

% Number of germs
if isscalar(germs)
    nGerms = germs;
    germs = [];
else
    nGerms = size(germs, 1);
end

% Number of iterations
nIter = 50;

verbose = false;

keepPaths = nargout > 1;

while length(varargin) > 1
    paramName = varargin{1};
    switch lower(paramName)
        case 'verbose'
            verbose = varargin{2};
        case 'niter'
            nIter = varargin{2};

        otherwise
            error(['Unknown parameter name: ' paramName]);
    end

    varargin(1:2) = [];
end

%% Initialisations

% bounding box of polygon
bbox = polygonBounds(poly);

% init germs if needed
if isempty(germs)
    germs = generatePointsInPoly(nGerms);
end
germIters = cell(nIter, 1);

%% Iteration of the Lloyd algorithm

for i = 1:nIter
     if verbose
        disp(sprintf('Iteration: %d/%d', i, nIter)); %#ok<DSPS>
    end

    if keepPaths
        germIters{i} = germs;
    end

    % Compute Clipped Voronoi diagram of germs
    if verbose
        disp('  compute Voronoi Diagram');
    end
    [n, e, f] = boundedVoronoi2d(bbox, germs);
    [n2, e2, f2] = clipMesh2dPolygon(n, e, f, poly); %#ok<ASGLU>

    % update the position of each germ
    if verbose
        disp('  compute centroids');
    end
    for iGerm = 1:nGerms
        polygon = n2(f2{iGerm}, :);
        germs(iGerm,:) = polygonCentroid(polygon);
    end

end

%% Evenutally compute germs trajectories

if nargout > 1
    % init
    germPaths = cell(nGerms, 1);
    path = zeros(nIter+1, 2);

    % Iteration on germs
    for i = 1:nGerms

        % create path corresponding to germ
        for j = 1:nIter
            pts = germIters{j};
            path(j,:) = pts(i,:);
        end
        path(nIter+1, :) = germs(i,:);

        germPaths{i} = path;
    end
end

function pts = generatePointsInPoly(nPts)
    % extreme coordinates
    xmin = bbox(1);  xmax = bbox(2);
    ymin = bbox(3);  ymax = bbox(4);

    % compute size of box
    dx = xmax - xmin;
    dy = ymax - ymin;

    % allocate memory for result
    pts = zeros(nPts, 2);

    % iterate until all points have been sampled within the polygon
    ind = (1:nPts)';
    while ~isempty(ind)
        NI = length(ind);
        x = rand(NI, 1) * dx + xmin;
        y = rand(NI, 1) * dy + ymin;
        pts(ind, :) = [x y];

        ind = ind(~polygonContains(poly, pts(ind, :)));
    end
end

end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function[varargout]=col2mat(varargin)
%COL2MAT  Expands 'column-appended' data into a matrix.
%
%   M=COL2MAT(C), where C is a column vector of data segments separated by
%   NANs, returns a matrix M in which each segment has its own column.
%
%   M will have enough rows to fit the longest data segment. Segments are
%   read into columns of M from the top down, leaving empty spaces at the
%   bottom which are filled in with NANs.
%
%   If the data is complex, then gaps are filled with NAN+SQRT(-1)*NAN;
%
%   [M1,M2,...,MN]=COL2MAT(C1,C2,...,CN), where the CN are input column
%   vectors, also works, as does [M1,M2,...,MN]=COL2MAT(MAT), where MAT is
%   a matrix of column vectors.  In both cases the locations of NaNs in the
%   first column is used as the reference for the others.
%
%   COL2MAT(C1,C2,...); with no output arguments overwrites the original
%   input variables.
%
%   COL2MAT, MAT2COL, and COLBREAKS together form a system for moving data
%   with segments of nonuniform length rapidly back and forth between a
%   column format and a padded-matrix format.
%
%   Note that while COL2CELL and CELL2COL work for arrays having multiple
%   columns, COL2MAT only works with column vectors.
%
%   See also MAT2COL, COLBREAKS, COL2CELL, CELL2COL.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2016 J.M. Lilly --- type 'help jlab_license' for details

%no loops!

if ischar(varargin{1})
    if strcmpi(varargin{1}(1:3),'--t')
        return
    end
end

if isempty(varargin{1})
   for i=1:nargout
       varargout{i}=[];
   end
   return
end

if nargin>1
	for i=1:nargin
        	eval(['data(:,' int2str(i) ')=varargin{' int2str(i) '};'])
	end
else
	data=varargin{1};
end

for ii=1:size(data,2)
	col=data(:,ii);
	if ii==1;
        %Account for potential missing NaNs at the end
        if ~isnan(col(end))
            col(end+1)=nan;
            data(end+1,:)=nan;
        end
		%use index into NANs and index derivative
		%to determine size of new matix
		nani=find(isnan(col));
		dnani=diff([0;nani]);
		nrows=max(dnani);
		ncols=length(nani);

		%determine index into top of each column ==a
		%and index into first NAN in each column ==b
		a=1+nrows*(0:ncols-1)';
        %size(dnani)
        %size(a)
		b=dnani+a;
		index=zeros(nrows*ncols,1);

		b=b(b<length(index));
		%mark all the numbers between each a and each b
		index(a)=1;
		index(b)=index(b)-1;%this matters (a may equal b)
		index=find(cumsum(index));
		if length(index)>length(col)
		   index=index(1:length(col));
		end
	end
	mat=nan*ones(nrows,ncols);
    mat(index)=col;
    if ~isreal(mat)
        vswap(mat,nan,nan+sqrt(-1)*nan);
    end
	varargout{ii}=mat;
end

if nargout>nargin
	varargout{int2str(size(data,2)+1)}=index;
end

if nargout==0
  eval(to_overwrite(nargin));
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [beta_gibbs,sigma_gibbs]=ndgibbs(It,Bu,beta0,omega0,X,Y,y,Bhat,n,T,q)

% function [beta_gibbs sigma_gibbs]=ndgibbs(It,Bu,beta0,omega0,X,Y,y,Bhat,n,T,q)
% performs the Gibbs algorithm 1.5.2 for the normal-diffuse prior, and returns draws from posterior distribution
% inputs:  - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - vector 'beta0': vector of prior values for beta (defined in 1.3.4)
%          - matrix 'omega0': prior covariance matrix for the VAR coefficients (defined in 1.3.8)
%          - matrix 'X': matrix of regressors for the VAR model (defined in 1.1.8)
%          - matrix 'Y': matrix of regressands for the VAR model (defined in 1.1.8)
%          - vector 'y': vectorised regressands for the VAR model (defined in 1.1.12)
%          - matrix 'Bhat': OLS VAR coefficients, in non vectorised form (defined in 1.1.9)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'T': number of sample time periods (defined p 7 of technical guide)
%          - integer 'q': total number of coefficients to estimate for the BVAR model (defined p 7 of technical guide)
% outputs: - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - matrix'sigma_gibbs': record of the gibbs sampler draws for the sigma matrix (vectorised)

% preliminary tasks

% invert omega0, as it will be used repeatedly during step 4
invomega0=diag(1./diag(omega0));

% set initial values for B (step 2); use OLS estimates
B=Bhat;

% create the progress bar
hbar = bear.parfor_progressbar(It,'Progress of the Gibbs sampler.');

% start iterations
for ii=1:It

% Step 3: at iteration ii, first draw sigma from IW, conditional on beta from previous iteration
% obtain first Shat, defined in (1.6.10)
Shat=(Y-X*B)'*(Y-X*B);
% Correct potential asymmetries due to rounding errors from Matlab
C=chol(bear.nspd(Shat));
Shat=C'*C;

% next draw from IW(Shat,T)
sigma=bear.iwdraw(Shat,T);

% step 4: with sigma drawn, continue iteration ii by drawing beta from a multivariate Normal, conditional on sigma obtained in current iteration
% first invert sigma
C=chol(bear.nspd(sigma));
invC=C\speye(n);
invsigma=invC*invC';

% then obtain the omegabar matrix
invomegabar=invomega0+kron(invsigma,X'*X);
C=chol(bear.nspd(invomegabar));
invC=C\speye(q);
omegabar=invC*invC';

% following, obtain betabar
betabar=omegabar*(invomega0*beta0+kron(invsigma,X')*y);

% draw from N(betabar,omegabar);
beta=betabar+chol(bear.nspd(omegabar),'lower')*mvnrnd(zeros(q,1),eye(q))';

% update matrix B with each draw
B=reshape(beta,size(B));

% update progress by one iteration
hbar.iterate(1);

% record the values if the number of burn-in iterations is exceeded
if ii>Bu
% values of vector beta
beta_gibbs(:,ii-Bu)=beta;
% values of sigma (in vectorized form)
sigma_gibbs(:,ii-Bu)=sigma(:);
% if current iteration is still a burn iteration, do not record the result
else
end

% go for next iteration
end
% close progress bar
close(hbar);

author: european-central-bank

repo: BEAR-toolbox


text:
function new_model = rotate_3D_model_around_z(model,phi)

r = [cosd(phi) -sind(phi) 0 ; sind(phi) cosd(phi) 0 ; 0 0 1];
new_model = transform_3d_rigid(model,r,[0;0;0]);

end

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function test_suite = test_intersectLinePolygon
%TESTINTERSECTLINEPOLYGON Test case for function intersectLinePolygon
%   output = testIntersectLinePolygon(input)
%
%   Example
%   testIntersectLinePolygon
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testSquare(testCase) %#ok<*DEFNU>
% test with a square and orthogonal lines, inside and outside

poly = [0 0;10 0;10 10;0 10];

lineH1 = [5 5 1 0];
testCase.assertEqual(2, size(intersectLinePolygon(lineH1, poly), 1));

lineH2 = [5 15 1 0];
testCase.assertEqual(0, size(intersectLinePolygon(lineH2, poly), 1));

lineV1 = [5 5 0 1];
testCase.assertEqual(2, size(intersectLinePolygon(lineV1, poly), 1));

lineV2 = [15 5 0 1];
testCase.assertEqual(0, size(intersectLinePolygon(lineV2, poly), 1));

function testClosedSquare(testCase)
% test when the polygon has same vertices at end and at bginning

poly = [0 0;10 0;10 10;0 10;0 0];

lineH1 = [5 5 1 0];
testCase.assertEqual(2, size(intersectLinePolygon(lineH1, poly), 1));

lineH2 = [5 15 1 0];
testCase.assertEqual(0, size(intersectLinePolygon(lineH2, poly), 1));

lineV1 = [5 5 0 1];
testCase.assertEqual(2, size(intersectLinePolygon(lineV1, poly), 1));

lineV2 = [15 5 0 1];
testCase.assertEqual(0, size(intersectLinePolygon(lineV2, poly), 1));

function testDiamond(testCase)

poly = [10 0;20 10;10 20;0 10];

% definition from 14/02/2022 do not intersect lower vertices
lineH0 = [0 0 3 0];
intersects = intersectLinePolygon(lineH0, poly);
testCase.assertTrue(isempty(intersects));

lineH1 = [10 10 3 0];
intersects = intersectLinePolygon(lineH1, poly);
target = [20 10];
testCase.assertTrue(ismember(target, intersects, 'rows'));
target = [0 10];
testCase.assertTrue(ismember(target, intersects, 'rows'));

[intersects, inds] = intersectLinePolygon(lineH1, poly);
testCase.assertEqual(size(intersects, 1), size(inds, 1));

lineH2 = [0 20 3 0];
intersects = intersectLinePolygon(lineH2, poly);
target = [10 20];
testCase.assertTrue(ismember(target, intersects, 'rows'));

[intersects, inds] = intersectLinePolygon(lineH2, poly);
testCase.assertEqual(size(intersects, 1), size(inds, 1));

% definition from 14/02/2022 intersects upper vertices twice
lineV0 = [0 0 0 3];
intersects = intersectLinePolygon(lineV0, poly);
testCase.assertEqual(size(intersects, 2), 2);

[intersects, inds] = intersectLinePolygon(lineV0, poly);
testCase.assertEqual(size(intersects, 1), size(inds, 1));

lineV1 = [10 10 0 3];
intersects = intersectLinePolygon(lineV1, poly);
testCase.assertEqual(size(intersects, 1), 2);

[intersects, inds] = intersectLinePolygon(lineV1, poly);
testCase.assertEqual(size(intersects, 1), size(inds, 1));

% definition from 14/02/2022 do not intersect lower vertices
lineV2 = [20 0 0 3];
intersects = intersectLinePolygon(lineV2, poly);
testCase.assertTrue(isempty(intersects));

[intersects, inds] = intersectLinePolygon(lineV2, poly);
testCase.assertEqual(size(intersects, 1), size(inds, 1));

function testMShape(testCase)
% a more complicated polygon, with 4 intersections

poly = [10 10;60 10;60 40;40 20;30 20;10 40];
line = [0 30 3 0];

[inters, inds] = intersectLinePolygon(line, poly);
expInters = [10 30;20 30;50 30;60 30];
for i = 1:4
    assertTrue(testCase, ismember(expInters(i,:), inters, 'rows'));
end
expInds = [6;5;3;2];
testCase.assertEqual(ismember(expInds, inds), true(4,1));

function testUniquePoints(testCase)
% Check that function returns unique results, even for vertex points

poly = [0 0;10 0;10 10;0 10];
line = [5 5 1 1];
intersects = intersectLinePolygon(line, poly);
testCase.assertEqual(2, size(intersects, 1), 'Wrong number of intersections');

function testGetEdgesIndices(testCase)

poly = [0 0;10 0;10 10;0 10];
line = [5 5 1 0];
[intersects, inds] = intersectLinePolygon(line, poly); %#ok<ASGLU>

testCase.assertEqual(2, length(inds));
testCase.assertEqual([4;2], inds);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function p = bst_cross(x,y,dim)
% BST_CROSS: Cross product between two sets, each set with three columns
%
% USAGE:  p = bst_cross(x,y)
%         p = bst_cross(x,y,dim)
%
% INPUT:
%     - x   : [1,3] double or single
%     - y   : [1,3] double or single
%     - dim : dimension along which to compute the cross product
%
% NOTE:
%     - Does exactly then same as the Matlab 'cross' function, but much faster.

% @=============================================================================
% 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.
% =============================================================================@

% Default: dim=1
if (nargin < 3) || isempty(dim)
    dim = 1;
end
% Check size
if (size(x,dim)~=3) || (size(y,dim)~=3)
    error(' Must have three columns ');
end
% Compute cross product
if (dim == 2)
    p = [x(:,2).*y(:,3) - x(:,3).*y(:,2),...
         x(:,3).*y(:,1) - x(:,1).*y(:,3),...
         x(:,1).*y(:,2) - x(:,2).*y(:,1)];
else
    p = [x(2,:).*y(3,:) - x(3,:).*y(2,:);...
         x(3,:).*y(1,:) - x(1,:).*y(3,:);...
         x(1,:).*y(2,:) - x(2,:).*y(1,:)];
end

author: brainstorm-tools


text:
% In a paper you'll often see a p-value and an N, based on either a T-test,
% a paired T-test or a one-way ANOVA.  To get an intuition how "good" a
% p-value is given the N, this script simulates some of those tests.

%% General parameter settings
nrSimulatedExperiments = 1000;
sampleSizes = [10 20 50];
testToSimulate = '2WAYRMMAIN';   % Pick one from 'TTEST','1WAY',2WAYRMMAIN 2WAYRMINT
nrSampleSizes= numel(sampleSizes);
scale =1;
bayesFactor = nan(nrSampleSizes,nrSimulatedExperiments);
bayesFactorBIC = nan(nrSampleSizes,nrSimulatedExperiments);
pValue = nan(nrSampleSizes,nrSimulatedExperiments);
stat =  nan(nrSampleSizes,nrSimulatedExperiments);
effectSize = 0.25;
tic;nn=0;

if ~exist('showTimeToCompletion.m','file')
    % Showtimetocompletion is a function outside the toolbox -do nothing
    showTimeToCompletion = @(x,y)(1); %NOP
end
switch upper(testToSimulate)
    case 'TTEST'
        %% Two-Sample T Test
        for i=1:nrSimulatedExperiments
            for j = 1:nrSampleSizes
                X = randn([sampleSizes(j) 1]);
                Y = randn([sampleSizes(j) 1])+effectSize;
                [bayesFactor(j,i),pValue(j,i),~,stats] = bf.ttest2(X,Y,'scale',scale);
                [bayesFactorBIC(j,i)] = bf.bfFromT(stats.tstat,stats.df); % Calculate BIC approximation
                stat(j,i) = stats.tstat;
                nn = showTimeToCompletion(((i-1)*nrSampleSizes+j)/(nrSimulatedExperiments*nrSampleSizes),nn);
            end
        end
    case '1WAY'
        %% 1-WAY ANOVA - main effect
        nrLevels = 4;
        cntr=0;
        for i=1:nrSimulatedExperiments
            for j = 1:nrSampleSizes
                X= repmat((0:nrLevels-1),[sampleSizes(j) 1]);
                X=X(:);
                y = randn([sampleSizes(j)*nrLevels 1])+X*effectSize;
                subject = repmat(1:sampleSizes(j),[nrLevels 1])';
                tbl = table(X,y,subject(:));
                [bayesFactor(j,i),lme] = bf.anova(tbl,'y~X','scale',scale);

                tmp =anova(lme);
                pValue(j,i) =tmp.pValue(2); % pValue of the Main factor
                stat(j,i) = tmp.FStat(2);
                df1 = nrLevels-1;
                df2 = nrLevels*(sampleSizes(j)-1);
                [bayesFactorBIC(j,i)] = bf.bfFromF(tmp.FStat(2),df1,df2,sampleSizes(j));
                nn = showTimeToCompletion(((i-1)*nrSampleSizes+j)/(nrSimulatedExperiments*nrSampleSizes),nn);                        end
        end
    case {'2WAYRMMAIN','2WAYRMINT'}
        %%
        levels = [2 3];
        nrFactors = numel(levels);
        for i=1:nrSimulatedExperiments
            for j = 1:nrSampleSizes
                fac1 = repmat((1:levels(1))',[1 levels(2)]);
                fac2 = repmat(1:levels(2),[levels(1) 1 ]);
                first = effectSize*randn([levels(1) 1]);
                second = effectSize*randn([1 levels(2)]);
                interaction = effectSize.*fac1.*fac2;
                singleSubject = repmat(first,[1 levels(2)]) + repmat(second,[levels(1) 1]) + interaction;
                response = reshape(repmat(singleSubject,[1 1 sampleSizes(j)])+ randn([levels sampleSizes(j)]),[prod(levels)*sampleSizes(j) 1]);
                fac1 = reshape(repmat(fac1,[1 1 sampleSizes(j)]),[prod(levels)*sampleSizes(j) 1]);
                fac2 = reshape(repmat(fac2,[1 1 sampleSizes(j)]),[prod(levels)*sampleSizes(j) 1]);
                subject = reshape(repmat((1:sampleSizes(j)),[prod(levels) 1 ]),[prod(levels)*sampleSizes(j) 1]);
                tbl = table(response,fac1,fac2,subject);
                switch upper(testToSimulate)
                    case '2WAYRMMAIN'
                        % Main effect Factor 1
                        full = bf.anova(tbl,'response~fac1+subject','scale',scale,'treatAsRandom',{'subject'});
                        restricted = bf.anova(tbl,'response~subject','scale',scale,'treatAsRandom',{'subject'});
                        bayesFactor(j,i) = full/restricted;
                        lm =fitlme(tbl,'response~fac1+(1|subject)');
                        tmp =anova(lm);
                        pValue(j,i) =tmp.pValue(2); % pValue of the Main factor
                        stat(j,i) = tmp.FStat(2);
                        df1 = levels(1)-1;
                        df2 = df1*(sampleSizes(j)-1);
                        bayesFactorBIC(j,i) = bf.bfFromF( stat(j,i),df1,df2,sampleSizes(j));

                    case '2WAYRMINT'
                        % Interaction effect
                        full = bf.anova(tbl,'response~fac1*fac2+subject','scale',scale,'treatAsRandom',{'subject'});
                        restricted = bf.anova(tbl,'response~fac1+fac2+subject','scale',scale,'treatAsRandom',{'subject'});
                        bayesFactor(j,i) = full/restricted;
                        lm =fitlme(tbl,'response~fac1*fac2+(1|subject)');
                        tmp =anova(lm);
                        pValue(j,i) =tmp.pValue(end); % pValue of the Main factor
                        stat(j,i) = tmp.FStat(end);
                        df1 = prod(levels-1);
                        % df error for repeated measures: total-within-subjects
                        df2 = df1*(sampleSizes(j)-1);
                        bayesFactorBIC(j,i) = bf.bfFromF(stat(j,i),df1,df2,sampleSizes(j));
                end
                nn = showTimeToCompletion(((i-1)*nrSampleSizes+j)/(nrSimulatedExperiments*nrSampleSizes),nn);
            end
            %             consistency =100*mean(sign(log10(bayesFactor(j,:)))==sign(log10(bayesFactorBIC(j,:))));
            %             sprintf('%.2f %.2f %.2f %.2f %.2f %%',effectSize(1),min(log10(bayesFactor(:,j))),nanmedian(log10(bayesFactor(:,j))),max(log10(bayesFactor(:,j))),consistency)
            %             sprintf('%.2f %.2f %.2f %.2f %.2f %%',effectSize(1),min(log10(bayesFactorBIC(:,j))),nanmedian(log10(bayesFactorBIC(:,j))),max(log10(bayesFactorBIC(:,j))),consistency)
        end

end

%% Generate a figure
% 1 - compare pValue with BF
% 2 - compare stat (F or T) with BF
% 3 -  compare BIC approximation with BF

figure;clf
minP = 1e-4;
maxP = 0.1;
h = nan(nrSampleSizes,1);
xtick = sort([10.^(-3:1:-1) 0.05]);
for j = 1:nrSampleSizes
    subplot(2,2,1);
    x = pValue(j,:,:);
    [x,ix]=sort(x(:));
    y = bayesFactor(j,:,:);
    y=y(ix);
    h(j)= plot(x,y(:),'.-');
    set(gca,'YScale','Log','YLim',[0 1000],'XScale','log','XLim',[minP maxP],'XTIck',xtick);
    hold on
    subplot(2,2,2);
    x = stat(j,:,:);
    [x,ix]=sort(x(:));
    y = bayesFactor(j,:,:);
    y=y(ix);
    h(j)= plot(x,y(:),'.-');
    set(gca,'YScale','Log','YLim',[0 1000]);%,'XScale','log','XLim',[minP maxP],'XTIck',xtick);
    hold on
end
subplot(2,2,1)
plot(xlim,[1 1],'k--'); text(maxP,1,'Equal','Color','k')
plot(xlim,[3 3],'r--'); text(maxP,3,'Barely','Color','r')
plot(xlim,[6 6],'g--');text(maxP,6,'Moderate','Color','g')
plot(xlim,[10 10],'b--');text(maxP,10,'Strong','Color','b')
ylabel 'Bayes Factor'
xlabel 'p-Value'
title (['Simulating: ' testToSimulate])
bfLim =ylim;  % Take the bf limits that go with the p-limits we imposed
legend(h,num2str(sampleSizes'))
subplot(2,2,2);
plot(xlim,[1 1],'k--'); text(maxP,1,'Equal','Color','k')
plot(xlim,[3 3],'r--'); text(maxP,3,'Barely','Color','r')
plot(xlim,[6 6],'g--');text(maxP,6,'Moderate','Color','g')
plot(xlim,[10 10],'b--');text(maxP,10,'Strong','Color','b')
ylabel 'Bayes Factor'
xlabel 'Test Statistic (T/F)'
legend(h,num2str(sampleSizes'))

subplot(2,2,3)
for j = 1:nrSampleSizes
    x = bayesFactor(j,:,:);
    y = bayesFactorBIC(j,:,:);
    plot(x(:),y(:),'.');
    hold on
end
axis square
axis equal
ylabel 'BIC approximation'
xlabel 'Bayes Factor'
set(gca,'XScale','log','YScale','Log','XLim',bfLim,'YLim',bfLim);
plot(bfLim,bfLim,'k--')
title 'Comparing BF with BIC Approximation'

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function odf = calcBinghamODF(ori,varargin)
% calculate ODF from individuel orientations via kernel density estimation
%
% *calcODF* is one of the core function of the MTEX toolbox.
% It estimates an ODF from a set of individual crystal orientations by
% <EBSD2ODF.html kernel density estimation>.
%
% Input
%  ori  - @orientation
%
% Output
%  odf - @SO3Fun
%
% See also
% ebsd_demo EBSD2odf EBSDSimulation_demo EBSD/load EBSD/calcKernel kernel/kernel

% maybe there is nothing to do
if isempty(ori), odf = ODF; return, end

% estimate Bingham parameters
[~,~,lambda,ev] = mean(ori,varargin{:});
kappa = evalkappa(lambda,varargin{:});

% set up Bingham ODF
odf = BinghamODF(kappa,ev,ori.CS,ori.SS);

author: mtex-toolbox

repo: mtex

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


text:
function [d pred] = shortest_paths(A,u,varargin)
% SHORTEST_PATHS Compute the weighted single source shortest path problem.
%
% [d pred] = shortest_paths(A,u) returns the distance (d) and the predecessor
% (pred) for each of the vertices along the shortest path from u to every
% other vertex in the graph.
%
% ... = shortest_paths(A,u,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   options.algname: the algorithm to use
%       [{'auto'} | 'dijkstra' | 'bellman_ford' | 'dag']
%   options.inf: the value to use for unreachable vertices
%       [double > 0 | {Inf}]
%   options.target: a special vertex that will stop the search when hit
%       [{'none'} | any vertex number besides the u]; target is ignored if
%       visitor is set.
%   options.visitor: a structure with visitor callbacks.  This option only
%       applies to dijkstra or bellman_ford algorithms.  See dijkstra_sp or
%       bellman_ford_sp for details on the visitors.
%   options.edge_weight: a double array over the edges with an edge
%       weight for each edge, see EDGE_INDEX and EXAMPLES/REWEIGHTED_GRAPHS
%       for information on how to use this option correctly
%       [{'matrix'} | length(nnz(A)) double vector]
%
% Note: if you need to compute shortest paths with 0 weight edges, you must
% use an edge_weight vector, see the examples for details.
%
% Note: 'auto' cannot be used with 'nocheck' = 1.  The 'auto' algorithm
% checks if the graph has negative edges and uses bellman_ford in that
% case, otherwise, it uses 'dijkstra'.  In the future, it may check if the
% graph is a dag and use 'dag'.
%
% Example:
%    load graphs/clr-25-2.mat
%    shortest_paths(A,1)
%    shortest_paths(A,1,struct('algname','bellman_ford'))
%
% See also DIJKSTRA_SP, BELLMAN_FORD_SP, DAG_SP

% David Gleich
% Copyright, Stanford University, 2006-2008

%% History
%  2006-04-19: Initial coding
%  2007-04-18: Added edge_weight option.
%  2007-04-19: Added target option.
%    Added additional error checks.
%  2007-07-12: Fixed edge_weight documentation
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end

options = struct('algname', 'auto', 'inf', Inf, 'edge_weight', 'matrix', ...
    'target', 'none');
options = merge_options(options,varargin{:});

% edge_weights is an indicator that is 1 if we are using edge_weights
% passed on the command line or 0 if we are using the matrix.
edge_weights = 0;
edge_weight_opt = 'matrix';

if strcmp(options.edge_weight, 'matrix')
    % do nothing if we are using the matrix weights
else
    edge_weights = 1;
    edge_weight_opt = options.edge_weight;
end

if strcmp(options.target,'none')
    target = 0; % a flag used to denote "no target" to the mex
elseif isa(options.target, 'double')
    target = options.target;
else
    error('matlab_bgl:invalidParameter', ...
        'options.target is not ''none'' or a vertex number.');
end

if check
    % check the values of the matrix
    check_matlab_bgl(A,struct('values',edge_weights ~= 1));

    if edge_weights && nnz(A) ~= length(edge_weight_opt)
        error('matlab_bgl:invalidParameter', 'the vector of edge weights must have length nnz(A)');
    end

    % set the algname
    if (strcmpi(options.algname, 'auto'))
        if edge_weights
            mv = min(edge_weights);
        else
            mv = min(min(A));
        end

        if (mv < 0)
            options.algname = 'bellman_ford';
        else
            options.algname = 'dijkstra';
        end
    else
        % check the data provided to match the algorithm
        if strcmpi(options.algname, 'dijkstra')
            if edge_weights
                mv = min(edge_weight_opt);
            else
                mv = min(min(A));
            end
            if mv < 0
                error('matlab_bgl:invalidParameter', ...
                    'dijkstra''s algorithm cannot be used with negative edge weights.');
            end
        end
    end

else
    if (strcmpi(options.algname, 'auto'))
        error('shortest_paths:invalidParameter', ...
            'algname auto is not compatible with no check');
    end
end

if options.inf < 0, error('options.inf must be larger than 0'); end

if trans, A = A'; end

if isfield(options,'visitor')
    [d pred] = matlab_bgl_sp_mex(A,u,target,lower(options.algname),options.inf,...
        edge_weight_opt, options.visitor);
else
    [d pred] = matlab_bgl_sp_mex(A,u,target,lower(options.algname),options.inf,...
        edge_weight_opt);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
clear, clc;

% This is an example for running the function mtLeastR
%
%  Problem:
%
%  min  1/2 || A x - y||^2 + rho * sum_j ||x^j||_q
%
%  x is grouped into k groups according to opts.ind
%  The indices of x_j in x is (ind(j)+1):ind(j+1)
%
% For detailed description of the function, please refer to the Manual.
%
%% ------------   History --------------------
% First version on August 10, 2009.
%
% September 5, 2009: adaptive line search is added
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/L1Lq;

m=1000;  n=100;     % the size of the data matrix
k=10;               % 10 tasks
ind=0:100:1000;     % the 1000 samples are from 10 tasks
randNum=1;          % a random number
q=2;                % the value of q in the L1/Lq regularization
rho=0.5;            % the regularization parameter

% ---------------------- generate random data ----------------------
randn('state',(randNum-1)*3+1);
A=randn(m,n);        % the data matrix

randn('state',(randNum-1)*3+2);
xOrin=randn(n,k);

randn('state',(randNum-1)*3+3);
noise=randn(m,1);

for i=1:k
    ind_i=(ind(i)+1):ind(i+1);
    y(ind_i,1)=A(ind_i,:)*xOrin(:,i)+...
        noise(ind_i,1)*0.01;
                     % the response
end

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% Termination
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=100;   % maximum number of iterations

% Normalization
opts.nFlag=0;       % without normalization

% Regularization
opts.rFlag=1;       % the input parameter 'rho' is a ratio in (0, 1)

% Group Property
opts.q=q;           % set the value for q
opts.ind=ind;       % set the group indices

%----------------------- Run the code mtLeastR -----------------------
fprintf('\n mFlag=0, lFlag=0 \n');
opts.mFlag=0;       % treating it as compositive function
opts.lFlag=0;       % Nemirovski's line search
tic;
[x1, funVal1, ValueL1]= mtLeastR(A, y, rho, opts);
toc;

opts.maxIter=200;

fprintf('\n mFlag=1, lFlag=0 \n');
opts.mFlag=1;       % smooth reformulation
opts.lFlag=0;       % Nemirovski's line search
opts.tFlag=2; opts.tol= funVal1(end);
tic;
[x2, funVal2, ValueL2]= mtLeastR(A, y, rho, opts);
toc;

fprintf('\n mFlag=1, lFlag=1 \n');
opts.mFlag=1;       % smooth reformulation
opts.lFlag=1;       % adaptive line search
opts.tFlag=2; opts.tol= funVal1(end);
tic;
[x3, funVal3, ValueL3]= mtLeastR(A, y, rho, opts);
toc;

figure;
plot(funVal1,'-r');
hold on;
plot(funVal2,'--b');
hold on;
plot(funVal3,':g');
legend('mFlag=0, lFlag=0', 'mFlag=1, lFlag=0', 'mFlag=1, lFlag=1');
xlabel('Iteration (i)');
ylabel('The objective function value');

% % --------------------- compute the pathwise solutions ----------------
% opts.fName='mtLeastR';    % set the function name to 'mtLeastR'
% Z=[0.9, 0.8, 0.5, 0.3];   % set the parameters
%
% % run the function pathSolutionLeast
% fprintf('\n Compute the pathwise solutions, please wait...');
% X=pathSolutionLeast(A, y, Z, opts);

author: zzstefan

repo: BrainNetClass


text:
function varargout = LoginGUI(varargin)
% LOGINGUI MATLAB code for LoginGUI.fig
%      LOGINGUI, by itself, creates a new LOGINGUI or raises the existing
%      singleton*.
%
%      H = LOGINGUI returns the handle to a new LOGINGUI or the handle to
%      the existing singleton*.
%
%      LOGINGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LOGINGUI.M with the given input arguments.
%
%      LOGINGUI('Property','Value',...) creates a new LOGINGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before LoginGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to LoginGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help LoginGUI

% Last Modified by GUIDE v2.5 30-May-2011 10:59:28

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @LoginGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @LoginGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [], ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% loginGUI(fcnHandle)
% fcnHandle has signature fcnHandle(username) and is executed when the
% Login button is pressed. If the window is killed, fcnHandle is not
% called.
function LoginGUI_OpeningFcn(hObject, eventdata, handles, varargin)
assert(numel(varargin)==1,'Expected exactly one input argument.');
handles.output = hObject;
handles.loginFcn = varargin{1};
guidata(hObject, handles);
centerfig(hObject);
set(handles.figure1,'WindowStyle','modal')

function varargout = LoginGUI_OutputFcn(hObject, eventdata, handles) %#ok<*INUSL>

% putting these two lines here rather than in OpeningFcn gets the GUI to
% open with the focus in etUsername
uicontrol(handles.etUsername);
uiwait(handles.figure1);

varargout{1} = [];

function etUsername_Callback(hObject, eventdata, handles) %#ok<*INUSD,*DEFNU>
% If user enters a username and hits return, move focus to pbLogin
c = get(handles.figure1,'CurrentCharacter');
if strcmp(c,sprintf('\r'))
    uicontrol(handles.pbLogin);
end

function edit2_Callback(hObject, eventdata, handles)

function pbLogin_Callback(hObject, eventdata, handles)
username = get(handles.etUsername,'String');
username = strtrim(username);
if isempty(username)
    errordlg('Please enter a username.','Missing Username','modal');
    return;
end
delete(handles.figure1);
handles.loginFcn(username);

function figure1_KeyPressFcn(hObject, eventdata, handles)
switch eventdata.Key
    case 'escape'
        delete(handles.figure1);
end

function pbLogin_KeyPressFcn(hObject, eventdata, handles)
switch eventdata.Key
    case 'return'
        pbLogin_Callback(hObject,[],handles);
    case 'escape'
        % this is not working
        set(handles.figure1,'CurrentObject',handles.figure1);
end

author: kristinbranson

repo: JAABA


text:
function tests = test_bdhizsqr
  tests = functiontests(localfunctions);
end

function verify_bd_full_svd(a, b, testCase)
    n = numel(a);
    A = full(spdiags([a [0; b]], [0 1], n, n));
    [UU, SS, VV] = svd(A);
    options.verbosity = 0;
    [U, S, V] = spx.fast.bdhizsqr_svd(a,b, options);
    verifyEqual(testCase, full(diag(S)), diag(SS), 'AbsTol', 1e-12);
    verifyEqual(testCase, spx.la.nonorthogonality(U), 0, 'AbsTol', 1e-12);
    verifyEqual(testCase, spx.la.nonorthogonality(V), 0, 'AbsTol', 1e-12);
    AA = U * S * V';
    verifyEqual(testCase, A, AA, 'AbsTol', 1e-12);
    verifyEqual(testCase, abs(U), abs(UU), 'AbsTol', 1e-12);
end

function verify_only_singular_values(a, b, testCase)
    n = numel(a);
    A = full(spdiags([a [0; b]], [0 1], n, n));
    SS = svd(A);
    S = spx.fast.bdhizsqr_svd(a,b);
    verifyEqual(testCase, S, SS, 'AbsTol', 1e-12);
end

function verify_u_last_row(a, b, testCase)
    n = numel(a);
    A = full(spdiags([a [0; b]], [0 1], n, n));
    [UU, SS, VV] = svd(A);
    options.u_rows = [n];
    [U, S] = spx.fast.bdhizsqr_svd(a,b, options);
    verifyEqual(testCase, S, diag(SS), 'AbsTol', 1e-12);
    U2 = UU(end, :);
    verifyEqual(testCase, abs(U), abs(U2), 'AbsTol', 1e-12);
end

function test_last_row(testCase)
    a = [1/2 1/4 1/8 1/16 1/32 1/64]'; b = [1 2 -1 -2 1]';
    verify_u_last_row(a, b, testCase);
    a = [1 2 3 3 2 1]'; b = [1 2 -1 -2 1]';
    verify_u_last_row(a, b, testCase);
    a = [1 2 3 4 5 6]'; b = [1 2 -1 -2 1]';
    verify_u_last_row(a, b, testCase);
    a = [1/2 1/4 1/8 1/16 1/32 1/64]'; b = [1 2 -1 -2 1]';
    verify_u_last_row(a, b, testCase);
end

function test_only_singular_values(testCase)
    a = [1 2 3 3 2 1]'; b = [1 2 -1 -2 1]';
    verify_only_singular_values(a, b, testCase);
    a = [1 2 3 4 5 6]'; b = [1 2 -1 -2 1]';
    verify_only_singular_values(a, b, testCase);
    a = [1/2 1/4 1/8 1/16 1/32 1/64]'; b = [1 2 -1 -2 1]';
    verify_only_singular_values(a, b, testCase);
end

function test_1(testCase)
    a = [1 2 3 3 2 1]'; b = [1 2 -1 -2 1]';
    verify_bd_full_svd(a, b, testCase);
end

function test_2(testCase)
    a = [1 2 3 4 5 6]'; b = [1 2 -1 -2 1]';
    verify_bd_full_svd(a, b, testCase);
end

function test_3(testCase)
    a = [1/2 1/4 1/8 1/16 1/32 1/64]'; b = [1 2 -1 -2 1]';
    verify_bd_full_svd(a, b, testCase);
end

function test_4(testCase)
    a = [1/64 1/32 1/16 1/8 1/4 1/2]'; b = [1 2 -1 -2 1]';
    verify_bd_full_svd(a, b, testCase);
end

function test_5(testCase)
    a = ones(6, 1); b = zeros(5, 1);
    verify_bd_full_svd(a, b, testCase);
end
function test_6(testCase)
    a = [1 1 1/100]'; b = [0  1]';
    verify_bd_full_svd(a, b, testCase);
end
function test_7(testCase)
    a = [-1 1 1/100]'; b = [0  -1]';
    verify_bd_full_svd(a, b, testCase);
end
function test_8(testCase)
    a = [-1 1/10 1/100]'; b = [.1 -.1]';
    verify_bd_full_svd(a, b, testCase);
end
function test_9(testCase)
    a = [1/100 -1/10 1]'; b = [-.1 .1]';
    verify_bd_full_svd(a, b, testCase);
end

function test_10(testCase)
    alpha = [2.0074 1.3912 0.9620 0.9705 1.1207 1.1491 1.0648 1.0474 0.9001 1.0274 1.0713 1.0241 1.0928 0.9850 1.0605 1.1760 0.9556 1.0401 1.2332 1.0288 0.9445 1.0836 1.0105 0.9715 1.1755 1.1439 0.9918 1.0790 1.0017 0.9292 1.0069 1.1838 0.9330 0.8743 0.9705 0.9932 0.9413 0.8130 0.8581 1.0076 0.9950 0.8767 0.9045 0.7637 0.8717 0.8697 0.8927 0.9939 0.9285 0.8069 0.8440 0.8879 0.8060 0.8034 0.8054 0.8393 0.8068 0.8077 0.8563 0.7771 0.8622 0.8232 0.8499 0.8687 0.8052 0.8193 0.0000]';
    beta =  [0.3983 1.0517 0.9578 1.1050 0.9405 0.9720 1.1747 1.0656 0.8632 1.1029 1.0015 0.9553 1.0921 1.0258 0.9547 1.0230 0.9264 1.0586 1.0231 0.8917 0.9671 0.9876 0.9214 0.9514 1.0351 1.0116 0.9881 0.9670 0.8618 0.9276 1.1199 1.0183 0.9144 0.8964 0.9943 0.9560 0.8559 0.9272 0.8591 1.0369 0.9838 0.8429 0.8651 0.8713 0.8399 0.9075 1.0259 0.9556 0.8768 0.8669 0.7960 0.8196 0.8519 0.8619 0.8613 0.7670 0.8661 0.8501 0.7928 0.8507 0.8567 0.8357 0.7723 0.8556 0.8403 0.8092]';
    n = numel(alpha);
    for i=2:6:n
        a = alpha(1:i);
        b = beta(1:i-1);
        verify_bd_full_svd(a, b, testCase);
    end
end

author: indigits

repo: sparse-plex


text:
% a simple demo of HOG template matching : Jianxiong

% you may need to compile some mex files in matlab: see features_compile.m

clear
close all

image = imread('demo_image.jpg');

feature = feature_pyramid(image);

%% pick a template
template_level = 1;
template_width = 8;
template_height= 6;
template_left  = 41;
tempalte_top   = 31;
% extract a tempalte
template = feature.feat{template_level}(tempalte_top:(tempalte_top+template_height-1),template_left:(template_left+template_width-1),:);
% normalize the tempalte to zero mean
% empirically better for some task, you can also try to disable this
template = template - mean(template(:));

% corresponding image
template_box = ([template_left tempalte_top template_left+size(template,2) tempalte_top+size(template,1)] - 1) * 8 / feature.scale(template_level) + [1 1 0 0];
template_box = round(template_box);
tempalte_image = image(template_box(2):template_box(4),template_box(1):template_box(3),:);

%% template matching
matchMatrix = [];

for level=1:length(feature.feat)
    % convolution
    match{level} = fconvblas(feature.feat{level}, {template}, 1, 1);
    match{level} = match{level}{1};

    % pick the top 5000 results
    [scores,indexes] = sort(match{level}(:),'descend');
    if numel(scores)>5000
        choosen = 1:5000;
        indexes = indexes(choosen);
        scores  = scores(choosen);
    end

    % obtain window pixel locations
    [indexes1 indexes2] = ind2sub(size(match{level}),indexes);
    bbs{level} = ([indexes2 indexes1 indexes2+size(template,2) indexes1+size(template,1)] - 1) * 8 / feature.scale(level) + repmat([1 1 0 0],numel(indexes),1);
    bbs{level} = [bbs{level} scores];

    % non maximal suppression to remove windows too close
    indexes = nmsMe(bbs{level}, 0.5); % pascal voc 0.5 criteria
    bbs{level} = bbs{level}(indexes,:);

    matchMatrix = [matchMatrix; bbs{level}];
end

% overall best match
indexes = nmsMe(matchMatrix, 0.5); % pascal voc 0.5 criteria
matchMatrix = matchMatrix(indexes,:);

%% visualization

visLevel = [];
for level = 1:3:length(feature.feat)
    if ~isempty(bbs{level})
        visLevel = [visLevel, level];
    end
end

subplot(length(visLevel)+1,3, 1);
showHOG(template);
title('the HOG tempalte');
subplot(length(visLevel)+1,3, 2);
imshow(tempalte_image);
title('image of the template');
subplot(length(visLevel)+1,3, 3);
imshow(image);
hold on
for b=1:min(5,size(matchMatrix))
    plot(matchMatrix(b,[1 3 3 1 1]),matchMatrix(b,[2 2 4 4 2]),'-y');
end
title('top 5 detection over all scale levels');

for l=1:length(visLevel)
    level = visLevel(l);
    subplot(length(visLevel)+1,3, 3*l+1);
    showHOG(feature.feat{level});
    title(sprintf('Level %d HOG ',level));

    subplot(length(visLevel)+1,3, 3*l+2);
    imagesc(match{level});
    axis image
    axis off
    title(sprintf('Level %d response ',level));

    subplot(length(visLevel)+1,3, 3*l+3);
    imshow(image);
    hold on
    for b=1:min(5,size(bbs{level},1)) % only visualization the top five results
        plot(bbs{level}(b,[1 3 3 1 1]),bbs{level}(b,[2 2 4 4 2]),'-y');
    end
    title(sprintf('Level %d top 5',level));
end

author: jianxiongxiao

repo: ProfXkit


text:
%Simulation test case
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Author: Jared Jacobs, jmjacobs@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 7/13/2010
classdef Simulation_KnowledgeBase_Test < TestCase
    methods
        function this = Simulation_KnowledgeBase_Test(name)
            this = this@TestCase(name);
        end

        function testInitializationFromKnowledgeBase(~)
            dbConnectionParameters = getConfig();
            database = edu.stanford.covert.db.MySQLDatabase(dbConnectionParameters);
            kbWID = edu.stanford.covert.cell.kb.KnowledgeBaseUtil.selectLatestKnowledgeBase(database);
            kb = edu.stanford.covert.cell.kb.KnowledgeBase(database, kbWID);

            sim = edu.stanford.covert.cell.sim.Simulation(kb.states, kb.processes);
            sim.initializeConstants(kb);
            sim.applyOptions('lengthSec', 10, 'verbosity', 0, 'seed', 1);
            fitter = edu.stanford.covert.cell.sim.util.FitConstants(sim, struct('method', 'heuristic', 'verbosity', 0));
            fitter.run();
            sim.run();

            assertElementsAlmostEqual(...
                1 / (9.0 * 3600) * log(2), sim.state('MetabolicReaction').growth(:,:,1),...
                'relative', 35e-2);

            database.close();
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function passed=moxunit_fieldtrip_runtests(varargin)
% run Fieldtrip tests with MOxUnit framework
%
% passed=moxunit_fieldtrip_runtests(...)
%
% Inputs:
%   'maxwalltime', w        Maximum duration as may be specified in each
%                           test; either numeric in seconds, or a string
%                           in format HH:MM:SS
%   'maxmem',m              Maximum memory usage as may be specified in
%                           each test; either numeric in bytes, or a string
%                           with a number and suffix 'kb', 'mb', 'gb', or
%                           'tb'
%   'upload',u              If 'true' or 'yes', then results are uploaded
%                           to the FieldTrip dashboard (default: true if
%                           the platform provides the 'weboptions' and
%                           'weboptions' functions)
%   'dependency',d          Do not run test if it matches a dependency in d
%   'xmloutput',f           Write JUnit-like XML file with test results to
%                           file f. This can be used with the shippable.com
%                           continuous integration service.
%   'exclude_if_prefix_equals_failed',e
%                           Exclude the test if the filename starts with
%                           'failed'.

%
% Output:
%   passed                  scalar boolean that is false if one or more
%                           tests failed or errored, and true otherwise.
%                           Tests that were skipped
%
% Examples:
%   % run all tests
%   moxunit_fieldtrip_runtests
%   %
%   % run tests with an expected execuation time of at most 10 minutes,
%   % maxmimum memory usage of 1 gigabyte, and with upload disabled.
%   moxunit_fieldtrip_runtests maxwalltime 00:10:00 maxmem 1gb upload no
%   %
%   % run all tests and store JUnit test results in file foo.xml. This file
%   % contains for each tests whether it passed, failed, errored, or was
%   % skipped; if it failed or errored, the full stack trace; if it was
%   % skipped, the reason why it was skipped.
%   moxunit_fieldtrip_runtests xmloutput foo.xml
%
% Notes:
% - This function aims to provide similar syntax as 'ft_test run', but
% provides an output
%
% #   For MOxUnit_fieldtrip's copyright information and license terms,   #
% #   see the COPYING file distributed with MOxUnit_fieldtrip.           #
%
% See also: ft_test

    passed=false; % fail unless all tests were successful

    % parse arguments
    opt=parse_args(varargin{:});

    % make a suite and add the tests
    suite=MOxUnitFieldTripTestSuite();
    suite=add_test_files(suite,opt.filelist);

    % show number of tests
    disp(suite);

    % apply filter args
    suite=filter(suite,opt.filter_args{:});

    % set output options
    verbosity=2;     % show detailed information
    output_stream=1; % standard output

    n_tests=countTestCases(suite);
    if n_tests==0
        fprintf(output_stream,['No test cases found - this is '...
                            'considered as a test suite failure\n']);
        return;
    end

    report=MOxUnitFieldTripTestReport(verbosity,output_stream);
    report=run(suite,report);

    if opt.upload
        sendToFieldtripDashboard(report);
    else
        fprintf(output_stream,'Not sending results to dashboard\n');
    end

    if ~isempty(opt.xmloutput)
        writeXML(report,opt.xmloutput);
        fprintf(output_stream,'XML test result written to ''%s''\n',...
                                                opt.xmloutput);
    end

    disp(report);
    passed=wasSuccessful(report);

function suite=add_test_files(suite,filelist)
% add tests specifed in filelist. If no tests are specified, then all
% tests in the 'test' directory are added
    if isempty(filelist)
        fieldtrip_test_dir=get_fieldtrip_test_dir();
        fieldtrip_test_file_pattern='^(test|failed).*\.m$';
        suite=addFromDirectory(suite,fieldtrip_test_dir,...
                                    fieldtrip_test_file_pattern);
    else
        % add tests from the FieldTrip test directory, so that file
        % names of tests do not need the full path to be added
        orig_pwd=pwd();
        pwd_resetter=onCleanup(@()cd(orig_pwd));
        cd(get_fieldtrip_test_dir());

        for k=1:numel(filelist)
            suite=addFromFile(suite,filelist{k});
        end
    end

function opt=parse_args(varargin)
    filter_args={};
    upload=[];
    filelist={};
    loadfile=[];
    xmloutput=[];

    pos=0;
    while true
        pos=pos+1;

        if pos>nargin
            break
        end

        key=varargin{pos};

        switch key
            case 'upload'
                [pos,value]=get_next_element(varargin,pos);
                upload=istrue(value);

            case 'loadfile'
                [pos,value]=get_next_element(varargin,pos);
                loadfile=istrue(value);

            case {'dependency','maxmem','maxwalltime',...
                        'exclude_if_prefix_equals_failed'}
                [pos,value]=get_next_element(varargin,pos);
                filter_args{end+1}=key;
                filter_args{end+1}=value;

            case 'xmloutput'
                [pos,value]=get_next_element(varargin,pos);
                xmloutput=value;

            otherwise
                filelist{end+1}=key;

        end
    end

    if isempty(loadfile)
        % 'loadfile' was not set as argument
        % set to true if at DCCN-like computer, false otherwise
        loadfile=exist(get_fieldtrip_data_dir(),'dir');
    end

    if isempty(upload)
        % TODO: decide whether upload should be disabled when not running
        % from DCCN-like computer

        % See if there are any upload methods available
        upload=~isempty(moxunit_fieldtrip_util_send_json());
    end

    filter_args=[filter_args,{'loadfile',loadfile}];

    opt=struct();
    opt.filter_args=filter_args;
    opt.upload=upload;
    opt.filelist=filelist;
    opt.xmloutput=xmloutput;

function [pos,value]=get_next_element(argcell,pos)

    if pos+1>numel(argcell)
        error('Missing value after key ''%s''', argcell{pos});
    end

    pos=pos+1;
    value=argcell{pos};

function result=ft_platform_supports_wrapper(arg)
    result=fieldtrip_function_wrapper(fullfile('utilities','private'),...
                                        'ft_platform_supports',arg);

function result=fieldtrip_function_wrapper(subdir,func_name,varargin)
    orig_pwd=pwd();
    cleaner=onCleanup(@()cd(orig_pwd));

    eval_dir=fullfile(get_fieldtrip_root_dir(),subdir);
    cd(eval_dir);

    func=str2func(func_name);
    result=func(varargin{:});

function fieldtrip_root_dir=get_fieldtrip_root_dir()
    fieldtrip_root_dir=fileparts(which('ft_defaults'));
    if isempty(fieldtrip_root_dir)
        error('Fieldtrip path is not set');
    end

function fieldtrip_test_dir=get_fieldtrip_test_dir()
    fieldtrip_root_dir=get_fieldtrip_root_dir();
    fieldtrip_test_dir=fullfile(fieldtrip_root_dir,'test');

function fieldtrip_data_dir=get_fieldtrip_data_dir()
    fieldtrip_data_dir='/home/common/matlab/fieldtrip/data/';
    fieldtrip_data_dir=fieldtrip_function_wrapper('test',...
                                'dccnpath',fieldtrip_data_dir);


text:
classdef MMF1 < PROBLEM
% <multi> <real> <multimodal>
% Multi-modal multi-objective test function

%------------------------------- Reference --------------------------------
% C. Yue, B. Qu, and J. Liang, A multi-objective particle swarm optimizer
% using ring topology for solving multimodal multiobjective Problems, IEEE
% Transactions on Evolutionary Computation, 2018, 22(5): 805-817.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties
        POS;    % Pareto optimal set for IGDX calculation
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            obj.D = 2;
            obj.lower    = [1,-1];
            obj.upper    = [3,1];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            PopObj(:,1) = abs(X(:,1)-2);
            PopObj(:,2) = 1-sqrt(PopObj(:,1))+2*(X(:,2)-sin(6*pi*PopObj(:,1)+pi)).^2;
        end
        %% Generate Pareto optimal solutions
        function R = GetOptimum(obj,N)
            % Generate points in Pareto optimal set
            obj.POS(:,1) = linspace(1,3,N)';
            obj.POS(:,2) = sin(6*pi*abs(obj.POS(:,1)-2)+pi);
            % Generate points on Pareto front
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R(:,1) = linspace(0,1,100)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case 'IGDX'
                    score = feval(metName,Population,obj.POS);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Display a population in the objective space
        function DrawObj(obj,Population)
            PopDec = Population.decs;
            temp   = PopDec(:,1)<=2;
            Draw(Population(temp).objs,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[1 .5 .5],'Markeredgecolor',[1 .2 .2],{'\it f\rm_1','\it f\rm_2',[]});
            Draw(Population(~temp).objs+0.1,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[.5 .5 1],'Markeredgecolor',[.2 .2 1]);
            Draw(obj.PF,'-','LineWidth',1,'Color',[1 .2 .2]);
            Draw(obj.PF+0.1,'-','LineWidth',1,'Color',[.2 .2 1]);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function quadrule_test105 ( )

%*****************************************************************************80
%
%% TEST105 tests JACOBI_EK_COMPUTE and JACOBI_SS_COMPUTE.
%
%  Discussion:
%
%    Compare with tabular values on page 178 of Stroud and Secrest.
%
%     In particular,
%
%             X              W
%
%     1  -0.9833999115   0.4615276287E-03
%     2  -0.9447138932   0.2732249104E-02
%     3  -0.8849310847   0.8045830455E-02
%    ..  .............   ................
%    19   0.9656375637   0.7613987785E-01
%    20   0.9934477866   0.3348337670E-01
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2011
%
%  Author:
%
%    John Burkardt
%
  order = 20;

  a = -1.0;
  b =  1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST105\n' );
  fprintf ( 1, '  JACOBI_EK_COMPUTE sets up Gauss-Jacobi quadrature;\n' );
  fprintf ( 1, '  JACOBI_SS_COMPUTE sets up Gauss-Jacobi quadrature;\n' );
  fprintf ( 1, '  Here, we simply compute a single rule and\n' );
  fprintf ( 1, '  print it, to check for accuracy.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The integration interval is [%f, %f]\n', a, b );

  alpha = 1.5;
  beta = 0.5;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  N = %d\n', order );
  fprintf ( 1, '  ALPHA = %f\n', alpha );
  fprintf ( 1, '  BETA =  %f\n', beta );

  [ x, w ] = jacobi_ek_compute ( order, alpha, beta );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  JACOBI_EK_COMPUTE:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I        X(I)            W(I)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : order
    fprintf ( 1, '  %4d  %14.8f  %14.8f\n', i, x(i), w(i) );
  end

  [ x, w ] = jacobi_ss_compute ( order, alpha, beta );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  JACOBI_SS_COMPUTE:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I        X(I)            W(I)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : order
    fprintf ( 1, '  %4d  %14.8f  %14.8f\n', i, x(i), w(i) );
  end

  return
end


text:
function [unquantized] = unQuantizeContinuous(values, min_v, max_v, num_bins)
    step_size = (max_v - min_v) / num_bins;
    unquantized = min_v + step_size/2 + (values-1) * step_size;

end

author: TadasBaltrusaitis

repo: OpenFace


text:
%go_calib_optim
%
%Main calibration function. Computes the intrinsic andextrinsic parameters.
%Runs as a script.
%
%INPUT: x_1,x_2,x_3,...: Feature locations on the images
%       X_1,X_2,X_3,...: Corresponding grid coordinates
%
%OUTPUT: fc: Camera focal length
%        cc: Principal point coordinates
%        alpha_c: Skew coefficient
%        kc: Distortion coefficients
%        KK: The camera matrix (containing fc and cc)
%        omc_1,omc_2,omc_3,...: 3D rotation vectors attached to the grid positions in space
%        Tc_1,Tc_2,Tc_3,...: 3D translation vectors attached to the grid positions in space
%        Rc_1,Rc_2,Rc_3,...: 3D rotation matrices corresponding to the omc vectors
%
%Method: Minimizes the pixel reprojection error in the least squares sense over the intrinsic
%        camera parameters, and the extrinsic parameters (3D locations of the grids in space)
%
%Note: If the intrinsic camera parameters (fc, cc, kc) do not exist before, they are initialized through
%      the function init_intrinsic_param.m. Otherwise, the variables in memory are used as initial guesses.
%
%Note: The row vector active_images consists of zeros and ones. To deactivate an image, set the
%      corresponding entry in the active_images vector to zero.
%
%VERY IMPORTANT: This function works for 2D and 3D calibration rigs, except for init_intrinsic_param.m
%that is so far implemented to work only with 2D rigs.
%In the future, a more general function will be there.
%For now, if using a 3D calibration rig, set quick_init to 1 for an easy initialization of the focal length

if ~exist('n_ima'),
   data_calib_no_read; % Load the images
   click_calib_no_read; % Extract the corners
end;

check_active_images;
check_extracted_images;
check_active_images;
desactivated_images = [];

recompute_extrinsic = (length(ind_active) < 100); % if there are too many images, do not spend time recomputing the extrinsic parameters twice..

%%% MAIN OPTIMIZATION CALL!!!!! (look into this function for the details of implementation)
go_calib_optim_iter;

if ~isempty(desactivated_images),
   param_list_save = param_list;
   fprintf(1,'\nNew optimization including the images that have been deactivated during the previous optimization.\n');
   active_images(desactivated_images) = ones(1,length(desactivated_images));
   desactivated_images = [];
   go_calib_optim_iter;
   if ~isempty(desactivated_images),
      fprintf(1,['List of images left desactivated: ' num2str(desactivated_images) '\n' ] );
   end;
   param_list = [param_list_save(:,1:end-1) param_list];
end;

repo: ekfmonoslam


text:
function [inface, outface]=innersurf(node,face,outface)
%
% outface=innersurf(node,face,outface)
%
% extract the interior triangles (shared by two enclosed compartments) of a complex surface
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%    node:  node coordinates
%    face:  surface triangle list
%    outface: (optional) the exterior triangle list, if not given, will
%           be computed using outersurf().
%
% output:
%    inface: the collection of interior triangles of the surface mesh
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin<3)
  outface=outersurf(node,face);
end

tf=ismember(sort(face,2),sort(outface,2),'rows');

inface=face(tf==0,:);

author: fangq


text:
function ica_input = compute_spatiotemporal_ica_input(pca_filters, pca_traces, mu)

	num_PCs = size(pca_traces,1);

	% Create concatenated input for spatio-temporal ICA
	if (mu == 1) % Pure temporal
	    ica_input = pca_traces;
	elseif (mu == 0) % Pure spatial
	    ica_input = pca_filters;
	else % Spatio-temporal
	    ica_input = [(1-mu)*pca_filters, mu*pca_traces];
	    for pc_idx = 1:num_PCs % Renormalize
	        ica_input_row = ica_input(pc_idx,:);
	        ica_input(pc_idx,:) = ica_input_row / sqrt(norm(ica_input_row));
	    end
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function pde = elli3DcircIntf3HM(am,ap,bm,bp,kappa,r1,r2,n1,n2)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,...
    'exactu1',@exactu1,'exactu2',@exactu2,'exactu3',@exactu3,...
    'um1',@um1,'um2',@um2,'um3',@um3,'up1',@up1,'up2',@up2,'up3',@up3,...
    'Dxu',@Dxu,'Dxum',@Dxum,'Dxup',@Dxup,'Dyu',@Dyu,...
    'Dyum',@Dyum,'Dyup',@Dyup,'Dzu',@Dzu,'Dzum',@Dzum,'Dzup',@Dzup,...
    'f1',@f1,'f2',@f2,'f3',@f3,...
    'fm1',@fm1,'fm2',@fm2,'fm3',@fm3,...
    'fp1',@fp1,'fp2',@fp2,'fp3',@fp3,...
    'A',@A,'Am',@Am,'Ap',@Ap,'one',@one,...
    'B',@B,'Bm',@Bm,'Bp',@Bp);

pde.am = am;
pde.ap = ap;
pde.bm = bm;
pde.bp = bp;
pde.kappa = kappa;
%% interface function
    function u = intf(x,y,z)
        u = (x.^2 + y.^2 + z.^2).^(1/2)/r1-1;
    end

%% exact solution
    function u = exactu1(x,y,z)
        u = um1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up1(x(id),y(id),z(id));
    end
    function u = exactu2(x,y,z)
        u = um2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up2(x(id),y(id),z(id));
    end
    function u = exactu3(x,y,z)
        u = um3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up3(x(id),y(id),z(id));
    end
    function u = um1(x,y,z)
        u = (x + n1*uker1(x,y,z).*(y-z))/am;
    end
    function u = um2(x,y,z)
        u = (y + n1*uker1(x,y,z).*(z-x))/am;
    end
    function u = um3(x,y,z)
        u = (z + n1*uker1(x,y,z).*(x-y))/am;
    end
    function u = up1(x,y,z)
        u = (x + n2*uker1(x,y,z).*uker2(x,y,z).*(y-z))/ap;
    end
    function u = up2(x,y,z)
        u = (y + n2*uker1(x,y,z).*uker2(x,y,z).*(z-x))/ap;
    end
    function u = up3(x,y,z)
        u = (z + n2*uker1(x,y,z).*uker2(x,y,z).*(x-y))/ap;
    end
    function u = uker1(x,y,z)
        u = r1^2 - (x.^2 + y.^2 + z.^2);
    end
    function u = uker2(x,y,z)
        u = r2^2 - (x.^2 + y.^2 + z.^2);
    end
%% Boundary Function
    function u = gD1(x,y,z)
        u = exactu1(x,y,z);
    end
    function u = gD2(x,y,z)
        u = exactu2(x,y,z);
    end
    function u = gD3(x,y,z)
        u = exactu3(x,y,z);
    end
%% Derivative of the exact solution
    function u = Dxu(x,y,z)
        u = Dxum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dxup(x(id),y(id),z(id));
    end
    function u = Dyu(x,y,z)
        u = Dyum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dyup(x(id),y(id),z(id));
    end
    function u = Dzu(x,y,z)
        u = Dzum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dzup(x(id),y(id),z(id));
    end
    function u = Dxum(x,y,z)
        u = -2*n1/am*(x.*y+x.*z-y.^2-z.^2 + uker1(x,y,z));
    end
    function u = Dyum(x,y,z)
        u = -2*n1/am*(y.*z+x.*y-z.^2-x.^2 + uker1(x,y,z));
    end
    function u = Dzum(x,y,z)
        u = -2*n1/am*(x.*z+y.*z-x.^2-y.^2 + uker1(x,y,z));
    end
    function u = Dxup(x,y,z)
        u = -2*n2*uker2(x,y,z)./ap.*(x.*y+x.*z-y.^2-z.^2 + uker1(x,y,z)) - ...
            2*n2*uker1(x,y,z).*(x.*y+x.*z - y.^2-z.^2)/ap;
    end
    function u = Dyup(x,y,z)
        u = -2*n2*uker2(x,y,z)./ap.*(y.*z+x.*y-z.^2-x.^2 + uker1(x,y,z)) - ...
            2*n2*uker1(x,y,z).*(y.*z+x.*y - x.^2-z.^2)/ap;
    end
    function u = Dzup(x,y,z)
        u = -2*n2*uker2(x,y,z)./ap.*(x.*z+y.*z-x.^2-y.^2 + uker1(x,y,z)) - ...
            2*n2*uker1(x,y,z).*(x.*z+y.*z - x.^2-y.^2)/ap;
    end

%% right hand side function
    function u = f1(x,y,z)
        u = fm1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp1(x(id),y(id),z(id));
    end
    function u = f2(x,y,z)
        u = fm2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp2(x(id),y(id),z(id));
    end
    function u = f3(x,y,z)
        u = fm3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp3(x(id),y(id),z(id));
    end

    function u = fm1(x,y,z)
        u = -10*n1*(z-y)- kappa*bm*um1(x,y,z);
    end
    function u = fm2(x,y,z)
        u = -10*n1*(x-z) - kappa*bm*um2(x,y,z);
    end
    function u = fm3(x,y,z)
        u = -10*n1*(y-x) - kappa*bm*um3(x,y,z);
    end
    function u = fp1(x,y,z)
        u = -10*n2*uker2(x,y,z).*(z-y) +8*n2*(z.^3-y.^3+z.*y.^2+z.*x.^2-y.*x.^2-y.*z.^2) -...
             10*n2*uker1(x,y,z).*(z-y) - kappa*bp*up1(x,y,z);
    end
    function u = fp2(x,y,z)
        u = -10*n2*uker2(x,y,z).*(x-z) +8*n2*(x.^3-z.^3+x.*z.^2+x.*y.^2-z.*y.^2-z.*x.^2) -...
             10*n2*uker1(x,y,z).*(x-z) - kappa*bp*up2(x,y,z);
    end
    function u = fp3(x,y,z)
        u = -10*n2*uker2(x,y,z).*(y-x) +8*n2*(y.^3-x.^3+y.*x.^2+y.*z.^2-x.*z.^2-x.*y.^2) -...
             10*n2*uker1(x,y,z).*(y-x) - kappa*bp*up3(x,y,z);
    end

%% Diffusion coefficient function
    function u = A(x,y,z)
        u = Am(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Ap(x(id),y(id),z(id));
    end
    function u = Am(x,y,z)
        u = am*ones(size(x));
    end
    function u = Ap(x,y,z)
        u = ap*ones(size(x));
    end

%% Mass coefficient function
    function u = B(x,y,z)
        u = Bm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Bp(x(id),y(id),z(id));
    end
    function u = Bm(x,y,z)
        u = bm*ones(size(x));
    end
    function u = Bp(x,y,z)
        u = bp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end

repo: ifem


text:
function trans = createRotation(varargin)
%CREATEROTATION Create the 3*3 matrix of a rotation.
%
%   TRANS = createRotation(THETA);
%   Returns the rotation corresponding to angle THETA (in radians)
%   The returned matrix has the form :
%   [cos(theta) -sin(theta)  0]
%   [sin(theta)  cos(theta)  0]
%   [0           0           1]
%
%   TRANS = createRotation(POINT, THETA);
%   TRANS = createRotation(X0, Y0, THETA);
%   Also specifies origin of rotation. The result is similar as performing
%   translation(-X0, -Y0), rotation(THETA), and translation(X0, Y0).
%
%   Example
%     % apply a rotation on a polygon
%     poly = [0 0; 30 0;30 10;10 10;10 20;0 20];
%     trans = createRotation([10 20], pi/6);
%     polyT = transformPoint(poly, trans);
%     % display the original and the rotated polygons
%     figure; hold on; axis equal; axis([-10 40 -10 40]);
%     drawPolygon(poly, 'k');
%     drawPolygon(polyT, 'b');
%
%   See also:
%   transforms2d, transformPoint, createRotation90, createTranslation
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 06/04/2004.
%

%   HISTORY
%   22/04/2009: rename as createRotation

% default values
cx = 0;
cy = 0;
theta = 0;

% get input values
if length(varargin)==1
    % only angle
    theta = varargin{1};
elseif length(varargin)==2
    % origin point (as array) and angle
    var = varargin{1};
    cx = var(1);
    cy = var(2);
    theta = varargin{2};
elseif length(varargin)==3
    % origin (x and y) and angle
    cx = varargin{1};
    cy = varargin{2};
    theta = varargin{3};
end

% compute coefs
cot = cos(theta);
sit = sin(theta);
tx =  cy*sit - cx*cot + cx;
ty = -cy*cot - cx*sit + cy;

% create transformation matrix
trans = [cot -sit tx; sit cot ty; 0 0 1];

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [X,Y,T,options,A,R2_pca,pca_opt,features] = preproc4hmm(X,Y,T,options)
% Prepare data to run TUDA
% 1. check parameters, including the type of classifier (regression is default)
% 2. Format X and Y accordingly to the classifier
% 3. Sets up state to be sequential , if asked
% 4. Preprocesses the data, including embedding

if length(size(X))==3 % 1st dim, time; 2nd dim, trials; 3rd dim, channels
    X = reshape(X,[size(X,1)*size(X,2), size(X,3)]);
end

N = length(T);
p = size(X,2);

if any(isnan(Y(:))), error('NaN found in Y'); end
if ~isfield(options,'Y_is_continuous') %ie not discrete
    options.Y_is_continuous = length(unique(Y(:))) > 50;
end

if size(X,1) ~= sum(T), error('Dimension of X not correct'); end
if size(Y,1) < size(Y,2); Y = Y'; end
if (size(Y,1) ~= sum(T)) && (size(Y,1) ~= length(T))
    error('Dimension of Y not correct');
end
q = size(Y,2);

if ~isfield(options,'K'), error('K needs to be specified'); end
if isfield(options,'downsample') && options.downsample
   error('Downsampling is not currently an option')
end
if isfield(options,'filter') && ~isempty(options.filter) && ~isfield(options,'Fs')
   error('You must specify options.Fs if you are to filter the data')
end

% options relative to the regression setting
if ~isfield(options,'classifier'), options.classifier = ''; end
% if empty, it is meant to be used with a continuous response
if ~isfield(options,'Nfeatures'), Nfeatures = p;
else, Nfeatures = options.Nfeatures; end
if ~isfield(options,'standardise'), standardise = 0;
else, standardise = options.standardise; end % otherwise it's done per trial
if ~isfield(options,'onpower'), onpower = 0;
else, onpower = options.onpower; end
if ~isfield(options,'embeddedlags'), embeddedlags = 0;
else, embeddedlags = options.embeddedlags; end
if ~isfield(options,'filter'), filter = [];
else, filter = options.filter; end
if ~isfield(options,'detrend'), detrend = 0;
else, detrend = options.detrend; end
% econ_embed saves memory at the expense of speed, when pca is applied
if ~isfield(options,'econ_embed'), econ_embed = 0;
else, econ_embed = options.econ_embed; end
if ~isfield(options,'parallel_trials'), parallel_trials = all(T==T(1)) & length(T)>1;
else, parallel_trials = options.parallel_trials; end
if ~isfield(options,'acrosstrial_constrained'), options.acrosstrial_constrained = 0; end

if ~isfield(options,'pca'), pca_opt = 0;
else, pca_opt = options.pca; end
if ~isfield(options,'A'), A = [];
else, A = options.A; end
if isfield(options,'downsample') && options.downsample~=0
    warning('Downsampling is not possible for TUDA')
end
if ~isfield(options,'encodemodel')
    options.encodemodel = false;
end
%options relative to classification models:
if ~isempty(options.classifier) || options.encodemodel
    if strcmp(options.classifier,'logistic')
        %set default options for logistic regression classification:
        options.distribution = 'logistic';
        demeanstim = false;
        %determine if multinomial or binomial:
        vals = unique(Y(:));
        if length(vals) == 2
            if all((vals == 0) | (vals == 1))
                Y = 2*(Y)-1;
            elseif q==1 && all((vals == 1) | (vals == 2))
                Y(Y==2) = -1;
            elseif any((vals ~= -1) & (vals ~= 1))
                error('Format of Y incorrect for classification tasks');
            end
        elseif length(vals) > 2 && q == 1
            % Y entered as categorical format
            fprintf(['\nFitting multinomial logistic classifier with classes ',int2str(vals'), '\n']);
            Y = convertToMultinomial(Y);
        elseif length(vals) == 3
            if any((vals ~= -1) & (vals ~= 1) & (vals ~= 0))
                error('Format of Y incorrect for classification tasks');
            end
        end
        options.logisticYdim=size(Y,2);
        if ~isfield(options,'balancedata')
            options.balancedata = 0;
        else
            options.balancedata = options.balancedata;
        end
        if ~isfield(options,'intercept'), options.intercept = 0; end
        if options.intercept
           X = [X,ones(size(X,1),1)];
        end
        options = rmfield(options,'intercept');
        if ~isfield(options,'sequential')
            options.sequential = 1;
        end
        if ~isfield(options,'inittype')
            if options.sequential ~= 0
                options.inittype = 'sequential';
            else
                options.inittype = 'HMM-MAR';
            end
        end
        add_noise = 0;
        if ~isfield(options,'cyc'), options.cyc = 1; end
   elseif strcmp(options.classifier,'LDA') || options.encodemodel
       % set default options for LDA model:
       options.distribution = 'Gaussian';
       demeanstim = false;
        if ~isfield(options,'intercept'), options.intercept = 1; end
        if options.intercept
            if size(Y,2)==1 && all((Y==0) + (Y==1))
                Y = [ones(size(Y,1),1),Y==1];
            else
                Y = [ones(size(Y,1),1),Y];
            end
            q = size(Y,2);
        end
        if ~isfield(options,'covtype')
            options.covtype = 'sharedfull';
        end
        options=rmfield(options,'intercept');
        if ~isfield(options,'sequential')
            options.sequential = 1;
        end
        if ~isfield(options,'inittype')
            if options.sequential ~= 0
                options.inittype = 'sequential';
            else
                options.inittype = 'HMM-MAR';
            end
        end
        options.add_noise = 0;
        add_noise = 0;
    elseif strcmp(options.classifier,'SVM') || strcmp(options.classifier,'SVM_rbf') || ...
            strcmp(options.classifier,'KNN') || ...
            strcmp(options.classifier,'decisiontree')
        add_noise = 0;
        demeanstim = false;
        options.sequential = 0;
    elseif strcmp(options.classifier,'regression')
        options.distribution = 'Gaussian';
        demeanstim = false;
        %determine if multinomial or binomial:
        vals = unique(Y(:));
        if length(vals) == 2 && q == 1
            if all((vals == 0) | (vals == 1))
                Y = 2*(Y)-1;
            elseif any((vals ~= -1) & (vals ~= 1))
                error('Format of Y incorrect for classification tasks');
            end
            c1 = sum(Y==-1); c2 = sum(Y==1); c12 = length(Y);
            Y(Y==-1) = - c12 / c1;  Y(Y==+1) = c12 / c2;
        elseif length(vals) > 2 && q == 1
            Ytmp = Y;
            Y = zeros(size(Y,1),length(vals));
            for jj = 1:length(vals), Y(Ytmp==vals(jj),jj) = 1; end
            q = length(vals);
        elseif q > 1
            if any((vals ~= 0) & (vals ~= 1))
                error('Format of Y incorrect for classification tasks');
            end
        end
        add_noise = ~options.Y_is_continuous;
        if ~isfield(options,'sequential')
            options.sequential = 1;
        end
        if ~isfield(options,'inittype')
            if options.sequential ~= 0
                options.inittype = 'sequential';
            else
                options.inittype = 'HMM-MAR';
            end
        end
    end

else % Standard regression problem

    options.distribution = 'Gaussian'; %default for all non-classification models
    demeanstim = true;
    if ~isfield(options,'add_noise'), add_noise = 0;
    else, add_noise = options.add_noise;
    end
    if ~isfield(options,'sequential')
        options.sequential = 1;
    end
    if ~isfield(options,'inittype')
        if options.sequential ~= 0
            options.inittype = 'fixedsequential';
        else
            options.inittype = 'HMM-MAR';
        end
    end
    if ~isfield(options,'intercept'), options.intercept = 0; end
    if options.intercept
       X = [X,ones(size(X,1),1)];
    end
    options = rmfield(options,'intercept');
    if ~isfield(options,'cyc'), options.cyc = 1; end
end

if ~isfield(options,'cyc'), options.cyc = 25; end
if ~isfield(options,'logisticYdim'), options.logisticYdim = 0; end

% Set up states to be a a sequence
if isfield(options,'sequential') && options.sequential ~= 0
    options.Pistructure = zeros(1,options.K);
    options.Pistructure(1:abs(options.sequential)) = 1;
    options.Pistructure = logical(options.Pistructure);
    options.Pstructure = logical(eye(options.K));
    for k = 1:abs(options.sequential)
        options.Pstructure = options.Pstructure + diag(ones(options.K-k,1),k);
    end
    if options.sequential < 0
        for k = 1:abs(options.sequential)
            options.Pstructure = options.Pstructure + diag(ones(options.K-k,1),-k);
        end
    end
end

% Options relative to constraints in the trans prob mat
if ~isfield(options,'K'), error('K was not specified'); end
if ~isfield(options,'Pstructure')
    options.Pstructure = true(options.K);
end
if ~isfield(options,'Pistructure')
    options.Pistructure = true(1,options.K);
end

options.parallel_trials = parallel_trials;
if ~isfield(options,'tudamonitoring'), options.tudamonitoring = 0; end
if ~isfield(options,'plotGamma'), options.plotGamma = 0; end

if parallel_trials && ~all(T==T(1))
    error('parallel_trials can be used only when all trials have equal length');
end
if options.tudamonitoring && ~all(T==T(1))
    error('tudamonitoring can be used only when all trials have equal length');
end

% options relative to the HMM
if ~isfield(options,'distribution'),options.distribution='Gaussian';end
if options.logisticYdim>0
    options.distribution = 'logistic';
end
if strcmp(options.distribution,'logistic')
    options.covtype = '';
end
if ~isfield(options,'covtype') && strcmp(options.distribution,'Gaussian')
    options.covtype = 'shareddiag';
end
if ~isfield(options,'inittype'), options.inittype = 'HMM-MAR'; end

options.order = 1;
options.zeromean = 1;
options.embeddedlags = 0; % it is done here
options.pca = 0; % it is done here
options.standardise = 0; % it is done here
options.onpower = 0; % it is done here
options.detrend = 0; % it is done here
options.filter = []; % it is done here
options.downsample = 0; % it is done here
options.dropstates = 0;

if isfield(options,'econ_embed'), options = rmfield(options,'econ_embed'); end
if isfield(options,'Nfeatures'), options = rmfield(options,'Nfeatures'); end
if isfield(options,'demeanstim'), options = rmfield(options,'demeanstim'); end
% Set a high prior for the initial probabilities because otherwise the
% model is biased to have the first time point of each trial to be assigned
% to just one state.
if ~isfield(options,'PriorWeightingPi'), options.PriorWeightingPi = length(T)/20; end
if ~isfield(options,'DirichletDiag'), options.DirichletDiag = 100; end

do_embedding = length(embeddedlags)>1;
do_pca = ~isempty(A) || length(pca_opt)>1 || (pca_opt>0 && pca_opt<(p*length(embeddedlags)));
do_pls = isfield(options,'pls');

if ~do_embedding && econ_embed
    econ_embed = 0;
end
if do_embedding && ~do_pca && econ_embed
    warning('It only makes sense to use econ_embed when using pca')
    econ_embed = 0;
end

emforw = max(max(embeddedlags),0);
emback = max(-min(embeddedlags),0);

if size(Y,1) == N % one value for the entire trial
    %Y = reshape(repmat(reshape(Y,[1 N q]),[ttrial 1 1]),[N*ttrial q]);
    Ytmp = Y;
    Y = zeros(sum(T),q);
    for n = 1:N
        Y(sum(T(1:n-1)) + (1:T(n)),:) = repmat(Ytmp(n,:),T(n),1);
    end; clear Ytmp
end

if q == 1 && length(unique(Y))==2
    if ismember(0,unique(Y)) || all(unique(Y)>0) && ~strcmp(options.distribution,'logistic')
        warning('Seems this is binary classification, transforming stimulus to have elements (-1,+1)')
        if islogical(Y);Y=1*Y;end
        v = unique(Y);
        Y(Y==v(1)) = -1; Y(Y==v(2)) = +1;
        if options.logisticYdim==0
            Y = Y - mean(Y);
        end
    end
end

if demeanstim
    % Demean stimulus
    Y = bsxfun(@minus,Y,mean(Y));
end
% Add noise, to avoid numerical problems
if add_noise > 0
    if add_noise == 1
        Y = Y + 1e-5 * randn(size(Y)) .* repmat(std(Y),size(Y,1),1);
    else
        Y = Y + add_noise * randn(size(Y)) .* repmat(std(Y),size(Y,1),1);
    end
end
% Standardise data
if standardise && N > 1
   warning(['You have set standardise=1, so each channel and trial will be standardized. ' ...
       'This will probably result in a loss of information in terms of how each stimulus is processed'])
   X = standardisedata(X,T,standardise);
end

% Filtering
if ~isempty(filter)
    data = filterdata(X,T,options.Fs,filter);
end
% Detrend data
if detrend
    X = detrenddata(X,T);
end

% adjust dimension of Y according to embedding
if do_embedding
    Ttmp = T-emforw-emback;
    Ytmp = Y;
    Y = zeros(sum(Ttmp),q);
    for n = 1:N
        Y(sum(Ttmp(1:n-1)) + (1:Ttmp(n)),:) = ...
            Ytmp(sum(T(1:n-1)) + (1+emforw:T(n)-emback) ,:);
    end; clear Ytmp Ttmp
end

% feature selection
if Nfeatures < p && Nfeatures > 0
    me = sum((repmat(mean(Y),size(Y,1),1) - Y).^2);
    C = zeros(p,1);
    for j=1:p
        if onpower
            x = rawsignal2power(X(:,j),T);
        else
            x = X(:,j);
        end
        if do_embedding
            x = embeddata(x,T,embeddedlags);
        end
        b = (x' * x) \ (x' * Y);
        e = sum((x * b - Y).^2);
        C(j) = sum(1 - e ./ me);
    end
    [~,features] = sort(C,1,'descend');
    features = features(1:Nfeatures);
    X = X(:,features);
    p = Nfeatures;
else
    features = 1:p;
end

% Hilbert envelope
if onpower
    X = rawsignal2power(X,T);
end

% do embedding + PCA
R2_pca = [];
if econ_embed

    % build gram matrix subject by subject
    for n = 1:N
        t = (1:T(n)) + sum(T(1:n-1));
        if do_embedding
            Xn = embeddata(X(t,:),T(n),embeddedlags);
        else
            Xn = X(t,:);
        end
        Xn = Xn - repmat(mean(Xn),size(Xn,1),1); % must center
        if n==1, C = zeros(size(Xn,2)); end
        C = C + Xn' * Xn;
    end

    % do SVD
    if isempty(A)
        [A,e,~] = svd(C);
        e = diag(e);
        e = cumsum(e)/sum(e);
        p = num_comp_pca(e,pca_opt);
        A = A(:,1:p);
        R2_pca = e(p);
    else
        R2_pca = []; p = size(A,2);
    end
    % eigendecompose subject by subject
    Xtmp = X; Ttmp = T;
    T = T-emforw-emback;
    X = zeros(sum(T),p);
    for n = 1:N
        t = (1:T(n)) + sum(T(1:n-1));
        ttmp = (1:Ttmp(n)) + sum(Ttmp(1:n-1));
        if do_embedding
            Xn = embeddata(Xtmp(ttmp,:),Ttmp(n),embeddedlags);
        else
            Xn = Xtmp(ttmp,:);
        end
        Xn = Xn - repmat(mean(Xn),size(Xn,1),1); % must center
        X(t,:) = Xn * A;
    end

else

    if do_embedding
        [X,T] = embeddata(X,T,embeddedlags);
        msg = '(embedded)';
    else
        msg = '';
    end
    if do_pca
        if isempty(A)
            [A,X,e] = pca(X);
            e = cumsum(e)/sum(e);
            p = num_comp_pca(e,pca_opt);
            R2_pca = e(p);
            X = X(:,1:p);
            A = A(:,1:p);
            fprintf('Working in PCA %s space, with %d components. \n',msg,p)
        else
            X = bsxfun(@minus,X,mean(X));
            X = X * A;
        end
    else
        R2_pca = 1;
    end
    if do_pls
        [X,A] = PLSdimreduce(X,Y,T,options.pls);
        options = rmfield(options,'pls');
    end
end

end

function ncomp = num_comp_pca(e,d)

if length(d)==1 && d<1
    ncomp = find(e>d,1);
elseif length(d)==1 && d>=1
    ncomp = d;
elseif length(d)==2 || (length(d)==3 && d(3)==1)
    ncomp = min(find(e>d(1),1),d(2));
elseif length(d)==3 && d(3)==2
    ncomp = max(find(e>d(1),1),d(2));
else
    error('pca parameters are wrongly specified')
end

end

author: OHBA-analysis

repo: HMM-MAR


text:
function har=arrow(x,y,z,vx,vy,vz,th,lr,c,hpar)
% har=arrow(x,y,z,vx,vy,vz,th,lr,c)
% har - cell array of handles to objects
% draw 3d arrow
% (x,y,z) - point of arrow begin
% (vx,vy,vz) - vector for arrow
% th - tip height
% lr - line radiuse
% c - color
% hpar - parent

% xyz=[x;y;z];
v=[vx;vy;vz];
vl=sqrt(v'*v); % lenght of arrow, each size will be normalazed to vl
if vl~=0
    %warning('length of arrow equal to 0');

    n=20; % discretization for circle in tip (tip is cone)
    dfi=2*pi/n;
    r=0.3*th; % radius in tip
    h=th; % height of tip
    v0=v*((vl-h))/vl; % rest part vector

    % draw tip of arrow

    % find some vector-perpendicular to v:
    if (vx==0)&&(vy==0)
        vp=[0;vz;-vy];
    else
        vp=[vy;-vx;0];
    end
    vp=r*vp/sqrt(vp'*vp);

    vpa=vp; % vectors for draw circle

    for fi=dfi:dfi:2*pi
        vcr=cross(v,vp);
        vcrn=vcr/sqrt(vcr'*vcr);
        vp=vp*cos(dfi)+r*sin(dfi)*vcrn;
        vpa=[vpa,vp];
    end

    X=x+v0(1)+vpa(1,:);
    Y=y+v0(2)+vpa(2,:);
    Z=z+v0(3)+vpa(3,:);

    % hold off
    % plot3(X,Y,Z);

    % add vertex of tip:
    X=[X;(x+vx)*ones(size(X))];
    Y=[Y;(y+vy)*ones(size(Y))];
    Z=[Z;(z+vz)*ones(size(Z))];

    hs=surf(X,Y,Z,'parent',hpar);
    har{1}=hs;
    set(hs,'FaceColor',c,'EdgeColor','none');
    har{2}=fill3(X(1,:),Y(1,:),Z(1,:),c,'parent',hpar);

    % draw line:
    vpal=lr*vpa/r;
    Xl=x+vpal(1,:);
    Yl=y+vpal(2,:);
    Zl=z+vpal(3,:);

    Xl=[Xl;x+v0(1)+vpal(1,:)];
    Yl=[Yl;y+v0(2)+vpal(2,:)];
    Zl=[Zl;z+v0(3)+vpal(3,:)];

    hs=surf(Xl,Yl,Zl,'parent',hpar);
    har{3}=hs;
    set(hs,'FaceColor',c,'EdgeColor','none');
    har{4}=fill3(Xl(1,:),Yl(1,:),Zl(1,:),c,'parent',hpar);
    har{5}=fill3(Xl(2,:),Yl(2,:),Zl(2,:),c,'parent',hpar);

else % zero length

    n=20; % discretization for circle in tip (tip is cone)
    dfi=2*pi/n;
    r=0.3*th; % radius in tip
    h=th; % height of tip

    % instead of tip of arrow:
    har{1}=spherei(x,y,z,r,c,1,hpar);
    crc=circle(0,0,r);
    har{2}=fill3(x+crc(:,1),y+crc(:,2),z+zeros(length(crc(:,1))),c,'parent',hpar,'EdgeColor','none');

    % insrtead of line:
    har{3}=spherei(x,y,z,r,c,1,hpar);
    har{4}=fill3(x+crc(:,1),y+crc(:,2),z+zeros(length(crc(:,1))),c,'parent',hpar,'EdgeColor','none');
    har{5}=fill3(x+crc(:,1),y+crc(:,2),z+zeros(length(crc(:,1))),c,'parent',hpar,'EdgeColor','none');
end

author: Sable

repo: mcbench-benchmarks


text:
function [P] = JacobiP(x,alpha,beta,N);

% function [P] = JacobiP(x,alpha,beta,N)
% Purpose: Evaluate Jacobi Polynomial of type (alpha,beta) > -1
%          (alpha+beta <> -1) at points x for order N and returns P[1:length(xp))]
% Note   : They are normalized to be orthonormal.

% Turn points into row if needed.
xp = x; dims = size(xp);
if (dims(2)==1) xp = xp'; end;

PL = zeros(N+1,length(xp));

% Initial values P_0(x) and P_1(x)
gamma0 = 2^(alpha+beta+1)/(alpha+beta+1)*gamma(alpha+1)*...
    gamma(beta+1)/gamma(alpha+beta+1);
PL(1,:) = 1.0/sqrt(gamma0);
if (N==0) P=PL'; return; end;
gamma1 = (alpha+1)*(beta+1)/(alpha+beta+3)*gamma0;
PL(2,:) = ((alpha+beta+2)*xp/2 + (alpha-beta)/2)/sqrt(gamma1);
if (N==1) P=PL(N+1,:)'; return; end;

% Repeat value in recurrence.
aold = 2/(2+alpha+beta)*sqrt((alpha+1)*(beta+1)/(alpha+beta+3));

% Forward recurrence using the symmetry of the recurrence.
for i=1:N-1
  h1 = 2*i+alpha+beta;
  anew = 2/(h1+2)*sqrt( (i+1)*(i+1+alpha+beta)*(i+1+alpha)*...
      (i+1+beta)/(h1+1)/(h1+3));
  bnew = - (alpha^2-beta^2)/h1/(h1+2);
  PL(i+2,:) = 1/anew*( -aold*PL(i,:) + (xp-bnew).*PL(i+1,:));
  aold =anew;
end;

P = PL(N+1,:)';
return

repo: Aero-matlab