text:
function [cfconds, cfshocks, cfblocks]=loadcfpan(endo,Units,panel,CFt,Fstartdate,Fenddate,Fperiods,pref)

% identify the number of endogenous variables
numendo=size(endo,1);
% identify the number of units
numunits=size(Units,1);
% initiate the cells
cfconds={};
cfshocks={};
cfblocks=[];

% recover the cfconds cell (for all types of conditional forecasts)

% load the data from Excel
[~, ~, strngs]=xlsread(pref.excelFile,'pan conditions');
% replace NaN entries by blanks
strngs(cellfun(@(x) any(isnan(x)),strngs))={[]};
% convert potential numeric entries into strings
strngs(cellfun(@isnumeric,strngs))=cellfun(@num2str,strngs(cellfun(@isnumeric,strngs)),'UniformOutput',0);
% identify the non-empty entries (pairs of rows and columns)
[nerows, neclmns]=find(~cellfun('isempty',strngs));
% count the number of such entries
neentries=size(nerows,1);
% all these entries contrain strings: fix them to correct potential user formatting errors
% loop over entries
for ii=1:neentries
    strngs{nerows(ii,1),neclmns(ii,1)}=bear.utils.fixstring(strngs{nerows(ii,1),neclmns(ii,1)});
end
% identify the rows coresponding respectively to the forecast start date and forecast end dates
temp=find(strcmp(strngs,Fstartdate));
% if empty, return an error
if isempty(temp)
    message=['conditional forecast error: forecast start date ' Fstartdate ' cannot be found. Please verify that the ''pan conditions'' sheet of the Excel data file is properly filled.'];
    error('bear:loadcfpan:FstartdateNotFound',message)
end
% otherwise, record
rows(1,1)=temp;
temp=find(strcmp(strngs,Fenddate));
% if empty, return an error
if isempty(temp)
    message=['conditional forecast error: forecast end date ' Fenddate ' cannot be found. Please verify that the ''pan conditions'' sheet of the Excel data file is properly filled.'];
    error('bear:loadcfpan:FenddateNotFound',message)
end
% otherwise, record
rows(2,1)=temp;

% now the treatment will differ according to the type of panel VAR model

% if the model is panel 2,3 or 4, different conditional forecasts will be estimated for each unit
% then cfconds becomes a cell with a number of pages equivalent to the number of units
if panel==2 || panel==3 || panel==4
    % create the cell cfconds
    cfconds=cell(Fperiods,numendo,numunits);
    % loop over units
    for ii=1:numunits
        % loop over endogenous
        for jj=1:numendo
            % find the column corresponding to the combination unit_variable
            [~,clmns]=find(strcmp(strngs,[Units{ii,1} '_' endo{jj,1}]));
            % if the string cannot be found, return an error
            if isempty(clmns)
                message=['conditional forecast error: endogenous variable ' endo{jj,1} ' cannot be found for unit ' Units{ii,1} '. Please verify that the ''pan conditions'' sheet of the Excel data file is properly filled.'];
                error('bear:loadcfpan:endoNotFound',message)
            end
            % else, loop over forecast periods
            for kk=1:Fperiods
                % fill the corresponding entry
                cfconds{kk,jj,ii}=str2double(strngs{rows(1,1)+kk-1,clmns});
            end
        end
    end

    % if the model is panel 5 or 6, a single model is estimated with all units, hence cfconds contains only one page
elseif panel==5 || panel==6
    % create the cell cfconds
    cfconds=cell(Fperiods,numendo*numunits);
    % initiate the column count
    clmncount=0;
    % loop over units
    for ii=1:numunits
        % loop over endogenous
        for jj=1:numendo
            % increment column count
            clmncount=clmncount+1;
            % find the column corresponding to the combination unit_variable
            [~,clmns]=find(strcmp(strngs,[Units{ii,1} '_' endo{jj,1}]));
            % if the string cannot be found, return an error
            if isempty(clmns)
                message=['conditional forecast error: endogenous variable ' endo{jj,1} ' cannot be found for unit ' Units{ii,1} '. Please verify that the ''pan conditions'' sheet of the Excel data file is properly filled.'];
                error('bear:loadcfpan:endoNotFound',message)
            end
            % else, loop over forecast periods
            for kk=1:Fperiods
                % fill the corresponding entry
                cfconds{kk,clmncount}=str2double(strngs{rows(1,1)+kk-1,clmns});
            end
        end
    end
end
% if no error is returned, record on Excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),strngs,'cf conditions','B2');
end

% recover the cfshocks cell (for shock-specific conditional forecasts only)

if CFt==2
    % load the data from Excel
    [~, ~, strngs]=xlsread(pref.excelFile,'pan shocks');
    % replace NaN entries by blanks
    strngs(cellfun(@(x) any(isnan(x)),strngs))={[]};
    % convert potential numeric entries into strings
    strngs(cellfun(@isnumeric,strngs))=cellfun(@num2str,strngs(cellfun(@isnumeric,strngs)),'UniformOutput',0);
    % identify the non-empty entries (pairs of rows and columns)
    [nerows, neclmns]=find(~cellfun('isempty',strngs));
    % count the number of such entries
    neentries=size(nerows,1);
    % all these entries contrain strings: fix them to correct potential user formatting errors
    % loop over entries
    for ii=1:neentries
        strngs{nerows(ii,1),neclmns(ii,1)}=bear.utils.fixstring(strngs{nerows(ii,1),neclmns(ii,1)});
    end
    % identify the rows coresponding respectively to the forecast start date and forecast end dates
    temp=find(strcmp(strngs,Fstartdate));
    % if empty, return an error
    if isempty(temp)
        message=['conditional forecast error: forecast start date ' Fstartdate ' cannot be found. Please verify that the ''pan shocks'' sheet of the Excel data file is properly filled.'];
        msgbox(message);
        error('programme termination: conditional forecast error');
    end
    % otherwise, record
    rows(1,1)=temp;
    temp=find(strcmp(strngs,Fenddate));
    % if empty, return an error
    if isempty(temp)
        message=['conditional forecast error: forecast end date ' Fenddate ' cannot be found. Please verify that the ''pan shocks'' sheet of the Excel data file is properly filled.'];
        msgbox(message);
        error('programme termination: conditional forecast error');
    end
    % otherwise, record
    rows(2,1)=temp;

    % now the treatment will differ according to the type of panel VAR model

    % if the model is panel 2,3 or 4, different conditional forecasts will be estimated for each unit
    % then cfshocks becomes a cell with a number of pages equivalent to the number of units
    if panel==2 || panel==3 || panel==4
        % create the cell cfshocks
        cfshocks=cell(Fperiods,numendo,numunits);
        % loop over units
        for ii=1:numunits
            % loop over endogenous
            for jj=1:numendo
                % find the column corresponding to the combination unit_variable
                [~,clmns]=find(strcmp(strngs,[Units{ii,1} '_' endo{jj,1}]));
                % if the string cannot be found, return an error
                if isempty(clmns)
                    message=['conditional forecast error: endogenous variable ' endo{jj,1} ' cannot be found for unit ' Units{ii,1} '. Please verify that the ''pan shocks'' sheet of the Excel data file is properly filled.'];
                    msgbox(message);
                    error('programme termination: conditional forecast error');
                end
                % else, loop over forecast periods
                for kk=1:Fperiods
                    % fill the corresponding entry
                    cfshocks{kk,jj,ii}=str2double(strngs{rows(1,1)+kk-1,clmns});
                end
            end
        end

        % if the model is panel 5 or 6, a single model is estimated with all units, hence cfshocks contains only one page
    elseif panel==5 || panel==6
        % create the cell cfshocks
        cfshocks=cell(Fperiods,numendo*numunits);
        % initiate the column count
        clmncount=0;
        % loop over units
        for ii=1:numunits
            % loop over endogenous
            for jj=1:numendo
                % increment column count
                clmncount=clmncount+1;
                % find the column corresponding to the combination unit_variable
                [~,clmns]=find(strcmp(strngs,[Units{ii,1} '_' endo{jj,1}]));
                % if the string cannot be found, return an error
                if isempty(clmns)
                    message=['conditional forecast error: endogenous variable ' endo{jj,1} ' cannot be found for unit ' Units{ii,1} '. Please verify that the ''pan shocks'' sheet of the Excel data file is properly filled.'];
                    msgbox(message);
                    error('programme termination: conditional forecast error');
                end
                % loop over forecast periods
                for kk=1:Fperiods
                    % fill the corresponding entry
                    cfshocks{kk,clmncount}=str2double(strngs{rows(1,1)+kk-1,clmns});
                end
            end
        end
    end

    % if no error is returned, record on Excel
    if pref.results==1
        bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),strngs,'cf shocks','B2');
    end
end

% recover the cfblocks matrix (for shock-specific conditional forecasts only)

if CFt==2
    % load the data from Excel
    [~, ~, strngs]=xlsread(pref.excelFile,'pan blocks');
    % replace NaN entries by blanks
    strngs(cellfun(@(x) any(isnan(x)),strngs))={[]};
    % convert potential numeric entries into strings
    strngs(cellfun(@isnumeric,strngs))=cellfun(@num2str,strngs(cellfun(@isnumeric,strngs)),'UniformOutput',0);
    % identify the non-empty entries (pairs of rows and columns)
    [nerows, neclmns]=find(~cellfun('isempty',strngs));
    % count the number of such entries
    neentries=size(nerows,1);
    % all these entries contrain strings: fix them to correct potential user formatting errors
    % loop over entries
    for ii=1:neentries
        strngs{nerows(ii,1),neclmns(ii,1)}=bear.utils.fixstring(strngs{nerows(ii,1),neclmns(ii,1)});
    end
    % identify the rows coresponding respectively to the forecast start date and forecast end dates
    % identify the rows coresponding respectively to the forecast start date and forecast end dates
    temp=find(strcmp(strngs,Fstartdate));
    % if empty, return an error
    if isempty(temp)
        message=['conditional forecast error: forecast start date ' Fstartdate ' cannot be found. Please verify that the ''pan shocks'' sheet of the Excel data file is properly filled.'];
        msgbox(message);
        error('programme termination: conditional forecast error');
    end
    % otherwise, record
    rows(1,1)=temp;
    temp=find(strcmp(strngs,Fenddate));
    % if empty, return an error
    if isempty(temp)
        message=['conditional forecast error: forecast end date ' Fenddate ' cannot be found. Please verify that the ''pan shocks'' sheet of the Excel data file is properly filled.'];
        msgbox(message);
        error('programme termination: conditional forecast error');
    end
    % otherwise, record
    rows(2,1)=temp;

    % now the treatment will differ according to the type of panel VAR model
    % if the model is panel 2,3 or 4, different conditional forecasts will be estimated for each unit
    % then cfblocks becomes a matrix with a number of pages equivalent to the number of units
    if panel==2 || panel==3 || panel==4
        % create the matrix cfshocks
        cfblocks=zeros(Fperiods,numendo,numunits);
        % loop over units
        for ii=1:numunits
            % loop over endogenous
            for jj=1:numendo
                % find the column corresponding to the combination unit_variable
                [~,clmns]=find(strcmp(strngs,[Units{ii,1} '_' endo{jj,1}]));
                % if the string cannot be found, return an error
                if isempty(clmns)
                    message=['conditional forecast error: endogenous variable ' endo{jj,1} ' cannot be found for unit ' Units{ii,1} '. Please verify that the ''pan blocks'' sheet of the Excel data file is properly filled.'];
                    msgbox(message);
                    error('programme termination: conditional forecast error');
                end
                % else, loop over forecast periods
                for kk=1:Fperiods
                    % recover the entry
                    temp=str2double(strngs{rows(1,1)+kk-1,clmns});
                    % if the entry is empty, ignore
                    if isempty(temp)
                        % if not empty, record in cfblocks
                    else
                        cfblocks(kk,jj,ii)=temp;
                    end
                end
            end
        end

        % if the model is panel 5 or 6, a single model is estimated with all units, hence cfblocks contains only one page
    elseif panel==5 || panel==6
        % create the cell cfblocks
        cfblocks=zeros(Fperiods,numendo*numunits);
        % initiate the column count
        clmncount=0;
        % loop over units
        for ii=1:numunits
            % loop over endogenous
            for jj=1:numendo
                % increment column count
                clmncount=clmncount+1;
                % find the column corresponding to the combination unit_variable
                [~,clmns]=find(strcmp(strngs,[Units{ii,1} '_' endo{jj,1}]));
                % if the string cannot be found, return an error
                if isempty(clmns)
                    message=['conditional forecast error: endogenous variable ' endo{jj,1} ' cannot be found for unit ' Units{ii,1} '. Please verify that the ''pan blocks'' sheet of the Excel data file is properly filled.'];
                    msgbox(message);
                    error('programme termination: conditional forecast error');
                end
                % else, loop over forecast periods
                for kk=1:Fperiods
                    % recover the entry
                    temp=str2double(strngs{rows(1,1)+kk-1,clmns});
                    % if the entry is empty, ignore
                    if isempty(temp)
                        % if not empty, record in cfblocks
                    else
                        cfblocks(kk,clmncount)=temp;
                    end
                end
            end
        end
    end

    % if no error is returned, record on Excel
    if pref.results==1
        bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),strngs,'cf blocks','B2');
    end
end

% final verification: check that cfconds is consistent with the other cells/matrices (i.e. cfshocks and cfblocks)
% the way this is done depends on the panel model
% if the panel model is 2, 3 or 4
if panel==2 || panel==3 || panel==4
    % if the type of conditional forecasts is shock-specific
    if CFt==2
        % loop over units
        for ii=1:numunits
            % check if there are conditions on unit ii
            temp1=cfconds(:,:,ii);
            nconds(ii,1)=numel(temp1(cellfun(@(x) any(~isempty(x)),temp1)));
            % if there are no conditions, ignore as there is nothing to check
            if nconds==0
                % if there are conditions
            else
                % identify the non-empty elements in cfconds
                [nerows1, neclmns1]=find(~cellfun('isempty',cfconds(:,:,ii)));
                % check that cfconds is consistent with cfshocks
                % identify the non-empty elements in cfshocks
                [nerows2, neclmns2]=find(~cellfun('isempty',cfshocks(:,:,ii)));
                % if the number of non-empty elements is not similar in the two cells, there is obvioulsy a problem: return an error
                if size(nerows1,1)~=size(nerows2,1)
                    message=['conditional forecast error: the conditions for unit ' Units{ii} ' seem to be inconsistent with the shocks. Please verify that the ''pan conditions'' and ''pan shocks'' sheets of the Excel data file are properly filled.'];
                    msgbox(message);
                    error('programme termination: conditional forecast error');
                    % else, if there a is a similar number of non-empty entries but they don't match, there is also a problem: return an error
                elseif (~isempty(find(nerows1-nerows2)) || ~isempty(find(neclmns1-neclmns2)))
                    message=['conditional forecast error: the conditions for unit ' Units{ii} ' seem to be inconsistent with the shocks. Please verify that the ''pan conditions'' and ''pan shocks'' sheets of the Excel data file are properly filled.'];
                    msgbox(message);
                    error('programme termination: conditional forecast error');
                end
                % similarly, check that cfconds is consistent with cfblocks
                % first turn cfblocks to cell
                temp=num2cell(cfblocks(:,:,ii));
                % switch zero entries to empty entries
                temp(cellfun(@(x) any(~(x)),temp))={[]};
                % identify the non-empty elements in cfblocks
                [nerows3, ~]=find(~cellfun('isempty',temp));
                % if the number of non-empty elements is not similar in the two cells, there is obvioulsy a problem: return an error
                if size(nerows1,1)~=size(nerows3,1)
                    message=['conditional forecast error: the conditions for unit ' Units{ii} ' seem to be inconsistent with the blocks. Please verify that the ''pan conditions'' and ''pan blocks'' sheets of the Excel data file are properly filled.'];
                    msgbox(message);
                    error('programme termination: conditional forecast error');
                    % else, if there a is a similar number of non-empty entries but they don't match, there is also a problem: return an error
                elseif (~isempty(find(nerows1-nerows2)) || ~isempty(find(neclmns1-neclmns2)))
                    message=['conditional forecast error: the conditions for unit ' Units{ii} ' seem to be inconsistent with the blocks. Please verify that the ''pan conditions'' and ''pan blocks'' sheets of the Excel data file are properly filled.'];
                    msgbox(message);
                    error('programme termination: conditional forecast error');
                end
            end
        end
    end

    %else if the panel model is 5 or 6
elseif panel==5 || panel==6
    % if the type of conditional forecasts is shock-specific
    if CFt==2
        % identify the non-empty elements in cfconds
        [nerows1, neclmns1]=find(~cellfun('isempty',cfconds));
        % check that cfconds is consistent with cfshocks
        % identify the non-empty elements in cfshocks
        [nerows2, neclmns2]=find(~cellfun('isempty',cfshocks));
        % if the number of non-empty elements is not similar in the two cells, there is obvioulsy a problem: return an error
        if size(nerows1,1)~=size(nerows2,1)
            message=['conditional forecast error: the conditions seem to be inconsistent with the shocks. Please verify that the ''pan conditions'' and ''pan shocks'' sheets of the Excel data file is properly filled.'];
            msgbox(message);
            error('programme termination: conditional forecast error');
            % else, if there a is a similar number of non-empty entries but they don't match, there is also a problem: return an error
        elseif (~isempty(find(nerows1-nerows2)) || ~isempty(find(neclmns1-neclmns2)))
            message=['conditional forecast error: the conditions seem to be inconsistent with the shocks. Please verify that the ''pan conditions'' and ''pan shocks'' sheets of the Excel data file is properly filled.'];
            msgbox(message);
            error('programme termination: conditional forecast error');
        end
        % similarly, check that cfconds is consistent with cfblocks
        % first turn cfblocks to cell
        temp=num2cell(cfblocks);
        % switch zero entries to empty entries
        temp(cellfun(@(x) any(~(x)),temp))={[]};
        % identify the non-empty elements in cfblocks
        [nerows3, neclmns3]=find(~cellfun('isempty',temp));
        % if the number of non-empty elements is not similar in the two cells, there is obvioulsy a problem: return an error
        if size(nerows1,1)~=size(nerows3,1)
            message=['conditional forecast error: the conditions seem to be inconsistent with the blocks. Please verify that the ''pan conditions'' and ''pan blocks'' sheets of the Excel data file is properly filled.'];
            msgbox(message);
            error('programme termination: conditional forecast error');
            % else, if there a is a similar number of non-empty entries but they don't match, there is also a problem: return an error
        elseif (~isempty(find(nerows1-nerows3)) || ~isempty(find(neclmns1-neclmns3)))
            message=['conditional forecast error: the conditions seem to be inconsistent with the blocks. Please verify that the ''pan conditions'' and ''pan blocks'' sheets of the Excel data file is properly filled.'];
            msgbox(message);
            error('programme termination: conditional forecast error');
        end
    end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function X = matrixReadFromFID(FID, varargin)

% MATRIXREADFROMFID Read a matrix from an FID.
% FORMAT
% DESC reads a matrix from an FID.
% ARG FID: the file ID to read the matrix from.
% RETURN X : the returned matrix read from the file.
%
% COPYRIGHT : Neil D. Lawrence, 2008
%
% SEEALSO : modelReadFromFID, doublematrixReadFromFID

% MLTOOLS

modelType = readStringFromFID(FID, 'type');
feval = str2func([modelType 'ReadFromFID']);
X = feval(FID, varargin{:});

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Euclidean projection on the nonnegative orthant
% Section 8.1.1, Boyd & Vandenberghe "Convex Optimization"
% Joelle Skaf - 10/07/05
%
% The projection of x0 on the proper cone K = R+^n is given by
%           minimize || x - x0 ||^2
%               s.t.    x >= 0
% It is also given by: P_K(x0)_k = max{x0_k,0}

% Input data
randn('seed',0);
n  = 10;
x0 = randn(n,1);

fprintf(1,'Computing the analytical solution...');

% Analytical solution
pk_x0 = max(x0,0);

fprintf(1,'Done! \n');

% Solution via CVX
fprintf(1,'Computing the solution via a QP...');

cvx_begin quiet
    variable x(n)
    minimize ( norm(x - x0) )
    x >= 0; %#ok
cvx_end

fprintf(1,'Done! \n');

% Verification
disp('-----------------------------------------------------------------');
disp('Verifying that the analytical solution and the solution obtained via QP are equal: ');
disp([pk_x0 x])
disp('They should be equal!');

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [A,twt]=loadrd3(fname,varargin)
% loads an rd3 file... (requires the .rad file to be present)
% syntax: [A,twt]=loadrd3('prof4',[starttrace],[tracecount],[options])
%
% options: 'stack'... return the mean of the profile
%
% Aslak Grinsted 2002
%

% load details from .RAD file

fname=strrep(lower(fname),'.rd3','');

numargcount=size(varargin,2);
tracecount=inf;
starttrace=1;
stack=0;

l=1;
numargcount=size(varargin,2); %how many numeric arguments are there before the settings?
while (l<=size(varargin,2))
   a=varargin{l};
   if ischar(a)
      v=[];
      if (l<numargcount), numargcount=l-1; end
      %if ((l+1)<=size(varargin,2)), v=varargin{l+1}; end
      switch lower(a)
      case {'stack','s'}
          stack=1;
      end
   end
   l=l+1;
end
switch numargcount
case 2
    starttrace=varargin{1};
    tracecount=varargin{2};
end

infos=rd3info(fname,'scans','samples','timewindow');
sz=infos{2};
timewindow=infos{3};

twt=((1:sz)'-1)*timewindow*1e-9/sz;

%loads the rd3 file
fid=fopen([fname '.rd3'],'r');
if (starttrace>1)
   fseek(fid,(starttrace-1)*sz*2,0);
end

if (stack==0)
    A=fread(fid,[sz tracecount],'short');
else
    blocksize=1+fix(1000000/sz); %number of traces to read at a time...

    A=zeros(sz,1);
    tracesread=0;
    while ((~feof(fid))&(tracesread<tracecount))
        B=fread(fid,[sz blocksize],'short');
        A=A+sum(B,2);
        tracesread=tracesread+size(B,2);
    end
    A=A/tracesread;
end
fclose(fid);
%A=reshape(A,sz,size(A,1)/sz);

author: Sable

repo: mcbench-benchmarks


text:
function test_suite = test_surficial_io_gifti
% tests for GIFTI input/output
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_gifti_dataset_io()
    if cosmo_skip_test_if_no_external('gifti')
        return;
    end

    ds=cosmo_synthetic_dataset('ntargets',2,'nchunks',1,...
                        'type','surface');
    ds.a.fdim.values{1}=1+[1 4 7 8 5 3];

    encodings=get_gifti_encodings();
    orders=get_indexing_orders();

    fn=cosmo_make_temp_filename('_tmp','.gii');
    cleaner=onCleanup(@()delete(fn));

    % only test subset
    ntest=2;

    n_encodings=numel(encodings);
    n_orders=numel(orders);

    rp=randperm(n_encodings*n_orders);
    test_idxs=rp(1:ntest);

    counter=0;
    for k=1:numel(encodings)
        encoding=encodings{k};
        for j=1:numel(orders)
            counter=counter+1;
            if ~any(test_idxs==counter)
                continue;
            end

            order=orders{j};

            cosmo_map2surface(ds,fn,'encoding',encoding);
            ds2=cosmo_surface_dataset(fn);
            assert_dataset_equal(ds,ds2);

            g_str=get_gifti_string(encoding, order);
            fid=fopen(fn,'w');
            fprintf(fid,'%s',g_str);
            fclose(fid);

            ds3=cosmo_surface_dataset(fn);
            assert_dataset_equal(ds,ds3);

            g=gifti(fn);
            assertElementsAlmostEqual(double(g.cdata)',ds.samples,...
                                            'absolute',1e-5);
            node_indices=ds.a.fdim.values{1}(ds.fa.node_indices);
            assertEqual(double(g.indices)',node_indices);

            g2=cosmo_map2surface(ds,'-gii');
            assertElementsAlmostEqual(g.cdata,g2.cdata,...
                                            'absolute',1e-5);
            assertEqual(g.indices,g2.indices);
        end
    end

function assert_dataset_equal(x,y)
    assertElementsAlmostEqual(x.samples,y.samples,'absolute',1e-5);
    assertEqual(sort(fieldnames(x)),sort(fieldnames(y)));
    assertEqual(x.fa,y.fa);
    assertEqual(x.a,y.a);

function string=get_gifti_string(encoding, indexing_order)
    header=get_header();
    body=get_data_arrays(encoding, indexing_order);
    footer=get_footer();

    string=[header body footer];

function string=get_header()
    string=sprintf(['<?xml version="1.0" encoding="UTF-8"?>\n'...
                '<!DOCTYPE GIFTI SYSTEM "http://gifti.projects.nitr'...
                        'c.org/gifti.dtd">\n'...
                '<GIFTI Version="1.0"  NumberOfDataArrays="3">\n'...
                '   <MetaData>\n'...
                '      <MD>\n'...
                '         <Name><![CDATA[UniqueID]]></Name>\n'...
                '         <Value><![CDATA[XYZ_Uj6iYZUAV3HKbGirxcn-c'...
                        'A]]></Value>\n'...
                '      </MD>\n'...
                '   </MetaData>\n'...
                '   <LabelTable/>\n'...
                ]);

function string=get_footer()
    string='</GIFTI>';

function encodings=get_gifti_encodings()
    encodings={'ASCII','Base64Binary','GZipBase64Binary'};

function orders=get_indexing_orders()
    orders={'ColumnMajorOrder','RowMajorOrder'};

function string=get_data_array_helper(data_str, encoding, data_type, ...
                                               indexing_order, intent)
    check_one('encoding',encoding,...
                    get_gifti_encodings());
    check_one('data_type',data_type,...
                    {'NIFTI_TYPE_INT32','NIFTI_TYPE_FLOAT32'});
    check_one('indexing_order',indexing_order,...
                    get_indexing_orders());

    string=sprintf([...
                '   <DataArray Intent="%s"\n'...
                '              DataType="%s"\n'...
                '              ArrayIndexingOrder="%s"\n'...
                '              Dimensionality="1"\n'...
                '              Dim0="6"\n'...
                '              Encoding="%s"\n'...
                '              Endian="LittleEndian"\n'...
                '              ExternalFileName=""\n'...
                '              ExternalFileOffset="">\n'...
                '      <MetaData>\n'...
                '         <MD>\n'...
                '            <Name><![CDATA[Name]]></Name>\n'...
                '            <Value><![CDATA[numeric]]></Value>\n'...
                '         </MD>\n'...
                '      </MetaData>\n'...
                '      <Data>\n%s'...
                '      </Data>\n'...
                '   </DataArray>\n'...
                ], intent, data_type, indexing_order, encoding, data_str);

function string=get_data_arrays(encoding, indexing_order)
    switch encoding
        case 'ASCII'
            data={sprintf([ '        2.031686\n'...
                            '        -3.684948\n'...
                            '        -1.050449\n'...
                            '        1.349442\n'...
                            '        -0.261723\n'...
                            '        -0.203955\n'...
                  ]),...
                  sprintf([ '        0.583806\n'...
                            '        1.723501\n'...
                            '        -1.326491\n'...
                            '        -0.397334\n'...
                            '        2.338657\n'...
                            '        0.482345\n'...
                  ]),...
                  sprintf([ '         1 \n'...
                            '         4  \n'...
                            '         7 \n'...
                            '         8 \n'...
                            '         5 \n'...
                            '         3 \n'...
                  ])};
        case 'Base64Binary'
            data={'JQcCQDDWa8AddYa/hLqsP48Ahr6U2VC+',...
                  'T3QVP66b3D91yqm/XW/Lvo6sFUDt9fY+',...
                  'AQAAAAQAAAAHAAAACAAAAAUAAAADAAAA',...
                  };
        case 'GZipBase64Binary'
            data={  'eAEBGADn/yUHAkAw1mvAHXWGv4S6rD+PAIa+lNlQvnPfCu4=',...
                    'eAEBGADn/090FT+um9w/dcqpv11vy76OrBVA7fX2PpcYDR0=',...
                    'eAEBGADn/wEAAAAEAAAABwAAAAgAAAAFAAAAAwAAAAGEAB0='...
                    };
    end

    narrays=numel(data);
    arrays_cell=cell(1,narrays);

    data_types={'NIFTI_TYPE_FLOAT32',...
                'NIFTI_TYPE_FLOAT32',...
                'NIFTI_TYPE_INT32',...
                };
    intents={'NIFTI_INTENT_NONE',...
                'NIFTI_INTENT_NONE',...
                'NIFTI_INTENT_NODE_INDEX'};

    assert(numel(data_types)==narrays);
    assert(numel(intents)==narrays);

    for k=1:narrays
        arrays_cell{k}=get_data_array_helper(data{k}, encoding, ...
                                            data_types{k}, ...
                                            indexing_order, intents{k});
    end

    string=cat(2,arrays_cell{:});

function check_one(label, value, allowed)
    if isempty(strmatch(value, allowed))
        error('illegal value for %s: %s', label, value);
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
% Plot the per-pixel KL divergence when one 2D distribution is used to
% approximate another.

% Zoya Bylinskii, April 2016
% linked to: "What do different evaluation metrics tell us about saliency models?"

% Note: if comparing multiple saliency maps, compute the resMap for each of
% the saliency maps using this function, but then normalize them jointly
% before plotting (to make range of values comparable across plots).

function [resMap,figtitle] = visualize_KL(salMap, fixMap, toplot)
% salMap is the saliency model prediction (distribution)
% fixMap is the ground truth fixation map (distribution)
% if only plotting one map, set toplot=1; but if comparing multiple maps,
% see note above

if nargin < 3
    toplot = 1;
end

figtitle = 'KL divergence';

% create the red color map, brighter red corresponds to higher KL
% divergence values
G = linspace(0,1,20)';
fpmap = horzcat(G,zeros(size(G)),zeros(size(G)));

% format, resize, and normalize maps
map1 = im2double(fixMap);
map2 = im2double(imresize(salMap, size(fixMap)));
map1 = map1/sum(map1(:));
map2 = map2/sum(map2(:));

% compute per-pixel KL divergence
resMap = map1 .* log(eps + map1./(map2+eps)); % map 1 is the target distribution, map 2 the prediction

% plot the visualization
if toplot
    figure;
    resMap_norm = resMap/max(resMap(:));
    imshow(resMap_norm);
    colormap(fpmap); colorbar;
    title(figtitle,'fontsize',14);
end

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function ms=zz_flood_fill_movie(r,tol)
global hi I
% this just make move
% flood fill, scan line algoritm
% ms=flood_fill(RR,GG,BB,r,tol)
% % ms - pixels numbers that flooded
% numbering: number=y+szy*(x-1), x y - pixels coordinaties
% szx szy - image size
% I - RGB image
% r- first point of selection
% tol - tolerance

clear F;
fc=1;

R=int32(I(:,:,1));
G=int32(I(:,:,2));
B=int32(I(:,:,3));

[szy szx]=size(R); % imae size

% stek, where seed will be stored:
stm=10000;
st=zeros(stm,2,'int32');
st(1,1)=r(1);
st(1,2)=r(2); % r - is start seed
stL=1; % stack length

hp=plot(r(1),r(2),'.r');
set(hp,'MarkerSize',10);

set(hp,'XData',st(1:stL,1),'YData',st(1:stL,2));
drawnow;

% found pixel store:
ms0m=1e6; % margin
ms0=zeros(ms0m,1,'int32'); % predifined array to increase speed
ms0L=0; % 0 points initially

% initial point color:
R0=R(r(2),r(1));
G0=G(r(2),r(1));
B0=B(r(2),r(1));

% to pixel number
tn=@(xx,yy) yy+szy*(xx-1);

while true;
    % get seed from stack:
    xt=st(stL,1);
    yt=st(stL,2);
    stL=stL-1;

    % line:
    % to right
    sku=false; % seed key, true if seed added, up
    skd=false; % same for down
    sku1=false; % this key is need to prewnt extra seed when move to left, up
    skd1=false; % same for left
    for xtt=xt:szx
        drawnow;
        F(fc)=getframe(gcf);
        fc=fc+1;

        if max(abs([(R(yt,xtt)-R0), (G(yt,xtt)-G0), (B(yt,xtt)-B0)]))<=tol
            % add pixel
            ms0L=ms0L+1;
            ms0(ms0L)=tn(xtt,yt);
            I(yt,xtt,1)=0;
            set(hi,'Cdata',I);
        else
            break;
        end

        % try to add seed up:
        if yt~=szy
            if max(abs([(R(yt+1,xtt)-R0), (G(yt+1,xtt)-G0), (B(yt+1,xtt)-B0)]))<=tol
                if ~sku
                    if all(tn(xtt,yt+1)~=ms0(1:ms0L)) % if free space
                        % add to stack
                        stL=stL+1;
                        st(stL,1)=xtt;
                        st(stL,2)=yt+1;
                        sku=true;
                    end
                end
            else
                sku=false;
            end
            if xtt==xt
                sku1=sku; % memorize, will be used when to left
            end
        end

        % try to add down
        if yt~=1
            if max(abs([(R(yt-1,xtt)-R0), (G(yt-1,xtt)-G0), (B(yt-1,xtt)-B0)]))<=tol
                if ~skd
                    if all(tn(xtt,yt-1)~=ms0(1:ms0L)) % if free space
                        % add to stack
                        stL=stL+1;
                        st(stL,1)=xtt;
                        st(stL,2)=yt-1;
                        skd=true;
                    end
                end
            else
                skd=false;
            end
            if xtt==xt
                skd1=skd; % memorize, will be used when to left
            end
        end
    end

    % to left
    %sku=false; % seed key, true if seed added
    %skd=false;
    sku=sku1;
    skd=skd1;
    if xt~=1
        for xtt=(xt-1):-1:1
            drawnow;
            F(fc)=getframe(gcf);
            fc=fc+1;

            if max(abs([(R(yt,xtt)-R0), (G(yt,xtt)-G0), (B(yt,xtt)-B0)]))<=tol
                % add pixel
                ms0L=ms0L+1;
                ms0(ms0L)=tn(xtt,yt);
                I(yt,xtt,1)=0;
                set(hi,'Cdata',I);
            else
                break;
            end

            % try to add seed up:
            if yt~=szy
                if max(abs([(R(yt+1,xtt)-R0), (G(yt+1,xtt)-G0), (B(yt+1,xtt)-B0)]))<=tol
                    if ~sku
                        if all(tn(xtt,yt+1)~=ms0(1:ms0L)) % if free space
                            % add to stack
                            stL=stL+1;
                            st(stL,1)=xtt;
                            st(stL,2)=yt+1;
                            sku=true;
                        end
                    end
                else
                    sku=false;
                end
            end

            % try to add down
            if yt~=1
                if max(abs([(R(yt-1,xtt)-R0), (G(yt-1,xtt)-G0), (B(yt-1,xtt)-B0)]))<=tol
                    if ~skd
                        if all(tn(xtt,yt-1)~=ms0(1:ms0L)) % if free space
                            % add to stack
                            stL=stL+1;
                            st(stL,1)=xtt;
                            st(stL,2)=yt-1;
                            skd=true;
                        end
                    end
                else
                    skd=false;
                end
            end
        end
    end

    set(hp,'XData',st(1:stL,1),'YData',st(1:stL,2));
    drawnow;
    %pause(1);

    if stL==0 % no more seed
        break; % stop
    end

end

ms=ms0(1:ms0L);

% save as avi file:
movie2avi(F(1:4:end),'flood_fill_movie','fps',20,'compression','Cinepak');

% save as animated gif
[im,map] = rgb2ind(F(20).cdata,256,'nodither');
im(1,1,1,20) = 0;
k=1;
for fc1=1:4:(fc-1)
    im(:,:,1,k) = rgb2ind(F(fc1).cdata,map,'nodither');
    k=k+1;
end
imwrite(im,map,'flood_fill_gif.gif','DelayTime',0,'LoopCount',inf) %g443800


text:
function obj = setFuncs(obj, funcs)
    % This function specify the functions (name) that calculates the
    % the function
    %
    % Parameters:
    % funcs:
    %
    % Required fields of Funcs:
    %   Func: a string of the function that computes the
    %   function value @type char
    %   Jac: a string of the function that computes the
    %   function Jacobian @type char
    %   JacStruct: a string of the function that computes the
    %   sparsity structure of function Jacobian @type char
    %   Hess: a string of the function that computes the
    %   function Hessian @type char
    %   HessStruct: a string of the function that computes the
    %   sparsity structure of function Hessian @type char
    % @type struct

    assert(isstruct(funcs),...
        'The input argument must be a struct.');

    assert(all(isfield(funcs, {'Func', 'Jac', 'JacStruct', 'Hess', 'HessStruct'})),...
        ['Required fields missing. \n',...
        'Following fields must be present: %s'], ...
        implode({'Func', 'Jac', 'JacStruct', 'Hess', 'HessStruct'},', '));

    obj.Funcs = funcs;

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function ids = database_ids( database, gt_set )
    index_file = fullfile(root_dir,'datasets',database, 'gt_sets',[gt_set '.txt']);
    fileID = fopen(index_file);
    ids = textscan(fileID, '%s');
    ids = ids{1};
    fclose(fileID);
end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
function call_price=european_call_div(S, K, r, sigma, time_to_maturity, dividend_times, dividend_amounts)

%--------------------------------------------------------------------------
%
% DESCRIPTION:
%
% European option price with an underlying stock paying a discrete-time
% dividend
%
%
% Reference:
%
% John Hull, "Options, Futures and other Derivative Securities",
% Prentice-Hall, second edition, 1993.
%
%--------------------------------------------------------------------------
%
% INPUTS:
%
%  S:                   spot price
%  K:                   exercice price
%  r:                   interest rate
%  y:                   continous payout
%  sigma:               volatility
%  time_to_maturity:    time to maturity
%  dividend_times:      periods of dividend payment
%  dividend_amounts:    amounts of dividend payment
%
%--------------------------------------------------------------------------
%
% OUTPUT:
%
% call_price: price of a call option
%
%--------------------------------------------------------------------------
%
% Author:  Paolo Z., February 2012
%
%--------------------------------------------------------------------------

adjusted_S = S;

for ( i=1:max(size(dividend_times)) )
    if (dividend_times(i)<=time_to_maturity)
        adjusted_S = adjusted_S-dividend_amounts(i)*exp(-r*dividend_times(i));
    end
end

call_price = bs_european_call(adjusted_S,K,r,sigma,time_to_maturity);


text:
function img2new = imSameHeight(img1, img2)
%
%
%        img2new = imSameHeight(img1, img2)
%
%
%        Input:
%           -img1: target size image
%           -img2: source image to adjust with the same height of image 1
%
%        Output:
%           -img2new: source image with the same height of img1
%
%     Copyright (C) 2011  Francesco Banterle
%
%     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/>.
%

[r1,c1,col1] = size(img1);
[r2,c2,col2] = size(img2);

c2new = round((c2*r1)/r2);
img2new = imresize(img2, [r1,c2new] , 'bilinear');

end

author: banterle


text:
function [mu,P,w,cost]=kMeanspp(z,K,maxIter)
%%KMEANSPP Run the k-means++ algorithm for clustering a set of more than K
%          data points into k clusters. The k-means++ algorithm is the
%          same as the general k-means algorithm, except a different
%          initialization routine is used.
%
%INPUTS: z A zDim X numPoints set of vectors that are to be clustered.
%        K The number of clusters to form. K>=numPoints.
%  maxIter The maximum number of iterations to perform for clustering. If
%          omitted, maxIter=1000;
%
%OUTPUTS: mu A zDim X K set of the cluster means.
%          P A zDim X zDim X K set of sample covariance matrices for each
%            of the k clusters, where mu(:,n) is the mean of the nth
%            cluster.
%          w A KX1 vector of weights such that w(n) is the fraction of the
%            total original points assigned to the nth cluster.
%       cost The cost of the k-means assignment. This is the sum of the
%            squared distances between the points assigned to a cluster and
%            the cluster mean.
%
%The k-means algorithm is a suboptimal algorithm that tries to find a set
%of k-means such that the sum of the squared distances from the points to
%the closest means is minimized. The implementation here uses a better
%initialization the k-means++ initialization from [1], which improves the
%probability that the algorithm will converge to a good solution.
%
%Note that the k-means algorithm used a randomized initialization, so one
%will not always get the same results when the function is run twice on the
%same data.
%
%REFERENCES:
%[1] D. Arthur and S. Vassilvitskii, "k-means++: The advantages of careful
%    seeding," in Proceedings of the Eighteenth Annual ACM-SIAM Symposium
%    on Discrete Algorithms, New Orleans, LA, Jan. 2007, pp. 1027-1035.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3)
   maxIter=1000;
end

%This uses the k-mean++ initalization algorithm on the given data.
%K is the number of clusters to form.

zDim=size(z,1);
numPoints=size(z,2);

mu=zeros(zDim,K);

%Use k-mean++ initialization.
%The first center is chosen randomly.
mu(:,1)=z(:,randi([1;numPoints]));
diff=mu(:,1)*ones(1,numPoints)-z;
nearestMuDist=sum(diff.*diff,1);%The squared distances from this point.

for k=2:K
    %The points are chosen with a probability equal to the ratio of the
    %squared distance to the sum of all distances.
    PMF=nearestMuDist/sum(nearestMuDist);%The PMF from which we will draw.
    CMF=cumsum(PMF);

    %We have to find the first element in the CMF that is greater than a
    %random draw; that will be the next center.
    idx=sum(CMF<rand(1))+1;
    %idx=randi([1;numPoints]);%For the regular k-means.
    %Set the center and update the nesrest center distances.
    mu(:,k)=z(:,idx);
    diff=mu(:,k)*ones(1,numPoints)-z;
    nearestMuDist=min(nearestMuDist,sum(diff.*diff,1));
end

%The k-means refinement.
muPrev=ones(zDim,K)*Inf;
numIter=0;
while(sum(sum(muPrev~=mu))~=0&&numIter<maxIter)
    %Calculate the distances between the points and the current set of
    %means.
    minCosts=ones(1,numPoints)*Inf;
    selClust=zeros(1,numPoints);
    for k=1:K
        diff=bsxfun(@minus,z,mu(:,k));%mu(:,k)*ones(1,numPoints)-z;
        dist=sum(diff.*diff,1);%The squared distances from this center.

        sel=dist<minCosts;
        minCosts(sel)=dist(sel);
        %Assign those points to this cluster if it has the lowest cost.
        selClust(sel)=k;
    end

    %Now, all of the points are assigned to clusters, and we can calculate
    %a new set of means.
    muPrev=mu;
    for k=1:K
        mu(:,k)=mean(z(:,selClust==k),2);
    end
    numIter=numIter+1;
end

%The clustering is done. We will now also return the weights and covariance
%estimates. The weights will be the fraction of the total points assigned
%to a given cluster. The covariance will be the sample covariance.
if(nargout>1)
    w=zeros(K,1);
    P=zeros(zDim,zDim,K);
    for k=1:K
        sel=selClust==k;
        w(k)=sum(sel);

        zSel=z(:,sel);

        diff=bsxfun(@minus,zSel,mu(:,k));
        P(:,:,k)=diff*diff'/w(k);
    end
    w=w/sum(w);

    cost=sum(minCosts)/numPoints;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


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

    n = length(y); % number of training examples

    % cost: J, this time with the penalty for the magnitude of theta
    J = -1/n * sum(...
        y      .* log(sigmoid(X * theta)) + ...
        (1 - y) .* log(1 - sigmoid(X * theta)) ...
    ) + lambda / (2 * n) * sum(theta(2:end) .* theta(2:end));

    % gradient: compute as the derivative of the cost function
    grad = 1/n * X' * ((sigmoid(X * theta)) - y) + theta * lambda / n;

    % we do not regularize the constant offset term,
    % undo the gradient penalization for the constant x_0 term
    grad(1) = grad(1) - lambda / n * theta(1);
end

author: worldveil

repo: coursera-ml


text:
function [ valLoss ] = HuberLoss( vX, paramDelta )
% ----------------------------------------------------------------------------------------------- %
% [ valLoss ] = HuberLoss( vX, paramDelta )
%   Applies the Huber Loss to the input vector.
% Input:
%   - vY            -   Input Vector.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - ballRadius    -   Ball Radius.
%                       Sets the Radius of the L1 Ball. For Unit L1 Ball
%                       set to 1.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: (0, inf).
% Output:
%   - vX            -   Output Vector.
%                       The projection of the Input Vector onto the L1
%                       Ball.
%                       Structure: Vector (Column).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  Huber Loss (Wikipedia) - https://en.wikipedia.org/wiki/Huber_loss.
% Remarks:
%   1.  a
% TODO:
%   1.  U.
% Release Notes:
%   -   1.0.000     21/03/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

numElements = length(vX);

valLoss = 0;

for ii = 1:numElements
    if(abs(vX(ii)) <= paramDelta)
        valLoss = valLoss + (0.5 * vX(ii) * vX(ii));
    else
        valLoss = valLoss + (paramDelta * (abs(vX(ii)) - (0.5 * paramDelta)));
    end

end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function window = mandelbrotViewer()
%mandelbrotViewer  view and explore the Mandelbrot set using a GPU
%
%   mandelbrotViewer() opens a MATLAB figure window showing the Mandelbrot
%   set. Use the usual zoom and pan controls from the figure window toolbar
%   to navigate around and explore, or click "animate" to see a pre-defined
%   path through the set. You can move back to the initial view at any time
%   by clicking the "reset" button or add the current view to the animation
%   list using "add".
%
%   The control panel can be hidden using the right-hand toolbar button.
%
%   A selector allows a choice of four ways to calculate each frame:
%
%   1. CPU: All calculations are performed by MATLAB on the host CPU. The
%   algorithm is fully vectorized and avoids indexing to give an efficient
%   calculation. Even so, this may take a few seconds to calculate each
%   frame.
%
%   2. GPU (simple): The same algorithm as (1) is used but the input
%   coordinates are switched to being stored on the GPU. This causes
%   MATLAB to operate on the resulting data array on the GPU with no other
%   code change. This gives some speedup at virtually no coding cost.
%
%   3. GPU Arrayfun: Now we change the code so that instead of many
%   operations running on the full data matrix, we now specify a single
%   "calculateElement" operation and run it for each element. Roughly
%   speaking, MATLAB translates one element into one thread on the GPU,
%   giving huge speedups for large arrays. Note that all code is still
%   written in MATLAB and the user needs no knowledge of how GPU kernels
%   are constructed and executed.
%
%   4. CUDAKernel: Taking things to the limit, we now hand-craft the
%   element-wise algorithm used in (3) in CUDA C++. The resulting kernel is
%   called from MATLAB using the CUDAKernel system, requiring the user to
%   specify the thread and block arrangements to use.
%
%   Note that version 3 gets us most of the speedup achieved by the
%   hand-crafted CUDA (version 4) but without any need to leave the comfort
%   of MATLAB!
%
%   See also:  gpuArray, mandelbrotViewerProcessElement

%   Copyright 2010-2011 The Mathworks, Inc.

% Check that we are running in R2011a or above and have a GPU
matlabVersionCheck();
gpuCheck();

% Define some global (to this file) data structures so that they can be
% used by all the helper functions.
data = createData();
gui = createGUI();
% Make sure the image is updated now that the window is onscreen
redraw();

% Return the window handle if requested
if nargout
    window = gui.Window;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Here are the four different Mandelbrot Set computations

    function logCount = computeMandelbrotCPU( xlim, numx, ylim, numy, maxIters )
        % Create a view of the Mandelbrot set using only the CPU.
        % This is the base-line version of the algorithm that is adapted to
        % run on the GPU in different ways in the functions below.

        % Create the input arrays
        escapeRadius2 = 400; % Square of escape radius
        x = linspace( xlim(1),  xlim(2), numx );
        y = linspace( ylim(1),  ylim(2), numy );
        [x0,y0] = meshgrid(x, y);
        count = zeros( size( x0 ) );
        z0 = complex( x0, y0 );
        z = z0;

        % Calculate
        for n = 0:maxIters
            inside = ((real(z).^2 + imag(z).^2) <= escapeRadius2);
            count = count + inside;
            z = inside.*(z.*z + z0) + (1-inside).*z;
        end
        magZ2 = real(z).^2 + imag(z).^2;
        logCount = log( count + 1 - log( log( max(magZ2,escapeRadius2) ) / 2 ) / log(2) );
    end % computeMandelbrotCPU

    function logCount = computeMandelbrotGPU( xlim, numx, ylim, numy, maxIters )
        % Compute using GPUArray overloads.
        % In this version the main calculation is exactly as it was for the
        % CPU version, we have simply changed the input grid to be on the
        % GPU. When MATLAB encounters GPU data it tries to run any
        % functions on the GPU. This provides a simple way to see if the
        % GPU helps without altering your code.

        % Setup the input grid on the GPU
        escapeRadius2 = 400; % Square of escape radius
        x = parallel.gpu.GPUArray.linspace( xlim(1),  xlim(2), numx );
        y = parallel.gpu.GPUArray.linspace( ylim(1),  ylim(2), numy );
        count = parallel.gpu.GPUArray.zeros( numy, numx );
        [x0,y0] = meshgrid(x, y);
        z0 = complex( x0, y0 );

        % Calculate
        z = z0;
        for n = 0:maxIters
            inside = ((real(z).^2 + imag(z).^2) <= escapeRadius2);
            count = count + inside;
            z = inside.*(z.*z + z0) + (1-inside).*z;
        end
        magZ2 = real(z).^2 + imag(z).^2;
        logCount = log( count + 1 - log( log( max(magZ2,escapeRadius2) ) / 2 ) / log(2) );

        % Gather the result back to the CPU
        logCount = gather( logCount );
    end % computeMandelbrotGPU

    function logCount = computeMandelbrotArrayFun( xlim, numx, ylim, numy, maxIters )
        % Compute using GPU arrayfun.
        % The second way in which MATLAB can use the GPU is by placing your
        % algorithm inside a helper function and calling it using ARRAYFUN
        % with some GPU data as input. The helper function is converted
        % into native GPU code (PTX) and each element of the input array is
        % processed in a separate GPU thread. The helper function must
        % operate only one scalars. Here we have taken the "calculate" code
        % above and put it in a helper "mandelbrotViewerProcessElement"
        % that will be converted into native code and run on the GPU.

        % Create the input arrays
        escapeRadius2 = 400;
        x = parallel.gpu.GPUArray.linspace( xlim(1),  xlim(2), numx );
        y = parallel.gpu.GPUArray.linspace( ylim(1),  ylim(2), numy );
        [x0,y0] = meshgrid(x, y);

        % Calculate
        [logCount] = arrayfun( @mandelbrotViewerProcessElement, x0, y0, ...
            escapeRadius2, maxIters );

        % Gather the result back to the CPU
        logCount = gather( logCount );
    end % computeMandelbrotCPU

    function logCount = computeMandelbrotCUDAKernel( xlim, numx, ylim, numy, maxIters )
        % Use pre-existing CUDA/C++ code.
        % The final way in which MATLAB can use the GPU is by calling some
        % hand-written CUDA code. The "CUDAKernel" interface allows the
        % function to be specified along with the number of threads and
        % blocks to use. This requires some knowledge of how GPUs work, but
        % does allow you to easily use existing CUDA kernels with MATLAB
        % data.

        % Create the input arrays
        escapeRadius = 20;
        x = parallel.gpu.GPUArray.linspace( xlim(1),  xlim(2), numx );
        y = parallel.gpu.GPUArray.linspace( ylim(1),  ylim(2), numy );
        [x0,y0] = meshgrid(x, y);

        % Make sure we have sufficient blocks to cover the whole array
        numElements = numel( x0 );
        data.Kernel.ThreadBlockSize = [data.Kernel.MaxThreadsPerBlock,1,1];
        data.Kernel.GridSize = [ceil(numElements/data.Kernel.MaxThreadsPerBlock),1];

        % Call the kernel
        logCount = parallel.gpu.GPUArray.zeros( size( x0 ) );
        logCount = feval( data.Kernel, logCount, ...
            x0, y0, ...
            escapeRadius, maxIters, numElements );
        logCount = gather( logCount );
    end % computeMandelbrotCUDAKernel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Everything else is a callback or helper function.

    function out = createData()
        out = struct( ...
            'MaxIterations', 5000, ...
            'OrigXLim', [-2 1], ...
            'OrigY', 0, ...
            'XLim', [-2 1], ...
            'Y', 0, ...
            'CalculationMethods', {{
            'CPU'
            'GPU (simple)'
            'GPU ArrayFun'
            'CUDAKernel'
            }}, ...
            'IsAnimating', false, ...
            'NextLocation', 1, ...
            'LastFrameTime', now(), ...
            'WindowPixelSize', [100 100], ...
            'SelectedCalculationMethod', 'CUDAKernel', ...
            'ControlsVisible', true, ...
            'WriteVideo', false, ...
            'VideoWriter',  [] );

        % open writer
        if out.WriteVideo
            out.VideoWriter = VideoWriter('out.avi');
            out.VideoWriter.FrameRate = 20;
            out.VideoWriter.Quality = 90;
            open( out.VideoWriter );
        end

        % Read the location list from a local file
        out.LocationList = readLocationList();

        % Load the CUDA kernel
        out.Kernel = loadKernel();
    end % createData

    function out = createGUI()
        % Create the GUI, storing handles in the global GUI structure
        out.Window = figure( ...
            'Name', 'Mandelbrot viewer v1.1', ...
            'NumberTitle', 'off', ...
            'HandleVisibility', 'off', ...
            'MenuBar', 'none', ...
            'ToolBar', 'figure', ...
            'Renderer', 'ZBuffer' ); % Can't use painters as colormaps are broken for >256 colors!
        out.MainAxes = axes( ...
            'Parent', out.Window, ...
            'Position', [0 0 1 1], ...
            'XLim', data.XLim, ...
            'YLim', [-1 1], ...
            'CLim', log([1 data.MaxIterations]), ...
            'XTick', [], 'YTick', [], ...
            'DataAspectRatio', [1 1 1] );
        out.Image = image( ...
            'XData', [0 1], ...
            'YData', [0 1], ...
            'XLimInclude', 'off', ...
            'YLimInclude', 'off', ...
            'CData', nan, ...
            'CDataMapping', 'Scaled', ...
            'HandleVisibility', 'off', ...
            'Parent', out.MainAxes );
        % Add a line so that zooming works. Strange but true.
        line( 'Parent', out.MainAxes, 'XData', [-2 2], 'YData', [-2 2], ...
            'Visible', 'off', ...
            'HitTest', 'off' );
        set( out.MainAxes, 'XLimMode', 'manual', 'YLimMode', 'manual', ...
            'CLim', [0 1] );
        colormap( out.MainAxes, jet2(1000) );

        out.ControlPanel = uipanel( ...
            'Parent', out.Window, ...
            'BackgroundColor', 'k', ...
            'Units', 'Pixels', ...
            'Position', [10 10 155 78] );

        % Some text for showing compute time
        out.ComputeText = uicontrol( ...
            'Style', 'Text', ...
            'String', 'Computed in 0ms', ...
            'BackgroundColor', 'k', ...
            'ForegroundColor', 'g', ...
            'FontSize', 7, ...
            'Parent', out.ControlPanel, ...
            'Position', [5 16 145 14] );
        out.FrameRateText = uicontrol( ...
            'Style', 'Text', ...
            'String', 'Displaying at 0fps', ...
            'BackgroundColor', 'k', ...
            'ForegroundColor', 'g', ...
            'FontSize', 7, ...
            'Parent', out.ControlPanel, ...
            'Position', [5 2 145 14] );

        % Create a drop-down for selecting the calculation method
        out.MethodSelector = uicontrol( ...
            'Style', 'PopupMenu', ...
            'String', data.CalculationMethods, ...
            'Value', numel( data.CalculationMethods ), ...
            'FontSize', 7, ...
            'Parent', out.ControlPanel, ...
            'BackgroundColor', 0.8*[1 1 1], ...
            'Position', [5 54 145 16], ...
            'Callback', @onCalculationMethodChanged );
        out.AddButton = uicontrol( ...
            'Style', 'ToggleButton', ...
            'String', 'Add', ...
            'FontSize', 7, ...
            'Parent', out.ControlPanel, ...
            'BackgroundColor', 0.6*[1 1 1], ...
            'Position', [5 30 45 16], ...
            'TooltipString', 'Add the current location to the animation list', ...
            'Callback', @onAddPressed );
        out.ResetButton = uicontrol( ...
            'Style', 'ToggleButton', ...
            'String', 'Reset', ...
            'FontSize', 7, ...
            'Parent', out.ControlPanel, ...
            'BackgroundColor', 0.6*[1 1 1], ...
            'Position', [55 30 45 16], ...
            'TooltipString', 'Reset the view to the top', ...
            'Callback', @onResetPressed );
        out.AnimateButton = uicontrol( ...
            'Style', 'ToggleButton', ...
            'String', 'Animate', ...
            'FontSize', 7, ...
            'Parent', out.ControlPanel, ...
            'BackgroundColor', 0.6*[1 1 1], ...
            'Position', [105 30 45 16], ...
            'TooltipString', 'Start/stop animating between stored locations', ...
            'Callback', @onPlayPressed );

        % Remove some things we don't want from the toolbar and add a
        % toggle to the toolbar to hide the controls
        tb = findall( out.Window, 'Type', 'uitoolbar' );
        delete( findall( tb, 'Tag', 'Standard.FileOpen' ) );
        delete( findall( tb, 'Tag', 'Standard.NewFigure' ) );
        delete( findall( tb, 'Tag', 'Standard.EditPlot' ) );
        delete( findall( tb, 'Tag', 'Exploration.Brushing' ) );
        delete( findall( tb, 'Tag', 'Exploration.DataCursor' ) );
        delete( findall( tb, 'Tag', 'Exploration.Rotate' ) );
        delete( findall( tb, 'Tag', 'DataManager.Linking' ) );
        delete( findall( tb, 'Tag', 'Plottools.PlottoolsOn' ) );
        delete( findall( tb, 'Tag', 'Plottools.PlottoolsOff' ) );
        out.AnimateToggle = uitoggletool( ...
            'Parent', tb, ...
            'CData', readIcon( 'icon_play.png' ), ...
            'TooltipString', 'Start/stop animating between stored locations', ...
            'State', 'off', ...
            'Separator', 'on', ...
            'ClickedCallback', @onPlayToolbarPressed );
        out.ShowControlsToggle = uitoggletool( ...
            'Parent', tb, ...
            'CData', readIcon( 'icon_mandelControls.png' ), ...
            'TooltipString', 'Show/hide the Mandelbrot control panel', ...
            'State', 'on', ...
            'ClickedCallback', @onControlsTogglePressed );

        % Add listeners so that we can redraw when the axes are moved
        axHandle = handle( out.MainAxes );
        out.Listeners = [
            handle.listener( axHandle, findprop( axHandle, 'YLim' ), 'PropertyPostSet', @onLimitsChanged )
            ]; %#ok<NBRAK>
        % Also redraw if resized
        set( out.Window, 'ResizeFcn', @onFigureResize, ...
            'CloseRequestFcn', @onFigureClose );
    end % createGUI

    function onLimitsChanged( ~, ~ )
        redraw();
    end % onLimitsChanged

    function onFigureResize( ~, ~ )
        % Change the axes limits to exactly fit the figure
        pos = get( gui.Window, 'Position' );
        xlim = get( gui.MainAxes, 'XLim' );
        ylim = get( gui.MainAxes, 'YLim' );
        delta_ylim = ( diff( xlim )*pos(4)/pos(3) - diff( ylim ) ) / 2;
        data.WindowPixelSize = pos(3:4);
        % Set the YLim to give the correct aspect. This will trigger a
        % redraw
        set( gui.MainAxes, 'YLim', ylim + delta_ylim*[-1 1] );
    end % onFigureResize

    function onFigureClose( ~, ~ )
        % Clear up
        data.IsAnimating = false;
        if data.WriteVideo
            close( data.VideoWriter );
        end
        delete( gui.Window );
    end % onFigureClose

    function onCalculationMethodChanged( ~, ~ )
        idx = get( gui.MethodSelector, 'Value' );
        data.SelectedCalculationMethod = data.CalculationMethods{idx};
        redraw();
    end % onCalculationMethodChanged

    function onAddPressed( ~, ~ )
        disp('Add')
        fprintf( 'Adding location [%1.15f, %1.15f], %1.15f\n', ...
            data.XLim, data.Y );
        idx = numel( data.LocationList ) + 1;
        data.LocationList(idx).XLim = data.XLim;
        data.LocationList(idx).Y = data.Y;
        writeLocationList( data.LocationList );
        % Release the button
        set( gui.AddButton, 'Value', 0 );
    end % onAddPressed

    function onResetPressed( ~, ~ )
        disp('Reset')
        fprintf( 'Leaving location [%1.15f, %1.15f], %1.15f\n', ...
            data.XLim, data.Y );
        pos = get( gui.Window, 'Position' );
        aspect = pos(4)/pos(3);
        ylim = diff( data.OrigXLim ) * aspect / 2 * [-1 1];
        set( gui.MainAxes, 'XLim', data.OrigXLim, 'YLim', data.OrigY + ylim );
        set( gui.ResetButton, 'Value', 0 );
    end % onResetPressed

    function onPlayPressed( ~, ~ )
        disp('Play')
        if get( gui.AnimateButton, 'Value' )==1
            updateAnimationControls(true);
            while ishandle(gui.AnimateButton) && (get( gui.AnimateButton, 'Value' )==1)
                newXLim = data.LocationList(data.NextLocation).XLim;
                newY    = data.LocationList(data.NextLocation).Y;
                animatedMove( newXLim, newY );
                if numel( data.LocationList )>1
                    % Choose a random location
                    thisLocation = data.NextLocation;
                    while data.NextLocation == thisLocation
                        data.NextLocation = randi( numel( data.LocationList ), 1 );
                    end
                    fprintf( 'Next location: %d\n', data.NextLocation )
                else
                    % Only one location, so stop
                    data.IsAnimating = false;
                    updateAnimationControls( false );
                end
            end
        else
            data.IsAnimating = false;
            updateAnimationControls( false );
        end
    end % onPlayPressed

    function onPlayToolbarPressed( ~, evt )
        ison = strcmpi(get( gui.AnimateToggle, 'State' ), 'on');
        updateAnimationControls( ison );
        onPlayPressed(gui.AnimateButton, evt);
    end % onPlayToolbarPressed

    function updateAnimationControls( isAnimating )
        if isAnimating
            set( gui.AnimateButton, 'Value', 1 );
            set( gui.AnimateToggle, 'State', 'on' );
        else
            set( gui.AnimateButton, 'Value', 0 );
            set( gui.AnimateToggle, 'State', 'off' );
        end
        drawnow();
    end % updateAnimationControls

    function onControlsTogglePressed( ~, ~ )
        % Toggle the control panel on and off
        disp('Toggle controls')
        pos = get( gui.ControlPanel, 'Position' );
        if strcmpi( get( gui.ShowControlsToggle, 'State' ), 'off' )
            % Turn it off (move offscreen)
            pos(1) = -pos(3)-10;
        else
            % Turn it on (move onscreen)
            pos(1) = 10;
        end
        set( gui.ControlPanel, 'Position', pos );

    end

    function animatedMove( targetXLim, targetY )
        % Form a zoom path between the two
        data.IsAnimating = true;
        if isequal( data.XLim, targetXLim ) && isequal( data.Y, targetY )
            data.IsAnimating = false;
            return;
        end

        % Perform a zoom and translate arc
        maxNumSteps = 1000;
        distTravelled = sqrt( (mean( data.XLim ) - mean( targetXLim )).^2 ...
            + (data.Y - targetY).^2 );
        adjustRatio = exp( -10*linspace(-3,3,maxNumSteps).^2 );
        adjustRatio = adjustRatio - min(adjustRatio);
        ratio = cumsum( adjustRatio ); ratio = ratio / ratio(end);
        minXPath = interp1( [0,1], [data.XLim(1),targetXLim(1)], ratio );
        maxXPath = interp1( [0,1], [data.XLim(2),targetXLim(2)], ratio );
        maxXRange = max( maxXPath - maxXPath );
        xlimAdjust = max(0, 0.3*distTravelled - maxXRange);
        minXPath = minXPath - xlimAdjust*adjustRatio;
        maxXPath = maxXPath + xlimAdjust*adjustRatio;

        % Cull the ends if there's negligable motion. This helps to keep
        % things smooth but without long periods of no apparant motion.
        tolerance = 0.001;
        xRange = maxXPath - minXPath;
        firstGood = find( (xRange > (1+tolerance)*xRange(1)) | (xRange < (1-tolerance)*xRange(1)), 1, 'first' );
        if ~isempty( firstGood ) && firstGood > 2
            toCull = 2:firstGood-1;
        else
            toCull = [];
        end
        lastGood = find( (xRange > (1+tolerance)*xRange(end)) | (xRange < (1-tolerance)*xRange(end)), 1, 'last' );
        if ~isempty( lastGood ) && lastGood < numel(xRange)-1
            toCull = [toCull, lastGood:numel(xRange)-1];
        end
        ratio(toCull) = [];
        minXPath(toCull) = [];
        maxXPath(toCull) = [];
        xRange(toCull) = [];

        if ~isempty( ratio )
            % Work out the aspect ratio
            pos = get( gui.Window, 'Position' );
            aspect = pos(4)/pos(3);

            YPath = interp1( [0,1], [data.Y,targetY], ratio );
            heightPath = aspect*xRange;
            minYPath = YPath - 0.5*heightPath;
            maxYPath = YPath + 0.5*heightPath;
            for ii=1:numel(ratio)
                % Setting the limits will cause a redraw
                set( gui.MainAxes, ...
                    'XLim', [minXPath(ii),maxXPath(ii)], ...
                    'YLim', [minYPath(ii),maxYPath(ii)] );
                if data.IsAnimating == false
                    break;
                end
            end
        end
        data.IsAnimating = false;
        % Do a final redraw at full res
        redraw();
    end % animatedMove

    function kernel = loadKernel()
        thisDir = fileparts( mfilename( 'fullpath' ) );
        baseName = 'mandelbrotViewerProcessElement';
        data.CUDAFile = fullfile( thisDir, [baseName,'.cu'] );
        ptxname = [baseName,'.',parallel.gpu.ptxext];
        data.PTXFile = fullfile( thisDir, ptxname );
        if exist( data.PTXFile, 'file' ) ~= 2
            close( gui.Window );
            error( 'mandelbrotViewer:MissingPTX', 'Could not find ''%s''. Please use NVCC to compile it', ptxname );
        end
        kernel = parallel.gpu.CUDAKernel( data.PTXFile, data.CUDAFile );
    end % loadKernel

    function redraw()
        % Protect against the window closing
        if ~ishandle(gui.MainAxes)
            return;
        end
        % To work out what to draw and at what resolution we need the axis
        % limits and pixel counts.
        xlim = get(gui.MainAxes,'XLim');
        ylim = get(gui.MainAxes,'YLim');
        data.XLim = xlim;
        data.Y = mean( ylim );
        imWidth = data.WindowPixelSize(1);
        imHeight = data.WindowPixelSize(2);
        if data.IsAnimating && (imWidth*imHeight>600000)
            % To speed up animations with large windows, subsample by 2
            imWidth = round(imWidth/2);
            imHeight = round(imHeight/2);
        end

        zoomLevel = imWidth / diff( xlim );
        maxIterations = min( data.MaxIterations, 200 + 0.1*sqrt(zoomLevel) );

        % Call the computation
        t = tic;
        switch( data.SelectedCalculationMethod )
            case 'CUDAKernel'
                logCount = computeMandelbrotCUDAKernel( xlim, imWidth, ...
                    ylim, imHeight, ...
                    maxIterations );

            case 'CPU'
                logCount = computeMandelbrotCPU( xlim, imWidth, ...
                    ylim, imHeight, ...
                    maxIterations );

            case 'GPU (simple)'
                logCount = computeMandelbrotGPU( xlim, imWidth, ...
                    ylim, imHeight, ...
                    maxIterations );

            case 'GPU ArrayFun'
                logCount = computeMandelbrotArrayFun( xlim, imWidth, ...
                    ylim, imHeight, ...
                    maxIterations );

            otherwise
                error( 'mandelbrotViewer:BadMethod', 'Unrecognised calculation method ''%s''', data.SelectedCalculationMethod );
        end

        computeTime = toc(t);

        minCount = min( logCount(:) );
        logCount = (logCount - minCount) ./ (log(maxIterations+1)-minCount);
        % Guard against a closed window
        if ~ishandle( gui.Image )
            return;
        end
        set( gui.Image, ...
            'XData', xlim, ...
            'YData', ylim, ...
            'CData', logCount );
        if data.ControlsVisible
            set( gui.ComputeText, 'String', sprintf( 'Computed in %dms', round(1000*computeTime) ) )

            % Capture the current time for frame-rate calculations
            thisFrameTime = now();
            framerate = 1 / (86400*(thisFrameTime - data.LastFrameTime)); % convert days to seconds
            set( gui.FrameRateText, 'String', sprintf( 'Displaying at %dfps', round(framerate) ) )
            data.LastFrameTime = thisFrameTime;

            % Force a redraw
            drawnow();
        end

        % Capture!
        if data.WriteVideo
            t0 = now();
            currFrame = getframe( gui.Window );
            writeVideo( data.VideoWriter, currFrame );
            % Also reset the frame time to exclude the video writing
            delta_t = now() - t0;
            data.LastFrameTime = data.LastFrameTime + delta_t;
        end

    end % redraw

    function locations = readLocationList()
        fid = fopen( 'locations.csv', 'rt' );
        if fid<0
            close( gui.Window );
            error( 'mandelbrotViewer:BadLocationRead', 'Could not open location list for reading: ''locations.csv''' );
        end

        locData = textscan( fid, '%f,%f,%f' );
        N = size( locData{1}, 1 );
        if N<1
            close( gui.Window );
            error( 'mandelbrotViewer:EmptyLocationFile', 'No locations found in: ''locations.csv''' );
        end
        locations = struct( ...
            'XLim', cell( N, 1 ), ...
            'Y', cell( N, 1 ) );
        for ii=1:N
            locations(ii).XLim = [locData{1}(ii), locData{2}(ii)];
            locations(ii).Y = locData{3}(ii);
        end

        fclose( fid );
    end % readLocationList

    function writeLocationList( locations )
        fid = fopen( 'locations.csv', 'wt' );
        if fid<0
            error( 'mandelbrotViewer:BadLocationWrite', 'Could not open location list for writing: ''locations.csv''' );
        end
        N = numel( locations );
        for ii=1:N
            fprintf( fid, '%1.15f,%1.15f,%1.15f\n', ...
                locations(ii).XLim(1), ...
                locations(ii).XLim(2), ...
                locations(ii).Y );
        end

        fclose( fid );
    end % writeLocationList

    function cdata = readIcon( filename )
        [cdata,~,alpha] = imread( filename );
        idx = find( ~alpha );
        page = size(cdata,1)*size(cdata,2);
        cdata = double( cdata ) / 255;
        cdata(idx) = nan;
        cdata(idx+page) = nan;
        cdata(idx+2*page) = nan;
    end % readIcon

    function matlabVersionCheck()
        % R2011a is v7.12
        majorMinor = sscanf( version, '%d.%d' );
        if (majorMinor(1)<7) || (majorMinor(1)==7 && majorMinor(2)<13)
            error( 'mandelbrotViewer:MATLABTooOld', 'mandelbrotViewer requires MATLAB R2011b or above.' );
        end
    end % matlabVersionCheck

    function gpuCheck()
        try
            d = gpuDevice();
        catch err
            error( 'mandelbrotViewer:NoGPU', 'mandelbrotViewer requires a GPU and none appear to be availble. Type "gpuDevice" for more information.' );
        end
        if ~d.DeviceSupported
            error( 'mandelbrotViewer:GPUNotSupported', 'The selected GPU is not supported. Type "gpuDevice" for more information.' );
        end
    end % matlabVersionCheck

end % mandelbrotViewer

repo: Aero-matlab


text:
% The COBRAToolbox: testElementalBalance.m
%
% Purpose:
%     - Tests computeMW functionality
%
% Authors:
%     - CI integration: Laurent Heirendt

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testElementalBalance'));
cd(fileDir);

% load the model and data
load('testElementalBalanceData.mat');

% run elmental balance with no optional functions
[MW, Ematrix] = computeMW(model);

% check that the molecular weights have been calculated properly
assert(isequal(MW, stdMW))

% check that the Ematrix has been computed correctly
assert(isequal(Ematrix, stdEmatrix))

% run computeMW with a specific met list
[MW, Ematrix, elements, knownWeights, unknownElements] = computeMW(model, model.mets(25:35), false);

% check that the molecular weights have been calculated properly
assert(isequal(MW, stdMW2))

% check that the Ematrix has been computed correctly
assert(isequal(Ematrix, stdEmatrix2))

% check the list of returned elements
assert(isequal(elements, {'C', 'N', 'O', 'H', 'P', 'Other'}))

% empty knownWeights and unknownElements (non-empty only if the genericFomrula
% flag is turned on)
assert(isempty(knownWeights) & isempty(unknownElements))

% check warning message
model = addMetabolite(model, 'uranium[c]', 'metFormula', 'U');
diary('testElementBalance_warning.txt')
computeMW(model, 'uranium[c]');
diary off
f = fopen('testElementBalance_warning.txt', 'r');
l = fgets(f);
text = '';
while ~isequal(l, -1)
    text = [text, l];
    l = fgets(f);
end
fclose(f);
delete('testElementBalance_warning.txt')
% the warning should contain the following:
re = regexp(text, 'formula \=');
assert(~isempty(re))
text = text(re:end);
re = regexp(text, '''U''');
assert(~isempty(re))
text = text(re:end);
re = regexp(text, 'comp \=');
assert(~isempty(re))
text = text(re:end);
re = regexp(text, 'U');
assert(~isempty(re))

% test genericFormula = true
model = addMetabolite(model, 'hypothetical1[c]', 'metFormula', '((CH2O)3H2O)2H2O');
nMets = numel(model.mets);
model = addMetabolite(model, 'hypothetical2[c]', 'metFormula', '(H2O)2CuSO4Element2');
model = addMetabolite(model, 'hypothetical3[c]', 'metFormula', 'C0H0');
model = addMetabolite(model, 'photon[c]', 'metFormula', 'Mass0');
[MW, Ematrix, elements, knownWeights, unknownElements] = computeMW(model, [], false, true);
% check weights
assert(all(isnan(MW) == isnan(stdMW3)) && max(abs(MW(~isnan(MW)) - stdMW3(~isnan(MW)))) < 1e-5)
% check elements and Ematrix
[yn, id] = ismember(elements, {'C', 'H', 'O', 'P', 'N', 'S', 'U', 'Cu', 'Element', 'Mass'});
assert(all(yn) & numel(elements) == 10)
assert(isequal(Ematrix, stdEmatrix3(:, id)))
% check weights for the known part
assert(all(knownWeights(1:nMets) == MW(1:nMets)))
assert(abs(knownWeights(nMets + 1) - 195.6392) < 1e-4)
assert(all(knownWeights(nMets + (2:3)) == 0))
% check unknownElements
assert(isequal(sort({'Element', 'Mass'}), sort(unknownElements)))

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function a = subsref(t,s)
%SUBSREF Subscripted reference for a ttensor.
%
%   Examples
%   core = tensor(rand(2,2,2));
%   X = TTENSOR(core, rand{4,2), rand(5,2),rand(3,2));
%   X.core %<-- returns core array
%   X.U %<-- returns a cell array of three matrices
%   X.U{1} %<-- returns the matrix corresponding to the first mode.
%
%   See also TTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

switch s(1).type
    case '.'
        switch s(1).subs
            case {'core','lambda'}
                a = tt_subsubsref(t.core,s);
            case {'U','u'}
                a = tt_subsubsref(t.u,s);
            otherwise
                error(['No such field: ', s.subs]);
        end
    case '()'
	error('Subsref with () not supported for ttensor.');
    case '{}'
	new_s(1).type = '.';
	new_s(1).subs = 'u';
	new_s(2:length(s)+1) = s;
	a = subsref(t, new_s);
     otherwise
        error('Invalid subsref.');
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ d ] = gsp_vec2mat( d,Nf )
%GSP_VEC2MAT vector to matrix transform
%   Usage: d  = gsp_vec2mat( d,Nf );
%
%   Input parameters:
%       d       : Data
%       Nf      : Number of filter
%
%   Ouput parameter
%       d       : Data
%
%   Reshape the data from the vector form to the matrix form

% TESTING: test_filter

[M,N] = size(d);

d = reshape(d,M/Nf,Nf,N);

end

repo: gspbox


text:
classdef SymbolicVoigtRotationMatrixGenerator < handle

    properties (Access = public)
        VoigtRotationMatrix
    end

    properties (Access = private)
        RotationMatrix
        Stress
        RotatedStress
    end

    methods (Access = public)

        function obj = SymbolicVoigtRotationMatrixGenerator()
            obj.createRotationMatrix();
            obj.createStressTensor();
            obj.computeRotatedStressTensor();
            obj.obtainVoigtRotationMatrix()
        end

    end

    methods (Access = private)

        function obj = createRotationMatrix(obj)
            theta    = obj.createRotationAngle();
            vect     = obj.createNormalVector();
            rotator = VectorRotator(theta,vect);
            obj.RotationMatrix = rotator.getRotationMatrix();
        end

        function u = createNormalVector(obj)
            u = sym('u',[3 1],'real');
        end

        function theta = createRotationAngle(obj)
            theta = sym('theta','real');
        end

        function createStressTensor(obj)
            obj.Stress = StressTensor();
            obj.Stress.tensor = sym('s',[3 3],'real');
            Tens = obj.Stress.tensor;
            Tens = obj.Stress.symmetrizeWithUpperDiagonal(Tens);
            obj.Stress.tensor = Tens;
            obj.Stress.transformTensor2Voigt();
        end

        function computeRotatedStressTensor(obj)
            R = obj.RotationMatrix;
            S = obj.Stress.tensor;
            RotS = R'*S*R;
            obj.RotatedStress = StressTensor();
            obj.RotatedStress.tensor = simplify(RotS);
            obj.RotatedStress.transformTensor2Voigt();
        end

        function obtainVoigtRotationMatrix(obj)
            RotStre = obj.RotatedStress.tensorVoigt;
            Stre    = obj.Stress.tensorVoigt;
            DimVoigt = length(RotStre);
            RotMatrix = sym(zeros(DimVoigt,DimVoigt));
            for iStress = 1:DimVoigt
                RS = RotStre(iStress);
                for jStress = 1:DimVoigt
                    S = Stre(jStress);
                    [TensorValue,~] = coeffs(RS,S);
                    TensorValue = obj.takeApropiateComponent(TensorValue);
                    RotMatrix(iStress,jStress) = TensorValue;
                end
            end
            obj.VoigtRotationMatrix = simplify(RotMatrix);
        end

        function value = takeApropiateComponent(obj,value)
            if isempty(value)
                value = [];
            else
                Dim = length(value);
                if Dim == 1
                    value = 0;
                else
                    value = value(1);
                end
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function DrawFaceOnImg( img, shape, out_loc, orig_bbox, visibilities)
%DRAWFACEONIMG Summary of this function goes here
%   Detailed explanation goes here
    f = figure('visible','off','Position', [100, 100, 600, 600]);

    if(nargin > 3)
        [height_img, width_img,~] = size(img);
        width = orig_bbox(3) - orig_bbox(1);
        height = orig_bbox(4) - orig_bbox(2);

        img_min_x = max(int32(orig_bbox(1) - width/3),1);
        img_max_x = min(int32(orig_bbox(3) + width/3),width_img);

        img_min_y = max(int32(orig_bbox(2) - height/3),1);
        img_max_y = min(int32(orig_bbox(4) + height/3),height_img);

        shape(:,1) = shape(:,1) - double(img_min_x);
        shape(:,2) = shape(:,2) - double(img_min_y);
        img = img(img_min_y:img_max_y, img_min_x:img_max_x, :);
    end

    % Downsize the image if it is very big
    if(size(img,1) > 500 || size(img,2) > 500)
        scale = min(500 / size(img,1), 500 / size(img,2));
        img = imresize(img, scale);
        shape = scale * shape;
    end
    if(nargin > 4)
        % valid points to draw (not to draw
        % occluded ones)
        v_points = visibilities;
    else
        v_points = true(size(shape,1),1);
    end

    if(max(img(:)) > 1)
        imshow(double(img)/255, 'Border', 'tight');
    else
        imshow(double(img), 'Border', 'tight');
    end
    hold on;

    plot(shape(v_points,1), shape(v_points,2),'.r','MarkerSize',20);
    plot(shape(v_points,1), shape(v_points,2),'.b','MarkerSize',10);
    print(f, '-djpeg', out_loc);
    close(f);
end

author: TadasBaltrusaitis

repo: OpenFace