text:
function pts = intersectPolylines(poly1, varargin)
%INTERSECTPOLYLINES Find the common points between 2 polylines
%
%   INTERS = intersectPolylines(POLY1, POLY2)
%   Returns the intersection points between two polylines. Each polyline is
%   defined by a N-by-2 array representing coordinates of its vertices:
%   [X1 Y1 ; X2 Y2 ; ... ; XN YN]
%   INTERS is a NP-by-2 array containing coordinates of intersection
%   points.
%
%   INTERS = intersectPolylines(POLY1)
%   Compute self-intersections of the polyline.
%
%   Example
%   % Compute intersection points between 2 simple polylines
%     poly1 = [20 10 ; 20 50 ; 60 50 ; 60 10];
%     poly2 = [10 40 ; 30 40 ; 30 60 ; 50 60 ; 50 40 ; 70 40];
%     pts = intersectPolylines(poly1, poly2);
%     figure; hold on;
%     drawPolyline(poly1, 'b');
%     drawPolyline(poly2, 'm');
%     drawPoint(pts);
%     axis([0 80 0 80]);
%
%   This function is largely based on the 'interX' function, found on the
%   FileExchange:
%   https://fr.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections
%
%   See also
%   polygons2d, polylineSelfIntersections, intersectLinePolygon
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2009-06-15,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

% The code is a slight rewritting of the interX function, consisting in
% avoiding argument transposition in the begining of the function. Comment
% of original submission are kept here for information.
%
%INTERX Intersection of curves
%   P = INTERX(L1,L2) returns the intersection points of two curves L1
%   and L2. The curves L1,L2 can be either closed or open and are described
%   by two-row-matrices, where each row contains its x- and y- coordinates.
%   The intersection of groups of curves (e.g. contour lines, multiply
%   connected regions etc) can also be computed by separating them with a
%   column of NaNs as for example
%
%         L  = [x11 x12 x13 ... NaN x21 x22 x23 ...;
%               y11 y12 y13 ... NaN y21 y22 y23 ...]
%
%   P has the same structure as L1 and L2, and its rows correspond to the
%   x- and y- coordinates of the intersection points of L1 and L2. If no
%   intersections are found, the returned P is empty.
%
%   P = INTERX(L1) returns the self-intersection points of L1. To keep
%   the code simple, the points at which the curve is tangent to itself are
%   not included. P = INTERX(L1,L1) returns all the points of the curve
%   together with any self-intersection points.
%
%   Example:
%       t = linspace(0,2*pi);
%       r1 = sin(4*t)+2;  x1 = r1.*cos(t); y1 = r1.*sin(t);
%       r2 = sin(8*t)+2;  x2 = r2.*cos(t); y2 = r2.*sin(t);
%       P = InterX([x1;y1],[x2;y2]);
%       plot(x1,y1,x2,y2,P(1,:),P(2,:),'ro')
%
%   Author : NS
%   Version: 3.0, 21 Sept. 2010
%
%   Two words about the algorithm: Most of the code is self-explanatory.
%   The only trick lies in the calculation of C1 and C2. To be brief, this
%   is essentially the two-dimensional analog of the condition that needs
%   to be satisfied by a function F(x) that has a zero in the interval
%   [a,b], namely
%           F(a)*F(b) <= 0
%   C1 and C2 exactly do this for each segment of curves 1 and 2
%   respectively. If this condition is satisfied simultaneously for two
%   segments then we know that they will cross at some point.
%   Each factor of the 'C' arrays is essentially a matrix containing
%   the numerators of the signed distances between points of one curve
%   and line segments of the other.

% Check number of inputs
narginchk(1, 2);

% Specific init depending on number of inputs
if nargin == 1
    % Compute self-intersections
    % -> Avoid the inclusion of common points
    poly2 = poly1;
    hF = @lt;
else
    % Compute intersections between distinct lines
    poly2 = varargin{1};
    hF = @le;
end

% Get coordinates of polyline vertices
x1 = poly1(:,1);
x2 = poly2(:,1)';
y1 = poly1(:,2);
y2 = poly2(:,2)';

% differentiate coordinate arrays
dx1 = diff(x1); dy1 = diff(y1);
dx2 = diff(x2); dy2 = diff(y2);

% Determine 'signed distances'
S1 = dx1 .* y1(1:end-1) - dy1 .* x1(1:end-1);
S2 = dx2 .* y2(1:end-1) - dy2 .* x2(1:end-1);

C1 = feval(hF, D(bsxfun(@times,dx1,y2) - bsxfun(@times,dy1,x2), S1), 0);
C2 = feval(hF, D((bsxfun(@times,y1,dx2) - bsxfun(@times,x1,dy2))', S2'), 0)';

% Obtain the segments where an intersection is expected
[i, j] = find(C1 & C2);

% Process case of no intersection
if isempty(i)
    pts = zeros(0, 2);
    return;
end

% Transpose and prepare for output
i=i'; dx2=dx2'; dy2=dy2'; S2 = S2';
L = dy2(j).*dx1(i) - dy1(i).*dx2(j);

% Avoid divisions by zero
i = i(L~=0);
j = j(L~=0);
L = L(L~=0);

% Solve system of eqs to get the common points
res = [dx2(j).*S1(i) - dx1(i).*S2(j), dy2(j).*S1(i) - dy1(i).*S2(j)] ./ [L L];
pts = unique(res, 'rows');

% Innre function computing a kind of cross-product
function u = D(x,y)
    u = bsxfun(@minus, x(:,1:end-1), y) .* bsxfun(@minus, x(:,2:end), y);
end

end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [E,J]=SynthMeasWatsonSHCylSingleRadTortGPD(x, protocol, fibredir, roots)
% Substrate: Impermeable cylinders with one radius in a homogeneous background.
% Orientation distribution: Watson's distribution with SH approximation
% Pulse sequence: Any
% Signal approximation: Gaussian phase distribution.
% Notes: This version estimates the hindered diffusivity from the free diffusivity
% and packing density using Szafer et al's tortuosity model for randomly
% packed cylinders.
%
% [E,J]=SynthMeasWatsonSHCylSingleRadTortGPD(x, protocol, fibredir, roots)
% returns the measurements E according to the model and the Jacobian J of the
% measurements with respect to the parameters.  The Jacobian does not
% include derivates with respect to the fibre direction.
%
% x is the list of model parameters in SI units:
% x(1) is the volume fraction of the intracellular space.
% x(2) is the free diffusivity of the material inside and outside the cylinders.
% x(3) is the radius of the cylinders.
% x(4) is the concentration parameter of the Watson's distribution.
%
% protocol is the object containing the acquisition protocol.
%
% fibredir is a unit vector along the symmetry axis of the Watson's
% distribution.  It must be in Cartesian coordinates [x y z]' with size [3 1].
%
% roots contains solutions to the Bessel function equation from function
% BesselJ_RootsCyl.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

f=x(1);
dPar=x(2);
% This version is for cylinders with regular packing.
%dPerp = dPar/((1 + f^(3/2))^2);
% This one is for randomly packed cylinders
dPerp = dPar*(1-f);
R=[x(3)];
kappa=x(4);

x_full = [f dPar dPerp R kappa];

% Call the model with no isotropic component to get the anisotropic component.
if(nargout == 1)
    E=SynthMeasWatsonSHCylSingleRadGPD(x_full, protocol, fibredir, roots);
else
    [E,J_full]=SynthMeasWatsonSHCylSingleRadGPD(x_full, protocol, fibredir, roots);
    J(:,1) = J_full(:,1) - J_full(:,3)*dPar;
    J(:,2) = J_full(:,2) + J_full(:,3)*(1-f);
    J(:,3) = J_full(:,4);
    J(:,4) = J_full(:,5);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function ConvertKlusta2Neurosuite(shank,basepath,basename,numCores)
% Converts .kwik/kwx files from Klusta into klusters-compatible
% fet,res,clu,spk files.  Works on a single shank of a recording, assumes a
% 16bit .dat and an .xml file is present in "basepath" (home folder) and
% that they are named basename.dat and basename.xml.  Also assumes that
% subdirectories in that basepath are made for each shank with names
% specified by numbers (ie 1,2,3,4..8).  In each shank folder should be
% .kwik and .kwx files made by klusta with names as follows:
% basename_sh[shankumber].kwik/kwx.  This
%
% Inputs:
% shank - the shank number (as a number, not text) equalling the name of
%           the folder under basepath with the data of interst.  Default = 1.
% basepath - directory path to the main recording folder with .dat and .xml
%           as well as shank folders made by makeProbeMapKlusta2.m (default is
%           current directory matlab is pointed to)
% basename - shared file name of .dat and .xml (default is last part of
%           current directory path, ie most immediate folder name)
% numCores  - double between 1 and inf that sets the number of cores for a
%             parloop to run (set 2-12 on SSD or RAID0 systems, set to 1
%             for data stored on a single hard drive

% Brendon Watson 2016
%
% Modified by Luke Sjulson, 2017-10

if ~exist('shank','var')
    shank = 1;
end
if ~exist('basepath','var');
    basepath = cd;
end
if ~exist('basename','var');
    [~,basename] = fileparts(cd);
elseif strcmp(basename, 'lfpfile')
    d = dir('*lfp');
    basename = d.name(1:end-4);
end
if ~exist('numCores','var');
    numCores = 1;
end

%% Start grabbing data
datpath = fullfile(basepath,[basename '.dat']);
tkwik = fullfile(basepath,num2str(shank),[basename '_sh' num2str(shank) '.kwik']);
tkwx = fullfile(basepath,num2str(shank),[basename '_sh' num2str(shank) '.kwx']);
% kwikinfo = h5info(tkwik,['/channel_groups/' num2str(shank) '/clusters/main']);
clu = h5read(tkwik,['/channel_groups/' num2str(shank) '/spikes/clusters/main']);
cluster_names = unique(clu);

try % different klusta versions have different names for this field
    totalch = h5readatt(tkwik,'/application_data/spikedetekt','nchannels');
catch
    totalch = h5readatt(tkwik,'/application_data/spikedetekt','n_channels');
end
sbefore = h5readatt(tkwik,'/application_data/spikedetekt','extract_s_before');
safter = h5readatt(tkwik,'/application_data/spikedetekt','extract_s_after');
channellist = h5readatt(tkwik,['/channel_groups/' num2str(shank)],'channel_order')+1;

%% find channels missing from spike groups becuase they were in "bad_channels.txt"
% this reflects what the xml and klusters expect but is not present
% will fill in data from missing channels with zeros
Par = LoadParameters(fullfile(basepath,[basename '.xml']));

spkgroupchannellist = Par.SpkGrps(shank).Channels + 1;
if length(spkgroupchannellist) > length(channellist)
    missingchannelidxs = find(spkgroupchannellist==setdiff(spkgroupchannellist,channellist));
else
    missingchannelidxs = [];
end

%% Getting spiketimes
spktimes = h5read(tkwik,['/channel_groups/' num2str(shank) '/spikes/time_samples']);
% spktimes = spktimes(1:10);

%% From Azahara: setting noise as 0 and MUA as 1
% when some files were corrupted I was able to still read them with this, not sure if all the corruption are like this one...
% Code that klustaviewa uses: 0 = noise, 1 = MUA, 2 = good
for ind = 1:length(cluster_names)
    kk=h5readatt(tkwik,['/channel_groups/' num2str(shank) '/clusters/main/' num2str(cluster_names(ind))],'cluster_group');
    cluster_group(ind) = kk(1,1);
    clear kk

    if cluster_group(ind) == 0 %NOISE
        clu(find(clu==cluster_names(ind))) = 0;
    elseif cluster_group(ind) == 1 %MUA
        clu(find(clu==cluster_names(ind))) = 1;
    end
end
clear ind cluster_group

%% spike extraction from dat
dat=memmapfile(datpath,'Format','int16');
tsampsperwave = (sbefore+safter);
ngroupchans = length(spkgroupchannellist);
valsperwave = tsampsperwave * ngroupchans;
wvforms_all = zeros(length(spktimes),valsperwave,'int16');
wvranges = zeros(length(spktimes),ngroupchans);
wvpowers = zeros(1,length(spktimes));

if numCores>1
    delete(gcp('nocreate')); parpool(numCores)
    parfor j=1:length(spktimes)
        % for j=1:length(spktimes)
        try
            byteList = [];
            for i = 1:length(channellist)
                ind = (double(spktimes(j))-sbefore).*totalch+1:(double(spktimes(j))+safter).*totalch;
                byteList = [byteList; ind(double(channellist(i)):totalch:end)];
            end
            wvforms = dat.data(byteList);
            %         w = dat.data((double(spktimes(j))-sbefore).*totalch+1:(double(spktimes(j))+safter).*totalch);
            %         wvforms=reshape(w,totalch,[]);
            %         %select needed channels
            %         wvforms = wvforms(channellist,:);
            % %         % detrend
            % %         wvforms = floor(detrend(double(wvforms)));
            for midx = 1:length(missingchannelidxs)%account for bad_channels that remained in spike groups in xml
                t = missingchannelidxs(midx);
                wvforms = cat(1,wvforms(1:t-1,:),zeros(1,size(wvforms,2)),wvforms(t:end,:));
            end
            % median subtract
            wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter);
        catch
            disp(['Error extracting spike at sample ' int2str(double(spktimes(j))) '. Saving as zeros']);
            disp(['Time range of that spike was: ' num2str(double(spktimes(j))-sbefore) ' to ' num2str(double(spktimes(j))+safter) ' samples'])
            wvforms = zeros(valsperwave,1);
        end
        wvforms_all(j,:)=wvforms(:);

        %some processing for fet file
        wvaswv = reshape(wvforms,ngroupchans,tsampsperwave)';
        wvranges(j,:) = range(wvaswv);
        wvpowers(j) = sum(sum(wvaswv.^2));

        %     lastpoint = tsampsperwave*ngroupchans*(j-1);
        wvforms_all(j,:)=wvforms(:);
        %     wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms;
        %     wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)')));
        if rem(j,50000) == 0
            disp([num2str(j) ' out of ' num2str(length(spktimes)) ' left'])
        end
    end

else % if using only one core, don't open a parpool
    for j=1:length(spktimes)
        % for j=1:length(spktimes)
        try
            byteList = [];
            for i = 1:length(channellist)
                ind = (double(spktimes(j))-sbefore).*totalch+1:(double(spktimes(j))+safter).*totalch;
                byteList = [byteList; ind(double(channellist(i)):totalch:end)];
            end
            wvforms = dat.data(byteList);
            %         w = dat.data((double(spktimes(j))-sbefore).*totalch+1:(double(spktimes(j))+safter).*totalch);
            %         wvforms=reshape(w,totalch,[]);
            %         %select needed channels
            %         wvforms = wvforms(channellist,:);
            % %         % detrend
            % %         wvforms = floor(detrend(double(wvforms)));
            for midx = 1:length(missingchannelidxs)%account for bad_channels that remained in spike groups in xml
                t = missingchannelidxs(midx);
                wvforms = cat(1,wvforms(1:t-1,:),zeros(1,size(wvforms,2)),wvforms(t:end,:));
            end
            % median subtract
            wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter);
        catch
            disp(['Error extracting spike at sample ' int2str(double(spktimes(j))) '. Saving as zeros']);
            disp(['Time range of that spike was: ' num2str(double(spktimes(j))-sbefore) ' to ' num2str(double(spktimes(j))+safter) ' samples'])
            wvforms = zeros(valsperwave,1);
        end
        wvforms_all(j,:)=wvforms(:);

        %some processing for fet file
        wvaswv = reshape(wvforms,ngroupchans,tsampsperwave)';
        wvranges(j,:) = range(wvaswv);
        wvpowers(j) = sum(sum(wvaswv.^2));

        %     lastpoint = tsampsperwave*ngroupchans*(j-1);
        wvforms_all(j,:)=wvforms(:);
        %     wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms;
        %     wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)')));
        if rem(j,50000) == 0
            disp([num2str(j) ' out of ' num2str(length(spktimes)) ' left'])
        end
    end

end

clear dat
wvforms_all = reshape(wvforms_all',size(wvforms_all,1)*size(wvforms_all,2),1);
wvranges = wvranges';

%% Spike features
fets = h5read(tkwx,['/channel_groups/' num2str(shank) '/features_masks']);
fets = double(squeeze(fets(1,:,:)));

% normalizing fets to be in a good numerical range for klusters (added by Luke, 2017-10)
f1 = size(fets,1);
f2 = size(fets,2);
fets = zscore(fets(:));
fets = reshape(fets, f1, f2);

fets(fets<-5) = -5; % trimming off everything >5 SDs from mean
fets(fets>5) = 5;

fets = fets.*1000;

% %mean activity per spike
% fetmeans = mean(fets,1);

%find first pcs, make means of those...
% featuresperspike = 4;%this is dumb, should read from file - not using it any more anyway.
% firstpcslist = 1:featuresperspike:size(fets,1);
% firstpcmeans = mean(fets(firstpcslist,:),1);

nfets = size(fets,1)+1;
% fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers,wvranges,double(spktimes')); % these are a bunch of extra features that many people don't find helpful.
fets = cat(1,fets,double(spktimes')); % just the PCA projections and timestamps
fets = fets';
% fets = cat(1,nfets,fets);

%% writing to clu, res, fet, spk
cluname = fullfile(basepath,[basename '.clu.' num2str(shank)]);
resname = fullfile(basepath,[basename '.res.' num2str(shank)]);
fetname = fullfile(basepath,[basename '.fet.' num2str(shank)]);
spkname = fullfile(basepath,[basename '.spk.' num2str(shank)]);

%clu
clu = cat(1,length(unique(clu)),double(clu));
fid=fopen(cluname,'w');
fprintf(fid,'%.0f\n',clu);
fclose(fid);
clear fid

%res
fid=fopen(resname,'w');
fprintf(fid,'%.0f\n',spktimes);
fclose(fid);
clear fid

%fet
SaveFetIn(fetname,fets);

%spk
fid=fopen(spkname,'w');
fwrite(fid,wvforms_all,'int16');
fclose(fid);
clear fid

disp(['Shank ' num2str(shank) ' done'])

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize)
    BufSize = inf;
end

nFeatures = size(Fet, 2);
formatstring = '%d';
for ii=2:nFeatures
    formatstring = [formatstring,'  %d'];
end
formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w');
fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)
    fprintf(outputfile,formatstring,round(Fet'));
else
    nBuf = floor(size(Fet,1)/BufSize)+1;

    for i=1:nBuf
        BufInd = [(i-1)*nBuf+1:min(i*nBuf,size(Fet,1))];
        fprintf(outputfile,formatstring,round(Fet(BufInd,:)'));
    end
end

fclose(outputfile);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [Jul1,Jul2]=TT2TCB(Jul1,Jul2,deltaT,clockLoc)
%%TT2TCB Convert from terrestrial time (TT) to barycentric coordinate
%        time (TCB) to an accuracy of nanoseconds (if deltaT is accurate)
%        using the routines from the International Astronomical Union's
%        library that do not require external ephemeris data.
%
%INPUTS: Jul1,Jul2 Two parts of a Julian date given in TT. The units of
%                  the date are days. The full date is the sum of both
%                  terms. The date is broken into two parts to provide
%                  more bits of precision. It does not matter how the date
%                  is split.
%           deltaT An optional parameter specifying the offset between TT
%                  and UT1 in seconds. If this parameter is omitted or an
%                  empty matrix is passed, then the value of the function
%                  getEOP will be used.
%         clockLoc An optional 3X1 vector specifying the location of the
%                  clock in WGS-84 ECEF Cartesian [x;y;z] coordinates with
%                  units of meters. Due to relativistic effects, clocks
%                  that are synchronized with respect to TT are not
%                  synchronized with respect to TCB. If this parameter is
%                  omitted, then a clock at the center of the Earth is
%                  used and the precision declines to microseconds.
%
%OUTPUTS:Jul1,Jul2 Two parts of a Julian date given in TCB.
%
%This function relies on a number of functions in the International
%Astronomical Union's Standards of Fundamental Astronomy library. The time
%is first converted to barycentric dynamical time (TDB) and then to
%barycentric coordinated time (TCB).
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%[Jul1,Jul2]=TT2TCB(Jul1,Jul2);
%or
%[Jul1,Jul2]=TT2TCB(Jul1,Jul2,deltaT);
%or
%[Jul1,Jul2]=TT2TCB(Jul1,Jul2,deltaT,clockLoc);
%
%Many temporal coordinate systems standards are compared in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report,
%    Naval Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016,
%    173 pages.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

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 Z = projectData(X, U, K)
%PROJECTDATA Computes the reduced data representation when projecting only
%on to the top k eigenvectors
%   Z = projectData(X, U, K) computes the projection of
%   the normalized inputs X into the reduced dimensional space spanned by
%   the first K columns of U. It returns the projected examples in Z.
%

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

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the projection of the data using only the top K
%               eigenvectors in U (first K columns).
%               For the i-th example X(i,:), the projection on to the k-th
%               eigenvector is given as follows:
%                    x = X(i, :)';
%                    projection_k = x' * U(:, k);
%

Ureduce = U(:, 1:K);

for i=1:size(X, 1)
    z = Ureduce' * X(i, :)';
    Z(i, :) = z;
end;

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

end

author: zsiciarz

repo: ml-coursera


text:
% FIT_CUBIC_BEZIER Fit a cubic bezier spline (G1 continuous) to an ordered list
% of input points in any dimension, according to "An algorithm for automatically
% fitting digitized curves" [Schneider 1990].
%
%  Inputs:
%    d  #d by dim list of points along a curve to be fit with a cubic bezier
%      spline (should probably be roughly uniformly spaced). If d(0)==d(end),
%      then will treat as a closed curve.
%    error  maximum squared distance allowed
%  Output:
%    cubics #cubics list of 4 by dim lists of cubic control points
%
% Example:
%   [P,~,C] = remove_duplicate_vertices(cell2mat(fit_cubic_bezier(X,1)),0);
%   C = reshape(C,4,[])';

author: alecjacobson

repo: gptoolbox


text:
function value = p06_param ( action, name, value )

%*****************************************************************************80
%
%% P06_PARAM gets or sets parameters for problem 6.
%
%  Discussion:
%
%    The parameter is named "M", and it represents the value of the exponent
%    in the integrand function:
%
%    Integral ( -oo < x < +oo ) x^m exp (-x*x) dx
%
%    M must be greater than -1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character ACTION, the action.
%    'S' to set the value,
%    'G' to get the value.
%
%    Input, character NAME, the parameter name.
%    'M', the exponent.
%
%    Input/output, integer VALUE, the parameter value.
%    If ACTION = 'S', then VALUE is an input quantity, and M is set to VALUE.
%    If ACTION = 'G', then VALUE is an output quantity, and VALUE is set to M.
%
  persistent m;

  if ( action == 'S' || action == 's' )

    if ( value <= -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P06_PARAM - Fatal error!\n' );
      fprintf ( 1, '  Parameter M must be greater than -1.\n' );
      error ( 'P06_PARAM - Fatal error!' );
    end

    m = value;

  elseif ( action == 'G' || action == 'g' )

    value = m;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P06_PARAM - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized value of ACTION = "%s".\n', action );
    error ( 'P06_PARAM - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function imsegs2 = bndinfo2imsegs(bndinfo, imsegs1)

tmpseg.segimage = bndinfo.wseg;
if isfield(bndinfo, 'imname')
    tmpseg.imname = bndinfo.imname;
end

stats = regionprops(bndinfo.wseg, 'Area');
tmpseg.npixels = vertcat(stats.Area);
tmpseg.nseg = numel(tmpseg.npixels);
tmpseg.imsize = bndinfo.imsize(1:2);

if exist('imsegs1', 'var') && isfield(imsegs1, 'labels')
    imsegs2 = APPtransferLabels(imsegs1, tmpseg);
else
    imsegs2 = tmpseg;
end

imsegs2 = orderfields(imsegs2);

author: Cloud-CV

repo: object-proposals


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

%%
clear; close all; clc;

%% Syntax
% |Nn=vec_normalize(N);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% icadefs() - function to read in a set of EEGLAB system-wide (i.e. lab-wide)
%             or working directory-wide constants and preferences. Change the
%             way these are defined in the master icadefs.m file (usually
%             in dir eeglab/functions/sigprocfunc) or make a custom copy of
%             the icadefs.m file in a project directory. Then, calling functions
%             that call icadefs from an EEGLAB session in that working directory
%             will read the local copy, which may set preferences different from
%             the system-wide copy.
%
% Author: Arnaud Delorme, Scott Makeig, SCCN/INC/UCSD, La Jolla, 05-20-97

% Copyright (C) 05-20-97 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% ----------------------------------------------------------------------
% ------ EEGLAB DEFINITION - YOU MAY CHANGE THE TEXT BELOW -------------
% ----------------------------------------------------------------------

EEGOPTION_PATH = ''; % if empty, the home folder of the current user is used
                     % Note that this may create problems under Windows
                     % when unicode characters are part of the user name
                     % In this case, enter the path name manually here.

YDIR  = 1;                  % positive potential up = 1; negative up = -1
                            % for most ERP plots

HZDIR = 'up';               % ascending freqs = 'up'; descending = 'down'
                            % (e.g., timef/newtimef frequency direction)

% Checking MATLAB version
tmpvers = version;
indp = find(tmpvers == '.');
if str2num(tmpvers(indp(1)+1)) >= 1, tmpvers = [ tmpvers(1:indp(1)) '0' tmpvers(indp(1)+1:end) ]; end;
indp = find(tmpvers == '.');
VERS = str2num(tmpvers(1:indp(2)-1));

% font size
tmpComputer   = computer;
tmpScreenSize = get(0, 'ScreenSize');
retinaDisplay = false;
if tmpScreenSize(3) == 1440 && ( tmpScreenSize(3) == 878 || tmpScreenSize(3) == 900 )
    retinaDisplay = true;
end;

% retinaDisplay = false; % uncoment this line if not retina display
if retinaDisplay && strcmpi(tmpComputer(1:3), 'MAC')
    W_MAIN = findobj('tag', 'EEGLAB');
    if isempty(W_MAIN)
        disp('Mac OSX retina display detected. If this is not the case uncoment line 58 of icadefs.m');
    end;
    GUI_FONTSIZE  = 18; % graphic interface font size
    AXES_FONTSIZE = 18; % Axis labels and legend font size
    TEXT_FONTSIZE = 18; % Miscellaneous font sizes
end

% Graph Definitions
DEFAULT_COLORMAP = 'jet';

if VERS < 8.04
    PLOT_LINEWIDTH   = 2;
    PLOT_LINEWIDTH_S = 1;

    % AXES FONTSIZE
    AXES_FONTSIZE   = 10;                % Axis labels and legend font size
    AXES_FONTSIZE_S = AXES_FONTSIZE - 2; % Axis labels and legend font size Small
    AXES_FONTSIZE_L = 16;                % Axis labels and legend font size Large

    % GUI FONTSIZE
    GUI_FONTSIZE    = 10;               % graphic interface font size
    GUI_FONTSIZE_S  = GUI_FONTSIZE - 2; % graphic interface font size Small
    GUI_FONTSIZE_L  = GUI_FONTSIZE + 2; % graphic interface font size Large

    % TEXT FONTSIZE
    TEXT_FONTSIZE = 10;                  % Miscellaneous font sizes
    TEXT_FONTSIZE_S = TEXT_FONTSIZE - 2; % Miscellaneous font sizes Small
    TEXT_FONTSIZE_L = TEXT_FONTSIZE + 2; % Miscellaneous font sizes Large

elseif VERS >= 8.04
    if strcmpi(tmpComputer(1:3), 'MAC')
        PLOT_LINEWIDTH   = 1;
        PLOT_LINEWIDTH_S = 0.5;

        % AXES FONTSIZE
        AXES_FONTSIZE   = 9;                 % Axis labels and legend font size
        AXES_FONTSIZE_S = AXES_FONTSIZE - 2; % Axis labels and legend font size Small
        AXES_FONTSIZE_L = 12.5;              % Axis labels and legend font size Large

        % GUI FONTSIZE
        GUI_FONTSIZE    = 12;                % graphic interface font size
        GUI_FONTSIZE_S  = GUI_FONTSIZE - 2; % graphic interface font size Small
        GUI_FONTSIZE_L  = GUI_FONTSIZE + 2; % graphic interface font size Large

        % TEXT FONTSIZE
        TEXT_FONTSIZE   = 12;                 % Miscellaneous font sizes
        TEXT_FONTSIZE_S = TEXT_FONTSIZE - 2; % Miscellaneous font sizes Small
        TEXT_FONTSIZE_L = TEXT_FONTSIZE + 4; % Miscellaneous font sizes Large
    else
        PLOT_LINEWIDTH   = 1;
        PLOT_LINEWIDTH_S = 0.5;

        % AXES FONTSIZE
        AXES_FONTSIZE   = 9;                 % Axis labels and legend font size
        AXES_FONTSIZE_S = AXES_FONTSIZE - 2; % Axis labels and legend font size Small
        AXES_FONTSIZE_L = 12.5;              % Axis labels and legend font size Large

        % GUI FONTSIZE
        GUI_FONTSIZE    = 10;                % graphic interface font size
        GUI_FONTSIZE_S  = GUI_FONTSIZE - 2; % graphic interface font size Small
        GUI_FONTSIZE_L  = GUI_FONTSIZE + 2; % graphic interface font size Large

        % TEXT FONTSIZE
        TEXT_FONTSIZE   = 10;                 % Miscellaneous font sizes
        TEXT_FONTSIZE_S = TEXT_FONTSIZE - 2; % Miscellaneous font sizes Small
        TEXT_FONTSIZE_L = TEXT_FONTSIZE + 4; % Miscellaneous font sizes Large
    end
end;

clear retinaDisplay tmpScreenSize tmpComputer tmpvers indp;

% the eeg_options.m file also countains additional options

% ----------------------------------------------------------------------
% ------------------------ END OF DEFINITIONS --------------------------
% ----------------------------------------------------------------------

% INSERT location of ica executable (LINUX ONLY) for binica.m below
eeglab_p = fileparts(which('eeglab'));
ICABINARY = fullfile(eeglab_p, 'functions', 'resources', 'ica_linux');

try
    set(0,'defaultaxesfontsize',AXES_FONTSIZE);
    set(0,'defaulttextfontsize',TEXT_FONTSIZE);
    set(0,'DefaultUicontrolFontSize',GUI_FONTSIZE);
catch
    % most likely Octave here
end;

TUTORIAL_URL = 'http://sccn.ucsd.edu/wiki/EEGLAB'; % online version
DEFAULT_SRATE = 256.0175;      % default local sampling rate (rarely used)
DEFAULT_TIMLIM = [-1000 2000]; % default local epoch limits (ms)

% Set EEGLAB figure and GUI colors
% --------------------------------
lowscreendepth = 0;
if ~exist('OCTAVE_VERSION')
    if get(0, 'screendepth') <=8 % if mono or 8-bit color
	lowscreendepth = 1;
    end;
end;
if lowscreendepth
    %fprintf('icadefs(): Setting display parameters for mono or 8-bit color\n');
    BACKCOLOR           = [1 1 1];    % Background figure color
    BACKEEGLABCOLOR     = [1 1 1];    % EEGLAB main window background
    GUIBUTTONCOLOR      = [1 1 1];    % Buttons colors in figures
    GUIPOPBUTTONCOLOR   = [1 1 1];    % Buttons colors in GUI windows
    GUIBACKCOLOR        = [1 1 1];    % GUI background color
    GUITEXTCOLOR        = [0 0 0];      % GUI foreground color for text
    PLUGINMENUCOLOR     = [.5 0 .5];  % plugin menu color

else % if full color screen
    BACKCOLOR           = [.93 .96 1];    % EEGLAB Background figure color
    BACKEEGLABCOLOR     = [.66 .76 1];    % EEGLAB main window background
    GUIBUTTONCOLOR      = BACKEEGLABCOLOR;% Buttons colors in figures
    GUIPOPBUTTONCOLOR   = BACKCOLOR;      % Buttons colors in GUI windows
    GUIBACKCOLOR        = BACKEEGLABCOLOR;% EEGLAB GUI background color <---------
    GUITEXTCOLOR        = [0 0 0.4];      % GUI foreground color for text
    PLUGINMENUCOLOR     = [.5 0 .5];      % plugin menu color
end;

% THE FOLLOWING PARAMETERS WILL BE DEPRECATED IN LATER VERSIONS
% -------------------------------------------------------------

SHRINKWARNING = 1;          % Warn user about the shrink factor in topoplot() (1/0)

MAXENVPLOTCHANS   = 264;  % maximum number of channels to plot in envproj.m
MAXPLOTDATACHANS  = 264;  % maximum number of channels to plot in dataplot.m
MAXPLOTDATAEPOCHS = 264;  % maximum number of epochs to plot in dataplot.m
MAXEEGPLOTCHANS   = 264;  % maximum number of channels to plot in eegplot.m
MAXTOPOPLOTCHANS  = 264;  % maximum number of channels to plot in topoplot.m

DEFAULT_ELOC  = 'chan.locs'; % default electrode location file for topoplot.m
DEFAULT_EPOCH = 10;       % default epoch width to plot in eegplot(s) (in sec)

SC  =  ['binica.sc'];           % Master .sc script file for binica.m
                                % MATLAB will use first such file found
                                % in its path of script directories.
                                % Copy to pwd to alter ICA defaults


text:
% "example_nufft1_reverse.m"
% This m-file is an example of applying the NUFFT method "in reverse,"
% meaning that you have a nonuniformly-sampled signal spectrum in hand
% and you want to try to "compute its inverse Fourier transform"
% to get uniform set of signal samples.  In the nomenclature of the 1999
% SIAM J Sci. Comput. paper by Nguyen and Liu, this is "Problem 1."
% (It is somewhat related to MRI reconstruction by gridding.)
% Before doing this, you should ask yourself if that is *really*
% what you want to do.  In most inverse problems, I think it is NOT
% what we should do, as argued in the T-SP paper on this method.
% Instead, we should solve the inverse problem *iteratively*, using
% the NUFFT method (and its adjoint) each iteration.
% Furthermore, the min-max optimality of the NUFFT method was only
% established for the "forward direction:" going from uniform signal
% samples to nonuniform frequency samples.
% (You can replace time and frequency in the above discussion.)
% Anyway, let's try it here and see how it works...

%
% synthesize some spectral data that is nonuniformly spaced
%
No = 2^7;				% number of frequency samples
rng(0)
om = 2*pi*sort(rand(No,1)-0.5);		% random frequency samples on [-pi,pi ]
%om = 2*pi*[-No/2:No/2-1]'/No;		% test with uniform samples

% a spectrum with periodic components, hopefully visible in other domain
X = inline('2 + 2*sin(om*10) + 4*cos(15*om) + 4*cos(20*om)', 'om');
Xm = X(om);

%
% go to the time domain the "exact" (slow) DTFT way.
% this implements essentially equation (1) in Nguyen and Liu (1999)
%
N1 = 2^6;		% # of time samples
n = [0:N1-1]'-N1/2;	% time sample locations
xd = dtft2_adj(Xm, [om 0*om], N1, 1, [N1/2 0]);

%
% now do it the fast NUFFT way
%
if 1 || ~isvar('st')	% create NUFFT structure
	J = 5;		% interpolation neighborhood
	K1 = N1*2;	% two-times oversampling
	st = nufft_init(om, N1, J, K1, N1/2, 'minmax:kb');
end

xn = nufft_adj(Xm, st);	% call ADJOINT to go "in reverse"

%
% compare slow exact to fast NUFFT
%
figure(1), clf, pl=220;
subplot(pl+1)
oo = [-200:200]'/400*2*pi;
plot(om, Xm, '.', oo, X(oo), '-')
xlabel \omega, ylabel X(\omega), title 'Spectrum and samples'

subplot(pl+2)
plot(n, real(xd), 'c.-', n, imag(xd), 'y.-')
xlabel n, ylabel x[n], title 'Exact "nonuniform FT"'

subplot(pl+3)
plot(n, real(xn), 'c.-', n, imag(xn), 'y.-')
xlabel n, ylabel x[n], title 'Fast NUFFT adjoint'

subplot(pl+4)
plot(n, real(xn-xd), 'c.-', n, imag(xn-xd), 'y.-')
xlabel n, ylabel x[n], title 'Approximation Error (very small!)'

%
% ok fine, so if you look at figure 1 you see that
% the NUFFT gives a great approximation to the "exact" formula
% but is the result *really* what you wanted?
% It seems that ideally the result should be 4 spikes with
% no background junk.  If we solved "Problem 5" in Nguyen and Liu,
% then we should get that!  This is how we have approached the
% MRI reconstruction problem, as described in ../mri.
%
% For iterative 2D version, see ../example/mri_example.m

return

%
% here is a different strategy: interpolate the nonuniform data
% onto a uniform grid, then simply take the inverse FFT.
% this didn't work so well because it is a lousy gridding method.
% todo: need to put a better gridding method here!
%
ok = [-N1/2:(N1/2-1)]'/N1*2*pi;
Xg = interp1(om, Xm, ok, 'linear', 'extrap');
xg = fftshift(ifft(fftshift(Xg)));

%xg = xg ./ nufft_sinc(n/N1).^2; % post-compensate?

figure(2), clf, pl=220;
subplot(pl+1)
plot(n, real(xg), 'c.-', n, imag(xg), 'y.-')
xlabel n, ylabel x[n], title 'Gridding "reconstruction"'

Xg = dtft1(xg, om, N1/2);
Xd = dtft1(xd/No, om, N1/2);

subplot(pl+2)
plot(om, real(Xg), '.', om, real(Xm), '-')
xlabel \omega, ylabel X(\omega), title ''

subplot(pl+3)
plot(om, real(Xd), '.', om, real(Xm), '-')
xlabel \omega, ylabel X(\omega), title ''

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% this code is inspired by VOCevaldet in the PASVAL VOC devkit
% Note: this function has been significantly optimized since ILSVRC2013
function [ap recall precision] = eval_detection(predict_file,gtruth_dir,meta_file,...
                                                eval_file,blacklist_file,optional_cache_file)
% Evaluate detection
% - predict_file: each line is a single predicted object in the
%   format
%    <image_id> <ILSVRC2014_DET_ID> <confidence> <xmin> <ymin> <xmax> <ymax>
% - gtruth_dir: a path to the directory of ground truth information,
%   e.g., ILSVRC2014_DET_bbox_val/
% - meta_file: information about the synsets
% - eval_file: list of images to evaluate on
% - blacklist_file: list of image/category pairs which aren't
%    considered in evaluation
% - optional_cache_file: to save the ground truth data and avoid
%    loading from scratch again

if nargin < 3
    meta_file = '../data/meta_det.mat';
end
if nargin < 4
    eval_file = '../data/det_lists/val.txt';
end
if nargin < 5
    blacklist_file = '../data/ILSVRC2014_det_validation_blacklist.txt';
end
if nargin < 6
    optional_cache_file = '';
end

defaultIOUthr = 0.5;
pixelTolerance = 10;

load(meta_file);
hash = make_hash(synsets);

bLoadXML = true;
if ~isempty(optional_cache_file) && exist(optional_cache_file,'file')
    fprintf('eval_detection :: loading cached ground truth\n');
    t = tic;
    load(optional_cache_file);
    fprintf('eval_detection :: loading cached ground truth took %0.1f seconds\n',toc(t));
    if exist('gt_obj_img_ids','var')
        fprintf(['eval_detection :: loaded cache from 2014, ' ...
                 'recomputing\n']);
    else
        bLoadXML = false;
    end
end
if bLoadXML
    fprintf('eval_detection :: loading ground truth\n');
    t = tic;

    [img_basenames gt_img_ids] = textread(eval_file,'%s %d');

    num_imgs = length(img_basenames);
    gt_obj_labels = cell(1,num_imgs);
    gt_obj_bboxes = cell(1,num_imgs);
    gt_obj_thr = cell(1,num_imgs);
    num_pos_per_class = [];
    tic
    for i=1:num_imgs
        if toc > 60
            fprintf('              :: on %i of %i\n',...
                    i,num_imgs);
            tic;
        end
        rec = VOCreadxml(sprintf('%s/%s.xml',gtruth_dir, ...
                                 img_basenames{i}));
        if ~isfield(rec.annotation,'object')
            continue;
        end
        for j=1:length(rec.annotation.object)
            obj = rec.annotation.object(j);
            c = get_class2node(hash, obj.name);
            gt_obj_labels{i}(j) = c;
            if length(num_pos_per_class) < c
                num_pos_per_class(c) = 1;
            else
                num_pos_per_class(c) = num_pos_per_class(c) + 1;
            end
            b = obj.bndbox;
            bb = str2double({b.xmin b.ymin b.xmax b.ymax});
            gt_obj_bboxes{i}(:,j) = bb;
        end

        bb = gt_obj_bboxes{i};
        gt_w = bb(4,:)-bb(2,:)+1;
        gt_h = bb(3,:)-bb(1,:)+1;
        thr = (gt_w.*gt_h)./((gt_w+pixelTolerance).*(gt_h+pixelTolerance));
        gt_obj_thr{i} = min(defaultIOUthr,thr);
    end
    fprintf('eval_detection :: loading ground truth took %0.1f seconds\n',toc(t));

    if ~isempty(optional_cache_file)
        fprintf('eval_detection :: saving cache in %s\n',optional_cache_file);
        save(optional_cache_file,'gt_img_ids','gt_obj_labels',...
             'gt_obj_bboxes','gt_obj_thr','num_pos_per_class');
    end
end

blacklist_img_id = [];
blacklist_label = [];
if ~isempty(blacklist_file) && exist(blacklist_file,'file')
    [blacklist_img_id wnid] = textread(blacklist_file,'%d %s');
    blacklist_label = zeros(length(wnid),1);
    for i=1:length(wnid)
        blacklist_label(i) = get_class2node(hash,wnid{i});
    end
    fprintf('eval_detection :: blacklisted %i image/object pairs\n',length(blacklist_label));
else
    fprintf('eval_detection :: no blacklist\n');
end

fprintf('eval_detection :: loading predictions\n');
t = tic;
[img_ids obj_labels obj_confs xmin ymin xmax ymax] = ...
        textread(predict_file,'%d %d %f %f %f %f %f');
obj_bboxes = [xmin ymin xmax ymax]';
fprintf('eval_detection :: loading predictions took %0.1f seconds\n',toc(t));

fprintf('eval_detection :: sorting predictions\n');
t = tic;
[img_ids ind] = sort(img_ids);
obj_confs = obj_confs(ind);
obj_labels = obj_labels(ind);
obj_bboxes = obj_bboxes(:,ind);

num_imgs = max(max(gt_img_ids),max(img_ids));
obj_labels_cell = cell(1,num_imgs);
obj_confs_cell = cell(1,num_imgs);
obj_bboxes_cell = cell(1,num_imgs);
start_i = 1;
id = img_ids(1);
tic
for i=1:length(img_ids)
    if toc > 60
        fprintf('               :: on %0.2fM of %0.2fM\n',...
                i/10^6,length(img_ids)/10^6);
        tic
    end
    if (i == length(img_ids)) || (img_ids(i+1) ~= id)
        % i is the last element of this group
        obj_labels_cell{id} = obj_labels(start_i:i)';
        obj_confs_cell{id} = obj_confs(start_i:i)';
        obj_bboxes_cell{id} = obj_bboxes(:,start_i:i);
        if i < length(img_ids)
            % start next group
            id = img_ids(i+1);
            start_i = i+1;
        end
    end
end

for i=1:num_imgs
    [obj_confs_cell{i} ind] = sort(obj_confs_cell{i},'descend');
    obj_labels_cell{i} = obj_labels_cell{i}(ind);
    obj_bboxes_cell{i} = obj_bboxes_cell{i}(:,ind);
end
tp_cell = cell(1,num_imgs);
fp_cell = cell(1,num_imgs);

fprintf('eval_detection :: sorting predictions took %0.1f seconds\n', ...
        toc(t));

fprintf('eval_detection :: accumulating\n');

num_classes = length(num_pos_per_class);

t = tic;
tic;
% iterate over images
for i=1:length(gt_img_ids)
    if toc > 60
        fprintf('               :: on %i of %i\n',...
                i,length(gt_img_ids));
        tic;
    end

    id = gt_img_ids(i);
    gt_labels = gt_obj_labels{i};
    gt_bboxes = gt_obj_bboxes{i};
    gt_thr = gt_obj_thr{i};
    num_gt_obj = length(gt_labels);
    gt_detected = zeros(1,num_gt_obj);

    bSameImg = blacklist_img_id == id;
    blacklisted_obj = blacklist_label(bSameImg);

    labels = obj_labels_cell{id};
    bboxes = obj_bboxes_cell{id};

    num_obj = length(labels);
    tp = zeros(1,num_obj);
    fp = zeros(1,num_obj);
    for j=1:num_obj
        if any(labels(j) == blacklisted_obj)
            continue; % just ignore this detection
        end
        bb = bboxes(:,j);
        ovmax = -inf;
        kmax = -1;
        for k=1:num_gt_obj
            if labels(j) ~= gt_labels(k)
                continue;
            end
            if gt_detected(k) > 0
                continue;
            end
            bbgt = gt_bboxes(:,k);
            bi=[max(bb(1),bbgt(1)) ; max(bb(2),bbgt(2)) ; min(bb(3),bbgt(3)) ; min(bb(4),bbgt(4))];
            iw=bi(3)-bi(1)+1;
            ih=bi(4)-bi(2)+1;
            if iw>0 & ih>0
                % compute overlap as area of intersection / area of union
                ua=(bb(3)-bb(1)+1)*(bb(4)-bb(2)+1)+...
                   (bbgt(3)-bbgt(1)+1)*(bbgt(4)-bbgt(2)+1)-...
                   iw*ih;
                ov=iw*ih/ua;
                % makes sure that this object is detected according
                % to its individual threshold
                if ov >= gt_thr(k) && ov > ovmax
                    ovmax=ov;
                    kmax=k;
                end
            end
        end
        if kmax > 0
            tp(j) = 1;
            gt_detected(kmax) = 1;
        else
            fp(j) = 1;
        end
    end

    % put back into global vector
    tp_cell{id} = tp;
    fp_cell{id} = fp;

    for k=1:num_gt_obj
        label = gt_labels(k);
        % remove blacklisted objects from consideration as positive examples
        if any(label == blacklisted_obj)
            num_pos_per_class(label) = num_pos_per_class(label)-1;
        end
    end
end
fprintf('eval_detection :: accumulating took %0.1f seconds\n', ...
        toc(t));

fprintf('eval_detection :: computing ap\n');
t = tic;
tp_all = [tp_cell{:}];
fp_all = [fp_cell{:}];
obj_labels = [obj_labels_cell{:}];
confs = [obj_confs_cell{:}];

[confs ind] = sort(confs,'descend');
tp_all = tp_all(ind);
fp_all = fp_all(ind);
obj_labels = obj_labels(ind);

for c=1:num_classes
    % compute precision/recall
    tp = cumsum(tp_all(obj_labels==c));
    fp = cumsum(fp_all(obj_labels==c));
    recall{c}=(tp/num_pos_per_class(c))';
    precision{c}=(tp./(fp+tp))';
    ap(c) =VOCap(recall{c},precision{c});
end

fprintf('eval_detection :: computing ap took %0.1f seconds\n', ...
        toc(t));

author: zhoubolei

repo: CAM

save_path: github-repos/MATLAB/zhoubolei-CAM


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2019.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtFfmHelmholtzBWneu.m                        |
%|    #    |   VERSION    : 0.61                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : Solve neumann scatering problem with          |
%|  `---'  |                Brackage-Werner formulation                   |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
N   = 1e3
tol = 1e-3
typ = 'P1'
gss = 3
X0  = [0 0 -1]

% Spherical mesh
sphere = mshSphere(N,1);
sigma  = dom(sphere,gss);
figure
plot(sphere)
axis equal

% Radiative mesh
square     = mshSquare(5*N,[5 5]);
square.vtx = [square.vtx(:,1) zeros(size(square.vtx,1),1) square.vtx(:,2)];
hold on
plot(square)

% Frequency adjusted to maximum esge size
stp = sphere.stp;
k   = 1/stp(2)
f   = (k*340)/(2*pi)

% Incident wave
PW         = @(X) exp(1i*k*X*X0');
gradxPW{1} = @(X) 1i*k*X0(1) .* PW(X);
gradxPW{2} = @(X) 1i*k*X0(2) .* PW(X);
gradxPW{3} = @(X) 1i*k*X0(3) .* PW(X);

% Incident wave representation
plot(sphere,real(PW(sphere.vtx)))
plot(square,real(PW(square.vtx)))
title('Incident wave')
xlabel('X');   ylabel('Y');   zlabel('Z');
hold off
view(0,10)
% camlight
% material dull
% lighting phong

%%% PREPARE OPERATORS
disp('~~~~~~~~~~~~~ PREPARE OPERATORS ~~~~~~~~~~~~~')

% Green kernel function --> G(x,y) = exp(ik|x-y|)/|x-y|
Gxy      = '[exp(ikr)/r]';
gradyGxy = {'grady[exp(ikr)/r]1','grady[exp(ikr)/r]2','grady[exp(ikr)/r]3'};

% Finite elements
u = fem(sphere,typ);
v = fem(sphere,typ);

% Coupling coeff
beta = 1i*k;

% Finite element mass matrix --> \int_Sx psi(x)' psi(x) dx
Id = integral(sigma,u,v);

% Finite element boundary operator -->
% k^2 * \int_Sx \int_Sy n.psi(x) G(x,y) n.psi(y) dx dy
% - \int_Sx \int_Sy nxgrad(psi(x)) G(x,y) nxgrad(psi(y)) dx dy
tic
Hr = 1/(4*pi) .* (k^2 * regularize(sigma,sigma,ntimes(u),'[1/r]',ntimes(v)) ...
    - regularize(sigma,sigma,nxgrad(u),'[1/r]',nxgrad(v)));
H  = 1/(4*pi) .* (k^2 * integral(sigma,sigma,ntimes(u),Gxy,k,ntimes(v),tol) ...
    - integral(sigma,sigma,nxgrad(u),Gxy,k,nxgrad(v),tol)) + Hr;
toc

% Finite element boundary operator --> \int_Sx \int_Sy psi(x)' dny G(x,y) psi(y) dx dy
tic
Dtr = 1/(4*pi) .* regularize(sigma,sigma,u,'grady[1/r]',ntimes(v)).';
Dt  = 1/(4*pi) .* integral(sigma,sigma,u,gradyGxy,k,ntimes(v),tol).' + Dtr;
toc

% Final operator [1i*k*beta*(-Id/2 + Dt) - H]
tic
LHS  = beta.*(-0.5*Id + Dt) - H;
toc

% Finite element incident wave trace --> \int_Sx psi(x) dnx(pw(x)) dx
RHS = - integral(sigma,ntimes(u),gradxPW);

%%% SOLVE LINEAR PROBLEM
disp('~~~~~~~~~~~~~ SOLVE LINEAR PROBLEM ~~~~~~~~~~~~~')

% Preconditionneur ILU
tic
[L,U] = ilu(beta.*(-0.5*Id + Dtr) - Hr);
toc

% Solve linear system : [H + 1i*k*beta*(Id/2 - Dt)] = -dnP0
tic
mu = mgcr(@(V) LHS*V,RHS,[],tol,100,L,U);
toc

% Jump for derivative
lambda = beta * mu;

%%% INFINITE SOLUTION
disp('~~~~~~~~~~~~~ INFINITE RADIATION ~~~~~~~~~~~~~')

% Plane waves direction
theta = 2*pi/1e3 .* (1:1e3)';
nu    = [sin(theta),zeros(size(theta)),cos(theta)];

% Green kernel function
Ginf      = '[exp(-ikxy)]';
gradxGinf = {'gradx[exp(-ikxy)]1','gradx[exp(-ikxy)]2','gradx[exp(-ikxy)]3'};

% Finite element infinite operators
Sinf = 1/(4*pi) .* integral(nu,sigma,Ginf,k,v,tol);
Dinf = 1/(4*pi) .* integral(nu,sigma,gradxGinf,k,ntimes(v),tol);

% Finite element radiation
sol = Sinf*lambda - Dinf*mu;

% Analytical solution
ref = sphereHelmholtz('inf','neu',1,k,nu);
norm(ref-sol,2)/norm(ref,2)
norm(ref-sol,'inf')/norm(ref,'inf')

% Graphical representation
figure
plot(theta,log(abs(sol)),'b',theta,log(abs(ref)),'--r')

%%% DOMAIN SOLUTION
disp('~~~~~~~~~~~~~ RADIATION ~~~~~~~~~~~~~')

% Finite element mass matrix --> \int_Sx psi(x)' psi(x) dx
Id = integral(sigma,u,v);

% Finite element boundary operator --> \int_Sx \int_Sy psi(x)' G(x,y) psi(y) dx dy
tic
Sbnd = 1/(4*pi) .* (integral(sigma,sigma,u,Gxy,k,v,tol) + ...
    regularize(sigma,sigma,u,'[1/r]',v));
toc

% Finite element boundary operator --> \int_Sx \int_Sy psi(x)' dny G(x,y) psi(y) dx dy
tic
Dbnd = 1/(4*pi) .* (integral(sigma,sigma,u,gradyGxy,k,ntimes(v),tol) + ...
    regularize(sigma,sigma,u,'grady[1/r]',ntimes(v)));
toc

% Boundary solution
Psca = Id\(Sbnd*lambda - (0.5*Id*mu + Dbnd*mu));
Pinc = PW(u.dof);
Pbnd = Pinc + Psca;

% Finite element radiative operator --> \int_Sy G(x,y) psi(y) dy
tic
Sdom = 1/(4*pi) .* (integral(square.vtx,sigma,Gxy,k,v,tol) + ...
    regularize(square.vtx,sigma,'[1/r]',v));
toc

% Finite element radiative operator --> \int_Sx \int_Sy psi(x)' grady(G(x,y)) ny.psi(y) dx dy
tic
Ddom = 1/(4*pi) .* ( integral(square.vtx,sigma,gradyGxy,k,ntimes(v),tol) + ...
    regularize(square.vtx,sigma,'grady[1/r]',ntimes(v)) );
toc

% Domain solution
Psca = Sdom*lambda - Ddom*mu;
Pinc = PW(square.vtx);
Pdom = Pinc + Psca;

% Annulation sphere interieure
r             = sqrt(sum(square.vtx.^2,2));
Pdom(r<=1.01) = Pinc(r<=1.01);

% Graphical representation
figure
plot(sphere,abs(Pbnd))
axis equal;
hold on
plot(square,abs(Pdom))
title('Total field solution')
colorbar
hold off
view(0,10)

%%% ANAYTICAL SOLUTIONS FOR COMPARISONS
% Analytical solution
Pbnd = sphereHelmholtz('dom','neu',1,k,1.001*sphere.vtx) + PW(sphere.vtx);
Pdom = sphereHelmholtz('dom','neu',1,k,square.vtx) + PW(square.vtx);

% Solution representation
figure
plot(sphere,abs(Pbnd))
axis equal;
hold on
plot(square,abs(Pdom))
title('Analytical solution')
colorbar
hold off
view(0,10)

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function y = normalizeAngle360(u)
% Normalize angle in degrees to [0 360]
y = mod(u, 360);
end

author: TUMFTM


text:
function sout = cvx_profile( flag )

% CVX_PROFILE	CVX-specific profiler control.
%    This is a function used for internal CVX development to help determine
%    performance limits within the CVX code itself, by turning off the profiler
%    when the solver is being called. End users will likely not find this
%    function to be useful.

global cvx___
cvx_global
s = cvx___.profile;
if nargin == 1,
    if isnumeric(flag) || islogical(flag),
        ns = double(flag) ~= 0;
    elseif ischar(flag) && size(flag,1) == 1,
        switch lower(flag),
            case 'true',
                ns = true;
            case 'false',
                ns = false;
            otherwise,
                error( 'String arugment must be ''true'' or ''false''.' );
        end
    else
        error( 'Argument must be a numeric scalar or a string.' );
    end
    if cvx___.profile ~= ns,
        cvx___.profile = ns;
        stat = profile('status');
        if ns && ~isempty( cvx___.problems ) && ~isequal( stat.ProfilerStatus, 'on' ),
            profile on
        elseif ~ns && isequal( stat.ProfilerStatus, 'on' ),
            profile off
        end
    end
end
if nargin == 0 || nargout > 0,
    sout = s;
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:
function [x, infos] = sparse_nmf(V, rank, in_options)
% Sparse nonnegative matrix factorization (sparseNMF)
%
% The problem of interest is defined as
%
%       min D(V||W*H) + lambda * sum(H(:)),
%       where
%       {V, W, H} > 0.
%
%       L1-based sparsity constraint on H.
%       Normalizes W column-wise.
%
% Given a non-negative matrix V, factorized non-negative matrices {W, H} are calculated.
%
%
% Inputs:
%       V           : (m x n) non-negative matrix to factorize
%       rank        : rank
%       in_options
%
%
% Output:
%       x           : non-negative matrix solution, i.e., x.W: (m x rank), x.H: (rank x n)
%       infos       : log information
%           epoch   : iteration nuber
%           cost    : objective function value
%           optgap  : optimality gap
%           time    : elapsed time
%           grad_calc_count : number of sampled data elements (gradient calculations)
%
% Reference:
%
%
% This file is part of NMFLibrary.
%
% Created by Patrik Hoyer, 2006 (and modified by Silja Polvi-Huttunen,
% University of Helsinki, Finland, 2014)
%
% Modified by H.Kasai on Jul. 23, 2018
%
% Change log:
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jul. 14, 2022 (Hiroyuki Kasai): Fixed algorithm.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options.lambda = 0;   % regularizer for sparsity
    local_options.cost  = 'euc'; % 'euc' or 'kl-div'

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize
    method_name = 'sparseNMF';
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, [], epoch, grad_calc_count, 0);
    % store additionally different cost
    reg_val = options.lambda*sum(sum(H));
    f_val_total = f_val + reg_val;
    infos.cost_reg = reg_val;
    infos.cost_total = f_val_total;

    if options.verbose > 1
        fprintf('sparseNMF: Epoch = 0000, cost = %.16e, optgap = %.4e\n', f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % update H with MU
        %H = (H.*(W'*(V./(W*H))))/(1+alpha);
        VC = V./(W*H + 1e-9);
        VC(V==0 & W*H==0) = 1+1e-9;
        H = (H.*(W'*VC))/(1+options.lambda);

        % update W by Lee and Seung's divergence step
        %W = W.*((V./(W*H))*H')./(ones(vdim,1)*sum(H'));
        VC = V./(W*H + 1e-9);
        VC(V==0 & W*H==0) = 1+1e-9;
        W = W.*(VC*H')./(ones(m,1)*sum(H,2)');

        % Liu, Zheng, and Lu add this normalization step
        W = W./(ones(m,1)*sum(W));

        % measure elapsed time
        elapsed_time = toc(start_time);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % update epoch
        epoch = epoch + 1;

        % store info
        [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, infos, epoch, grad_calc_count, elapsed_time);
        % store additionally different cost
        reg_val = options.lambda*sum(sum(H));
        f_val_total = f_val + reg_val;
        infos.cost_reg = [infos.cost_reg reg_val];
        infos.cost_total = [infos.cost_total f_val_total];

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.H = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function val=intSinPow(u,n)
%%INTSINPOWER Evaluate the integral of sin(u)^n du.  A definite integral
%             can be evaluated, or an indefinite integral (with a
%             particular additive constant).
%
%INPUTS: u A 2XN (for definite integral) or a 1XN (for indefinite
%          integrals) set of N points. For definite integrals, u(1,:) are
%          the real lower bounds and u(2,:) are the real upper bounds. For
%          indefinite integrals, the integral is evaluated at the points in
%          u. The values in u should be real.
%        n The positive integer exponent of the sine term.
%
%OUTPUTS: val The 1XN set of values of the integral of sin(u)^n.
%
%This function simply implements the recursion that arises from integration
%by parts from basic calculus, as are given in the tables in the back of
%[1].
%
%REFERENCES:
%[1] J. Stewart, Calculus, 7th ed. Belmont, CA: Brooks/Cole, 2012.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(u,1);

if(isempty(u))
   val=[];
   return;
end

if(numDim==1)%An indefinite integral
    val=indefIntSinPow(u,n);
else%A definite integral
    val=indefIntSinPow(u(2,:),n)-indefIntSinPow(u(1,:),n);
end
end

function val=indefIntSinPow(u,n)

%Special cases for n=0 and n=1.
if(n==0)
    val=u;
    return;
elseif(n==1)
    val=-cos(u);
    return;
end

%If here, n>=1
sinVal=sin(u);
cosVal=cos(u);

val=0;
if(mod(n,2)==0)%If n is even
    endVal=2;
else
    endVal=3;
end

coeffProd=1;
for k=n:-2:endVal
    val=val-coeffProd*(1/k)*sinVal.^(k-1).*cosVal;
    coeffProd=coeffProd*(k-1)/k;
end

if(endVal==2)
    %The final integral is over sin^0
    val=val+coeffProd*u;
else
    %The final integral is over sin^1
    val=val-coeffProd*cosVal;
end

end

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


text:
function z = daxpy ( n, sa, x, incx, y, incy )

%*****************************************************************************80
%
%% DAXPY adds a constant times one vector to another.
%
%  Modified:
%
%    25 February 2004
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Lawson, Hanson, Kincaid, Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real SA, the multiplier.
%
%    Input, real X(*), the vector to be scaled and added to Y.
%
%    Input, integer INCX, the increment between successive entries of X.
%
%    Input, real Y(*), the vector to which a multiple of X is to be added.
%
%    Input, integer INCY, the increment between successive entries of Y.
%
%    Output, real Z(*), the vector Y(*) + SA * X(*).
%
  if ( n <= 0 )

    z = [];

  elseif ( sa == 0.0 )

    ny = length ( y );
    z(1:ny) = y(1:ny);

  elseif ( incx == 1 & incy == 1 )

    z(1:n) = y(1:n) + sa * x(1:n);

  else

    ny = length ( y );
    z(1:ny) = y(1:ny);

    if ( 0 <= incx )
      ix = 1;
    else
      ix = ( - n + 1 ) * incx + 1;
    end if

    if ( 0 <= incy  )
      iy = 1;
    else
      iy = ( - n + 1 ) * incy + 1;
    end

    for i = 1 : n
      z(iy) = z(iy) + sa * x(ix);
      ix = ix + incx;
      iy = iy + incy;
    end

  end

  return
end


text:
%% FUNCTION LSSMTC
%   permute labels of L2 match L1 as good as possible
%
%% LICENSE
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%   Copyright (C) 2011 - 2012 Quanquan Gu, Jiayu Zhou, and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 17, 2012.

function [newL2] = bestMap(L1,L2)
%bestmap: permute labels of L2 match L1 as good as possible
%   [newL2] = bestMap(L1,L2);

%===========
L1 = L1(:);
L2 = L2(:);
if size(L1) ~= size(L2)
    error('size(L1) must == size(L2)');
end
L1 = L1 - min(L1) + 1;      %   min (L1) <- 1;
L2 = L2 - min(L2) + 1;      %   min (L2) <- 1;
%===========    make bipartition graph  ============
nClass = max(max(L1), max(L2));
G = zeros(nClass);
for i=1:nClass
    for j=1:nClass
        G(i,j) = length(find(L1 == i & L2 == j));
    end
end
%===========    assign with hungarian method    ======
[c,t] = hungarian(-G);
newL2 = zeros(nClass,1);
for i=1:nClass
    newL2(L2 == i) = c(i);
end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function mfcc = mfccFeature(config, x)
%MFCCFEATURE    MFCC feature extraction.
%   Format: mfcc = mfccFeature(config, x)
%   Inputs:
%       config: Configuration for MFCC features.
%       x:      Input signal.
%   Output:
%       mfcc:   MFCC features, one column per frame.

    X = abs(spectrogram(config, x));
    frames = size(X,2);
    mfcc = zeros(config.totalDims, frames);

    mfcc(1:config.coefs, :) = config.D * log(config.M * X);
    if (config.energy == true)
        energy = zeros(1,frames);
        for f = 1:frames
            start = (f-1) * config.frameShift + 1;
            finish = start - 1 + config.frameLen;
            frame = x(start:finish) .* config.window;
            energy(f) = sum(frame.^2);
        end
        mfcc(config.staticDims, :) = log(energy);
    end

    % Cepstral mean subtraction or normalization
    switch upper(config.normalization)
        case 'CMS'
            for u = 1:config.staticDims
                mfcc(u,:) = (mfcc(u,:) - mean(mfcc(u,:)));
            end
        case 'CMN'
            for u = 1:config.staticDims
                mfcc(u,:) = (mfcc(u,:) - mean(mfcc(u,:))) / std(mfcc(u,:));
            end
    end

    % Differentials
    zeroCol = zeros(config.staticDims, 1);
    diff = mfcc(1:config.staticDims, :);
    for d = 1:config.diffs
        diff = [diff(:, 2:end), zeroCol] - [zeroCol, diff(:, 1:end-1)];
        mfcc(d * config.staticDims + 1 : (d+1) * config.staticDims, :) = diff;
    end
end

author: MaigoAkisame


text:
% Computes a rotation of the PC factors consistent with a recursive assumption with the policy instrument
%
% Syntax:
%
% Fr = facrot(F,Ffast,Fslow)
%
% where:      F:      Unrestricted PC estimates (from all the dataset)
%             Ffast:  Factors assumed to be fast moving (e.g. policy instrument)
%             Fslow:  Proxy of the slow moving factors
%
%
% Bernanke, Boivin and Eliasz (2002)
% 12/17/02

function Fr = favar_facrot(F,Ffast,Fslow)

k1=size(Ffast,2);

b=bear.favar_olssvd(F,[ones(size(Ffast,1),1) Ffast Fslow]);
Fr = F - Ffast*b(2:k1+1,:);

author: european-central-bank

repo: BEAR-toolbox


text:
classdef prtClassMilVbDpGmm < prtClass
    % prtClassMilVbDpGmm VBDPGMM MIL Classifier
    %
    %  A variational bayes classifier for multiple instance data.  There
    %  are a number of parameters available for pruning/training.
    %  Performnace should (in theory) be somewhat invariant to these.  See
    %  the text of Manandhar et al., for descriptions.
    %
    %  You may need to reduce the dimensionality of the data being used via
    %  PCA, PLS, or some other pre-processing prior to using
    %  prtClassMilVbDpGmm.  We also recommend making sure the data is
    %  approximately zero-mean, unit-variance in each of the columns.
    %
    %  dsTrain = prtDataGenMultipleInstance;
    %  dsTest = prtDataGenMultipleInstance;
    %  class = prtClassMilVbDpGmm;
    %  class = class.train(dsTrain);
    %  yOutTrain = class.run(dsTrain);
    %  yOutTest = class.run(dsTest);
    %
    %  [pfTrain,pdTrain] = prtScoreRoc(yOutTrain);
    %  [pfTest,pdTest] = prtScoreRoc(yOutTest);
    %  subplot(1,1,1);
    %  plot(pfTrain,pdTrain,pfTest,pdTest);
    %

    properties (SetAccess=private)
        name = 'MilVbDpGmm'
        nameAbbreviation = 'MilVbDpGmm'
        isNativeMary = false;  % False
    end

    properties
        maxVbIter = 300;

        K0 = 5;
        K1 = 5;
        gamma_0_1_m0 = 1;
        gamma_0_2_m0 = .001;
        gamma_0_1_m1 = 1;
        gamma_0_2_m1 = .001;
        alpha = [10 10];

        rvH1
        rvH0
    end

    properties (Hidden)

        initialFit = [];
        %Dependent on "d"
        v_0_m
        v_1_m
        beta_1_m
        beta_0_m

        rho_0_1
        rho_0_0
        Phi_0_1
        Phi_0_0
    end

	methods
		function self = prtClassMilVbDpGmm(varargin)

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

            self.classTrain = 'prtDataSetClassMultipleInstance';
            self.classRun = 'prtDataSetClassMultipleInstance';
            self.classRunRetained = false;
        end

    end

    methods (Access=protected, Hidden = true)
		function self = trainAction(self,dsMil)

            milStruct = dsMil.data;
            x = cat(1,milStruct.data);

            bagLabels = dsMil.targets;
            expandedBagLabels = dsMil.expandedTargets;

            bagIndices = dsMil.getBagInds;
            uniqueBagIndices = unique(bagIndices);

            d = size(x,2);

            self.v_0_m = d+2;
            self.v_1_m = d+2;
            self.beta_1_m = d;
            self.beta_0_m = d;
            self.rho_0_0 = zeros(1,d);
            self.rho_0_1 = zeros(1,d);

            self.Phi_0_0 = eye(d)*d;
            self.Phi_0_1 = eye(d)*d;

%             params = struct('K0',5,'K1',5,'gamma_0_1_m0',1,'gamma_0_2_m0',.001,'gamma_0_1_m1',1,'gamma_0_2_m1',.001,...
%                 'v_0_m',d+2,'v_1_m',d+2,'beta_1_m',d,'beta_0_m',d, ...
%                 'rho_0_1',zeros(1,d),'rho_0_0',zeros(1,d),'Phi_0_1',eye(d)*d,'Phi_0_0',eye(d)*d,'alpha',[4 1]*10000);

            [rvH0,rvH1,xH0,xH1] = initialize(self,{milStruct.data},bagLabels);
            %self.initialFit = self.calculateDataLogLikelihood(self,xH0,xH1);
            ll0 = sum(rvH0.logPdf(xH0));
            ll1 = rvH0.pdf(xH1).*self.alpha(1)/sum(self.alpha) + rvH1.pdf(xH1).*self.alpha(2)/sum(self.alpha);
            self.initialFit = ll0 + sum(log(ll1));

            [p0,pp] = rvH0.pdf(x);
            phi_0_i = bsxfun(@rdivide,pp,sum(pp,2));

            [p1,pp] = rvH1.pdf(x);
            phi_1_i = bsxfun(@rdivide,pp,sum(pp,2));
            phi_1_i(expandedBagLabels == 0,:) = 0;
            p1(expandedBagLabels ==0) = 0;

            phi_1_M = p1./(p0 + p1);
            phi_0_M = 1-phi_1_M;

            locEps = 1e-6;

            for vbIter = 1:self.maxVbIter
                %Equations 67 and 68; note, these are different from 20-21 (!)
                %     gamma_i_1_m0 = sum(bsxfun(@times,phi_0_M,phi_0_i)) + self.gamma_0_1_m0;
                %     gamma_i_2_m0 = cat(2,0,cumsum(gamma_i_1_m0(1:end-1),2)) + self.gamma_0_2_m0;
                %
                %     gamma_i_1_m1 = sum(bsxfun(@times,phi_1_M,phi_1_i)) + self.gamma_0_1_m1;
                %     gamma_i_2_m1 = cat(2,0,cumsum(gamma_i_1_m1(1:end-1),2)) + self.gamma_0_2_m1;

                %Equations 20-21
                gamma_i_1_m0 = sum(bsxfun(@times,phi_0_M,phi_0_i)) + self.gamma_0_1_m0;
                gamma_i_2_m0 = fliplr(cumsum(fliplr(gamma_i_1_m0))) - gamma_i_1_m0 + self.gamma_0_2_m0;

                gamma_i_1_m1 = sum(bsxfun(@times,phi_1_M,phi_1_i)) + self.gamma_0_1_m1;
                gamma_i_2_m1 = fliplr(cumsum(fliplr(gamma_i_1_m1))) - gamma_i_1_m1 + self.gamma_0_2_m1;

                %Update priors... H0 Gaussian; Equations 77-79 (H0)
                respMat = bsxfun(@times,phi_0_M,phi_0_i);
                nBar0 = sum(respMat)+locEps;  %avoid 1./0 problems
                for cluster = 1:self.K0
                    mu_i_0(cluster,:) = 1/nBar0(cluster)*sum(bsxfun(@times,respMat(:,cluster),x));

                    xx = bsxfun(@times,sqrt(respMat(:,cluster)),bsxfun(@minus,x,mu_i_0(cluster,:)));
                    c = 1/nBar0(cluster)*xx'*xx;
                    cov_i_0(cluster,:) = c(:);
                end

                %Update priors... H1 Gaussian; Equations 77-79
                respMat = bsxfun(@times,phi_1_M,phi_1_i);
                nBar1 = sum(respMat)+locEps;
                for cluster = 1:self.K1
                    mu_i_1(cluster,:) = 1/nBar1(cluster)*sum(bsxfun(@times,respMat(:,cluster),x));

                    xx = bsxfun(@times,sqrt(respMat(:,cluster)),bsxfun(@minus,x,mu_i_1(cluster,:)));
                    c = 1/nBar1(cluster)*xx'*xx;
                    cov_i_1(cluster,:) = c(:);
                end

                %Updating hyper-parameters; eqns 81-84
                nu_0 = self.v_0_m + nBar0;
                nu_1 = self.v_1_m + nBar1;

                beta_0 = self.beta_0_m + nBar0;
                beta_1 = self.beta_1_m + nBar1;

                for cluster = 1:self.K0
                    rho_0(cluster,:) = self.beta_0_m*self.rho_0_0 + nBar0(cluster)*mu_i_0(cluster,:);
                    rho_0(cluster,:) = rho_0(cluster,:)./beta_0(cluster);

                    %See equation 84; Note: this is different from Equation 84 - 84 is
                    %wrong
                    %covPart = nBar0(cluster) * reshape(cov_i_0(cluster,:),d,d) * ((mu_i_0(cluster,:) - rho_0(cluster,:))'*(mu_i_0(cluster,:) - rho_0(cluster,:)));
                    covPart = nBar0(cluster) * self.beta_0_m * ((mu_i_0(cluster,:) - self.rho_0_0)'*(mu_i_0(cluster,:) - self.rho_0_0));
                    covPart = covPart ./ beta_0(cluster);

                    phiCov = self.Phi_0_0 + nBar0(cluster)*reshape(cov_i_0(cluster,:),d,d) + covPart;
                    Phi_0(cluster,:) = phiCov(:)';
                end

                for cluster = 1:self.K1
                    rho_1(cluster,:) = self.beta_1_m*self.rho_0_1 + nBar1(cluster)*mu_i_1(cluster,:);
                    rho_1(cluster,:) = rho_1(cluster,:)./beta_1(cluster);

                    %See equation 84; Note: this is different from Equation 84 - 84 is
                    %wrong
                    %covPart = nBar1(cluster) * reshape(cov_i_1(cluster,:),d,d) * ((mu_i_1(cluster,:) - rho_1(cluster,:))'*(mu_i_1(cluster,:) - rho_1(cluster,:)));
                    %         covPart = nBar1(cluster) * reshape(cov_i_1(cluster,:),d,d) * ((mu_i_1(cluster,:) - self.rho_0_0)'*(mu_i_1(cluster,:) - self.rho_0_0));
                    covPart = nBar1(cluster) * self.beta_1_m * ((mu_i_1(cluster,:) - self.rho_0_1)'*(mu_i_1(cluster,:) - self.rho_0_1));
                    covPart = covPart ./ beta_1(cluster);

                    phiCov = self.Phi_0_1 + nBar1(cluster)*reshape(cov_i_1(cluster,:),d,d) + covPart;
                    Phi_1(cluster,:) = phiCov(:)';
                end

                % <log_pi>  Equation 102
                matGamma = psi(gamma_i_1_m0) - psi(gamma_i_1_m0 + gamma_i_2_m0);
                matGamma2 = psi(gamma_i_2_m0) - psi(gamma_i_1_m0 + gamma_i_2_m0);
                log_pi_0 = matGamma + cat(2,0,cumsum(matGamma2(1:end-1)));

                matGamma = psi(gamma_i_1_m1) - psi(gamma_i_1_m1 + gamma_i_2_m1);
                matGamma2 = psi(gamma_i_2_m1) - psi(gamma_i_1_m1 + gamma_i_2_m1);
                log_pi_1 = matGamma + cat(2,0,cumsum(matGamma2(1:end-1)));

                % <Log \eta_m>  Equation 105
                nBarM0 = self.alpha(1) +  sum(phi_0_M(expandedBagLabels == 1,:));
                nBarM1 = self.alpha(2) +  sum(phi_1_M(expandedBagLabels == 1,:));
                log_eta = psi([nBarM0,nBarM1]) - psi(nBarM1+nBarM0);

                % <log det ||>
                for cluster = 1:self.K0
                    temp_det0 = sum(psi((nu_0(cluster) + 1 - (1:d))/2));
                    temp_det0 = temp_det0 + d*log(2);
                    log_det_0(cluster) = temp_det0 + log(det(reshape(Phi_0(cluster,:),d,d)^-1));
                end
                for cluster = 1:self.K1
                    temp_det1 = sum(psi((nu_1(cluster) + 1 - (1:d))/2));
                    temp_det1 = temp_det1 + d*log(2);
                    log_det_1(cluster) = temp_det1 + log(det(reshape(Phi_1(cluster,:),d,d)^-1));
                end

                %Equatioin 113
                for cluster = 1:self.K0
                    c = reshape(Phi_0(cluster,:),d,d);
                    %         v = prtUtilCalcDiagXcInvXT(bsxfun(@minus,x,mu_i_0(cluster,:)),c);
                    v = prtUtilCalcDiagXcInvXT(bsxfun(@minus,x,rho_0(cluster,:)),c);
                    inner_0_n(:,cluster) = d/beta_0(cluster) + nu_0(cluster).*v;
                end

                for cluster = 1:self.K1
                    c = reshape(Phi_1(cluster,:),d,d);
                    %         v = prtUtilCalcDiagXcInvXT(bsxfun(@minus,x,mu_i_1(cluster,:)),c);
                    v = prtUtilCalcDiagXcInvXT(bsxfun(@minus,x,rho_1(cluster,:)),c);
                    inner_1_n(:,cluster) = d/beta_1(cluster) + nu_1(cluster).*v;
                end

                for cluster = 1:self.K0
                    eLogPdf0(:,cluster) = -d*log(2*pi) + log_det_0(cluster) - inner_0_n(:,cluster);
                end
                eLogPdf0 = eLogPdf0*1/2;

                for cluster = 1:self.K1
                    eLogPdf1(:,cluster) = -d*log(2*pi) + log_det_1(cluster) - inner_1_n(:,cluster);
                end
                eLogPdf1 = eLogPdf1*1/2;

                phiHat0 = bsxfun(@plus,eLogPdf0,log_pi_0);
                phiHat1 = bsxfun(@plus,eLogPdf1,log_pi_1);
                phiHat1(expandedBagLabels == 0,:) = -inf;

                phi_1_i = bsxfun(@rdivide,exp(phiHat1),sum(exp(phiHat1),2));
                phi_0_i = bsxfun(@rdivide,exp(phiHat0),sum(exp(phiHat0),2));

                phi_1_i(exp(phiHat1) == 0) = 0;
                phi_0_i(exp(phiHat0) == 0) = 0;

                %
                n1 = (log(sum(exp(phiHat1),2))+log_eta(2));
                n0 = (log(sum(exp(phiHat0),2))+log_eta(1));
                phi_1_M = exp(n1)./(exp(n1)+exp(n0));
                phi_1_M(expandedBagLabels == 0) = 0;

                phi_0_M = 1-phi_1_M;

                for i = 1:size(rho_0,1)
                    c = Phi_0(i,:)./nu_0(i);
                    c = reshape(c,d,d);
                    if vbIter == 1
                        mm0(i) = prtRvMvn('mu',rho_0(i,:),'sigma',c);
                    else
                        mm0(i).mu = rho_0(i,:);
                        mm0(i).sigma = c;
                    end
                end
                for i = 1:size(rho_1,1)
                    c = Phi_1(i,:)./nu_1(i);
                    c = reshape(c,d,d);
                    if vbIter == 1
                        mm1(i) = prtRvMvn('mu',rho_1(i,:),'sigma',c);
                    else
                        mm1(i).mu = rho_1(i,:);
                        mm1(i).sigma = c;
                    end
                end

                pi0 = exp(log_pi_0); pi0 = pi0./sum(pi0);
                rv0 = prtRvGmm('nComponents',size(mu_i_0,1),'mixingProportions',pi0,'components',mm0);

                pi1 = exp(log_pi_1); pi1 = pi1./sum(pi1);
                rv1 = prtRvGmm('nComponents',size(mu_i_1,1),'mixingProportions',pi1,'components',mm1);

                self.rvH1 = rv1;
                self.rvH0 = rv0;

                drawnow;
                if ~mod(vbIter,100);
                    disp(vbIter)
                    subplot(2,4,1);
                    stem(exp(log_pi_0));
                    subplot(2,4,2);
                    stem(exp(log_eta))
                    subplot(2,4,3);
                    stem(exp(log_pi_1));
                    fprintf('exp(log_pi_0): %.2f\n',exp(log_pi_0))
                    fprintf('exp(log_pi_1): %.2f\n',exp(log_pi_1))

                    subplot(2,4,4);
                    self.isTrained = true;
                    yOut = self.run(dsMil);
                    prtScoreRoc(yOut);

                    subplot(2,4,5:8);
                    if d > 1
                        plot(x(expandedBagLabels == 0,1),x(expandedBagLabels == 0,2),'b.');
                        hold on;
                        plot(x(expandedBagLabels == 1,1),x(expandedBagLabels == 1,2),'r.');

                        for i = 1:size(rho_0,1)
                            c = Phi_0(i,:)./nu_0(i);
                            c = reshape(c,d,d);
                            if d > 1
                                %                                 plotMvnEllipse(rho_0(i,1:2),c(1:2,1:2),1);
                            end
                            mm0(i) = prtRvMvn('mu',rho_0(i,:),'sigma',c);
                        end
                        for i = 1:size(rho_1,1)
                            c = Phi_1(i,:)./nu_1(i);
                            c = reshape(c,d,d);
                            if d > 1
                                %                                 plotMvnEllipse(rho_1(i,1:2),c(1:2,1:2),1);
                            end
                            mm1(i) = prtRvMvn('mu',rho_1(i,:),'sigma',c);
                        end
                        hold on; h = plot(rho_0(:,1),rho_0(:,2),'ko'); set(h,'MarkerFaceColor','k');
                        hold on; h = plot(rho_1(:,1),rho_1(:,2),'go'); set(h,'MarkerFaceColor','g');

                        hold off;
                    end
                    drawnow;
                end

            end
        end

        function yOut = runAction(self,dsMil)

            for n = 1:dsMil.nObservations
                milStruct = dsMil.data(n);
                data = milStruct.data;

                h1 = self.rvH1.logPdf(data);
                h0 = self.rvH0.logPdf(data);
                l0 = h0'-prtUtilSumExp([h0';h1']);

                h0LogProbabilityBag = sum(l0);
                %exp(h0LogProbabilityBag)

                y(n,1) = 1-exp(h0LogProbabilityBag);
            end

            yOut = prtDataSetClass(y,dsMil.targets);
        end

        function [h0Mix,h1Mix,xH0,xH1] = initialize(self,xBag,bagLabels)

            %H0 Kmeans
            xH0 = cat(1,xBag{bagLabels == 0});
            idx0 = kmeans(xH0,self.K0,'Replicates',1,'EmptyAction','singleton','MaxIter',100);
            uIds = unique(idx0);

            for idx = 1:length(uIds)
                rvH0struct(idx).mean = mean(xH0(idx0 == idx,:),1);

                try
                    rvH0struct(idx).cov = cov(xH0(idx0 == idx,:));
                    chol(rvH0struct(idx).cov);
                    assert(size(xH0,2) == size(rvH0struct(idx).cov,2));
                catch ME
                    rvH0struct(idx).cov = diag(var(xH0(idx0 == idx,:))) + eye(size(xH0,2));
                end

                rvH0(idx) = prtRvMvn('mu',rvH0struct(idx).mean,'sigma',rvH0struct(idx).cov);
                pi(idx) = sum(idx0 == idx)./length(idx0);

            end
            pi = pi./sum(pi);
            h0Mix = prtRvGmm('nComponents',length(rvH0),'mixingProportions',pi,'components',rvH0);

            %Least likely H1
            h1Bags = xBag(bagLabels == 1);
            for h1Ind = 1:length(h1Bags)
                ll = h0Mix.logPdf(h1Bags{h1Ind});
                [~,ind] = min(ll);
                xH1(h1Ind,:) = h1Bags{h1Ind}(ind,:);
            end

            %H1 Kmeans
            idx1 = kmeans(xH1,self.K1,'Replicates',1,'EmptyAction','singleton','MaxIter',100);
            uIds = unique(idx1);

            pi = [];
            for idx = 1:length(uIds)
                rvH1struct(idx).mean = mean(xH1(idx1 == idx,:),1);

                try
                    rvH1struct(idx).cov = cov(xH1(idx0 == idx,:));
                    chol(rvH1struct(idx).cov);
                catch ME
                    rvH1struct(idx).cov = diag(var(xH1)) + eye(size(xH1,2));
                end

                rvH1(idx) = prtRvMvn('mu',rvH1struct(idx).mean,'sigma',rvH1struct(idx).cov);
                pi(idx) = sum(idx1 == idx)./length(idx1);

            end
            pi = pi./sum(pi);
            h1Mix = prtRvGmm('nComponents',length(rvH1),'mixingProportions',pi,'components',rvH1);
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function minimise( varargin )

%MINIMISE Specifiies a convex (or affine) objective to be maximized.

cvx_problem = evalin( 'caller', 'cvx_problem', '[]' );
if isempty( cvx_problem ) || ~isa( cvx_problem, 'cvxprob' ),
    error( 'A cvx problem does not exist in this scope.' );
end
if ~isempty( cvx_problem.direction ),
    if isequal( cvx_problem.direction, 'find' ),
        error( 'Objective functions cannot be added to sets.' );
    else
        error( 'An objective function has already been supplied.' );
    end
end

if nargin < 1,
    error( 'Objective expression missing.' );
elseif iscellstr( varargin ),
    arg = evalin( 'caller', sprintf( '%s ', varargin{:} ) );
elseif nargin > 1,
    error( 'Too many input arguments.' );
else
    arg = varargin{1};
end

if ~isa( arg, 'cvx' ) && ~isa( arg, 'double' ) && ~isa( arg, 'sparse' ),
    error( 'Cannot accept an objective of type ''%s''.', class( arg ) );
end
persistent remap
if isempty( remap ),
    remap = cvx_remap( 'convex', 'log-convex' );
end
vx = remap( cvx_classify( arg ) );
if ~all( vx ),
    error( 'Disciplined convex programming error:\n   Cannot minimise a(n) %s expression.', cvx_class(arg(vx==0),false,true) );
end

newobj( cvx_problem, 'minimize', arg );

% 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.


text:
function t_opf_dc_bpmpd(quiet)
%T_OPF_DC_BPMPD  Tests for DC optimal power flow using BPMPD_MEX solver.

%   MATPOWER
%   Copyright (c) 2004-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

num_tests = 43;

t_begin(num_tests, quiet);

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

casefile = 't_case9_opf';
if quiet
    verbose = 0;
else
    verbose = 0;
end

t0 = 'DC OPF (BPMPD): ';
mpopt = mpoption('out.all', 0, 'verbose', verbose);
mpopt = mpoption(mpopt, 'opf.dc.solver', 'BPMPD');

%% run DC OPF
if have_feature('bpmpd')
    %% set up indices
    ib_data     = [1:BUS_AREA BASE_KV:VMIN];
    ib_voltage  = [VM VA];
    ib_lam      = [LAM_P LAM_Q];
    ib_mu       = [MU_VMAX MU_VMIN];
    ig_data     = [GEN_BUS QMAX QMIN MBASE:APF];
    ig_disp     = [PG QG VG];
    ig_mu       = (MU_PMAX:MU_QMIN);
    ibr_data    = (1:ANGMAX);
    ibr_flow    = (PF:QT);
    ibr_mu      = [MU_SF MU_ST];
    ibr_angmu   = [MU_ANGMIN MU_ANGMAX];

    %% get solved DC power flow case from MAT-file
    load soln9_dcopf;       %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF
    t = t0;
    [baseMVA, bus, gen, gencost, branch, f, success, et] = rundcopf(casefile, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with angle difference limits  -----
    t = [t0 'w/angle diff lims : '];
    mpc = loadcase(casefile);
    mpc.branch(4, ANGMAX) = 3;
    mpc.branch(7, ANGMIN) = -4.5;
    r = rundcopf(mpc, mpopt);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(   f, 6456.7213, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,PG        ),    [99.98497;89.35133;125.66371], 4, [t 'gen dispatch']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    e = zeros(size(branch, 1), 1);
    e(4) = 297.83776;
    e(7) = -26.94788;
    t_is(branch(:,MU_ANGMAX )-branch(:,MU_ANGMIN ), e, 4, [t 'branch ang diff mu']);

    t = [t0 'w/ignored angle diff lims : '];
    mpopt1 = mpoption(mpopt, 'opf.ignore_angle_lim', 1);
    r = rundcopf(mpc, mpopt1);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  run OPF with extra linear user constraints & costs  -----
    %% two new z variables
    %%      0 <= z1, P3 - P1 <= z1
    %%      0 <= z2, P3 - P2 <= z2
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[10;12;13;12;11;14],[-1;1;-1;1;-1;-1],2,14);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [13;14], [1;1], 2, 14);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 1 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 5, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 5, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 4, [t 'user costs']);

    %% with A and N sized for AC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[19;21;25;21;20;26],[-1;1;-1;1;-1;-1],2,26);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [25;26], [1;1], 2, 26);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 2 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 5, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 5, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 4, [t 'user costs']);

    t = [t0 'infeasible : '];
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1], [10;11], [1;1], 1, 14);   %% Pg1 + Pg2
    mpc.u = Inf;
    mpc.l = 600;
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(~success, [t 'no success']);

    %% OPF with all buses isolated
    t = [t0 'all buses isolated : '];
    mpc = loadcase(casefile);
    mpc.bus(:, BUS_TYPE) = NONE;
    try
        r = rundcopf(mpc, mpopt);
        t_is(r.success, 0, 12, [t 'success = 0']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
    end
else
    t_skip(num_tests, 'BPMPD_MEX not available');
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function varargout = deal(X, simplifyFlag)
%DEAL   Deal for chebmatrices.
%   [A,B,C,...] = DEAL(X) assigns to the variables A,B,C... the entries in the
%   chebmatrix X by row-major ordering so that A=X(1,1), B=X(1,2), C=X(1,3),
%   etc. If the number of outputs is less than the number of entries in X,
%   only those entries are assigned. By default, outputs are simplified (i.e.
%   they are passed through chebfun.simplify) before they are returned.
%
%   [A,B,C,...] = DEAL(X, SIMPLIFYFLAG), with SIMPLIFYFLAG = FALSE skips the
%   default simplification step.
%
%   DEAL is especially useful for the solution of differential equations,
%   wherein one may write, for example,
%
%       N = chebop(@(x,u,v) [diff(u) + v; diff(v) - u], [0,1]);
%       N.bc = @(x,u,v) [u(0) - 1; v(0)];
%       [u,v] = deal(N \ 0);  % solution components of a system
%                             % of differential equations
%
%   However, users do not normally need this syntax as it is
%   simpler to write  [u,v] = N \ 0 ; .

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

if ( nargin < 2 )
    % By default, call simplify:
    simplifyFlag = true;
end

if ( nargout > prod(size(X)) )
    error('CHEBFUN:CHEBMATRIX:deal:tooManyOutputs', ...
        ['Number of outputs should be less than or equal ' ...
         'to the number of elements in input chebmatrix.'])
end

X = X.';
blocks = X.blocks(:);

if ( simplifyFlag )
    for blockCounter = 1:length(blocks)
        if ( isa(blocks{blockCounter}, 'chebfun') )
            blocks{blockCounter} = simplify(blocks{blockCounter});
        end
    end
end

varargout = blocks(1:nargout);

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function h = hline(y, varargin)
%HLINE Easy plotting of a horizonal line.
% Usage:
%   hline
%   hline(y)
%   hline(ax, _)
%   hline(_, ...) % plot args (e.g., linespec)
%   h = hline(_)
%
% Args:
%   y: y-value
%
% See also: vline

if nargin < 1; y = []; end
ax = [];
if isax(y)
    ax = y;
    y = [];
    if nargin > 1
        arg1 = varargin{1};
        if isnumeric(arg1)
            y = arg1;
            varargin(1) = [];
        end
    end
end

if isempty(ax)
    ax = gca;
end

if isempty(y)
    y = mean(ylim(ax));
end

state = ax.NextPlot;
ax.NextPlot = 'add';

X = xlim(ax) .* ones(numel(y),1);
Y = [1 1] .* y(:);
if isvector(y)
    X = [X NaN(size(X,1),1)]';
    Y = [Y NaN(size(Y,1),1)]';
end

h = plot(ax, X(:), Y(:), varargin{:});

ax.NextPlot = state;

if nargout < 1; clear h; end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap