text:
function [C,T]=hungarian(A)
% HUNGARIAN - Solve the assignment problem using the Hungarian method.
%
% Usage: >> [C,T]=hungarian(A)
%
% Inputs:
%   A - a square (correlation/distance/cost) matrix.
%   C - the optimal assignment (closest row/col pairs)
%   T - the total (minimized) cost of the optimal assignment.
%
% Author: Niclas Borlin, 1996

% Adapted from the FORTRAN IV code in Carpaneto and Toth, "Algorithm 548:
% Solution of the assignment problem [H]", ACM Transactions on
% Mathematical Software, 6(1):104-111, 1980.
%
% v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.
%                 Department of Computing Science, Ume University,
%                 Sweden.
%
% NOTE: A substantial effort was put into this code. If you use it for a
% publication or otherwise, please include an acknowledgement and notify
% me by email. /Niclas

[m,n]=size(A);

if (m~=n)
    error('HUNGARIAN: Cost matrix must be square!');
end

% Save original cost matrix.
orig=A;

% Reduce matrix.
A=hminired(A);

% Do an initial assignment.
[A,C,U]=hminiass(A);

% Repeat while we have unassigned rows.
while (U(n+1))
    % Start with no path, no unchecked zeros, and no unexplored rows.
    LR=zeros(1,n);
    LC=zeros(1,n);
    CH=zeros(1,n);
    RH=[zeros(1,n) -1];

    % No labelled columns.
    SLC=[];

    % Start path in first unassigned row.
    r=U(n+1);
    % Mark row with end-of-path label.
    LR(r)=-1;
    % Insert row first in labelled row set.
    SLR=r;

    % Repeat until we manage to find an assignable zero.
    while (1)
        % If there are free zeros in row r
        if (A(r,n+1)~=0)
            % ...get column of first free zero.
            l=-A(r,n+1);

            % If there are more free zeros in row r and row r in not
            % yet marked as unexplored..
            if (A(r,l)~=0 && RH(r)==0)
                % Insert row r first in unexplored list.
                RH(r)=RH(n+1);
                RH(n+1)=r;

                % Mark in which column the next unexplored zero in this row
                % is.
                CH(r)=-A(r,l);
            end
        else
            % If all rows are explored..
            if (RH(n+1)<=0)
                % Reduce matrix.
                [A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR);
            end

            % Re-start with first unexplored row.
            r=RH(n+1);
            % Get column of next free zero in row r.
            l=CH(r);
            % Advance "column of next free zero".
            CH(r)=-A(r,l);
            % If this zero is last in the list..
            if (A(r,l)==0)
                % ...remove row r from unexplored list.
                RH(n+1)=RH(r);
                RH(r)=0;
            end
        end

        % While the column l is labelled, i.e. in path.
        while (LC(l)~=0)
            % If row r is explored..
            if (RH(r)==0)
                % If all rows are explored..
                if (RH(n+1)<=0)
                    % Reduce cost matrix.
                    [A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR);
                end

                % Re-start with first unexplored row.
                r=RH(n+1);
            end

            % Get column of next free zero in row r.
            l=CH(r);

            % Advance "column of next free zero".
            CH(r)=-A(r,l);

            % If this zero is last in list..
            if(A(r,l)==0)
                % ...remove row r from unexplored list.
                RH(n+1)=RH(r);
                RH(r)=0;
            end
        end

        % If the column found is unassigned..
        if (C(l)==0)
            % Flip all zeros along the path in LR,LC.
            [A,C,U]=hmflip(A,C,LC,LR,U,l,r);
            % ...and exit to continue with next unassigned row.
            break;
        else
            % ...else add zero to path.

            % Label column l with row r.
            LC(l)=r;

            % Add l to the set of labelled columns.
            SLC=[SLC l];

            % Continue with the row assigned to column l.
            r=C(l);

            % Label row r with column l.
            LR(r)=l;

            % Add r to the set of labelled rows.
            SLR=[SLR r];
        end
    end
end

% Calculate the total cost.
T=sum(orig(logical(sparse(C,1:size(orig,2),1))));

function A=hminired(A)
%HMINIRED Initial reduction of cost matrix for the Hungarian method.
%
%B=assredin(A)
%A - the unreduced cost matris.
%B - the reduced cost matrix with linked zeros in each row.

% v1.0  96-06-13. Niclas Borlin, niclas@cs.umu.se.

[m,n]=size(A);

% Subtract column-minimum values from each column.
colMin=min(A);
A=A-colMin(ones(n,1),:);

% Subtract row-minimum values from each row.
rowMin=min(A')';
A=A-rowMin(:,ones(1,n));

% Get positions of all zeros.
[i,j]=find(A==0);

% Extend A to give room for row zero list header column.
A(1,n+1)=0;
for k=1:n
    % Get all column in this row.
    cols=j(k==i)';
    % Insert pointers in matrix.
    A(k,[n+1 cols])=[-cols 0];
end

function [A,C,U]=hminiass(A)
%HMINIASS Initial assignment of the Hungarian method.
%
%[B,C,U]=hminiass(A)
%A - the reduced cost matrix.
%B - the reduced cost matrix, with assigned zeros removed from lists.
%C - a vector. C(J)=I means row I is assigned to column J,
%              i.e. there is an assigned zero in position I,J.
%U - a vector with a linked list of unassigned rows.

% v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.

[n,np1]=size(A);

% Initialize return vectors.
C=zeros(1,n);
U=zeros(1,n+1);

% Initialize last/next zero "pointers".
LZ=zeros(1,n);
NZ=zeros(1,n);

for i=1:n
    % Set j to first unassigned zero in row i.
	lj=n+1;
	j=-A(i,lj);

    % Repeat until we have no more zeros (j==0) or we find a zero
	% in an unassigned column (c(j)==0).

	while (C(j)~=0)
		% Advance lj and j in zero list.
		lj=j;
		j=-A(i,lj);

		% Stop if we hit end of list.
		if (j==0)
			break;
		end
	end

	if (j~=0)
		% We found a zero in an unassigned column.

		% Assign row i to column j.
		C(j)=i;

		% Remove A(i,j) from unassigned zero list.
		A(i,lj)=A(i,j);

		% Update next/last unassigned zero pointers.
		NZ(i)=-A(i,j);
		LZ(i)=lj;

		% Indicate A(i,j) is an assigned zero.
		A(i,j)=0;
	else
		% We found no zero in an unassigned column.

		% Check all zeros in this row.

		lj=n+1;
		j=-A(i,lj);

		% Check all zeros in this row for a suitable zero in another row.
		while (j~=0)
			% Check the in the row assigned to this column.
			r=C(j);

			% Pick up last/next pointers.
			lm=LZ(r);
			m=NZ(r);

			% Check all unchecked zeros in free list of this row.
			while (m~=0)
				% Stop if we find an unassigned column.
				if (C(m)==0)
					break;
				end

				% Advance one step in list.
				lm=m;
				m=-A(r,lm);
			end

			if (m==0)
				% We failed on row r. Continue with next zero on row i.
				lj=j;
				j=-A(i,lj);
			else
				% We found a zero in an unassigned column.

				% Replace zero at (r,m) in unassigned list with zero at (r,j)
				A(r,lm)=-j;
				A(r,j)=A(r,m);

				% Update last/next pointers in row r.
				NZ(r)=-A(r,m);
				LZ(r)=j;

				% Mark A(r,m) as an assigned zero in the matrix . . .
				A(r,m)=0;

				% ...and in the assignment vector.
				C(m)=r;

				% Remove A(i,j) from unassigned list.
				A(i,lj)=A(i,j);

				% Update last/next pointers in row r.
				NZ(i)=-A(i,j);
				LZ(i)=lj;

				% Mark A(r,m) as an assigned zero in the matrix . . .
				A(i,j)=0;

				% ...and in the assignment vector.
				C(j)=i;

				% Stop search.
				break;
			end
		end
	end
end

% Create vector with list of unassigned rows.

% Mark all rows have assignment.
r=zeros(1,n);
rows=C(C~=0);
r(rows)=rows;
empty=find(r==0);

% Create vector with linked list of unassigned rows.
U=zeros(1,n+1);
U([n+1 empty])=[empty 0];

function [A,C,U]=hmflip(A,C,LC,LR,U,l,r)
%HMFLIP Flip assignment state of all zeros along a path.
%
%[A,C,U]=hmflip(A,C,LC,LR,U,l,r)
%Input:
%A   - the cost matrix.
%C   - the assignment vector.
%LC  - the column label vector.
%LR  - the row label vector.
%U   - the
%r,l - position of last zero in path.
%Output:
%A   - updated cost matrix.
%C   - updated assignment vector.
%U   - updated unassigned row list vector.

% v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.

n=size(A,1);

while (1)
    % Move assignment in column l to row r.
    C(l)=r;

    % Find zero to be removed from zero list..

    % Find zero before this.
    m=find(A(r,:)==-l);

    % Link past this zero.
    A(r,m)=A(r,l);

    A(r,l)=0;

    % If this was the first zero of the path..
    if (LR(r)<0)
        % remove row from unassigned row list and return.
        U(n+1)=U(r);
        U(r)=0;
        return;
    else

        % Move back in this row along the path and get column of next zero.
        l=LR(r);

        % Insert zero at (r,l) first in zero list.
        A(r,l)=A(r,n+1);
        A(r,n+1)=-l;

        % Continue back along the column to get row of next zero in path.
        r=LC(l);
    end
end

function [A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR)
%HMREDUCE Reduce parts of cost matrix in the Hungerian method.
%
%[A,CH,RH]=hmreduce(A,CH,RH,LC,LR,SLC,SLR)
%Input:
%A   - Cost matrix.
%CH  - vector of column of 'next zeros' in each row.
%RH  - vector with list of unexplored rows.
%LC  - column labels.
%RC  - row labels.
%SLC - set of column labels.
%SLR - set of row labels.
%
%Output:
%A   - Reduced cost matrix.
%CH  - Updated vector of 'next zeros' in each row.
%RH  - Updated vector of unexplored rows.

% v1.0  96-06-14. Niclas Borlin, niclas@cs.umu.se.

n=size(A,1);

% Find which rows are covered, i.e. unlabelled.
coveredRows=LR==0;

% Find which columns are covered, i.e. labelled.
coveredCols=LC~=0;

r=find(~coveredRows);
c=find(~coveredCols);

% Get minimum of uncovered elements.
m=min(min(A(r,c)));

% Subtract minimum from all uncovered elements.
A(r,c)=A(r,c)-m;

% Check all uncovered columns..
for j=c
    % ...and uncovered rows in path order..
    for i=SLR
        % If this is a (new) zero..
        if (A(i,j)==0)
            % If the row is not in unexplored list..
            if (RH(i)==0)
                % ...insert it first in unexplored list.
                RH(i)=RH(n+1);
                RH(n+1)=i;
                % Mark this zero as "next free" in this row.
                CH(i)=j;
            end
            % Find last unassigned zero on row I.
            row=A(i,:);
            colsInList=-row(row<0);
            if (length(colsInList)==0)
                % No zeros in the list.
                l=n+1;
            else
                l=colsInList(row(colsInList)==0);
            end
            % Append this zero to end of list.
            A(i,l)=-j;
        end
    end
end

% Add minimum to all doubly covered elements.
r=find(coveredRows);
c=find(coveredCols);

% Take care of the zeros we will remove.
[i,j]=find(A(r,c)<=0);

i=r(i);
j=c(j);

for k=1:length(i)
    % Find zero before this in this row.
    lj=find(A(i(k),:)==-j(k));
    % Link past it.
    A(i(k),lj)=A(i(k),j(k));
    % Mark it as assigned.
    A(i(k),j(k))=0;
end

A(r,c)=A(r,c)+m;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [OUT, W] = CrossWeightAverage(DATA,WEIGHT)
% =======================================================================
% Compute the weighted average of DATA, computing the weights form the
% matrix WEIGHT. Note that if there are NaNs in DATA the weights are
% rescaled. WEIGHTS cannot contain NaN
% =======================================================================
% [average, weight] = CrossWeightAverage(DATA,WEIGHT)
% -----------------------------------------------------------------------
% INPUT
%	- DATA: panel of time series DATA(T,N) of T observation x N series
%       (NaN are accepted and weights rescaled)
%   - WEIGHT: time series of WEIGHTS (TxN). You can also provide a (1xN)
%       matrix in this case the code would automatically compute fixed
%       weights
% -----------------------------------------------------------------------
% OUTPUT
%	- average : weighted average vector (T x 1) [double]
%	- weight  : matrix of weights (T x N) [double]
% -----------------------------------------------------------------------
% EXAMPLE 1: foreign variables as in GVAR (example with fixed weights)
%   x = [1 NaN  3 ;
%         2   3  4 ;
%         3   4  5 ;
%         4   5  6];
%    w = [0  .4 .6];
%    [average, weight] = CrossWeightAverage(x,w)
% -----------------------------------------------------------------------
% EXAMPLE 2: weighted average (with time-varying weights)
%    x = [ 1 NaN   3 ;
%          2   3   4 ;
%          3   4   5 ;
%          4   5   6];
%    w = [.4  .2  .4 ;
%         .3  .4  .3 ;
%         .4  .3  .3 ;
%         .5  .3  .2];
%    [average, weight] = CrossWeightAverage(x,w)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

%% Preliminaries
% Check inputs
[r1, c1] = size(DATA);
[r2, c2] = size(WEIGHT);
if r2==1 && c2==c1
    % Create a matrix of fixed weights
    aux = [];
    for ii=1:r1
        aux = [aux; WEIGHT];
    end
    WEIGHT = aux;
end

% Re-Check inputs
[r1, c1] = size(DATA);
[r2, c2] = size(WEIGHT);
if r1~=r2 || c1~=c2
    disp(' ')
    disp('The matrices of data and weights are not conformable.')
    return
end

%% Cross section weihgted average
[nobs, nvars] = size(DATA);
nans = isnan(DATA);
WEIGHT(nans) = NaN;

OUT(1:nobs,1) = NaN;
W(1:nobs,1:nvars) = NaN;

for ii=1:nobs
    aux_weight = WEIGHT(ii,:)./nansum(WEIGHT(ii,:));
    W(ii,:) = aux_weight;
    aux_weight(nans(ii,:)) = 0;
    DATA(ii,nans(ii,:)) = 0;
    if sum(nans(ii,:))==nvars
        OUT(ii,1) = NaN;
    else
        OUT(ii,1) = DATA(ii,:) * aux_weight';
    end
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function y=rdct(x,n,a,b)
%RDCT     Discrete cosine transform of real data Y=(X,N,A,B)
% Data is truncated/padded to length N.
%
% This routine is equivalent to multiplying by the matrix
%
%   rdct(eye(n)) = diag([sqrt(2)*B/A repmat(2/A,1,n-1)]) * cos((0:n-1)'*(0.5:n)*pi/n)
%
% Default values of the scaling factors are A=sqrt(2N) and B=1 which
% results in an orthogonal matrix. Other common values are A=1 or N and/or B=1 or sqrt(2).
% If b~=1 then the columns are no longer orthogonal.
%
% see IRDCT for the inverse transform

% BUG: in line 51 we should do chopping after transform and not before

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: rdct.m,v 1.6 2007/05/04 07:01:39 dmb Exp $
%
%   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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fl=size(x,1)==1;
if fl x=x(:); end
[m,k]=size(x);
if nargin<2 n=m;
end
if nargin<4 b=1;
    if nargin<3 a=sqrt(2*n);
    end
    end
if n>m x=[x; zeros(n-m,k)];
elseif n<m x(n+1:m,:)=[];
end

x=[x(1:2:n,:); x(2*fix(n/2):-2:2,:)];
z=[sqrt(2) 2*exp((-0.5i*pi/n)*(1:n-1))].';
y=real(fft(x).*z(:,ones(1,k)))/a;
y(1,:)=y(1,:)*b;
if fl y=y.'; end

author: decouples


text:
% Script: call_llforecast.m
% -----------------------
% Script to input parameters for an Omori law fit
% Calculation is done in function plot_fitloglike_a2
% The function works on the catalog newt2!!!
% J.Woessner
% last update: 10.07.03

report_this_filefun(mfilename('fullpath'));
% Get input parameters
prompt  = {,'Enter length of learning period (days)','Enter forecast period (days):','Enter number of bootstraps:'};
title   = 'Parameters ';
lines= 1;
def     = {'50','50','100'};
answer  = inputdlg(prompt,title,lines,def);
time = str2double(answer{1});
timef = str2double(answer{2});
bootloops = str2double(answer{3});

% maepi is the mainshock of sequence
plot_fitloglike_a2(newt2,time,timef,bootloops,maepi);

author: CelsoReyes


text:
function A = set(A, varargin)
%
% Set psfMatrix properties and return the updated object.
%
%     A = set(A, ...)
%
%  This funtion accepts an psfMatrix object, A, and a variable list of
%  property name/value pairs and returns the modified object.
%  Valid properties are:
%    'psf', 'matdata', 'type', 'boundary', 'transpose'
%

%  J. Nagy  5/2/01

property_argin = varargin;

while length( property_argin ) >= 2
  prop = property_argin{1};
  val = property_argin{2};
  property_argin = property_argin(3:end);

  switch prop
  case 'psf'
    A.psf = val;
  case 'matdata'
    A.matdata = val;
  case 'type'
    A.type = val;
  case 'boundary'
    A.boundary = val;
  case 'transpose'
    A.transpose = val;
  otherwise
    error('Valid psfMatrix properties: psf, matdata, type, boundary, transpose')
  end
end

repo: IRtools


text:
function [ sicd_meta_current ] = sicd_update_meta_0_5( sicd_meta_0_5 )
%SICD_UPDATE_META_0_5 Update a SICD metadata structure from version 0.5 to
%current version (whatever that may be)
%
% Author: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

sicd_meta_1_0 = sicd_meta_0_5;

% Add RadarCollection.TxPolarization, now required, but optional prior to version 1.0
if isfield(sicd_meta_0_5, 'RadarCollection') && ...
        ~isfield(sicd_meta_0_5.RadarCollection, 'TxPolarization') && ...
        isfield(sicd_meta_0_5.RadarCollection, 'RcvChannels') && ...
        isfield(sicd_meta_0_5.RadarCollection.RcvChannels, 'ChanParameters')
    if numel(sicd_meta_0_5.RadarCollection.RcvChannels.ChanParameters)==1
        sicd_meta_1_0.RadarCollection.TxPolarization = sicd_meta_0_5.RadarCollection.RcvChannels.ChanParameters.TxRcvPolarization(1);
    else
        sicd_meta_1_0.RadarCollection.TxPolarization = 'SEQUENCE';
        tx_pols = {sicd_meta_0_5.RadarCollection.RcvChannels.ChanParameters.TxRcvPolarization};
        tx_pols = unique(cellfun(@(x) x(1), tx_pols));
        for i = 1:numel(tx_pols)
            sicd_meta_1_0.RadarCollection.TxSequence.TxStep(i).TxPolarization = tx_pols(i);
        end
        % Note: If there are multiple waveforms and multiple polarizations,
        % there is no deconfliction done here.
    end
end

% RadarCollection.Area.Corner was optional in version 0.5, but required in
% version 1.0.  Fortunately, Corner is easily derived from Plane.
if isfield(sicd_meta_0_5, 'RadarCollection') && ...
        isfield(sicd_meta_0_5.RadarCollection, 'Area') && ...
        ~isfield(sicd_meta_0_5.RadarCollection.Area, 'Corner') && ...
        isfield(sicd_meta_0_5.RadarCollection.Area, 'Plane')
    try % If Plane substructure is misformed, this may fail
        plane = sicd_meta_0_5.RadarCollection.Area.Plane; % For concise notation
        ref_pt = [plane.RefPt.ECF.X plane.RefPt.ECF.Y plane.RefPt.ECF.Z];
        x_uvect = [plane.XDir.UVectECF.X plane.XDir.UVectECF.Y plane.XDir.UVectECF.Z];
        y_uvect = [plane.YDir.UVectECF.X plane.YDir.UVectECF.Y plane.YDir.UVectECF.Z];
        x_offsets = [plane.XDir.FirstLine plane.XDir.FirstLine ...
            plane.XDir.NumLines plane.XDir.NumLines];
        y_offsets = [plane.YDir.FirstSample plane.YDir.NumSamples ...
            plane.YDir.NumSamples plane.YDir.FirstSample];
        for i = 1:4
            acp = ref_pt + x_uvect * plane.XDir.LineSpacing * double(x_offsets(i) - plane.RefPt.Line) + ...
                y_uvect * plane.YDir.SampleSpacing * double(y_offsets(i) - plane.RefPt.Sample);
            sicd_meta_1_0.RadarCollection.Area.Corner.ACP(i) = ...
                cell2struct(num2cell(ecf_to_geodetic(acp)),{'Lat','Lon','HAE'});
        end
    end
end

% PolarizationHVAnglePoly no longer a valid field in version 1.0.
if isfield(sicd_meta_0_5, 'RadarCollection') && ...
        isfield(sicd_meta_0_5.RadarCollection, 'PolarizationHVAnglePoly')
    sicd_meta_1_0.RadarCollection = rmfield(sicd_meta_1_0.RadarCollection, 'PolarizationHVAnglePoly');
end

% Antenna.Tx/Rcv/TwoWay.HPBW no longer a valid field in version 1.0.
if isfield(sicd_meta_0_5, 'Antenna')
    if isfield(sicd_meta_0_5.Antenna, 'Tx') && ...
            isfield(sicd_meta_0_5.Antenna.Tx, 'HPBW')
        sicd_meta_1_0.Antenna.Tx = rmfield(sicd_meta_0_5.Antenna.Tx, 'HPBW');
    end
    if isfield(sicd_meta_0_5.Antenna, 'Rcv') && ...
            isfield(sicd_meta_0_5.Antenna.Rcv, 'HPBW')
        sicd_meta_1_0.Antenna.Rcv = rmfield(sicd_meta_0_5.Antenna.Rcv, 'HPBW');
    end
    if isfield(sicd_meta_0_5.Antenna, 'TwoWay') && ...
            isfield(sicd_meta_0_5.Antenna.TwoWay, 'HPBW')
        sicd_meta_1_0.Antenna.TwoWay = rmfield(sicd_meta_0_5.Antenna.TwoWay, 'HPBW');
    end
end

% NoiseLevel got its own substructure between SICD 0.5 and SICD 1.0
if isfield(sicd_meta_0_5, 'Radiometric') && ...
        isfield(sicd_meta_0_5.Radiometric, 'NoisePoly')
    sicd_meta_1_0.Radiometric.NoiseLevel.NoisePoly = ...
            sicd_meta_0_5.Radiometric.NoisePoly;
    sicd_meta_1_0.Radiometric = rmfield(sicd_meta_1_0.Radiometric, 'NoisePoly');
    if isfield(sicd_meta_0_5.Radiometric, 'NoiseLevelType')
        sicd_meta_1_0.Radiometric.NoiseLevel.NoiseLevelType = ...
            sicd_meta_0_5.Radiometric.NoiseLevelType;
    else
        % Even if NoiseLevelType wasn't given, we know that relative noise
        % levels should be 1 at SCP.
        if abs(sicd_meta_0_5.Radiometric.NoisePoly(1)-1)<eps
            sicd_meta_1_0.Radiometric.NoiseLevel.NoiseLevelType = 'RELATIVE';
        else
            sicd_meta_1_0.Radiometric.NoiseLevel.NoiseLevelType = 'ABSOLUTE';
        end
    end
end
if isfield(sicd_meta_0_5, 'Radiometric') && ...
        isfield(sicd_meta_0_5.Radiometric, 'NoiseLevelType')
    sicd_meta_1_0.Radiometric = rmfield(sicd_meta_1_0.Radiometric, 'NoiseLevelType');
end

% MatchInfo
if isfield(sicd_meta_0_5, 'MatchInfo')
    sicd_meta_1_0.MatchInfo = struct(); % Clear this out so we can reconstruct it
    if isfield(sicd_meta_0_5.MatchInfo.Collect, 'MatchType') % MatchType was optional field in 0.5
        types = unique({sicd_meta_0_5.MatchInfo.Collect.MatchType});
    else
        types = {''}; % TypeID (equivalent of MatchType) required in 1.0
    end
    sicd_meta_1_0.MatchInfo.NumMatchTypes = numel(types);
    sicd_meta_1_0.MatchInfo.MatchType = struct();
    for i=1:numel(sicd_meta_1_0.MatchInfo.NumMatchTypes)
        if isfield(sicd_meta_0_5.MatchInfo.Collect, 'MatchType')
            ind = find(strcmp({sicd_meta_0_5.MatchInfo.Collect.MatchType},types{i}));
        else
            ind = 1:numel(sicd_meta_0_5.MatchInfo.Collect);
        end
        sicd_meta_1_0.MatchInfo.MatchType(i).TypeID = strtrim(types{i});
        sicd_meta_1_0.MatchInfo.MatchType(i).NumMatchCollections = ...
            numel(ind) - 1; % 0.5 included current instance as one of the collections
        for j = 1:numel(ind)
            if isfield(sicd_meta_0_5.MatchInfo.Collect(ind(j)), 'Parameter')
                if iscell(sicd_meta_0_5.MatchInfo.Collect(ind(j)).Parameter) % Multiple parameters
                    current_index = find(cellfun(@(x) strcmpi(strtrim(x.name),'CURRENT_INSTANCE'), ...
                        sicd_meta_0_5.MatchInfo.Collect(ind(j)).Parameter),1);
                    if ~isempty(current_index)
                        current_index = str2double(sicd_meta_0_5.MatchInfo.Collect(ind(j)).Parameter{current_index}.value);
                    end
                elseif strcmpi(strtrim(sicd_meta_0_5.MatchInfo.Collect(ind(j)).Parameter.name), 'CURRENT_INSTANCE') % Single parameter
                    current_index = str2double(sicd_meta_0_5.MatchInfo.Collect(ind(j)).Parameter.value);
                else
                    current_index = [];
                end
            end
            if ~isempty(current_index)
                sicd_meta_1_0.MatchInfo.MatchType(i).CurrentIndex = current_index;
            else
                sicd_meta_1_0.MatchInfo.MatchType(i).MatchCollection(j).CoreName = ...
                    sicd_meta_0_5.MatchInfo.Collect(ind(j)).CoreName;
                if isfield(sicd_meta_0_5.MatchInfo.Collect(ind(j)), 'Parameter')
                    sicd_meta_1_0.MatchInfo.MatchType(i).MatchCollection(j).Parameter = ...
                        sicd_meta_0_5.MatchInfo.Collect(ind(j)).Parameter;
                end
            end
        end
    end
end

% Add AzimAng and LayoverAng to SCPCOA
sicd_meta_1_0 = derived_sicd_fields(sicd_meta_1_0);

% We are now updated to version 1.0.  Now do rest of updates.
% sicd_meta_current = sicd_update_meta_1_0(sicd_meta_1_0); % If we needed
% to to more updates past version 1.0, they would go here
sicd_meta_current = sicd_meta_1_0; % No further updates

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function plotGamma(as,bs)
% as = [1 1.5 2];
% b = 1;
% bs = b*ones(1,length(as));
% figure;
[styles, colors, symbols] = plotColors;
legendStr = cell(1, length(as));
%styles = {'k:',  'b--', 'r-'};
% textstr={'Iter.=20', 'Iter.=22', 'Iter.=25', 'Iter.=33'};
for i=1:length(as)
    a = as(i); b = bs(i);
    xs = linspace(max(0, a/b - 5* sqrt(a/(b^2))),a/b + 5*sqrt(a/(b^2)), 50);
    model.a = a;
    model.b = b;
    ps = exp(gammaLogprob(model, xs));
    style = [styles{i}, colors(i), symbols(i)];
    plot(xs , ps, style, 'linewidth', 1.5, 'markersize', 1.2);
%     plot(xs , ps, styles{i}, 'color', colors(i), 'linewidth', 1.5);
%     text(a/b,exp(gammaLogprob(model, a/b)), sprintf('a=%2.1f,b=%2.1f', a, b));
    hold on
    legendStr{i} = sprintf('a=%2.1f,\nb=%2.1f', a, b);
    %     axis tight;
end
% legend(legendStr, 'fontsize', 10);
% title('Gamma distributions');
hold off;
end

function [styles, colors, symbols, str] =  plotColors()
% Use plot(x,y,str{i}) to print in i'th style
% We use colors and linestyles, not markers

% This file is from pmtk3.googlecode.com

colors =  ['b' 'r' 'k' 'g' 'c' 'y' 'm' ...
    'r' 'b' 'k' 'g' 'c' 'y' 'm'];
symbols = ['o' 'x' '*' '>' '<' '^' 'v' ...
    '+' 'p' 'h' 's' 'd' 'o' 'x'];
styles = {'-', ':', '-.', '--', '-', ':', '-.', '--', ...
    '-', ':', '-.', '--', '-', ':', '-.', '--'};

for i=1:length(colors)
    %str{i} = sprintf('-%s%s', colors(i), symbols(i));
    str{i} = sprintf('%s%s', colors(i), styles{i});
end

end

function logp = gammaLogprob(arg1, arg2, arg3)
% logp(i) = log p(X(i) | a, b)
% logp = gammaLogprob(model, X); OR logp = gammaLogprob(a, b, X);
% a is the shape,
% b is the rate, i.e. 1/scale

% This file is from pmtk3.googlecode.com

if isstruct(arg1)
    model = arg1;
    X = arg2;
    a = model.a;
    b = model.b;
else
    a = arg1;
    b = arg2;
    X = arg3;
end

logZ = gammaln(a) - a.*log(b);
logp = (a-1).*log(X) - b.*X - logZ;
logp = logp(:);

end

author: andrewssobral


text:
function [varargout] = funname(varargin)

% WRITE_CTF_SHM writes metainformation and data as a packet to shared memory.
% This function can be used for real-time processing of data while it is
% being acquired.
%
% Use as
%   write_ctf_shm(msgType, msgId, sampleNumber, numSamples, numChannels, data);
%
% See also READ_CTF_SHM

% Copyright (C) 2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: write_ctf_shm.m 2885 2011-02-16 09:41:58Z roboos $

% remember the original working directory
pwdir = pwd;

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc  = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

try
  % try to compile the mex file on the fly
  warning('trying to compile MEX file from %s', mexsrc);
  cd(mexdir);
  mex(mexsrc);
  cd(pwdir);
  success = true;

catch
  % compilation failed
  disp(lasterr);
  error('could not locate MEX file for %s', mexname);
  cd(pwdir);
  success = false;
end

if success
  % execute the mex file that was juist created
  funname   = mfilename;
  funhandle = str2func(funname);
  [varargout{1:nargout}] = funhandle(varargin{:});
end

author: PatternRecognition

repo: OpenBMI


text:
% eegplot2event() - convert EEGPLOT rejections into events
%                   compatible with the eeg_eegrej function for rejecting
%                   continuous portions of datasets.
%
% Usage:
%   >> [events] = eegplot2event( eegplotrej, type, colorin, colorout );
%
% Inputs:
%   eegplotrej - EEGPLOT output (TMPREJ; see eegplot for more details)
%   type       - type of the event. Default -1.
%   colorin    - only extract rejection of specific colors (here a n x 3
%                array must be given). Default: extract all rejections.
%   colorout   - do not extract rejection of specified colors.
%
% Outputs:
%   events     - array of events compatible with the eeg_eegrej function
%                for rejecting continuous portions of datasets.
%
% Example:
% NEWEEG = eeg_eegrej(EEG,eegplot2event(TMPREJ, -1));
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: eegplot(), eeg_multieegplot(), eegplot2trial(), eeglab()

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% 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 should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function events = eegplot2event( TMPREJ, type, color, colorout )

if nargin < 1
   help eegplot2event;
   return;
end;

if nargin < 2
   type = -1;
end;

% only take into account specific colors
% --------------------------------------
if exist('color') == 1
	for index = 1:size(color,1)
   		tmpcol1 = TMPREJ(:,3) + 255*TMPREJ(:,4) + 255*255*TMPREJ(:,5);
   		tmpcol2 = color(index,1)+255*color(index,2)+255*255*color(index,3);
   		I = find( tmpcol1 == tmpcol2);
   		if isempty(I)
      		fprintf('Warning: color [%d %d %d] not found\n', ...
      			color(index,1), color(index,2), color(index,3));
   		end;
   		TMPREJ = TMPREJ(I,:);
   	end;
end;

% remove other specific colors
% ----------------------------
if exist('colorout') == 1
	for index = 1:size(colorout,1)
   		tmpcol1 = TMPREJ(:,3) + 255*TMPREJ(:,4) + 255*255*TMPREJ(:,5);
   		tmpcol2 = colorout(index,1)+255*colorout(index,2)+255*255*colorout(index,3);
   		I = find( tmpcol1 ~= tmpcol2);
   		TMPREJ = TMPREJ(I,:);
   	end;
end;

events = [];
if ~isempty(TMPREJ)
   events = TMPREJ(:,1:5);
   events = [ type*ones(size(events,1), 1) ones(size(events,1), 1) round(events(:,1:2)) events(:,3:5)];
end;

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [pcacoef] =computePca(namefile,ncoef)
% [pcacoef] =computePca(namefile,ncoef)

r=load(namefile);
if isfield(r,'peval')
    if isfield(r.peval, 'data_path')
        if isfield(r.peval,'data_file')
            d=load([r.peval.data_path '/' r.peval.data_file]);
        else
            d=load(r.peval.data_path); % there was some confusion what is path and waht is file
        end
    else
        error('Can not find the location of the data.')
    end
else
    d=load('dpixc'); % Data stored directlu in the directory
end

pcacoef = pca(reshape(d.dpixc,r.peval.nx*r.peval.ny,r.peval.nt),ncoef);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%

X = [ones(m, 1) X];

for j=1:m,

    %first layer propagation
    z2 = sigmoid(X(j,:) * Theta1');

    %bias to hidden layer
    z2 = [1 z2];

    %hidden layer propagation and getting max (candidate)
    [trash,p(j)] = max(sigmoid(z2 * Theta2'));
end;

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

end

author: vkosuri

repo: CourseraMachineLearning


text:
function dM = Bloch(t, M, Param, Pulse)
%BLOCH Bloch-McConnell ordinary differential equations.
%   dM = Bloch(t, M, Param, Pulse)
%
%   --args--
%   t: Function handle variable, represents the time.
%   M : magnetization vector [Mxf, Myf, Mzf, Mzr]
%   Param: Tissue parameter structure.
%
%          -fields-
%          M0f: Free pool equilibrium magnetization.
%          M0r: Restricted pool equilibrium magnetization. (M0r = F*M0f)
%          R1f: Free pool longitudinal relaxation rate.
%          R1r: Restricted pool longitudinal relaxation rate.
%          R2f: Free pool transverse relaxation rate.
%          kf: Free pool magnetization exchange rate.
%          kr: Restricted pool magnetization exchange rate. (= kf/F)
%          G: Lineshape value of the restricted pool (generated by
%             computeG.m)
%
%   Pulse: RF Pulse structure (generated by GetPulse.m).
%
%   Reference: R. Mark Henkelman, Xuemei Huang, Qing-San Xiang, G. J.
%               Stanisz, Scott D. Swanson, Michael J. Bronskill.
%               Quantitative Interpretation of Magnetization Transfer,Mag.
%               Res. in Med., 29, 759-766, Eqs. 1-6, (1993)
%
%
%   See also COMPUTEG, GETPULSE, BLOCH, BLOCHNOMT.
%

dM = zeros(4,1);

if (nargin < 4)
    omega  = 0;
    omega2 = 0;
    delta  = 0;
else
    omega  = Pulse.omega(t);
    omega2 = Pulse.omega2(t);
    delta  = Pulse.delta;
end

W = pi*Param.G.*omega2;

if isfield(Param,'T2f')
    R2f = 1./Param.T2f;
else
    R2f = Param.R2f;
end

if isfield(Param,'F')
    kf = Param.kr*Param.F;
    M0r = Param.M0f*Param.F;
else
    kf = Param.kf;
    M0r = Param.M0r;
end

dM(1) = -R2f*M(1) - 2*pi*delta*M(2);
dM(2) = -R2f*M(2) + 2*pi*delta*M(1) + omega*M(3);
dM(3) =  Param.R1f*(Param.M0f-M(3)) - kf*M(3) + Param.kr*M(4) - omega*M(2);
dM(4) =  Param.R1r*(M0r-M(4)) + kf*M(3) - Param.kr*M(4) - W*M(4);

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function CPD = deterministic_CPD(bnet, self, fname, pfail)
% DETERMINISTIC_CPD Make a tabular CPD representing a (noisy) deterministic function
%
% CPD = deterministic_CPD(bnet, self, fname)
% This calls feval(fname, pvals) for each possible vector of parent values.
% e.g., suppose there are 2 ternary parents, then pvals =
%  [1 1], [2 1], [3 1],   [1 2], [2 2], [3 2],   [1 3], [2 3], [3 3]
% If v = feval(fname, pvals(i)), then
%  CPD(x | parents=pvals(i)) = 1 if x==v, and = 0 if x<>v
% e.g., suppose X4 = X2 AND (NOT X3). Then
%    bnet.CPD{4} = deterministic_CPD(bnet, 4, inline('((x(1)-1) & ~(x(2)-1)) + 1'));
% Note that x(1) refers pvals(1) = X2, and x(2) refers to pvals(2)=X3
% See also boolean_CPD.
%
% CPD = deterministic_CPD(bnet, self, fname, pfail)
% will put probability mass 1-pfail on f(parents), and distribute pfail over the other values.
% This is useful for simulating noisy deterministic functions.
% If pfail is omitted, it is set to 0.
%

if nargin==0
  % This occurs if we are trying to load an object from a file.
  CPD = tabular_CPD(bnet, self);
  return;
elseif isa(bnet, 'deterministic_CPD')
  % This might occur if we are copying an object.
  CPD = bnet;
  return;
end

if nargin < 4, pfail = 0; end

ps = parents(bnet.dag, self);
ns = bnet.node_sizes;
psizes = ns(ps);
self_size = ns(self);

psucc = 1-pfail;

CPT = zeros(prod(psizes), self_size);
pvals = zeros(1, length(ps));
for i=1:prod(psizes)
  pvals = ind2subv(psizes, i);
  x = feval(fname, pvals);
  %fprintf('%d ', [pvals x]); fprintf('\n');
  if psucc == 1
    CPT(i, x) = 1;
  else
    CPT(i, x) = psucc;
    rest = mysetdiff(1:self_size, x);
    CPT(i, rest) = pfail/length(rest);
  end
end
CPT = reshape(CPT, [psizes self_size]);

CPD = tabular_CPD(bnet, self, 'CPT',CPT, 'clamped',1);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function pred = knntest(Xtest,class_t,X,class,K,dist_type,pret_type)

% prediction of new samples with calculated model
%
% pred = knnpred(Xtest,X,class,K,dist_type,pret_type)
%
% ------------ INPUT ---------------------------------------------------
% Xtest:        dataset to be predicted [n_test x p] n objects, p variables
% X:            training data matrix (n x p)
% class:        training class vector (n x 1)
% K:            number of neighbors
% dist_type:    'euclidean' Euclidean distance
%               'mahalanobis' Mahalanobis distance
%               'cityblock' City Block metric
%               'minkowski' Minkowski metric
%               'sm' Sokal-Michener
%               'jt' Jaccard Tanimoto
%               'gle' Gleason-Dice
%               'ct4' Consonni-Todeschini
%               'ac' Austin-Colwell
% pret_type:    'cent' cenering
%               'scal' variance scaling
%               'auto' for autoscaling (centering + variance scaling)
%               'rang' range scaling (0-1)
%               'fp'   fingerprints
%
% ------------ OUTPUT --------------------------------------------------
% pred is a structure conyaining
% class_pred    predicted class vector [n_test x 1]
% neighbors     list of k neighbors for each predicted sample [n_test x k]
%
% version 1.0 - september 2009
% Davide Ballabio
% Milano Chemometrics and QSAR Research Group
% www.disat.unimib.it/chm

% version 2.0 - February 2012
% Kamel Mansouri
% Milano Chemometrics and QSAR Research Group
% www.disat.unimib.it/chm

% data check
if length(class)~=size(X,1)
    disp('the class input should be for the training set')
    %class_tr=input('class tr');
    %class=evalin(WS,);
    %keyboard
end

[n,p] = size(Xtest);
[X_scal_train,param] = data_pretreatment(X,pret_type);
X_scal = test_pretreatment(Xtest,param);
Xd = [X_scal;X_scal_train];
% D = pdist(Xd,model.set.dist_type);
% D = squareform(D);
D = knn_calc_dist(X_scal_train,X_scal,dist_type,pret_type);
neighbors = zeros(n,K);
for i=1:n
    D_in = D(i,:);
    [d_tmp,n_tmp] = sort(D_in);
    neighbors(i,:) = n_tmp(1:K);
    d_neighbors = d_tmp(1:K);
    class_calc(i) = knnclass(class(neighbors(i,:)),d_neighbors,max(class),K);
    dc(i,:)=d_neighbors;
end

pred.neighbors  = neighbors;
pred.class_pred = class_calc';
pred.class_param = calc_class_param(pred.class_pred ,class_t);
pred.D=D;
pred.dc=dc;

author: kmansouri

repo: OPERA

save_path: github-repos/MATLAB/kmansouri-OPERA


text:
function w=getwords(s)
% w=getwords(s)
%
% Returns a cell array with the individual words in the string 's'

s = upper2space(s);
w = strread(s, '%s', 'delimiter', '-_., ');
w = strrep(w, '(', '');
w = strrep(w, ')', '');
w = strrep(w, '''', '');
w = strrep(w, '"', '');
w = strrep(w, ';', '');
w = strrep(w, '&', '');
w = strrep(w, '?', '');
w = strrep(w, '!', '');

author: CSAILVision

repo: LabelMeToolbox


text:
%demo for affine_fit
%Author: Adrien Leygue
%Date: December 3 2014
close all
clear all
figure;
%generate points that lie approximately in the Z=0 plane
N = 10;
[X,Y] = meshgrid(linspace(0,1,N));
XYZ_1 = [X(:) Y(:) 0.05*randn(N^2,1)];
plot3(XYZ_1(:,1),XYZ_1(:,2),XYZ_1(:,3),'r.');
hold on
%compute the normal to the plane and a point that belongs to the plane
[n_1,~,p_1] = affine_fit(XYZ_1);

%generate points that lie approximately in the Z=X plane
%the normal vector is
n_2_exact = [-sqrt(2)/2 0 sqrt(2)/2];
N = 12;
[X,Y] = meshgrid(linspace(0,1,N));
XYZ_2 = [X(:) Y(:) X(:)] + bsxfun(@times,0.05*randn(N^2,1),n_2_exact);
plot3(XYZ_2(:,1),XYZ_2(:,2),XYZ_2(:,3),'b.');

%compute the normal to the plane and a point that belongs to the plane
[n_2,V_2,p_2] = affine_fit(XYZ_2);

%plot the two points p_1 and p_2
plot3(p_1(1),p_1(2),p_1(3),'ro','markersize',15,'markerfacecolor','red');
plot3(p_2(1),p_2(2),p_2(3),'bo','markersize',15,'markerfacecolor','blue');

%plot the normal vector
quiver3(p_1(1),p_1(2),p_1(3),n_1(1)/3,n_1(2)/3,n_1(3)/3,'r','linewidth',2)
h = quiver3(p_2(1),p_2(2),p_2(3),n_2(1)/3,n_2(2)/3,n_2(3)/3,'b','linewidth',2)

%plot the two adjusted planes
[X,Y] = meshgrid(linspace(0,1,3));

%first plane
surf(X,Y, - (n_1(1)/n_1(3)*X+n_1(2)/n_1(3)*Y-dot(n_1,p_1)/n_1(3)),'facecolor','red','facealpha',0.5);

%second plane
%NB: if the plane is vertical the above method cannot be used, one should
%use the secont output of affine_fit which contains a base of the plane.
%this is illustrated below
%S1 and S2 are the coordinates of the plane points in the basis made of the
%columns ov V_2
[S1,S2] = meshgrid([-1 0 1]);
%generate the pont coordinates
X = p_2(1)+[S1(:) S2(:)]*V_2(1,:)';
Y = p_2(2)+[S1(:) S2(:)]*V_2(2,:)';
Z = p_2(3)+[S1(:) S2(:)]*V_2(3,:)';
%plot the plane
surf(reshape(X,3,3),reshape(Y,3,3),reshape(Z,3,3),'facecolor','blue','facealpha',0.5);

xlabel('x');
ylabel('y');
zlabel('z');
axis equal
%compute the angle between the planes in [0 90] degrees
angle = acosd(dot(n_1,n_2));
if angle>90
    angle = 180-angle;
end
angle

author: zouchuhang

repo: LayoutNet


text:
function model = load_svm_model(filename)
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

load(filename, 'weights', 'bias');
model = {weights; bias};
end

author: gidariss

repo: mrcnn-object-detection


text:
function plotroute(city, route, current_distance, temperature)
% PLOTROUTE
% PLOTROUTE(city, route, current_distance, temperature) plots the route and
% display current temperautre and distance.

global h;
cycle = route([1:end, 1]);
% update route
set(h,'Xdata',[city(cycle).long],'Ydata',[city(cycle).lat]);

% display current temperature and total distance
xlabel(sprintf('T = %6.1f        Total Distance = %6.1f', ...
                    temperature,                  current_distance));
drawnow

author: HuangCongQing


text:
%% Demo of gaimc - 'Graph Algorithms In Matlab Code'
% Matlab includes great algorithms to work with sparse matrices but does
% provide a reasonable set of algorithms to work with sparse matrices as
% graph data structures.  My other project -- MatlabBGL -- provides a
% high-performance solution to this problem by directly interfacing the
% Matlab sparse matrix data structure with the Boost Graph Library.
% That library, however, suffers from enormous complication because
% it must be compiled for each platform.  The Boost Graph Library
% heavily uses advanced C++ features that impair easy portability
% between platforms.  In contrast, the gaimc library is implemented in
% pure Matlab code, making it completely portable.
%
% The cost of the portability for this library is a 2-4x slowdown in
% the runtime of the algorithms as well as significantly fewer
% algorithms to choose from.

%% Sparse matrices as graphs
% To store the connectivity structure of the graph, gaimc uses the
% adjacency matrix of a graph.
%
% A graph is represented by a set of vertices and a set of
% edges between the vertices.  Often, we write $G = (V,E)$ to denote the
% graph, the set of vertices, and the set of edges, respectively.  In
% gaimc, like in my other package MatlabBGL, we represent graphs with their
% adjacency matrices.  This representation is handy in Matlab because
% Matlab is rather efficient at working with the large sparse matrices that
% typically arise as adjacency matrices.
%
% To convert from $G=(V,E)$ to an adjacency matrix, we identify each vertex
% with a row of the matrix via a bijective map.  The adjacency matrix is
% then a $|V| \times |V|$ matrix called A.  The entry A(i,j) = 1 for
% any edge between in $E$ and 0 otherwise.  Let's look at an example.

load_gaimc_graph('bfs_example');
graph_draw(A,xy,'labels',labels)
full(A)
labels'

%%
% This output means that vertex 'r' is row 1, vertex 's' is row 2 and
% because A(1,2) = 1, then there is an edge between them, just like in the
% picture.
%
% One funny property is that A(2,1) = 1 too!  So we actually have to store
% each edge twice in the adjacency matrix.  This might seem wasteful, but
% its hard to avoid as I've learned while working on graph algorithms.  So
% don't worry about it!  It also makes the generalization to directed
% graphs (below) easy.
%
% For more information about the adjacency matrix representation of a
% graph, see a standard book on graph algorithms.
%

%% Weighted and directed graphs
% Our previous case handled the situation for undirected graphs only.  To
% encode weighted and directed graphs, we use weighted and non-symmetric
% adjacency matrices.
%
% For a weighted matrix, A(i,j) = distance between i and j for most of the
% algorithms in gaimc.  But A(i,j) = 0 means there is no edge, and so
% sometimes things can get a little tricky to get what you want.
%
% For a directed graph, just set A(i,j) ~= A(j,i).  The adjacency matrix
% won't be symmetric, but that's what you want!
%
% To understand more, explore the examples or read up on adjacency matrices
% in graph theory books.

%% Loading helper
% To make loading our sample graphs easy, gaimc defines it's own function
% to load graphs.

load_gaimc_graph('dfs_example'); % loads one of our example graphs
whos

%%
% This helps make our examples work regardless of where the current
% directory lies.

%% Search algorithms
% The two standard graph search algorithms are depth first search and
% breadth first search.  This library implements both.

%%
% Load the example matrix from the Boost Graph Library
load_gaimc_graph('dfs_example');
figure(1); graph_draw(A,xy,'labels',labels);

%%
% Run a depth first search.  The output records the distance to the other
% vertices, except where the vertices are not reachable starting from the
% first node A.
d=dfs(A,1)

%%
% From this example, we see that vertices a-f are reachable from vertex a,
% but that verice g-i are not reachable.  Given the of the edges, this
% makes sense.

%%
% Let's look at breadth first search too, using a different example.
load_gaimc_graph('bfs_example');
figure(1); clf; graph_draw(A,xy,'labels',labels);

%%
% The breadth first search algorithm records the distance from the starting
% vertex to each vertex it visits in breadth first order.  This means it
% visits all the vertices in order of their distance from the starting
% vertex.  The d output records the distance, and the dt output records the
% step of the algorithm when the breadth first search saw the node.
[d dt] = bfs(A,2);
% draw the graph where the label is the "discovery time" of the vertex.
figure(1); clf; graph_draw(A,xy,'labels',num2str(dt));

%%
% Notice how the algorithm visits all vertices one edge away from the start
% vertex (0) before visiting those two edges away.

%% Shortest paths
% In the previous two examples, the distance between vertices was
% equivalent to the number of edges.  Some graphs, however, have specific
% weights, such as the graph of flights between airports.  We can use this
% information to build information about the _shortest path_ between two
% nodes in a network.

% Find the minimum travel time between Los Angeles (LAX) and
% Rochester Minnesota (RST).
load_gaimc_graph('airports')
A = -A; % fix funny encoding of airport data
lax=247; rst=355;

[d pred] = dijkstra(A,lax); % find all the shorest paths from Los Angeles.

fprintf('Minimum time: %g\n',d(rst));
% Print the path
fprintf('Path:\n');
path =[]; u = rst; while (u ~= lax) path=[u path]; u=pred(u); end
fprintf('%s',labels{lax});
for i=path; fprintf(' --> %s', labels{i}); end, fprintf('\n');

%% Minimum spanning trees
% A minimum spanning tree is a set of edges from a graph that ...
%
% This demo requires the mapping toolbox for maximum effect, but we'll do
% okay without it.

%%
% Our data comes from a graph Brendan Frey prepared for his affinity
% propagation clustering tool.  For 456 cities in the US, we have the mean
% travel time between airports in those cities, along with their latitude
% and longitude.
load_gaimc_graph('airports')

%%
% For some reason, the data is stored with the negative travel time between
% cities.  (I believe this is so that closer cities have larger edges
% between them.)  But for a minimum spanning tree, we want the actual
% travel time between cities.
A = -A;

%%
% Now, we just call MST and look at the result.
% T = mst_prim(A);
% This command means we can't run the demo, so it's commented out.

%%
% Oops, travel time isn't symmetric!  Let's just pick the longest possible
% time.
A = max(A,A');
T = mst_prim(A);
sum(sum(T))/2 % total travel time in tree

%%
% Well, the total weight isn't that helpful, let's _look_ at the data
% instead.
clf;
gplot(T,xy);

%%
% Hey!  That looks like the US!  You can see regional airports and get some
% sense of the overall connectivity.

%% Connected components
% The connected components of a network determine which parts of the
% network are reachable from other parts.  One of your first questions
% about any network should generally be: is it connected?
%
% There are two types of connected components: components and strongly
% connected components.  gaimc only implements an algorithm for the latter
% case, but that's okay!  It turns out it computes exactly the right thing
% for connected components as well.  The difference only occurs when the
% graph is undirected vs. directed.

load_gaimc_graph('dfs_example')
graph_draw(A,xy)

%%
% This picture shows there are 3 strongly connected components and 2
% connected components

% get the number of strongly connected components
max(scomponents(A))

%%

% get the number of connected components
max(scomponents(A|A'))  % we make the graph symmetric first by "or"ing each entry

%%
% Let's look at the vertices in the strongly connected components
cc = scomponents(A)
%%
% The output tells us that vertices 1,2,3,5,6 are in one strong component,
% vertex 4 is it's own strong component, and vertices 7,8,9 are in another
% one.  Remember that a strong component is all the vertices mutually
% reachable from a given vertex.  If you start at vertex 4, you can't get
% anywhere else!  That's why it is in a different component than vertices
% 1,2,3,5,6.

%%
% We also have a largest_component function that makes it easy to just get
% the largest connected component.
clf;
[Acc,f] = largest_component(A);
graph_draw(Acc,xy(f,:))
%%
% The filter variable f, tells us which vertices in the original graph made
% it into the largest strong component.  We can just apply that filter to
% the coordinates xy and reuse them for drawing the graph!

%% Statistics
% Graph statistics are just measures that indicate a property of the graph
% at every vertex, or at every edges.  Arguably the simplest graph
% statistic would be the average vertex degree.  Because such statistics
% are easy to compute with the adjaceny matrix in Matlab, they do not have
% special functions in gaimc.

%%
% Load a road network to use for statistical computations
load_gaimc_graph('minnesota');
gplot(A,xy);

%%
% Average vertex degree
d = sum(A,2);
mean(d)

%%
% So the average number of roads at any intersection is 2.5.  My guess is
% that many roads have artificial intersections in the graph structure that
% do not correspond to real intersections.  Try validating that hypothesis
% using the library!

%%
% Average clustering coefficients
ccfs = clustercoeffs(A);
mean(ccfs)
% The average clustering coefficient is a measure of the edge density
% throughout the graph.  A small value indicates that the network has few
% edges and they are well distributed throughout the graph.

%%
% Average core numbers
cn = corenums(A);
mean(cn)

%% Efficient repetition
% Every time a gaimc function runs, it converts the adjacency matrix into a
% set of compressed sparse row arrays.  These arrays yield efficient access
% to the edges of the graph starting at a particular vertex.  For many
% function calls on the same graph, this conversion process slows the
% algorithms.  Hence, gaimc also accepts pre-converted input, in which case
% it skips it's conversion.
%
% Let's demonstrate how this works by calling Dijkstra's algorithm to
% compute the shortest paths between all vertices in the graph.  The
% Floyd Warshall algorithm computes these same quantities more efficiently,
% but that would just be one more algorithm to implement and maintain.

%%
% Load and convert the graph.
load_gaimc_graph('all_shortest_paths_example');
A = spfun(@(x) x-min(min(A))+1,A); % remove the negative edges
As = convert_sparse(A);
%%
% Now, we'll run Dijkstra's algorithm for every vertex and save the result
% On my 2GHz laptop, this takes 0.000485 seconds.
n = size(A,1);
D = zeros(n,n);
tic
for i=1:n
    D(i,:) = dijkstra(As,i);
end
toc
%%
% Let's try it without the conversion to see if we can notice the
% difference in speed.
% On my 2GHz laptop, this takes 0.001392 seconds.
D2 = zeros(n,n);
tic
for i=1:n
    D2(i,:) = dijkstra(A,i);
end
toc
%%
% And just to check, let's make sure the output is the same.
isequal(D,D2)

author: ckczzj

repo: CHAN

save_path: github-repos/MATLAB/ckczzj-CHAN


text:
function C= cprintf(fmt, varargin)
%CPRINTF - Print array into cell array of string
%
%Synopsis:
%  cprintf(FORMAT, A, ...)
%
%Arguments
%  FORMAT - String that describes the format of the output as in sprintf
%
%Returns:
%  C - Cell array of strings
%
%Examples:
%  cprintf('chan #%d', 1:8)
%  cprintf('log(%d)= %f', 1:5, log(1:5))
%  cprintf('candidate #%d: %s', [1:3]', {'abraham','bebraham','zebraham'})

if isempty(varargin),
  C= {};
  return;
end

len= cellfun(@numel, varargin);
if any(diff(len)),
  error('all arguments must have the same number of elements');
end

N= len(1);
C= cell(N, 1);
Nv= length(varargin);
args= cell(1, Nv);
for n= 1:N,
  for m= 1:Nv,
    v= varargin{m};
    if iscell(v),
      args{m}= v{n};
    else
      args{m}= v(n);
    end
  end
  C{n}= sprintf(fmt, args{:});
end
C= reshape(C, size(varargin{1}));

author: bbci


text:
section_initializeVariablesAndSettings;

%% read images
targetFile = {'tomo.img', 'lena.bmp', 'bridge.img'};   % images
imageChoice = 1;
targetFile1 = strcat(targetFolder, '\', char(targetFile(imageChoice)));
image1 = readImage(targetFile1, 1);

imageChoice = 2;
targetFile2 = strcat(targetFolder, '\', char(targetFile(imageChoice)));
image2 = readImage(targetFile2, 1);

%% process
image1InGray = setImageToGray( image1 );
image2InGray = setImageToGray( image2 );
% adjust the size
newSize = [350 350];
image1InGray = imresize(image1InGray, newSize);
image2InGray = imresize(image2InGray, newSize);

% perform fft2, obtain magnitude and phase information
[image1InGray_fft_amplitudeOnly, image1InGray_fft_phaseOnly] = FFT_image( image1InGray );
[image2InGray_fft_amplitudeOnly, image2InGray_fft_phaseOnly] = FFT_image( image2InGray );

% interleaved the 2 FFT2 results
image1_amplitudeOnlyWithImage2_phaseOnly = image1InGray_fft_amplitudeOnly.*image2InGray_fft_phaseOnly;
image2_amplitudeOnlyWithImage1_phaseOnly = image2InGray_fft_amplitudeOnly.*image1InGray_fft_phaseOnly;
% take only the real parts
image1_amplitudeWithImage2_phase_realOnly = real(ifft2(image1_amplitudeOnlyWithImage2_phaseOnly));
image2_amplitudeWithImage1_phase_realOnly = real(ifft2(image2_amplitudeOnlyWithImage1_phaseOnly));

% recover image 1
imageInGray_InSpatialDomain = IFFT_image(image1InGray_fft_amplitudeOnly, image1InGray_fft_phaseOnly);

% display the FFT2 interleaved
section_plotResults;
% results commentary
section_resultsReporting;

author: Sable

repo: mcbench-benchmarks


text:
% DEMO3C - Chaining files using ncml and finding attributes using 'value4key'

echo('on')
o2_ = 'Oxygen';
time_ = 'esecs';

%% ---- Create the full pathname for an ncml file

scriptname = which('demo3c');  %% find the absolute path of a local file
fName = strrep(scriptname,'.m','.ncml') ; % create a path to the associated ncml file

if ~exist(fName,'file')
  error(strcat(fName,' does not exist.'));
end

%% ---- Access the OpenDAP datasets through an ncml file

nc = ncdataset(fName); % read the local ncml file
t = [nc.time(time_)];
o2 = [nc.data(o2_)];

%% ---- Find the units of Oxygen and label the display
attr = nc.attributes(o2_);
units = value4key(attr, 'units'); % Retrieve the units value
name = value4key(attr, 'long_name'); % Retrieve the long_name value

%% ---- Plot the data
figure;
plot(t, o2);...
datetick('x');...
ylabel([name ' [' units ']']);...
title({'M1 Mooring in Monterey Bay',nc.location},'interpreter','none');...
grid;...
shg
echo('off')

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function B = barymat(y, x, w, s, r, doFlip)
%BARYMAT  Barycentric Interpolation Matrix.
%   BARYMAT(Y, X, W), where Y is a column vector of length M and X and W are
%   column and row vectors of length N, respectively, returns the M*N matrix
%   which interpolates data from the grid X to the grid Y using the 2nd-kind
%   barycentric interpolation formula with barycentric weights W. If W is not
%   supplied it is assumed to be the weights for polynomial interpolation at a
%   2nd-kind Chebyshev grid: W(j) = (-1)^j, W([1, N]) = 0.5*W([1, N]).
%
%   BARYMAT(Y, X, W, S, R) is the same, where S = acos(Y) and R = acos(X). The
%   purpose of this is that Y(j) - X(k) can be more accurately computed in this
%   'theta space'. This is sometimes referred to as the 'trig trick' in spectral
%   collocation. BARYMAT(Y, X, W, S, R, 1) also performs the 'flipping trick',
%   which takes advantage of the fact that the smaller entries in R and S can be
%   computed more accurately. Note that X and Y should be symmetric about zero
%   for this work, and it is assumed that S and R are sorted in descending
%   order.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( isempty(y) || isempty(x) )
    % Nothing to do here!
    B = [];
    return
end

if ( size(x, 2) > 1 || size(y, 2) > 1 )
    error('CHEBFUN:barymat:dims', 'Inputs should be column vectors.');
end

% Obtain lengths of inputs:
N = length(x);
M = length(y);

% Nothing to do here!
if ( M == N && all(x == y) )
    B = eye(N);
    return
end

% Default to the Chebyshev barycentric weights:
if ( nargin < 3 || isempty(w) )
    w = ones(1, N);
    w(2:2:end) = -1;
    w([1, N]) = 0.5*w([1, N]);
else
    % Ensure w is a row vector:
    w = reshape(w, 1, N);
end

% Repmat(Y-X'):
if ( nargin < 5 )
    B = bsxfun(@minus, y, x.');
else
    % Use the 'trig trick' that y-x = cos(s)-cos(r) = 2*sin((s+r)/2)*sin((r-s)/2).
    B = 2*bsxfun(@(r, s) sin((s+r)/2).* sin((r-s)/2), r.', s);
end

% Construct the matrix:
if ( M >= 500 && N >= 1000 )            % <-- Experimentally determined.
    % Testing shows BSXFUN is faster in this regime
    B = bsxfun(@rdivide, w, B);         % w(k)/(y(j)-x(k))
    B = bsxfun(@rdivide, B, sum(B, 2)); % Normalisation.
else
    % Else use FOR loops
    for k = 1:N
        B(:,k) = w(k)./B(:,k);          % w(k)/(y(j)-x(k))
    end
    c = 1./sum(B, 2);                   % Normalisation.
    for j = 1:M
        B(j,:) = B(j,:)*c(j);
    end
end

% Where points coincide there will be division by zeros (as with bary.m).
% Replace these entries with the identity:
B(isnan(B)) = 1;

% Flipping trick:
if ( nargin > 5 && doFlip )
    ii = logical(rot90(tril(ones(M, N)), 2));
    ii = fliplr(ii);
    rot90D = rot90(B, 2);
    B(ii) = rot90D(ii);
    B(isnan(B)) = 1;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function p_eff = gp_peff(gp, x, y, varargin);
%GP_PEFF  The effective number of parameters in GP model with focus
%         on latent variables
%
%  Description
%    P_EFF = EP_PEFF(GP, X, Y) Takes the Gaussian process structure
%    GP, training inputs X and training outputs and returns the
%    effective number of parameters as defined by Spiegelhalter et
%    al. (2002).
%
%    NOTE! The effective number of parameters is evaluated with
%    focus on latent variable f. This means that the parameters th
%    (parameters of covariance function and likelihood) are
%    considered fixed. (See Spiegelhalter et al (2002) for
%    discussion on the parameters in focus in Bayesian model).
%    Thus, the returned p_eff tells the effective number of latent
%    variables. This statistics is important for example when
%    assessing the goodness of Laplace or EP approximation in case
%    of non-Gaussian likelihood (See Vanhatalo et al for
%    discussion).
%
%    If you want to evaluate the effective number of parameters
%    with focus on parameters, see GP_DIC.
%
%    The effective number of parameters is approximated as follows:
%        p_eff = n - trace( K\C ),
%
%    where K is the prior covariance matrix and C the posterior
%    covariance matrix. This approximation is introduced by
%    Spiegelhalter et al. (2002) in equation (16). If the
%    likelihood is non-Gaussian and gp.latent_method is either
%    Laplace or EP, then C is the Laplace or EP approximation for
%    the posterior covariance.
%
%  References:
%    Spiegelhalter, Best, Carlin and van der Linde (2002).
%    Bayesian measures of model complexity and fit. J. R.
%    Statist. Soc. B, 64(4):583-639.
%
%    Vanhatalo, J., Pietil�inen V. and Vehtari, A. (2010).
%    Approximate inference for disease mapping with sparse
%    Gaussian processes. Statistics in Medicine, 29(15):1580-1607.
%
%  See also
%    GP_DIC, DEMO_MODELASSESMENT1
%
%
% Copyright (c) 2009-2010 Jarno Vanhatalo

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

  ip=inputParser;
  ip.FunctionName = 'GP_PEFF';
  ip.addRequired('gp',@isstruct);
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.parse(gp, x, y, varargin{:});
  z=ip.Results.z;

  tn = size(x,1);

    if isfield(gp.lik.fh,'trcov')
      % a Gaussian likelihood

        switch gp.type
          case 'FULL'

            [K, C] = gp_trcov(gp, x);
            L = chol(C);
            p_eff = trace( L\(L'\K) );

          case 'FIC'
            u = gp.X_u;
            m = size(u,1);
            % Turn the inducing vector on right direction
            if size(u,2) ~= size(x,2)
                u=u';
            end
            % Calculate some help matrices
            [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
            K_fu = gp_cov(gp, x, u);   % f x u
            K_uu = gp_trcov(gp, u);     % u x u, noiseles covariance K_uu
            Luu = chol(K_uu)';

            % Evaluate Lambda (La) for specific model
            % Q_ff = K_fu*inv(K_uu)*K_fu'
            % Here we need only the diag(Q_ff), which is evaluated below
            B=Luu\(K_fu');
            Qv_ff=sum(B.^2)';
            Lav = Kv_ff-Qv_ff;
            Lav2 = Cv_ff-Qv_ff;

            iLaKfu = zeros(size(K_fu));  % f x u,
            n = size(x,1);
            for i=1:n
                iLaKfu(i,:) = K_fu(i,:)./Lav2(i);  % f x u
            end
            A = K_uu+K_fu'*iLaKfu;
            A = (A+A')./2;

            L = iLaKfu/chol(A);

            p_eff = sum(Lav./Lav2) + sum(sum( repmat(Lav2,1,m).\B'.*B' - L.*(L'*B'*B)' - L.*(L'.*repmat(Lav',m,1))', 2));

% $$$             % Check the result using full matrices
% $$$             C = B'*B + diag(Lav2);
% $$$             K = B'*B + diag(Lav);
% $$$             L = chol(C);
% $$$             p_eff = trace( L\(L'\K) );

          case {'PIC' 'PIC_BLOCK'}
            u = gp.X_u;
            ind = gp.tr_index;
            if size(u,2) ~= size(x,2)
                u=u';
            end

            % Calculate some help matrices
            [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
            K_fu = gp_cov(gp, x, u);         % f x u
            K_uu = gp_trcov(gp, u);    % u x u, noiseles covariance K_uu
            Luu = chol(K_uu)';

            % Evaluate the Lambda (La) for specific model
            % Q_ff = K_fu*inv(K_uu)*K_fu'
            % Here we need only the diag(Q_ff), which is evaluated below
            B=Luu\K_fu';
            iLaKfu = zeros(size(K_fu));  % f x u
            for i=1:length(ind)
                Qbl_ff = B(:,ind{i})'*B(:,ind{i});
                [Kbl_ff, Cbl_ff] = gp_trcov(gp, x(ind{i},:));
                La{i} = Kbl_ff - Qbl_ff;
                La2{i} = Cbl_ff - Qbl_ff;
                iLaKfu(ind{i},:) = La2{i}\K_fu(ind{i},:);
            end
            A = K_uu+K_fu'*iLaKfu;
            A = (A+A')./2;            % Ensure symmetry
            L = iLaKfu/chol(A);

            p_eff = sum(sum(- L.*(L'*B'*B)',2));
            for i=1:length(ind)
                LLa2 = chol(La2{i});
                p_eff = p_eff + trace(LLa2\(LLa2'\La{i})) + trace( LLa2\(LLa2'\B(:,ind{i})'*B(:,ind{i})) - L(ind{i},:)*L(ind{i},:)'*La{i} );
            end

          case 'CS+FIC'
            u = gp.X_u;
            m = size(u,1);
            % Turn the inducing vector on right direction
            if size(u,2) ~= size(x,2)
                u=u';
            end

            % Indexes to all non-compact support and compact support covariances.
            cf1 = [];
            cf2 = [];

            ncf = length(gp.cf);
            % Loop through all covariance functions
            for i = 1:ncf
                % Non-CS covariances
                if ~isfield(gp.cf{i},'cs')
                    cf1 = [cf1 i];
                    % CS-covariances
                else
                    cf2 = [cf2 i];
                end
            end

            [Kv_ff, Cv_ff] = gp_trvar(gp, x, cf1);  % f x 1  vector
            K_fu = gp_cov(gp,x,u,cf1);         % f x u
            K_uu = gp_trcov(gp,u,cf1);    % u x u, noiseles covariance K_uu
            K_uu = (K_uu+K_uu')./2;     % ensure the symmetry of K_uu

            Kcs = gp_trcov(gp, x, cf2);
            Luu = chol(K_uu)';
            B=Luu\(K_fu');       % u x f
            Qv_ff=sum(B.^2)';
            Lav = Cv_ff-Qv_ff;   % f x 1, Vector of diagonal elements
            Lav2 = Kv_ff-Qv_ff;   % f x 1, Vector of diagonal elements

            La = sparse(1:tn,1:tn,Lav,tn,tn) + Kcs;
            La2 = sparse(1:tn,1:tn,Lav2,tn,tn) + Kcs;

            iLaKfu = La\K_fu;
            A = K_uu+K_fu'*iLaKfu;
            A = (A+A')./2;     % Ensure symmetry
            L = iLaKfu/chol(A);

            LLa2 = L'*La2;
            LaB = La\B';
            LBB = L'*B'*B;

            VD = ldlchol(La);
            spiLa = spinv(VD,1);

            p_eff = sum( sum(LaB.*B',2) + sum(spiLa.*La2,2) - sum(L.*LBB',2) - sum(L.*LLa2',2) );

% $$$
% $$$             C = B'*B + Kcs + diag(Lav);
% $$$             K = B'*B + Kcs + diag(Lav2);
% $$$
% $$$             L = chol(C);
% $$$             p_eff = trace( L\(L'\K) );
        end

    else
      % ============================
      % A non Gaussian likelihood
      % ============================

        switch gp.type
          case 'FULL'
            switch gp.latent_method
              case 'EP'
                %[e, edata, eprior, tautilde, nutilde, L] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
                [e, edata, eprior, p] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
                [tautilde, L] = deal(p.tautilde, p.L);

                % The prior variance
                K=gp_trcov(gp,x);

                if all(tautilde > 0) && ~isequal(gp.latent_opt.optim_method, 'robust-EP')
                    sqrttautilde = sqrt(tautilde);
                    Stildesqroot = sparse(1:tn, 1:tn, sqrttautilde, tn, tn);

                    if issparse(L)
                        p_eff = trace(Stildesqroot*ldlsolve(L, Stildesqroot*K));
                    else
                        p_eff = trace(Stildesqroot*(L'\(L\Stildesqroot)*K));
                    end
                else
                    C = L*L';
                    L = chol(K);
                    Varf = tn - trace(L\(L'\C));
                end

              case 'Laplace'
                %[e, edata, eprior, f, L] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
                [e, edata, eprior, p] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
                [f, L] = deal(p.f, p.L);

                W = -gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);

                % Evaluate the prior variance
                K = gp_trcov(gp,x);

                if W >= 0
                    if issparse(K) && issparse(L)

                        sqrtW = sparse(1:tn, 1:tn, sqrt(W), tn, tn);
                        sqrtWK = sqrtW*K;
                        p_eff = trace(sqrtW*ldlsolve(L,sqrtWK));
                    else
                        W = diag(W);
                        p_eff = trace(sqrt(W)*(L'\(L\(sqrt(W)*K))));
                    end
                else
                    C = L*L';
                    L = chol(K);
                    p_eff = tn - trace(L\(L'\C));
                end
            end

          case 'FIC'
            u = gp.X_u;
            m = size(u,1);
            % Calculate some help matrices
            [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
            K_fu = gp_cov(gp, x, u);   % f x u
            K_uu = gp_trcov(gp, u);     % u x u, noiseles covariance K_uu
            Luu = chol(K_uu)';
            B=Luu\(K_fu');
            Qv_ff=sum(B.^2)';
            Lav = Kv_ff-Qv_ff;

            switch gp.latent_method
              case 'EP'
                %[e, edata, eprior, tautilde, nutilde, L, La, b] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
                %[e, edata, eprior, p] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
                [L, La] = deal(p.L, p.La2);

                k = gp_trvar(gp,x);

                p_eff = sum( sum( (repmat(1./La,1,m).*B').*B',2) ) - sum(sum(L.*((L'*B')*B)',2));
                Lav = k - sum(B.^2)';
                p_eff = p_eff + sum(Lav./La) - sum(sum(L.*L,2).*Lav);

% $$$                 C = diag(1./La) - L*L';
% $$$                 K = B'*B + diag(k - sum(B.^2)');
% $$$
% $$$                 p_eff = trace(C*K);

              case 'Laplace'
                %[e, edata, eprior, f, L, a, La2] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
                [e, edata, eprior, p] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
                [f, La2] = deal(p.f, p.La2);

                W = -gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);
                La = W.*Lav;
                Lahat = 1 + La;
                sqrtW = sqrt(W);
                B = (repmat(sqrtW,1,m).*K_fu);

                % Components for (I + W^(1/2)*(Qff + La2)*W^(1/2))^(-1) = Lahat^(-1) - L2*L2'
                B2 = repmat(Lahat,1,m).\B;
                A2 = K_uu + B'*B2; A2=(A2+A2)/2;
                L2 = B2/chol(A2);

                BB=Luu\(K_fu');
                BB2=B/Luu';

                p_eff = sum(W./Lahat.*Lav);
                p_eff = p_eff + sum(sqrtW .* (sum((repmat(Lahat,1,m).\BB2).*BB',2)...
                                             - sum(L2.*(L2'.*repmat(sqrtW'.*La2',m,1))',2)...
                                             - sum(L2.*(L2'*BB2*BB)',2)) );

% $$$                 K = BB'*BB + diag(Lav);
% $$$                 sW = diag(sqrt(W));
% $$$                 B = eye(size(K)) + sW*K*sW;
% $$$                 L = chol(B)';
% $$$
% $$$                 W = diag(W);
% $$$                 p_eff = trace(sqrt(W)*(L'\(L\(sqrt(W)*K))));
            end

          case {'PIC' 'PIC_BLOCK'}
            u = gp.X_u;
            K_fu = gp_cov(gp, x, u);
            K_uu = gp_trcov(gp, u);
            K_uu = (K_uu+K_uu')./2;

            m = size(u,1);
            ind = gp.tr_index;
            Luu = chol(K_uu)';
            B=Luu\(K_fu');

            switch gp.latent_method
              case 'EP'

                %[e, edata, eprior, tautilde, nutilde, L, La, b] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
                [e, edata, eprior, p] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
                [L, La] = deal(p.L, p.La2);

                p_eff = - sum(sum(L.*((L'*B')*B)',2));

                for i=1:length(ind)
                    La1 = gp_trcov(gp, x(ind{i},:)) - B(:,ind{i})'*B(:,ind{i});
                    p_eff = p_eff + trace(La{i}\B(:,ind{i})'*B(:,ind{i}));
                    p_eff = p_eff + trace(La{i}\La1);
                    p_eff = p_eff - trace(L(ind{i},:)*L(ind{i},:)'*La1);
                end

              case 'Laplace'
                %[e, edata, eprior, f, L, a, La2] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
                [e, edata, eprior, p] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
                [f, La2] = deal(p.f, p.La2);

                iKuuKuf = K_uu\K_fu';

                % Evaluate the variance

                W = -gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);
                sqrtW = sqrt(W);

                % Components for (I + W^(1/2)*(Qff + La2)*W^(1/2))^(-1) = Lahat^(-1) - L2*L2'
                for i=1:length(ind)
                    La{i} = diag(sqrtW(ind{i}))*La2{i}*diag(sqrtW(ind{i}));
                    Lahat{i} = eye(size(La{i})) + La{i};
                    LLahat{i} = chol(Lahat{i});
                end
                sKfu = (repmat(sqrt(W),1,m).*K_fu);
                for i=1:length(ind)
                    iLasKfu(ind{i},:) = Lahat{i}\sKfu(ind{i},:);
                end
                A2 = K_uu + sKfu'*iLasKfu; A2=(A2+A2)/2;
                L2 = iLasKfu/chol(A2);

                p_eff = -sum(sqrtW.*sum(L2.*(L2'*(repmat(sqrtW,1,m).*B')*B)',2));
                for i=1:length(ind)
                    dsqrtW = diag(sqrtW(ind{i}));
                   p_eff = p_eff + trace( dsqrtW*(LLahat{i}\(LLahat{i}'\(dsqrtW*La2{i}))));
                   p_eff = p_eff + trace( dsqrtW*(LLahat{i}\(LLahat{i}'\(dsqrtW*B(:,ind{i})'*B(:,ind{i})))));
                   p_eff = p_eff - trace(dsqrtW*L2(ind{i},:)*L2(ind{i},:)'*dsqrtW*La2{i});
                end

% $$$                 K = B'*B;
% $$$                 C = -L2*L2';
% $$$                 for i=1:length(ind)
% $$$                    K(ind{i},ind{i}) =  K(ind{i},ind{i}) + La2{i};
% $$$                    C(ind{i},ind{i}) =  C(ind{i},ind{i}) + inv(Lahat{i});
% $$$                 end
% $$$
% $$$                 p_eff = trace (diag(sqrtW)*C*diag(sqrtW)*K) ;

            end

          case 'CS+FIC'
            u = gp.X_u;
            [n,nin]=size(x);
            % Indexes to all non-compact support and compact support covariances.
            cf1 = [];
            cf2 = [];

            ncf = length(gp.cf);
            % Loop through all covariance functions
            for i = 1:ncf
                % Non-CS covariances
                if ~isfield(gp.cf{i},'cs')
                    cf1 = [cf1 i];
                    % CS-covariances
                else
                    cf2 = [cf2 i];
                end
            end

            K_fu = gp_cov(gp,x,u,cf1);         % f x u
            K_uu = gp_trcov(gp,u,cf1);    % u x u, noiseles covariance K_uu
            K_uu = (K_uu+K_uu')./2;     % ensure the symmetry of K_uu

            Kcs = gp_trcov(gp, x, cf2);
            Luu = chol(K_uu)';
            B=Luu\(K_fu');

            switch gp.latent_method
              case 'EP'

                %[e, edata, eprior, tautilde, nutilde, L, La, b] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
                [e, edata, eprior, p] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
                [L, La] = deal(p.L, p.La2);

                k = gp_trvar(gp,x,cf1);
                Lav = k - sum(B.^2)';
                La1 = Kcs + sparse(1:n, 1:n, Lav, n, n);

                issparse(La)
                VD = ldlchol(La);
                siLa = spinv(VD,1);
                p_eff = sum(sum(ldlsolve(VD,B').*B')) + sum(sum(siLa.*La1)) - sum(sum(L.*((L'*B')*B)',2));
                p_eff = p_eff - sum(sum(L.*(L'*La1)',2));

% $$$                 C = inv(La) - L*L';
% $$$                 K = B'*B + diag(k - sum(B.^2)') + Kcs;
% $$$
% $$$                 p_eff = trace(C*K);

              case 'Laplace'
                %[e, edata, eprior, f, L, a, La2] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
                [e, edata, eprior, p] = gpla_e(gp_pak(gp), gp, x, y, 'z', z);
                [f, La2] = deal(p.f, p.La2);

                W = -gp.lik.fh.llg2(gp.lik, y, f, 'latent', z);
                sqrtW = sparse(1:tn,1:tn,sqrt(W),tn,tn);
                Lahat = sparse(1:tn,1:tn,1,tn,tn) + sqrtW*La2*sqrtW;
                B = sqrtW*K_fu;

                % Components for (I + W^(1/2)*(Qff + La2)*W^(1/2))^(-1) = Lahat^(-1) - L2*L2'
                B2 = Lahat\B;
                A2 = K_uu + B'*B2; A2=(A2+A2)/2;
                L2 = B2/chol(A2);

                BB = Luu\(K_fu');
                BB2 = B/Luu';

                VD = ldlchol(Lahat);
                spiLahat = spinv(VD,1);

                p_eff = sum(sum(sqrtW*spiLahat*sqrtW.*La2,2)) + sum(sqrtW * (sum(ldlsolve(VD, BB2).*BB',2) - sum(L2.*(L2'*sqrtW*La2)',2) - sum(L2.*(L2'*BB2*BB)',2)) );

            end

        end
    end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function test_partial_basis_1d(N, M, K)

if nargin < 1
    N = 2048;
end

if nargin < 2
    M = floor(N/8);
end

if nargin < 3
    K = floor(M/5);
end

% prepare demo signal
X = spx.data.synthetic.SparseSignalGenerator()

end

author: indigits

repo: sparse-plex


text:
function w = ymdf_to_weekday_julian2 ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_WEEKDAY_JULIAN2 returns the weekday of a Julian YMDF date.
%
%  Discussion:
%
%    This routine computes the day of the week from the date in
%    the Julian calendar, that is, the calendar in force before the
%    Gregorian calendar, in which every fourth year was a leap year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm A,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 307-308.
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, integer W, the day of the week of the date.
%    The days are numbered from Sunday through Saturday, 1 through 7.
%
  m_prime = mod ( 9 + m, 12 );
  q = floor ( m_prime / 10 );
  z = floor ( ( 13 * m_prime + 2 ) / 5 );
  t = 28 * m_prime + z + d - 365 * q + 59;

  c = i4_wrap ( t, 1, 7 );

  y_prime = y - q;
  v = floor ( y / 4 ) - floor ( y_prime / 4 );
  p = y + floor ( y / 4 ) + 4 - v;
  n = 7 - mod ( p, 7 );

  w = i4_wrap ( c + 1 - n, 1, 7 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function F = diff(d, varargin)
%DIFF      Differentiation operator.
%   This function is deprecated. Use OPERATORBLOCK.DIFF.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

F = linop( operatorBlock.diff(double(d), varargin{:}) );

end


text:
classdef mp_idx_manager < handle
%MP_IDX_MANAGER  MATPOWER Index Manager abstract class
%
%   A MATPOWER Index Manager object can be used to manage the indexing of
%   various named and indexed blocks of various set types, such as variables,
%   constraints, etc. This class helps keep track of the ordering and
%   indexing of the various blocks as they are added to the object.
%
%   The types of named sets to be managed by the class are defined by the
%   DEF_SET_TYPES method, which assigns a struct to the 'set_types' field.
%
%   Properties:
%       set_types   - a struct defined by DEF_SET_TYPES method
%       userdata    - a struct containing arbitrary data added by the user
%
%   Private Methods
%       add_named_set - Adds named subset of a particular type to the
%           object.
%
%       def_set_types - (must be implemented in the subclass) Returns a
%           struct defining the various set types, where the key is the
%           set type name, which must also be declared as a property in
%           the object's class, and the value is a string name used for
%           display purposes.
%
%           E.g.
%               function obj = def_set_types(obj)
%                   obj.set_types = struct(...
%                           'var', 'variable', ...
%                           'lin', 'linear constraint' ...
%                       );
%               end
%
%       init_set_types - Initializes the structures needed to track the
%           ordering and indexing of each set type and can be overridden
%           to initialize any additional data to be stored with each block
%           of each set type. Ideally, this would be called at the end
%           of the MP_IDX_MANAGER constructor, but this is not possible
%           in Octave 5.2 and earlier due to a bug related to altering
%           fields of an object not yet fully constructed.
%           Fixed in Octave v6.x: https://savannah.gnu.org/bugs/?52614
%           The workaround is to be sure that this method is called by
%           some subclass method(s) after object construction, but before
%           any other use (e.g. in display() and add_*() methods).
%
%       valid_named_set_type - Returns a label for the given named set
%           type if valid, empty otherwise.
%
%   Public Methods
%       add_<SET-TYPE> - (must be implemented in the subclass) Used to
%           add a block of the given set type to the object, using the
%           private ADD_NAMED_SET method internally followed by the
%           handling of any data specific to that set type.
%
%       copy - Makes a shallow copy of the object.
%
%       describe_idx - Identifies indices of a given set type.
%           E.g. variable 361 corresponds to Pg(68)
%
%       display - (must be implemented in the subclass) Displays the
%           object (called automatically when you omit the semicolon at
%           the command-line).
%
%       display_set - Prints to screen the indexing details for the
%           specified set type. Intended to be called by DISPLAY method.
%
%       get_idx - Returns index structure(s) for specified set type(s),
%           with starting/ending indices and number of elements for
%           each named (and optionally indexed) block.
%
%       get_userdata - Retreives values of user data stored in the object.
%
%       get - Return the value of any individual field.
%
%       getN - Returns the number of elements of any given set type,
%           optionally for a single named block.
%
%       init_indexed_name - Initializes the dimensions for an indexed
%           named set.
%
%       params_<SET-TYPE> - (must be implemented in the subclass)
%           Returns set-type-specific data for a given type.
%
%   The following is the structure of the data in the object, using a set
%   type named 'var' for illustration. Each field of .idx or .data is a
%   struct whose field names are the names of the corresponding blocks of
%   elements of that type (e.g. variables, constraints, etc.). They are
%   found in order in the corresponding .order field. The description next
%   to these fields gives the meaning of the value for each named sub-field.
%   E.g. obj.var.data.v0.Pg contains a vector of initial values for the 'Pg'
%   block of variables.
%
%   obj
%       .var        - data for 'var' set type, e.g. variable sets that
%                     make up the full optimization variable x
%           .idx
%               .i1 - starting index within x
%               .iN - ending index within x
%               .N  - number of elements in this variable set
%           .N      - total number of elements in x
%           .NS     - number of variable sets or named blocks
%           .data   - additional set-type-specific data for each block
%           .order  - struct array of names/indices for variable
%                     blocks in the order they appear in x
%               .name   - name of the block, e.g. Pg
%               .idx    - indices for name, {2,3} => Pg(2,3)
%       .userdata   - any user defined data
%           .(user defined fields)

%   MP-Opt-Model
%   Copyright (c) 2008-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%    es = struct();

    properties
        userdata = struct();
        set_types
    end     %% properties

    methods
        %% constructor
        function obj = mp_idx_manager(s)
            if nargin > 0
                if isa(s, 'mp_idx_manager')
                    %% this copy constructor will not be inheritable under
                    %% Octave until the fix has been included for:
                    %%      https://savannah.gnu.org/bugs/?52614
                    if have_feature('octave')
                        s1 = warning('query', 'Octave:classdef-to-struct');
                        warning('off', 'Octave:classdef-to-struct');
                    end
                    props = fieldnames(s);
                    if have_feature('octave')
                        warning(s1.state, 'Octave:classdef-to-struct');
                    end
                    for k = 1:length(props)
                        obj.(props{k}) = s.(props{k});
                    end
                elseif isstruct(s)
                    props = fieldnames(obj);
                    for k = 1:length(props)
                        if isfield(s, props{k})
                            obj.(props{k}) = s.(props{k});
                        end
                    end
                else
                    error('@mp_idx_manager/mp_idx_manager: input must be an ''mp_idx_manager'' object or a struct');
                end
            end

            obj.def_set_types();

            %% The INIT_SET_TYPES() method should ideally be (1) called here
            %% in the MP_IDX_MANAGER constructor and (2) skipped if constructed
            %% from an existing object, neither of which work in Octave 5.2
            %% and earlier due to inheritance issues in constructors:
            %%   https://savannah.gnu.org/bugs/?52614
            %%
            %% WORKAROUND: In some subclass method(s), call INIT_SET_TYPES()
            %%             automatically as needed after construction, but
            %%             before use, checking first to ensure that it
            %%             hasn't already been called.
%             if isempty(obj.????)        %% skip if already initialized (e.g.
%                 obj.init_set_types();   %% constructed from existing object)
%             end
        end

        function obj = init_set_types(obj)
            %% base data struct for each type
            es = struct();
            ds = struct( ...
                'idx', struct( ...
                    'i1', es, ...
                    'iN', es, ...
                    'N', es ), ...
                'N', 0, ...
                'NS', 0, ...
                'order', struct( ...
                    'name', [], ...
                    'idx', [] ), ...
                'data', es );

            %% initialize each (set_type) field with base data structure
            for f = fieldnames(obj.set_types)'
                obj.(f{1}) = ds;
            end
        end

        function new_obj = copy(obj)
            %% make shallow copy of object
            new_obj = eval(class(obj));  %% create new object
            if have_feature('octave')
                s1 = warning('query', 'Octave:classdef-to-struct');
                warning('off', 'Octave:classdef-to-struct');
            end
            props = fieldnames(obj);
            if have_feature('octave')
                warning(s1.state, 'Octave:classdef-to-struct');
            end
            for k = 1:length(props)
                new_obj.(props{k}) = obj.(props{k});
            end
        end

        function display_set(obj, stype, sname)
            if nargin < 3
                sname = obj.set_types.(stype);
            end
            st = obj.(stype);    %% data for set type of interest
            if st.NS
                fmt = '%-26s %6s %8s %8s %8s\n';
                fprintf(fmt, sname, 'name', 'i1', 'iN', 'N');
                fprintf(fmt, repmat('=', 1, length(sname)), '------', '-----', '-----', '------');
                idx = st.idx;
                fmt = '%10d:%22s %8d %8d %8d\n';
                for k = 1:st.NS
                    name = st.order(k).name;
                    if isempty(st.order(k).idx)
                        fprintf(fmt, k, name, idx.i1.(name), idx.iN.(name), idx.N.(name));
                    else
                        vsidx = st.order(k).idx;
                        str = '%d'; for m = 2:length(vsidx), str = [str ',%d']; end
                        s = substruct('.', name, '()', vsidx);
                        nname = sprintf(['%s(' str, ')'], name, vsidx{:});
                        fprintf(fmt, k, nname, ...
                                subsref(idx.i1, s), subsref(idx.iN, s), subsref(idx.N, s));
                    end
                end
                fmt = sprintf('%%10d = %%s.NS%%%dd = %%s.N\\n\\n', 35-length(stype));
                fprintf(fmt, st.NS, stype, st.N, stype);
            else
                fprintf('%-26s  :  <none>\n', sname);
            end
        end

        obj = add_named_set(obj, set_type, name, idx, N, varargin)

        label = describe_idx(obj, set_type, idxs)

        varargout = get_idx(obj, varargin)

        rv = get_userdata(obj, name)

        val = get(obj, varargin)

        N = getN(obj, set_type, name, idx)

        obj = init_indexed_name(obj, set_type, name, dim_list)

        str = valid_named_set_type(obj, set_type)
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower