text:
function [voitype]=dicomrt_checkvoitype(VOI)
% dicomrt_checkvoi(VOI)
%
% Check if the VOI has its own original format (set of 2D contours defined in the Z axis)
% or if it has been modified to contain outlines defined also in the X and Y planes.
%
% See also: dicomrt_checkinput, dicomrt_loadvoi
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

VOI_pointer=VOI{1,1};
VOI_array=VOI{2,1};

nvois=size(VOI_array,1);

for i=1:nvois
    n3d(i)=isempty(strfind(VOI_array{i,1},'3D'));
end

if find(n3d)
    voitype='2D';
else
    voitype='3D';
end


text:
%ISOLATE_AXES Isolate the specified axes in a figure on their own
%
% Examples:
%   fh = isolate_axes(ah)
%   fh = isolate_axes(ah, vis)
%
% This function will create a new figure containing the axes/uipanels
% specified, and also their associated legends and colorbars. The objects
% specified must all be in the same figure, but they will generally only be
% a subset of the objects in the figure.
%
% IN:
%    ah - An array of axes and uipanel handles, which must come from the
%         same figure.
%    vis - A boolean indicating whether the new figure should be visible.
%          Default: false.
%
% OUT:
%    fh - The handle of the created figure.

% Copyright (C) Oliver Woodford 2011-2013

% Thank you to Rosella Blatt for reporting a bug to do with axes in GUIs
% 16/3/2012 Moved copyfig to its own function. Thanks to Bob Fratantonio
% for pointing out that the function is also used in export_fig.m.
% 12/12/12 - Add support for isolating uipanels. Thanks to michael for
% suggesting it.
% 08/10/13 - Bug fix to allchildren suggested by Will Grant (many thanks!).
% 05/12/13 - Bug fix to axes having different units. Thanks to Remington
% Reid for reporting the issue.

function fh = isolate_axes(ah, vis)
% Make sure we have an array of handles
if ~all(ishandle(ah))
    error('ah must be an array of handles');
end
% Check that the handles are all for axes or uipanels, and are all in the same figure
fh = ancestor(ah(1), 'figure');
nAx = numel(ah);
for a = 1:nAx
    if ~ismember(get(ah(a), 'Type'), {'axes', 'uipanel'})
        error('All handles must be axes or uipanel handles.');
    end
    if ~isequal(ancestor(ah(a), 'figure'), fh)
        error('Axes must all come from the same figure.');
    end
end
% Tag the objects so we can find them in the copy
old_tag = get(ah, 'Tag');
if nAx == 1
    old_tag = {old_tag};
end
set(ah, 'Tag', 'ObjectToCopy');
% Create a new figure exactly the same as the old one
fh = copyfig(fh); %copyobj(fh, 0);
if nargin < 2 || ~vis
    set(fh, 'Visible', 'off');
end
% Reset the object tags
for a = 1:nAx
    set(ah(a), 'Tag', old_tag{a});
end
% Find the objects to save
ah = findall(fh, 'Tag', 'ObjectToCopy');
if numel(ah) ~= nAx
    close(fh);
    error('Incorrect number of objects found.');
end
% Set the axes tags to what they should be
for a = 1:nAx
    set(ah(a), 'Tag', old_tag{a});
end
% Keep any legends and colorbars which overlap the subplots
lh = findall(fh, 'Type', 'axes', '-and', {'Tag', 'legend', '-or', 'Tag', 'Colorbar'});
nLeg = numel(lh);
if nLeg > 0
    set([ah(:); lh(:)], 'Units', 'normalized');
    ax_pos = get(ah, 'OuterPosition');
    if nAx > 1
        ax_pos = cell2mat(ax_pos(:));
    end
    ax_pos(:,3:4) = ax_pos(:,3:4) + ax_pos(:,1:2);
    leg_pos = get(lh, 'OuterPosition');
    if nLeg > 1;
        leg_pos = cell2mat(leg_pos);
    end
    leg_pos(:,3:4) = leg_pos(:,3:4) + leg_pos(:,1:2);
    ax_pos = shiftdim(ax_pos, -1);
    % Overlap test
    M = bsxfun(@lt, leg_pos(:,1), ax_pos(:,:,3)) & ...
        bsxfun(@lt, leg_pos(:,2), ax_pos(:,:,4)) & ...
        bsxfun(@gt, leg_pos(:,3), ax_pos(:,:,1)) & ...
        bsxfun(@gt, leg_pos(:,4), ax_pos(:,:,2));
    ah = [ah; lh(any(M, 2))];
end
% Get all the objects in the figure
axs = findall(fh);
% Delete everything except for the input objects and associated items
delete(axs(~ismember(axs, [ah; allchildren(ah); allancestors(ah)])));
end

function ah = allchildren(ah)
ah = findall(ah);
if iscell(ah)
    ah = cell2mat(ah);
end
ah = ah(:);
end

function ph = allancestors(ah)
ph = [];
for a = 1:numel(ah)
    h = get(ah(a), 'parent');
    while h ~= 0
        ph = [ph; h];
        h = get(h, 'parent');
    end
end
end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [clout] = mask_princomp(clusters,varargin)
% :Usage:
% ::
%
%     [clusters] = mask_princomp(clusters,[behavioral score vector],[corr flag],[plotflag],[saveflag])
%
% This function is just like cluster_princomp, except that it works on the SET
% of activations in all clusters, rather than within each cluster.
%
% clusters is structure of clusters from tor_extract_rois.m
%
% behavioral vector is row vector of behavioral or other scores to correlate
%
% corr flag:  1 = work on correlations among voxels, 2 = work on covariance
%
% plotflag:   1 = yes, 0 = no.  plots.
%
% try this to test the program on random data:
% ::
%
%    cl(1).all_data = randn(23,30);cl(1).numVox = 30;cl = cluster_princomp(cl,EXPT.behavior,1,1);
%    cl(1).all_data(:,1:10) = cl(1).all_data(:,1:10) + 10; cl = cluster_princomp(cl,EXPT.behavior,1,1);
%    cl(1).all_data(:,25:30) = cl(1).all_data(:,25:30) + repmat((EXPT.behavior .* 3)',1,6);
%    cl(1).all_data(:,21:24) = cl(1).all_data(:,21:24) + repmat((1:23)',1,4);
%    cl = cluster_princomp(cl,EXPT.behavior,1,1);
%
% mean-center everything now:
% ::
%
%    cl.PCA = []; cl.all_data - cl.all_data - repmat(mean(cl.all_data),size(cl.all_data,1),1);
%    cl = cluster_princomp(cl,EXPT.behavior,1,1);
%
% add another correlated group:
% ::
%
%    cl.all_data(:,1:5) = cl.all_data(:,1:5) + repmat(rand(23,1)*5,1,5);
%    cl = cluster_princomp(cl,EXPT.behavior,1,1);
%
% if component scores are used and correlated with behavior, this means that the subjects
% tend to show the behavioral effect who also show the pattern associated with comp. x.
% this may mean high on a number of voxels, or high on some and low on others.
% the weights may be used to interpret what the components mean, and this can be done
% graphically.
%
% t-tests on component scores have ambiguous interpretations, because a high t-score
% may indicate negative values or close-to-zero values on some voxels.
% a component could have the interpretation, "high on this component means high on V1
% and low on V2."
%
% classifying voxels is done using cluster analysis (hierarchical, centroid linkage)
% on the voxels (observations) using the PCA weights (eigenvectors) as variables.
% This lets the clustering algorithm work in the reduced variable space with dimensionality
% equal to the number of components.
% The max number of clusters is restricted based on the gradient of the eigenvalues in the PCA
% maxclusters = 1 + the number of eigenvalues with gradient at least 20% of the initial drop
% from 1 to 2 eigenvalues.
%
% Requires clustering library in Matlab.
% Robust option also uses the robust PCA algorithm RAPCA,
% created by:
% Hubert, M., Rousseeuw, P.J., Verboven, S. (2002),
%  "A fast method for robust principal components with applications to chemometrics", by Mia Hubert, Peter J. Rousseeuw,
%  Chemometrics and Intelligent Laboratory Systems, 60, 101-111.
%

corrflag = 1; plotflag = 1; robustflag = 1; saveflag = 0;
if length(varargin) > 1, corrflag = varargin{2};, end
if length(varargin) > 2, plotflag = varargin{3};, end
if length(varargin) > 3, saveflag = varargin{4};, end

% get all the voxels and data together

CLU = clusters2CLU(clusters);
CLU.all_data = cat(2,clusters.all_data);

clear clusters;
clusters(1) = CLU;
i = 1;
subclusters{1} = []; clusters(i).PCA = [];

disp('Output saved in mask_princomp.out')
diary mask_princomp.out

a = CLU.all_data;

    if size(a,2) > 2    % must have 3 voxels to try clustering

    % -------------------------------------------------------------------------------
    % * All the real work is done here.  Compute pc's and clustering
    % -------------------------------------------------------------------------------

    if ~robustflag
        [clusters(i).PCA.pcomps,clusters(i).PCA.weights,clusters(i).PCA.eigval,clusters(i).PCA.class] = pc(a,corrflag);

        % automatically pick number of clusters, based on gradient in eigenvalues
        g = abs(gradient(clusters(i).PCA.eigval));
        maxclusters = sum(g > g(1).*.2) + 1;
    else

        % pick number of clusters by hand
        out = rapca(a);
        clusters(i).PCA.pcomps = out.T;
        clusters(i).PCA.weights = out.P;
        clusters(i).PCA.eigval = out.L;
        maxclusters = length(out.L);

        if saveflag,
            try mkdir mask_pca_results, catch, end
            saveas(gcf,'mask_pca_results/robust_pca_diagnostic','fig'), close
            saveas(gcf,'mask_pca_results/robust_eigenvalues','fig'), close
        end

    end

    clusters(i).PCA.class = docluster(a,maxclusters,plotflag);

    grps = unique(clusters(i).PCA.class(clusters(i).PCA.class~=0)); % values are component of origin
    for j = 1:length(grps),
        clusters(i).PCA.avgs(:,j) = mean(a(:,find(clusters(i).PCA.class==grps(j))),2);,
        freq(j) = sum(clusters(i).PCA.class == grps(j));
    end

    disp([clusters(i).title ', ' num2str(clusters(i).numVox) ' voxels: ' num2str(size(clusters(i).PCA.pcomps,2)) ' components'])
    fprintf(1,'\tMean\tEigval\tcorrel\t')

    % -------------------------------------------------------------------------------
    % * display each component and correlation with behavior
    % -------------------------------------------------------------------------------

    for j = 1:size(clusters(i).PCA.pcomps,2)

        %[H,P,CI,STATS] = TTEST(clusters(i).PCA.pcomps(:,j),0,.05,0);
        % skip the t-test.  t-tests on component scores don't make a lot of sense.
        fprintf(1,'\n\t%3.3f\t%3.3f\t',mean(clusters(i).PCA.pcomps(:,j)),clusters(i).PCA.eigval(j))

        if length(varargin) > 0
            if ~isempty(varargin{1})
                co = corrcoef(clusters(i).PCA.pcomps(:,j),varargin{1});
                co = co(1,2);
                fprintf(1,'%3.3f\t',co)
            end
        end

    end
    fprintf(1,'\n')

    % -------------------------------------------------------------------------------
    % * display classification info
    % -------------------------------------------------------------------------------
    disp(['Classified into ' num2str(max(clusters(i).PCA.class)) ' groups:'])
    fprintf(1,'\tClass\tVoxels\tMean\tt\tp\tcorrect. p\tcorrel\t')

    % for each component, test mean value and correlation with behavior
    for j = 1:length(grps)

        [H,P,CI,STATS] = TTEST(clusters(i).PCA.avgs(:,j),0,.05,0);
        fprintf(1,'\n\t%3.0f\t%3.0f\t%3.3f\t%3.3f\t%3.3f\t%3.3f\t',j,freq(j),mean(clusters(i).PCA.avgs(:,j)),STATS.tstat,P,P .* size(clusters(i).PCA.avgs,2))

        if length(varargin) > 0
            if ~isempty(varargin{1})
                co = corrcoef(clusters(i).PCA.avgs(:,j),varargin{1});
                co = co(1,2);
                fprintf(1,'%3.3f\t',co)
            end
        end

    end
    fprintf(1,'\n')

    % -------------------------------------------------------------------------------
    % * Plot, if requested
    % -------------------------------------------------------------------------------

    if plotflag,
        figure('Color','w'), subplot(1,3,1), imagesc(a), title(['Cl ' num2str(i) ': Data']), xlabel('Voxels'),ylabel('Subjects')
        subplot(1,3,2), imagesc(clusters(i).PCA.weights'), title(['Weights (eigenvectors)']), xlabel('Voxels'),ylabel('Eigenvectors')
        subplot(1,3,3), imagesc(clusters(i).PCA.pcomps), title(['Component scores (predictions)']), xlabel('Voxels'),ylabel('Subjects')

        if saveflag,
            try mkdir mask_pca_results, catch, end
            disp('Images saved in mask_pca_results folder')
            saveas(gcf,'mask_pca_results/PCA1','fig'),
        end

        a = [clusters(i).PCA.class' a']; a=sortrows(a,1); a = a(:,2:end)';
        figure;subplot 131; imagesc(a),title(['Cl ' num2str(i) ':Data sorted by class']), xlabel('Class'),ylabel('Subjects'),
        xlab = [sort(clusters(i).PCA.class(clusters(i).PCA.class~=0)) clusters(i).PCA.class(clusters(i).PCA.class==0)];
        set(gca,'XTick',1:length(clusters(i).PCA.class)); set(gca,'XTickLabel',xlab)
        subplot 132; imagesc(clusters(i).PCA.avgs),title('Class averages'), xlabel('Class'),ylabel('Subjects'),
        subplot 133; if length(varargin) > 0, if ~isempty(varargin{1}), imagesc(varargin{1}'), title('Behavior'),end,end

        if saveflag, saveas(gcf,'mask_pca_results/PCA_with_clustering','fig'),end
    end

    else
        disp(['Cluster ' num2str(i) ' has less than 3 voxels.'])
        clusters(i).PCA.class = ones(1,clusters(i).numVox);
        clusters(i).PCA.avgs = clusters(i).timeseries;
        grps = 1;
        if ~isfield(clusters,'correl'), clusters(1).correl = [];, end
    end

    diary off

    % -------------------------------------------------------------------------------
    % * separate into subclusters, based on class membership
    % -------------------------------------------------------------------------------

    for j = 1:length(grps)

        wh = find(clusters(i).PCA.class == grps(j));

        if length(varargin) > 0
            if ~isempty(varargin{1})
                co = corrcoef(clusters(i).PCA.avgs(:,j),varargin{1});
                subc(j).correl = co(1,2);
            end
        end

        CLUcomp = clusters(1);
        CLUcomp.XYZmm = CLUcomp.XYZmm(:,wh);
        CLUcomp.XYZ = CLUcomp.XYZ(:,wh);
        CLUcomp.Z = CLUcomp.Z(:,wh);
        CLUcomp.all_data = CLUcomp.all_data(:,wh);
        CLUcomp.numVox = size(CLUcomp.all_data,2);

        clout{j} = tor_extract_rois([],CLUcomp,CLUcomp);
        if plotflag
            montage_clusters([],clout{j})
            if saveflag, saveas(gcf,['mask_pca_results/montage_subset' num2str(j)],'fig'),end
        end

    end

end

return

function [b,v,d,class] = pc(a,corrflag)
% a is original matrix, b is principal components, v is eigenvectors
% (weights on columns, which = weights on voxels)
% class is classification of voxels into groups based on component loadings

if corrflag, [v,d]=eig(corrcoef(a));, else, [v,d]=eig(cov(a));,end
b = (pinv(v) * a')' ./ repmat((diag(d)').^.5,size(a,1),1);
% i made this up: think of rptating each subject's scores (in cols of a')
% by the rotation matrix pinv(v), and normalizing by the sqrt of the eigenvalues
% pinv(v) and v are rotation matrices because det = 1, no shearing or dilation
%
% this appears to work to give scores as well
% both methods (above,below) are scaled versions of the splus factor scores
% the problem is that doing it two different ways in splus flips the signs
% of some components and not others (gui vs cmd line).

%X = a; R = corrcoef(a); A = v * (d^.5); B = inv(R) * A;
%scores = X * B;
% X is data, A is factor loading matrix, B is factor score coeff matrix
% this method, from the text, and the one giving b above produce identical results

A = v * (d^.5);

b = fliplr(b); v = fliplr(v); A = fliplr(A); %scores = fliplr(scores);

num = min(10,sum(diag(d) >= 1));
b = b(:,1:num); v = v(:,1:num); A = A(:,1:num);
origd = diag(d);
d = diag(d)'; d= fliplr(d); d = d(1:num);

if num == 0, warning('No eigenvalues above 1!');, origd, class = [];

else
    % classify each voxel into a group based on loading
    % use A, which re-introduces the comp variance, because we
    % want relationships with more variance to count more.
    % This just doesn't work so hot.  See docluster, below.

    %wh = A' == repmat(max(A'),size(A,2),1);
    %for i = 1:size(wh,2), tmp = find(wh(:,i));, class(i) = tmp(1); end
    %class(max(A') < .3) = 0;

end

%figure;plot(b,'r'),hold on;plot(a,'k'), hold on; plot(mean(a,2),'g--'),legend({'eig' 'orig' 'avg'})

return

function class = docluster(a,maxclusters,doplot)

    if size(a,2) > 300,
        Y = pdist_by_parts(a');
    else
        Y = pdist(a','euclid');     % transpose so the voxels are observations, eigenvectors the variables
    end
    Z = linkage(Y,'centroid');
    class = cluster(Z,maxclusters)';

    if doplot,
        dendrogram(Z,0); title('Dendrogram for clustering')
    end

return

function pd =  pdist_by_parts(x)

disp(['Computing distances between observations'])

indx = 1;
vindx = 2;

for i = 1:size(x,1)
    indx = 1;   % faster this way
    for j = vindx:size(x,1)
        pd{i}(indx) = (sum((x(i,:) - x(j,:)).^2)).^.5;
        indx = indx+1;
    end
    vindx = vindx+1;
    if mod(i,100)==0,fprintf(1,'.'),end
end

pd = cat(2,pd{:});
fprintf(1,'done\n')

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function import_wizard( varargin )

if check_option(varargin,'EBSD')
  select = 2;
elseif check_option(varargin,'ODF')
  select = 3;
elseif check_option(varargin,'tensor')
  select = 4;
else
  select = 1;
end

api.presetDataImport    = select;
api.isPublishMode       = getMTEXpref('generatingHelpMode');

if check_option(varargin,'test')
  test_wizard_export(api,varargin{:});
  return
end

api = WizardDataApi(api);
api = WizardEmptyGUIApi(api);
api = WizardProgressApi(api,@pageImportData);
api.initProgressApi();

function test_wizard_export(api,varargin)

api = WizardDataApi(api);

% testcases

switch lower(get_option(varargin,'test',''))
  case '3d'

    api.loadDataFiles(1,fullfile(mtexDataPath,'EBSD','3dData','S00.ang'))
    api.loadDataFiles(1,fullfile(mtexDataPath,'EBSD','3dData','S01.ang'))
    api.loadDataFiles(1,fullfile(mtexDataPath,'EBSD','3dData','S02.ang'))
    api.Export.setOptEBSD('is3d',true);
    api.Export.setOptEBSD('Z',num2cell(1:3));

  case 'ebsd'

    api.loadDataFiles(1,fullfile(mtexDataPath,'EBSD','data.ctf'))

  case 'odf'

    api.loadDataFiles(2,fullfile(mtexDataPath,'ODF','odf.mtex'))

  case 'pf'

    api.loadDataFiles(4,fullfile(mtexDataPath,'PoleFigure','dubna','Q(10-11)(01-11)_amp.cnv'));
    api.loadDataFiles(4,fullfile(mtexDataPath,'PoleFigure','dubna','Q(02-21)_amp.cnv'));
    api.loadDataFiles(4,fullfile(mtexDataPath,'PoleFigure','dubna','Q(11-20)_amp.cnv'));

    data = api.getData();

    data{1} = set(data{1},'CS',crystalSymmetry('trigonal'));

    api.setData(data);

  case 'bg'

    api.loadDataFiles(4,fullfile(mtexDataPath,'PoleFigure','dubna','Q(02-21)_amp.cnv'));
    api.loadDataFiles(5,fullfile(mtexDataPath,'PoleFigure','dubna','Q(11-20)_amp.cnv'));

  case 'def'

    api.loadDataFiles(4,fullfile(mtexDataPath,'PoleFigure','dubna','Q(02-21)_amp.cnv'));
    api.loadDataFiles(5,fullfile(mtexDataPath,'PoleFigure','dubna','Q(11-20)_amp.cnv'));
    api.loadDataFiles(6,fullfile(mtexDataPath,'PoleFigure','dubna','Q(02-21)_amp.cnv'));
    api.loadDataFiles(7,fullfile(mtexDataPath,'PoleFigure','dubna','Q(11-20)_amp.cnv'));

  case 'lcp'

    api.loadDataFiles(4,fullfile(mtexDataPath,'PoleFigure','BearTex','Test_2.XPa'));
    api.loadDataFiles(4,fullfile(mtexDataPath,'PoleFigure','beartex.XPa'));

  otherwise

    api.loadDataFiles(4,fullfile(mtexDataPath,'PoleFigure','dubna','Q(02-21)_amp.cnv'));

    api = WizardEmptyGUIApi(api);
    api = WizardProgressApi(api,@pageImportData);
    api.initProgressApi();

    for j = 1:7
      pause(.5)
      api.nextCallback();
      snapnow;
    end

    close(api.hFigure);

    return
    % end

    % api = WizardProgressApi(api,@(a) pageCS(a,3));
    % api = WizardProgressApi(api,@pageKernel);
    % api = WizardProgressApi(api,@page3dEBSD);
    % api = WizardProgressApi(api,@pageFinish);

end

WizardFinish(api)

% api = WizardProgressApi(api,@(a) pageCS(a,3));
% api = WizardProgressApi(api,@pageKernel);
% api = WizardProgressApi(api,@page3dEBSD);
% api = WizardProgressApi(api,@pageFinish);

author: mtex-toolbox

repo: mtex

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


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the RUBBERBAND-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Rubberband()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  32;        % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny =  32;        % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
Ly = 1.0;        % Length of Eulerian Grid in y-Direction

% Immersed Structure Geometric / Dynamic Parameters %
N = 2*Nx;        % Number of Lagrangian Pts. (2x resolution of Eulerian grid)
a = 0.2;         % Length of semi-major axis.
b = 0.2;         % Length of semi-minor axis.
ds_Rest = 0;     % Resting length of springs
struct_name = 'rubberband'; % Name for .vertex, .spring, etc files.

% Call function to construct geometry
[xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(N,a,b);

% Plot Geometry to test
plot(xLag,yLag,'r-'); hold on;
plot(xLag,yLag,'*'); hold on;
xlabel('x'); ylabel('y');
axis square;

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .spring file!
%k_Spring = 1e7;
%print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name);

% Prints .beam file!
%k_Beam = 0.5; C = 0.0;
%print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name);

% Prints .target file!
k_Target = 1e7;
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints Vertex points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called rubberband.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C);
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called rubberband.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            else
                %Case s=N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 1,   k_Spring, ds_Rest);
            end
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(N,rmin,rmax)

% The immsersed structure is an ellipse %
t = 2*pi/N;
for i=1:N

    xLag(i) = 0.25 + rmax * cos( 2*pi/N*(i-1) );
    yLag(i) = 0.25 + rmin * sin( 2*pi/N*(i-1) );
end

author: nickabattista


text:
% Test file for chebtech/simplify.m

function pass = test_simplify(pref)

% Get preferences:
if ( nargin < 1 )
    pref = chebtech.techPref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
x = 2 * rand(100, 1) - 1;

% Tolerance for passing to simplify:
simptol = 1e-6;

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    %%
    % Test pathological inputs.

    % Empty CHEBTECH objects should be left alone.
    f = testclass.make();
    g = simplify(f);
    pass(n, 1) = isequal(f, g);

    % Unhappy CHEBTECH objects should be left alone.
    f = testclass.make(@(x) sqrt(x));
    g = simplify(f);
    pass(n, 2) = ~f.ishappy && isequal(f, g);

    %%
    % Test for a scalar-valued function:

    f = testclass.make(@(x) sin(100*(x + 0.1)));
    g = simplify(f, simptol);
    pass(n, 3) = abs(g.coeffs(end)) ~= 0;
    pass(n, 4) = length(g) < length(f);
    pass(n, 5) = norm(feval(f, x) - feval(g, x), inf) < 1e2*simptol*vscale(f);

    %%
    % Lengths of simplifications should be invariant under scaling:

    f1 = 1e-8*f;
    g1 = simplify(f1, simptol);
    pass(n, 6) = all(abs(g1.coeffs) >= eps*vscale(g1));
    pass(n, 7) = length(g1) == length(g);

    f2 = 1e8*f;
    g2 = simplify(f2, simptol);
    pass(n, 8) = all(abs(g2.coeffs) >= eps*vscale(g2));
    pass(n, 9) = length(g2) == length(g);

    %%
    % Test for an array-valued function:

    f = testclass.make(@(x) [sin(100*(x + 0.1)) cos(100*(x + 0.1)) exp(x)]);
    g = simplify(f, simptol);
    pass(n, 10) = any(abs(g.coeffs(1, :)) ~= 0);
    pass(n, 11) = length(g) < length(f);
    pass(n, 12) = all(norm(feval(f, x) - feval(g, x), inf) < ...
        10*max(simptol.*vscale(f)));

    %%
    % Try a contrived example which will leave a length 1 CHEBTECH:

    f = testclass.make(@(x) sin(100*(x + 0.1)));
    g = simplify(f, 1e20);
    pass(n, 13) = length(g) == 1;

    % Check that a long identically-zero CHEBTECH simplifies correctly:
    f = testclass.make(@(x) 0*x, [], struct('fixedLength', 8));
    g = simplify(f);
    pass(n, 14) = iszero(g) && (length(g) == 1);

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [metabolite_structure] = addMetFormulaCharge(metabolite_structure,startSearch,endSearch)
% This function uses getInchiString2ChargedFormula.m to calculate charge
% and neutral formula.
%
% INPUT
% metabolite_structure  metabolite structure
% startSearch           specify where the search should start in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
% endSearch             specify where the search should end in the
%                       metabolite structure. Must be numeric (optional, default: all metabolites
%                       in the structure will be search for)
% OUTPUT
% metabolite_structure  updated metabolite structure
%
% Ines Thiele 09/21
annotationSource = 'Calculated using metaboAnnotator and inchiString';
annotationType = 'automatic';

F = fieldnames(metabolite_structure);
if ~exist('startSearch','var')
    startSearch = 1;
end
if ~exist('endSearch','var')
    endSearch = length(F);
end

for i = startSearch : endSearch
    if (isempty(metabolite_structure.(F{i}).chargedFormula) || ~isempty(find(isnan(metabolite_structure.(F{i}).chargedFormula)))) && ...
            ~isempty(metabolite_structure.(F{i}).inchiString) && isempty(find(isnan(metabolite_structure.(F{i}).inchiString)))
        % compute charged formula for each entry from inchiString
        inchiString = metabolite_structure.(F{i}).inchiString;
        [metFormulaNeutral,metFormulaCharged,metCharge] = getInchiString2ChargedFormula({metabolite_structure.(F{i}).VMHId},cellstr(inchiString));
        if ~isempty(metFormulaCharged)
            metabolite_structure.(F{i}).chargedFormula = metFormulaCharged;
            metabolite_structure.(F{i}).chargedFormula_source = [annotationSource,':',annotationType,':',datestr(now)];
            metabolite_structure.(F{i}).charge = metCharge;
            metabolite_structure.(F{i}).charge_source = [annotationSource,':',annotationType,':',datestr(now)];
        end
        metabolite_structure.(F{i}).neutralFormula = metFormulaNeutral;
        metabolite_structure.(F{i}).neutralFormula_source = [annotationSource,':',annotationType,':',datestr(now)];
    end
end

author: opencobra

repo: cobratoolbox


text:
function submit(partId, webSubmit)
%SUBMIT Submit your code and output to the pgm-class servers
%   SUBMIT() will connect to the pgm-class server and submit your solution
%   There is no penalty for submitting, so go ahead and try this!
%
%   If this function does not work for you, use the web-submission mechanism.
%   Call the submitWeb function, which will produce a file for the part you
%   wish to submit. Then, submit the file to the class servers using the
%   "Web Submission" button on the Programming Assignments page on the course
%   website.
%
%   webSubmit is a boolean variable that specifies whether to prepare
%   a file for web-submission (if webSubmit = 1) or to submit directly
%   to the server (if webSubmit = 0, or is not specified); you should call
%   submitWeb if you want to do a web-submission.
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

fprintf('======================================================\n');
fprintf('                 Honor Code Statement:\n');
fprintf('------------------------------------------------------\n');
fprintf('The Honor Code is an undertaking of the students,\n');
fprintf('individually and collectively, that they will not\n');
fprintf('give or receive unpermitted aid in class work and\n');
fprintf('will do their share and take an active part in seeing\n');
fprintf('to it that others as well as themselves uphold the\n');
fprintf('spirit and letter of the Honor Code.\n\n');
fprintf('Sharing, collaboration, or looking at any code related\n');
fprintf('to Programming Assignments that is not your own is\n');
fprintf('considered a violation of the Honor Code.\n');
fprintf('------------------------------------------------------\n');
fprintf('By submitting the Programming Assignment, I\n');
fprintf('acknowledge and accept the Honor Code.\n');
fprintf('======================================================\n\n');

  fprintf('==\n== [pgm-class] Submitting Solutions | Programming Assignment %s\n==\n', ...
	  homework_id());

  if ~exist('partId', 'var') || isempty(partId)
    partId = promptPart();
  end

  if ~exist('webSubmit', 'var') || isempty(webSubmit)
    webSubmit = 0; % submit directly by default
  end

   % Check valid partId
  partNames = validParts();
  if ~isValidPartId(partId)
    fprintf('!! Invalid assignment part selected.\n');
    fprintf('!! Expected an integer from 1 to %d.\n', numel(partNames) + 1);
    fprintf('!! Submission Cancelled\n');
    return
  end

  if ~exist('pgm_login_data.mat','file')
    [login password] = loginPrompt();
    save('pgm_login_data.mat','login','password');
  else
    load('pgm_login_data.mat');
    [login password] = quickLogin(login,password);
    save('pgm_login_data.mat','login','password');
  end

  if isempty(login)
    fprintf('!! Submission Cancelled\n');
    return
  end

  fprintf('\n== Connecting to pgm-class ... ');
  if exist('OCTAVE_VERSION')
    fflush(stdout);
  end

  % Setup submit list
  if partId == numel(partNames) + 1
    submitParts = 1:numel(partNames);
  else
    submitParts = [partId];
  end

  for s = 1:numel(submitParts)
    % Submit this part
    partId = submitParts(s);

    for thisPartId = subParts(partId)
      if (~webSubmit) % submit directly to server

        [login, ch, signature, auxstring] = getChallenge(login, thisPartId);
        if isempty(login) || isempty(ch) || isempty(signature)
          % Some error occured, error string in first return element.
          fprintf('\n!! Error: %s\n\n', login);
          return
        end

        % Attempt Submission with Challenge
        ch_resp = challengeResponse(login, password, ch);

        [result, str] = submitSolution(login, ch_resp, thisPartId, ...
               output(thisPartId, auxstring), source(partId), signature);

        if (~isTest(thisPartId))
          partName = partNames{partId};
        else
          partName = [partNames{partId} ' (test)'];
        end

        fprintf('\n== [pgm-class] Submitted Assignment %s - Part %d - %s\n', ...
          homework_id(), partId, partName);
        fprintf('== %s\n', strtrim(str));

        if exist('OCTAVE_VERSION')
          fflush(stdout);
        end

      else % make web submission files

        [result] = submitSolutionWeb(login, thisPartId, output(thisPartId), ...
                            source(partId));
        result = base64encode(result);

        if (~isTest(thisPartId))
          partType = 'sample';
        else
          partType = 'test';
        end

        fprintf('\nSave as submission file [submit_pa%s_part%d_%s.txt (enter to accept default)]:', ...
          homework_id(), partId, partType);
        saveAsFile = input('', 's');
        if (isempty(saveAsFile))
          saveAsFile = sprintf('submit_pa%s_part%d_%s.txt', homework_id(), partId, partType);
        end

        fid = fopen(saveAsFile, 'w');
        if (fid)
          fwrite(fid, result);
          fclose(fid);
          fprintf('\nSaved your solutions to %s.\n\n', saveAsFile);
          fprintf(['You can now submit your solutions using the ' ...
             'Web Submission button\non the Assignments page\n']);
        else
          fprintf('Unable to save to %s\n\n', saveAsFile);
          fprintf(['You can create a submission file by saving the \n' ...
             'following text in a file: (press enter to continue)\n\n']);
          pause;
          fprintf(result);
        end

      end
    end
  end
end

% ====================== SERVER CONFIGURATION ===========================

% ***************** REMOVE -staging WHEN YOU DEPLOY *********************
function url = site_url()
  url = 'http://www.coursera.org/pgm';
end

function url = challenge_url()

  url = [site_url() '/assignment/challenge'];
end

function url = submit_url()
  url = [site_url() '/assignment/submit'];
end

% ================== CONFIGURABLES FOR EACH HOMEWORK ==================
function id = homework_id()
  id = '8';
end

function [partNames] = validParts()
  partNames = { 'FitGaussianParameters', ...
		'FitLinearGaussianParameters', ...
		'ComputeLogLikelihood', ...
        'LearnCPDsGivenGraph', ...
        'ClassifyDataset', ...
        'LearnGraphStructure', ...
        'LearnGraphAndCPDs'
	      };
end

function srcs = sources()
% Separated by part
  srcs = { { 'FitGaussianParameters.m' }, ... %1
	   { 'FitLinearGaussianParameters.m' }, ... %3
	   { 'ComputeLogLikelihood.m' }, ... %5
	   { 'LearnCPDsGivenGraph.m' }, ... %7
	   { 'ClassifyDataset.m' }, ... %9
	   { 'LearnGraphStructure.m' }, ... %11
	   { 'LearnGraphAndCPDs.m' } %13
	 };
end

% defines the shown part to back-end part mappings
function parts = subParts(partId)
  first = 2 * (partId - 1) + 1;
  parts = [first, first + 1];
end

% specifies which parts are test parts
function result = isTest(partId)
  if (mod(partId, 2) == 0)
      result = true;
  else
      result = false;
  end
end

function out = output(partId, auxstring)

  if partId == 1
    load 'PA8SampleCases.mat';
    [i j] = FitGaussianParameters(exampleINPUT.t1a1);
    out = [SerializeFloat(i) ' ' SerializeFloat(j) ' '];

  elseif partId == 2
    load 'submit_input.mat';
    [i j] = FitGaussianParameters(INPUT.t1a1);
    out = [SerializeFloat(i) ' ' SerializeFloat(j) ' '];
    clear 'INPUT';

  elseif partId == 3

    load 'PA8SampleCases.mat';
    [i j] = FitLinearGaussianParameters(exampleINPUT.t2a1, exampleINPUT.t2a2);
    out = '';
    for idx = 1:length(i)
        out = [out ' ' SerializeFloat(i(idx))];
    end
    out = [out ' ' num2str(j)];
  elseif partId == 4

    load 'submit_input.mat';
    [i j] = FitLinearGaussianParameters(INPUT.t2a1, INPUT.t2a2);
    out = '';
    for idx = 1:length(i)
        out = [out ' ' SerializeFloat(i(idx))];
    end
    out = [out ' ' num2str(j)];
    clear 'INPUT';

  elseif partId == 5

    load 'PA8SampleCases.mat';
    i = ComputeLogLikelihood(exampleINPUT.t3a1, exampleINPUT.t3a2, exampleINPUT.t3a3);
    out = SerializeFloat(i);

  elseif partId == 6

    load 'submit_input.mat';
    i = ComputeLogLikelihood(INPUT.t3a1, INPUT.t3a2, INPUT.t3a3);
    out = SerializeFloat(i);
    clear 'INPUT';

  elseif partId == 7

    load 'PA8SampleCases.mat';
    [P L] = LearnCPDsGivenGraph(exampleINPUT.t4a1, exampleINPUT.t4a2, exampleINPUT.t4a3);
    tmp = [[P.c] [P.clg.sigma_x] [P.clg.sigma_y] [P.clg.sigma_angle]];
    out = '';
    for idx = 1:length(tmp)
        out = [out ' ' SerializeFloat(tmp(idx))];
    end
    out = [out ' ' SerializeFloat(L)];

  elseif partId == 8

    load 'submit_input.mat';
    [P L] = LearnCPDsGivenGraph(INPUT.t4a1, INPUT.t4a2, INPUT.t4a3);
    tmp = [[P.c] [P.clg.sigma_x] [P.clg.sigma_y] [P.clg.sigma_angle]];
    out = '';
    for idx = 1:length(tmp)
        out = [out ' ' SerializeFloat(tmp(idx))];
    end
    out = [out ' ' SerializeFloat(L)];
    clear 'INPUT';

  elseif partId == 9
    load 'PA8SampleCases.mat';
    i = ClassifyDataset(exampleINPUT.t5a1, exampleINPUT.t5a2, exampleINPUT.t5a3, exampleINPUT.t5a4);
    out = SerializeFloat(i);

  elseif partId == 10

    load 'submit_input.mat';
    i = ClassifyDataset(INPUT.t5a1, INPUT.t5a2, INPUT.t5a3, INPUT.t5a4);
    out = SerializeFloat(i);
    clear 'INPUT';

  elseif partId == 11
    load 'PA8SampleCases.mat';
    [A W] = LearnGraphStructure(exampleINPUT.t6a1);
    out = '';
    for idx = 1:numel(A)
        out = [out ' ' num2str(A(idx))];
    end
  elseif partId == 12
    load 'submit_input.mat';
    [A W] = LearnGraphStructure(INPUT.t6a1);
    out = '';
    for idx = 1:numel(A)
        out = [out ' ' num2str(A(idx))];
    end
    clear 'INPUT';

  elseif partId == 13
    load 'PA8SampleCases.mat';
    [P G L] = LearnGraphAndCPDs(exampleINPUT.t7a1, exampleINPUT.t7a2);
    tmp = [[P.c] [P.clg.sigma_x] [P.clg.sigma_y] [P.clg.sigma_angle]];
    out = '';
    for idx = 1:numel(tmp)
        out = [out ' ' SerializeFloat(tmp(idx))];
    end
    for idx = 1:numel(G)
        out = [out ' ' SerializeFloat(G(idx))];
    end
    out = [out ' ' SerializeFloat(L)];
  elseif partId == 14
    load 'submit_input.mat';
    [P G L] = LearnGraphAndCPDs(INPUT.t7a1, INPUT.t7a2);
    tmp = [[P.c] [P.clg.sigma_x] [P.clg.sigma_y] [P.clg.sigma_angle]];
    out = '';
    for idx = 1:length(tmp)
        out = [out ' ' SerializeFloat(tmp(idx))];
    end
    for idx = 1:numel(G)
        out = [out ' ' SerializeFloat(G(idx))];
    end
    out = [out ' ' SerializeFloat(L)];
    clear 'INPUT';
  end

  % end of output function.

end

function out = SerializeFloat( x )
  out = sprintf('%.4f', x);
end

function out = SerializeFactorsFg(F)
% Serializes a factor struct array into the .fg format for libDAI
% http://cs.ru.nl/~jorism/libDAI/doc/fileformats.html
%
% To avoid incompatibilities with EOL markers, make sure you write the
% string to a file using the appropriate file type ('wt' for windows, 'w'
% for unix)

  lines = cell(5*numel(F) + 1, 1);

  lines{1} = sprintf('%d\n', numel(F));
  lineIdx = 2;
  for i = 1:numel(F)
    lines{lineIdx} = sprintf('\n%d\n', numel(F(i).var));
    lineIdx = lineIdx + 1;

    lines{lineIdx} = sprintf('%s\n', num2str(F(i).var(:)')); % ensure that we put in a row vector
    lineIdx = lineIdx + 1;

    lines{lineIdx} = sprintf('%s\n', num2str(F(i).card(:)')); % ensure that we put in a row vector
    lineIdx = lineIdx + 1;

    lines{lineIdx} = sprintf('%d\n', numel(F(i).val));
    lineIdx = lineIdx + 1;

    % Internal storage of factor vals is already in the same indexing order
    % as what libDAI expects, so we don't need to convert the indices.
    vals = [0:(numel(F(i).val) - 1); F(i).val(:)'];
    lines{lineIdx} = sprintf('%d %0.8g\n', vals);
    lineIdx = lineIdx + 1;
  end

  out = sprintf('%s', lines{:});

end

function out = SerializeMEUOptimization(meu, optdr)
  optdr = SortFactorVars(optdr);
  optdr_part = SerializeFactorsFg(optdr);
  out = sprintf('%s\n%.4f\n', optdr_part, meu);
end

function f = SortAllFactors(factors)

  for i = 1:length(factors)
    factors(i) = SortFactorVars(factors(i));
  end

  varMat = vertcat(factors(:).var);
  [unused, order] = sortrows(varMat);

  f = factors(order);

end

function G = SortFactorVars(F)

  [sortedVars, order] = sort(F.var);
  G.var = sortedVars;

  G.card = F.card(order);
  G.val = zeros(numel(F.val), 1);

  assignmentsInF = IndexToAssignment(1:numel(F.val), F.card);
  assignmentsInG = assignmentsInF(:,order);
  G.val(AssignmentToIndex(assignmentsInG, G.card)) = F.val;

end

% ========================= CHALLENGE HELPERS =========================

function src = source(partId)
  src = '';
  src_files = sources();
  if partId <= numel(src_files)
      flist = src_files{partId};
      for i = 1:numel(flist)
          fid = fopen(flist{i});
          if (fid == -1)
            error('Error opening %s (is it missing?)', flist{i});
          end
          line = fgets(fid);
          while ischar(line)
            src = [src line];
            line = fgets(fid);
          end
          fclose(fid);
          src = [src '||||||||'];
      end
  end
end

function ret = isValidPartId(partId)
  partNames = validParts();
  ret = (~isempty(partId)) && (partId >= 1) && (partId <= numel(partNames) + 1);
end

function partId = promptPart()
  fprintf('== Select which part(s) to submit:\n', ...
          homework_id());
  partNames = validParts();
  srcFiles = sources();
  for i = 1:numel(partNames)
    fprintf('==   %d) %s [', i, partNames{i});
    fprintf(' %s ', srcFiles{i}{:});
    fprintf(']\n');
  end
  fprintf('==   %d) All of the above \n==\nEnter your choice [1-%d]: ', ...
          numel(partNames) + 1, numel(partNames) + 1);
  selPart = input('', 's');
  partId = str2num(selPart);
  if ~isValidPartId(partId)
    partId = -1;
  end
end

function [email,ch,signature,auxstring] = getChallenge(email, part)
  str = urlread(challenge_url(), 'post', {'email_address', email, 'assignment_part_sid', [homework_id() '-' num2str(part)], 'response_encoding', 'delim'});

  str = strtrim(str);
  r = struct;
  while(numel(str) > 0)
    [f, str] = strtok (str, '|');
    [v, str] = strtok (str, '|');
    r = setfield(r, f, v);
  end

  email = getfield(r, 'email_address');
  ch = getfield(r, 'challenge_key');
  signature = getfield(r, 'state');
  auxstring = getfield(r, 'challenge_aux_data');
end

function [result, str] = submitSolutionWeb(email, part, output, source)

  result = ['{"assignment_part_sid":"' base64encode([homework_id() '-' num2str(part)], '') '",' ...
            '"email_address":"' base64encode(email, '') '",' ...
            '"submission":"' base64encode(output, '') '",' ...
            '"submission_aux":"' base64encode(source, '') '"' ...
            '}'];
  str = 'Web-submission';
end

function [result, str] = submitSolution(email, ch_resp, part, output, ...
                                        source, signature)

  params = {'assignment_part_sid', [homework_id() '-' num2str(part)], ...
            'email_address', email, ...
            'submission', base64encode(output, ''), ...
            'submission_aux', base64encode(source, ''), ...
            'challenge_response', ch_resp, ...
            'state', signature};

  str = urlread(submit_url(), 'post', params);

  % Parse str to read for success / failure
  result = 0;

end

% =========================== LOGIN HELPERS ===========================

function [login password] = loginPrompt()
  % Prompt for password
  [login password] = basicPrompt();

  if isempty(login) || isempty(password)
    login = []; password = [];
  end
end

function [login password] = basicPrompt()
  login = input('Login (Email address): ', 's');
  password = input('Submission Password (from Assignments page): ', 's');
end

function [login password] = quickLogin(login,password)
  cont_token = input(['You are currently logged in as ' login '.\nIs this you? (y/n - type n to reenter password)'],'s');
  if(isempty(cont_token) || cont_token(1)=='Y'||cont_token(1)=='y')
    return;
  else
    [login password] = loginPrompt();
  end
end

function [str] = challengeResponse(email, passwd, challenge)
  str = sha1([challenge passwd]);
end

% =============================== SHA-1 ================================

function hash = sha1(str)

  % Initialize variables
  h0 = uint32(1732584193);
  h1 = uint32(4023233417);
  h2 = uint32(2562383102);
  h3 = uint32(271733878);
  h4 = uint32(3285377520);

  % Convert to word array
  strlen = numel(str);

  % Break string into chars and append the bit 1 to the message
  mC = [double(str) 128];
  mC = [mC zeros(1, 4-mod(numel(mC), 4), 'uint8')];

  numB = strlen * 8;
  if exist('idivide')
    numC = idivide(uint32(numB + 65), 512, 'ceil');
  else
    numC = ceil(double(numB + 65)/512);
  end
  numW = numC * 16;
  mW = zeros(numW, 1, 'uint32');

  idx = 1;
  for i = 1:4:strlen + 1
    mW(idx) = bitor(bitor(bitor( ...
                  bitshift(uint32(mC(i)), 24), ...
                  bitshift(uint32(mC(i+1)), 16)), ...
                  bitshift(uint32(mC(i+2)), 8)), ...
                  uint32(mC(i+3)));
    idx = idx + 1;
  end

  % Append length of message
  mW(numW - 1) = uint32(bitshift(uint64(numB), -32));
  mW(numW) = uint32(bitshift(bitshift(uint64(numB), 32), -32));

  % Process the message in successive 512-bit chs
  for cId = 1 : double(numC)
    cSt = (cId - 1) * 16 + 1;
    cEnd = cId * 16;
    ch = mW(cSt : cEnd);

    % Extend the sixteen 32-bit words into eighty 32-bit words
    for j = 17 : 80
      ch(j) = ch(j - 3);
      ch(j) = bitxor(ch(j), ch(j - 8));
      ch(j) = bitxor(ch(j), ch(j - 14));
      ch(j) = bitxor(ch(j), ch(j - 16));
      ch(j) = bitrotate(ch(j), 1);
    end

    % Initialize hash value for this ch
    a = h0;
    b = h1;
    c = h2;
    d = h3;
    e = h4;

    % Main loop
    for i = 1 : 80
      if(i >= 1 && i <= 20)
        f = bitor(bitand(b, c), bitand(bitcmp(b), d));
        k = uint32(1518500249);
      elseif(i >= 21 && i <= 40)
        f = bitxor(bitxor(b, c), d);
        k = uint32(1859775393);
      elseif(i >= 41 && i <= 60)
        f = bitor(bitor(bitand(b, c), bitand(b, d)), bitand(c, d));
        k = uint32(2400959708);
      elseif(i >= 61 && i <= 80)
        f = bitxor(bitxor(b, c), d);
        k = uint32(3395469782);
      end

      t = bitrotate(a, 5);
      t = bitadd(t, f);
      t = bitadd(t, e);
      t = bitadd(t, k);
      t = bitadd(t, ch(i));
      e = d;
      d = c;
      c = bitrotate(b, 30);
      b = a;
      a = t;

    end
    h0 = bitadd(h0, a);
    h1 = bitadd(h1, b);
    h2 = bitadd(h2, c);
    h3 = bitadd(h3, d);
    h4 = bitadd(h4, e);

  end

  hash = reshape(dec2hex(double([h0 h1 h2 h3 h4]), 8)', [1 40]);

  hash = lower(hash);

end

function ret = bitadd(iA, iB)
  ret = double(iA) + double(iB);
  ret = bitset(ret, 33, 0);
  ret = uint32(ret);
end

function ret = bitrotate(iA, places)
  t = bitshift(iA, places - 32);
  ret = bitshift(iA, places);
  ret = bitor(ret, t);
end

% =========================== Base64 Encoder ============================
% Thanks to Peter John Acklam
%

function y = base64encode(x, eol)
%BASE64ENCODE Perform base64 encoding on a string.
%
%   BASE64ENCODE(STR, EOL) encode the given string STR.  EOL is the line ending
%   sequence to use; it is optional and defaults to '\n' (ASCII decimal 10).
%   The returned encoded string is broken into lines of no more than 76
%   characters each, and each line will end with EOL unless it is empty.  Let
%   EOL be empty if you do not want the encoded string broken into lines.
%
%   STR and EOL don't have to be strings (i.e., char arrays).  The only
%   requirement is that they are vectors containing values in the range 0-255.
%
%   This function may be used to encode strings into the Base64 encoding
%   specified in RFC 2045 - MIME (Multipurpose Internet Mail Extensions).  The
%   Base64 encoding is designed to represent arbitrary sequences of octets in a
%   form that need not be humanly readable.  A 65-character subset
%   ([A-Za-z0-9+/=]) of US-ASCII is used, enabling 6 bits to be represented per
%   printable character.
%
%   Examples
%   --------
%
%   If you want to encode a large file, you should encode it in chunks that are
%   a multiple of 57 bytes.  This ensures that the base64 lines line up and
%   that you do not end up with padding in the middle.  57 bytes of data fills
%   one complete base64 line (76 == 57*4/3):
%
%   If ifid and ofid are two file identifiers opened for reading and writing,
%   respectively, then you can base64 encode the data with
%
%      while ~feof(ifid)
%         fwrite(ofid, base64encode(fread(ifid, 60*57)));
%      end
%
%   or, if you have enough memory,
%
%      fwrite(ofid, base64encode(fread(ifid)));
%
%   See also BASE64DECODE.

%   Author:      Peter John Acklam
%   Time-stamp:  2004-02-03 21:36:56 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   if isnumeric(x)
      x = num2str(x);
   end

   % make sure we have the EOL value
   if nargin < 2
      eol = sprintf('\n');
   else
      if sum(size(eol) > 1) > 1
         error('EOL must be a vector.');
      end
      if any(eol(:) > 255)
         error('EOL can not contain values larger than 255.');
      end
   end

   if sum(size(x) > 1) > 1
      error('STR must be a vector.');
   end

   x   = uint8(x);
   eol = uint8(eol);

   ndbytes = length(x);                 % number of decoded bytes
   nchunks = ceil(ndbytes / 3);         % number of chunks/groups
   nebytes = 4 * nchunks;               % number of encoded bytes

   % add padding if necessary, to make the length of x a multiple of 3
   if rem(ndbytes, 3)
      x(end+1 : 3*nchunks) = 0;
   end

   x = reshape(x, [3, nchunks]);        % reshape the data
   y = repmat(uint8(0), 4, nchunks);    % for the encoded data

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Split up every 3 bytes into 4 pieces
   %
   %    aaaaaabb bbbbcccc ccdddddd
   %
   % to form
   %
   %    00aaaaaa 00bbbbbb 00cccccc 00dddddd
   %
   y(1,:) = bitshift(x(1,:), -2);                  % 6 highest bits of x(1,:)

   y(2,:) = bitshift(bitand(x(1,:), 3), 4);        % 2 lowest bits of x(1,:)
   y(2,:) = bitor(y(2,:), bitshift(x(2,:), -4));   % 4 highest bits of x(2,:)

   y(3,:) = bitshift(bitand(x(2,:), 15), 2);       % 4 lowest bits of x(2,:)
   y(3,:) = bitor(y(3,:), bitshift(x(3,:), -6));   % 2 highest bits of x(3,:)

   y(4,:) = bitand(x(3,:), 63);                    % 6 lowest bits of x(3,:)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Now perform the following mapping
   %
   %   0  - 25  ->  A-Z
   %   26 - 51  ->  a-z
   %   52 - 61  ->  0-9
   %   62       ->  +
   %   63       ->  /
   %
   % We could use a mapping vector like
   %
   %   ['A':'Z', 'a':'z', '0':'9', '+/']
   %
   % but that would require an index vector of class double.
   %
   z = repmat(uint8(0), size(y));
   i =           y <= 25;  z(i) = 'A'      + double(y(i));
   i = 26 <= y & y <= 51;  z(i) = 'a' - 26 + double(y(i));
   i = 52 <= y & y <= 61;  z(i) = '0' - 52 + double(y(i));
   i =           y == 62;  z(i) = '+';
   i =           y == 63;  z(i) = '/';
   y = z;

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Add padding if necessary.
   %
   npbytes = 3 * nchunks - ndbytes;     % number of padding bytes
   if npbytes
      y(end-npbytes+1 : end) = '=';     % '=' is used for padding
   end

   if isempty(eol)

      % reshape to a row vector
      y = reshape(y, [1, nebytes]);

   else

      nlines = ceil(nebytes / 76);      % number of lines
      neolbytes = length(eol);          % number of bytes in eol string

      % pad data so it becomes a multiple of 76 elements
      y = [y(:) ; zeros(76 * nlines - numel(y), 1)];
      y(nebytes + 1 : 76 * nlines) = 0;
      y = reshape(y, 76, nlines);

      % insert eol strings
      eol = eol(:);
      y(end + 1 : end + neolbytes, :) = eol(:, ones(1, nlines));

      % remove padding, but keep the last eol string
      m = nebytes + neolbytes * (nlines - 1);
      n = (76+neolbytes)*nlines - neolbytes;
      y(m+1 : n) = '';

      % extract and reshape to row vector
      y = reshape(y, 1, m+neolbytes);

   end

   % output is a character array
   y = char(y);

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function x = lprec1(c, d, h, g, extmod)
% LPREC1   One-level Laplacian pyramid reconstruction
%
%	x = lprec1(c, d, h, g)
%
% Input:
%   c:      coarse signal at half size
%   d:      detail signal at full size
%   h, g:   two biorthogonal 1-D lowpass filters
%   extmod: [optional] extension mode (default is 'per')
%
% Output:
%   x:      reconstructed signal
%
% Note:     This uses a new reconstruction method by Do and Vetterli,
%           "Framming pyramids", IEEE Trans. on Sig Proc., Sep. 2003.
%
% See also:	LPDEC1

if ~exist('extmod', 'var')
    extmod = 'per';
end

nd = ndims(c);

% First, filter and downsample the detail image
r = d;
for dim = 1:nd
    r = filtdn(r, h, dim, extmod, 0);
end

% Then subtract the result from the coarse signal
p = c - r;

% Even size filter needs to be adjusted to obtain perfect reconstruction
adjust = mod(length(g) + 1, 2);

% Then upsample and filter
for dim = 1:nd
    p = upfilt(p, g, dim, extmod, adjust);
end

% Final combination
x = p + d;

author: Sable

repo: mcbench-benchmarks


text:
function value = indexn0 ( n, i_min, i, i_max )

%*****************************************************************************80
%
%% INDEXN0 indexes an N-dimensional array by rows, with zero base.
%
%  Discussion:
%
%    Entries of the array are indexed starting at entry
%      ( I_MIN(1), I_MIN(2),...,I_MIN(N) ),
%    and increasing the last index up to I_MAX(N),
%    then the next-to-last and so on.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of indices.
%
%    Input, integer I_MIN(N), the minimum indices.
%
%    Input, integer I(N), the indices.
%
%    Input, integer I_MAX(N), for maximum indices.
%
%    Output, integer VALUE, the index of element I.
%
  index_min = 0;

  value = ( i(1) - i_min(1) );

  for j = 2 : n
    value = value * ( i_max(j) + 1 - i_min(j) ) + ( i(j) - i_min(j) );
  end
  value = value + index_min;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% op_CSImakeK0fid.m
% Jamie Near & Sneha Senthil, Sunnybrook Research Insitute 2021.
%
% USAGE:
% out=op_CSImakeK0fid(in,start,interval);
%
% DESCRIPTION:
% Extract points with repeating periodic interval from a spatio-spectral
% MRSI signal.  This is useful for Rosette MRSI data, when we want to extract
% the FIDs from the centre of k-space.
%
% INPUTS:
% in        = input MRSI data in matlab structure format.
% start     = The index of the first point to extract.
% interval  = The interval between subsequently extracted points.
%
% OUTPUTS:
% out       = Extracted FID.

function out=op_CSImakeK0fid(in,start,interval);

%extract the desired points;
fids=in.data(start:interval:end,:,:,:);

%recalculate the size;
sz=size(fids);

%remake the time vector:
t=in.adcTime(start:interval:end);

%dwelltime
dwelltime=t(2)-t(1);

%FILLING IN DATA STRUCTURE
out=in;
out.fids=fids;
out.sz=sz;
out.t=t;
out.dwelltime=dwelltime;
out.sw=1/dwelltime;

%FILLING IN THE FLAGS
out.flags=in.flags;
out.flags.writtentostruct=1;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
%PROCM Fixed-cell mapping to execute arbitray Matlab function on objects
%      Low-level routine
%
%    [OUT1,OUT2, ..] = A*PROCM(COMMAND,PAR1,PAR2,....)
%    [OUT1,OUT2, ..] = {A,B}*PROCM(COMMAND,PAR1,PAR2,....)
%
% INPUT
%    A               Dataset, datafile, cell array or double
%    B               Dataset, datafile, cell array or double
%    COMMAND         String with function name or PRTools nmapping
%    PAR1,...        Optional parameters for COMMAND
%
% OUTPUT
%    OUT1, ...       Cell arrays with results from COMMAND
%
% DESCRIPTION
% This routine executes a given COMMAND over a set of objects. These can
% be supplied by a dataset, a datafile, a cell array or a double array.
% Rows are interpreted as objects. For each object N stored in A the
% routine COMMAND is processed by
%
%    [OUT1{N}, ...] = COMMAND(OBJECT_N,PAR1,PAR2,....)
%
% For diadic operations the two inputs should be combined in a cell: {A,B}.
% They are processed by
%
%    [OUT1{N}, ...] = COMMAND(OBJECT_A_N,OBJECT_B_N,PAR1,PAR2,....)
%
% This is a low-level routine that feeds the unpacked objects in datafiles
% and dataset (doubles) to COMMAND. Use FILTM for a higher level of
% processing. Use MAPM to handle entire matrices, datasets or datafiles
% instead of processing object by object.
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, IM2OBJ, FEAT2OBJ, DATA2IM, FILTM, FILTM, MAPM

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function varargout = procm(varargin)

  if nargin == 0
    error('No command found')
  end
  varargin = shiftargin(varargin,'char');
  if isempty(varargin{1})
    w = prmapping(mfilename,'fixed_cell',varargin(2:end));
    varargout{1} = setname(w,'DataProc');
    return
  end

  % now we have: procm(data,command,pars)
  % data might be dataset, datafile, cell array, double, structure array

  nout = max(1,nargout);
  varargout = cell(1,nout); % output list
  argout = cell(1,nout);    % temp cell array for output per object
  pars = cell(1,nargin-2);     % input pars
  [a,command,pars{:}] = deal(varargin{:}); % use nice names
%   if numel(pars) > 1 % allow argument list as cell array or not.
%    pars = {pars}
%   end
  if iscell(a) && all(size(a)==[1,2])      % diadic operation
    [a,b] = deal(a{:});
    diadic = true;
  else
    diadic = false;
  end
  m = size(a,1);
  for j=1:nout              % output space
    varargout{j} = cell(m,1);
  end
  t = sprintf(['Processing %i objects by ' command ': '],m);
  prwaitbar(m,t)
  for i=1:m                   % loop over objects
    prwaitbar(m,i,[t num2str(i)]);
    if diadic                 % call command
      [argout{:}] = feval(command,g(a,i),g(b,i),pars{:});
    else
      x = g(a,i);
      [argout{:}] = feval(command,g(a,i),pars{:});
    end
    for j=1:nout           % store results
      varargout{j}{i} = argout{j};
    end
  end
  prwaitbar(0);

return

function x = g(a,i)
  if iscell(a)              % get cell data
    x =a{i};
  elseif isa(a,'prdataset') % datasets and datafiles
    x = +a(i,:);            % get their data
  else                      % all other datatypes
    x = a(i,:);
  end
return


text:
function tests = LineFeatureTest
    tests = functiontests(localfunctions)
    clc
end

function setupOnce(tc)
    im = testpattern('squares', 256, 256, 128);
    im = irotate(im, -0.3);
    edges = icanny(im);

    tc.TestData.edges = edges;
end

function teardownOnce(tc)
    close all
end

function constructor_test(tc)

    % create Hough object
    h = Hough(tc.TestData.edges);
    tc.verifyClass(h, 'Hough');

    s = char(h);
    tc.verifyClass(s, 'char');

    h.display;
    h.show();
end

function simple_test(tc)
    % check accumulators are correct size
    im = [0 0 0; 0 1 0; 0 0 0];
    h = Hough(im, 'nbins', [6 5]);
    tc.verifyEqual( size(h.A), [5 6]);

    h = Hough(im, 'nbins', 4);
    tc.verifyEqual( size(h.A), [5 4]);

    % simple example with 2 edge points
    edge = [
        0 0 0
        0 1 1
        0 0 0];
    out = [
        0     0     0     0
        2     0     0     0
        0     2     0     0
        0     0     2     1
        0     0     0     1
     ];
    h = Hough(edge, 'nbins', 4)
    tc.verifyEqual(h.A, out);

     % test vote weighting
     h = Hough(edge*2, 'nbins', 4);
     tc.verifyEqual( h.A, out*2);

     h = Hough(edge*2, 'nbins', 4, 'equal');
     tc.verifyEqual( h.A, out);

     % test xy input mode
     h = Hough([2 2; 2 3]', 'nbins', 4, 'points', [3 3]);
     tc.verifyEqual( h.A, out);

     % test xy input mode with weights
     h = Hough([2 2 5; 2 3 5]', 'nbins', 4, 'points', [3 3]);
     tc.verifyEqual( h.A, out*5);

     % test point, rather than image, mode
%      h=Hough([2 3; 2 2], 'points', 'nbins', 4)
%      tc.verifyEqual( h.A, out);
end

function square_test(tc)
    edges = tc.TestData.edges;

    h = Hough(edges);

    % create LineFeature objects
    lines = h.lines();
    tc.verifyClass(lines, 'LineFeature');

    tc.verifyEqual(numel(lines), 11);
    tc.verifyTrue( isempty(lines.length) );
    tc.verifyEqual( length(lines.theta), 11);
    tc.verifyEqual( length(lines.rho), 11);
    tc.verifyEqual( length(lines.strength), 11);

    s = char(lines);
    tc.verifyClass(s, 'char');
    tc.verifyEqual( size(s,1), 11);

end

function show_test(tc)
    h = Hough(tc.TestData.edges);

    clf
    h.show();
    a = gca;
    tc.verifyNotEmpty(a.Children);
    tc.verifyNumElements(a.Children, 1);
    tc.verifyMatches(a.Children.Type, 'image')

end

function plot_test(tc)
    h = Hough(tc.TestData.edges);
    clf
    h.plot('b');

    a = gca;
    tc.verifyNotEmpty(a.Children);
    tc.verifyNumElements(a.Children, 11);
    tc.verifyTrue(all( strcmp({a.Children.Type}, 'line') ) )
end

function square_suppress_test(tc)
    edges = tc.TestData.edges;

    % create new Hough and lines with non-local maxima suppression
    h = Hough(edges, 'suppress', 5);
    tc.verifyClass(h, 'Hough');

    lines = h.lines();
    tc.verifyClass(lines, 'LineFeature');
    tc.verifyEqual(numel(lines), 4);

    s = char(lines);
    tc.verifyClass(s, 'char');
    tc.verifyEqual(size(s,1), 4);

    % test line support
    lines = lines.seglength(edges);
    tc.verifyClass(lines, 'LineFeature');
    tc.verifyEqual(numel(lines), 4);

    tc.verifyEqual( numel(lines), 4);
    tc.verifyEqual( length(lines.theta), 4);
    tc.verifyEqual( length(lines.rho), 4);
    tc.verifyEqual( length(lines.length), 4);
    tc.verifyEqual( length(lines.strength), 4);
end

author: petercorke

repo: machinevision-toolbox-matlab


text:
% P. Perona and W. T. Freeman, "A factorization approach to grouping",
% In H. Burkardt and B. Neumann, editors, Proc ECCV, pages 655-670, 1998.

% Yair Weiss, "Segmentation Using Eigenvectors: A Unifying View", UC
% Berkeley.

% Asad Ali
% GIK Institute of Engineering Sciences & Technology, Pakistan
% Email: asad_82@yahoo.com

% CONCEPT: Clustering algorithm based on thresholding of the first eigenvector of
% the affinity matrix
clear all;
close all;

THRESHOLD = 0.28; % change the threshold depending on the data.
% generate the data
data = GenerateData(1);
% break the block matrix structure
data = data(3:14,:);
figure,plot(data(:,1), data(:,2),'r+'),title('Original Data Points'); grid on;shg

affinity = CalculateAffinity(data);
figure,imshow(affinity,[]),title('Affinity Matrix')

[eigVectors,eigValues] = eig(affinity);

% first eigenvector is the one which has the highest eigenvalue
sz = size(eigVectors);
firstEig = eigVectors(:,sz(1,2));

% plot the eigen vector corresponding to the largest eigen value
[xx1,yy1,val1] = find(firstEig > THRESHOLD);
[xx2,yy2,val2] = find(firstEig <= THRESHOLD);
figure,
hold on;
plot(data(xx1,1),data(xx1,2),'g*')
plot(data(xx2,1),data(xx2,2),'b*')
hold off;
title('Clustering Results after thresholding the First Eigenvector');
grid on;shg


text:
classdef HaloArriveDepartManifoldEnum < matlab.mixin.SetGet
    %HaloArriveDepartManifoldEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        Arrive('Arrival Manifold');
        Depart('Departure Manifold');
    end

    properties
        name char = '';
    end

    methods
        function obj = HaloArriveDepartManifoldEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function [listBoxStr, enums] = getListBoxStr()
            m = enumeration('HaloArriveDepartManifoldEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
            enums = m;
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('HaloArriveDepartManifoldEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('HaloArriveDepartManifoldEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [x,fval,exitflag,info] = multiSolveOpti(optObj,user_x0,ndivs,penalty,solveAll)
%MULTISOLVESOPTI an OPTI object
%
%   Called By opti multisolve

%   Very basic and naive implementation of a multi-start solver, using
%   pretty much an exhaustive search. However it has been very useful on
%   problems where an initial guess is hard to find!

%   Copyright (C) 2013 Jonathan Currie (I2C2)

%Number of 'best points' to keep during search
nbestpts = 10;
%Number of points to use for phase 2
nphase2pts = ceil(nbestpts/2);
%Number of divisions for phase 2 grid
ndivs2 = 2; %over written below if user specifies a vector
%Number of points to use for 'around x0' checking
nx0pts = 5;
%'around x0' radius factor
nx0rad = 0.05;

%If penalty is empty, set default
if(isempty(penalty)), penalty = 1e4; end

%Get display level
verb = dispLevel(optObj.opts.display);
%Build a new OPTI object for use here, disabling output as required (note sense = 1 to prevent double inversion)
Opt = opti(optObj,'sense',1,'options',optiset(optObj.opts,'display','off','warnings','none','derivCheck','off'));

%Allocate input args
prob = Opt.prob;
opts = Opt.opts;
%Initialization
where = '';
x_best = []; fval_best = inf; ef_best = []; info_best = []; i_best = 1;
titer = 0; tfeval = 0;

%If divs is empty, attempt to solve for a version that results in approx 100 sols, min 2 divisions
ndec = length(prob.lb);
if(isempty(ndivs))
    if(ndec==1)
        ndivs = 25;
    else
        ndivs = max(ceil(10^(2/ndec)),2);
    end
elseif(length(ndivs) == 2)
    ndivs2 = ndivs(2);
    if(ndivs2 < 2), nphase2pts = 0; end
    ndivs = ndivs(1);
end
%Determine number of points we are going to check
total_pts = ndivs^ndec;
%For search problems, limit nop based on total_pts
nbestpts = min(nbestpts,total_pts);

%Check division input
if(ndivs(1) == 1), error('You must have more than one division to use this function'); end
if(ndivs(1) > 1e4), error('The maximum number of divisions is 1e4'); end

%Check we have finite bounds
if(isempty(prob.lb) || any(isinf(prob.lb))), error('Multisolve currently requires finite lower bounds on all variables'); end
if(isempty(prob.ub) || any(isinf(prob.ub))), error('Multisolve currently requires finite upper bounds on all variables'); end

%Start Timer
t = tic;
%Build bound checking vectors
cvec = cell(ndec,1);
for i = 1:ndec
    cvec{i} = linspace(prob.lb(i),prob.ub(i),ndivs+1);
end

%Display header
if(verb)
    fprintf('------------------------------------------------------\n');
    fprintf('OPTI Multi-Start Solver using %s (beta)\n',upper(opts.solver));
end

%Construct bound selection matrix for phase 1
lbnd = zeros(total_pts,ndec); lbnd2 = [];
ubnd = zeros(total_pts,ndec); ubnd2 = [];
col = ones(1,ndec+1);
for i = 1:total_pts
    for j = 1:ndec
        lbnd(i,j) = cvec{j}(col(j));
        ubnd(i,j) = cvec{j}(col(j)+1);
    end
    col(1) = col(1) + 1;
    ind = col > ndivs;
    while(any(ind))
        idx = find(ind);
        col(idx) = 1;
        col(idx+1) = col(idx+1) + 1;
        ind = col > ndivs;
    end
end

if(verb)
    if(~solveAll)
        fprintf('\n Phase 1: Exhaustively Searching %d bounded regions (%d divisions per variable):\n',total_pts,ndivs);
    else
        fprintf('\n Exhaustively Solving %d bounded regions (%d divisions per variable):\n',total_pts,ndivs);
    end
end

%Begin searching bounded regions
x_search = cell(nbestpts,1); x_search2 = [];
lb_search = cell(nbestpts,1); ub_search = cell(nbestpts,1);
f_search = Inf(nbestpts,1); f_search2 = [];
for i = 1:total_pts
    %Assign Bounds for this iter
    lb = lbnd(i,:)'; ub = ubnd(i,:)';
    %Make x0
    x0 = (ub-lb)./2 + lb;

    %Search Every Point
    if(~solveAll)
        [x_search,f_search,lb_search,ub_search] = searchPoint(prob,x0,lb,ub,penalty,verb,i,total_pts,x_search,f_search,lb_search,ub_search);
    %Else Solve Every Point
    else
        [x_best,fval_best,ef_best,info_best,i_best,where,titer,tfeval] = solvePoint(Opt,lb,ub,x0,verb,i,total_pts,titer,tfeval,x_best,fval_best,ef_best,info_best,i_best,where,'Exhaustive Solve');
    end
end

%If we are not solving every problem, now solve all best points within their bounds, save best solution then enter phase 2
if(~solveAll)
    if(verb), fprintf('\n Solving Problem in the Best %d Phase 1 Regions:\n',nbestpts); end
    for i = 1:length(x_search)
        %Solve the problem within the search region
        [x_best,fval_best,ef_best,info_best,i_best,where,titer,tfeval] = solvePoint(Opt,lb_search{i},ub_search{i},x_search{i},verb,i,nbestpts,titer,tfeval,x_best,fval_best,ef_best,info_best,i_best,where,'Phase 1');
    end

    %Enter Phase 2, take noPhase2 best boxes and divide into divs2 divisions
    noPhase2 = min(nphase2pts,length(lb_search));
    if(noPhase2)
        total_pts2 = ndivs2^ndec;
        if(verb)
            fprintf('\n Phase 2: Exhaustively Searching %d Sub-Regions Within the %d Best Phase 1 Regions using %d Divisions...\n',total_pts2*noPhase2,noPhase2,ndivs2);
        end
        cvec2 = cell(ndec,noPhase2);
        for i = 1:noPhase2
            for j = 1:ndec
                cvec2{j,i} = linspace(lb_search{i}(j),ub_search{i}(j),ndivs2+1);
            end
        end
        lbnd2 = zeros(total_pts2*noPhase2,ndec);
        ubnd2 = zeros(total_pts2*noPhase2,ndec);
        for p = 1:noPhase2
            col = ones(1,ndec+1);
            for i = 1:total_pts2
                for j = 1:ndec
                    lbnd2(i + total_pts2*(p-1),j) = cvec2{j,p}(col(j));
                    ubnd2(i + total_pts2*(p-1),j) = cvec2{j,p}(col(j)+1);
                end
                col(1) = col(1) + 1;
                ind = col > ndivs2;
                while(any(ind))
                    idx = find(ind);
                    col(idx) = 1;
                    col(idx+1) = col(idx+1) + 1;
                    ind = col > ndivs2;
                end
            end
        end
        %Begin searching bounded regions
        x_search2 = cell(nbestpts,1); lb_search2 = cell(nbestpts,1); ub_search2 = cell(nbestpts,1); f_search2 = Inf(nbestpts,1);
        titer = 0; tfeval = 0; t = tic;
        for i = 1:total_pts2*noPhase2
            %Assign Bounds for this iter
            lb = lbnd2(i,:)'; ub = ubnd2(i,:)';
            %Make x0
            x0 = (ub-lb)./2 + lb;
            %Search Point
            [x_search2,f_search2,lb_search2,ub_search2] = searchPoint(prob,x0,lb,ub,penalty,verb,i,total_pts2*noPhase2,x_search2,f_search2,lb_search2,ub_search2);
        end

        if(verb), fprintf('\n Solving Problem in the Best %d Phase 2 Regions:\n',nbestpts); end
        for i = 1:length(x_search2)
            %Solve the problem within the search region
            [x_best,fval_best,ef_best,info_best,i_best,where,titer,tfeval] = solvePoint(Opt,lb_search2{i},ub_search2{i},x_search2{i},verb,i,length(x_search2),titer,tfeval,x_best,fval_best,ef_best,info_best,i_best,where,'Phase 2');
        end

        %Solve Problem within min and max bounds of best phase 2 divisions
        llb = Inf(ndec,1); gub = -Inf(ndec,1);
        for i = 1:length(lb_search2)
            for j = 1:ndec
                llb(j) = min([llb(j) lb_search2{i}(j) ub_search2{i}(j)]);
                gub(j) = max([gub(j) lb_search2{i}(j) ub_search2{i}(j)]);
            end
        end
        if(verb), fprintf('\n Solving Problem within Bounds of the Best Phase 2 Regions\n'); end
        [x_best,fval_best,ef_best,info_best,i_best,where,titer,tfeval] = solvePoint(Opt,llb,gub,(gub-llb)./2 + llb,verb,1,1,titer,tfeval,x_best,fval_best,ef_best,info_best,i_best,where,'Phase 2 Bounds');
    end
end

%Solve from user supplied start point
if(~isempty(user_x0))
    if(verb), fprintf('\n Phase 3: Searching Problem around User Supplied Start Point:\n'); end
    %Create search cells
    x_searchx0 = cell(nx0pts+1,1);
    f_searchx0 = Inf(nx0pts+1,1);
    %Variance of search points about user x0
    var = (prob.ub-prob.lb).*nx0rad;
    %Search user point
    [x_searchx0,f_searchx0] = searchPoint(prob,user_x0,[],[],penalty,verb,1,nx0pts+1,x_searchx0,f_searchx0);
    %Search random points
    for i = 2:nx0pts+1
        %Create random search points within bounds
        x0 = randnbnd(user_x0,var,prob.lb,prob.ub);
        %Search Point
        [x_searchx0,f_searchx0] = searchPoint(prob,x0,[],[],penalty,verb,i,nx0pts+1,x_searchx0,f_searchx0);
    end
    if(verb), fprintf('\n Solving Problem from the Best 3 Phase 3 Points:\n'); end
    %Solve Problem
    [x_best,fval_best,ef_best,info_best,i_best,where,titer,tfeval] = solvePoint(Opt,prob.lb,prob.ub,user_x0,verb,1,3,titer,tfeval,x_best,fval_best,ef_best,info_best,i_best,where,'User x0');
    for i = 1:2
        [x_best,fval_best,ef_best,info_best,i_best,where,titer,tfeval] = solvePoint(Opt,prob.lb,prob.ub,user_x0,verb,i+1,3,titer,tfeval,x_best,fval_best,ef_best,info_best,i_best,where,'Point Around User x0');
    end
else
    x_searchx0 = [];
    f_searchx0 = [];
end

%Check we have a solution so far
if(isempty(x_best))
    error('OPTI Multi-Solve did not find a feasible/successful solution.\n\nPlease increase the number of search points via the third argument to multisolve()%s','.');
end

%Check best solution in original problem
if(verb)
    if(~solveAll)
        fprintf('\nSolving Original Problem from the Best Point found...\n');
    else
        fprintf('\nRe-solving Original Problem from the Best Point found...\n');
    end
end
%Now Resolve from the best point, given original problem
[x,fval,exitflag,info] = solveOpti(optObj,x_best);
%Sum Iterations + FuncEvals
titer = incField(info,'Iterations',titer);
tfeval = incField(info,'FuncEvals',tfeval);
%Print Results
if(verb)
    status = getStatusMsg(exitflag);
    if(~isfield(info,'Iterations'))
        fprintf('Final Solution: Fval %12.5g, Status %s\n',fval,status);
    else
        fprintf('Final Solution: Fval %12.5g, Iter %4d, Status %s\n',fval,info.Iterations,status);
    end
end
%Check if final solution is actually better than best solution
if(fval_best < fval)
    x = x_best;
    fval = fval_best;
    exitflag = ef_best;
    info = info_best;
    if(verb)
        fprintf('\nFinal Solution is worse than best found, returning best found (%s, Run %d)\n',where,i_best);
    end
end

%Modify Info Status
info.Time = toc(t);
if(titer), info.Iterations = titer; end
if(tfeval), info.FuncEvals = tfeval; end

%Save Search Area
optObj.prob.multi.lbnd = lbnd;
optObj.prob.multi.ubnd = ubnd;
optObj.prob.multi.lbnd2 = lbnd2;
optObj.prob.multi.ubnd2 = ubnd2;
optObj.prob.multi.x_search = x_search;
optObj.prob.multi.f_search = f_search;
optObj.prob.multi.x_search2 = x_search2;
optObj.prob.multi.f_search2 = f_search2;
optObj.prob.multi.x_searchx0 = x_searchx0;
optObj.prob.multi.f_searchx0 = f_searchx0;

if(verb)
    fprintf('------------------------------------------------------\n');
end

%Search a Point, print as required
function [x_search,f_search,lb_search,ub_search] = searchPoint(prob,x0,lb,ub,penalty,verb,i,total_pts,x_search,f_search,lb_search,ub_search)
if(nargin < 12), ub_search = []; end
if(nargin < 11), lb_search = []; end
%Evaluate point
[fpen,f,c] = evalPoint(prob,x0,penalty);
%Decide whether we save this point
idx = fpen < f_search;
bstr = 1;
if(any(idx))
    %Find first index we are better than, insert, and shift remainder
    fidx = find(idx); fidx = fidx(1);
    %Insert
    f_search = [f_search(1:fidx-1);fpen;f_search(fidx:end-1)];
    x_search = [x_search(1:fidx-1);x0;x_search(fidx:end-1)];
    if(~isempty(lb_search))
        lb_search = [lb_search(1:fidx-1);lb;lb_search(fidx:end-1)];
        ub_search = [ub_search(1:fidx-1);ub;ub_search(fidx:end-1)];
    end
    %Indicate new solution found
    bstr = 2;
end
if(verb)
    if(~isempty(c))
        fprintf(bstr,'Point %3d of %3d: FPenalty %12.5g, Fval %12.5g, ConViol %12.5g\n',i,total_pts,fpen,f,c);
    else
        fprintf(bstr,'Point %3d of %3d: Fval %12.5g\n',i,total_pts,f);
    end
end

%Solve a Point, print as required
function [x_best,fval_best,ef_best,info_best,i_best,where,titer,tfeval] = solvePoint(Opt,lb,ub,x0,verb,i,total_pts,titer,tfeval,x_best,fval_best,ef_best,info_best,i_best,where,phase)
%Save new bounds into OPTI object
Opt.prob.lb = lb; Opt.prob.ub = ub;
Opt.nlprob.lb = lb; Opt.nlprob.ub = ub;
Opt.nlprob.options.lb = lb; Opt.nlprob.options.ub = ub;
bstr = 1;
%Solve Problem
try
    %Solve Problem
    [xi,fi,ei,ii] = solveOpti(Opt,x0);
    %Save results into object
    Opt.sol = xi; Opt.obj = fi; Opt.ef = ei; Opt.info = ii;
    %If best, save it
    if(fi < fval_best)
        %Check point is actually feasible
        if(checkOptiSol(Opt,1e-4)) %slightly relaxed tol
            x_best = xi; fval_best = fi; ef_best = ei; info_best = ii; i_best = i; bstr = 2; where = phase;
        end
    end
    %Sum Iterations + FuncEvals
    titer = incField(ii,'Iterations',titer);
    tfeval = incField(ii,'FuncEvals',tfeval);
catch ME
    fi = NaN; ei = NaN; ii = NaN;
    optiwarn('OPTI:SolverError',ME.message);
end
%Display Iteration if requested
if(verb)
    status = getStatusMsg(ei);
    if(~isfield(ii,'Iterations'))
        fprintf(bstr,'Run %3d of %3d: Fval %12.5g, Status %s\n',i,total_pts,fi,status);
    else
        fprintf(bstr,'Run %3d of %3d: Fval %12.5g, Iter %4d, Status %s\n',i,total_pts,fi,ii.Iterations,status);
    end
end

function [fpen,f,c] = evalPoint(prob,x0,penalty)
%Evaluate point
f = prob.objective(x0);
if(~isempty(prob.constraints))
    c = sum(prob.constraints(x0));
    fpen = f + c*penalty;
else
    c = [];
    fpen = f;
end

function status = getStatusMsg(ef)
switch(ef)
    case 1
        status = 'OK';
    case 0
        status = 'Exceeded It/Time';
    case -1
        status = 'Infeasible';
    otherwise
        status = 'Error';
end

function val = incField(strct,field,val)
if(isfield(strct,field) && isnumeric(strct.(field)))
    val = val + strct.(field);
end

function x0 = randnbnd(mean,var,lb,ub)
ndec = length(lb);
%Clamp x0 outside bounds
il = mean < lb; mean(il) = lb(il);
iu = mean > ub; mean(iu) = ub(iu);
%Check if avg sitting on a bound, if so try move by var in
for i = 1:ndec
    if(mean(i) == lb(i))
        mean(i) = min(mean(i) + var(i),ub(i));
    elseif(mean(i) == ub(i))
        mean(i) = max(mean(i) - var(i),lb(i));
    end
end
%Now generate random numbers until all fall within bounds, or 100 iter
idx = true(ndec,1); iter = 1; x0 = zeros(ndec,1);
while(any(idx) && iter < 100)
    xnew = mean(idx) + var(idx).*randn(sum(idx),1);
    x0(idx) = xnew;
    idx = x0 < lb | x0 > ub;
    iter = iter + 1;
end
if(iter==100)
    %re clamp
    il = x0 < lb; x0(il) = lb(il);
    iu = x0 > ub; x0(iu) = ub(iu);
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function cg_test01 ( n )

%*****************************************************************************80
%
%% CG_TEST01 tests R8GE_CG for a full storage matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the size of the system.
%
  if ( nargin < 1 )
    n = 10;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CG_TEST01\n' );
  fprintf ( 1, '  Test R8GE_CG, applying CG to a full storage matrix.\n' );

  seed = 123456789;
%
%  Choose a random positive definite symmetric matrix A.
%
  [ a, seed ] = pds_random ( n, seed );
%
%  Choose a random solution.
%
  [ x1, seed ] = r8vec_uniform_01 ( n, seed );
%
%  Compute the corresponding right hand side.
%
  b = r8ge_mv ( n, n, a, x1 );
%
%  Call the CG routine.
%
  x2 = ones ( n, 1 );
  x2 = r8ge_cg ( n, a, b, x2 );
%
%  Compute the residual.
%
  r = r8ge_resid ( n, n, a, x2, b );
  r_norm = r8vec_norm ( n, r );
%
%  Compute the error.
%
  e_norm = r8vec_diff_norm ( n, x1, x2 );
%
%  Report.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of variables N = %d\n', n );
  fprintf ( 1, '  Norm of residual ||Ax-b|| = %g\n', r_norm );
  fprintf ( 1, '  Norm of error ||x1-x2|| = %g\n', e_norm );

  return
end


text:
function r = intrand(p, m, n);
%INTRAND Random matrices from uniform integer distribution.
%
%   R = INTRAND(P) returns a matrix of random numbers chosen
%   from the uniform integer distribution with parameter P. The
%   minimum and maximum output integer numbers are P(1) and P(2).
%   R = INTRAND(P,M,N) returns an M by N matrix.
%

% Copyright (c) 2000 Aki Vehtari

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

if nargin < 1 | any(size([p(:)])~=[2 1])
  error('Incorrect range assignment in INTRAND');
end
if p(1)<p(2)
  rmin=ceil(p(1));
  rrange=floor(p(2))-rmin+1;
else
  rmin=ceil(p(2));
  rrange=floor(p(1))-rmin+1;
end

if nargin < 2
  r=floor(rand.*rrange+rmin);
else
  if nargin < 3
    n=m;
  end
  r=floor(rand([m n]).*rrange+rmin);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function alpha = backtracking_line_search_prox(problem, w, alpha, rho)
% This is a proximal backtracking algorithm.
%
% Inputs:
%       problem     function (cost/grad/hess)
%       w      current point
%       alpha       current stepsize
%       rho         shrink constant (<1)
% Output:
%       alpha       new stepsize
%
% References:
%
%
% This file is part of SGDLibrary.
%
% Created by H.Kasai on Nov. 19, 2018.
%
%   F(x) := f(x) + g(x);

    w0 = w;

    %% f0
    % calculate f0
%     F0 = problem.calculate_cost(w0);
%     reg0 = problem.calculate_reg(w0);
%     g0 = problem.lambda * reg0;
%     f0 = F0 - g0;
    f0 = problem.differentiable_cost(w0);

    % calculate g0
    grad0 = problem.full_grad(w0);

    % w = w - alpha * grads0
    w_out = w0 - alpha * grad0;

    % prox
    w_out = problem.prox(w_out, alpha);

    %% f1
    fk = problem.differentiable_cost(w_out);

    diff = w_out - w0;

    while fk > f0 + grad0'*diff + 1/(2*alpha) * (diff'*diff)
        alpha = rho * alpha;

        % w = w - alpha * grads0
        w_out = w0 - alpha * grad0;

        % prox
        w_out = problem.prox(w_out, alpha);

        %% fk
        fk = problem.differentiable_cost(w_out);
        diff = w_out - w0;

    end

%     while fk >= f0 + grad_f(x)'*z + (0.5/step)*norm(z,2)^2
%         lambda = rho * step;
%         z = prox(x - lambda*grad_f(x),lambda*g);
%     end

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function plot_recall_all_datasets()
  % Plot the average recall of PASCAL, ImageNet, and COCO for all methods.
  %
  % This function requires the proposals to already be saved to disk. It will
  % compute a matching between ground truth and proposals (if the result is not
  % yet found on disk) and then plot recall curves. The plots are saved to
  % figures/.

  methods = get_method_configs();
  methods([14 16 19:25]) = [];

  pascal = {methods.best_voc07_candidates_file};
  imagenet = {methods.best_imagenet_candidates_file};
  coco = {methods.best_coco14_candidates_file};

  valid_methods = check_for_data(cat(2, pascal(:), imagenet(:), coco(:)));
  if ~any(valid_methods)
    fprintf('no methods have the recall matchings computed, run the plot scripts for each of the datasets, to get the required data\n');
    return
  end
  methods = methods(valid_methods);
  pascal = {methods.best_voc07_candidates_file};
  imagenet = {methods.best_imagenet_candidates_file};
  coco = {methods.best_coco14_candidates_file};

%   plot_legend(methods);
%   printpdf('figures/imagenet_recall_legend.pdf');

  fprintf('\n\n');
  for i = 1:numel(methods)
    fprintf('Valid data for plots: %s\n', methods(i).name);
  end

  dataset_labels = {'PASCAL 2007', 'ImageNet 2013', 'COCO 2014'};
  AR_plots(methods, cat(2, pascal(:), imagenet(:), coco(:)), dataset_labels);
end

function AR_plots(methods, iou_files, dataset_labels)
  num_candidates = 1000;

  n = size(iou_files, 1);
  num_datasets = size(iou_files, 2);

  for i = 1:n
    ARs = zeros(1, num_datasets);
    for dataset_i = 1:num_datasets, iou_file = iou_files{i,dataset_i};
      try
        data = load(iou_file);

        num_experiments = numel(data.best_candidates);
        exp_idx = find([data.best_candidates.candidates_threshold] == num_candidates);
        assert(numel(exp_idx) == 1);
        experiment = data.best_candidates(exp_idx);
        [~, ~, ARs(dataset_i)] = compute_average_recall(experiment.best_candidates.iou);
      catch
        ARs(dataset_i) = nan;
      end
    end

    line_style = '-';
    if methods(i).is_baseline
      line_style = '--';
    end
    if ~isempty(methods(i).line_style)
      line_style = methods(i).line_style;
    end
    plot(1:num_datasets, ARs, '.', 'MarkerSize', 20, 'Color', methods(i).color, 'LineWidth', 1.5, 'LineStyle', line_style);
    hold on;
  end
  xlim([0.5, num_datasets + 0.5]);
  ylim([0.1, 0.6]);
  ylabel('average recall');
  hold off; grid on;
  set(gca,'XTick',(1:num_datasets),'XTickLabel',dataset_labels);
  set(gca,'XTickLabelRotation',25);
  set(gca, 'box','off');
  set(gca,'YTick', 0.1:0.1:0.6);

  hei = 10;
  wid = 8;
  set(gcf, 'Units','centimeters', 'Position',[0 0 wid hei]);
  set(gcf, 'PaperPositionMode','auto');
  printpdf('figures/average_recall_all_datasets.pdf');
end

function valid = check_for_data(candidate_files)
  valid = true(size(candidate_files));
  for i = 1:size(candidate_files,1)
    for j = 1:size(candidate_files,2)
      try
        ld = load(candidate_files{i,j});
      catch
        fprintf('file missing: %s\n', candidate_files{i,j});
        valid(i,j) = false;
      end
    end
  end
  valid = all(valid,2);
end

author: hosang

repo: detection-proposals


text:
function [GS, pp, info] = slgmm(X, varargin)
%SLGMM Learns Gaussian Mixture model from samples
%
% $ Syntax $
%   - GS = slgmm(X, ...)
%   - [GS, pp] = slgmm(X, ...)
%   - [GS, pp, info] = slgmm(X, ...)
%
% $ Arguments $
%   - GS:       The Gaussian model struct with mixture weights
%   - pp:       the posteriori of samples to component models
%   - info:     the information of learning process
%   - X:        the sample matrix
%
% $ Description $
%   - GS = slgmm(X, ...) learns Gaussian mixture model from samples
%     in matrix X. The following properties can be specified:
%     \*
%     \t    Table. Properties of Learning Gausssian Mixture Model
%     \h    name         &      description                      \\
%          'method'      & The method using for GMM learning. Currently,
%                          there are only one method available: 'EM',
%                          default = 'EM'.
%          'K'           & The number of initial components in the mixture
%                          model (default = 3)                     \\
%          'update'      & The way of updating (default = 'pass'):
%                          1. 'pass'     Pass-wise update;
%                          2. 'comp'     Component-wise update     \\
%          'cyclecn'     & The ratio of components to be updated in each
%                          cycle for 'comp' update scheme. default = 1.
%          'maxiter'     & The maximum number of iterations
%                          (default = 100) \\
%          'tol'         & The maximum tolerance of posteriori error when
%                          the iteration terminates, (default = 1e-6) \\
%          'verbose'     & whether to display information while iteration,
%                          default = true                              \\
%          'annthres'    & The threshold of annealing in FJ algorithm,
%                          default = 0 (unit = average mixture weight) \\
%          'weights'     & The weights of the samples, default = [],
%                          indicating non-weighted. Weights should be given
%                          by a 1 x n row vector. \\
%          'varform'     & The form of variance: 'univar'|'diagvar'|'covar',
%                          default = 'covar';          \\
%          'sharevar'    & Whether the variance is shared by all models,
%                          default = false;            \\
%          'invparams'   & The parameters for computing inverse of variance
%                          default = {}                \\
%
% $ Remarks $
%   - In current version, the component-wise update scheme is only
%     supported when sharevar is false.
%
% $ History $
%   - Created by Dahua Lin, on Aug 28, 2006
%

%% parse and verify input arguments

if ~isnumeric(X) || ndims(X) ~= 2
    error('sltoolbox:invalidarg', ...
        'X should be a 2D numeric matrix');
end

opts.method = 'EM';
opts.K = 3;
opts.update = 'pass';
opts.cyclecn = 1;
opts.maxiter = 100;
opts.tol = 1e-6;
opts.verbose = true;
opts.annthres = 0;
opts.weights = [];
opts.varform = 'covar';
opts.sharevar = false;
opts.invparams = {};
opts = slparseprops(opts, varargin{:});

n = size(X, 2);

%% initialization by random clustering

if opts.K > n
    error('sltoolbox:sizoverflow', ...
        'The K is larger than the number of samples');
end

init_centerinds = randsample(n, opts.K);
init_centers = X(:, init_centerinds);
initc = slclassify_eucnn(init_centers, X);

%% perform learning based on FMM

estfunctor = {@gmm_estfunc, opts.varform, opts.sharevar, opts.invparams};
evalfunctor = {@gmm_evalfunc};

[GS, cw, pp, info] = slfmm(X, n, estfunctor, evalfunctor, ...
    'method', opts.method, ...
    'update', opts.update, ...
    'cyclecn', opts.cyclecn, ...
    'iter', {'maxiter', opts.maxiter}, ...
    'tol', opts.tol, ...
    'verbose', opts.verbose, ...
    'initc', initc, ...
    'annthres', opts.annthres, ...
    'weights', opts.weights, ...
    'estmode', 'innermul', ...
    'condpmode', 'log', ...
    'manifunc', @gmm_manifunc);

%% post actions

GS.mixweights = cw;

%% Core slot functions

function GS = gmm_estfunc(GS, X, n, W, selinds, varform, sharevar, invparams)

slignorevars(n);

if isempty(selinds)
    GS = slgaussest(X, 'weights', W, ...
        'varform', varform, 'sharevar', sharevar, ...
        'compinv', true, 'invparams', invparams);
else
    if sharevar
        error('sltoolbox:rterror', ...
            'The selective updating scheme is only supported when sharevar is false');
    end

    Wsel = W(selinds, :);
    GSu = slgaussest(X, 'weights', Wsel, ...
        'varform', varform, 'sharevar', sharevar, ...
        'compinv', true, 'invparams', invparams);

    nu = length(selinds);
    switch varform
        case 'univar'
            for idx = 1 : nu
                iu = selinds(idx);
                GS.means(:, iu) = GSu.means(:,idx);
                GS.vars(iu) = GSu.vars(idx);
                GS.invvars(iu) = GSu.invvars(idx);
            end
        case 'diagvar'
            for idx = 1 : nu
                iu = selinds(idx);
                GS.means(:, iu) = GSu.means(:, idx);
                GS.vars(:, iu) = GSu.vars(:, idx);
                GS.invvars(:, iu) = GSu.invvars(:, idx);
            end
        case 'covar'
            for idx = 1 : nu
                iu = selinds(idx);
                GS.means(:, iu) = GSu.means(:, idx);
                GS.covs(:,:, iu) = GSu.covs(:,:, idx);
                GS.invcovs(:,:, iu) = GSu.invcovs(:,:, idx);
            end
    end

end

function condp = gmm_evalfunc(GS, X, n)

slignorevars(n);
condp = slgausspdf(GS, X, 'output', 'log');

function GS = gmm_manifunc(GS0, op, varargin)

switch op
    case 'select'

        selinds = varargin{1};

        tyi = slgausstype(GS0);

        GS.dim = GS0.dim;
        GS.nmodels = length(selinds);
        GS.means = GS0.means(:, selinds);

        switch tyi.varform
            case {'univar', 'diagvar'}
                if tyi.sharevar
                    GS.vars = GS0.vars;
                    GS.invvars = GS0.invvars;
                else
                    GS.vars = GS0.vars(:, selinds);
                    GS.invvars = GS0.invvars(:, selinds);
                end
            case 'covar'
                if tyi.sharevar
                    GS.covs = GS0.covs;
                    GS.invcovs = GS0.invcovs;
                else
                    GS.covs = GS0.covs(:,:,selinds);
                    GS.invcovs = GS0.invcovs(:,:,selinds);
                end
        end

    otherwise
        error('sltoolbox:invalidarg', ...
            'Unsupported manipulation option for GMM: %s', op);
end

author: lmthang


text:
function   r = assortativity(CIJ,flag)
%ASSORTATIVITY      Assortativity coefficient
%
%   r = assortativity(CIJ,flag);
%
%   The assortativity coefficient is a correlation coefficient between the
%   degrees of all nodes on two opposite ends of a link. A positive
%   assortativity coefficient indicates that nodes tend to link to other
%   nodes with the same or similar degree.
%
%   Inputs:     CIJ,        binary directed/undirected connection matrix
%               flag,       1 = directed graph; 0 = non-directed graph
%
%   Outputs:    r,          assortativity
%
%   Notes: The function accepts weighted networks, but all connection
%   weights are ignored. The main diagonal should be empty.
%
%   Reference:  Newman (2002) Phys Rev Lett 89:208701.
%
%
%   Olaf Sporns, Indiana University, 2007/2008
%   Vassilis Tsiaras, University of Crete, 2009

if (flag==0)
    [deg] = degrees_und(CIJ);
    [i,j] = find(triu(CIJ,1)>0);
    K = length(i);
    for k=1:K
        degi(k) = deg(i(k));
        degj(k) = deg(j(k));
    end;
end;
if (flag==1)
    [id,od,deg] = degrees_dir(CIJ);
    [i,j] = find(CIJ>0);
    K = length(i);
    for k=1:K
        degi(k) = deg(i(k));
        degj(k) = deg(j(k));
    end;
end;

% compute assortativity
r = (sum(degi.*degj)/K - (sum(0.5*(degi+degj))/K)^2)/(sum(0.5*(degi.^2+degj.^2))/K - (sum(0.5*(degi+degj))/K)^2);