text:
function [m,v,w,g,f,pp,gg]=gaussmix(x,c,l,m0,v0,w0,wx)
%GAUSSMIX fits a gaussian mixture pdf to a set of data observations [m,v,w,g,f]=(x,c,l,m0,v0,w0,wx)
%
% Usage:
%    (1) [m,v,w]=gaussmix(x,[],[],k);    % create GMM with k mixtures and diagonal covariances
%    (2) [m,v,w]=gaussmix(x,[],[],k,'v);    % create GMM with k mixtures and full covariances
%
% Inputs: n data values, k mixtures, p parameters, l loops
%
%     X(n,p)   Input data vectors, one per row.
%     c(1)     Minimum variance of normalized data (Use [] to take default value of 1/n^2)
%     L        The integer portion of l gives a maximum loop count. The fractional portion gives
%              an optional stopping threshold. Iteration will cease if the increase in
%              log likelihood density per data point is less than this value. Thus l=10.001 will
%              stop after 10 iterations or when the increase in log likelihood falls below
%              0.001.
%              As a special case, if L=0, then the first three outputs are omitted.
%              Use [] to take default value of 100.0001
%     M0       Number of mixtures required (or initial mixture means - see below)
%     V0       Initialization mode:
%                'f'    Initialize with K randomly selected data points [default]
%                'p'    Initialize with centroids and variances of random partitions
%                'k'    k-means algorithm ('kf' and 'kp' determine initialization)
%                'h'    k-harmonic means algorithm ('hf' and 'hp' determine initialization) [default]
%                's'    do not scale data during initialization to have equal variances
%                'm'    M0 contains the initial centres
%                'v'    full covariance matrices
%              Mode 'hf' [the default] generally gives the best results but 'f' is faster and often OK
%     W0(n,1)  Data point weights
%
%     Alternatively, initial values for M0, V0 and W0 can be given  explicitly:
%
%     M0(k,p)  Initial mixture means, one row per mixture.
%     V0(k,p)  Initial mixture variances, one row per mixture.
%      or V0(p,p,k)  one full-covariance matrix per mixture
%     W0(k,1)  Initial mixture weights, one per mixture. The weights should sum to unity.
%     WX(n,1)  Data point weights
%
% Outputs: (Note that M, V and W are omitted if L==0)
%
%     M(k,p)   Mixture means, one row per mixture. (omitted if L==0)
%     V(k,p)   Mixture variances, one row per mixture. (omitted if L==0)
%       or V(p,p,k) if full covariance matrices in use (i.e. either 'v' option or V0(p,p,k) specified)
%     W(k,1)   Mixture weights, one per mixture. The weights will sum to unity. (omitted if L==0)
%     G       Average log probability of the input data points.
%     F        Fisher's Discriminant measures how well the data divides into classes.
%              It is the ratio of the between-mixture variance to the average mixture variance: a
%              high value means the classes (mixtures) are well separated.
%     PP(n,1)  Log probability of each data point
%     GG(l+1,1) Average log probabilities at the beginning of each iteration and at the end
%
% The fitting procedure uses one of several initialization methods to create an initial guess
% for the mixture centres and then uses the EM (expectation-maximization) algorithm to refine
% the guess. Although the EM algorithm is deterministic, the initialization procedures use
% random numbers and so the routine will not give identical answers if you call it multiple
% times with the same input data.

%  Bugs/Suggestions
%     (1) Allow processing in chunks by outputting/reinputting an array of sufficient statistics
%     (2) Other initialization options:
%              'l'    LBG algorithm
%              'm'    Move-means (dog-rabbit) algorithm
%     (3) Allow updating of weights-only, not means/variances

%      Copyright (C) Mike Brookes 2000-2009
%      Version: $Id: gaussmix.m 7784 2016-04-15 11:09:50Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n,p]=size(x);
wn=ones(n,1);
mx0=sum(x,1)/n;         % calculate mean and variance of input data in each dimension
vx0=sum(x.^2,1)/n-mx0.^2;
sx0=sqrt(vx0);
sx0(sx0==0)=1;      % do not divide by zero when scaling
scaled=0;           % data is not yet scaled
memsize=voicebox('memsize');    % set memory size to use
if isempty(c)
    c=1/n^2;
else
    c=c(1);         % just to prevent legacy code failing
end
fulliv=0;           % initial variance is not full
if isempty(l)
    l=100+1e-4;         % max loop count + stopping threshold
end
if nargin<5 || isempty(v0) || ischar(v0)             % no initial values specified for m0, v0, w0
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %  No initialvalues given, so we must use k-means or equivalent
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if nargin<6
        if nargin<5 || isempty(v0)
            v0='hf';                 % default initialization mode: hf
        end
        wx=wn;                      % no data point weights
    else
        wx=w0(:);                   % data point weights
    end
    if any(v0=='m')
        k=size(m0,1);
    else
        k=m0;
    end
    fv=any(v0=='v');                % full covariance matrices requested
    if n<=k                         % each data point can have its own mixture
        xs=(x-mx0(wn,:))./sx0(wn,:);          % scale the data
        m=xs(mod((1:k)-1,n)+1,:);   % just include all points several times
        v=zeros(k,p);               % will be set to floor later
        w=zeros(k,1);
        w(1:n)=1/n;
        if l>0
            l=0.1;                  % no point in iterating
        end
    else                            % more points than mixtures
        if any(v0=='s')
            xs=x;                   % do not scale data during initialization
        else
            xs=(x-mx0(wn,:))./sx0(wn,:);  % else scale now
            if any(v0=='m')
                m=(m0-mx0(ones(k,1),:))./sx0(ones(k,1),:);  % scale specified means as well
            end
        end
        w=repmat(1/k,k,1);                      % all mixtures equally likely
        if any(v0=='k')                         % k-means initialization
            if any(v0=='m')
                [m,e,j]=v_kmeans(xs,k,m);
            elseif any(v0=='p')
                [m,e,j]=v_kmeans(xs,k,'p');
            else
                [m,e,j]=v_kmeans(xs,k,'f');
            end
        elseif any(v0=='h')                     % k-harmonic means initialization
            if any(v0=='m')
                [m,e,j]=kmeanhar(xs,k,[],4,m);
            else
                if any(v0=='p')
                    [m,e,j]=kmeanhar(xs,k,[],4,'p');
                else
                    [m,e,j]=kmeanhar(xs,k,[],4,'f');
                end
            end
        elseif any(v0=='p')                     % Initialize using a random partition
            j=ceil(rand(n,1)*k);                % allocate to random clusters
            j(rnsubset(k,n))=1:k;               % but force at least one point per cluster
            for i=1:k
                m(i,:)=mean(xs(j==i,:),1);
            end
        else
            if any(v0=='m')
                m=m0;                           % use specified centres
            else
                m=xs(rnsubset(k,n),:);          % Forgy initialization: sample k centres without replacement [default]
            end
            [e,j]=v_kmeans(xs,k,m,0);             % find out the cluster allocation
        end
        if any(v0=='s')
            xs=(x-mx0(wn,:))./sx0(wn,:);      % scale data now if not done previously
        end
        v=zeros(k,p);                   % diagonal covariances
        w=zeros(k,1);
        for i=1:k
            ni=sum(j==i);               % number assigned to this centre
            w(i)=(ni+1)/(n+k);          % weight of this mixture
            if ni
                v(i,:)=sum((xs(j==i,:)-repmat(m(i,:),ni,1)).^2,1)/ni;
            else
                v(i,:)=zeros(1,p);
            end
        end
    end
else
    %%%%%%%%%%%%%%%%%%%%%%%%
    % use initial values given as input parameters
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [k,p]=size(m0);
    xs=(x-mx0(wn,:))./sx0(wn,:);          % scale the data
    m=(m0-mx0(ones(k,1),:))./sx0(ones(k,1),:);          % and the means
    v=v0;
    w=w0;
    fv=ndims(v)>2 || size(v,1)>k;                       % full covariance matrix is supplied
    if fv
        mk=eye(p)==0;                                    % off-diagonal elements
        fulliv=any(v(repmat(mk,[1 1 k]))~=0);            % check if any are non-zero
        if ~fulliv
            v=reshape(v(repmat(~mk,[1 1 k])),p,k)'./repmat(sx0.^2,k,1);   % just pick out and scale the diagonal elements for now
        else
            v=v./repmat(sx0'*sx0,[1 1 k]);              % scale the full covariance matrix
        end
    end
    if nargin<7
        wx=wn;              % no data point weights
    end
end
if length(wx)~=n
    error('%d datapoints but %d weights',n,length(wx));
end
lsx=sum(log(sx0));
xsw=xs.*repmat(wx,1,p); % weighted data points
nwt=sum(wx);        % number of data points counting duplicates
if ~fulliv          % initializing with diagonal covariance
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Diagonal Covariance matrices  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    v=max(v,c);         % apply the lower bound
    xs2=xs.^2.*repmat(wx,1,p);          % square and weight the data for variance calculations

    % If data size is large then do calculations in chunks

    nb=min(n,max(1,floor(memsize/(8*p*k))));    % chunk size for testing data points
    nl=ceil(n/nb);                  % number of chunks
    jx0=n-(nl-1)*nb;                % size of first chunk

    im=repmat(1:k,1,nb); im=im(:);
    th=(l-floor(l))*n;
    sd=(nargout > 3*(l~=0)); % = 1 if we are outputting log likelihood values
    lp=floor(l)+sd;   % extra loop needed to calculate final G value

    lpx=zeros(1,n);             % log probability of each data point
    wk=ones(k,1);
    wp=ones(1,p);
    wnb=ones(1,nb);
    wnj=ones(1,jx0);

    % EM loop

    g=0;                           % dummy initial value for comparison
    gg=zeros(lp+1,1);
    ss=sd;                       % initialize stopping count (0 or 1)
    for j=1:lp
        g1=g;                    % save previous log likelihood (2*pi factor omitted)
        m1=m;                       % save previous means, variances and weights
        v1=v;
        w1=w;
        vi=-0.5*v.^(-1);                % data-independent scale factor in exponent
        lvm=log(w)-0.5*sum(log(v),2);   % log of external scale factor (excluding -0.5*p*log(2pi) term)

        % first do partial chunk (of length jx0)

        jx=jx0;
        ii=1:jx;                        % indices of data points in this chunk
        kk=repmat(ii,k,1);              % kk(jx,k): one row per data point, one column per mixture
        km=repmat(1:k,1,jx);            % km(jx,k): one row per data point, one column per mixture
        py=reshape(sum((xs(kk(:),:)-m(km(:),:)).^2.*vi(km(:),:),2),k,jx)+lvm(:,wnj); % py(k,jx) pdf of each point with each mixture
        mx=max(py,[],1);                % mx(1,jx) find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
        ps=sum(px,1);                   % total normalized likelihood of each data point
        px=px./ps(wk,:);                % relative mixture probabilities for each data point (columns sum to 1)
        lpx(ii)=log(ps)+mx;
        pk=px*wx(ii);                       % pk(k,1) effective number of data points for each mixture (could be zero due to underflow)
        sx=px*xsw(ii,:);
        sx2=px*xs2(ii,:);
        for il=2:nl                     % process the data points in chunks
            ix=jx+1;
            jx=jx+nb;                   % increment upper limit
            ii=ix:jx;                   % indices of data points in this chunk
            kk=repmat(ii,k,1);
            py=reshape(sum((xs(kk(:),:)-m(im,:)).^2.*vi(im,:),2),k,nb)+lvm(:,wnb);
            mx=max(py,[],1);            % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:));        % find normalized probability of each mixture for each datapoint
            ps=sum(px,1);               % total normalized likelihood of each data point
            px=px./ps(wk,:);            % relative mixture probabilities for each data point (columns sum to 1)
            lpx(ii)=log(ps)+mx;
            pk=pk+px*wx(ii);                % pk(k,1) effective number of data points for each mixture (could be zero due to underflow)
            sx=sx+px*xsw(ii,:);
            sx2=sx2+px*xs2(ii,:);
        end
        g=lpx*wx;                       % total log probability summed over all data points
        gg(j)=g;                        % save log prob at each iteration
        w=pk/nwt;                       % normalize to get the weights
        if pk                           % if all elements of pk are non-zero
            m=sx./pk(:,wp);             % calculate mixture means
            v=sx2./pk(:,wp);            % and variances
        else
            wm=pk==0;                   % mask indicating mixtures with zero weights
            nz=sum(wm);              	% number of zero-weight mixtures
            [vv,mk]=sort(lpx);          % find the lowest probability data points
            m=zeros(k,p);               % initialize means and variances to zero (variances are floored later)
            v=m;
            m(wm,:)=xs(mk(1:nz),:); 	% set zero-weight mixture means to worst-fitted data points
            w(wm)=1/n;               	% set these weights non-zero
            w=w*n/(n+nz);            	% normalize so the weights sum to unity
            wm=~wm;                 	% mask for non-zero weights
            m(wm,:)=sx(wm,:)./pk(wm,wp);  % recalculate means and variances for mixtures with a non-zero weight
            v(wm,:)=sx2(wm,:)./pk(wm,wp);
        end
        v=max(v-m.^2,c);                % apply floor to variances
        if g-g1<=th && j>1
            if ~ss, break; end  %  stop
            ss=ss-1;       % stop next time
        end

    end
    if sd && ~fv  % we need to calculate the final probabilities
        pp=lpx'-0.5*p*log(2*pi)-lsx;   % log of total probability of each data point
        gg=gg(1:j)/n-0.5*p*log(2*pi)-lsx;    % average log prob at each iteration
        g=gg(end);
        %     gg' % *** DEBUG ***
        m=m1;       % back up to previous iteration
        v=v1;
        w=w1;
        mm=sum(m,1)/k;
        f=(m(:)'*m(:)-k*mm(:)'*mm(:))/sum(v(:));
    end
    if ~fv
        m=m.*sx0(ones(k,1),:)+mx0(ones(k,1),:);	% unscale means
        v=v.*repmat(sx0.^2,k,1);                % and variances
    else
        v1=v;
        v=zeros(p,p,k);
        mk=eye(p)==1;                           % mask for diagonal elements
        v(repmat(mk,[1 1 k]))=v1';              % set from v1
    end
end
if fv              % check if full covariance matrices were requested
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Full Covariance matrices  %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    pl=p*(p+1)/2;
    lix=1:p^2;
    cix=repmat(1:p,p,1);
    rix=cix';
    lix(cix>rix)=[];                                        % index of lower triangular elements
    cix=cix(lix);                                           % index of lower triangular columns
    rix=rix(lix);                                           % index of lower triangular rows
    dix=find(rix==cix);
    lixi=zeros(p,p);
    lixi(lix)=1:pl;
    lixi=lixi';
    lixi(lix)=1:pl;                                        % reverse index to build full matrices
    v=reshape(v,p^2,k);
    v=v(lix,:)';                                            % lower triangular in rows

    % If data size is large then do calculations in chunks

    nb=min(n,max(1,floor(memsize/(24*p*k))));    % chunk size for testing data points
    nl=ceil(n/nb);                  % number of chunks
    jx0=n-(nl-1)*nb;                % size of first chunk
    %
    th=(l-floor(l))*n;
    sd=(nargout > 3*(l~=0)); % = 1 if we are outputting log likelihood values
    lp=floor(l)+sd;   % extra loop needed to calculate final G value
    %
    lpx=zeros(1,n);             % log probability of each data point
    wk=ones(k,1);
    wp=ones(1,p);
    wpl=ones(1,pl);             % 1 index for lower triangular matrix
    wnb=ones(1,nb);
    wnj=ones(1,jx0);

    % EM loop

    g=0;                        % dummy initial value for comparison
    gg=zeros(lp+1,1);
    ss=sd;                      % initialize stopping count (0 or 1)
    vi=zeros(p*k,p);            % stack of k inverse cov matrices each size p*p
    vim=zeros(p*k,1);       	% stack of k vectors of the form inv(v)*m
    mtk=vim;                  	% stack of k vectors of the form m
    lvm=zeros(k,1);
    wpk=repmat((1:p)',k,1);
    for j=1:lp
        g1=g;               	% save previous log likelihood (2*pi factor omitted)
        m1=m;                	% save previous means, variances and weights
        v1=v;
        w1=w;
        for ik=1:k

            % these lines added for debugging only
            %             vk=reshape(v(k,lixi),p,p);
            %             condk(ik)=cond(vk);
            %%%%%%%%%%%%%%%%%%%%
            [uvk,dvk]=eig(reshape(v(ik,lixi),p,p));	% convert lower triangular to full and find eigenvalues
            dvk=max(diag(dvk),c);                	% apply variance floor to eigenvalues
            vik=-0.5*uvk*diag(dvk.^(-1))*uvk';      % calculate inverse
            vi((ik-1)*p+(1:p),:)=vik;               % vi contains all mixture inverses stacked on top of each other
            vim((ik-1)*p+(1:p))=vik*m(ik,:)';       % vim contains vi*m for all mixtures stacked on top of each other
            mtk((ik-1)*p+(1:p))=m(ik,:)';           % mtk contains all mixture means stacked on top of each other
            lvm(ik)=log(w(ik))-0.5*sum(log(dvk));       % vm contains the weighted sqrt of det(vi) for each mixture
        end
        %
        %         % first do partial chunk
        %
        jx=jx0;
        ii=1:jx;
        xii=xs(ii,:).';
        py=reshape(sum(reshape((vi*xii-vim(:,wnj)).*(xii(wpk,:)-mtk(:,wnj)),p,jx*k),1),k,jx)+lvm(:,wnj);
        mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
        px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
        ps=sum(px,1);                   % total normalized likelihood of each data point
        px=px./ps(wk,:);                % relative mixture probabilities for each data point (columns sum to 1)
        lpx(ii)=log(ps)+mx;
        pk=px*wx(ii);                       % effective number of data points for each mixture (could be zero due to underflow)
        sx=px*xsw(ii,:);
        sx2=px*(xsw(ii,rix).*xs(ii,cix));	% accumulator for variance calculation (lower tri cov matrix as a row)
        for il=2:nl
            ix=jx+1;
            jx=jx+nb;        % increment upper limit
            ii=ix:jx;
            xii=xs(ii,:).';
            py=reshape(sum(reshape((vi*xii-vim(:,wnb)).*(xii(wpk,:)-mtk(:,wnb)),p,nb*k),1),k,nb)+lvm(:,wnb);
            mx=max(py,[],1);                % find normalizing factor for each data point to prevent underflow when using exp()
            px=exp(py-mx(wk,:));            % find normalized probability of each mixture for each datapoint
            ps=sum(px,1);                   % total normalized likelihood of each data point
            px=px./ps(wk,:);                % relative mixture probabilities for each data point (columns sum to 1)
            lpx(ii)=log(ps)+mx;
            pk=pk+px*wx(ii);                    % effective number of data points for each mixture (could be zero due to underflow)
            sx=sx+px*xsw(ii,:);             % accumulator for mean calculation
            sx2=sx2+px*(xsw(ii,rix).*xs(ii,cix));	% accumulator for variance calculation
        end
        g=lpx*wx;                   % total log probability summed over all data points
        gg(j)=g;                    % save convergence history
        w=pk/nwt;               	% w(k,1) normalize to get the column of weights
        if pk                       % if all elements of pk are non-zero
            m=sx./pk(:,wp);         % find mean and mean square
            v=sx2./pk(:,wpl);
        else
            wm=pk==0;                       % mask indicating mixtures with zero weights
            nz=sum(wm);                  % number of zero-weight mixtures
            [vv,mk]=sort(lpx);             % find the lowest probability data points
            m=zeros(k,p);                   % initialize means and variances to zero (variances are floored later)
            v=zeros(k,pl);
            m(wm,:)=xs(mk(1:nz),:);                % set zero-weight mixture means to worst-fitted data points
            w(wm)=1/n;                      % set these weights non-zero
            w=w*n/(n+nz);                   % normalize so the weights sum to unity
            wm=~wm;                         % mask for non-zero weights
            m(wm,:)=sx(wm,:)./pk(wm,wp);  % recalculate means and variances for mixtures with a non-zero weight
            v(wm,:)=sx2(wm,:)./pk(wm,wpl);
        end
        v=v-m(:,cix).*m(:,rix);                 % subtract off mean squared
        if g-g1<=th && j>1
            if ~ss, break; end  %  stop
            ss=ss-1;       % stop next time
        end
    end
    if sd  % we need to calculate the final probabilities
        pp=lpx'-0.5*p*log(2*pi)-lsx;   % log of total probability of each data point
        gg=gg(1:j)/nwt-0.5*p*log(2*pi)-lsx;    % average log prob at each iteration
        g=gg(end);
        %             gg' % *** DEBUG ONLY ***
        m=m1;                                           % back up to previous iteration
        v=zeros(p,p,k);                                 % reserve spave for k full covariance matrices
        trv=0;                                          % sum of variance matrix traces
        for ik=1:k                                      % loop for each mixture to apply variance floor
            [uvk,dvk]=eig(reshape(v1(ik,lixi),p,p));	% convert lower triangular to full and find eigenvectors
            dvk=max(diag(dvk),c);                       % apply variance floor to eigenvalues
            v(:,:,ik)=uvk*diag(dvk)*uvk';               % reconstitute full matrix
            trv=trv+sum(dvk);                           % add trace to the sum
        end
        w=w1;
        mm=sum(m,1)/k;
        f=(m(:)'*m(:)-k*mm(:)'*mm(:))/trv;
    else
        v1=v;                                           % lower triangular form
        v=zeros(p,p,k);                                 % reserve spave for k full covariance matrices
        for ik=1:k                                      % loop for each mixture to apply variance floor
            [uvk,dvk,]=eig(reshape(v1(ik,lixi),p,p));	% convert lower triangular to full and find eigenvectors
            dvk=max(diag(dvk),c);                       % apply variance floor
            v(:,:,ik)=uvk*diag(dvk)*uvk';               % reconstitute full matrix
        end
    end
    m=m.*sx0(ones(k,1),:)+mx0(ones(k,1),:);  % unscale means
    v=v.*repmat(sx0'*sx0,[1 1 k]);
end
if l==0         % suppress the first three output arguments if l==0
    m=g;
    v=f;
    w=pp;
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [nodes, edges] = readGraph(fileName)
%READGRAPH Read a graph from a text file.
%
%   [NODES EDGES] = readGraph(FILENAME)
%
%   Example
%     % create a basic graph, save it to a file, and read it again
%     nodes = [10 10;20 10;10 20;20 20;27 15];
%     edges = [1 2;1 3;2 4;2 5;3 4;4 5];
%     writeGraph(nodes, edges, 'simpleGraph.txt');
%     [n2 e2] = readGraph('simpleGraph.txt');
%     figure; drawGraph(n2, e2); axis equal; axis([0 40 0 30]);
%
%   See also
%     writeGraph
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2014-01-21, using Matlab 7.9.0.529 (R2009b)
% Copyright 2014-2022 INRA - Cepia Software Platform

%% Open file and read header

% open file for reading in text mode
f = fopen(fileName, 'rt');
if f == -1
    error(['could not open file for reading: ' fileName]);
end

% check header
line = fgetl(f);
if ~ischar(line)
    error(['can not read graph from empty file: ' fileName]);
end
if ~strcmpi(strtrim(line(2:end)), 'graph')
    error(['Wrong header line in file:' fileName]);
end

%% read node section

% check sub-header
line = fgetl(f);
if ~strcmpi(strtrim(line(2:end)), 'nodes')
    error(['Could not interpret node section in file:' fileName]);
end

% read the number of nodes
line = fgetl(f);
[nNodesStr, line] = strtok(line);
nNodes = str2double(nNodesStr);

% read number of dimension, assumes 2 by default if not specified
nDims = 2;
if ~isempty(line)
    nDims = str2double(strtok(line));
end

% read node coordinates
[nodes, nRead] = fscanf(f, '%g', [nDims nNodes]);
assert(nRead == nNodes * nDims, ...
    'Could not read all node info in file %s', fileName);
nodes = nodes';

% terminate the reading of current line
fgetl(f);

%% read edge section

% check sub-header
line = fgetl(f);
if ~strcmpi(strtrim(line(2:end)), 'edges')
    error(['Could not interpret edge section in file:' fileName]);
end

% read the number of nodes
line = fgetl(f);
nEdgesStr = strtok(line);
nEdges = str2double(nEdgesStr);

% read node indices of each edge
[edges, nRead] = fscanf(f, '%d %d\n', [2 nEdges]);
assert(nRead == nEdges * 2, ...
    'Could not read all edge info in file %s', fileName);
edges = edges';

%% Clean up

fclose(f);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function nonlocalarr = blockmatch(patchmat,neighborindarr,neighbornumarr,selfindarr,para)
%BLOCKMATCH Block matching of patches with non-local similarity.
%   nonlocalarr=BLOCKMATCH(patchmat,neighborindarr,neighbornumarr,
%   selfindarr,para) returns the indexes of the patches with non-local
%   similarity.
%   See slao WNNM_IMDENOISE.
gridsize = length(neighbornumarr);

nonlocalarr = zeros([para.patchnum,gridsize],'uint32');
for i = 1:gridsize
    patch = patchmat(:,selfindarr(i)); % key patch
    neighbors = patchmat(:,neighborindarr(1:neighbornumarr(i),i)); % all neighbors
    distance = sum(bsxfun(@minus,neighbors,patch).^2); % ell_2 distance
    [~,ascendind] = sort(distance); % sort distance in an ascending order
    % indexes of the most similar (shortest distance) para.patchnum of neighbors
    nonlocalarr(:,i) = neighborindarr(ascendind(1:para.patchnum),i);
end

end

repo: DeSCI


text:
% extract_data.m
% --------------
% This program extracts the raw signal data from the CD included in the
% Cumming/Wong book.  The data on CD are in CEOS format.
% It is assumed that the run parameters are stored in CD_run_params.mat
% in the current directory.
% Run "specify_run_parameters.m" first to create this file.
% -------------------------------------------------------------------------
% Created :   Nov 01, 2004  by Kaan Ersahin
% Modified:   Nov 22, 2004  by Kaan Ersahin
% Modified:   Dec 3,  2004  by Ian Cumming
%              - changed function to m-file
%              - fixed the parameter file name
%              - added more radar parameters
% -------------------------------------------------------------------------

clear,   home,   format compact

%  Load the input parameters from a matlab data file
load CD_run_params

disp ' '
disp '------------------------------------------------'
disp ' UBC RRSG - CEOS Reader for RADARSAT-1 RAW DATA'
disp '------------------------------------------------'
disp ' '

% -------------------------------------------------------------------------
% Quantize the range line limits and block size, if necessary
% -------------------------------------------------------------------------

% Move the first range line to the beginning of an 8-line block
first_rg_line = 8 * ( ceil(first_rg_line / 8) - 1 ) + 1;

% Make 'Nrg_lines_blk' a multiple of 8, to get complete the 8-line blocks
Nrg_lines_blk = 8 * ceil(Nrg_lines_blk / 8);

% Find the number of complete blocks required to cover the area of interest
Nblocks = ceil(Nrg_lines / Nrg_lines_blk);

% Make 'Nrg_lines' a multiple of 'Nblocks', to get complete blocks
Nrg_lines = Nrg_lines_blk * Nblocks;

% =========================================================================
% These values are specific to the data set, DO NOT CHANGE for this CD
% =========================================================================

length_replica  =  2880;         % Total length (I&Q) of replica record
tot_Nrg_cells   =  9288;         % Total number of range cells per line
tot_Nrg_lines   = 19432;         % Total number of range lines (records)
first_replica   =     7;         % First record that contains the replica
PRF             = 1256.98;       % Pulse Reputation Frequency (Hz)
Fr              = 32.317e+6;     % Radar sampling rate (Hz)
f0              = 5.300e+9;      % Radar center frequency (Hz)
c               = 2.9979e+8;     % Speed of light (m/s)
R0              = 0.0065956*c/2; % Slant range of first radar sample (m)
Nrepl           = 1349;          % No. of valid samples in the replica
Kr              = 0.72135e+12;   % FM rate of radar pulse (Hz/s)
Tr              = 41.75e-6;      % Chirp duration (s)

% -------------------------------------------------------------------------
% Save parameters in a MAT file that can be used by subsequent programs

if     SaveV6,   save -v6 CD_run_params
else,            save     CD_run_params,      end
% -------------------------------------------------------------------------

fprintf('Total number of range lines      : %5d \n', tot_Nrg_lines )
fprintf('Total number of range cells      : %5d \n', tot_Nrg_cells )
fprintf('First range cell to be extracted : %5d \n', first_rg_cell )
fprintf('First range line to be extracted : %5d \n', first_rg_line )
fprintf('Number of range cells            : %5d \n', Nrg_cells )
fprintf('Number of range lines            : %5d \n', Nrg_lines )
fprintf('Number of range lines per block  : %5d \n', Nrg_lines_blk )
fprintf('Number of blocks                 : %5d \n', Nblocks )
disp ' '
disp '------------------------------------------------'

% -------------------------------------------------------------------------
%  Check the dimensions of the selected area to be processed.
% -------------------------------------------------------------------------

if (first_rg_line <= 0) | ((first_rg_line + Nrg_lines - 1) > tot_Nrg_lines)
    disp ' ',  disp '*****************************************************'
    disp ' ERROR: Check the limits of the range lines !', beep, return
end
if (first_rg_cell <= 0) | ((first_rg_cell + Nrg_cells - 1) > tot_Nrg_cells)
    disp ' ',  disp '*****************************************************'
    disp ' ERROR: Check the limits of the range cells !', beep, return
end

% -------------------------------------------------------------------------
% EXTRACT DATA from the area of interest and write data files
% -------------------------------------------------------------------------

for blk = 1 : Nblocks
    % find the first range line of block 'blk'
    start_line_blk = Nrg_lines_blk * (blk-1) + first_rg_line;
    fprintf('\nExtracting block number%3.0f,  RC1 =%5.0f,  RL1 =%6.0f\n',...
        blk, first_rg_cell, start_line_blk )

    % create the output file name for block 'blk'
    output_file_pre = strcat(output_path,output_prefix,'_',num2str(blk));

    % Call 'read_ceos_raw' function to extract the data for block 'blk'
    read_CEOS_raw( output_file_pre, start_line_blk, blk );

end  % of the 'blk' for loop

beep,   pause(0.3),   beep

author: denkywu

repo: SAR-Synthetic-Aperture-Radar


text:
function OUT = localtopography(DEM,varargin)

%LOCALTOPOGRAPHY Local topography
%
% Syntax
%
%     H = localtopography(DEM)
%     H = localtopography(DEM,radius)
%     H = localtopography(DEM,radius,pn,pv,...)
%
% Description
%
%     localtopography quantifies local relief, e.g. the elevation range
%     within a specific radius. localtopography may take a while to
%     evaluate large DEMs with large kernels. You may speed up calculations
%     by adjusting the parameter 'N' for calculating the max, min or range
%     filter.
%
%     localtopography uses symmetric boundary padding. NaNs
%     are inpainted by nearest neighbor interpolation prior to the
%     calculation (affects only 'mean', 'median', 'prctile', 'std'). For
%     'max', 'min' and 'range', localtopography adopts the padding behavior
%     of the function imdilate and imerode.
%
%
% Input arguments
%
%     DEM    digital elevation model (GRIDobj)
%     radius radius of the moving window filter in map units. The default
%            value is 5000 (m).
%
% Parameter Name/Value (pn,pv) pairs
%
%     'type'   'range' (default), 'max', 'min', 'mean', 'median',
%              'prctile', 'std' (standard deviation)
%     'prc'    scalar between 0 and 100 [%]. Only applicable if 'prctile'
%              is chosen as 'type'.
%     'N'      enable speed improvement for 'max', 'min' and 'range'.
%              N must be 0, 4, 6,or 8. When N is greater than 0, the
%              disk-shaped structuring element is approximated by a sequence
%              of N periodic-line structuring elements. When N equals 0,
%              no approximation is used, and the structuring element members
%              consist of all pixels whose centers are no greater than R
%              away from the origin. If N is not specified, the default
%              value is 0.
%
% Output arguments
%
%     H      local topography grid (GRIDobj)
%
% Example
%
%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
%     H = localtopography(DEM,500);
%     imageschs(DEM,H)
%
% See also: IMDILATE, IMERODE
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 28. January, 2013

narginchk(1,inf)

p = inputParser;
p.FunctionName = 'GRIDobj/localtopography';
expectedTypes = {'range','max','min','mean','median','prctile','std'};

addRequired(p,'DEM',@(x) issparse(x) || isa(x,'GRIDobj'));
addOptional(p,'radius',5000,@(x) isscalar(x) && x>DEM.cellsize);

addParamValue(p,'type','range',@(x) ischar(validatestring(x,expectedTypes)));
addParamValue(p,'N',0,@(x) isscalar(x) && ismember(x,[0 4 6 8]));
addParamValue(p,'thin',1,@(x) x>0.1 && x<=1);
addParamValue(p,'prc',90,@(x) x>0 && x<100);

parse(p,DEM,varargin{:});

dem = DEM.Z;
cs  = DEM.cellsize;

% any nans
INAN = isnan(dem);
flaginan = any(INAN(:));

% structuring element
radiuspx = ceil(p.Results.radius/cs);
SE = strel('disk',radiuspx,p.Results.N);

switch p.Results.type
    case 'max'
        % Maximum filter
        if flaginan;
            dem(INAN) = -inf;
        end
        H = imdilate(dem,SE);
    case 'min'
        % Minimum filter
        if flaginan;
            dem(INAN) = inf;
        end
        H = imerode(dem,SE);
    case 'range'
        if flaginan;
            dem(INAN) = -inf;
        end
        H1 = imdilate(dem,SE);
        if flaginan;
            dem(INAN) = inf;
        end
        H2 = imerode(dem,SE);
        H  = H1-H2;
    case {'mean','average'}
        if flaginan;
            [~,L] = bwdist(~INAN,'e');
            dem = dem(L);
        end
        H   = fspecial('disk',radiuspx);
        H   = imfilter(dem,H,'symmetric','same','conv');
    case 'median'
        if flaginan;
            [~,L] = bwdist(~INAN,'e');
            dem = dem(L);
        end
        H   = getnhood(SE);
        n   = round(sum(H(:))/2);
        H   = ordfilt2(dem,n,H,'symmetric');

    case 'prctile'
        if flaginan;
            [~,L] = bwdist(~INAN,'e');
            dem = dem(L);
        end
        H   = getnhood(SE);
        n   = round(sum(H(:))*p.Results.prc/100);
        H   = ordfilt2(dem,n,H,'symmetric');

    case 'std'
        if flaginan;
            [~,L] = bwdist(~INAN,'e');
            dem = dem(L);
        end
        H   = getnhood(SE);
        H   = stdfilt(dem,H);
end

% Handle nans
if flaginan;
    H(INAN) = nan;
end

% prepare output
OUT = DEM;
OUT.Z = H;
OUT.name = ['local topography (' p.Results.type ')'];

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function Y = subs(S,old,new)
    % Symbolic substitution.  Also used to evaluate expressions numerically.
    %   SUBS(S,OLD,NEW) replaces OLD with NEW in the symbolic expression S.
    %   OLD is a symbolic variable, a string representing a variable name, or
    %   a symbolic expression. NEW is a symbolic or numeric variable
    %   or expression.  That is, SUBS(S,OLD,NEW) evaluates S at OLD = NEW.
    %
    %   SUBS(S,VALUES), where VALUES is a STRUCT, replaces the symbolic
    %   variables in S which are field names in VALUES by the corresponding
    %   entries of the struct.
    %
    %
    %   If OLD and NEW are vectors or cell arrays of the same size, each element
    %   of OLD is replaced by the corresponding element of NEW.  If S and OLD
    %   are scalars and NEW is a vector or cell array, the scalars are expanded
    %   to produce an array result.  If NEW is a cell array of numeric matrices,
    %   the substitutions are performed elementwise.
    %
    %   Examples:
    %     Single input:
    %       Suppose
    %          y = exp(-a*t)*C1
    %       After that, set a = 980 and C1 = 3 in the workspace.
    %       Then the statement
    %          subs(y)
    %       produces
    %          ans = 3*exp(-980*t)
    %
    %     Single Substitution:
    %       subs(a+b,a,4) returns 4+b.
    %       subs(a*b^2, a*b, 5) returns 5*b.
    %
    %     Multiple Substitutions:
    %       subs(cos(a)+sin(b),{a,b},[sym('alpha'),2]) or
    %       subs(cos(a)+sin(b),{a,b},{sym('alpha'),2}) returns
    %       cos(alpha)+sin(2)
    %
    %     Scalar Expansion Case:
    %       subs(exp(a*t),'a',-magic(2)) returns
    %
    %       [   exp(-t), exp(-3*t)]
    %       [ exp(-4*t), exp(-2*t)]
    %
    %     Multiple Scalar Expansion:
    %       subs(x*y,{x,y},{[0 1;-1 0],[1 -1;-2 1]}) returns
    %         0  -1
    %         2   0
    %
    %   See also SYM/SUBEXPR, SYM/VPA, SYM/DOUBLE.

    % Convert inputs to SymExpression
    S = SymExpression(S);

    narginchk(2,3);

    if nargin == 2
        assert(isstruct(old),'SymExpression:invalidInputArgs',...
            'The second argument must be a struct if there are two input arguments.');
        srule = SymExpression(old);

        sstr = ['ReplaceAll[' S.s ',' srule.s ']'];

        % create a new object with the evaluated string
        Y = SymExpression(sstr);
    elseif nargin == 3

        if iscell(old) && iscell(new)
            old_s = cellfun(@(x)SymExpression(x),old,'UniformOutput',false);
            new_s = cellfun(@(x)SymExpression(x),new,'UniformOutput',false);

            rarray = cellfun(@(x,y)[x.s '->' y.s], old_s, new_s ,'UniformOutput',false);
            srule = SymExpression(['<| ', ...
                implode(rarray,', '), ...
                ' |>']);
        elseif ischar(old)
            old_s = SymExpression(old);
            new_s = SymExpression(new);
            srule = SymExpression(['<| ', ...
                old_s.s, '->', new_s.s, ...
                ' |>']);
        elseif isa(old,'SymExpression')
            siz_o = size(old);
            siz_n = size(new);

            assert(prod(siz_o) == prod(siz_n),...
                'The sizes of the second and third argument must be the same.');
            old_s = flatten(SymExpression(old));
            new_s = flatten(SymExpression(new));
            rarray = arrayfun(@(x,y)['First@' x.s '->First@' y.s], old_s, new_s ,'UniformOutput',false);
            srule = SymExpression(['<| ', ...
                implode(rarray,', '), ...
                ' |>']);
        end

        sstr = ['ReplaceAll[' S.s ',' srule.s ']'];

        % create a new object with the evaluated string
        Y = SymExpression(sstr);
    end

author: ayonga

repo: frost-dev

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


text:
function [grade,X,SR,improve,AppSet] = LocalSearch2(Problem,X,SR,improve,AppSet)
% Local Search 2 of MTS

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

    if ~improve
        SR = SR / 2;
        if all(SR<1e-8)
            SR = (Problem.upper-Problem.lower).*(rand(1,length(SR))/10+0.4);
        end
    end
    improve = false;
    grade   = 0;
    for l = 1 : length(SR)
        chosen = rand(1,length(SR)) < 1/4;
        old_X  = X;
        dec    = X.dec;
        dec(chosen) = dec(chosen) + SR(chosen).*(rand(1,sum(chosen))*2-1);
        X      = Problem.Evaluation(dec);
        [grade,improve,AppSet] = Grading(X,old_X,grade,improve,AppSet,Problem.N);
        if all(old_X.obj<=X.obj)
            dec = old_X.dec;
            dec(chosen) = dec(chosen) - 0.5*SR(chosen).*(rand(1,sum(chosen))*2-1);
            X   = Problem.Evaluation(dec);
            [grade,improve,AppSet] = Grading(X,old_X,grade,improve,AppSet,Problem.N);
            if all(old_X.obj<=X.obj)
                X = old_X;
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%%*******************************************************************
%% schurmat_qblk: compute schur matrix corresponding to SOCP blocks.
%%
%% HKM direction: output = schur + Ax*Ae' + Ae*Ax' - Ad*Ad'
%% NT  direction: output = schur + Ae*Ae' - Ad*Ad'
%%
%% where schur = A*D*A', and Ad is the modification to ADA'
%% so that the latter is positive definite.
%%
%% [schur,UU,EE] = schurmat_qblk(blk,At,schur,UU,EE,p,dd,ee,xx);
%%
%% UU: stores the dense columns of Ax, Ae, Ad, and possibly
%%     those of A*D^{1/2}. It has the form UU = [Ax Ae Ad].
%% EE: stores the assocaited (2,2) block matrix when the
%%     output matrix is expressed as an augmented matrix.
%%     It has the form EE = [0 -lam 0; -lam 0 0; 0 0 I].
%%
%% options = 0, HKM
%%         = 1, NT
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 10 Jun 02
%%*******************************************************************

   function [schur,UU,VV] = schurmat_qblk(blk,At,schur,UU,VV,p,dd,ee,xx);

   if (nargin == 9); options = 0; else; options = 1; end;

   pblk = blk(p,:); n = sum(pblk{2});  numblk = length(pblk{2});

   Ae = qprod(pblk,At{p}',ee{p});
   if (options == 0)
      Ax = qprod(pblk,At{p}',xx{p});
   end;
   decolidx = checkdense(Ae);
   ddsch = dd{p};
   if ~isempty(decolidx);
      spcolidx = setdiff([1:numblk],decolidx);
      s = 1 + [0 cumsum(pblk{2})];
      idx = s(decolidx);
      tmp = zeros(n,1);
      tmp(idx) = sqrt(2*abs(ddsch(idx)));
      Ad = qprod(pblk,At{p}',tmp);
      ddsch(idx) = abs(ddsch(idx));
      if (options == 0)
         UU = [UU Ax(:,decolidx) Ae(:,decolidx)  Ad];
         VV = [VV Ae(:,decolidx) Ax(:,decolidx) -Ad];
         Ax = Ax(:,spcolidx); Ae = Ae(:,spcolidx);
         schur = schur + Ax*Ae'+ Ae*Ax';
      else
         UU = [UU Ae(:,decolidx)  Ad];
         VV = [VV Ae(:,decolidx) -Ad];
         Ae = Ae(:,spcolidx);
         schur = schur + Ae*Ae';
      end
   else
      if (options == 0)
         schur = schur + Ax*Ae'+ Ae*Ax';
      else
         schur = schur + Ae*Ae';
      end
   end
   decolidx = checkdense(At{p}');
   if ~isempty(decolidx);
      len = length(decolidx);
      tmp = (spdiags(sqrt(abs(ddsch(decolidx))),0,len,len)*At{p}(decolidx,:))';
      UU = [UU tmp];
      VV = [VV tmp*spdiags(sign(ddsch(decolidx)),0,len,len)];
      ddsch(decolidx) = zeros(len,1);
   end
   schur = schur + At{p}' *spdiags(ddsch,0,n,n) *At{p};
%%*******************************************************************

author: yu-jiang

repo: radpbook

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


text:
% SBPLOT - create axes in arbitrary subplot grid positions and sizes
%
% Usage:  >> axis_handle = sbplot(v,h,index)
%         >> axis_handle = sbplot(v,h,[index1 index2])
%         >> axis_handle = sbplot(v,h,[index1 index2],axprop,..)
%         >> axis_handle = sbplot(v,h,[index1 index2],'ax',handle,axprop,..)
%
% Inputs:
%   v,h    - Integers giving the vertical and horizontal ranks of the tiling.
%   index  - Either a single subplot index, in which case the command
%            is equivalent to subplot, or a two-element vector giving
%            the indices of two corners of the SBPLOT area according
%            to SUBPLOT convention (e.g., left-to-right, top-to-bottom).
%   axprop - Any axes property(s), e.g., >> sbplot(3,3,3,'color','w')
%   handle - Following keyword 'ax', sbplot tiles the given axes handle
%            instead of the whole figure
%
% Output:
%  axis_handle - matlab axis handle
%
% Note:
%   sbplot is essentially the same as the subplot command except that
%   sbplot axes may span multiple tiles. Also, SBPLOT will not erase
%   underlying axes.
%
%  Examples: >> sbplot(3,3,6);plot(rand(1,10),'g');
%            >> sbplot(3,3,[7 2]);plot(rand(1,10),'r');
%            >> sbplot(8,7,47);plot(rand(1,10),'b');
%
% Authors: Colin Humphries, Arnaud Delorme & Scott Makeig, SCCN/INC/UCSD, La Jolla, June, 1998

% Copyright (C) June 1998, Colin Humphries & Scott Makeig, SCCN/INC/UCSD,
% scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% reformatted by Scott Makeig, 6/10/98
% 12/22/00 test nargin<3 -sm
% 01/21/01 added (recursive) axes option 'ax' -sm
% 01-25-02 reformated help & licence -ad

function [out] = sbplot(m,n,gridpos,varargin) % varargin is std. matlab arg list

if nargin<3
   error('            requires >=3 arguments');
end

if nargin>3 && strcmp(varargin{1},'ax')
   pos = get(varargin{2},'Position'); %  sbplot(3,1,[2 3]) -> 0.4111 0.1100 0.4939 0.815
   varargin = {varargin{3:end}};
else
   pos = get(gcf,'DefaultAxesPosition'); %  [0.1300    0.1100    0.7750    0.815]
end

Xpad = pos(1);      % lower-left distance from left side of figure
Ypad = pos(2);      % lower-left distance from bottom of figure

Xlen = pos(3);      % axes width
Ylen = pos(4);      % axes height

if n == 2
  xspace = Xlen*0.27/(n-0.27);      % xspace between axes as per subplot
else
  xspace = (0.9*Xlen)*0.27/(n-0.9*0.27);
end

if m == 2
  yspace = Ylen*0.27/(m-0.27);      % yspace between axes as per subplot
else                                % WHY Xlen (.775) instead of Ylen (.815) ??
  yspace = (0.9*Ylen)*0.27/(m-0.9*0.27);
end

xlength = (Xlen-xspace*(n-1))/n;  % axes width
ylength = (Ylen-yspace*(m-1))/m;  % axes height

% Convert tile indices to grid positions
if length(gridpos) == 1
  xgridpos(1) = mod(gridpos,n);   % grid position
  if xgridpos(1) == 0
    xgridpos(1) = n;
  end
  xgridpos(2) = 1;                   % grid length
  ygridpos(1) = m-ceil(gridpos/n)+1; % grid position
  ygridpos(2) = 1;                   % grid length
else
  xgridpos(1) = mod(gridpos(1),n);
  if xgridpos(1) == 0
    xgridpos(1) = n;
  end
  tmp = mod(gridpos(2),n);
  if tmp == 0
    tmp = n;
  end
  if tmp > xgridpos(1)
    xgridpos(2) = tmp-xgridpos(1)+1;
  else
    xgridpos(2) = xgridpos(1)-tmp+1;
    xgridpos(1) = tmp;
  end

  ygridpos(1) = m-ceil(gridpos(1)/n)+1;
  tmp = m-ceil(gridpos(2)/n)+1;
  if tmp > ygridpos(1)
    ygridpos(2) = tmp-ygridpos(1)+1;
  else
    ygridpos(2) = ygridpos(1)-tmp+1;
    ygridpos(1) = tmp;
  end
end

% Calculate axes coordinates
position(1) = Xpad+xspace*(xgridpos(1)-1)+xlength*(xgridpos(1)-1);
position(2) = Ypad+yspace*(ygridpos(1)-1)+ylength*(ygridpos(1)-1)-0.03;
position(3) = xspace*(xgridpos(2)-1)+xlength*xgridpos(2);
position(4) = yspace*(ygridpos(2)-1)+ylength*ygridpos(2);

% Create new axes
ax = axes('Position',position,varargin{:});

% Output axes handle
if nargout > 0
  out = ax;
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function P = subGFM(PopObj,Center,R,FrontNo)
% PF modeling for each subregion

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

    K     = size(Center,1);
    [N,M] = size(PopObj);

    % Normalize the population
    fmin = min(PopObj,[],1);
    fmax = max(PopObj,[],1);
    Obj  = (PopObj-repmat(fmin,N,1))./repmat(fmax-fmin,N,1);

    % PF modeling
    if K == 1
        P = GFM(Obj(FrontNo==1,:));
    else
        P = ones(K,M);

        % Allocation
        transformation = Allocation(Obj,Center,R);
        subFirstFront  = false(N,1);

        % Non-dominated sorting od each subregion
        for i = 1 : K
            current = find(transformation == i);
            if  ~isempty(current)
                [FNo,MFNo] = NDSort(PopObj(current,:),length(current));
                subFirstFront(current(FNo<MFNo|FNo==1)) = true;
            end
        end

        FTransformation = transformation(subFirstFront);
        PopObj          = PopObj(subFirstFront,:);
        RemainObj       = Obj(subFirstFront,:);

        % GFM of each subregion
        if size(PopObj,1) > M
           for i = 1 : K
               current = find(FTransformation==i);
               if ~isempty(current)
                    if length(current) < M+1
                        [~,sDis] = sort(pdist2(RemainObj ,Center(i,:)));
                        current  = sDis(1:M+1);
                    end
                    p  = GFM(Obj(current,:));
                    P(i,:) = p;
               end
           end
        end

    end

end

function P = GFM(X)
% Generic front modeling

    [N,M] = size(X);
    X     = max(X,1e-12);
    P     = ones(1,M);
    lamda = 1;
    E     = sum(X.^repmat(P,N,1),2) - 1;
    MSE   = mean(E.^2);
    for epoch = 1 : 1000
        % Calculate the Jacobian matrix
        J = X.^repmat(P,N,1).*log(X);
        % Update the value of each weight
        while true
            Delta  = -(J'*J+lamda*eye(size(J,2)))^-1*J'*E;
            newP   = P + Delta(1:end)';
            newE   = sum(X.^repmat(newP,N,1),2) - 1;
            newMSE = mean(newE.^2);
            if newMSE < MSE && all(newP>1e-3)
                P     = newP;
                E     = newE;
                MSE    = newMSE;
                lamda = lamda/1.08;
                break;
            elseif lamda > 1e8
                return;
            else
                lamda = lamda*1.08;
            end
        end
    end
end


text:
% Test file for trigtech/sign.m

function pass = test_sign(pref)

if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

% Test a positive function:
F = @(x) sin(pi*x) + 2;
f = testclass.make(@(x) F(x), [], pref);
h = sign(f);
pass(1) = normest(h - 1) < 10*eps;

% Test a negative function:
f2 = testclass.make(@(x) -F(x), [], pref);
h = sign(f2);
pass(2) = normest(h + 1) < 10*eps;

% Test a complex-valued function:
F = @(x) exp(1i*pi*x);
f = testclass.make(@(x) F(x), [], pref);
h = sign(f);
pass(3) = normest(h - f) < 10*eps;

% Test a complex array-valued function:
xx = linspace(-.95, .97);
F = @(x) [(2+sin(pi*x)).*exp(1i*pi*x), -(2+sin(pi*x)).*exp(1i*pi*x), 2+sin(pi*x)];
f = testclass.make(@(x) F(x), [], pref);
ff = feval(f, xx);
gg = ff./abs(ff);
h = sign(f);
hh = feval(h, xx);
pass(4) = norm(hh - gg, inf) < 10*eps;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Helper scriprt for running the DC Motor Simulation
close all % figures
disp('Running Simulation ...')
tic
plotMotorSim([Voltage, Inertia],0,'Starting Simulation...');

% perform simulation
for i = 1:noRuns
    V = Voltage(i);     % input voltage         (V)
    Jd = Inertia(i);    % load inertia          (g-cm^2)
    sim('MaxonDCMotor.mdl',[0 20]);
    time{i} = tout;
    velocity{i} = vel;
    plotMotorSim([Voltage, Inertia],i,['Run ',num2str(i), ' of ',num2str(noRuns)])
end
clear NbBlocks blk m powerguiblock vel tout t s simData
plotMotorSim([Voltage, Inertia],i,'Finished Simulation')
toc

author: Sable

repo: mcbench-benchmarks


text: Mac OS X        	2 RMATFMATL

repo: Aero-matlab


text:
filename='Cantilever_hexahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.1;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
Network torchvision.models.squeezenet {
Layer Conv2d-1 {
Type: CONV
Stride { X: 2, Y: 2 }
Dimensions { K: 96, C: 3, R: 7, S: 7, Y: 224, X: 224 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-2 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 16, C: 96, R: 1, S: 1, Y: 54, X: 54 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-3 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 64, C: 16, R: 1, S: 1, Y: 54, X: 54 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-4 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 64, C: 16, R: 3, S: 3, Y: 54, X: 54 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-5 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 16, C: 128, R: 1, S: 1, Y: 54, X: 54 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-6 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 64, C: 16, R: 1, S: 1, Y: 54, X: 54 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-7 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 64, C: 16, R: 3, S: 3, Y: 54, X: 54 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-8 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 32, C: 128, R: 1, S: 1, Y: 54, X: 54 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-9 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 128, C: 32, R: 1, S: 1, Y: 54, X: 54 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-10 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 128, C: 32, R: 3, S: 3, Y: 54, X: 54 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-11 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 32, C: 256, R: 1, S: 1, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-12 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 128, C: 32, R: 1, S: 1, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-13 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 128, C: 32, R: 3, S: 3, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-14 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 48, C: 256, R: 1, S: 1, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-15 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 192, C: 48, R: 1, S: 1, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-16 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 192, C: 48, R: 3, S: 3, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-17 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 48, C: 384, R: 1, S: 1, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-18 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 192, C: 48, R: 1, S: 1, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-19 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 192, C: 48, R: 3, S: 3, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-20 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 64, C: 384, R: 1, S: 1, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-21 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 256, C: 64, R: 1, S: 1, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-22 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 256, C: 64, R: 3, S: 3, Y: 27, X: 27 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-23 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 64, C: 512, R: 1, S: 1, Y: 13, X: 13 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-24 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 256, C: 64, R: 1, S: 1, Y: 13, X: 13 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-25 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 256, C: 64, R: 3, S: 3, Y: 13, X: 13 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
Layer Conv2d-26 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1000, C: 512, R: 1, S: 1, Y: 13, X: 13 }
Dataflow {
        SpatialMap(1,1) K;
        TemporalMap(64,64) C;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        Cluster(64, P);
        SpatialMap(1,1) C;
        TemporalMap(Sz(R),1) Y;
        TemporalMap(Sz(S),1) X;
        TemporalMap(Sz(R),Sz(R)) R;
        TemporalMap(Sz(S),Sz(S)) S;
}
}
}

author: maestro-project

repo: maestro


text:
function [uGal] = cmps22uGal(cmps2)
% Convert acceleration from centimeters per square centimeter to microgalileos
% Chad A. Greene 2012
uGal = cmps2*1e+6;
end


text:
function doseProfileFigure(command, varargin)
%"doseProfileFigure"
%   Given a start point and end point (startPt, endPt) as [x,y,z]
%   coordinates, create a figure containing a plot doses along the line
%   connecting startPt to endPt.  This figure has controls to activate
%   plots of any or all existing CT scans or doses along the same line, but
%   defaults to the dose passed in initialDoseV and initialScanV.
%
%JRA 1/11/05
%AI  80/04/16  Added scrolling/pagination to handle large number of scans and dose distributions.
%Usage:
%   function doseProfileFigure('init', initialDoseV, initialScanV);
%   function doseProfileFigure('new_points', startPt, endPt);
%   function doseProfileFigure('refresh');
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC
global stateS
indexS = planC{end};
%persistent drawing;

initialDoseV = [];
initialScanV = [];

if isnumeric(command) && length(command) == 3
    startPt = command;
    endPt = varargin{1};
    command = 'init';
elseif ischar(command) && strcmpi(command, 'init')
    startPt = [];
    endPt = [];
    if nargin > 2
        initialDoseV = varargin{1};
        initialScanV = varargin{2};
    end
elseif ischar(command)
else
    error('Invalid call to doseProfileFigure. Help doseProfileFigure for details.');
end

switch upper(command)
    case 'INIT'
        stateS.doseDiffScale = 'Abs';
        drawing = 0;
        %oldFigs = findobj('tag', 'CERR_DoseLineProfile');
        if isfield(stateS.handle,'doseProfileFigure') && ...
                ~isempty(ishandle(stateS.handle.doseProfileFigure)) && ...
                ishandle(stateS.handle.doseProfileFigure)
            delete(stateS.handle.doseProfileFigure);
        end

        units = 'pixels';
        screenSize = get(0,'ScreenSize');
        % APA: Commented
        % w = 500; h = 500;
        % APA: Commented
        w = 510; h = 650;

        %Initial size of figure in pixels. Figure scales fairly well.
        hFig = figure('name', 'Dose Line Profile', 'units', units, 'position',[(screenSize(3)-w) 35 w h], 'MenuBar', 'none', 'NumberTitle', 'off', 'resize', 'off', 'Tag', 'CERR_DoseLineProfile', 'DeleteFcn', 'doseProfileFigure(''CLOSE'')');
        stateS.handle.doseProfileFigure = hFig;

        nDoses = length(planC{indexS.dose});
        nScans = length(planC{indexS.scan});

        %nScans = min(nScans,20); % 20 scans max allowed. Add pagination in future.

        %Draw frames.
        yFrameSize = 110;
        yStart = 100;
        hDoseFrame = uipanel(hFig, 'units', units, 'position', [10 10 235 yFrameSize]);
        hScanFrame = uipanel(hFig, 'units', units, 'position', [255 10 235 yFrameSize]);
        uicontrol(hFig, 'style', 'text', 'units', units, 'position', [10 yFrameSize-10 50 20], 'string', 'Doses', 'fontweight', 'bold');
        uicontrol(hFig, 'style', 'text', 'units', units, 'position', [255 yFrameSize-10 50 20], 'string', 'Scans', 'fontweight', 'bold');

        %AI: Changed
        %Create dose & scan uicontrols.
        doseUI = uitable(hFig,'RowName',[],'ColumnName',[],'Position',[20 20 215 80],'visible','off','Tag','doseUI');
        scanUI = uitable(hFig,'RowName',[],'ColumnName',[],'Position',[265 20 215 80],'visible','off','tag','scanUI');

        %Get dose list
        colors = stateS.optS.colorOrder;
        outString = cell(nDoses+1,1);
        outString{1} = 'Pick Reference Dose';
        chkvaldoseC = cell(nDoses,1);
        chkvaldoseC(:) = {false};
        doseStringC = cell(nDoses,1);
        colorDose3M = zeros(nDoses,3);
        for i=1:nDoses
            if ismember(i, initialDoseV)
                chkvaldoseC{i} = true;
            else
                chkvaldoseC{i} = false;
            end
            outString{i+1} = planC{indexS.dose}(i).fractionGroupID;
            colorDose3M(i,:) = getColor(i, colors, 'loop');
            if ~isLocal(planC{indexS.dose}(i).doseArray)
                doseStringC{i} = [num2str(i) '. ' planC{indexS.dose}(i).fractionGroupID,'(R)'];
            else
                doseStringC{i} = [num2str(i) '. ' planC{indexS.dose}(i).fractionGroupID];
            end
        end

        %Get scan list
        chkvalscanC = cell(nScans,1);
        chkvalscanC(:) = {false};
        scanStringC = cell(nScans,1);
        colorScan3M = zeros(nDoses,3);
        for i=1:nScans
            if ismember(i, initialScanV);
                chkvalscanC{i} = true;
            else
                chkvalscanC{i} = false;
            end

            colorScan3M(i,:) = getColor(i+nDoses, colors, 'loop');
            if ~isLocal(planC{indexS.scan}(i).scanArray)
                scanStringC{i} = [num2str(i) '. ' planC{indexS.scan}(i).scanType,'(R)'];
            else
                scanStringC{i} = [num2str(i) '. ' planC{indexS.scan}(i).scanType];
            end
        end

        %Display dose uicontrols.
        doseDataC = cat(2,chkvaldoseC,doseStringC);
        set(doseUI, 'data', doseDataC, 'visible', 'on','BackgroundColor',colorDose3M,'cellSelectionCallback',@objCheck);

        %Display scan uicontrols.
        scanDataC = cat(2,chkvalscanC,scanStringC);
        set(scanUI, 'data', scanDataC,'visible', 'on','BackgroundColor',colorScan3M,'cellSelectionCallback',@objCheck);

        %AI: End changed

        %Store UI handles and other info for later use.
        ud.doseUI = doseUI;
        ud.baseDose = [];
        ud.scanUI = scanUI;
        ud.nDoses = nDoses;
        ud.nScans = nScans;
        ud.startPt = startPt;
        ud.endPt = endPt;
        ud.hDoseProfPlotV = [];

        % APA: Commented
        %         %Create axes.
        %         ud.doseaxis = axes('units', 'pixels', 'position', [50 yStart+100 400 h-yStart-150], 'parent', hFig, 'tickdir', 'out');
        %         xlabel('Distance along line');
        %         ylabel('Dose Value');
        %         ud.scanaxis = axes('units', 'pixels', 'position', [50 yStart+100 400 h-yStart-150], 'parent', hFig, 'tickdir', 'out', 'Color','none','YAxisLocation','right', 'xtick', []);
        %         ylabel('Scan Value');
        % APA: Commented ends

        % APA: try
        %Create axes.
        % ud.doseaxis = axes('units', 'pixels', 'position', [50 yStart+100+170 400 h-yStart-330], 'parent', hFig, 'tickdir', 'out','fontSize',8);
        ud.doseaxis = axes('units', 'pixels', 'position', [50 yStart+270+25 400 h-yStart-300], 'parent', hFig, 'tickdir', 'out','fontSize',8);
        % xlabel('\bfDistance along line','fontSize',10);
        xlabel({});
        ylabel('\bfDose Value','fontSize',10);
        % ud.scanaxis = axes('units', 'pixels', 'position', [50 yStart+100+170 400 h-yStart-330], 'parent', hFig, 'tickdir', 'out', 'Color','none','YAxisLocation','right', 'xtick', [],'fontSize',8);
        ud.scanaxis = axes('units', 'pixels', 'position', [50 yStart+270+25 400 h-yStart-300], 'parent', hFig, 'tickdir', 'out', 'Color','none','YAxisLocation','right', 'xtick', [],'fontSize',8);
        ylabel('\bfScan Value','fontSize',10);

        % APA: try ends

        % DK Radio button for selecting absolute or normalized scale
        clrBg = get(hFig,'color');

        radioToolTip = 'Choose Absolute(ABS) or Normalized (NOR) Scale for Dose Difference';
        h = uibuttongroup('parent',hFig,'units', 'pixels', 'visible','off','Position',[60 yStart+240 100 25]);
        u0 = uicontrol('Style','radiobutton','String','Abs',...
            'position',[5 5 40 15],'parent',h,'HandleVisibility','off');
        u1 = uicontrol('Style','radiobutton','String','Rel',...
            'position',[50 5 40 15],'parent',h,'HandleVisibility','off');
        set(h,'selectionchangedfcn',@d_ct_profRadio);
        set(h,'SelectedObject',u0);  % No selection
        set(h,'Visible','on')

        % DK Drop down list to select the primary dose
        uicontrol('parent',hFig,'Style','Text','Units','Pixels','Position',[170 yStart+230 120 30],...
            'String','Select Base Dose','backgroundcolor',clrBg);

        refDoseToolTip = 'Select the reference/Base Dose from which all doses will be subtracted';

        ud.refDoseSelect = uicontrol('parent',hFig,'Style','popupmenu', 'Units','Pixels','Position',[300 yStart+220 100 40],'Tag','refDoseSelect',...
            'callback', 'doseProfileFigure(''selectrefDose'')','String',outString,'TooltipString',refDoseToolTip);

        % APA: Create  difference axis
        % ud.diffaxis = axes('units', 'pixels', 'position', [50 yStart+90 400 yStart+100+5], 'parent', hFig, 'tickdir', 'out','fontSize',8);
        ud.diffaxis = axes('units', 'pixels', 'position', [50 yStart+100 400 120], 'parent', hFig, 'tickdir', 'out','fontSize',8);
        xlabel('\bfDistance along line','fontSize',10);
        ylabel('\bfDifference','fontSize',10);
        ud.htext = [];

        % APA: Create  slider axis
        ud.slideraxis = axes('units', 'pixels', 'position', [40 yStart+220 420 15], 'parent', hFig,'visible','off');
        %         udS.xL = patch([0.025 0 0], [0.5 0 1], 'w', 'edgecolor', 'k', 'buttondownfcn', 'doseProfileFigure(''RANGERCLICKED'', ''xL'')', 'userdata', ud.slideraxis, 'parent', ud.slideraxis, 'erasemode', 'xor');
        %         udS.xR = patch([0.975 1 1], [0.5 0 1], 'w', 'edgecolor', 'k', 'buttondownfcn', 'doseProfileFigure(''RANGERCLICKED'', ''xR'')', 'userdata', ud.slideraxis, 'parent', ud.slideraxis, 'erasemode', 'xor');
        udS.xL = patch([0.025 0.025 0], [0 1 1], 'w', 'edgecolor', 'k', 'buttondownfcn', 'doseProfileFigure(''RANGERCLICKED'', ''xL'')', 'userdata', ud.slideraxis, 'parent', ud.slideraxis);
        udS.xR = patch([0.975 0.975 1], [0 1 1], 'w', 'edgecolor', 'k', 'buttondownfcn', 'doseProfileFigure(''RANGERCLICKED'', ''xR'')', 'userdata', ud.slideraxis, 'parent', ud.slideraxis);

        % Display initial text on difference axis
        grid(ud.diffaxis,'off')
        set(ud.diffaxis,'xLim',[0 1]);
        set(ud.diffaxis,'yLim',[0 1]);
        htextInit = text(0.1,0.5,'Select More than one dose to display difference','parent', ud.diffaxis,'HorizontalAlignment','left','FontWeight','bold');
        ud.htextInit = htextInit;
        set(hFig, 'userdata',ud);

        %Save userdata.
        set(hFig, 'userdata', ud);
        set(ud.slideraxis, 'userdata', udS);

    case 'SELECTREFDOSE'

        %hFig = findobj('tag', 'CERR_DoseLineProfile');
        hFig = stateS.handle.doseProfileFigure;

        ud = get(hFig, 'userdata');

        %doseUI  = ud.doseUI;

        value = get(ud.refDoseSelect,'value')-1;

        if value ~=0
            ud.baseDose = value;
            %set(doseUI(ud.baseDose), 'FontWeight','bold');
        end

        set(hFig,'Userdata',ud);

    case 'RANGERCLICKED'
        %hFig = findobj('tag', 'CERR_DoseLineProfile');
        hFig = stateS.handle.doseProfileFigure;
        ud    = get(hFig, 'userdata');
        udS   = get(ud.slideraxis, 'userdata');

        %Suspend figure while motion is going on, restore later.
        udS.UISTATE = uisuspend(hFig);

        %Determine which triangle is clicked.
        switch upper(varargin{1})
            case 'XL'
                udS.handle.movingObject = udS.xL;
                RangeMax = get(udS.xR, 'xData');
                udS.motionbounds = [0.025 RangeMax(1)];
            case 'XR'
                udS.handle.movingObject = udS.xR;
                RangeMin = get(udS.xL, 'xData');
                udS.motionbounds = [RangeMin(1) 0.975];
        end
        %Prepare object/figure for motion.
        %set(udS.handle.movingObject, 'erasemode', 'xor');
        set(hFig, 'windowbuttonmotionfcn', 'doseProfileFigure(''INDICATORMOVING'')')
        set(hFig, 'windowbuttonupfcn', 'doseProfileFigure(''MOTIONDONE'');')
        set(ud.slideraxis, 'userdata', udS);

    case 'INDICATORMOVING'
        %hFig = findobj('tag', 'CERR_DoseLineProfile');
        hFig = stateS.handle.doseProfileFigure;
        ud    = get(hFig, 'userdata');
        hAxis = ud.slideraxis;
        udS   = get(ud.slideraxis, 'userdata');
        hObj  = udS.handle.movingObject;
        xData = get(hObj, 'xData');
        hAxis = ud.slideraxis;
        cp = get(hAxis, 'currentpoint');
        mB = udS.motionbounds;
        if cp(2,1) > udS.motionbounds(2)
            cp(2,1) = udS.motionbounds(2);
        elseif cp(2,1) < udS.motionbounds(1)
            cp(2,1) = udS.motionbounds(1);
        end
        delta = cp(2,1) - xData(1);
        %Move both the object and text string to new position.
        set(hObj, 'xData', xData+delta);
        doseProfileFigure('refresh');

    case 'MOTIONDONE'
        %hFig = findobj('tag', 'CERR_DoseLineProfile');
        hFig = stateS.handle.doseProfileFigure;
        ud    = get(hFig, 'userdata');
        udS   = get(ud.slideraxis, 'userdata');

        %Restore the figure to pre-motion state.
        uirestore(udS.UISTATE);

        %Match position of arrow with rounded text value.
        %         hObj  = udS.handle.movingObject;
        %         xData = get(hObj, 'xData');
        %         cP    = str2num(get(ud.associatedText, 'string'));
        %         delta = cP - xData(1);
        %         set(hObj, 'xData', xData+delta);

        %set(udS.handle.movingObject, 'erasemode', 'xor');
        doseProfileFigure('refresh');

    %AI : Removed 'DOSE_CHECK','SCAN_CHECK'

    case 'NEW_POINTS'
        %New points specified for the profile line.
        %hFig = findobj('tag', 'CERR_DoseLineProfile');
        hFig = stateS.handle.doseProfileFigure;
        ud = get(hFig, 'userdata');
        ud.startPt = varargin{1};
        ud.endPt = varargin{2};
        set(hFig, 'userdata', ud);

    case 'REFRESH'
        %Make sure only the most recent calls are drawn.
        %global CERR_doseProfileLastCallTime;
        %thisCallTime = now;
        %if isempty(CERR_doseProfileLastCallTime)
        %    CERR_doseProfileLastCallTime = thisCallTime;
        %else
        %    if thisCallTime > CERR_doseProfileLastCallTime
        %        CERR_doseProfileLastCallTime = thisCallTime;
        %    end
        %end
        %clear CERR_doseProfileLastCallTime;

        %Make dose profile figure the foreground.
        %hFig = findobj('tag', 'CERR_DoseLineProfile');
        hFig = stateS.handle.doseProfileFigure;
        if isempty(hFig)
            return;
        end
        ud = get(hFig, 'userdata');
        udS   = get(ud.slideraxis, 'userdata');
        xDataL = get(udS.xL,'xData');
        x1 = xDataL(1);
        xDataR = get(udS.xR,'xData');
        x2 = xDataR(1);

        delete(ud.htext)
        ud.htext = [];

        if ~isequal(gcf,hFig)
            figure(hFig);
        end

        %delete(findobj(hFig, 'Tag', 'CERR_DOSEPROFILE_PLOT'));
        delete(ud.hDoseProfPlotV)

        nSamples = stateS.optS.numDoseProfileSamples;
        ptsDel = 1:nSamples;
        ptsL = nSamples/0.95*(x1-0.025);
        ptsR = nSamples/0.95*(x2-0.025);
        ptsDel = unique([1:floor(ptsL) ceil(ptsR+1e-6):nSamples]);

        %Get UD stored variables.
        doseUI  = ud.doseUI;
        baseDose = ud.baseDose;
        scanUI  = ud.scanUI;
        nDoses  = ud.nDoses;
        nScans  = ud.nScans;
        startPt = ud.startPt;
        endPt   = ud.endPt;

        if isempty(startPt) || isempty(endPt)
            return;
        end

        colors = stateS.optS.colorOrder;
        %AI: Changed
        drawDoses = get(doseUI,'Data');
        drawDoses = [drawDoses{:,1}];
        drawScans = get(scanUI,'Data');
        drawScans = [drawScans{:,1}];
        %AI: End changed

        drawDoses = find(drawDoses);
        drawScans = find(drawScans);

        %Get x,y,z coords of samples and distance from first pt.
        xV = linspace(startPt(1), endPt(1), nSamples);
        yV = linspace(startPt(2), endPt(2), nSamples);
        zV = linspace(startPt(3), endPt(3), nSamples);
        %distV = sqrt(sepsq(startPt', [xV;yV;zV]));
        len = sqrt(sum((startPt - endPt).^2));
        distV = 0:len/(nSamples-1):len;
        %max(abs(distV - dist1V))

        for i=1:length(drawDoses)
            set(ud.doseaxis, 'nextplot', 'add', 'tickdir', 'out');
            doseNum = drawDoses(i);
            transM = getTransM('dose', doseNum, planC);
            if ~isempty(transM)
                [xV1, yV1, zV1] = applyTransM(inv(transM), xV, yV, zV);
            else
                xV1 = xV; yV1 = yV; zV1 = zV;
            end
            dV{i} = getDoseAt(doseNum, xV1, yV1, zV1, planC);
            dV{i}(ptsDel) = NaN;
        end

        for i=1:length(drawScans)
            set(ud.scanaxis, 'nextplot', 'add', 'tickdir', 'out');
            scanNum = drawScans(i);
            transM = getTransM('scan', scanNum, planC);
            if ~isempty(transM)
                [xV1, yV1, zV1] = applyTransM(inv(transM), xV, yV, zV);
            else
                xV1 = xV; yV1 = yV; zV1 = zV;
            end
            scanV = getScanAt(scanNum, xV1, yV1, zV1, planC);
            offset = planC{indexS.scan}(scanNum).scanInfo(1).CTOffset;
            sV{i} = single(scanV) - offset;
            sV{i}(ptsDel) = NaN;
        end

        %persistent CERR_doseProfileLastCallTime
        hDoseProfPlotV = [];
        %if isequal(CERR_doseProfileLastCallTime, thisCallTime)
            for i=1:length(drawDoses)
                doseNum = drawDoses(i);
                hDoseProfPlotV(end+1) = plot(distV, dV{i}, 'parent', ud.doseaxis, ...
                    'color', getColor(doseNum, colors, 'loop'),...
                    'Tag', 'CERR_DOSEPROFILE_PLOT','linewidth',2);
            end
            axis(ud.doseaxis,'tight')
            for i=1:length(drawScans)
                scanNum = drawScans(i);
                hDoseProfPlotV(end+1) = plot(distV, sV{i}, 'parent', ud.scanaxis,...
                    'color', getColor(scanNum+nDoses, colors, 'loop'),...
                    'Tag', 'CERR_DOSEPROFILE_PLOT','linewidth',2);
            end
            axis(ud.scanaxis,'tight')
            ud.hDoseProfPlotV = hDoseProfPlotV;

            if exist('dV','var') && length(dV)>1 && ~isempty(ud.baseDose)
                delete(ud.htextInit)
                ud.htextInit = [];
                set(ud.diffaxis,'XTickLabelMode','auto','YTickLabelMode','auto')
                set(ud.diffaxis,'nextplot','add')
                otherDoses = drawDoses;
                baseDoseInd = find(otherDoses==baseDose);
                otherDosesInd = 1:length(otherDoses);
                otherDosesInd(baseDoseInd) = [];
                otherDoses(baseDoseInd) = [];

                for i=1:length(otherDosesInd)
                    ddV{i} = dV{otherDosesInd(i)}-dV{baseDoseInd};
                    if max(ddV{i})== 0
                        nddV{i} = ddV{i};
                    else
                        nddV{i} = ddV{i}/max(dV{baseDoseInd});
                    end
                end

                for i=1:length(ddV)
                    ddoseNum = otherDoses(i);
                    if strcmpi(stateS.doseDiffScale, 'Abs')
                        hDoseProfPlotV(end+1) = plot(distV, ddV{i}, 'parent', ud.diffaxis, 'color', ...
                            getColor(ddoseNum, colors, 'loop'),...
                            'Tag', 'CERR_DOSEPROFILE_PLOT','linewidth',2);
                    elseif strcmpi(stateS.doseDiffScale, 'Rel')
                        hDoseProfPlotV(end+1) = plot(distV, nddV{i}, 'parent', ud.diffaxis, 'color', ...
                            getColor(ddoseNum, colors, 'loop'),...
                            'Tag', 'CERR_DOSEPROFILE_PLOT','linewidth',2);
                    end

                end
                axis(ud.diffaxis,'tight')
                grid(ud.diffaxis,'on')
                ax = axis(ud.diffaxis);
                %text1 = get(doseUI(baseDose),'String');
                %dispText = [text1,' - other doses'];
                dispText = ['Dose - Ref'];
                htext = text(ax(2)-(ax(2)-ax(1))*0.03,ax(4)-(ax(4)-ax(3))*0.1,dispText,'parent', ud.diffaxis,'HorizontalAlignment','right','FontWeight','bold');
                ud.htext = htext;
                set(hFig, 'userdata',ud);
                xlabel(ud.diffaxis,'\bfDistance along line','fontSize',10);
                ylabel(ud.diffaxis,'\bfDifference','fontSize',10);
            else
                grid(ud.diffaxis,'off')
                set(ud.diffaxis,'xLim',[0 1]);
                set(ud.diffaxis,'yLim',[0 1]);
                set(ud.diffaxis,'XTickLabelMode','manual','YTickLabelMode','manual','XTickLabel',[],'YTickLabel',[]);
                delete(ud.htextInit)
                htextInit = text(0.1,0.5,'Select More than one dose to display difference','parent', ud.diffaxis,'HorizontalAlignment','left','FontWeight','bold');
                ud.htextInit = htextInit;
                ud.htext =[];
                %ud.hDoseProfPlotV = [];
                set(hFig, 'userdata',ud);
            end

        %end

    case 'CLOSE'
        if stateS.doseProfileState
            sliceCallBack('TOGGLEDOSEPROFILE');
        end
end

% AI: Added
function objCheck(hObj,hEvent)
%Dose/scan has been checked/unchecked
if ~isempty(hEvent.Indices)
hFig = get(hObj,'Parent');
if isempty(hFig)
    return
end
ud = get(hFig,'UserData');
delete(ud.htext)

if strcmp(get(hObj,'Tag'), 'doseUI')
    objUI  = ud.doseUI;
else
    objUI  = ud.scanUI;
end

selectedRowsV = hEvent.Indices(:,1);
objDataC = get(objUI,'Data');
drawObjV = get(objUI,'Data');
drawObjV = [drawObjV{:,1}];
drawObjV(selectedRowsV) = ~drawObjV(selectedRowsV);
objDataC(:,1) = num2cell(drawObjV);
set(objUI,'Data',objDataC);
%objUI.Data = objDataC;

if strcmp(get(hObj,'Tag'), 'doseUI')
    ud.doseUI = objUI;
else
    ud.scanUI = objUI;
end
ud.htext = [];
set(hFig, 'userdata',ud);
doseProfileFigure('refresh');
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [IX,IX1,IX2] = FindLinearIdx(x,y,lon,lat)
% given points in vectors x,y (np x 1) find their linear indices IX (np x 1)
% from a matrix of x-locations X and y-locations Y of both size nx columns and ny rows.
% lon and lat must be matrices created by ndgrid.
% kjr, 20171210 chl, und.
% wjp, 20201120 making sure dx and dy can be different
%      20200325 implementing method for irregular grid, cleaning up
ny = size(lon,1);
nx = size(lon,2);
np = numel(x);

if nx == 1 && ny == 1
   IX = 1; IX1 = 1; IX2 = 1;
   return
end

% make sure entry points are column vectors
x = x(:);
y = y(:);

% Get the grid spacing in x and y directions
% (trying both directions so could be meshgrid or ndgrid format)
dx  = diff(lon(:,1));
dy  = diff(lat(1,:));
dx_cutoff = 0.1/111e3; % approx 10 cm
if (max(dx) - min(dx)) > dx_cutoff || (max(dy) - min(dy)) > dx_cutoff
    % % IRREGULAR GRID (SLOWER)

    % convert ndgrid to vector
    lonV = reshape(lon(:,1),1,[]);
    % repeat the vector along the length of x
    LON = repmat(lonV,np,1);
    % find the closest lon to each point of x
    [~,IX1] = min(abs(x - LON),[],2);

    % convert ndgrid to vector
    latV = lat(1,:);
    % repeat the vector along the length of y
    LAT = repmat(latV,np,1);
    % find the closest lon to each point of y
    [~,IX2] = min(abs(y - LAT),[],2);

else
    % % REGULAR GRID (FASTER)
    dx = dx(1); dy = dy(1);

    IX1 = (x-lon(1,1))/dx + 1;
    IX2 = (y-lat(1,1))/dy + 1;

    IX1 = round(IX1);
    IX2 = round(IX2);

    IX1 = max([IX1,ones(np,1)],[],2);
    IX1 = min([IX1,ny*ones(np,1)],[],2);

    IX2 = max([IX2,ones(np,1)],[],2);
    IX2 = min([IX2,nx*ones(np,1)],[],2);

end

IX = sub2ind([ny,nx],IX1,IX2);

end

author: CHLNDDEV


text:
function [paramtransformsettings, names] = nddisimKernExtractParamTransformSettings(kern)

% NDDISIMKERNEXTRACTPARAMTRANSFORMSETTINGS Extract parameter transform settings from the NDDISIM kernel structure.
% FORMAT
% DESC Extract parameters from the single input motif kernel structure
% into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC Extract parameters and their names from the single input
% motif kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel matrix, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings containing parameter names.
%
% SEEALSO disimKernParamInit, disimKernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Antti Honkela, 2007-2009
%
% COPYRIGHT : Jaakko Peltonen, 2011
%
% KERN

paramtransformsettings = {kern.transforms(1).transformsettings,kern.transforms(2).transformsettings,kern.transforms(3).transformsettings,kern.transforms(4).transformsettings,kern.transforms(5).transformsettings};

if nargout > 1
  names = {'inverse width', 'di_variance', 'decay', 'variance', 'delay'};
end

%fprintf(1, 'disimKern parameters physical values:\n');
%params
%names

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [L,S] = gosus(M,opts)
  regularizer = 0.1;
  param.superpixel.slicParam = [10, regularizer; 20, regularizer; 40,regularizer; 80, regularizer;];
  %  param.superpixel.slicParam = [5,regularizer;  10, regularizer; 20, regularizer;  40,regularizer ; ];
  param.maxIter = 200;
  param.tol = 1e-4;
  param.lambda = 0.33;    %   bestAccu = 0.9762, lambda = 0.33; threshold = 8e-6;
  param.rank = 2; % subspace dimension
  param.eta = 1e-3;   % stepsize for subspace updating
  %param.sampleSize = 20; % number of sample frames to learn the background
  %param.trainSize = 24; % sample from the fist 240 frames to learn the background
  %param.startIndex = 1;   %  which frame to start
  %param.randomStart = true;
  param.admm.I = speye(size(M,1)*3);

  %U = orth(randn(size(M,1)*3, param.rank));
  [U, ~] = svds([M(:,1:25);M(:,1:25);M(:,1:25)], param.rank);

  for i = 1:size(M,2)
    disp(i);
    I = M(:,i);
    im = im2uint8(reshape(I,[opts.rows opts.cols]));
    im = repmat(im,[1 1 3]);

    param.admm.G = getGroupSuperColor27(im, param);
    param.admm.Z = sparse(size(M,1)*3, size(param.admm.G, 2), 0);
    param.admm.Y = param.admm.Z;

    v = double(im(:))/255.0;

    [x, w] = solveWXADMM(U, v, param);

    pL = U*w; % low-rank
    pS = pL-v; % sparse
    %E = v;

    imgL = reshape(pL, size(im));
    imgS = reshape(pS, size(im));

    vL = imgL(:,:,1);
    vS = imgS(:,:,1);

    vL = vL(:);
    vS = vS(:);

    L(:,i) = vL;
    S(:,i) = vS;

    %L = U*w; % low-rank
    %E = v;

    %fg = L-v;
    %fg(abs(x)<=8e-10) = 0;
    %fg(abs(x)>8e-10) = 1;

    %subplot(1,3,1), imshow(im,[]);
    %subplot(1,3,2), imshow(reshape(L, size(im)),[]);
    %subplot(1,3,3), imshow(reshape(fg, size(im)),[]);

    residual = param.lambda*(U*w + x - v);
    U = updateSubspace(U, residual, w, param);

    %pause(.01);
  end
end

author: andrewssobral

repo: lrslibrary