text:
%% Machine Learning Online Class
%  Exercise 8 | Anomaly Detection and Collaborative Filtering
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     estimateGaussian.m
%     selectThreshold.m
%     cofiCostFunc.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% =============== Part 1: Loading movie ratings dataset ================
%  You will start by loading the movie ratings dataset to understand the
%  structure of the data.
%
fprintf('Loading movie ratings dataset.\n\n');

%  Load data
load ('ex8_movies.mat');

%  Y is a 1682x943 matrix, containing ratings (1-5) of 1682 movies on
%  943 users
%
%  R is a 1682x943 matrix, where R(i,j) = 1 if and only if user j gave a
%  rating to movie i

%  From the matrix, we can compute statistics like average rating.
fprintf('Average rating for movie 1 (Toy Story): %f / 5\n\n', ...
        mean(Y(1, R(1, :))));

%  We can "visualize" the ratings matrix by plotting it with imagesc
imagesc(Y);
ylabel('Movies');
xlabel('Users');

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

%% ============ Part 2: Collaborative Filtering Cost Function ===========
%  You will now implement the cost function for collaborative filtering.
%  To help you debug your cost function, we have included set of weights
%  that we trained on that. Specifically, you should complete the code in
%  cofiCostFunc.m to return J.

%  Load pre-trained weights (X, Theta, num_users, num_movies, num_features)
load ('ex8_movieParams.mat');

%  Reduce the data set size so that this runs faster
num_users = 4; num_movies = 5; num_features = 3;
X = X(1:num_movies, 1:num_features);
Theta = Theta(1:num_users, 1:num_features);
Y = Y(1:num_movies, 1:num_users);
R = R(1:num_movies, 1:num_users);

%  Evaluate cost function
J = cofiCostFunc([X(:) ; Theta(:)], Y, R, num_users, num_movies, ...
               num_features, 0);

fprintf(['Cost at loaded parameters: %f '...
         '\n(this value should be about 22.22)\n'], J);

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

%% ============== Part 3: Collaborative Filtering Gradient ==============
%  Once your cost function matches up with ours, you should now implement
%  the collaborative filtering gradient function. Specifically, you should
%  complete the code in cofiCostFunc.m to return the grad argument.
%
fprintf('\nChecking Gradients (without regularization) ... \n');

%  Check gradients by running checkNNGradients
checkCostFunction;

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

%% ========= Part 4: Collaborative Filtering Cost Regularization ========
%  Now, you should implement regularization for the cost function for
%  collaborative filtering. You can implement it by adding the cost of
%  regularization to the original cost computation.
%

%  Evaluate cost function
J = cofiCostFunc([X(:) ; Theta(:)], Y, R, num_users, num_movies, ...
               num_features, 1.5);

fprintf(['Cost at loaded parameters (lambda = 1.5): %f '...
         '\n(this value should be about 31.34)\n'], J);

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

%% ======= Part 5: Collaborative Filtering Gradient Regularization ======
%  Once your cost matches up with ours, you should proceed to implement
%  regularization for the gradient.
%

%
fprintf('\nChecking Gradients (with regularization) ... \n');

%  Check gradients by running checkNNGradients
checkCostFunction(1.5);

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

%% ============== Part 6: Entering ratings for a new user ===============
%  Before we will train the collaborative filtering model, we will first
%  add ratings that correspond to a new user that we just observed. This
%  part of the code will also allow you to put in your own ratings for the
%  movies in our dataset!
%
movieList = loadMovieList();

%  Initialize my ratings
my_ratings = zeros(1682, 1);

% Check the file movie_idx.txt for id of each movie in our dataset
% For example, Toy Story (1995) has ID 1, so to rate it "4", you can set
my_ratings(1) = 4;

% Or suppose did not enjoy Silence of the Lambs (1991), you can set
my_ratings(98) = 2;

% We have selected a few movies we liked / did not like and the ratings we
% gave are as follows:
my_ratings(7) = 3;
my_ratings(12)= 5;
my_ratings(54) = 4;
my_ratings(64)= 5;
my_ratings(66)= 3;
my_ratings(69) = 5;
my_ratings(183) = 4;
my_ratings(226) = 5;
my_ratings(355)= 5;

fprintf('\n\nNew user ratings:\n');
for i = 1:length(my_ratings)
    if my_ratings(i) > 0
        fprintf('Rated %d for %s\n', my_ratings(i), ...
                 movieList{i});
    end
end

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

%% ================== Part 7: Learning Movie Ratings ====================
%  Now, you will train the collaborative filtering model on a movie rating
%  dataset of 1682 movies and 943 users
%

fprintf('\nTraining collaborative filtering...\n');

%  Load data
load('ex8_movies.mat');

%  Y is a 1682x943 matrix, containing ratings (1-5) of 1682 movies by
%  943 users
%
%  R is a 1682x943 matrix, where R(i,j) = 1 if and only if user j gave a
%  rating to movie i

%  Add our own ratings to the data matrix
Y = [my_ratings Y];
R = [(my_ratings ~= 0) R];

%  Normalize Ratings
[Ynorm, Ymean] = normalizeRatings(Y, R);

%  Useful Values
num_users = size(Y, 2);
num_movies = size(Y, 1);
num_features = 10;

% Set Initial Parameters (Theta, X)
X = randn(num_movies, num_features);
Theta = randn(num_users, num_features);

initial_parameters = [X(:); Theta(:)];

% Set options for fmincg
options = optimset('GradObj', 'on', 'MaxIter', 100);

% Set Regularization
lambda = 10;
theta = fmincg (@(t)(cofiCostFunc(t, Ynorm, R, num_users, num_movies, ...
                                num_features, lambda)), ...
                initial_parameters, options);

% Unfold the returned theta back into U and W
X = reshape(theta(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(theta(num_movies*num_features+1:end), ...
                num_users, num_features);

fprintf('Recommender system learning completed.\n');

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

%% ================== Part 8: Recommendation for you ====================
%  After training the model, you can now make recommendations by computing
%  the predictions matrix.
%

p = X * Theta';
my_predictions = p(:,1) + Ymean;

movieList = loadMovieList();

[r, ix] = sort(my_predictions, 'descend');
fprintf('\nTop recommendations for you:\n');
for i=1:10
    j = ix(i);
    fprintf('Predicting rating %.1f for movie %s\n', my_predictions(j), ...
            movieList{j});
end

fprintf('\n\nOriginal ratings provided:\n');
for i = 1:length(my_ratings)
    if my_ratings(i) > 0
        fprintf('Rated %d for %s\n', my_ratings(i), ...
                 movieList{i});
    end
end

author: Ayatans

repo: Machine-Learning-homework


text:
function [L,S] = RPCA(X)
[n1,n2] = size(X);
mu = n1*n2/(4*sum(abs(X(:))));
lambda = 1/sqrt(max(n1,n2));
thresh = 1e-7*norm(X,'fro');

S = zeros(size(X));
Y = zeros(size(X));
count = 0;
while((norm(X-L-S,'fro')>thresh)&&(count<1000))
    L = SVT(X-S+(1/mu)*Y,1/mu);
    S = shrink(X-L+(1/mu)*Y,lambda/mu);
    Y = Y + mu*(X-L-S);
    count = count + 1
end

author: dynamicslab


text:
function ni = niftiRead(fileName, volumesToLoad)
% Matlab wrapper to call the mex readFileNifti
%
% *******************************************************
% NOTE:  As of 2017 Matlab has a "niftiread" function in its
%        distribution.
%
% NOTE: JAN 2020 GLU & JW: this version uses an external tool called nii_util.m
%       taken from the toolset https://github.com/xiangruili/dicm2nii
%       It will be used only if it detects that the file to read is a Nifti-2,
%       otherwise it will use the old code
% *******************************************************
%
%   niftiImage = niftiRead(fileName,volumesToLoad)
%
% Reads a NIFTI image and populates a structure that should be the
% NIFTI 1 standard.  We expect that the filename has an extension of
% either .nii or .nii.gz
%
% If volumesToLoad is not included in the arguments, all the data
% are returned. If volumesToLoad is empty ([]) returns only the header
%
% Web Resources
%  web('http://nifti.nimh.nih.gov/nifti-1/','-browser')
%
% See also:  niftiCreate, niftiGetStruct
%
% Example:
%  niFile = fullfile(mrvDataRootPath,'mrQ','T1_lsqnabs_SEIR_cfm.nii.gz');
%  ni = niftiRead(niFile);
%
% Copyright, Vista Team Stanford, 2011

% This normally calls the mex file for your system
if ~exist('fileName','var') || isempty(fileName)
    % Return the default structure.  Equivalent to niftiCreate
    % ni = niftiRead;
    ni = readFileNifti;
elseif ischar(fileName) && exist(fileName,'file')
    % fileName is a string and the file exists.
    % For some reason, the volumeToLoad is not yet implemented.
    % We should just implement it here, by reading the whole
    % thing and only returning the relevant volumes.  I think
    % that is represented by the 4th dimension, but I should ask
    % someone who knows.
    if exist('volumesToLoad','var')
        % ni = niftiRead('foo.nii.gz',1:20);
        % We let readFileNifti complain about not implemented for
        % now.
        %
        % GLU EDIT: throw same error as the mex file
        % ni = readFileNifti(fileName,volumesToLoad);
        error('(volume selection) Not yet implemented!')
    else
        % Example:
        %     ni = niftiRead('foo.nii.gz');
        % Select the reader depending on nifti version
        % We need to read only the header first, to save time.
        % Matlab's niftiinfo will be a nice option in the future, it was
        % introduced in 2017a.
        hdr = nii_tool('hdr',fileName);
        if isfield(hdr,'version')
            if hdr.version==2
                % New call to the nii_tool wrapper
                ni = nii_tool_read_wrapper(fileName);
            else
                % Old call to the mex file, for backward compatibility
                ni = readFileNifti(fileName);
            end
        else
            % Old call to the mex file, for backward compatibility
            ni = readFileNifti(fileName);
        end

    end
else
    % Make sure the the file includes the .nii.gz extensions
    % Really, someone should
    [p,n,e] = fileparts(fileName);
    if isempty(e), fileNameExtended = [p,filesep,n,'.nii.gz'];
    else
        % The can be file.nii or file.nii.gz
        if strcmp(e,'.gz') || strcmp(e,'.nii')
            fileNameExtended = fileName;
        else
            warning('Unexpected file name extension %s\n',e);
        end
    end
    if exist(fileNameExtended,'file')
        % HERE WE MODIFY THE OLD VERSION
        hdr = nii_tool('hdr',fileName);
        if isfield(hdr,'version')
            if hdr.version==2
                % New call to the nii_tool wrapper
                ni = nii_tool_read_wrapper(fileName);
            else
                % Old call to the mex file, for backward compatibility
                ni = readFileNifti(fileName);
            end
        else
            % Old call to the mex file, for backward compatibility
            ni = readFileNifti(fileName);
        end
    else
        error('Cannot find the file %s\n',fileNameExtended);
    end
end

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Copyright 2010-2013 The MathWorks, Inc.
example_sel2html('wave_list.txt')

author: Sable

repo: mcbench-benchmarks


text:
%% HELP_ellipsoidFit_centered
% Below is a demonstration of the features of the |ellipsoidFit_centered| function

%% Syntax
% |[M,ellipStretch,R,MU]=ellipsoidFit_centered(X,MU);|

%% Description
% The |ellipsoidFit_centered| function fits an ellipsoid to data when the
% ellipsoid centre is known. If the centre is not provided the mean of the
% input point set will be assumed to be the centre.

%% Examples

%%
clear; close all; clc;

%%
% Plot settings

figColor='w';
figColorDef='white';
fontSize=11;

%% Example: Using |ellipsoidFit_centered| to fit an ellipsoid to a point cloud with known centre

%%
% Simulating an ellipsoid with known directions

% Ellipsoid axis stretch factors
ellipStretchTrue=[pi 2 0.5];
MU_true=[1 6 pi];

% Create ellipsoid patch data
[F,V,~]=geoSphere(3,1);
v=V(:,1);
FX=mean(v(F),2);
logicKeep=FX>0;
F=F(logicKeep,:);
indKeep=unique(F(:));
indFix=nan(size(V,1),1);
indFix(indKeep)=1:numel(indKeep);
V=V(indKeep,:);
F=indFix(F);
V=V.*ellipStretchTrue(ones(size(V,1),1),:);

%Create Euler angles to set directions
E=[0.25*pi 0.25*pi -0.25*pi];
[R_true,~]=euler2DCM(E); %The true directions for X, Y and Z axis
V=(R_true*V')'; %Rotate polyhedron

V=V+MU_true(ones(size(V,1),1),:); %Centre points around mean

%Add noise
n_std=0.2;  %Standard deviation
Vn=V+n_std.*randn(size(V));

%%
% This is the true axis system
R_true

%%
% These are the true stretch factors
ellipStretchTrue

%%

[M,ellipsStretchFit,R_fit,MU]=ellipsoidFit_centered(Vn,MU_true);

%%
% This is the fitted axis system. The system axes should be colinear with
% the true axes but can be oposite in direction.
R_fit=R_fit(1:3,1:3)

%%
% These are the fitted stretch factors
ellipsStretchFit

%%
% Building a fitted (clean) ellipsoid for visualization

%Create sphere
[F_fit,V_fit,~]=geoSphere(4,1);

%Transforming sphere to ellipsoid
V_fit_t=V_fit;
V_fit_t(:,end+1)=1;
V_fit_t=(M*V_fit_t')'; %Rotate polyhedron
V_fit=V_fit_t(:,1:end-1);

%%
% Visualizing results

cFigure; hold on;
title('The true (green) and fitted ellipsoid (red) and axis directions (solid, transparant respectively)','FontSize',fontSize);
plotV(Vn,'k.','MarkerSize',15);
gpatch(F,V,'gw','k',1);
gpatch(F_fit,V_fit,'rw','none',0.2);
quiverTriad(MU,R_fit,7,[],0.2);
quiverTriad(MU,R_true,7,[],1);
axisGeom;
camlight('headlight');
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef prtRv < prtAction
    % prtRv Base class for all prt random variables
    %
    %   This is an abstract class from which all prt random variables
    %   inherit. It can not be instantiated. prtRv contains the following
    %   properties:
    %
    %   name           - Name of the random variable.
    %   userData       - Structure for holding additional related to the
    %                    random variable.
    %   nDimensions    - Number of dimensions of the vector space
    %                    represented by the random variable.
    %
    %   The prtRv class has the following methods
    %
    %   plotPdf - Plot the pdf of the random variable
    %   plotCdf - Plot the cdf of the random variable
    %
    %   The prtRv class has the following methods, most of which are
    %   overloaded. If a method is not overloaded, it is because it is not
    %   possible to implement the functionality.
    %
    %   pdf - Output the pdf of the random variable evaluated at the points
    %         specified
    %
    %   logPdf - Output the log-pdf of the random variable evaluated at the
    %            points specified (for many distributions, this can be
    %            calculated more easily than simply log(pdf(R,X))
    %
    %   cdf - Output the cdf of the random variable evaluated at the
    %         points specified
    %
    %   draw - Draw samples from the random variable
    %
    %   mle - Perform maximum likelihood estimation of the objects parameters
    %         using the specified data
    %
    %   See also: prtRvMvn, prtRvGmm, prtRvMultinomial, prtRvUniform,
    %   prtRvUniformImproper, prtRvVq

    properties
        plotOptions  = prtRv.initializePlotOptions()
    end
    properties (Abstract = true, Hidden = true, Dependent = true)
        nDimensions % The number of dimensions
    end
    methods (Abstract, Hidden=true)
        [bool, reasonStr] = isValid(R)
    end
    methods
        % These functions are default "error" functions incase a
        % subclass has not specified these standard methods.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function vals = pdf(R,X) %#ok
            % PDF Output the pdf of the random variable evaluated at the points specified
            %
            % pdf = RV.pdf(X) returns the value of the pdf of the prtRv
            % object evaluated at X. X must be an N x nDims matrix, where N
            % is the number of locations to evaluate the pdf, and nDims is
            % the same as the number of dimensions, nDimensions, of the
            % prtRv object RV.
            missingMethodError(R,'pdf');
        end

        function vals = logPdf(R,X)
            % LOGPDF Output the log pdf of the random variable evaluated at the points specified
            %
            % logpdf = RV.logpdf(X) returns the logarithm of value of the
            % pdf of the prtRv object evaluated at X. X must be an N x
            % nDims matrix, where N is the number of locations to evaluate
            % the pdf, and nDims is the same as the number of dimensions,
            % nDimensions, of the prtRv object RV.
            vals = log(R.pdf(X));
        end

        function vals = cdf(R,X) %#ok
            % CDF Output the cdf of the random variable evaluated at the points specified
            %
            % cdf = RV.cdf(X) returns the value of the cdf of the prtRv
            % object evaluated at X. X must be an N x nDims matrix, where
            % N is the number of locations to evaluate the pdf, and nDims
            % is the same as the number of dimensions, nDimensions, of the
            % prtRv object RV.

            missingMethodError(R,'cdf');
        end

        function vals = draw(R,N) %#ok
            % DRAW  Draw random samples from the distribution described by the prtRv object
            %
            % VAL = RV.draw(N) generates N random samples drawn from the
            % distribution described by the prtRv object RV. VAL will be a
            % N x nDimensions vector, where nDimensions is the number of
            % dimensions of RV.

            missingMethodError(R,'draw');
        end

        function vals = mle(R,X) %#ok
            % MLE Compute the maximum likelihood estimate
            %
            % RV = RV.mle(X) computes the maximum likelihood estimate based
            % the data X. X should be nObservations x nDimensions.
            missingMethodError(R,'mle');
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % These functions are default plotting functions
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function varargout = plotPdf(R,varargin)
            % PLOTPDF  Plot the pdf of the prtRv
            %
            % R.plotPdf()
            % R.plotPdf([xMin yMin (zMin) xMax yMax (zMax)]);

            varargout = {};
            [isValid, reasonStr] = R.isValid;
            if R.isPlottable

                if nargin > 1 % Calculate appropriate limits from covariance
                    plotLims = varargin{1};
                else
                    plotLims = plotLimits(R);
                end
                [linGrid, gridSize] = prtPlotUtilGenerateGrid(plotLims(1:2:end), plotLims(2:2:end), R.plotOptions.nSamplesPerDim);

                imageHandle = prtPlotUtilPlotGriddedEvaledFunction(R.pdf(linGrid), linGrid, gridSize, R.plotOptions.colorMapFunction(R.plotOptions.nColorMapSamples),varargin{2:end});

                if nargout
                    varargout = {imageHandle};
                end
            else
                if isValid
                    error('prt:prtRv:plot','This RV object cannont be plotted because it has too many dimensions.')
                else
                    error('prt:prtRv:plot','This RV object cannot be plotted. This RV is not yet valid %s.',reasonStr);
                end
            end
        end
        function varargout = plotLogPdf(R,varargin)
            % plotLogPdf  Plot the pdf of the prtRv
            %
            % R.plotLogPdf()
            % R.plotLogPdf([xMin yMin (zMin) xMax yMax (zMax)]);

            varargout = {};
            [isValid, reasonStr] = R.isValid;
            if R.isPlottable

                if nargin > 1 % Calculate appropriate limits from covariance
                    plotLims = varargin{1};
                else
                    plotLims = plotLimits(R);
                end

                [linGrid, gridSize] = prtPlotUtilGenerateGrid(plotLims(1:2:end), plotLims(2:2:end), R.plotOptions.nSamplesPerDim);

                imageHandle = prtPlotUtilPlotGriddedEvaledFunction(R.logPdf(linGrid), linGrid, gridSize, R.plotOptions.colorMapFunction(R.plotOptions.nColorMapSamples),varargin{2:end});

                if nargout
                    varargout = {imageHandle};
                end
            else
                if isValid
                    error('prt:prtRv:plot','This RV object cannont be plotted because it has too many dimensions.')
                else
                    error('prt:prtRv:plot','This RV object cannot be plotted. This RV is not yet valid %s.',reasonStr);
                end
            end
        end

        function varargout = plotCdf(R,varargin)
            %PLOTCDF Plot the cdf of the prtRv
            % R.plotCdf()
            % R.plotCdf([xMin yMin (zMin) xMax yMax (zMax)]);

            varargout = {};
            [isValid, reasonStr] = R.isValid;
            if R.isPlottable

                if nargin > 1 % Calculate appropriate limits from covariance
                    plotLims = varargin{1};
                else
                    plotLims = plotLimits(R);
                end

                [linGrid,gridSize] = prtPlotUtilGenerateGrid(plotLims(1:2:end), plotLims(2:2:end), R.plotOptions.nSamplesPerDim);

                imageHandle = prtPlotUtilPlotGriddedEvaledFunction(R.cdf(linGrid), linGrid, gridSize, R.plotOptions.colorMapFunction(R.plotOptions.nColorMapSamples));

                if nargout
                    varargout = {imageHandle};
                end
            else
                if isValid
                    error('prt:prtRv:plot','This RV object cannont be plotted because it has too many dimensions for plotting.')
                else
                    error('prt:prtRv:plot','This RV object cannot be plotted. This RV is not yet valid %s.',reasonStr);
                end
            end
        end
    end

    methods (Hidden = true)
        % As a default we just use the monte carlo limits.
        % If a sub-class implements plotLimits ezPdfPlot or ezCdfPlot
        % should use that.
        function limits = plotLimits(R)
            limits = monteCarloPlotLimits(R);
        end

        function val = isPlottable(R)
            val = ~isempty(R.nDimensions) && R.nDimensions < 4 && R.isValid;
        end
    end

    methods (Access = protected, Hidden = true)
        function Obj = trainAction(Obj, DataSet)
            Obj = Obj.mle(DataSet);
        end

        function DataSet = runAction(Obj, DataSet)
            DataSet = DataSet.setObservations(Obj.logPdf(DataSet));
        end

        function xOut = runActionFast(Obj, xIn, ds) %#ok<INUSD>
            xOut = Obj.logPdf(xIn);
        end
    end

    methods (Access = 'private', Hidden = true)
        function missingMethodError(R,methodName) %#ok<MANU>
            error('The method %s is not defined for this prtRv object',methodName);
        end
    end
    methods (Access = 'private',Hidden = true,Static = true)
        function plotOptions = initializePlotOptions()
            if ~isdeployed
                plotOptions = prtOptionsGet('prtOptionsRvPlot');
            else
                plotOptions = prtOptions.prtOptionsRvPlot;
            end
        end
    end

    methods (Access = 'protected', Hidden = true)
        function X = dataInputParse(R,X) %#ok<MANU>
            % dataInputParse - Parse inputs for RVs that only require the
            % data. Since most RVs need only the observations (X). This
            % function allows RVs to operate on prtDataSetStandard()s and
            % a data matrix.

            if isnumeric(X) || islogical(X)
                % Quick exit from this ifelse so we don't call isa
                % which can be slow
            elseif isa(X,'prtDataSetStandard')
                X = X.getObservations();
            else
                error('prt:prtRv','Input to mle() must be a matrix of data or a prtDataSet.');
            end

        end

        function R = constructorInputParse(R,varargin)

            nIn = length(varargin);

            % Quick Exit for the zero input constructor
            if nIn == 0
                return
            elseif mod(nIn,2)
                error('prt:prtRv:constructorInputParse','Inputs must be supplied by as string/value pairs');
            end

            R = prtUtilAssignStringValuePairs(R,varargin{:});

        end

        function val = monteCarloCovariance(R,nSamples)
            % Calculates the sample covariance of drawn data
            if nargin < 2 || isempty(nSamples)
                nSamples = 1e3;
            end
            val = cov(draw(R,nSamples));
        end

        function val = monteCarloMean(R,nSamples)
            % Calculates the sample mean of drawn data
            if nargin < 2 || isempty(nSamples)
                nSamples = 1e3;
            end
            val = mean(draw(R,nSamples));
        end

        function limits = monteCarloPlotLimits(R,nStds,nSamples)
            % Calculate plotting limits
            if nargin < 2 || isempty(nStds)
                nStds = 2;
            end
            if nargin < 3 || isempty(nSamples)
                nSamples = []; % Let the defaults from monteCarloMean and
                               % monteCarloCovariance decide.
            end

            mu = monteCarloMean(R,nSamples);
            C = monteCarloCovariance(R,nSamples);

            minX = min(mu, [], 1)' - nStds*sqrt(diag(C));
            maxX = max(mu, [], 1)' + nStds*sqrt(diag(C));

            limits = zeros(1,2*length(minX));
            limits(1:2:end) = minX;
            limits(2:2:end) = maxX;
        end
    end
    methods (Hidden = true)
        function S = toStructure(obj)
            % toStructure(obj)
            % This default prtRv method adds all properties defined in
            % the class of obj into the structure, that are:
            %   GetAccess: public
            %   Hidden: false
            % other prtRvs (that are properties, contained in cells,
            %   or fields of structures) are also converted to structures.

            MetaInfo = meta.class.fromName(class(obj));

            propNames = {};
            for iProperty = 1:length(MetaInfo.Properties)
                if isequal(MetaInfo.Properties{iProperty}.DefiningClass,MetaInfo) && strcmpi(MetaInfo.Properties{iProperty}.GetAccess,'public') && ~MetaInfo.Properties{iProperty}.Hidden
                    propNames{end+1} = MetaInfo.Properties{iProperty}.Name; %#ok<AGROW>
                end
            end

            S.class = 'prtRv';
            S.prtRvType = class(obj);
            for iProp = 1:length(propNames)
                cProp = obj.(propNames{iProp});
                for icProp = 1:length(cProp) % Allow for arrays of objects
                    cOut = prtUtilFindPrtActionsAndConvertToStructures(cProp(icProp));
                    if icProp == 1
                        cVal = repmat(cOut,size(cProp));
                    else
                        cVal(icProp) = cOut;
                    end
                end
                S.(propNames{iProp}) = cVal;
            end
            S.userData = obj.userData;
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function K = spm_kron(A,B)
% Kronecker tensor product with sparse outputs
% FORMAT K = spm_kron(A,B)
%        K = spm_kron(A)
%
%   KRON(X,Y) is the Kronecker tensor product of X and Y.
%   The result is a large matrix formed by taking all possible
%   products between the elements of X and those of Y.   For
%   example, if X is 2 by 3, then KRON(X,Y) is
%
%      [ X(1,1)*Y  X(1,2)*Y  X(1,3)*Y
%        X(2,1)*Y  X(2,2)*Y  X(2,3)*Y ]
%
%   Class support for inputs X,Y:
%      float: double, single
%
% When called with a single cell array input, the tensor product
% is formed recursively

%   Previous versions by Paul Fackler, North Carolina State,
%   and Jordan Rosenthal, Georgia Tech.
%   Copyright 1984-2004 The MathWorks, Inc.
%   $Revision: 7760 $ $Date: 2004/06/25 18:52:18 $
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_kron.m 7760 2019-12-29 17:45:58Z karl $

% Deal with cell arrays
%--------------------------------------------------------------------------
if iscell(A)
    K = 1;
    for i = 1:numel(A)
        K = spm_kron(A{i},K);
    end
    return
end

% Kronecker tensor product
%--------------------------------------------------------------------------
[ma,na] = size(A);
[mb,nb] = size(B);

[ia,ja,sa] = find(A);
[ib,jb,sb] = find(B);
ia = ia(:); ja = ja(:); sa = sa(:);
ib = ib(:); jb = jb(:); sb = sb(:);
ka = ones(size(sa));
kb = ones(size(sb));
t  = mb*(ia-1)';
ik = t(kb,:)+ib(:,ka);
t  = nb*(ja-1)';
jk = t(kb,:)+jb(:,ka);
K  = sparse(ik,jk,sb*sa.',ma*mb,na*nb);

author: spm


text:
function [opt_radius,opt_interval,result,r,alldt,nummod,numreal,sigma] = calc_optim(a,rmin,rstep,rmax,tmin,tstep,Nmin,bootloops,maepi,choice)

    % function [opt_radius,opt_interval,result,r,alldt,nummod,numreal,sigma] = calc_optim(a,rmin,rstep,rmax,tmin,tstep,Nmin,bootloops,maepi,choice);
    % -------------------------------------------------------------------------
    %
    % optimisation of free parameters looking for seismic quiescences preceding large aftershocks
    %
    % Input parameters:
    %   a               Earthquake catalog (has to be complete in magnitude!)
    %   rmin,rstep,rmax Radius [degrees]
    %   tmin,tstep      Time [days]
    %   Nmin            minimum earthquake number in catalog
    %   bootloops       number of bootstrap loops
    %   maepi           mainshock (->ZMAP)
    %   choice          Nr. of M5+ aftershock
    %
    % Output parameters:
    %   opt_radius      radius around large aftershock
    %   opt_interval    time interval before large aftershock
    %   result          matrix containing relative seismic quiescences
    %   r, alldt        vectors with radii, forecast intervals
    %   numreal,nummod  observed/modeled nr. of aftershocks in interval
    %   sigma           uncertainty of forecast
    %
    % Samuel Neukomm
    % last update 25.02.2004

report_this_filefun(mfilename('fullpath'));
    % get mainshock / calculate delay times / cut catalogue at mainshock
    [m_main, main] = max(a.Magnitude);
    if size(a,2) == 9
        date_matlab = datenum(a.Date.Year,a.Date.Month,a.Date.Day,a.Date.Hour,a.Date.Minute,zeros(size(a,1),1));
    else
        date_matlab = datenum(a.Date.Year,a.Date.Month,a.Date.Day,a.Date.Hour,a.Date.Minute,a(:,10));
    end
    date_main = date_matlab(main);
    time_aftershock = date_matlab-date_main;
    l = time_aftershock(:) > 0;
    tas = time_aftershock(l);
    eqcatalogue = a.subset(l);

    % choose one of large aftershocks
    l = eqcatalogue(:,6) >= 5;
    largeas = eqcatalogue(l,:);
    largetime = tas(l);
    eq = largeas(choice,:);
    eqt = largetime(choice);

    % get optimum parameters
    dist_lon = eqcatalogue(:,1)-eq(1);
    dist_lat = eqcatalogue(:,2)-eq(2);
    dist_dep = km2deg(eqcatalogue(:,7)-eq(7),6371);
    result = []; i = 1;
    numreal = []; nummod = []; sigma = [];
    for r = rmin:rstep:rmax % range of collection radius
        disp(num2str(r/rmax))
        l = (dist_lon.^2+dist_lat.^2+dist_dep.^2).^0.5 <= r;
        gpi = eqcatalogue(l,:);
        time_as = tas(l);
        l = time_as < eqt;
        gpi = gpi(l,:); % sub-catalogue to be analysed
        time_as = time_as(l); % corresponding delay times

        if length(time_as) >= Nmin
            dt = tmin; % interval before large aftershock
            j = 1; alldt = [];
            while dt < eqt % Changed eqt/2 to eqt
                [rc,realnum,modnum,sig] = calc_optrc(gpi,time_as,eqt-dt,eqt,bootloops,maepi); % determine significance
                if isnan(rc)==0
                    result(i,j) = rc;
                    numreal(i,j) = realnum;
                    nummod(i,j) = modnum;
                    sigma(i,j) = sig;
                else
                    result(i,j) = 0;
                    numreal(i,j) = 0;
                    nummod(i,j) = 0;
                    sigma(i,j) = 0;
                end
                alldt = [alldt; dt];
                dt = dt+tstep;
                j = j+1;
            end
        end
        i = i+1;
    end
    r = (rmin:rstep:rmax)';

    if isempty(result) == 0  &&  max(max(result)) > 0
        % find minimum
        [dum, j] = min(min(result));
        [dum, i] = min(result(:,j));
        opt_radius = r(i);
        opt_interval = alldt(j);

        % contourplot of results
        figure_w_normalized_uicontrolunits('Numbertitle','off','Name','Optimum quiescence')
        lim = ceil(abs(dum));
        pcolor(alldt,r,result)
        set(gca,'pos',[0.15 0.15 0.7 0.7])
        shading interp
        caxis([-lim lim]);
        load fourcolors
        j = [j(256:-1:1,:) ];
        colormap(j)
        set(gcf,'renderer','zbuffer');
        set(gca,'TickDir','out')
        h5 = colorbar;
        set(h5,'Tickdir','out','pos',[0.8 0.3 0.02 0.4]);
        hold on
        plot(opt_interval,opt_radius','*','Markersize',8,'linewidth',2,'color','k');
        string =['Radius = ' num2str(opt_radius) ' deg; Interval = ' num2str(opt_interval) ' days; RC = ' num2str(dum,3) ];
        title(string)
        string = ['Time interval before t1 = ' num2str(round(100*eqt)/100) ' [days]'];
        xlabel(string)
        ylabel('Radius [deg]');
        set(gcf,'color','w');

        l = (dist_lon.^2+dist_lat.^2+dist_dep.^2).^0.5 <= opt_radius;
        optcat = eqcatalogue(l,:);
        time_as = tas(l);
        l = time_as < eqt;
        optcat = optcat(l,:);
        plot_optfit(optcat,eqt-opt_interval,opt_interval,bootloops,maepi);
    else
        opt_radius = NaN;
        opt_interval = NaN;
        alldt = NaN;
    end

author: CelsoReyes


text:
% [INPUT]
% p = A vector of floats [0,Inf) of length t representing the prices.
% bw = An integer [21,252] representing the dimension of each rolling window.
% w = An integer [500,Inf) representing the number of sweeps (optional, default=1000).
% c = A float (0,Inf) representing the starting coefficient value (optional, default=0.01).
% s2 = A float (0,Inf) representing the starting variance of innovations (optional, default=0.0004).
%
% [OUTPUT]
% ris = A column vector of floats (0,Inf) of length t representing the Roll implicit spread.

function ris = roll_implicit_spread(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('p',@(x)validateattributes(x,{'double'},{'real' 'finite' 'nonnegative' 'vector' 'nonempty'}));
        ip.addRequired('bw',@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 21 '<=' 252 'scalar'}));
        ip.addOptional('w',1000,@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 500}));
        ip.addOptional('c',0.01,@(x)validateattributes(x,{'double'},{'real' 'finite' 'positive'}));
        ip.addOptional('s2',0.0004,@(x)validateattributes(x,{'double'},{'real' 'finite' 'positive'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    p = validate_input(ipr.p);
    bw = ipr.bw;
    w = ipr.w;
    c = ipr.c;
    s2 = ipr.s2;

    nargoutchk(1,1);

    ris = roll_implicit_spread_internal(p,bw,w,c,s2);

end

function ris = roll_implicit_spread_internal(p,bw,w,c,s2)

    up = isempty(getCurrentTask());

    windows = extract_rolling_windows(log(max(1e-6,p)),bw);
    ris = zeros(numel(windows),1);

    if (up)
        parfor i = 1:numel(windows)
            ris(i) = gibbs_sampler(windows{i},w,c,s2);
        end
    else
        for i = 1:numel(windows)
            ris(i) = gibbs_sampler(windows{i},w,c,s2);
        end
    end

    alpha = 2 / (bw + 1);

    ris = [ris(1); filter(alpha,[1 (alpha - 1)],ris(2:end),(1 - alpha) * ris(1))];

end

function g = gibbs_sampler(p,w,c,s2)

    dp = diff(p);
    q = [1; sign(dp)];

    for i = 1:w
        dq = diff(q);

        d = 1 + ((1 / s2) * (dq.' * dq));
        mu = linsolve(d,(1 / s2) * (dq.' * dp));
        rho = inv(d);
        c = mvnrnd_truncated(mu,rho);

        u = dp - (c .* dq);
        alpha = 1e-12 + (numel(u) / 2);
        beta = 1e-12 + (sum(u.^2) / 2);
        s2 = 1 / ((1 / beta) * gamrnd(alpha,1));

        q = perform_draw(p,q,c,s2);
    end

    g = 2 * c;

end

function q = perform_draw(p,q,c,s2)

    t = numel(p);
    q_nnz = q ~= 0;
    s22 = s2 * 2;

    m = mod((1:t).',2);
    r = rand(t,1);

    q = [q q];
    p = [p p];

    for i = 1:2
        o = (m == (i - 1)) & q_nnz;

        if (~any(o))
            continue;
        end

        q(o,1) = 1;
        q(o,2) = -1;

        u = diff(p - (c * q));

        s = u.^2 ./ s22;
        s_sum = [s; [0 0]] + [[0 0]; s];
        s_sum = s_sum(o,:);

        odds_log = diff(s_sum,1,2);
        in_range = odds_log < 500;
        odds = exp(in_range .* odds_log);

        buy = odds ./ (1 + odds);
        buy = (in_range .* buy) + ~in_range;

        if (i == 1)
            q(o,:) = repmat(1 - (2 * (r(o) > buy)),1,2);
        else
            q(o,1) = 1 - (2 * (r(o) > buy));
        end
    end

    q = q(:,1);

end

function r = mvnrnd_truncated(mu,rho)

    f = sqrt(rho);
    low = -mu / f;

    if (low > 6)
        eta = low + (100 * eps());
    else
        plow = normcdf(low);

        p = plow + (rand() * (1 - plow));

        if (p == 1)
            eta = low + (100 * eps());
        else
            eta = norminv(p);
        end
    end

    r = mu + (f * eta);

end

function p = validate_input(p)

    p = p(:);

    if (numel(p) < 5)
        error('The value of ''p'' is invalid. Expected input to be a vector containing at least 5 elements.');
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function [fProbability, fAICc] = calc_loglikelihood_nochange2(mCat1, mCat2)
% function [fProbability, fAICc] = calc_loglikelihood_nochange2(mCat1, mCat2);
% ----------------------------------------------------------------------------------------------
% Calculate log-likelihood estimation of NO CHANGE between the two periods
%
% Incoming variable
% mCat1 : EQ catalog period 1 (Catalog to be modified)
% mCat2 : EQ catalog period 2 (Observed catalog)
%
% Outgoing variable
% fProbability : log-likelihood probabilty
% fAICc         : Corrected Akaike Information Criterion
%
% Author: J. Woessner, woessner@seismo.ifg,.ethz.ch
% last update: 26.11.02

fBinning = 0.1;

% Determine exact time period
fPeriod1 = max(mCat1(:,3)) - min(mCat1(:,3));
fPeriod2 = max(mCat2(:,3)) - min(mCat2(:,3));

% Initialize values
fMinMag = min([min(mCat1(:,6)) min(mCat2(:,6))]);
fMaxMag = max([max(mCat1(:,6)) max(mCat2(:,6))]);

try
    %% Calculate model for best fiting Mc
    [mResult, fMls, fMc, fMu, fSigma, mDatPred, vPredBest, fBvalue] = calc_McCdfnormal(mCat1, fBinning);
    vPredFMD = mDatPred(:,1)'./fPeriod1;
    vMags = mDatPred(:,2)';
    % FMD to be modeled
    [vObsFMD,vBin2] = hist(mCat2(:,6),min(vMags):0.1:max(vMags));
    vObsFMD = ceil(vObsFMD./fPeriod2);
    % Calculate the likelihoods for both models
    vProb_ = calc_log10poisspdf2(vObsFMD', vPredFMD');
    % Sum the probabilities
    fProbability = (-1) * sum(vProb_);

    nDegFree = 0; % degree of freedom
    n_samples = length(mCat2(:,6));
    %% Corrected Akaike Information Criterion (AICc)
    fAICc = -2*(-fProbability)+2*nDegFree+2*nDegFree*(nDegFree+1)/(n_samples-nDegFree-1);
catch
    fAICc = nan;
    fProbability = nan;
end


text:
function [t,x_new,f_new,g_new,funEvals,H] = ArmijoBacktrack(...
    x,t,d,f,fr,g,gtd,c1,LS,tolX,debug,doPlot,saveHessianComp,funObj,varargin)
%
% Backtracking linesearch to satisfy Armijo condition
%
% Inputs:
%   x: starting location
%   t: initial step size
%   d: descent direction
%   f: function value at starting location
%   fr: reference function value (usually funObj(x))
%   gtd: directional derivative at starting location
%   c1: sufficient decrease parameter
%   debug: display debugging information
%   LS: type of interpolation
%   tolX: minimum allowable step length
%   doPlot: do a graphical display of interpolation
%   funObj: objective function
%   varargin: parameters of objective function
%
% Outputs:
%   t: step length
%   f_new: function value at x+t*d
%   g_new: gradient value at x+t*d
%   funEvals: number function evaluations performed by line search
%   H: Hessian at initial guess (only computed if requested

% Evaluate the Objective and Gradient at the Initial Step
if nargout == 6
    [f_new,g_new,H] = feval(funObj, x + t*d, varargin{:});
else
    [f_new,g_new] = feval(funObj, x + t*d, varargin{:});
end
funEvals = 1;

while f_new > fr + c1*t*gtd || ~isLegal(f_new)

    temp = t;
    if LS == 0 || ~isLegal(f_new)
        % Backtrack w/ fixed backtracking rate
        if debug
            fprintf('Fixed BT\n');
        end
        t = 0.5*t;
    elseif LS == 2 && isLegal(g_new)
        % Backtracking w/ cubic interpolation w/ derivative
        if debug
            fprintf('Grad-Cubic BT\n');
        end
        t = polyinterp([0 f gtd; t f_new g_new'*d],doPlot);
    elseif funEvals < 2 || ~isLegal(f_prev)
        % Backtracking w/ quadratic interpolation (no derivative at new point)
        if debug
            fprintf('Quad BT\n');
        end
        t = polyinterp([0 f gtd; t f_new sqrt(-1)],doPlot);
    else%if LS == 1
        % Backtracking w/ cubic interpolation (no derivatives at new points)
        if debug
            fprintf('Cubic BT\n');
        end
        t = polyinterp([0 f gtd; t f_new sqrt(-1); t_prev f_prev sqrt(-1)],doPlot);
    end

    % Adjust if change in t is too small/large

    if t < temp*1e-3
        if debug
            fprintf('Interpolated Value Too Small, Adjusting\n');
        end
        t = temp*1e-3;
    elseif t > temp*0.6
        if debug
            fprintf('Interpolated Value Too Large, Adjusting\n');
        end
        t = temp*0.6;
    end

    f_prev = f_new;
    t_prev = temp;
    if ~saveHessianComp && nargout == 6
        [f_new,g_new,H] = feval(funObj, x + t*d, varargin{:});
    else
        [f_new,g_new] = feval(funObj, x + t*d, varargin{:});
    end
    funEvals = funEvals+1;

    % Check whether step size has become too small
    if sum(abs(t*d)) <= tolX
        if debug
            fprintf('Backtracking Line Search Failed\n');
        end
        t = 0;
        f_new = f;
        g_new = g;
        break;
    end
end

% Evaluate Hessian at new point
if nargout == 6 && funEvals > 1 && saveHessianComp
    [f_new,g_new,H] = feval(funObj, x + t*d, varargin{:});
    funEvals = funEvals+1;
end

x_new = x + t*d;

end

author: singaxiong

repo: SignalGraph


text:
% IMSHAPES Generation of images representing geometric shapes
% Version 1.1 01-Jul-2011 .
%
%   This module contains several function for creating synthetic images of
%   common shapes, in 2D (ellipses, rectangles...), or in 3D (balls and
%   ellipsoids, cuboids, torus...).
%   Such discrete representation can be used as phantom to test and / or
%   validate algorithms of binary images measurements.
%
%   Requires the 'geom2d', the 'geom3d', and the 'imFilters' toolboxes.
%
% Shapes 2D:
%   discretePoints         - Discretize a set of points
%   discreteDisc           - Discretize a disc
%   discreteEllipse        - Discretize a planar ellipse
%   discreteSquare         - Discretize a planar square
%   discreteRectangle      - Discretize a planar rectangle
%   discreteCapsule        - Discretize a planar capsule
%   discretePolygon        - Discretize a planar polygon
%   discretePolyline       - Discretize a planar polyline
%   discreteHalfPlane      - Discretize a half plane
%   discreteParabola       - Discretize a planar parabola
%   discreteEgg            - Discretize a planar egg
%   discreteStarfish       - Discretize a starfish curve
%   discreteTrefoil        - Discretize a trefoil curve
%   discreteCurve          - Discretize a planar curve
%
% Shapes 3D:
%   discreteBall           - Discretize a 3D Ball
%   discreteHalfBall       - discretize a 3D half-ball
%   discreteEllipsoid      - discretize a 3D ellipsoid
%   discreteCuboid         - discretize a 3D cuboid
%   discreteCube           - discretize a 3D cube
%   discreteTorus          - Discretize a 3D Torus
%   discreteCylinder       - Discretize a 3D cylinder
%   discreteCapsule3d      - Create binary image of a 3D capsule
%   discreteReuleauxRevol  - Discretize the revolution of a Reuleaux triangle
%   discreteSphereEighth   - Discretize a 3D sphere eighth
%
% Tessellations:
%   imPointsInfluenceZones - Maps influence zones of a set of 2D/3D points
%   imvoronoi2d            - Generate a 2D voronoi image from a set of points
%   imvoronoi3d            - generate a 3D voronoi image from a set of points
%   dilatedVoronoi         - Simulate a 'thick' voronoi tesselation
%   imAWVoronoi            - generate Additively Weighted Voronoi Diagram image
%   imPowerDiagram         - Power diagramm of a set of points
%
%
% Author: David Legland
% e-mail: david.legland@inra.fr
% http://github.com/mattools/matImage
% Copyright 2009 INRA - Cepia Software Platform.

help('imShapes');

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% problem 11 - graph of r[n] and of r[n+1]

% r[n]
 n=-3:6;
 u=(n>=0);
 r=n.*u;

 stem(n,r)
 title('Unit ramp sequence r[n]')
 ylim([-.1 6.1])

% r[n+1]
 figure
 n=-3:6;
 u1=(n>=-1);
 r=(n+1).*u1;

 stem(n,r)
 title('Unit ramp sequence r[n+1]')
 ylim([-.1 7.1])


text:
function [imVol, mmPerVox, gradDirs, bVals, xformToCannonical] = dtiAverageRawData(firstDataFile, gradsFile)
% Intelligent averaging of diffusion-weighted images.
%
% [imVol, mmPerVox, gradDirs, bVals, xformToCannonical] = dtiAverageRawData(firstDataFile, gradsFile);
%
% Loads a whole directory of raw image files and combines them according to
% the gradient direction information in gradsFile. All repeated
% measurements are averaged. If coreg is true (default), it will also try
% to align each b>0 image to the b=0 image using a 12-parameter
% mutual-information algorithm (from SPM2). This should remove much of the
% eddy-current distortion.
%
% Also returns the xform that will rotate the images into our cannonical
% orientation (as close to axial as possible- see
% computeCannonicalXformFromIfile).
%
% TODO:
% * coreg should also do motion correction. To implement this, we need to
% do the coregistration step before we combine the repeated measurements.
% Of course, we should also correct the gradDirs for non-trivial motion.
%
%
% HISTORY:
% 2005.05.25: RFD (bob@sirl.stanford.edu) wrote it.
%

if(~exist('gradsFile','var')|isempty(gradsFile))
    % We could also try to infer the file based on info in the dicom
    % header.
    if(isunix)
      default = '/usr/local/dti/diffusion_grads/dwepi.13.grads';
    else
      default = pwd;
    end
    [f,p] = uigetfile({'*.grads';'*.*'},'Select dwepi grads file...',default);
    if(isequal(f,0)||isequal(p,0))
        error('User cancelled.');
    end
    gradsFile = fullfile(p,f);
end
if(~exist('firstDataFile','var')|isempty(firstDataFile))
    [f,p] = uigetfile({'*.dcm';'*.*'},'Select the first data file...');
    if(isequal(f,0)||isequal(p,0))
        error('User cancelled.');
    end
    firstDataFile = fullfile(p,f);
end

% Flag to turn on MI-based coregistration (eddy-current correction)
if(~exist('coreg','var')|isempty(coreg))
    coreg = false;
end

gradDirs = dlmread(gradsFile);
scannerToPhysXform = computeXformFromIfile(firstDataFile);
[xformToCannonical, baseName, mmPerVox, imDim, notes, sliceSkip] = computeCannonicalXformFromIfile(firstDataFile);
% *** TODO: may want to include im_hdr.scanspacing!
%mmPerVox(3) = mmPerVox(3)+sliceSkip;

% *** TODO: We assume that all non-zero b-vals are the same
hdr = dicominfo(firstDataFile);
bVals = hdr.Private_0019_10b0;

% HACK! TO deal with our interleaved data:
imDim(3) = imDim(3)*2;

[dataDir,f,ext] = fileparts(firstDataFile);
allFiles = dir(fullfile(dataDir, ['*' ext]));
nDirs = size(gradDirs,1);
nSlices = imDim(3);
filesPerRep = nDirs*nSlices;
nReps = length(allFiles)/filesPerRep;
if(nReps~=round(nReps))
    error('Non-integer number of repeats. Something is funky- refusing to go on...');
end

% TO DO:
% * motion correction/outlier rejection
%
disp('Loading all the data files...');
imVol = zeros([imDim nDirs]);
for(d=1:nDirs)
  for(n=1:nReps)
    fn = [[1:nSlices]+(d-1)*nSlices+((n-1)*filesPerRep)];
    for s = 1:nSlices;
      curName = sprintf('I%04d.dcm', fn(s));
      tmpImg = readRawImage(fullfile(dataDir,curName), 0, nPix(1:2), 'b');
      imVol(:,:,s,d) = imVol(:,:,s,d)+double(tmpImg);
    end
  end
  %tmpIm = makeCubeIfiles(dataDir, nPix(1:2), fileNums);
end
imVol = imVol./nReps;

disp('Merging duplicate directions...');
% Merge duplicates
[gDirs,inds] = unique(gradDirs,'rows');
% undo unique's sorting
gDirs = gradDirs(sort(inds),:);
nDirs = size(gDirs,1);
junk= [];
for(ii=1:nDirs)
  dups = gDirs(ii,1)==gradDirs(:,1)&gDirs(ii,2)==gradDirs(:,2)& ...
         gDirs(ii,3)==gradDirs(:,3);
  dups = find(dups);
  if(length(dups)>1)
    imVol(:,:,:,dups(1)) = mean(imVol(:,:,:,dups),4);
    junk = [junk dups(2:end)];
  end
end
imVol(:,:,:,junk) = [];

if(coreg)
  disp('Coregistering...');
  % NOTE: The following MI-based image registration tries to
  % coregister each direction map to the b0. It doen't seem to do
  % anything bad, but isn't particularly effective at removing the
  % eddy-current distortion. I guess it's not a 3d affine xform?
  % But, if we ran it on the individual images before averaging, it
  % would probably be a good way to remove motion.
  b0 = imVol(:,:,:,1);
  b0 = uint8(mrAnatHistogramClip(b0, 0.4, 0.99)*255+0.5);
  h = mrvWaitbar(0,'Performing eddy current correction...');
  for(d=2:nDirs)
    xf{d} = mrAnatRegister(imVol(:,:,:,d), b0);
    tmpIm = mrAnatResliceSpm(imVol(:,:,:,d),xf{d},[],[1 1 1],[7 7 7 0 0 0],0);
    tmpIm(isnan(tmpIm)) = 0;
    tmpIm(tmpIm<0) = 0;
    imVol(:,:,:,d) = tmpIm;
    mrvWaitbar(d/nDirs,h);
  end
  close(h);
end
imVol = int16(round(imVol));
return;

% Showing a movie (to see the image mis-alignment)
for(d=1:nDirs)
  f = squeeze(imVol(:,:,20,d));
  f = uint8(f./max(f(:)).*255+0.5);
  m(d) = im2frame(f,gray(256));
end
figure(99);
image(f); colormap(gray(256));
truesize; axis off
movie(m,1,5);


text:
function files=cosmo_dir(varargin)
% list files recursively in a directory
%
% files=cosmo_dir([directory][,file_pattern])
%
% Inputs:
%   directory      optional parent directory to look for files. If omitted
%                  the current directory is used.
%   file_pattern   optional pattern to match file names. The wildcards '*'
%                  (zero or more of any character) and '?' (any single
%                  character) can be used. If omitted, '*' is used, meaning
%                  that all files are listed. A file pattern may be
%                  prefixed by a directory.
% Output:
%   files          Px1 struct (if P files were found in the directory or
%                  recursively in any (sub^X)-directorie) with fields
%                  .name, .date, .bytes, .isdir, .datenum. Unlike the
%                  built-in 'dir' function, .name contains the path and
%                  filename.
%
% Notes:
%  - the output is similar to built-in 'dir', except that .name contains
%    the path of the files as well.
%  - to list files in a directory but not those in its subdirectories, use
%    the built-in 'dir' function.
%  - this function does not return '.' (current directory) or '..' (parent
%    directory.
%
% Examples:
%  % list recursively all files in the current directory
%  d=cosmo_dir();
%
%  % Assuming that the directory 'my_dir' exists, list recursively all
%  % files the directory contains
%  d=cosmo_dir('my_dir');
%
%  % Assuming that the directory 'my_file' does not exist, return
%  % either a struct with a single entry (if a file named 'my_file'
%  % exists), or an emtpy struct (if no such file exists)
%  d=cosmo_dir('my_file');%
%
%  % list recursively all files in the current directory for which the name
%  % ends with '.jpg'
%  d=cosmo_dir('*.jpg')
%
%  % list recursively all files in the directory 'my_dir' with extension
%  % '.jpg'
%  d=cosmo_dir('my_dir', '*.jpg')
%
%  % as above but with a single argument; list recursively all files in
%  % the directory 'my_dir' with extension '.jpg'
%  d=cosmo_dir('my_dir/*.jpg')
%
%  % list recursively all files in my_dir that are two characters long,
%  % followed by the extension '.jpg'
%  d=cosmo_dir('my_dir', '??.jpg')
%
%  % list recursively all files in my_dir that start with 'a', followed by
%  % any character, followed by 'b', followed by any number of characters,
%  % followed by '.jpg'
%  d=cosmo_dir('my_dir', 'a?b*.jpg')
%
% See also: dir
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    [root_dir,file_pat]=process_arguments(varargin{:});

    file_re=['^' ... % start of the string
             regexptranslate('wildcard',file_pat) ...
             '$'];   % end of the string

    files=find_files_recursively(root_dir,file_re);

function [root_dir,file_pat]=process_arguments(varargin)
    narg=numel(varargin);
    if narg>2
        error('need at most two arguments');
    end

    if ~iscellstr(varargin)
        error('arguments must be strings');
    end

    root_dir='';
    file_pat='*';

    if narg>=1
        arg1=varargin{1};
        if isdir(arg1)
            root_dir=arg1;
            if narg>=2
                file_pat=varargin{2};
            end
        else
            if narg>=2
                error('Directory not found: ''%s''', arg1);
            end
            [root_dir,nm,ext]=fileparts(arg1);
            file_pat=[nm ext];
        end
    end

function res=find_files_recursively(root_dir,file_re)
    if isempty(root_dir)
        dir_arg={};
    else
        dir_arg={root_dir};
    end

    d=dir(dir_arg{:});
    n=numel(d);

    res_cell=cell(n,1);
    for k=1:n
        d_k=d(k);
        fn=d_k.name;

        % do not return current and parent directorys
        ignore_fn=strcmp(fn,'.') || strcmp(fn,'..');

        if ~ignore_fn
            path_fn=fullfile(root_dir, fn);

            if isdir(path_fn)
                % recursive call
                res=find_files_recursively(path_fn,file_re);

            elseif ~isempty(regexp(fn,file_re,'once'))
                % file matches the pattern
                d_k.name=path_fn;
                res=d_k;

            else
                % ignore file
                continue;
            end

            if ~isempty(res)
                res_cell{k}=res;
            end
        end
    end

    keep_msk=~cellfun(@isempty,res_cell);

    if any(keep_msk)
        % at least one file found
        res_keep=res_cell(keep_msk);
        res=cat(1,res_keep{:});
    else
        % no files found, return empty struct
        if cosmo_wtf('is_octave')
            labels={'name','date','bytes','isdir',...
                        'datenum','statinfo'};
        else
            labels={'name','date','bytes','isdir',...
                        'datenum'};
        end
        n_labels=numel(labels);
        res=cell2struct(cell(n_labels,0),labels);
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
clc;
clear;
%全文原理介绍见：https://zhuanlan.zhihu.com/p/57967971
%————————————————————————————————————————————————————————%
%q1:ifft点数难道不是应该等于子载波数吗？子载波数与ifft点数的关系？
%a:ifft点数等于子载波数
%q2：对矩阵进行fft？
%a:y可以是一向量或矩阵，若y为向量，则Y是y的FFT，并且与y具有相同的长度。若y为一矩阵，则Y是对矩阵的每一列向量进行FFT。
%q3：怎么对ofdm信号上变频
%————————————————————————————————————————————————————————%

%% 参数设置

N_sc=52;      %系统子载波数（不包括直流载波）、number of subcarrierA
N_fft=64;            % FFT 长度
N_cp=16;             % 循环前缀长度、Cyclic prefix
N_symbo=N_fft+N_cp;        % 1个完整OFDM符号长度
N_c=53;             % 包含直流载波的总的子载波数、number of carriers
M=4;               %4PSK调制
SNR=0:1:25;         %仿真信噪比
N_frm=10;            % 每种信噪比下的仿真帧数、frame
Nd=6;               % 每帧包含的OFDM符号数
P_f_inter=6;      %导频间隔
data_station=[];    %导频位置
L=7;                %卷积码约束长度
tblen=6*L;          %Viterbi译码器回溯深度
stage = 3;          % m序列的阶数
ptap1 = [1 3];      % m序列的寄存器连接方式
regi1 = [1 1 1];    % m序列的寄存器初始值

%% 基带数据数据产生
P_data=randi([0 1],1,N_sc*Nd*N_frm);

%% 信道编码（卷积码、或交织器）
%卷积码：前向纠错非线性码
%交织：使突发错误最大限度的分散化
trellis = poly2trellis(7,[133 171]);       %(2,1,7)卷积编码
code_data=convenc(P_data,trellis);

%% qpsk调制
data_temp1= reshape(code_data,log2(M),[])';             %以每组2比特进行分组，M=4
data_temp2= bi2de(data_temp1);                             %二进制转化为十进制
modu_data=pskmod(data_temp2,M,pi/M);              % 4PSK调制
% figure(1);
scatterplot(modu_data),grid;                  %星座图(也可以取实部用plot函数)

%% 扩频
%————————————————————————————————————————————————————————%
%扩频通信信号所占有的频带宽度远大于所传信息必需的最小带宽
%根据香农定理，扩频通信就是用宽带传输技术来换取信噪比上的好处，这就是扩频通信的基本思想和理论依据。
%扩频就是将一系列正交的码字与基带调制信号内积
%扩频后数字频率变成了原来的m倍。码片数量 = 2（符号数）* m（扩频系数）
%————————————————————————————————————————————————————————%

code = mseq(stage,ptap1,regi1,N_sc);     % 扩频码的生成
code = code * 2 - 1;         %将1、0变换为1、-1
modu_data=reshape(modu_data,N_sc,length(modu_data)/N_sc);
spread_data = spread(modu_data,code);        % 扩频
spread_data=reshape(spread_data,[],1);

%% 插入导频
P_f=3+3*1i;                       %Pilot frequency
P_f_station=[1:P_f_inter:N_fft];%导频位置（导频位置很重要，why?）
pilot_num=length(P_f_station);%导频数量

for img=1:N_fft                        %数据位置
    if mod(img,P_f_inter)~=1          %mod(a,b)就是求的是a除以b的余数
        data_station=[data_station,img];
    end
end
data_row=length(data_station);
data_col=ceil(length(spread_data)/data_row);

pilot_seq=ones(pilot_num,data_col)*P_f;%将导频放入矩阵
data=zeros(N_fft,data_col);%预设整个矩阵
data(P_f_station(1:end),:)=pilot_seq;%对pilot_seq按行取

if data_row*data_col>length(spread_data)
    data2=[spread_data;zeros(data_row*data_col-length(spread_data),1)];%将数据矩阵补齐，补0是虚载频~
end;

%% 串并转换
data_seq=reshape(data2,data_row,data_col);
data(data_station(1:end),:)=data_seq;%将导频与数据合并

%% IFFT
ifft_data=ifft(data);

%% 插入保护间隔、循环前缀
Tx_cd=[ifft_data(N_fft-N_cp+1:end,:);ifft_data];%把ifft的末尾N_cp个数补充到最前面

%% 并串转换
Tx_data=reshape(Tx_cd,[],1);%由于传输需要

%% 信道（通过多经瑞利信道、或信号经过AWGN信道）
 Ber=zeros(1,length(SNR));
 Ber2=zeros(1,length(SNR));
for jj=1:length(SNR)
    rx_channel=awgn(Tx_data,SNR(jj),'measured');%添加高斯白噪声

%% 串并转换
    Rx_data1=reshape(rx_channel,N_fft+N_cp,[]);

%% 去掉保护间隔、循环前缀
    Rx_data2=Rx_data1(N_cp+1:end,:);

%% FFT
    fft_data=fft(Rx_data2);

%% 信道估计与插值（均衡）
    data3=fft_data(1:N_fft,:);
    Rx_pilot=data3(P_f_station(1:end),:); %接收到的导频
    h=Rx_pilot./pilot_seq;
    H=interp1( P_f_station(1:end)',h,data_station(1:end)','linear','extrap');%分段线性插值：插值点处函数值由连接其最邻近的两侧点的线性函数预测。对超出已知点集的插值点用指定插值方法计算函数值

%% 信道校正
    data_aftereq=data3(data_station(1:end),:)./H;
%% 并串转换
    data_aftereq=reshape(data_aftereq,[],1);
    data_aftereq=data_aftereq(1:length(spread_data));
    data_aftereq=reshape(data_aftereq,N_sc,length(data_aftereq)/N_sc);

%% 解扩
    demspread_data = despread(data_aftereq,code);       % 数据解扩

%% QPSK解调
    demodulation_data=pskdemod(demspread_data,M,pi/M);
    De_data1 = reshape(demodulation_data,[],1);
    De_data2 = de2bi(De_data1);
    De_Bit = reshape(De_data2',1,[]);

%% （解交织）
%% 信道译码（维特比译码）
    trellis = poly2trellis(7,[133 171]);
    rx_c_de = vitdec(De_Bit,trellis,tblen,'trunc','hard');   %硬判决

%% 计算误码率
    [err,Ber2(jj)] = biterr(De_Bit(1:length(code_data)),code_data);%译码前的误码率
    [err, Ber(jj)] = biterr(rx_c_de(1:length(P_data)),P_data);%译码后的误码率

end
 figure(2);
 semilogy(SNR,Ber2,'b-s');
 hold on;
 semilogy(SNR,Ber,'r-o');
 hold on;
 legend('4PSK调制、卷积码译码前（有扩频）','4PSK调制、卷积码译码后（有扩频）');
 hold on;
 xlabel('SNR');
 ylabel('BER');
 title('AWGN信道下误比特率曲线');

 figure(3)
 subplot(2,1,1);
 x=0:1:30;
 stem(x,P_data(1:31));
 ylabel('amplitude');
 title('发送数据（以前30个数据为例)');
 legend('4PSK调制、卷积译码、有扩频');

 subplot(2,1,2);
 x=0:1:30;
 stem(x,rx_c_de(1:31));
 ylabel('amplitude');
 title('接收数据(以前30个数据为例)');
 legend('4PSK调制、卷积译码、有扩频');

repo: OFDM


text:
function [J, grad] = linearRegCostFunction(X, y, theta, lambda)
%LINEARREGCOSTFUNCTION Compute cost and gradient for regularized linear
%regression with multiple variables
%   [J, grad] = LINEARREGCOSTFUNCTION(X, y, theta, lambda) computes the
%   cost of using theta as the parameter for linear regression to fit the
%   data points in X and y. Returns the cost in J and the gradient in grad

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

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

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost and gradient of regularized linear
%               regression for a particular choice of theta.
%
%               You should set J to the cost and grad to the gradient.
%

H = X*theta;
J = sum((H - y).^2) / (2 * m) + lambda*sum(theta(2:end).^2) / (2 * m);

grad = X'*(H - y) / m + lambda*[0;theta(2:end)] / m;

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

grad = grad(:);

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function D = eig3(A)
% function D = eig3(A)
%
% Compute in one shot the eigen-values of multiples (3 x 3) matrices
%
% INPUT:
%   A: (3 x 3 x n) array
% OUTPUT:
%   D: (3 x n). EIG3 returns in D(:,k) three eigen-values of A(:,:,k)
%
% See also: CardanRoots, eig2, eig
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% History:
%     Original 20-May-2010

if size(A,1) ~= 3 || size(A,2) ~= 3
    error('A must be [3x3xn] array');
end

A = reshape(A, 9, []).';

P3 = 1;
% Trace
P2 = -(A(:,1)+A(:,5)+A(:,9));
% Principal minors
M11 = A(:,5).*A(:,9) - A(:,8).*A(:,6);
M22 = A(:,9).*A(:,1) - A(:,3).*A(:,7);
M33 = A(:,1).*A(:,5) - A(:,4).*A(:,2);
P1 = (M11 + M22 + M33);
% Determinant
P0 = - A(:,1).*M11 ...
     + A(:,4).*(A(:,2).*A(:,9)-A(:,8).*A(:,3)) ...
     - A(:,7).*(A(:,2).*A(:,6)-A(:,5).*A(:,3));

D = CardanRoots(P3, P2, P1, P0).';

end


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

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

% Last Modified by GUIDE v2.5 16-Jan-2009 09:06:46

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

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

% --- Executes just before VAudio is made visible.
function VAudio_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to VAudio (see VARARGIN)

% Choose default command line output for VAudio
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes VAudio wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = VAudio_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on selection change in popupmenu1.
function varargout = popupmenu1_Callback(hObject, eventdata, handles, varargin)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns popupmenu1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from popupmenu1

gdt = get(handles.popupmenu1,'value')

switch gdt
    case 1
         y=audiorecorder(8000,8,2);
         set(handles.text3,'string','SampleRate=8000, SampleSize=8, Stereo')
         handles.data=y;
    guidata(handles.popupmenu1,handles)
    case 2
        y=audiorecorder(8000,16,2);
        set(handles.text3,'string','SampleRate=8000, SampleSize=16, Stereo')
        handles.data=y;
    guidata(handles.popupmenu1,handles)
    case 3
        y=audiorecorder(11025,8,2);
        set(handles.text3,'string','SampleRate=11025, SampleSize=8, Stereo')
        handles.data=y;
    guidata(handles.popupmenu1,handles)
    case 4
        y=audiorecorder(11025,16,2);
        set(handles.text3,'string','SampleRate=11025, SampleSize=16, Stereo')
        handles.data=y;
    guidata(handles.popupmenu1,handles)
    case 5
        y=audiorecorder(22050,8,2);
        set(handles.text3,'string','SampleRate=22050, SampleSize=8, Stereo')
        handles.data=y;
    guidata(handles.popupmenu1,handles)
    case 6
        y=audiorecorder(22050,16,2);
        set(handles.text3,'string','SampleRate=22050, SampleSize=16, Stereo')
        handles.data=y;
    guidata(handles.popupmenu1,handles)
    case 7
        y=audiorecorder(44100,8,2);
        set(handles.text3,'string','SampleRate=44100, SampleSize=8, Stereo')
        handles.data=y;
    guidata(handles.popupmenu1,handles)
    case 8
        y=audiorecorder(44100,16,2);
        set(handles.text3,'string','SampleRate=44100, SampleSize=16, Stereo')
        handles.data=y;
    guidata(handles.popupmenu1,handles)
end

% --- Executes during object creation, after setting all properties.
function popupmenu1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to popupmenu1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in togglebutton1.
function varargout = togglebutton1_Callback(h, eventdata, handles, varargin)
% hObject    handle to togglebutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of togglebutton1

if (get(handles.togglebutton1,'value')) == get(handles.togglebutton1,'Max')

     y=audiorecorder(8000,16,2);
% y=handles.data;
% y.samplerate=za
% set(handles.text3,'string',za)
    record(y);
    set(handles.text3,'string','Now Recording...PLease speak into mic')
     handles.data=y;
     guidata(handles.togglebutton1,handles)
elseif (get(handles.togglebutton1,'value')) == get(handles.togglebutton1,'Min')
    y=handles.data;
    stop(y);
    set(handles.text3,'string','Recording Stopped')
    pd=getaudiodata(y,'int16');
%      t=0:0.001;100;
    plot(pd,'k')
    grid on;
%     axis([0 10])
    xlabel('Time (Second)');
    ylabel('Gain (Decibel)');
end

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% % % y=handles.data;
% % % if (get(handles.togglebutton1,'value')) == get(handles.togglebutton1,'Max') && (get(handles.togglebutton2,'value')) == get(handles.togglebutton2,'Min')
% % %     set(handles.text3,'string','Cannot Play, Recording is in progress...')
% % % elseif (get(handles.togglebutton1,'value')) == get(handles.togglebutton1,'Max') && (get(handles.togglebutton2,'value')) == get(handles.togglebutton2,'Max')
% % %     play(y);
% % %     set(handles.text3,'string','Now Playing...')
% % % elseif (get(handles.togglebutton1,'value')) == get(handles.togglebutton1,'Min')
% % %     play(y);
% % %     set(handles.text3,'string','Now Playing...')
% % % end

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

sure

% --- Executes on button press in togglebutton2.
function varargout = togglebutton2_Callback(h, eventdata, handles, varargin)
% hObject    handle to togglebutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of togglebutton2

if (get(handles.togglebutton2,'value')) == get(handles.togglebutton2,'Max')
    y=handles.data;
    pause(y);
    set(handles.text3,'string','Recording Paused')
elseif (get(handles.togglebutton2,'value')) == get(handles.togglebutton2,'Min')
    y=handles.data;
    resume(y);
    set(handles.text3,'string','Recording Resumed')
end

% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% % % sure3

% --------------------------------------------------------------------
function Untitled_5_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

[F] = uigetfile('*.fig','Select any .fig file');
open(F);

% --------------------------------------------------------------------
function Untitled_4_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

close('VAudio');
close('sure');

% --------------------------------------------------------------------
function Untitled_6_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Untitled_1_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Untitled_2_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Untitled_3_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Untitled_7_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

y=audiorecorder(8000,16,2);
% y=handles.data;
% y.samplerate=za
% set(handles.text3,'string',za)
record(y);
set(handles.text3,'string','Now Recording...PLease speak into mic')
handles.data=y;
guidata(handles.togglebutton1,handles)

% --------------------------------------------------------------------
function Untitled_8_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_8 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

y=handles.data;
stop(y);
set(handles.text3,'string','Recording Stopped')
pd=getaudiodata(y,'int16');
%      t=0:0.001;100;
plot(pd,'k')
grid on;
%     axis([0 10])
xlabel('Time (Second)');
ylabel('Gain (Decibel)');

% --------------------------------------------------------------------
function Untitled_9_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_9 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

y=handles.data;
pause(y);
set(handles.text3,'string','Recording Paused')

% --------------------------------------------------------------------
function Untitled_10_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_10 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

y=handles.data;
resume(y);
set(handles.text3,'string','Recording Resumed')

% --------------------------------------------------------------------
function Untitled_11_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_11 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Untitled_12_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_12 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

open heelp.txt

% --------------------------------------------------------------------
function Untitled_13_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_13 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Untitled_16_Callback(hObject, eventdata, handles)
% hObject    handle to Untitled_16 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

sure


text:
function [ kn, fn, wn ] = r4_nor_setup ( )

%*****************************************************************************80
%
%% R4_NOR_SETUP sets data needed by R4_NOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2013
%
%  Author:
%
%    Original C version by George Marsaglia, Wai Wan Tsang
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    George Marsaglia, Wai Wan Tsang,
%    The Ziggurat Method for Generating Random Variables,
%    Journal of Statistical Software,
%    Volume 5, Number 8, October 2000, seven pages.
%
%  Parameters:
%
%    Output, uint32 KN(128), data needed by R4_NOR.
%
%    Output, real FN(128), WN(128), data needed by R4_NOR.
%
  kn = zeros ( 128, 'uint32' );
  fn = zeros ( 128 );
  wn = zeros ( 128 );

  m1 = 2147483648.0;
  vn = 9.91256303526217E-03;

  dn = 3.442619855899;
  tn = 3.442619855899;

  q = vn / exp ( - 0.5 * dn * dn );
  kn(1) = uint32 ( ( dn / q ) * m1 );
  kn(2) = 0;

  wn(1) = q / m1;
  wn(128) = dn / m1;

  fn(1) = 1.0;
  fn(128) = exp ( - 0.5 * dn * dn );

  for i = 127 : -1 : 2
    dn = sqrt ( - 2.0 * log ( vn / dn + exp ( - 0.5 * dn * dn ) ) );
    kn(i+1) = uint32 ( ( dn / tn ) * m1 );
    tn = dn;
    fn(i) = exp ( - 0.5 * dn * dn );
    wn(i) = dn / m1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function D = driving_function_mono_localwfs_sbl(x0,xs,src,f,conf)
%DRIVING_FUNCTION_MONO_LOCALWFS_SBL driving signal for local WFS using spatial
%bandwidth limitation
%
%   Usage: D = driving_function_mono_localwfs_sbl(x0,xs,src,f,conf)
%
%   Input parameters:
%       x0          - position and direction of the secondary source / m [nx7]
%       xs          - position of point source or direction of plane
%                     wave / m [1x3]
%       src         - source type of the virtual source
%                         'pw' - plane wave
%                         'ps' - point source
%       f           - frequency of the monochromatic source / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function [nx1]
%
%   See also: sound_field_mono, sound_field_mono_localwfs_sbl

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 5;
nargmax = 5;
narginchk(nargmin,nargmax);
isargsecondarysource(x0);
isargxs(xs);
isargchar(src);
isargpositivescalar(f);
isargstruct(conf);

%% ===== Computation ==========================================================
switch src
case 'ps'
    D = driving_function_mono_localwfs_sbl_ps(x0,xs,f,conf);
case 'pw'
    xs = xs./norm(xs);
    D = driving_function_mono_localwfs_sbl_pw(x0,xs,f,conf);
otherwise
    error('%s: %s is not a known source type.',upper(mfilename),src);
end

author: sfstoolbox

repo: sfs-matlab


text:
% run this script AFTER using talairach space utility (tsu.m) to get clusters
% this copies the orignial tsu stuff in ihb_TalSpaceView3D.m, but in the workspace

%------------------------------------------------------------------------------
% Get cluster information
%------------------------------------------------------------------------------
clusters = getappdata(hFigMain, 'clusters');
indClusterToView = getappdata(hFigMain, 'indClusterToView');
cl = clusters(indClusterToView);
iv = cl.hThreshold;

%------------------------------------------------------------------------------
% Figure does not exist - create new one
%------------------------------------------------------------------------------
	hFig3D = figure(...
		'Visible', 'off',...
		'NumberTitle', 'off',...
		'MenuBar', 'none',...
		'Name', '3D view of cluster',...
        'Renderer', ihbdfl_renderer,...
		'Tag', 'ihb_TalSpace3D_fig'...
		);

	%------------------------------------------------------------------------------
    % Create and/or adjust light position
	%------------------------------------------------------------------------------
    hLight = findobj('Type', 'light', 'Tag', 'ihb_lightTalSpace3D');
    cameraPos = get(hAxis3D, 'CameraPosition');
    if isempty(hLight)
	    hLight = light('Tag', 'ihb_lightTalSpace3D', 'Position', cameraPos);
    else
        set(hLight, 'Position', cameraPos);
    end
	%------------------------------------------------------------------------------
    % set callback to light follow the camera
	%------------------------------------------------------------------------------
	set(hFig3D, 'WindowButtonUpFcn', 'ihb_View3DButtonUpFcn');
	view([-150,32]);    % set initial camera position
	lighting phong;     % set lighting algorithm
	rotate3d on;        % enable 3D rotation

    %------------------------------------------------------------------------------
    % Smooth volume data
	%------------------------------------------------------------------------------
    if ihbdfl_bApplySmoothing == 1
        switch ihbdfl_smooth_filter
        case 'box'
            V = smooth3(cl.vTal, 'box', ihbdfl_smooth_size);
        case 'gaussian'
            V = smooth3(cl.vTal, 'gaussian', ihbdfl_smooth_size, ihbdfl_smooth_sd');
        end
    else
	    V = cl.vTal;
    end

    %------------------------------------------------------------------------------
    % Prepare volume data and create isosurface patch object
	%------------------------------------------------------------------------------
    % V is the volume of the cluster.

    FV = isosurface(cl.xTal, cl.yTal, cl.zTal, V, cl.hThreshold);
	hPatch = patch(FV);
	isonormals(cl.xTal, cl.yTal, cl.zTal, cl.vTal, hPatch);
	set(hPatch, 'Tag', 'ihb_surfaceCluster', 'FaceColor', ihbdfl_color_surface, 'EdgeColor', 'none', 'FaceAlpha', ihbdfl_transparency);
    set(hFig3D, 'Visible', 'on');

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function cut = knapsack_cheapcut(row)
% Dedicated easy cut for sum ai xi >= 0
% which appears effective on some models
if all(row(2:end)>=0) & row(1)<=0
    a = row(2:end);
    b = -row(1);
    [val,loc] = sort(a,'ascend');
    g = cumsum(val);
    m = max(find(g < -row(1)));
    ok = setdiff(loc,loc(1:m));
    b_cut = -1;
    a_cut = spalloc(1,length(val),0);
    a_cut(ok) = 1;
    cut = [b_cut a_cut];
else
    cut = [];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef prtOutlierRemovalNonFinite < prtOutlierRemoval
    % prtOutlierRemovalNonFinite  Removes missing data from a prtDataSet
    %
    %   OutRemove = prtOutlierRemovalNonFinite creates a pre-processing
    %   object that removes all missing data as outliers. Missing data are
    %   respresented as NaN's in the original data set.
    %
    %   A prtOutlierRemovalNonFinite object also inherits all properties and
    %   functions from the prtOutlierRemoval class.  For more information
    %   on how to control the behaviour of outlier removal objects, see the
    %   help for prtOutlierRemoval.
    %
    %   Example:
    %
    %   dataSet = prtDataGenUnimodal;               % Load a data Set
    %   outlier = prtDataSetClass([NaN NaN; inf inf],[1;1]);     % Create and insert
    %   dataSet = catObservations(dataSet,outlier); % an Outlier
    %
    %   % Create the prtOutlierRemoval object
    %   nonFiniteRemove = prtOutlierRemovalNonFinite('runMode','removeObservation');
    %
    %   nonFiniteRemove = nonFiniteRemove.train(dataSet);    % Train and run
    %   dataSetNew = nonFiniteRemove.run(dataSet);
    %
    %   See Also:  prtOutlierRemoval,
    %   prtOutlierRemovalMissingData ,prtOutlierRemovalNstd

    properties (SetAccess=private)
        % Required by prtAction
        name = 'Non Finite Data Outlier Removal';  % NonFinite Data Outlier Removal
        nameAbbreviation = 'NonFiniteDataRemove'   % NonFiniteDataRemove
    end

    methods

        % Allow for string, value pairs
        function Obj = prtOutlierRemovalNonFinite(varargin)

            %Need to check this with the string - setting the string in
            %prtOutlierRemoval should change this value...
            Obj.isCrossValidateValid = true;

            Obj = prtUtilAssignStringValuePairs(Obj,varargin{:});
        end
    end

    methods (Access = protected, Hidden = true)

        function Obj = trainAction(Obj,~)
            %Nothing to do
        end

        function outlierIndices = calculateOutlierIndices(Obj,DataSet) %#ok<MANU> This is OK - don't make it static.
            outlierIndices = ~isfinite(DataSet.getObservations);
        end

    end

end


text:
classdef TestContourFitting
    %TestContourFitting

    methods (Static)
        function test_contourSampling
            c = get_contour();

            % cell array
            c2 = cv.ContourFitting.contourSampling(c, 100);
            validateattributes(c2, {'cell'}, {'numel',100});

            % Nx2
            c = cat(1, c{:});
            c2 = cv.ContourFitting.contourSampling(c, 100);
            validateattributes(c2, {'numeric'}, {'size',[100 2]});

            % Nx1x2
            c = permute(c, [1 3 2]);
            c2 = cv.ContourFitting.contourSampling(c, 100);
            validateattributes(c2, {'numeric'}, {'size',[100 1 2]});
        end

        function test_fourierDescriptor
            c = get_contour();

            % cell array
            fd = cv.ContourFitting.fourierDescriptor(c, 'NumFD',30);
            validateattributes(fd, {'numeric'}, {});
            assert(length(fd)==30 && size(fd,3)==2);

            % Nx2
            c = cat(1, c{:});
            fd = cv.ContourFitting.fourierDescriptor(c, 'NumFD',30);
            validateattributes(fd, {'numeric'}, {});
            assert(length(fd)==30 && size(fd,3)==2);

            % Nx1x2
            c = permute(c, [1 3 2]);
            fd = cv.ContourFitting.fourierDescriptor(c, 'NumFD',30);
            validateattributes(fd, {'numeric'}, {});
            assert(length(fd)==30 && size(fd,3)==2);
        end

        function test_transformFD
            c = get_contour();
            t = [0 pi/2 1 10 20];  % [alpha, phi, s, Tx, Ty]

            % cell array
            c2 = cv.ContourFitting.transformFD(c, t, 'FD',false);
            validateattributes(c2, {'cell'}, {'nonempty'});

            % Nx2
            c = cat(1, c{:});
            c2 = cv.ContourFitting.transformFD(c, t, 'FD',false);
            validateattributes(c2, {'cell'}, {'nonempty'});

            % Nx1x2
            c = permute(c, [1 3 2]);
            c2 = cv.ContourFitting.transformFD(c, t, 'FD',false);
            validateattributes(c2, {'cell'}, {'nonempty'});
        end

        function test_transformFD_2
            c = get_contour();
            t = [0 pi/2 1 10 20];  % [alpha, phi, s, Tx, Ty]

            n = cv.getOptimalDFTSize(numel(c));
            cc = cv.ContourFitting.contourSampling(c, n);
            fd = cv.ContourFitting.fourierDescriptor(cc);

            fd2 = cv.ContourFitting.transformFD(fd, t, 'FD',true);
            validateattributes(fd2, {'numeric'}, {'nonempty'});
            assert(length(fd2)==length(fd) && size(fd2,3)==2);
        end

        function test_estimateTransformation
            c1 = get_contour();
            c2 = cellfun(@(pt) pt+[10 20], c1, 'UniformOutput',false);

            obj = cv.ContourFitting();
            [t, d] = obj.estimateTransformation(c1, c2, 'FD',false);
            validateattributes(t, {'double'}, {'vector', 'numel',5});
            validateattributes(d, {'numeric'}, {'scalar', 'nonnegative'});
        end

        function test_estimateTransformation_2
            c1 = get_contour();
            c2 = cellfun(@(pt) pt+[10 20], c1, 'UniformOutput',false);

            n = cv.getOptimalDFTSize(numel(c1));
            cc1 = cv.ContourFitting.contourSampling(c1, n);
            cc2 = cv.ContourFitting.contourSampling(c2, n);
            fd1 = cv.ContourFitting.fourierDescriptor(cc1);
            fd2 = cv.ContourFitting.fourierDescriptor(cc2);
            fd1 = reshape(fd1, [], 1, 2);
            fd2 = reshape(fd2, [], 1, 2);

            obj = cv.ContourFitting();
            [t, d] = obj.estimateTransformation(fd1, fd2, 'FD',true);
            validateattributes(t, {'double'}, {'vector', 'numel',5});
            validateattributes(d, {'numeric'}, {'scalar', 'nonnegative'});
        end
    end

end

function [c, img] = get_contour()
    im = fullfile(mexopencv.root(),'test','shape06.png');
    img = cv.imread(im, 'Grayscale',true);
    c = cv.findContours(img, 'Mode','List', 'Method','None');
    [~,idx] = max(cellfun(@numel,c));
    c = c{idx};
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function fileGrouper = DMImportRecursive(filenameOrRootDirectory, dicomLibrary, reporting)
    % DMImportRecursive. Loads metadata from a series of DICOM files
    %
    %     Syntax
    %     ------
    %
    %         file_grouper = DMLoadMetaDMImportRecursivedataFromDicomFiles(filenameOrRootDirectory, dicomLibrary, reporting)
    %
    %             file_grouper    a DMFileGrouper object containing the
    %                             metadata grouped into coherent sequences of images
    %
    %             filenameOrRootDirectory specify the location of the DICOM
    %                             files.
    %
    %             dicomLibrary    (Optional) An object implementing
    %                             DMDicomLibraryInterface, used to parse
    %                             the Dicom files. If no object is provided
    %                             then the default DMDicomLibrary is used
    %
    %             reporting       A CoreReporting or implementor of the same interface,
    %                             for error and progress reporting. Create a CoreReporting
    %                             with no arguments to hide all reporting. If no
    %                             reporting object is specified then a default
    %                             reporting object with progress dialog is
    %                             created
    %
    %
    %     Licence
    %     -------
    %     Part of DicoMat. https://github.com/tomdoel/dicomat
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the BSD 3-Clause license. Please see the file LICENSE for details.
    %

    % Create a reporting object if none was provided
    if nargin < 3 || isempty(reporting)
        reporting = CoreReportingDefault;
    end

    % Create a library object if none was provided
    if nargin < 2 || isempty(dicomLibrary)
        dicomLibrary = DMDicomLibrary.getLibrary;
    end

    % Show a progress dialog
    reporting.ShowProgress('Importing data');
    reporting.UpdateProgressValue(0);

    % The file grouper performs the sorting of image metadata
    fileGrouper = DMFileGrouper;

    dictionary = DMDicomDictionary.EssentialDictionaryWithoutPixelData;

    [importFolder, filenameOnly] = CoreDiskUtilities.GetFullFileParts(filenameOrRootDirectory);

    % Find out file type. 0=does not exist; 2=file; 7=directory
    exist_result = exist(filenameOrRootDirectory, 'file');

    if exist_result == 7
        % For a directory, we import every file
        ImportDirectoryRecursive(fileGrouper, dicomLibrary, importFolder, dictionary, reporting);

    elseif exist_result == 2
        % For a single file, we check if it is Dicom.
        if dicomLibrary.isdicom(filenameOrRootDirectory)
            % For a Dicom file we import the whole folder
            ImportDirectoryRecursive(fileGrouper, dicomLibrary, importFolder, dictionary, reporting);
        else
            % If this is not a Dicom image we exclude it from the set and warn the user
            reporting.ShowWarning('DMImportRecursive:NotADicomFile', ['The file ' filenameOrRootDirectory ' is not a DICOM file and will be removed from this series.']);
        end

    else
        reporting.Error('DMImportRecursive:FileDoesNotExist', ['The file or directory ' filenameOrRootDirectory ' does not exist.']);
    end

    reporting.CompleteProgress;
end

function ImportDirectoryRecursive(fileGrouper, dicomLibrary, directory, dictionary, reporting)
    directories_to_do = CoreStack(directory);

    while ~directories_to_do.IsEmpty
        current_dir = directories_to_do.Pop;

        reporting.UpdateProgressMessage(['Importing data from: ' current_dir]);

        if reporting.HasBeenCancelled
            reporting.Error(CoreReporting.CancelErrorId, 'User cancelled');
        end

        ImportFilesInDirectory(fileGrouper, dicomLibrary, current_dir, dictionary, reporting);
        next_dirs = CoreDiskUtilities.GetListOfDirectories(current_dir);
        for dir_to_add = next_dirs
            new_dir = fullfile(current_dir, dir_to_add{1});
            directories_to_do.Push(new_dir);
        end
    end
end

function ImportFilesInDirectory(fileGrouper, dicomLibrary, directory, dictionary, reporting)
    all_filenames = CoreTextUtilities.SortFilenames(CoreDiskUtilities.GetDirectoryFileList(directory, '*'));
    for filename = all_filenames
        fullFileName = fullfile(directory, filename{1});
        if ~strcmp(fullFileName, 'DICOMDIR')
            try
                if dicomLibrary.isdicom(fullFileName)
                    fileGrouper.AddFile(dicomLibrary.dicominfo(fullFileName, dictionary));
                else
                    % If this is not a Dicom image we exclude it from the set and warn the user
                    reporting.ShowWarning('DMImportRecursive:NotADicomFile', ['DMLoadMetadataFromDicomFiles: The file ' fullFileName ' is not a DICOM file and will be removed from this series.']);
                end
            catch ex
                reporting.ShowWarning('DMImportRecursive:ImportFileFailed', ['Failed to import file ' fullFileName ' due to error: ' ex.message]);
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit