text:
% evaluate simulated data - ICA and NMF
% V ~ WH
% V -> N_pix x N_t
% W -> N_pix x N_comp - xth pixel of the ith components
% H -> N_copm x N_t - contribution of the i-th component in the time t

function separNMFICA(sep0, offset0, path_data, path_res, prename, niter, savethis, initval, sep_how)

if ~exist('initval', 'var')
    initval = 0;
end

if ~exist('sep_how', 'var')
    sep_how = 'in';
end
ncomp = 3; %number of components to be separated + background

for rr = 1: length(offset0)
    for ll=1 : length(sep0)
        namedir = [prename num2str(100*sep0(ll)) 'offset_' num2str(offset0(rr))];
        cd ([path_data namedir])
        for mm=1:niter
            %reads the first...
            namefile = [namedir '-iter_' num2str(mm)];
            load ([namefile '.mat'])
            %cat the folowing,,,
            p.catitervec=[1];
            [dpixc, dveccr, blinkmat, p] = catsimul(namedir, p.catitervec);
            if sum(sep_how == 'i')>0 %ICA
                [icasig{mm}, A{mm}, W{mm}] = fastica (dveccr, 'numOfIC', ncomp, 'g', 'tanh');
                icapixICA{mm} = reshape(A{mm},p.nx, p.ny, ncomp);
            end

            if sum(sep_how == 'n')>0 %NMF
                if initval
                    % background estimation:
                    %[out, bg(mm), bg_im]=backgroundoffset(dpixc);
                    [out, bg(mm), bg_im]=backgroundoffset(dpixc, 'no', 5, 20, 8); %empirical values...

                    dvec_bg = ones(p.nx*p.ny, 1);
                    %                     dvec_bg = p.offset*ones(1, p.nx*p.ny); %changed for
                    %                     offset 10...

                    dvec_ind = squeeze(reshape(double(array2im(dpixc_ind)), p.nx*p.ny, 1, 2)); % vectors of resized images
%                     sum_dvec_ind = sum(dvec_ind, 1);
%                     dvec_ind = dvec_ind./repmat(sum_dvec_ind, p.nx*p.ny,1); %normlaized

                    %                     winit = [f*dvec_ind'; dvec_bg];       %original 'true' points + background
% % %                     winittmp = [dvec_ind, dvec_bg];
                    winittmp = [rand(size(dvec_ind)), dvec_bg];
                    sumw = sum(winittmp,1);
                    winit = winittmp./repmat(sumw, p.nx*p.ny, 1); %normalized to 1
                    f = mean(dveccr(:)-bg(mm))/mean(mean(winit(:, 1:2))); %ration of the data/psf
                    %                     winit = [rand(ncomp, p.nx*p.ny); dvec_bg];
                    blinkmatrand = rand(ncomp-1, p.Nt); %uniform random;
                    hinit = [f*blinkmatrand; bg(mm)*sumw(ncomp)*ones(1, p.Nt)];             %random weights will be assigned to firts two and bg fixed
% % %                     hinit = [f*blinkmat./repmat(mean(blinkmat,2),1,size(blinkmat,2)); bg(mm)*sumw(ncomp)*ones(1, p.Nt)];             %random weights will be assigned to firts two and bg fixed

                    %                    [w{mm},h{mm}, wtrace{mm},wtrace{mm}]=nmf_test(double(dveccr'),ncomp+1,1,hinit,winit, [3], [3]);
                    %                     [w{mm},h{mm}, wtrace,htrace,ddiv{mm}]=nmf_testconvD(double(dveccr'),ncomp+1,1,hinit,winit, [3], [3]);
                    [c,w{mm},h{mm}, X1,X2, dhr, minXr, miXvalr, xcovpeakr, meanabsxcr,mhdr, p]=nmf_S41(double(dveccr),ncomp,1,winit,hinit, [3], [3],p);

                    hinit = [f*blinkmat./repmat(mean(blinkmat,2),1,size(blinkmat,2)); bg(mm)*sumw(ncomp)*ones(1, p.Nt)];             %random weights will be assigned to firts two and bg fixed
                    [c,w{mm},h{mm}, X1,X2, dht, minXt, miXvalt, xcovpeakt, meanabsxct,mhdt, p]=nmf_S41(double(dveccr),ncomp,1,winit,hinit, [3], [3],p);
                    qqq=[];
                else
                    [w{mm},h{mm}]=nmf(double(dveccr'),ncomp,1);
                end
                icapixNMF{mm} = reshape(w{mm},p.nx,p.ny,ncomp);
            end

        end

        p.path_data = path_data;
        p.path_res = path_res;

        if savethis == 1
            fprintf('saving data \n');
            if ~(strcmp(path_res, path_data)) %not identical
                mkdir ([path_res namedir]);
                cd ([path_res namedir]);
            end
            %             save (p.namedir)
            save ([namedir '_separMulti'],'p', 'X1','X2','dhr','minXr', 'miXvalr','dht','minXt', 'miXvalt', 'xcovpeakr', 'meanabsxcr','mhdr', 'xcovpeakt', 'meanabsxct','mhdt')
            writedata([],[],p,[namedir '_param'])
        end
    end
end

fprintf('\n')

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [B_trn, B_tst] = compressSpH(X, SpHparam)

% input:
%          X: n*d, n is the number of database samples
%          SpHparam:
%              SpHparam.nbits---encoding length
%              SpHparam.centers---spherical centers
%              SpHparam.radii---spherical radii
% output:
%          B_trn: compacted binary code of training samples
%          B_tst: compacted binary code of test samples

ntrain =SpHparam.ntrain;
xData = X;
centers = SpHparam.centers;
radii = SpHparam.radii;

% compute distances from centers
dData = distMat( xData , centers );

% compute binary codes for data points
th = repmat( radii' , size(dData , 1) , 1);
bData = zeros( size(dData) );
bData( dData <= th ) = 1;
bData = compactbit(bData);

B_trn = bData(1: ntrain, :);
B_tst = bData(ntrain+1:end, :);

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
% STD_PROPPLOT - Command line function to plot component cluster
%                  properties for a STUDY set.
%                  Displays mean cluster scalp map, ERP, ERSP;
%                  dipole model, spectrum, and ITC in one figure
%                  per cluster. Only measures computed during
%                  pre-clustering (by POP_PRECLUST or STD_PRECLUST)
%                  are plotted. Called by POP_CLUSTEDIT.
%                  Leaves the plotted grand mean cluster measures
%                  in STUDY.cluster for quick replotting.
% Usage:
%             >> [STUDY] = std_propplot(STUDY, ALLEEG, clusters);
% Inputs:
%   STUDY      - STUDY set including some or all EEG datasets in ALLEEG.
%   ALLEEG     - vector of EEG dataset structures including the datasets
%                in the STUDY. Yypically created using LOAD_ALLEEG.
%
% Optional inputs:
%   clusters   - [numeric vector | 'all'] -> cluster numbers to plot.
%                Else 'all' -> make plots for all clusters in the STUDY
%                {default: 'all'}.
% Outputs:
%   STUDY      - the input STUDY set structure modified with the plotted
%                cluster mean properties to allow quick replotting (unless
%                cluster means already existed in the STUDY).
% Example:
%              % Plot mean properties of Cluster 5 in one figure.
%              >> [STUDY] = std_propplot(STUDY,ALLEEG, 5);
%
% See also:  POP_CLUSTEDIT
%
% Authors:  Arnaud Delorme, Hilit Serby, Scott Makeig, SCCN/INC/UCSD, 2005

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

function STUDY = std_propplot(STUDY, ALLEEG,  varargin)

icadefs; % read EEGLAB defaults
warningon = 0;

if iscell(varargin{1}) % channel plotting
    chans = varargin{1};
    for k = 1: length(chans);
        figure
        orient tall
        set(gcf,'Color', BACKCOLOR);
        subplot(2,2,1),
        try,
            STUDY = std_erpplot(STUDY,ALLEEG, 'channels', chans(k), 'mode', 'together', 'plotmode', 'condensed' );
            erpAxisHangle = gca;
        catch
            axis off; text(0.5, 0.5, 'No ERP information', 'horizontalalignment', 'center');
            warningon = 1;
        end
        subplot(2,2,2),
        try,
            [STUDY] = std_erspplot(STUDY,ALLEEG, 'channels', chans(k), 'mode', 'together', 'plotmode', 'condensed' );
        catch,
            axis off; text(0.5, 0.5, 'No ERSP information', 'horizontalalignment', 'center');
            warningon = 1;
        end
        subplot(2,2,3),
        try,
            STUDY = std_specplot(STUDY,ALLEEG, 'channels', chans(k), 'mode', 'together', 'plotmode', 'condensed');
        catch
            axis off; text(0.5, 0.5, 'No spectral information', 'horizontalalignment', 'center');
            warningon = 1;
        end
        subplot(2,2,4),
        try,
            [STUDY] = std_itcplot(STUDY,ALLEEG, 'channels', chans(k), 'mode', 'together', 'plotmode', 'condensed' );
        catch,
            axis off; text(0.5, 0.5, 'No ITC information', 'horizontalalignment', 'center');
            warningon = 1;
        end

        %subplot('position', [0.77 0.16 0.15 0.28]),
        maintitle = ['Channel ''' chans{k} ''' average properties' ];
        a = textsc(maintitle, 'title');
        set(a, 'fontweight', 'bold');

        if warningon
            disp('Some properties could not be plotted. To plot these properties, first');
            disp('include them in pre-clustering. There, specify 0 dimensions if you do');
            disp('now want a property (scalp map, ERSP, etc...) to be included');
            disp('in the clustering procedure. See the clustering tutorial.');
        end

    end  % Finished all conditions
    return;
end;

% Set default values
cls = 1:length(STUDY.cluster); % plot all clusters in STUDY

if length(varargin) > 0
    if length(varargin) == 1, varargin{2} = varargin{1}; end; % backward compatibility
    if isnumeric(varargin{2})
        cls = varargin{2};
    elseif ischar(varargin{2}) && strcmpi(varargin{2}, 'all')
        cls = 1:length(STUDY.cluster);
    else
        error('cluster input should be either a vector of cluster indices or the keyword ''all''.');
    end
end

len = length(cls);
% Plot clusters mean properties
for k = 1: len
    if k == 1
        try
            % optional 'CreateCancelBtn', 'delete(gcbf); error(''USER ABORT'');',
            h_wait = waitbar(0,['Computing cluster properties ...'], 'Color', BACKEEGLABCOLOR,'position', [300, 200, 300, 48]);
        catch % for Matlab 5.3
            h_wait = waitbar(0,['Computing cluster properties ...'],'position', [300, 200, 300, 48]);
        end
    end
    warningon = 0;
    figure
    orient tall
    set(gcf,'Color', BACKCOLOR);
    subplot(2,3,1),
    try,
        STUDY = std_topoplot(STUDY,ALLEEG, 'clusters', cls(k), 'mode', 'together', 'figure', 'off');
    catch
        axis off; text(0.5, 0.5, 'No scalp map information', 'horizontalalignment', 'center');
        warningon = 1;
    end
    waitbar(k/(len*6),h_wait)
    subplot(2,3,2),
    try,
        STUDY = std_erpplot(STUDY,ALLEEG, 'clusters', cls(k), 'mode', 'together', 'plotmode', 'condensed' );
        erpAxisHangle = gca;
    catch
        axis off; text(0.5, 0.5, 'No ERP information', 'horizontalalignment', 'center');
        warningon = 1;
    end
    waitbar((k*2)/(len*6),h_wait)
    subplot(2,3,3),
    try,
        [STUDY] = std_erspplot(STUDY,ALLEEG, 'clusters', cls(k), 'mode', 'together', 'plotmode', 'condensed' );
    catch,
        axis off; text(0.5, 0.5, 'No ERSP information', 'horizontalalignment', 'center');
        warningon = 1;
    end
    waitbar((k*3)/(len*6),h_wait)
    axes('unit', 'normalized', 'position', [0.1 0.16 0.2 0.28]); %subplot(2,3,4),
    try,
        STUDY = std_dipplot(STUDY,ALLEEG, 'clusters', cls(k), 'mode', 'apart', 'figure', 'off'); set(gcf,'Color', BACKCOLOR);
    catch
        axis off; text(0.5, 0.5, 'No dipole information', 'horizontalalignment', 'center');
        warningon = 1;
    end
    waitbar((k*4)/(len*6),h_wait)
    subplot(2,3,5),
    try,
        STUDY = std_specplot(STUDY,ALLEEG, 'clusters', cls(k), 'mode', 'together', 'plotmode', 'condensed');
    catch
        axis off; text(0.5, 0.5, 'No spectral information', 'horizontalalignment', 'center');
        warningon = 1;
    end
    waitbar((k*5)/(len*6),h_wait)
    subplot(2,3,6),
    try,
        [STUDY] = std_itcplot(STUDY,ALLEEG, 'clusters', cls(k), 'mode', 'together', 'plotmode', 'condensed' );
    catch,
        axis off; text(0.5, 0.5, 'No ITC information', 'horizontalalignment', 'center');
        warningon = 1;
    end
    waitbar((k*6)/(len*6),h_wait);
    %subplot('position', [0.77 0.16 0.15 0.28]),
    maintitle = ['Cluster '''  STUDY.cluster(cls(k)).name ''' mean properties (' num2str(length(STUDY.cluster(cls(k)).comps)) ' comps).' ];
    a = textsc(maintitle, 'title');
    set(a, 'fontweight', 'bold');
    set(gcf,'name',maintitle);
    if warningon
        disp('Some properties could not be plotted. To plot these properties, first');
        disp('include them in pre-clustering. There, specify 0 dimensions if you do');
        disp('not want a property (scalp map, ERSP, etc...) to be included');
        disp('in the clustering procedure. See the clustering tutorial.');
    end

end  % Finished all conditions

if ishandle(erpAxisHangle) % make sure it is a valid graphics handle
    legend(erpAxisHangle,'off');

    set(erpAxisHangle,'YTickLabelMode','auto');
    set(erpAxisHangle,'YTickMode','auto');
end
delete(h_wait)

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function watchdog

% WATCHDOG will trigger an exit() if the controller disappears or if the allowed time elapsed
%
% To enable the watchdog you should do
%   watchdog(controllerid, time)
% and to disable it again
%   clear watchdog

% Copyright (C) 2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

warning('could not locate mex file');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function UNew = diffusionDirichlet3D(varargin);
% diffusionDirichlet3D: solve diffusion registraion in 3D with Dirichlet
%        boundary conditions
%
%
% author: Nathan D. Cahill
% email: nathan.cahill@rit.edu
% affiliation: Rochester Institute of Technology
% date: January 2014
% licence: GNU GPL v3
%
% Copyright Nathan D. Cahill
% Code available from https://github.com/tomdoel/npReg
%
%

% parse input arguments
[U,F,PixSize,M,N,P,RegularizerFactor] = parse_inputs(varargin{:});

% add displacement vectors to multiple of force field
FNew = F/RegularizerFactor;

% compute sine transform of new force field
FS = discreteSineTransform(FNew,M,N,P);

% construct images of coordinates scaled by pi/(N or M or P)
[a,b,c] = ndgrid(pi*(0:(M-1))/(M-1),pi*(0:(N-1))/(N-1),pi*(0:(P-1))/(P-1));

% construct LHS factor
LHSfactor = 2*cos(a) + 2*cos(b) + 2*cos(c) - 6;

% if gamma is zero, set origin term to 1, as DC term does not matter
LHSfactor(1,1,1) = 1;

% solve for FFT of U
US = cat(4,FS(:,:,:,1)./LHSfactor,FS(:,:,:,2)./LHSfactor,FS(:,:,:,3)./LHSfactor);

% perform inverse DST
UNew = discreteSineTransform(US,M,N,P);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FS = discreteSineTransform(F,M,N,P);
% compute discrete sine transform of 3-D vector field

% initialize resulting array
FS = F;

% first perform sine transform down columns
len = 2*M-2; ind = 1:M;
for p=1:P
    for n=1:N
        s = fft(FS(:,n,p,:),len,1);
        FS(:,n,p,:) = imag(s(ind,:,:,:));
    end
end
FS = sqrt(2/(M-1))*FS;

% next perform sine transform across rows
len = 2*N-2; ind = 1:N;
for p=1:P
    for m=1:M
        s = fft(FS(m,:,p,:),len,2);
        FS(m,:,p,:) = imag(s(:,ind,:,:));
    end
end
FS = sqrt(2/(N-1))*FS;

% finally perform sine transform across pages
len = 2*P-2; ind = 1:P;
for n=1:N
    for m=1:M
        s = fft(FS(m,n,:,:),len,3);
        FS(m,n,:,:) = imag(s(:,:,ind,:));
    end
end
FS = sqrt(2/(P-1))*FS;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U,F,PixSize,M,N,P,RegularizerFactor] = parse_inputs(varargin);

% get displacement field and check size
U = varargin{1};
F = varargin{2};
PixSize = varargin{4}(1:3);
M = varargin{5};
N = varargin{6};
P = varargin{7};
RegularizerFactor = varargin{10};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: tomdoel

repo: pulmonarytoolkit


text:
function str = optiMLVer()

v = ver('matlab');
if(v.Release(1) == '(')
    str = v.Release(2:end-1);
else
    str = v.Release;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% -   Minus.
%    X - Y subtracts matrix Y from X. X and Y must have compatible sizes. In
%    the simplest cases, they can be the same size or one can be a scalar.
%    Two inputs have compatible sizes if, for every dimension, the dimension
%    sizes of the inputs are either the same or one of them is 1.
%
%    C = MINUS(A,B) is called for the syntax 'A - B' when A or B is an
%    object.
%
%    Reference page in Doc Center
%       doc minus
%
%    Other functions named minus
%
%       calendarDuration/minus    fints/minus       sym/minus
%       cell/minus                gpuArray/minus    tall/minus
%       codistributed/minus       LagOp/minus       timeseries/minus
%       datetime/minus            mtree/minus       ts/minus
%       duration/minus
%

author: jmaih


text:
% Chapter 2 - Nonlinear Discrete Dynamical Systems.
% Program 2g - Computing the Lyapunov Exponents of the Henon map.
% Copyright Birkhauser 2013. Stephen Lynch.

itermax=500;
a=1.2;b=0.4;x=0;y=0;
vec1=[1;0];vec2=[0;1];
for i=1:itermax
    x1=1-a*x^2+y;y1=b*x;
    x=x1;y=y1;
    J=[-2*a*x 1;b 0];
    vec1=J*vec1;
    vec2=J*vec2;
    dotprod1=dot(vec1,vec1);
    dotprod2=dot(vec1,vec2);
    vec2=vec2-(dotprod2/dotprod1)*vec1;
    lengthv1=sqrt(dotprod1);
    area=vec1(1)*vec2(2)-vec1(2)*vec2(1);
    h1=log(lengthv1)/i;
    h2=log(area)/i-h1;
end
fprintf('h1= %12.10f\n',h1)
fprintf('h2= %12.10f\n',h2)

% End of Program 2g.

author: Sable

repo: mcbench-benchmarks


text:
function [X,r] = imstack2vectors(S,MASK)
%IMSTACK2VECTORS Extracts vectors from an image stack.
%   [X,R] = imstack2vectors(S,MASK) extracts vectors from S, which is an
%   M-by-N-by-n stack array of n registered images of size M-by-N each
%   (see Fig. 13.32 in DIPUM3). The extracted vectors are arranged as
%   the rows of matrix X. Input MASK is an M-by-N logical or numeric
%   image with nonzero values in the locations where elements of S are
%   to be used in forming X and 0s in locations to be ignored. The
%   number of row vectors in X is equal to the number of nonzero
%   elements of MASK. If MASK is omitted, all M*N locations are used in
%   forming X.  A simple way to obtain MASK interactively is to use
%   function roipoly. Finally, R is a column vector that contains the
%   linear indices of the locations of the vectors extracted from S.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Stack size.
[M,N,n] = size(S);

% Default.
if nargin == 1
   MASK = true(M,N);
else
   MASK = MASK ~= 0;
end

% Find the linear indices of the 1-valued elements in MASK. Each
% element of r identifies the location in the M-by-N array of the
% vector extracted from S.
r = find(MASK);

% The remaining code finds X.

% First reshape S into X by turning each set of n values along the
% third dimension of S so that it becomes a row of X. The order is
% from top to bottom along the first column, the second column, and
% so on.
Q = M*N;
X = reshape(S,Q,n);

% Now reshape MASK so that it corresponds to the right locations
% vertically along the elements of X.
MASK = reshape(MASK,Q,1);

% Keep the rows of X at locations where MASK is not 0.
X = X(MASK,:);

author: dipum

repo: dipum-toolbox

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


text:
classdef TestApproxPolyDP
    %TestApproxPolyDP

    methods (Static)
        function test_1
            curve = {[0 0], [1 1], [2 2], [3 3], [4 4]};
            approxCurve = cv.approxPolyDP(curve);
            assert(iscell(approxCurve));
            if ~isempty(approxCurve)
                assert(all(cellfun(@numel, approxCurve) == 2));
            end
        end

        function test_2
            curve = single([0 0; 10 0; 10 10; 5 4]);
            approxCurve = cv.approxPolyDP(curve, 'Epsilon',5, 'Closed',true);
            assert(isnumeric(approxCurve) && isfloat(approxCurve));
            if ~isempty(approxCurve)
                assert(ismatrix(approxCurve) && size(approxCurve,2)==2);
            end
        end

        function test_3
            curve = int32([0 0; 1 1; 2 2; 3 3; 4 4]);
            approxCurve = cv.approxPolyDP(curve);
            assert(isnumeric(approxCurve) && isinteger(approxCurve));
        end

        function test_error_argnum
            try
                cv.approxPolyDP();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function test_online_nmf(varargin)
% Test file for stochastic/online NMF.
%
% Created by H.Kasai on Apr. 17, 2018

    if nargin < 1
        clc;
        clear;
        close all;
        rng('default')

        m = 500;
        n = 100;
        V = rand(m,n);
        rank = 20;
        options = [];
        options.verbose = 1;
        max_epoch = 100;
        options.max_epoch = max_epoch;
        health_check_mode = false;
        calc_sol = 1;
        outlier_rho = 0.0;
        lambda = 1;
    else
        V = varargin{1};
        rank = varargin{2};
        options = varargin{3};
        max_epoch = options.max_epoch;
        calc_sol = 0;
        health_check_mode = true;
    end

    max_h_repeat    = 5;
    svrmu_inneriter = 5;

    % initialize factor matrices
    [x_init, ~] = generate_init_factors(V, rank, []);
    options.x_init = x_init;

    %% calc solution
    if calc_sol

        clear options_sol;
        options_sol.max_epoch = max_epoch*100;
        options_sol.x_init = x_init;
        options_sol.verbose = 0;
        options_sol.verbose = 1;

        if ~outlier_rho
            fprintf('Calculating f_opt by HALS ...\n');
            options_sol.alg = 'hals';
            [w_sol, infos_sol] = als_nmf(V, rank, options_sol);
        else
            fprintf('Calculating f_opt by R-NMF ...\n');
            options_sol.lambda = lambda;
            [w_sol, infos_sol] = robust_mu_nmf(V, rank, options_sol);
        end

        f_opt = infos_sol.cost(end);
        fprintf('Done.. f_opt: %.16e\n', f_opt);
    else
        f_opt = -Inf;
    end

    %% online_mu_nmf
    [w_onmf, infos_onmf] = online_mu_nmf(V, rank, options);

    %% acc_online_mu_nmf
    options.rep_mode = 'fix';
    options.w_repeat = 1;
    options.h_repeat = max_h_repeat;
    [w_onmf_acc, infos_onmf_acc] = acc_online_mu_nmf(V, rank, options);

    options.rep_mode = 'adaptive';
    options.w_repeat = 1;
    options.h_repeat = max_h_repeat;
    [w_onmf_acc, infos_onmf_acc] = acc_online_mu_nmf(V, rank, options);

    %% inf
    [w_inmf, infos_inmf] = incremental_mu_nmf(V, rank, options);

    options.online = 1; % online mode
    [w_online_inmf, infos_online_inmf] = incremental_mu_nmf(V, rank, options);

    %% asag_mu
    [w_asag_mu_nmf, infos_asag_mu_nmf] = asag_mu_nmf(V, rank, options);

    %% smu
    [w_smu_nmf, infos_smu_nmf] = smu_nmf(V, rank, options);

    options.accel = true;
    options.rep_mode = 'fix';
    options.w_repeat = 1;
    options.h_repeat = max_h_repeat;
    [w_smu_nmf_acc, infos_smu_nmf_acc] = smu_nmf(V, rank, options);
    options.accel = false;

    options.ls = false;
    options.rep_mode = 'adaptive';
    options.w_repeat = 1;
    options.h_repeat = max_h_repeat;
    [w_smu_nmf_acc, infos_smu_nmf_acc] = smu_nmf(V, rank, options);
    options.accel = 0;
    options.ls = 0;

    %% spg
    [w_spg_nmf, infos_spg_nmf] = spg_nmf(V, rank, options);

    options.accel = 1;
    options.h_repeat = max_h_repeat;
    [w_spg_nmf, infos_spg_nmf] = spg_nmf(V, rank, options);

    options.accel = 0;
    %options.h_repeat = max_h_repeat;
    options.W_sub_mode = 'Nesterov';
    [w_spg_nmf, infos_spg_nmf] = spg_nmf(V, rank, options);

    options.accel = 0;
    options.ls = 1;
    [w_spg_nmf, infos_spg_nmf] = spg_nmf(V, rank, options);

    %% svrmu
    options.x_init.R = zeros(size(V)); % enfore non-robust mode
    options.fast_calc = 0;
    options.repeat_inneriter = svrmu_inneriter;
    options.accel = 0;
    options.ls = 0;
    options.ratio = 1; % 1: original 1<: adaptive
    options.robust = false;
    options.max_epoch = floor(max_epoch / (options.repeat_inneriter + 1));
    [w_svrmu_nmf, infos_svrmu_nmf] = svrmu_nmf(V, rank, options);

    options.accel = 1;
    options.rep_mode = 'fix';
    [w_svrmu_acc_nmf, infos_svrmu_acc_nmf] = svrmu_nmf(V, rank, options);

    options.rep_mode = 'adaptive';
    options.accel = 0;
    [w_svrmu_acc_nmf, infos_svrmu_acc_nmf] = svrmu_nmf(V, rank, options);

    options.accel = 1;
    [w_svrmu_acc_nmf, infos_svrmu_acc_nmf] = svrmu_nmf(V, rank, options);

    options.fast_calc = 0;
    options.accel = 1;
    options.rep_mode = 'adaptive';
    options.ratio = 0.3; % 1: original 1<: adaptive
    options.robust = false;
    [w_svrmu_acc_adaptive_nmf, infos_svrmu_acc_adaptive_nmf] = svrmu_nmf(V, rank, options);

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [w,s]=fram2wav(x,tt,mode)
%FRAM2WAV  converts frame values to a continuous waveform [W]=(X,TT,MODE)
%  Inputs:
%          x(nf,p)      is the input signal: one row per frame
%	       tt(nf,3)     specifies the frames. Each row has the form [start_sample end_sample flag]
%                       where flag = 1 for the start of a new spurt.
%                       If tt(:,3) is omitted, a new spurt will be started whenever there is a gap
%                       of more than one between the end of one frame and the beginning or the next.
%                       A new spurt is automatically started if x() = NaN.
%          mode         consists of one or more of the following letters:
%                          z for zero-order hold interpolation (i.e. constant within each frame)
%                          l for linear interpolation within each spurt [default]
%
% Outputs:
%          w(n,p)       contains the interpolated waveforms. Their length is n = tt(nf,2)
%          s(ns,2)      gives the starting and ending sample numbers of each spurt (excluding NaN spurts)
%
%    This routine converts frame-based values to continuous waveforms by performing
%    a chosen method of interpolation. Interpolation is restarted at the beginning of each spurt.

%    Bugs/Suggestions
%      (1)   Additional mode option for cubic interpolation
%      (2)   Additional mode option for interpolation in log domain
%      (3)   Additional mode option for x values being
%            frame averages rather than mid-frame values.

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: fram2wav.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3
    mode='l';
end
[nf,m]=size(x);
n=round(tt(end,2));
w=repmat(NaN,n,m);
nt=size(tt,2);
ix1=ceil(tt(:,1)); % start of frame sample
ix2=floor(tt(:,2)); % end of frame sample

% determine the start and end of spurts

if nt>2
    ty=tt(:,3)>0;   % frame type set by user
else
    ty=zeros(nf,1);
    ty(2:end)=ix1(2:end)>ix2(1:end-1)+1;    % new spurt whenever a gap
end
ty(1)=1;           % first frame always starts a spurt
ty(isnan(x))=1;    % NaN always ends previous spurt
ty(1+find(isnan(x(1:end-1))))=1; % NaN always forces a new spurt
ty=double(ty);
ty(1:end-1)=ty(1:end-1)+2*ty(2:end);
ty(end)=ty(end)+2;   % last frame always ends a spurtw=repmat(NaN,n,m);  % initialize output to all NaN
nx=ix2-ix1+1;

if any(mode=='z')   % zero-order hold
    for i=1:nf
        if nx(i)
            w(ix1(i):ix2(i),:)=repmat(x(i,:),nx(i),1);
        end
    end
else   % linear interpolation is the default
    ttm=(tt(:,1)+tt(:,2))/2;    % mid point of frame
    ixm=floor(ttm); % end of first half of frame
    for i=1:nf
        if i==176
            i
        end
        if nx(i)
            tyi=ty(i);
            if tyi==3    % use a zero order hold
                w(ix1(i):ix2(i),:)=repmat(x(i,:),nx(i),1);
            else
                nxm=ixm(i)-ix1(i)+1;
                if nxm
                    if tyi==1
                        grad=(x(i+1,:)-x(i,:))/(ttm(i+1)-ttm(i));
                    else
                        grad=(x(i,:)-x(i-1,:))/(ttm(i)-ttm(i-1));
                    end
                    w(ix1(i):ixm(i),:)=repmat(x(i,:),nxm,1)+((ix1(i):ixm(i))'-ttm(i))*grad;
                end
                if nx(i)>nxm
                    if tyi==2
                        grad=(x(i,:)-x(i-1,:))/(ttm(i)-ttm(i-1));
                    else
                        grad=(x(i+1,:)-x(i,:))/(ttm(i+1)-ttm(i));
                    end
                    w(ixm(i)+1:ix2(i),:)=repmat(x(i,:),ix2(i)-ixm(i),1)+((ixm(i)+1:ix2(i))'-ttm(i))*grad;
                end
            end
        end
    end
end

% now sort out the start and end spurt positions

ty(isnan(x))=0;    % Don't count NaN spurts
s=repmat(ix1(bitand(ty,1)>0),1,2);
s(:,2)=ix2(bitand(ty,2)>0);
if ~nargout
    tw=(1:n)';
    for i=size(s,1):-1:2
        j=s(i,1);   % start of new spurt
        tw=[tw(1:j-1); tw(j); tw(j:end)];
        w=[w(1:j-1); NaN; w(j:end)];        % insert a NaN to force a plotting break
    end
    plot(tt(:,1:2)',repmat(x(:)',2,1),'r-+',tw,w,'b-');
end

author: jtkim-kaist

repo: VAD

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


text:
function [Vx,Vy,reliab]=opticalFlow( I1, I2, varargin )
% Coarse-to-fine optical flow using Lucas&Kanade or Horn&Schunck.
%
% Implemented 'type' of optical flow estimation:
%  LK: http://en.wikipedia.org/wiki/Lucas-Kanade_method
%  HS: http://en.wikipedia.org/wiki/Horn-Schunck_method
% LK is a local, fast method (the implementation is fully vectorized).
% HS is a global, slower method (an SSE implementation is provided).
%
% Common parameters for LK and HS: 'smooth' determines smoothing prior to
% flow computation and can make flow estimation more robust. 'resample' can
% be used to downsample an image for faster but lower quality results, e.g.
% resample=.5 makes flow computation about 4x faster. LK: 'radius' controls
% integration window size (and smoothness of flow). HS: 'alpha' controls
% tradeoff between data and smoothness term (and smoothness of flow) and
% 'nIter' determines number of gradient decent steps.
%
% USAGE
%  [Vx,Vy,reliab] = opticalFlow( I1, I2, pFlow )
%
% INPUTS
%  I1, I2   - input images to calculate flow between
%  pFlow    - parameters (struct or name/value pairs)
%   .type       - ['LK'] may be either 'LK' or 'HS'
%   .smooth     - [1] smoothing radius for triangle filter (may be 0)
%   .resample   - [1] resampling amount (must be a power of 2)
%   .radius     - [5] integration radius for weighted window [LK only]
%   .alpha      - [1] smoothness constraint [HS only]
%   .nIter      - [250] number of iterations [HS only]
%
% OUTPUTS
%  Vx, Vy   - x,y components of flow  [Vx>0->right, Vy>0->down]
%  reliab   - reliability of flow in given window [LK only]
%
% EXAMPLE - compute LK flow on test images
%  load opticalFlowTest;
%  [Vx,Vy]=opticalFlow(I1,I2,'smooth',1,'radius',10,'type','LK');
%  figure(1); im(I1); figure(2); im(I2);
%  figure(3); im([Vx Vy]); colormap jet;
%
% EXAMPLE - rectify I1 to I2 using computed flow
%  load opticalFlowTest;
%  [Vx,Vy]=opticalFlow(I1,I2,'smooth',1,'radius',10,'type','LK');
%  I1=imtransform2(I1,[],'vs',-Vx,'us',-Vy,'pad','replicate');
%  figure(1); im(I1); figure(2); im(I2);
%
% EXAMPLE - compare LK and HS flow
%  load opticalFlowTest;
%  prm={'smooth',1,'radius',10,'alpha',20,'nIter',200,'type'};
%  tic, [Vx1,Vy1]=opticalFlow(I1,I2,prm{:},'LK'); toc
%  tic, [Vx2,Vy2]=opticalFlow(I1,I2,prm{:},'HS'); toc
%  figure(1); im([Vx1 Vy1; Vx2 Vy2]); colormap jet;
%
% See also convTri, imtransform2
%
% Piotr's Image&Video Toolbox      Version 3.02
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% get default parameters and do error checking
dfs={'type','LK','smooth',1,'resample',1,'radius',5,'alpha',1,'nIter',250};
[type,smooth,resample,radius,alpha,nIter]=getPrmDflt(varargin,dfs,1);
assert(any(strcmp(type,{'LK','HS'}))); useLk=strcmp(type,'LK');
if( ~ismatrix(I1) || ~ismatrix(I2) || any(size(I1)~=size(I2)) )
  error('Input images must be 2D and have same dimensions.'); end

% run optical flow in coarse to fine fashion
if(~isa(I1,'single')), I1=single(I1); I2=single(I2); end
[h,w]=size(I1); nScales=floor(log2(min(h,w)))-2;
for s=1:nScales + round(log2(resample))
  % get current scale and I1s and I2s at given scale
  scale=2^(nScales-s); h1=round(h/scale); w1=round(w/scale);
  if( scale==1 ), I1s=I1; I2s=I2; else
    I1s=imResample(I1,[h1 w1]); I2s=imResample(I2,[h1 w1]); end
  % initialize Vx,Vy or upsample from previous scale
  if(s==1), Vx=zeros(h1,w1,'single'); Vy=Vx; else r=sqrt(h1*w1/numel(Vx));
    Vx=imResample(Vx,[h1 w1])*r; Vy=imResample(Vy,[h1 w1])*r; end
  % transform I1s according to current estimate of Vx and Vy
  if(s), I1s=imtransform2(I1s,[],'pad','replciate','vs',-Vx,'us',-Vy); end
  % smooth images
  I1s=convTri(I1s,smooth); I2s=convTri(I2s,smooth);
  % run optical flow on current scale
  if( useLk ), [Vx1,Vy1,reliab]=opticalFlowLk(I1s,I2s,radius);
  else [Vx1,Vy1]=opticalFlowHs(I1s,I2s,alpha,nIter); reliab=[]; end
  Vx=Vx+Vx1; Vy=Vy+Vy1;
end
if(s~=nScales), r=sqrt(h*w/numel(Vx));
  Vx=imResample(Vx,[h w])*r; Vy=imResample(Vy,[h w])*r; end

end

function [Vx,Vy,reliab] = opticalFlowLk( I1, I2, radius  )
% Compute elements of A'A and also of A'b
radius=min(radius,floor(min(size(I1,1),size(I1,2))/2)-1);
[Ix,Iy]=gradient2(I1); It=I2-I1; AAxy=convTri(Ix.*Iy,radius);
AAxx=convTri(Ix.^2,radius)+1e-5; ABxt=convTri(-Ix.*It,radius);
AAyy=convTri(Iy.^2,radius)+1e-5; AByt=convTri(-Iy.*It,radius);
% Find determinant and trace of A'A
AAdet=AAxx.*AAyy-AAxy.^2; AAdeti=1./AAdet; AAtr=AAxx+AAyy;
% Compute components of velocity vectors (A'A)^-1 * A'b
Vx = AAdeti .* ( AAyy.*ABxt - AAxy.*AByt);
Vy = AAdeti .* (-AAxy.*ABxt + AAxx.*AByt);
% Check for ill conditioned second moment matrices
reliab = 0.5*AAtr - 0.5*sqrt(AAtr.^2-4*AAdet);
end

function [Vx,Vy] = opticalFlowHs( I1, I2, alpha, nIter )
% compute derivatives (averaging over 2x2 neighborhoods)
A00=shift(I1,0,0); A10=shift(I1,1,0);
A01=shift(I1,0,1); A11=shift(I1,1,1);
B00=shift(I2,0,0); B10=shift(I2,1,0);
B01=shift(I2,0,1); B11=shift(I2,1,1);
Ex=0.25*((A01+B01+A11+B11)-(A00+B00+A10+B10));
Ey=0.25*((A10+B10+A11+B11)-(A00+B00+A01+B01));
Et=0.25*((B00+B10+B01+B11)-(A00+A10+A01+A11));
Ex([1 end],:)=0; Ex(:,[1 end])=0;
Ey([1 end],:)=0; Ey(:,[1 end])=0;
Et([1 end],:)=0; Et(:,[1 end])=0;
Z=1./(alpha*alpha + Ex.*Ex + Ey.*Ey);
% iterate updating Ux and Vx in each iter
if( 1 )
  [Vx,Vy]=opticalFlowHsMex(Ex,Ey,Et,Z,nIter);
  Vx=Vx(2:end-1,2:end-1); Vy=Vy(2:end-1,2:end-1);
else
  Vx=zeros(size(I1),'single'); Vy=Vx;
  for i = 1:nIter
    Mx=.25*(shift(Vx,-1,0)+shift(Vx,1,0)+shift(Vx,0,-1)+shift(Vx,0,1));
    My=.25*(shift(Vy,-1,0)+shift(Vy,1,0)+shift(Vy,0,-1)+shift(Vy,0,1));
    m=(Ex.*Mx+Ey.*My+Et).*Z; Vx=Mx-Ex.*m; Vy=My-Ey.*m;
    Vx=Vx(2:end-1,2:end-1); Vy=Vy(2:end-1,2:end-1);
  end
end
end

function J = shift( I, y, x )
% shift I by -1<=x,y<=1 pixels
[h,w]=size(I); J=zeros(h+2,w+2,'single');
J(2-y:end-1-y,2-x:end-1-x)=I;
end

author: zouchuhang

repo: LayoutNet


text:
function B = rowNormalize(A)
% DESCRIPTION: normalize each row of the matrix
A_n_r = sqrt(sum(A.^2,2,'omitnan'));
B = A./repmat(A_n_r,1,size(A,2));
end

author: goGPS-Project


text:
 function [err, sn, kernel] = nufft1_error(om, N1, J1, K1, kernel, sn)
%function [err, sn, kernel] = nufft1_error(om, N1, J1, K1, kernel, sn)
%
% Compute worst-case error for each input frequency for 1D NUFFT
% using specified (inline) `ad hoc' kernel function (e.g., gaussian).
% This is worst-case for a unit-norm signal of length N1.
% in
%	om	[M,1]	digital frequency omega in radians
%	N1		signal length
%	J1		# of neighbors used per frequency location
%	K1		FFT size (should be > N1)
%	kernel		inline kernel function, args (k,J)
%				(or choose from built-in examples - see code!)
%	sn		optional scaling factors (otherwise do-no-harm)
% out
%	err	[M,1]	worst-case error over unit-norm signals
%	sn	[N,1]	scaling factors
%
% examples for kernel:
% linear:		inline('(1 - abs(k/(J/2))) .* (abs(k) < J/2)', 'k', 'J')
% truncated diric:	inline('sinc(k) .* (abs(k) < J/2)', 'k', 'J')
%
% Copyright 2001-12-7, Jeff Fessler, The University of Michigan

% if no arguments, give an example, comparing triangular kernel to min-max
if nargin < 4
	help(mfilename)
	N = 2^7; K = 2*N; gam = 2*pi/K;
	Jlist = [2:10]';
	om = gam * linspace(0,1,101);

	err.linear	= zeros(size(Jlist));
	err.minmaxu	= zeros(size(Jlist));
	err.minmax2	= zeros(size(Jlist));
	err.minmaxo	= zeros(size(Jlist));
	err.minmaxk	= zeros(size(Jlist));	% kaiser sn's
	err.gauss_zn	= zeros(size(Jlist));
	err.gauss_ft	= zeros(size(Jlist));
	err.kaiser	= zeros(size(Jlist));

	for ii=1:length(Jlist)
		J = Jlist(ii);
		printf('J=%d', J)
		err.linear(ii) = ...
			max(nufft1_error(om, N, J, K, 'linear'));
		err.minmaxu(ii) = ...
			max(nufft1_error(om, N, J, K, 'minmax,uniform'));
		err.minmax2(ii) = ...
			max(nufft1_error(om, N, J, K, 'minmax,best,L=2'));
		err.minmaxo(ii) = ...
			max(nufft1_error(om, N, J, K, 'minmax,best'));
		err.gauss_zn(ii) = ...
			max(nufft1_error(om, N, J, K, 'gauss'));
		err.gauss_ft(ii) = ...
			max(nufft1_error(om, N, J, K, 'gauss', 'ft'));
		[tmp, sn] = nufft1_error(om, N, J, K, 'kaiser', 'ft');
		err.kaiser(ii) = max(tmp);
		err.minmaxk(ii) = ...
			max(nufft1_err_mm(om, N, J, K, 'qr', sn));
	end
	clf, semilogy(Jlist, err.linear, 'g-x', ...
		Jlist, err.minmaxu, 'c-+', ...
		Jlist, err.gauss_zn, 'b-*', ...
		Jlist, err.gauss_ft, 'b-o', ...
		Jlist, err.minmax2, 'r-^', ...
		Jlist, err.kaiser, 'm->', ...
		Jlist, err.minmaxo, 'y-<', ...
		Jlist, err.minmaxk, 'w-o'), axis tight
	xlabel J, ylabel 'worst-case error'
	legend('linear', 'min-max, uniform', ...
		'gaussian (zn)', 'gaussian (FT)', ...
		'min-max, best L=2', 'kaiser', 'min-max, optimized', ...
		'min-max, kaiser s', 3)
	clear err
return
end

%
% kernel selection
%
if ischar(kernel)

	% min-max interpolators (uniform, best, etc., see nufft1_err_mm.m)
	if strncmp(kernel, 'minmax,', 7)
		type = kernel(8:end);	% uniform or best or best,L=2 etc.
		[err, sn] = nufft1_err_mm(om, N1, J1, K1, 'qr', type);
		return

	% cos^3-tapered dirichlet
	elseif streq(kernel, 'cos3diric')
		kernel = 'diric(2*pi*k/J, J) .* cos((2*pi*k/J)/2).^3';

	% Dirichlet (truncated)
	elseif streq(kernel, 'diric')
		kernel = 'nufft_diric(k,%d,%d,1) .* (abs(k) < J/2)';
		kernel = sprintf(kernel, N1, N1);

	% gaussian (truncated) with previously numerically-optimized width
	elseif streq(kernel, 'gauss')
		if isvar('sn') & ischar(sn) & streq(sn, 'ft')
			stype = 'ft';
		else
			stype = 'zn';
		end
		[dummy, kernel, kernel_ft] = nufft_best_gauss(J1, K1/N1, stype);
		kernel_ft = inline(kernel_ft, 't');

	% kaiser-bessel with previously numerically-optimized shape
	elseif streq(kernel, 'kaiser')
		[kernel, kb_a, kb_m] = kaiser_bessel('string', J1, 'best', 0, K1/N1);
		kernel_ft = kaiser_bessel_ft('inline', J1, kb_a, kb_m, 1);

	% linear interpolation via triangular function (possibly "wide"!)
	elseif streq(kernel, 'linear')
		kernel = '(1 - abs(k/(J/2))) .* (abs(k) < J/2)';

	else
		error(sprintf('unknown kernel "%s"', kernel))
	end

	kernel = inline(kernel, 'k', 'J');

elseif ~streq('inline', class(kernel))
	error 'need inline kernel'
end

gam = 2*pi/K1;

if 0
	%	plot interpolator
	k = linspace(-J/2-1,J/2+1,101);
	clf, subplot(221), plot(k, kernel(k, J))
	xlabel k, ylabel kernel(k), axis tight, grid
end

%
% Compute scaling factors using the "do no harm" strategy.
% This may not be optimal; analytical FT could also be reasonable.
%
if ~isvar('sn')
	sn = 1 ./ nufft_interp_zn(0, N1, J1, K1, kernel);	% [N]
%	sn = 1 ./ mean(nufft_interp_zn([0 1/2], N1, J1, K1, kernel), 2); % alt

elseif isa(sn, 'inline')
	n = [0:(N1-1)]'-(N1-1)/2;
	sn = 1 ./ sn(n/K1);		% [N]

% trick to use Gaussian FT scaling factors
elseif ischar(sn) & streq(sn, 'ft') & isvar('kernel_ft')
	n = [0:(N1-1)]'-(N1-1)/2;
	sn = 1 ./ kernel_ft(n/K1);		% [N]
	if 0 & J1 > 2
		sn_zn = 1 ./ nufft_interp_zn(0, N1, J1, K1, kernel);	% [N]
		clf, plot(n, [sn reale(sn_zn)])
		keyboard
	end

% trick to use Kaiser-Bessel FT scaling factors
elseif isstruct(sn) & streq(sn.type, 'kaiser')
	n = [0:(N1-1)]'-(N1-1)/2;
	sn = 1 ./ kaiser_bessel_ft(n/K1, J1, sn.alpha, sn.m, 1);

else
	error 'unsupport scaling factors type'
end

%
% interpolator worst-case error for each frequency (scaled by 1/sqrt(N)),
% from equations (46)-(47) in Fessler&Sutton NUFFT paper, T-SP, Feb. 2003
%
zn = nufft_interp_zn(om/gam, N1, J1, K1, kernel);		% [N,M]
err = sqrt(mean(abs(spdiag(sn) * zn - 1).^2, 1))';	% [M]

author: thomaskuestner


text:
function a = hilbert_fun ( m, n )

%*****************************************************************************80
%
%% HILBERT_FUN assembles a copy of the Hilbert matrix using SPMD workers.
%
%  Discussion:
%
%    The matrix is assembled in blocks.  Each SPMD worker builds one
%    block, and the client worker assembles the pieces into a single
%    array.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix
%    to be generated.
%
%    Output, real A(M,N), the Hilbert matrix.
%
  spmd
    ilo = 1;
    ihi = m;
    jlo = floor ( ( ( labindex ( ) - 1 ) * n ) / numlabs ( ) ) + 1;
    jhi = floor ( (   labindex ( )       * n ) / numlabs ( ) );
    a_block = hilbert_block ( ilo, jlo, ihi, jhi );
  end
%
%  The client collects the stripes into a single matrix.
%
  a = [ a_block{:} ];

  return
end
function a = hilbert_block ( ilo, jlo, ihi, jhi )

%*****************************************************************************80
%
%% HILBERT_BLOCK returns a block of the Hilbert matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ILO, JLO, IHI, JHI, the row and column of the first
%    and last entries in the block.
%
%    Output, real A(1:IHI+1-ILO,1:JHI+1-JLO), the block of the Hilbert matrix.
%
  a = zeros ( ihi + 1 - ilo, jhi + 1 - jlo );

  i = ilo:ihi;

  for j = jlo : jhi
    a(i+1-ilo,j+1-jlo) = i / j;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [baseMVA, bus, gen, branch, areas, gencost] = case57
%CASE57    Power flow data for IEEE 57 bus test case.
%   Please see 'help caseformat' for details on the case file format.
%   This data was converted from IEEE Common Data Format
%   (ieee57cdf.txt) on 20-Sep-2004 by cdf2matp, rev. 1.11
%   See end of file for warnings generated during conversion.
%
%   Converted from IEEE CDF file from:
%       http://www.ee.washington.edu/research/pstca/
%
%   Manually modified Qmax, Qmin on generator 1 to 200, -140, respectively.
%
%  08/25/93 UW ARCHIVE           100.0  1961 W IEEE 57 Bus Test Case

%   MATPOWER
%   $Id: case57.m,v 1.5 2004/09/21 01:47:48 ray Exp $

%%-----  Power Flow Data  -----%%
%% system MVA base
baseMVA = 100;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
bus = [
	1	3	55	17	0	0	1	1.04	0	0	1	1.06	0.94;
	2	2	3	88	0	0	1	1.01	-1.18	0	1	1.06	0.94;
	3	2	41	21	0	0	1	0.985	-5.97	0	1	1.06	0.94;
	4	1	0	0	0	0	1	0.981	-7.32	0	1	1.06	0.94;
	5	1	13	4	0	0	1	0.976	-8.52	0	1	1.06	0.94;
	6	2	75	2	0	0	1	0.98	-8.65	0	1	1.06	0.94;
	7	1	0	0	0	0	1	0.984	-7.58	0	1	1.06	0.94;
	8	2	150	22	0	0	1	1.005	-4.45	0	1	1.06	0.94;
	9	2	121	26	0	0	1	0.98	-9.56	0	1	1.06	0.94;
	10	1	5	2	0	0	1	0.986	-11.43	0	1	1.06	0.94;
	11	1	0	0	0	0	1	0.974	-10.17	0	1	1.06	0.94;
	12	2	377	24	0	0	1	1.015	-10.46	0	1	1.06	0.94;
	13	1	18	2.3	0	0	1	0.979	-9.79	0	1	1.06	0.94;
	14	1	10.5	5.3	0	0	1	0.97	-9.33	0	1	1.06	0.94;
	15	1	22	5	0	0	1	0.988	-7.18	0	1	1.06	0.94;
	16	1	43	3	0	0	1	1.013	-8.85	0	1	1.06	0.94;
	17	1	42	8	0	0	1	1.017	-5.39	0	1	1.06	0.94;
	18	1	27.2	9.8	0	10	1	1.001	-11.71	0	1	1.06	0.94;
	19	1	3.3	0.6	0	0	1	0.97	-13.2	0	1	1.06	0.94;
	20	1	2.3	1	0	0	1	0.964	-13.41	0	1	1.06	0.94;
	21	1	0	0	0	0	1	1.008	-12.89	0	1	1.06	0.94;
	22	1	0	0	0	0	1	1.01	-12.84	0	1	1.06	0.94;
	23	1	6.3	2.1	0	0	1	1.008	-12.91	0	1	1.06	0.94;
	24	1	0	0	0	0	1	0.999	-13.25	0	1	1.06	0.94;
	25	1	6.3	3.2	0	5.9	1	0.982	-18.13	0	1	1.06	0.94;
	26	1	0	0	0	0	1	0.959	-12.95	0	1	1.06	0.94;
	27	1	9.3	0.5	0	0	1	0.982	-11.48	0	1	1.06	0.94;
	28	1	4.6	2.3	0	0	1	0.997	-10.45	0	1	1.06	0.94;
	29	1	17	2.6	0	0	1	1.01	-9.75	0	1	1.06	0.94;
	30	1	3.6	1.8	0	0	1	0.962	-18.68	0	1	1.06	0.94;
	31	1	5.8	2.9	0	0	1	0.936	-19.34	0	1	1.06	0.94;
	32	1	1.6	0.8	0	0	1	0.949	-18.46	0	1	1.06	0.94;
	33	1	3.8	1.9	0	0	1	0.947	-18.5	0	1	1.06	0.94;
	34	1	0	0	0	0	1	0.959	-14.1	0	1	1.06	0.94;
	35	1	6	3	0	0	1	0.966	-13.86	0	1	1.06	0.94;
	36	1	0	0	0	0	1	0.976	-13.59	0	1	1.06	0.94;
	37	1	0	0	0	0	1	0.985	-13.41	0	1	1.06	0.94;
	38	1	14	7	0	0	1	1.013	-12.71	0	1	1.06	0.94;
	39	1	0	0	0	0	1	0.983	-13.46	0	1	1.06	0.94;
	40	1	0	0	0	0	1	0.973	-13.62	0	1	1.06	0.94;
	41	1	6.3	3	0	0	1	0.996	-14.05	0	1	1.06	0.94;
	42	1	7.1	4.4	0	0	1	0.966	-15.5	0	1	1.06	0.94;
	43	1	2	1	0	0	1	1.01	-11.33	0	1	1.06	0.94;
	44	1	12	1.8	0	0	1	1.017	-11.86	0	1	1.06	0.94;
	45	1	0	0	0	0	1	1.036	-9.25	0	1	1.06	0.94;
	46	1	0	0	0	0	1	1.05	-11.89	0	1	1.06	0.94;
	47	1	29.7	11.6	0	0	1	1.033	-12.49	0	1	1.06	0.94;
	48	1	0	0	0	0	1	1.027	-12.59	0	1	1.06	0.94;
	49	1	18	8.5	0	0	1	1.036	-12.92	0	1	1.06	0.94;
	50	1	21	10.5	0	0	1	1.023	-13.39	0	1	1.06	0.94;
	51	1	18	5.3	0	0	1	1.052	-12.52	0	1	1.06	0.94;
	52	1	4.9	2.2	0	0	1	0.98	-11.47	0	1	1.06	0.94;
	53	1	20	10	0	6.3	1	0.971	-12.23	0	1	1.06	0.94;
	54	1	4.1	1.4	0	0	1	0.996	-11.69	0	1	1.06	0.94;
	55	1	6.8	3.4	0	0	1	1.031	-10.78	0	1	1.06	0.94;
	56	1	7.6	2.2	0	0	1	0.968	-16.04	0	1	1.06	0.94;
	57	1	6.7	2	0	0	1	0.965	-16.56	0	1	1.06	0.94;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin
gen = [
	1	128.9	-16.1	200	-140	1.04	100	1	575.88	0;
	2	0	-0.8	50	-17	1.01	100	1	100	0;
	3	40	-1	60	-10	0.985	100	1	140	0;
	6	0	0.8	25	-8	0.98	100	1	100	0;
	8	450	62.1	200	-140	1.005	100	1	550	0;
	9	0	2.2	9	-3	0.98	100	1	100	0;
	12	310	128.5	155	-150	1.015	100	1	410	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status
branch = [
	1	2	0.0083	0.028	0.129	9900	0	0	0	0	1;
	2	3	0.0298	0.085	0.0818	9900	0	0	0	0	1;
	3	4	0.0112	0.0366	0.038	9900	0	0	0	0	1;
	4	5	0.0625	0.132	0.0258	9900	0	0	0	0	1;
	4	6	0.043	0.148	0.0348	9900	0	0	0	0	1;
	6	7	0.02	0.102	0.0276	9900	0	0	0	0	1;
	6	8	0.0339	0.173	0.047	9900	0	0	0	0	1;
	8	9	0.0099	0.0505	0.0548	9900	0	0	0	0	1;
	9	10	0.0369	0.1679	0.044	9900	0	0	0	0	1;
	9	11	0.0258	0.0848	0.0218	9900	0	0	0	0	1;
	9	12	0.0648	0.295	0.0772	9900	0	0	0	0	1;
	9	13	0.0481	0.158	0.0406	9900	0	0	0	0	1;
	13	14	0.0132	0.0434	0.011	9900	0	0	0	0	1;
	13	15	0.0269	0.0869	0.023	9900	0	0	0	0	1;
	1	15	0.0178	0.091	0.0988	9900	0	0	0	0	1;
	1	16	0.0454	0.206	0.0546	9900	0	0	0	0	1;
	1	17	0.0238	0.108	0.0286	9900	0	0	0	0	1;
	3	15	0.0162	0.053	0.0544	9900	0	0	0	0	1;
	4	18	0	0.555	0	9900	0	0	0.97	0	1;
	4	18	0	0.43	0	9900	0	0	0.978	0	1;
	5	6	0.0302	0.0641	0.0124	9900	0	0	0	0	1;
	7	8	0.0139	0.0712	0.0194	9900	0	0	0	0	1;
	10	12	0.0277	0.1262	0.0328	9900	0	0	0	0	1;
	11	13	0.0223	0.0732	0.0188	9900	0	0	0	0	1;
	12	13	0.0178	0.058	0.0604	9900	0	0	0	0	1;
	12	16	0.018	0.0813	0.0216	9900	0	0	0	0	1;
	12	17	0.0397	0.179	0.0476	9900	0	0	0	0	1;
	14	15	0.0171	0.0547	0.0148	9900	0	0	0	0	1;
	18	19	0.461	0.685	0	9900	0	0	0	0	1;
	19	20	0.283	0.434	0	9900	0	0	0	0	1;
	21	20	0	0.7767	0	9900	0	0	1.043	0	1;
	21	22	0.0736	0.117	0	9900	0	0	0	0	1;
	22	23	0.0099	0.0152	0	9900	0	0	0	0	1;
	23	24	0.166	0.256	0.0084	9900	0	0	0	0	1;
	24	25	0	1.182	0	9900	0	0	1	0	1;
	24	25	0	1.23	0	9900	0	0	1	0	1;
	24	26	0	0.0473	0	9900	0	0	1.043	0	1;
	26	27	0.165	0.254	0	9900	0	0	0	0	1;
	27	28	0.0618	0.0954	0	9900	0	0	0	0	1;
	28	29	0.0418	0.0587	0	9900	0	0	0	0	1;
	7	29	0	0.0648	0	9900	0	0	0.967	0	1;
	25	30	0.135	0.202	0	9900	0	0	0	0	1;
	30	31	0.326	0.497	0	9900	0	0	0	0	1;
	31	32	0.507	0.755	0	9900	0	0	0	0	1;
	32	33	0.0392	0.036	0	9900	0	0	0	0	1;
	34	32	0	0.953	0	9900	0	0	0.975	0	1;
	34	35	0.052	0.078	0.0032	9900	0	0	0	0	1;
	35	36	0.043	0.0537	0.0016	9900	0	0	0	0	1;
	36	37	0.029	0.0366	0	9900	0	0	0	0	1;
	37	38	0.0651	0.1009	0.002	9900	0	0	0	0	1;
	37	39	0.0239	0.0379	0	9900	0	0	0	0	1;
	36	40	0.03	0.0466	0	9900	0	0	0	0	1;
	22	38	0.0192	0.0295	0	9900	0	0	0	0	1;
	11	41	0	0.749	0	9900	0	0	0.955	0	1;
	41	42	0.207	0.352	0	9900	0	0	0	0	1;
	41	43	0	0.412	0	9900	0	0	0	0	1;
	38	44	0.0289	0.0585	0.002	9900	0	0	0	0	1;
	15	45	0	0.1042	0	9900	0	0	0.955	0	1;
	14	46	0	0.0735	0	9900	0	0	0.9	0	1;
	46	47	0.023	0.068	0.0032	9900	0	0	0	0	1;
	47	48	0.0182	0.0233	0	9900	0	0	0	0	1;
	48	49	0.0834	0.129	0.0048	9900	0	0	0	0	1;
	49	50	0.0801	0.128	0	9900	0	0	0	0	1;
	50	51	0.1386	0.22	0	9900	0	0	0	0	1;
	10	51	0	0.0712	0	9900	0	0	0.93	0	1;
	13	49	0	0.191	0	9900	0	0	0.895	0	1;
	29	52	0.1442	0.187	0	9900	0	0	0	0	1;
	52	53	0.0762	0.0984	0	9900	0	0	0	0	1;
	53	54	0.1878	0.232	0	9900	0	0	0	0	1;
	54	55	0.1732	0.2265	0	9900	0	0	0	0	1;
	11	43	0	0.153	0	9900	0	0	0.958	0	1;
	44	45	0.0624	0.1242	0.004	9900	0	0	0	0	1;
	40	56	0	1.195	0	9900	0	0	0.958	0	1;
	56	41	0.553	0.549	0	9900	0	0	0	0	1;
	56	42	0.2125	0.354	0	9900	0	0	0	0	1;
	39	57	0	1.355	0	9900	0	0	0.98	0	1;
	57	56	0.174	0.26	0	9900	0	0	0	0	1;
	38	49	0.115	0.177	0.003	9900	0	0	0	0	1;
	38	48	0.0312	0.0482	0	9900	0	0	0	0	1;
	9	55	0	0.1205	0	9900	0	0	0.94	0	1;
];

%%-----  OPF Data  -----%%
%% area data
areas = [
	1	1;
];

%% generator cost data
%	1	startup	shutdown	n	x0	y0	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
gencost = [
	2	0	0	3	0.0775795	20	0;
	2	0	0	3	0.01	40	0;
	2	0	0	3	0.25	20	0;
	2	0	0	3	0.01	40	0;
	2	0	0	3	0.0222222	20	0;
	2	0	0	3	0.01	40	0;
	2	0	0	3	0.0322581	20	0;
];

return;

% Warnings from cdf2matp conversion:
%
% ***** Qmax = Qmin at generator at bus    1 (Qmax set to Qmin + 10)
% ***** area data conversion not yet implemented (creating dummy area data)
% ***** Insufficient generation, setting Pmax at slack bus (bus 1) to 575.88
% ***** MVA limit of branch 1 - 2 not given, set to 9900
% ***** MVA limit of branch 2 - 3 not given, set to 9900
% ***** MVA limit of branch 3 - 4 not given, set to 9900
% ***** MVA limit of branch 4 - 5 not given, set to 9900
% ***** MVA limit of branch 4 - 6 not given, set to 9900
% ***** MVA limit of branch 6 - 7 not given, set to 9900
% ***** MVA limit of branch 6 - 8 not given, set to 9900
% ***** MVA limit of branch 8 - 9 not given, set to 9900
% ***** MVA limit of branch 9 - 10 not given, set to 9900
% ***** MVA limit of branch 9 - 11 not given, set to 9900
% ***** MVA limit of branch 9 - 12 not given, set to 9900
% ***** MVA limit of branch 9 - 13 not given, set to 9900
% ***** MVA limit of branch 13 - 14 not given, set to 9900
% ***** MVA limit of branch 13 - 15 not given, set to 9900
% ***** MVA limit of branch 1 - 15 not given, set to 9900
% ***** MVA limit of branch 1 - 16 not given, set to 9900
% ***** MVA limit of branch 1 - 17 not given, set to 9900
% ***** MVA limit of branch 3 - 15 not given, set to 9900
% ***** MVA limit of branch 4 - 18 not given, set to 9900
% ***** MVA limit of branch 4 - 18 not given, set to 9900
% ***** MVA limit of branch 5 - 6 not given, set to 9900
% ***** MVA limit of branch 7 - 8 not given, set to 9900
% ***** MVA limit of branch 10 - 12 not given, set to 9900
% ***** MVA limit of branch 11 - 13 not given, set to 9900
% ***** MVA limit of branch 12 - 13 not given, set to 9900
% ***** MVA limit of branch 12 - 16 not given, set to 9900
% ***** MVA limit of branch 12 - 17 not given, set to 9900
% ***** MVA limit of branch 14 - 15 not given, set to 9900
% ***** MVA limit of branch 18 - 19 not given, set to 9900
% ***** MVA limit of branch 19 - 20 not given, set to 9900
% ***** MVA limit of branch 21 - 20 not given, set to 9900
% ***** MVA limit of branch 21 - 22 not given, set to 9900
% ***** MVA limit of branch 22 - 23 not given, set to 9900
% ***** MVA limit of branch 23 - 24 not given, set to 9900
% ***** MVA limit of branch 24 - 25 not given, set to 9900
% ***** MVA limit of branch 24 - 25 not given, set to 9900
% ***** MVA limit of branch 24 - 26 not given, set to 9900
% ***** MVA limit of branch 26 - 27 not given, set to 9900
% ***** MVA limit of branch 27 - 28 not given, set to 9900
% ***** MVA limit of branch 28 - 29 not given, set to 9900
% ***** MVA limit of branch 7 - 29 not given, set to 9900
% ***** MVA limit of branch 25 - 30 not given, set to 9900
% ***** MVA limit of branch 30 - 31 not given, set to 9900
% ***** MVA limit of branch 31 - 32 not given, set to 9900
% ***** MVA limit of branch 32 - 33 not given, set to 9900
% ***** MVA limit of branch 34 - 32 not given, set to 9900
% ***** MVA limit of branch 34 - 35 not given, set to 9900
% ***** MVA limit of branch 35 - 36 not given, set to 9900
% ***** MVA limit of branch 36 - 37 not given, set to 9900
% ***** MVA limit of branch 37 - 38 not given, set to 9900
% ***** MVA limit of branch 37 - 39 not given, set to 9900
% ***** MVA limit of branch 36 - 40 not given, set to 9900
% ***** MVA limit of branch 22 - 38 not given, set to 9900
% ***** MVA limit of branch 11 - 41 not given, set to 9900
% ***** MVA limit of branch 41 - 42 not given, set to 9900
% ***** MVA limit of branch 41 - 43 not given, set to 9900
% ***** MVA limit of branch 38 - 44 not given, set to 9900
% ***** MVA limit of branch 15 - 45 not given, set to 9900
% ***** MVA limit of branch 14 - 46 not given, set to 9900
% ***** MVA limit of branch 46 - 47 not given, set to 9900
% ***** MVA limit of branch 47 - 48 not given, set to 9900
% ***** MVA limit of branch 48 - 49 not given, set to 9900
% ***** MVA limit of branch 49 - 50 not given, set to 9900
% ***** MVA limit of branch 50 - 51 not given, set to 9900
% ***** MVA limit of branch 10 - 51 not given, set to 9900
% ***** MVA limit of branch 13 - 49 not given, set to 9900
% ***** MVA limit of branch 29 - 52 not given, set to 9900
% ***** MVA limit of branch 52 - 53 not given, set to 9900
% ***** MVA limit of branch 53 - 54 not given, set to 9900
% ***** MVA limit of branch 54 - 55 not given, set to 9900
% ***** MVA limit of branch 11 - 43 not given, set to 9900
% ***** MVA limit of branch 44 - 45 not given, set to 9900
% ***** MVA limit of branch 40 - 56 not given, set to 9900
% ***** MVA limit of branch 56 - 41 not given, set to 9900
% ***** MVA limit of branch 56 - 42 not given, set to 9900
% ***** MVA limit of branch 39 - 57 not given, set to 9900
% ***** MVA limit of branch 57 - 56 not given, set to 9900
% ***** MVA limit of branch 38 - 49 not given, set to 9900
% ***** MVA limit of branch 38 - 48 not given, set to 9900
% ***** MVA limit of branch 9 - 55 not given, set to 9900


text:
function sample = SamplePose(P, G, label)

% sample from the distribution specified by P and G,
% label == 0: class label unknown
% label == k: class label = k; k=1,2,3, ...
%
% Copyright (C) Daphne Koller, Stanford Univerity, 2012

sample = zeros(10,3);

if label == 0
    k = SampleMultinomial(P.c); % sample label using P.c
else
    k = label;
end

visited = zeros(10,1);

while sum(visited) < 10

    for i=1:10 % for each body part

        par = G(i,2); % parent body part, if exists

        if G(i,1) == 0 % parametrized by (2),(3),(4)

            muy = P.clg(i).mu_y(k);
            mux = P.clg(i).mu_x(k);
            muangle = P.clg(i).mu_angle(k);
            sample(i,1) = SampleGaussian(muy, P.clg(i).sigma_y(k));
            sample(i,2) = SampleGaussian(mux, P.clg(i).sigma_x(k));
            sample(i,3) = SampleGaussian(muangle, P.clg(i).sigma_angle(k));
            visited(i) = 1;

        elseif G(i,1) == 1 % parametrized by (5),(6),(7)

            if visited(par)

                muy = P.clg(i).theta(k,1) + ...
                    P.clg(i).theta(k,2) * sample(par,1) + ...
                    P.clg(i).theta(k,3) * sample(par,2) + ...
                    P.clg(i).theta(k,4) * sample(par,3);
                mux = P.clg(i).theta(k,5) + ...
                    P.clg(i).theta(k,6) * sample(par,1) + ...
                    P.clg(i).theta(k,7) * sample(par,2) + ...
                    P.clg(i).theta(k,8) * sample(par,3);
                muangle = P.clg(i).theta(k,9) + ...
                    P.clg(i).theta(k,10) * sample(par,1) + ...
                    P.clg(i).theta(k,11) * sample(par,2) + ...
                    P.clg(i).theta(k,12) * sample(par,3);
                sample(i,1) = SampleGaussian(muy, P.clg(i).sigma_y(k));
                sample(i,2) = SampleGaussian(mux, P.clg(i).sigma_x(k));
                sample(i,3) = SampleGaussian(muangle, P.clg(i).sigma_angle(k));
                visited(i) = 1;
            end

        elseif G(i,1) == 2 % parametrized by (8),(9),(10)

            if visited(par)

                muy = P.clg(i).gamma(k,1) + ...
                    P.clg(i).gamma(k,2) * sample(par,1) + ...
                    P.clg(i).gamma(k,3) * sample(par,2) + ...
                    P.clg(i).gamma(k,4) * sin(sample(par,3)) + ...
                    P.clg(i).gamma(k,5) * cos(sample(par,3));
                mux = P.clg(i).gamma(k,6) + ...
                    P.clg(i).gamma(k,7) * sample(par,1) + ...
                    P.clg(i).gamma(k,8) * sample(par,2) + ...
                    P.clg(i).gamma(k,9) * sin(sample(par,3)) + ...
                    P.clg(i).gamma(k,10) * cos(sample(par,3));
                muangle = P.clg(i).gamma(k,11) + ...
                    P.clg(i).gamma(k,12) * sample(par,1) + ...
                    P.clg(i).gamma(k,13) * sample(par,2) + ...
                    P.clg(i).gamma(k,14) * sample(par,3);
                sample(i,1) = SampleGaussian(muy, P.clg(i).sigma_y(k));
                sample(i,2) = SampleGaussian(mux, P.clg(i).sigma_x(k));
                sample(i,3) = SampleGaussian(muangle, P.clg(i).sigma_angle(k));
                visited(i) = 1;
            end
        end

    end
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [symbol,labelx,labely] = sectionLabels(type)

switch lower(type)

  case 'alpha'

    symbol = '\alpha';
    labelx = '$\gamma$';
    labely = '$\beta$';

  case 'gamma'

    symbol = '\gamma';
    labelx = '$\alpha$';
    labely = '$\beta$';

  case 'phi1'

    symbol = '\varphi_1';
    labelx = '$\varphi_2$';
    labely = '$\Phi$';

  case 'phi2'

    symbol = '\varphi_2';
    labelx = '$\varphi_1$';
    labely = '$\Phi$';

  case 'sigma'

    symbol = '\sigma';
    labelx = '$\sigma_2$';
    labely = '$\Phi$';

  case 'omega'

    symbol = '\omega';
    labelx = '$\theta$';
    labely = '$\rho$';

  case 'axisangle'

    symbol = '\theta';
    labelx = '$\r$';
    labely = '$\Theta$';

end

author: mtex-toolbox

repo: mtex

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


text:
classdef Polyhedron

  properties
    A;
    b;
    Aeq = [];
    beq = [];
    vertices;
    has_vertices = false;
  end

  methods
    function obj = Polyhedron(A, b, Aeq, beq)
      if nargin < 3
        Aeq = [];
      end
      if nargin < 4
        beq = [];
      end
      obj.A = A;
      obj.b = b;
      obj.Aeq = Aeq;
      obj.beq = beq;
    end

    function vertices = getVertices(obj)
      obj = obj.reduce();
      if ~obj.has_vertices
        if exist('cddmex', 'file')
          H = struct('A', [obj.Aeq; obj.A], 'B', [obj.beq; obj.b], 'lin', (1:size(obj.Aeq, 1))');
          V = cddmex('extreme', H);
          obj.vertices = V.V';
        else
          obj.vertices = iris.thirdParty.polytopes.lcon2vert(obj.A, obj.b, obj.Aeq, obj.beq)';
        end
        obj.has_vertices = true;
      end
      vertices = obj.vertices;
    end

    function reduced_poly = reduce(obj)
      % Find a minimal representation of the polyhedron
      if ~exist('cddmex', 'file')
        error('IRIS:MissingDependency', 'This function requires the cddmex tool. The easiest way to get it is using tbxmanager: http://www.tbxmanager.com/');
      end
      H = struct('A', [obj.Aeq; obj.A], 'B', [obj.beq; obj.b], 'lin', (1:size(obj.Aeq, 1))');
      Hred = cddmex('reduce_h', H);
      assert(isempty(Hred.lin), 'as far as I know, Hred.lin should always be empty. That is, the reduced polyhedron should not contain equality constraints. -rdeits');
      reduced_poly = iris.Polyhedron(Hred.A, Hred.B);
    end

    function plotVertices(obj, varargin)
      vertices = obj.getVertices();
      if ~isempty(vertices)
        if size(vertices, 2) > 2
          k = convhull(vertices(1,:), vertices(2,:));
        else
          k = [1:size(vertices, 2), 1];
        end
        plot(vertices(1,k), vertices(2,k), varargin{:});
      end
    end

    function drawLCMGL(obj, lcmgl)
      lcmgl.glBegin(lcmgl.LCMGL_LINES);
      vertices = obj.getVertices();
      if size(vertices, 2) > 2
        k = convhull(vertices(1,:), vertices(2,:));
      else
        k = [1:size(vertices, 2), 1];
      end
      for j = 1:length(k)-1
        lcmgl.glVertex3d(vertices(1,k(j)), vertices(2,k(j)), vertices(3,k(j)));
        lcmgl.glVertex3d(vertices(1,k(j+1)), vertices(2,k(j+1)), vertices(3,k(j+1)));
      end
      lcmgl.glEnd();
    end

    function obj = normalize(obj)
      n = zeros(size(obj.A, 1), 1);
      for j = 1:size(obj.A, 1)
        n(j) = norm(obj.A(j,:));
        obj.A(j,:) = obj.A(j,:) / n(j);
        obj.b(j) = obj.b(j) / n(j);
      end
      obj.A = obj.A(n > 0, :);
      obj.b = obj.b(n > 0);

      n = zeros(size(obj.Aeq, 1), 1);
      for j = 1:size(obj.Aeq, 1)
        n(j) = norm(obj.Aeq(j,:));
        obj.Aeq(j,:) = obj.Aeq(j,:) / n(j);
        obj.beq(j) = obj.beq(j) / n(j);
      end
      obj.Aeq = obj.Aeq(n > 0, :);
      obj.beq = obj.beq(n > 0);
    end

  end

  methods(Static)
    function obj = fromVertices(vertices)
      [A, b] = iris.thirdParty.polytopes.vert2lcon(vertices');
      obj = iris.Polyhedron(A, b);
    end

    function obj = from2DVertices(vertices)
      assert(size(vertices, 1) == 2);
      x = vertices(1,:);
      y = vertices(2,:);
      k = convhull(x,y, 'simplify', true);
      A = [(y(k(2:end)) - y(k(1:end-1)))', (x(k(1:end-1)) - x(k(2:end)))'];
      b = sum(A' .* [x(k(1:end-1)); y(k(1:end-1))], 1)';
      obj = iris.Polyhedron(A, b);
    end

    function obj = from2DVerticesAndPlane(vertices, normal, v)
      assert(size(vertices, 1) == 2);
      % normal' * [x;y;z] = v;
      obj = iris.Polyhedron.from2DVertices(vertices);
      obj.Aeq = reshape(normal, 1, []);
      obj.beq = v;
      obj.A = [obj.A, zeros(size(obj.A, 1), 1)];
    end

    function obj = fromBounds(lb, ub)
      % create a polyhedron representing a bounding box in n dimensions
      dim = length(lb);
      assert(length(lb) == length(ub));
      A = [eye(dim); -eye(dim)];
      b = [reshape(ub,[],1); reshape(-lb,[],1)];
      obj = iris.Polyhedron(A, b);
    end
  end
end

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function i4block_print_test ( )

%*****************************************************************************80
%
%% I4BLOCK_PRINT_TEST tests I4BLOCK_PRINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2012
%
%  Author:
%
%    John Burkardt
%
  l = 4;
  m = 3;
  n = 2;
  xvec = [ ...
        1,  2,  3,   4,  1, ...
        4,  9, 16,   1,  8, ...
       27, 64,  2,   4,  6, ...
        8,  2,  8,  18, 32, ...
        2, 16, 54, 128 ];

  x = reshape ( xvec, l, m, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4BLOCK_PRINT_TEST\n' );
  fprintf ( 1, '  I4BLOCK_PRINT prints an I4BLOCK.\n' );

  i4block_print ( l, m, n, x, '  The 3D array:' )

  return
end


text:
function data = read_1_pinnacle_file(filename,display)
%
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%
data = [];
if ~exist('display','var')
	display = 1;
end

if display == 1
	h = waitbar(0,sprintf('Reading file "%s" ...',filename));
end

fid = fopen(filename);
if fid < 0
	fprintf('File "%s" cannot be opened.\n',filename);
    if display == 1
        close(h);
    end
	return;
end

%try
	while 1
		[fieldname,val] = read_1_field_pinnacle(fid);
		if isempty(fieldname) && isempty(val)
			break;
		else
			data = add_1_field_to_struct(data,fieldname,val);
		end
	end
	fclose(fid);
% catch
% 	fclose(fid);
% 	disp(lasterr);
% end

if display == 1
	close(h);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [B]=ceven(C)
%these function choice element of vector with even indices of column
k=1;
for i=2:2:length(C)
    B(1,k)=C(1,i);
    k=k+1;
end


text:
% Compute the log intensity for the inverse link function g(f) = 1/(1+exp(-f)).
%
% The function is used in GLM likelihoods such as likPoisson, likGamma, likBeta
% and likInvGauss.
%
% Copyright (c) by Hannes Nickisch, 2013-10-16.

function varargout = glm_invlink_logit(f)
  varargout = cell(nargout, 1);  % allocate the right number of output arguments
  [varargout{:}] = glm_invlink_logistic(f);
  if nargout>0
    elg = exp(varargout{1});
    varargout{1} = f - elg;
    if nargout>1
      dlg = varargout{2};
      varargout{2} = 1 - elg.*dlg;
      if nargout>2
        d2lg = varargout{3};
        varargout{3} = -elg.*(dlg.^2+d2lg);
        if nargout>3
          d3lg = varargout{4};
          varargout{4} = -elg.*(dlg.^3+3*d2lg.*dlg+d3lg);
        end
      end
    end
  end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function ODG = PQnNetB (MOV)
% Neural net to get the final ODG

% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:27:44 $

persistent amin amax wx wxb wy wyb bmin bmax I J CLIPMOV
global PQopt

if (isempty (amin))
    I = length (MOV);
    if (I == 11)
        [amin, amax, wx, wxb, wy, wyb, bmin, bmax] = NNetPar ('Basic');
    else
        [amin, amax, wx, wxb, wy, wyb, bmin, bmax] = NNetPar ('Advanced');
    end
    [I, J] = size (wx);
end

sigmoid = inline ('1 / (1 + exp(-x))');

% Scale the MOV's
Nclip = 0;
MOVx = zeros (1, I);
for (i = 0:I-1)
    MOVx(i+1) = (MOV(i+1) - amin(i+1)) / (amax(i+1) - amin(i+1));
    if (~ isempty (PQopt) & PQopt.ClipMOV ~= 0)
        if (MOVx(i+1) < 0)
            MOVx(i+1) = 0;
            Nclip = Nclip + 1;
        elseif (MOVx(i+1) > 1)
            MOVx(i+1) = 1;
            Nclip = Nclip + 1;
        end
    end
end
if (Nclip > 0)
    fprintf ('>>> %d MOVs clipped\n', Nclip);
end

% Neural network
DI = wyb;
for (j = 0:J-1)
    arg = wxb(j+1);
    for (i = 0:I-1)
        arg = arg + wx(i+1,j+1) * MOVx(i+1);
    end
    DI = DI + wy(j+1) * sigmoid (arg);
end

ODG = bmin + (bmax - bmin) * sigmoid (DI);

function [amin, amax, wx, wxb, wy, wyb, bmin, bmax] = NNetPar (Version)

if (strcmp (Version, 'Basic'))
    amin = ...
        [393.916656,  361.965332,  -24.045116,    1.110661,  -0.206623, ...
           0.074318,    1.113683,    0.950345,    0.029985,   0.000101, ...
           0];
    amax = ...
        [921,         881.131226,   16.212030,  107.137772,   2.886017, ...
          13.933351,   63.257874, 1145.018555,   14.819740,   1,        ...
           1];
    wx = ...
        [ [ -0.502657,  0.436333,   1.219602 ];
          [  4.307481,  3.246017,   1.123743 ];
          [  4.984241, -2.211189,  -0.192096 ];
          [  0.051056, -1.762424,   4.331315 ];
          [  2.321580,  1.789971,  -0.754560 ];
          [ -5.303901, -3.452257, -10.814982 ];
          [  2.730991, -6.111805,   1.519223 ];
          [  0.624950, -1.331523,  -5.955151 ];
          [  3.102889,  0.871260,  -5.922878 ];
          [ -1.051468, -0.939882,  -0.142913 ];
          [ -1.804679, -0.503610,  -0.620456 ] ];
    wxb = ...
        [ -2.518254,  0.654841, -2.207228 ];
    wy = ...
        [ -3.817048,  4.107138,  4.629582 ];
    wyb = -0.307594;
    bmin = -3.98;
    bmax = 0.22;
else
    amin = ...
        [   13.298751,  0.041073,  -25.018791,  0.061560,   0.024523 ];
    amax = ...
        [ 2166.5,      13.24326,    13.46708,  10.226771,  14.224874 ];
    wx = ...
        [ [  21.211773, -39.913052, -1.382553, -14.545348,  -0.320899 ];
          [  -8.981803,  19.956049,  0.935389,  -1.686586,  -3.238586 ];
          [   1.633830,  -2.877505, -7.442935,   5.606502,  -1.783120 ];
          [   6.103821,  19.587435, -0.240284,   1.088213,  -0.511314 ];
          [  11.556344,   3.892028,  9.720441,  -3.287205, -11.031250 ] ];
    wxb = ...
        [   1.330890,   2.686103,    2.096598, -1.327851,   3.087055 ];
    wy = ...
        [ -4.696996,   -3.289959,    7.004782,  6.651897,   4.009144 ];
    wyb = -1.360308;
    bmin = -3.98;
    bmax = 0.22;
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
if ~exist('fig', 'dir')
    mkdir('fig');
end
fig1;
fig2;
fig3;
fig4;
fig5;

author: zhoupc


text:
function [iV,jV,kV]=find3d(mask3M)
%"find3d"
%   This is the 3D equivalent of the builtin command, find.
%   It returns the i,j,k indices of non-zero entries.
%
%JOD, 16 Nov 98
%JOD, 26 Feb 03, bugfix.
%JRA, 26 Feb 04, new algorithm, also implements the new fastind2sub.
%
%Usage:
%   [iV,jV,kV]=find3d(mask3M);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indV = find(mask3M(:));
[iV,jV,kV] = fastind2sub(size(mask3M), indV);
iV = iV';
jV = jV';
kV = kV';


text:
%% Copyright (C) 2016 Lagu
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

function r = octave_array_to_python(a)
  t = '';

  [numrows, numcols] = size (a);

  if (numrows == 0)
    r = mat2str([]);
    return;

  elseif (numrows > 1)
    t = strcat('[', t);
    t = strcat(t, octave_array_to_python(a(1, :)));

    for i = 2:numrows
      t = strcat(t, ', ');
      t = strcat(t, octave_array_to_python(a(i, :)));
    end

    t = strcat(t, ']');

  else
    t = strcat('[', t);
    t = strcat(t, mat2str(a(1, 1)));

    for i = 2:numcols
      t = strcat(t, ', ');
      t = strcat(t, mat2str(a(1, i)));
    end

    t = strcat(t, ']');

  end

  r = t;

repo: octsympy


text:
function h = p10_fh ( p, varargin )

%*****************************************************************************80
%
%% P10_FH returns a mesh size function for problem 10.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    06 February 2006
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real P(NP,ND), the point coordinates.
%
%    Input, VARARGIN, room for extra arguments.
%
%    Output, real H(NP,1), the mesh size function.
%
  np = size ( p, 1 );
  h = ones ( np, 1 );

  return
end


text:
function figure_handle = PTKShow2DSlicesInOneFigure(viewer_panel_handle, orientation, skip_sices, reporting)
    % PTKShowAll2DSlicesInOneFigure. Creates a figure showing every slice from a
    %     MimViewerPanel object displayed in a grid
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    if nargin < 4
        reporting = CoreReportingDefault;
    end

    % Number of plots along x axis
    num_plots_x = 3;

    % Maximum number of plots along y axis
    max_plots_y = 4;

    figure_handle = figure;
    axis off;
    hold on;
    gap = 0.01;

    number_slices = viewer_panel_handle.BackgroundImage.ImageSize(orientation);

    % Ignore slices at beginning and end
    number_slices = number_slices - 2*skip_sices;

    max_images_on_page = num_plots_x*max_plots_y;
    slice_spacing = ceil(number_slices/max_images_on_page);
    number_slices_shown = floor(number_slices/slice_spacing);

    num_plots_y  = ceil(number_slices_shown/num_plots_x);

    panel_slice_number = [1, 1, 1];
    viewer_panel_handle.Orientation = orientation;
    figure_num = 0;

    for slice_num = skip_sices : slice_spacing : number_slices + skip_sices

        figure_num = figure_num + 1;

        % Make viewer display this slice
        panel_slice_number(orientation) = slice_num;
        viewer_panel_handle.SliceNumber = panel_slice_number;

        % Capture the frame (including overlay)
        frame = viewer_panel_handle.Capture;

        yc = floor((figure_num - 1)/num_plots_x);
        xc = mod(figure_num - 1, num_plots_x);
        y_pos = 1 - yc/num_plots_y;
        y_size =  1/num_plots_y;
        x_pos = xc/num_plots_x;
        x_size =  1/num_plots_x;
        im = frame2im(frame);

        % Make figure active
        figure(figure_handle);
        pos = [x_pos + gap, y_pos - y_size + gap, x_size - 2*gap, y_size - 2*gap];
        axes_handle = axes('Position', pos, 'Visible', 'off', 'DataAspectRatio', [1 1 1]);
        image(im);
        set(axes_handle, 'DataAspectRatio', [1 1 1])
        hold on;
        axis off;
        set(axes_handle, 'Units', 'normalized', 'OuterPosition', pos, 'Position', pos);
        set(axes_handle, 'Units', 'normalized', 'OuterPosition', pos, 'Position', pos);
    end

    % Choose a figure size that fils the vertical screen size, but with the
    % correct ratio
    display_dimensions = CoreSystemUtilities.GetMonitorDimensions;

    % Adjust to take into account toolbars etc.
    display_height = display_dimensions(2) - 300;
    proportional_width = display_height*(size(im,2)/size(im,1))*(num_plots_x/num_plots_y)*(1 - 2*gap*num_plots_y)/(1 - 2*gap*num_plots_x);

    set(figure_handle, 'Units', 'pixels', 'Position', [0 0 proportional_width, display_height]);
    set(figure_handle, 'PaperPositionMode', 'auto');
end

author: tomdoel


load '..\uavW.mat'
i1 = uavW;
clear uavW;


text:
function disp( x, prefix )
if nargin < 2,
    prefix = '';
end
disp( [ prefix, 'cvx tuple object: ' ] );
prefix = [ prefix, '   ' ];
do_disp( x.value_, {}, prefix, prefix, '' );

function do_disp( x, f, fprefix, prefix, suffix )
switch class( x ),
    case 'struct',
        do_disp( struct2cell(x), fieldnames(x), fprefix, prefix, suffix );
    case 'cell',
        fprefix = [ fprefix, '{ ' ];
        prefix  = [ prefix, '  ' ];
        nsuffix = '';
        kend = numel( x );
        for k = 1 : kend,
            if k == kend, nsuffix = [ ' }', suffix ]; end
            if ~isempty( f ),
                fprefix = sprintf( '%s%s: ', fprefix, f{k} );
            end
            do_disp( x{k}, {}, fprefix, prefix, nsuffix );
            fprefix = prefix;
        end
    case 'cvx',
        dual = getdual( x );
        if ~isempty( dual ),
            suffix = sprintf( ' (dual: %s)%s', dual, suffix );
        end
        disp( [ fprefix, cvx_class( x, true, true ), ' ', type( x ), suffix ] );
    case 'double',
        fprintf( 1, '%s%g%s\n', fprefix, x, suffix );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
dtDirs = {'/biac3/wandell4/data/reading_longitude/dti_adults/ak090724/dti40trilin/'...
    '/biac3/wandell4/data/reading_longitude/dti_adults/am090121/dti06trilinrt/'...
    '/biac3/wandell4/data/reading_longitude/dti_adults/rfd080930/dti06trilinrt/'};

BPFDirs = {'/biac3/wandell5/data/relaxometry/ak090721/trilin_ak_2mms'...
    '/biac3/wandell5/data/relaxometry/am_090126/trilin2mml'...
    '/biac3/wandell5/data/relaxometry/rb_090716/trilin_rb_2mms'};

sub={'ak' 'am' 'rfd'}

for j=1:1%2 %1.5T or 3T anatomy
    for i=2:2%length(sub),
        name3T=sub{i};
        name15T=[sub{i} '1_5'];
        path3=['/biac3/wandell5/data/relaxometry/100405HLF3T/anal/' name3T 'lin_2mms1/'];
        path15=['/biac3/wandell5/data/relaxometry/100405HLF3T/anal/' name15T 'lin_2mms'];%Nbs1

        if j==1,
            path=path15;
        elseif j==2
            path=path3;
        end;
        %%% load maps  %%%
        hlf=niftiRead([path '/HLF_F.nii.gz']);
        t1=niftiRead([path '/T1_LFit_F.nii.gz']);
        T1hWf=niftiRead([path '/T1fh_F.nii.gz']);
        wf=niftiRead([path '/Wf_F.nii.gz']);
        dt = dtiLoadDt6(fullfile(dtDirs{i},'dt6'));
        bpf=niftiRead(fullfile(BPFDirs{i},'f.nii.gz'));
       Class=niftiRead([path15 '/T1_class2mm.nii.gz']);

%        t1=t1.data;
%        hlf=hlf.data;
%        wf=wf.data;
%        T1hWf=T1hWf.data;
%        bpf=bpf.data;
%
      [eigVec,eigVal] = dtiEig(dt.dt6);
    [fa,md,rd] = dtiComputeFA(eigVal);
    [cl, cp, cs] = dtiComputeWestinShapes(eigVal);

%       we will make wm masks for the maps

maskAnat=find(Class.data==3 | Class.data==4);

if(~all(dt.xformToAcpc(:)==Class.qto_xyz(:)))
    maskR= mrAnatResliceSpm(double(Class.data),inv(Class.qto_xyz),dt.bb,dt.mmPerVoxel,1);
    maskDif=find(maskR==3 | maskR==4);
    clear maskR;

else
    maskDif=maskAnat;
end;

%compute means andsd's for each measure within white matter mask
fa_std=std(fa(maskDif(~isnan(fa(maskDif)))));
rd_std=std(rd(maskDif(~isnan(rd(maskDif)))));
md_std=std(md(maskDif(~isnan(md(maskDif)))));
cl_std=std(cl(maskDif(~isnan(cl(maskDif)))));

hlf_std=std(hlf.data(maskAnat(~isnan(hlf.data(maskAnat)))));
wf_std=std(wf.data(maskAnat(~isnan(wf.data(maskAnat)))));
t1_std=std(t1.data(maskAnat(~isnan(t1.data(maskAnat)))));
T1hWf_std=std(T1hWf.data(maskAnat(~isnan(T1hWf.data(maskAnat)))));
bpf_std=std(bpf.data(maskAnat(~isnan(bpf.data(maskAnat)))));

fa_mean=mean(fa(maskDif(~isnan(fa(maskDif)))));
rd_mean=mean(rd(maskDif(~isnan(rd(maskDif)))));
md_mean=mean(md(maskDif(~isnan(md(maskDif)))));
cl_mean=mean(fa(maskDif(~isnan(cl(maskDif)))));

hlf_mean=mean(hlf.data(maskAnat(~isnan(hlf.data(maskAnat)))));
wf_mean=mean(wf.data(maskAnat(~isnan(wf.data(maskAnat)))));
t1_mean=mean(t1.data(maskAnat(~isnan(t1.data(maskAnat)))));
T1hWf_mean=mean(T1hWf.data(maskAnat(~isnan(T1hWf.data(maskAnat)))));
bpf_mean=mean(bpf.data(maskAnat(~isnan(bpf.data(maskAnat)))));

%load morigroups and rois
%fg will contain all mori groups
fbDir=[dtDirs{i} 'fibers/MoriGroups'];

fg=dtiReadFibers(fbDir)

roiDir=[dtDirs{i} 'ROIs'];

fgName(3)={'Left Cortico-Spinal'};
roi1{3}=dtiReadRoi(fullfile(roiDir,'CST_roi1_L'));
roi2{3}=dtiReadRoi(fullfile(roiDir,'CST_roi2_L'));
fgName(4)={'Right Cortico-Spinal'};
roi1{4}=dtiReadRoi(fullfile(roiDir,'CST_roi1_R'));
roi2{4}=dtiReadRoi(fullfile(roiDir,'CST_roi2_R'));
fgName(11)={'Left Inferior Fronto-Occ'};
roi1{11}=dtiReadRoi(fullfile(roiDir,'IFO_roi1_L'));
roi2{11}=dtiReadRoi(fullfile(roiDir,'IFO_roi2_L'));
fgName(12)={'Right Inferior Frontal-Occ'};
roi1{12}=dtiReadRoi(fullfile(roiDir,'IFO_roi1_R'));
roi2{12}=dtiReadRoi(fullfile(roiDir,'IFO_roi2_R'));
fgName(13)={'Left Inferior Longitude'};
roi1{13}=dtiReadRoi(fullfile(roiDir,'ILF_roi1_L'));
roi2{13}=dtiReadRoi(fullfile(roiDir,'ILF_roi2_L'));
fgName(14)={'Right Inferior Longitude'};
roi1{14}=dtiReadRoi(fullfile(roiDir,'ILF_roi1_R'));
roi2{14}=dtiReadRoi(fullfile(roiDir,'ILF_roi2_R'))
fgName(19)={'Left Arcuate'};
roi1{19}=dtiReadRoi(fullfile(roiDir,'SLF_roi1_L'));
roi2{19}=dtiReadRoi(fullfile(roiDir,'SLFt_roi2_L'));
fgName(20)={'Right Arcuate'};
roi1{20}=dtiReadRoi(fullfile(roiDir,'SLF_roi1_R'));
roi2{20}=dtiReadRoi(fullfile(roiDir,'SLFt_roi2_R'));
%compute properties on tract trajectory for the specific fiber goups of
%interest

groups=[3  11  13  19 4 12 14 20];
for ii=groups;
    %compute fa and md along the trajectory of the fiber group as a
    %weighted average of fa and md at each fiber
    [fa_ md_ rd_ d cl_]=dtiComputeDiffusionPropertiesAlongFG(fg(ii), dt, roi1{ii}, roi2{ii}, 30);
    %clip fiber group between 2 rois to obtain core segment
    fgClipped = dtiClipFiberGroupToROIs(fg(ii),roi1{ii},roi2{ii});
    %compute quantitative measures along trajectory of fiber group
    [hlf_, SuperFiber, weightsNormalized] =dtiFiberGroupPropertyWeightedAverage(fgClipped, hlf, 30,'image');
    [wf_, SuperFiber, weightsNormalized] =dtiFiberGroupPropertyWeightedAverage(fgClipped, wf, 30,'image');
    [t1_, SuperFiber, weightsNormalized] =dtiFiberGroupPropertyWeightedAverage(fgClipped, t1, 30,'image');
    [T1hWf_, SuperFiber, weightsNormalized] =dtiFiberGroupPropertyWeightedAverage(fgClipped, T1hWf, 30,'image');
    [bpf_, SuperFiber, weightsNormalized] =dtiFiberGroupPropertyWeightedAverage(fgClipped, bpf, 30,'image');

    %Standardize measures based on variance within white matter mask
    faZ=((fa_-fa_mean)./fa_std);
    rdZ=((rd_-rd_mean)./rd_std);
    mdZ=((md_-md_mean)./md_std);
    clZ=((cl_-cl_mean)./cl_std);

    hlfZ=(hlf_-hlf_mean)./hlf_std;
    T1hWfZ=(T1hWf_-T1hWf_mean)./T1hWf_std;
    t1Z=(t1_-t1_mean)./t1_std;
    wfZ=(wf_-wf_mean)./wf_std;
    bpfZ=(bpf_-bpf_mean)./bpf_std;

    %Make plot for fiber group ii
    figure(ii);
     subplot(2,4,find(groups==ii));

       subplot(2,1,1);
   plot(horzcat(faZ,rdZ,mdZ,clZ),'LineWidth',2);axis([0 30 -2 2]);title([sub{i} 'Dif' fgName{ii}])
         legend('fa', 'rd', 'md', 'cl');ylabel('Z Score');
   subplot(2,1,2);
    plot(horzcat(hlfZ,T1hWfZ,t1Z,wfZ,bpfZ),'LineWidth',2);axis([0 30 -2 2]);title([sub{i} 'Anat' fgName{ii}])
   %if find(groups==ii)==1 | find(groups==ii)==5
        legend( 'hlf', 'T1hWf', 't1', 'wf' ,'bpf'); ylabel('Z Score');
   % end

end

    end

end;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function pred = nnrpred(Xtest,X,y,K,dist_type,pret_type)

% prediction of new samples with knn regression model
%
% pred = nnrpred(Xtest,X,y,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)
% y:            training y vector (n x 1)
% K:            number of neighbors
% dist_type:    'euclidean' Euclidean distance
%               'mahalanobis' Mahalanobis distance
%               'cityblock' City Block metric
%               'minkowski' Minkowski metric
% 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
% y_pred            predicted y vector [n_test x 1]
% y_pred_weighted   predicted wighted y vector [n_test x 1]
% neighbors         list of k neighbors for each predicted sample [n_test x k]
%
% version 2.0 - may 2012
% Davide Ballabio
% Milano Chemometrics and QSAR Research Group
% www.disat.unimib.it/chm

[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 = nnr_calc_dist(X_scal_train,X_scal,dist_type,pret_type);
neighbors = zeros(n,K);
yc=zeros(1,n);
yc_weighted=zeros(1,n);
w=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);
    [yc(i),yc_weighted(i),w(i,:)] = nnrcalcy(y(neighbors(i,:)),d_neighbors,K);
    dc(i,:)=d_neighbors;
end

pred.neighbors  = neighbors;
pred.y_pred = yc';
pred.y_pred_weighted = yc_weighted';
pred.D=D;
pred.dc=dc;
pred.w=w;

author: kmansouri

repo: OPERA

save_path: github-repos/MATLAB/kmansouri-OPERA


text:
function u = ref_to_koorn ( r )

%*****************************************************************************80
%
%% REF_TO_KOORN maps points from the reference to Koornwinder's triangle.
%
%  Discussion:
%
%    The reference triangle has vertices:
%
%      ( -1, -1/sqrt(3) )
%      ( +1, -1/sqrt(3) )
%      (  0, +2/sqrt(3) )
%
%    Koornwinder's triangle has vertices:
%
%      ( -1, -1 )
%      ( +1, -1 )
%      ( -1, +1 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R(2), the coordinates of a point in the
%    reference triangle.
%
%    Output, real U(2), the coordinates of the point in
%    the Koornwinder triangle.
%
  a10 = -1.0 / 3.0;
  a11 =  1.0;
  a12 = -1.0 / sqrt ( 3.0 );

  a20 = - 1.0 / 3.0;
  a21 =   0.0;
  a22 =   2.0 * sqrt ( 3.0 ) / 3.0;

  u(1) = a10 + r(1) + a12 * r(2);
  u(2) = a20 +        a22 * r(2);

  return
end


text:
function [gradISE,gradJhr,gradJrr]=GaussMixISEMeanGrads(w1,mu1,P1,w2,mu2,P2,PDFVals12,PDFVals22)
%%GAUSSMIXISEMEANGRADS Compute the gradient of the non-normalized
%               integrated squared error (ISE) between two Gaussian mixture
%               PDFs with respect to the mean vectors of the individual
%               components of the second Gaussian mixture. This gradient
%               can arise when optimizing over the ISE with respect to the
%               means.
%
%INPUTS: w1 The N1X1 or 1XN1 vector of weights for the first Gaussian
%           mixture. All w1>0 and sum(w1)=1.
%       mu1 The xDimXN1 set of mean vectors for the first Gaussian mixture
%           distribution.
%        P1 The xDimXxDimXN1 set of positive definite covariance matrices
%           for the first Gaussian mixture distirbution.
% w2, mu2, P2, The length N2, xDimXN2, and xDimXxDimXN2 set of weights,
%           mean vectors and positive-definite covariance matrices for the
%           second Gaussian mixture distribution.
% PDFVals12 A matrix such that the value in element (i,j) is
%           N(mu1(:,i);mu2(:,j),P1(:,:,i)+P2(:,:,j)), where N indicates the
%           multivariate Gaussian PDF evaluated at the first argument with
%           the second and third arguments being the mean and covarince
%           matrix. This parameter is returned by computeGaussMixISE.
% PDFVals22 A matrix such that the value in element (i,j) is
%           N(mu2(:,i);mu2(:,j),P2(:,:,i)+P2(:,:,j)). This parameter is
%           returned by computeGaussMixISE.
%
%OUTPUTS: gradISE The xDimXn2 set of derivatives of the ISE with respect to
%                 the elements of q (the square roots of w2).
% gradJhr, gradJrr In Chapter 3 of [1], the ISE is expressed in terms of
%                 Jhr and Jrr terms. These are the xDimXn2 gradients of
%                 those terms.
%
%Formule for gradJhr and gradJrr are Equation 3.34 in Section 3.3.3.2 of
%[1]. They relate to the ISE via Equation 3.20. See the function
%computeGaussMixISE to compute the ISE.
%
%EXAMPLE 1:
%In this example with a scalar random variable, we verify that the gradient
%obtained from this function is consistent with numerical differentiation.
% w1=[0.03,0.18,0.12,0.19,0.02,0.16,0.06,0.1,0.08,0.06];
% n1=length(w1);
% mu1=[1.45,2.20,0.67,0.48,1.49,0.91,1.01,1.42,2.77,0.89];
% P1=[0.0487,0.0305,0.1171,0.0174,0.0295,0.0102, 0.0323, 0.0380, 0.0115, 0.0679];
% P1=reshape(P1,[1,1,n1]);
%
% %The second PDF is the first with the five least-weight components deleted.
% w2=[0.18,0.12,0.19,0.16,0.1,0.08];
% w2=w2/sum(w2);
% n2=length(w2);
% mu2=[2.20,0.67,0.48,0.91,1.42,2.77];
% P2=[0.0305,0.1171,0.0174,0.0102,0.0380,0.0115];
% P2=reshape(P2,[1,1,n2]);
%
% [ISEVal,PDFVals12,PDFVals22]=computeGaussMixISE(w1,mu1,P1,w2,mu2,P2);
%
% epsVal=1e-8;
% gradISENumDiff=zeros(1,n2);
% for j=1:n2
%     muCur=mu2;
%     muCur(1,j)=muCur(1,j)+epsVal;
% 	ISEValCur=computeGaussMixISE(w1,mu1,P1,w2,muCur,P2);
%     gradISENumDiff(1,j)=(ISEValCur-ISEVal)/epsVal;
% end
% gradISE=GaussMixISEMeanGrads(w1,mu1,P1,w2,mu2,P2,PDFVals12,PDFVals22);
% RelErr=max(abs((gradISENumDiff-gradISE)./gradISENumDiff))
%The relative error will be about 5.4042e-6, which indicates good numeric
%agreement.
%
%EXAMPLE 2:
%In this example with a bivariate distribution, we verify that the gradient
%obtained from this function is consistent with numerical differentiation.
% w1=[0.25;0.5;0.25];
% mu1=zeros(2,2);
% mu1(:,1)=[1;-1];
% mu1(:,2)=[-1;1];
% mu1(:,3)=[0;0];
% P1=zeros(2,2,2);
% P1(:,:,1)=[4/9,  14/45;
%            14/45,4/9];
% P1(:,:,2)=[4/9, 0;
%            0, 4/9];
% P1(:,:,3)=[2/9, -1/9;
%           -1/9, 3/9];
%
% %The second distribution just throws out the first component.
% w2=w1(2:3);
% w2=w2/sum(w2);
% mu2=mu1(:,2:3);
% P2=P1(:,:,2:3);
%
% n2=length(w2);
% [ISEVal,PDFVals12,PDFVals22]=computeGaussMixISE(w1,mu1,P1,w2,mu2,P2);
%
% numDim=size(mu1,1);
% epsVal=1e-8;
% gradISENumDiff=zeros(numDim,n2);
% for curDim=1:numDim
%     for j=1:n2
%         muCur=mu2;
%         muCur(curDim,j)=muCur(curDim,j)+epsVal;
%         ISEValCur=computeGaussMixISE(w1,mu1,P1,w2,muCur,P2);
%         gradISENumDiff(curDim,j)=(ISEValCur-ISEVal)/epsVal;
%     end
% end
% gradISE=GaussMixISEMeanGrads(w1,mu1,P1,w2,mu2,P2,PDFVals12,PDFVals22);
% RelErr=max(max(abs((gradISENumDiff-gradISE)./gradISENumDiff)))
%The relative error will be about 1.2876e-7, which indicates good numeric
%agreement.
%
%REFERENCES:
%[1] J. L. Williams, "Gaussian mixture reduction for tracking multiple
%    maneuvering targets in clutter," Master's thesis, Air Force Institute
%    of Technology, Mar. 2003. [Online].
%    Available: http://www.dtic.mil/srch/doc?collection=t3&id=ADA415317
%
%May 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

Nh=length(w1);
Nr=length(w2);
xDim=size(mu1,1);

%The formulae for the gradient of Jhr and Jrr are given in Equation 3.34 of
%Section 3.3.3.2 of [1].
gradJhr=zeros(xDim,Nr);
gradJrr=zeros(xDim,Nr);

for j=1:Nr
    for i=1:Nh
        gradJhr(:,j)=gradJhr(:,j)+(P1(:,:,i)+P2(:,:,j))\(mu2(:,j)-mu1(:,i))*PDFVals12(i,j)*w1(i);
    end
    gradJhr(:,j)=-w2(j)*gradJhr(:,j);
end

for j=1:Nr
    for i=1:Nr
        gradJrr(:,j)=gradJrr(:,j)+(P2(:,:,i)+P2(:,:,j))\(mu2(:,j)-mu2(:,i))*PDFVals22(i,j)*w2(i);
    end
    gradJrr(:,j)=-2*w2(j)*gradJrr(:,j);
end

gradISE=gradJrr-2*gradJhr;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function str = stripzeros(strin)
     %STRIPZEROS Strip trailing zeros, leaving one digit right of decimal point.
     % Remove trailing zeros while leaving at least one digit to the right of
     % the decimal place.

     %   Copyright 2010 The MathWorks, Inc.

     str = strin;
     n = regexp(str,'\.0*$');
    if ~isempty(n)
        % There is nothing but zeros to the right of the decimal place;
        % the value in n is the index of the decimal place itself.
        % Remove all trailing zeros except for the first one.
        str(n+2:end) = [];
    else
        % There is a non-zero digit to the right of the decimal place.
        m = regexp(str,'0*$');
        if ~isempty(m)
            % There are trailing zeros, and the value in m is the index of
            % the first trailing zero. Remove them all.
            str(m:end) = [];
        end
    end

author: Arrowstar

repo: ksptot