text:
% POLYGONS Manipulation of planar polygons and polylines
% Version 1.6 21-Mar-2011 .
%
%   The 'polygons' module contains functions operating on shapes composed
%   of a vertex list, like polygons or polylines.
%
%   We call 'polyline' the curve defined by a series of vertices.
%   A polyline can be either closed or open, depending on whether the last
%   vertex is connected to the first one or not. This can be given as an
%   option is some functions in the module.
%   A 'polygon' is the planar domain delimited by a closed polyline. We
%   sometimes want to consider 'complex polygons', whose boundary is
%   composed of several disjoint domains. The domain defined by a single
%   closed polyline is called 'simple polygon'.
%   We call 'curve' a polyline with many vertices, such that the polyline
%   can be considered as a discrete approximation of a "real" curve.
%
%   A simple polygon or polyline is represented by a N-by-2 array, each row
%   of the array representing the coordinates of a vertex.
%   Simple polygons are assumed to be closed, so there is no need to repeat
%   the first vertex at the end.
%   As both polygons and polylines can be represented by a list of vertex
%   coordinates, some functions also consider the vertex list itself. Such
%   functions are prefixed by 'pointSet'. Also, many functions prefixed by
%   'polygon' or 'polyline' works also on the other type of shape.
%
%   For multiple-connected polygons, the different connected boundaries are
%   separated by a row [NaN NaN].
%
%   For some functions, the orientation of the polygon can be relevant: CCW
%   stands for 'Conter-Clockwise' (positive orientation), CW stands for
%   'Clockwise'.
%
%   Polylines are parametrized in the following way:
%   * the i-th vertex is located at position i-1
%   * points of the i-th edge have positions ranging linearly from i-1 to i
%   The parametrization domain for an open polyline is from 0 to Nv-1, and
%   from 0 to Nv for a closed polyline (positions 0 and Nv correspond to
%   the same point).
%
%   Example:
%   % Simple polygon:
%   P1 = [1 1;2 1;2 2;1 2];
%   drawPolygon(P1);
%   axis([0 5 0 5]);
%   % Multiple polygon:
%   P2 = [10 10;40 10;40 40;10 40;NaN NaN;20 20;20 30;30 30;30 20];
%   figure;drawPolygon(P2); axis([0 50 0 50]);
%
%
% Polylines
%   polylinePoint             - Extract a point from a polyline
%   polylineLength            - Return length of a polyline given as a list of points
%   polylineCentroid          - Compute centroid of a curve defined by a series of points
%   polylineSubcurve          - Extract a portion of a polyline
%   resamplePolyline          - Distribute N points equally spaced on a polyline
%   resamplePolylineByLength  - Resample a polyline with a fixed sampling step
%   reversePolyline           - Reverse a polyline, by iterating vertices from the end
%   isPointOnPolyline         - Test if a point belongs to a polyline
%   projPointOnPolyline       - Compute position of a point projected on a polyline
%   distancePointPolyline     - Compute shortest distance between a point and a polyline
%   distancePolylines         - Compute the shortest distance between 2 polylines
%   intersectLinePolyline     - Intersection points between a line and a polyline
%   intersectPolylines        - Find the common points between 2 polylines
%   polylineSelfIntersections - Find self-intersection points of a polyline
%   simplifyPolyline          - Douglas-Peucker simplification of a polyline
%   smoothPolyline            - Smooth a polyline using local averaging
%   removeMultipleVertices    - Remove multiple vertices of a polygon or polyline
%
% Polygon basic manipulation
%   polygonPoint              - Extract a point from a polygon
%   polygonSubcurve           - Extract a portion of a polygon
%   polygonEdges              - Return the edges of a simple or multiple polygon
%   reversePolygon            - Reverse a polygon, by iterating vertices from the end
%   smoothPolygon             - Smooth a polygon using local averaging
%   simplifyPolygon           - Douglas-Peucker simplification of a polygon
%   projPointOnPolygon        - Compute position of a point projected on a polygon
%   splitPolygons             - Convert a NaN separated polygon list to a cell array of polygons
%   polygonLoops              - Divide a possibly self-intersecting polygon into a set of simple loops
%
% Polygon clipping and intersections
%   intersectLinePolygon      - Intersection points between a line and a polygon
%   intersectRayPolygon       - Intersection points between a ray and a polygon
%   intersectEdgePolygon      - Intersection point of an edge with a polygon
%   polygonSelfIntersections  - Find self-intersection points of a polygon
%   clipPolygon               - Clip a polygon with a rectangular box
%   clipPolygonHP             - Clip a polygon with a Half-plane defined by a directed line
%
% Point Sets
%   pointSetsAverage          - Compute the average of several point sets
%   minimumCaliperDiameter    - Minimum caliper diameter of a set of points
%   findPoint                 - Find index of a point in an set from its coordinates
%   convexHull                - Convex hull of a set of points
%   randomPointInPolygon      - Generate random point(s) in a polygon
%
% Measures on Polygons
%   isPointInPolygon          - Test if a point is located inside a polygon
%   polygonContains           - Test if a point is contained in a multiply connected polygon
%   polygonCentroid           - Compute the centroid (center of mass) of a polygon
%   polygonArea               - Compute the signed area of a polygon
%   polygonInertiaEllipse     - Compute ellipse with same inertia moments as polygon
%   polygonSecondAreaMoments  - Compute second-order area moments of a polygon
%   polygonLength             - Perimeter of a polygon
%   polygonNormalAngle        - Compute the normal angle at a vertex of the polygon
%   polygonBounds             - Compute the bounding box of a polygon
%   polygonOuterNormal        - Outer normal vector for a given vertex(ices)
%   distancePointPolygon      - Shortest distance between a point and a polygon
%   distancePolygons          - Compute the shortest distance between 2 polygons
%   distancePolygonsNoCross   - Compute the shortest distance between 2 polygons
%   polygonSignature          - Polar signature of a polygon (polar distance to origin)
%   signatureToPolygon        - Reconstruct a polygon from its polar signature
%
% More complex operations on polygons
%   resamplePolygon           - Distribute N points equally spaced on a polygon
%   resamplePolygonByLength   - Resample a polygon with a fixed sampling step
%   densifyPolygon            - Add several points on each edge of the polygon
%   expandPolygon             - Expand a polygon by a given (signed) distance
%   triangulatePolygon        - Compute a triangulation of the polygon
%   polygonSymmetryAxis       - Try to identify symmetry axis of polygon
%   medialAxisConvex          - Compute medial axis of a convex polygon
%
% Curves (polylines with lot of vertices)
%   parametrize               - Parametrization of a polyline, based on edges lengths
%   curvature                 - Estimate curvature of a polyline defined by points
%   cart2geod                 - Convert cartesian coordinates to geodesic coord.
%   geod2cart                 - Convert geodesic coordinates to cartesian coord.
%   curveMoment               - Compute inertia moment of a 2D curve
%   curveCMoment              - Compute centered inertia moment of a 2D curve
%   curveCSMoment             - Compute centered scaled moment of a 2D curve
%
% Functions from stochastic geometry
%   steinerPoint              - Compute steiner point (weighted centroid) of a polygon
%   steinerPolygon            - Create a Steiner polygon from a set of vectors
%   supportFunction           - Compute support function of a polygon
%   convexification           - Compute the convexification of a polygon
%
% Input, Output and conversions
%   polygonToRow              - Convert polygon coordinates to a row vector
%   rowToPolygon              - Create a polygon from a row vector
%   contourMatrixToPolylines  - Converts a contour matrix array into a polyline set
%   readPolygonSet            - Read a set of simple polygons stored in a file
%   writePolygonSet           - Write a set of simple polygons into a file
%
% Drawing functions
%   drawPolyline              - Draw a polyline specified by a list of points
%   drawPolygon               - Draw a polygon specified by a list of points
%   fillPolygon               - Fill a polygon specified by a list of points
%   drawVertices              - Draw the vertices of a polygon or polyline
%
%
%   Credits:
%   * function intersectPolylines uses the 'interX' contribution from "NS"
%       (file exchange 22441, called 'curve-intersections')
%
% -----
% Author: David Legland
% e-mail: david.legland@inra.fr
% created the  07/11/2005.
% Project homepage: http://github.com/mattools/matGeom
% http://www.pfl-cepia.inra.fr/index.php?page=geom2d
% Copyright INRA - Cepia Software Platform.

help('Contents');

%% Requires further development

%% Others...

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [dose_ratio] = dicomrt_doseratio(doseone,dosetwo)
% dicomrt_doseratio(doseone,dosetwo)
%
% Calculate dose ratio between two 3D matrices
%
% doseone and dosetwo can be rtplan and/or monte carlo 3D dose distributions.
%
% NOTE:
% Warnings for divisions by zero are switched off during the call to this function.
% Infinite values are set to nan (not-a-number).
%
% Example:
%
% [doseratio]=dicomrt_doseratio(doseone,dosetwo)
%
% returns in doseratio the ratio: dosetwo/doseone.
%
% See also: dicomrt_dosediff
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check case and set-up some parameters and variables
[doseone_dose_temp,type_doseone_dose,labeld1,PatientPosition]=dicomrt_checkinput(doseone);
[dosetwo_dose_temp,type_dosetwo_dose,labeld1,PatientPosition]=dicomrt_checkinput(dosetwo);

doseone_dose=dicomrt_varfilter(doseone_dose_temp);
dosetwo_dose=dicomrt_varfilter(dosetwo_dose_temp);

% Perform ratio: swith off and back on warnings
warning off MATLAB:divideByZero;
dose_ratio=dosetwo_dose./doseone_dose;
warning on MATLAB:divideByZero;

% Discard infinite values
dose_ratio(find(isinf(dose_ratio)==1))=nan;

% Restore original variable format
[dose_ratio]=dicomrt_restorevarformat(dosetwo,dose_ratio);

% Label Plan and update time of creation
if iscell(dose_ratio)==1
    dose_ratio{1,1}{1}.RTPlanLabel=[dose_ratio{1,1}{1}.RTPlanLabel,'-DRATIO'];
    dose_ratio{1,1}{1}.RTPlanDate=date;
    time=fix(clock);
    creationtime=[num2str(time(4)),':',num2str(time(5))];
    dose_ratio{1,1}{1}.RTPlanTime=creationtime;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function plotThermoModelExtractStats(model, activeInactiveRxn, rxnWeights, presentAbsentMet, metWeights, thermoModelMetBool, thermoModelRxnBool)

nMet=length(thermoModelMetBool);
nRxn=length(thermoModelRxnBool);

if 0
    figure
    H1 = subplot(2,1,1);
    plot(abs(solution.v),g0orig,'.')
    xlim(H1,[-1/max(abs(solution.v)),max(abs(solution.v))])
    xlabel('Reaction flux magnitude')
    ylabel('Reaction weight')
    H2 = subplot(2,1,2);
    plot(abs(solution.d),model.h0,'.')
    xlim(H2,[-1/max(abs(solution.d)),max(abs(solution.d))])
    xlabel('Metabolite production rate')
    ylabel('Metabolite weight')
else

    if 0
        T = false(3,nRxn);
        T(1,:) = abs(solution.v)>=param.epsilon;
        T(2,:) = abs(solution.v)>=param.epsilon | abs(solution.v)<param.epsilon;
        T(2,:) = abs(solution.v)<param.epsilon;
        P = false(3,nRxn);
        P(1,:) = g0orig < 0;
        P(2,:) = g0orig == 0;
        P(3,:) = g0orig > 0;
        plotconfusion(T,P,'Reaction activity')
    else

        if 0
            T = cell(nRxn,1);
            T(abs(solution.v)>=param.epsilon & g0orig ~= 0)={'Active'};
            T(abs(solution.v)<param.epsilon & g0orig ~= 0)={'Inactive'};
            T(g0orig == 0)={'Unspecified'};
            P = cell(nRxn,1);
            P(g0orig < 0)={'Active'};
            P(g0orig == 0)={'Unspecified'};
            P(g0orig > 0)={'Inactive'};
            %             usage:
            %             % Find the confusionmat matrix first
            %             [C,order] = confusionmat(TestResults,PredictionResults)
            %             % Then, plot
            %             plotConfMat(C) plots the confmat with integers 1 to n as class labels
            %             plotConfMat(C, order) plots the confmat with the specified labels
            C = confusionmat(T,P,'order',{'Active';'Inactive';'Unspecified'});
            plotConfMat(C,{'Active';'Inactive';'Unspecified'});
        else

            if any(rxnWeights~=0)
                figure;
                hist(rxnWeights)
                title('Reaction weights')
                figure
                P = cell(nRxn,1);
                P(thermoModelRxnBool==1)={'Active'};
                P(thermoModelRxnBool==0)={'Inactive'};

                T = cell(nRxn,1);
                if 1
                    T(rxnWeights < 0)={'Active'};
                    T(rxnWeights > 0)={'Inactive'};
                    T=T(rxnWeights~=0);
                    P=P(rxnWeights~=0);
                else
                    T(activeInactiveRxn==1)={'Active'};
                    T(activeInactiveRxn==-1)={'Inactive'};
                    T=T(activeInactiveRxn~=0);
                    P=P(activeInactiveRxn~=0);
                end

                %             usage:
                %             % Find the confusionmat matrix first
                %             [C,order] = confusionmat(TestResults,PredictionResults)
                %             % Then, plot
                %             plotConfMat(C) plots the confmat with integers 1 to n as class labels
                %             plotConfMat(C, order) plots the confmat with the specified labels

                %    CM = CONFUSIONMAT(G,GHAT) returns the confusion matrix CM determined

                %    by the known group labels G and the predicted group labels GHAT.
                if any(rxnWeights == 0) && 0
                    labels = {'Active';'Inactive';'Unspecified'};
                else
                    labels = {'Active';'Inactive'};
                end
                C = confusionmat(P,T,'order',labels);
                plotConfMat(C,labels);
                accuracy = sum(diag(C),1)/sum(sum(C,1));
                title(['Reaction confusion matrix, accuracy = ' num2str(accuracy)])
            end

            if any(metWeights~=0)
                figure;
                hist(metWeights)
                title('Metabolite weights')

                figure;
                P = cell(nMet,1);
                P(thermoModelMetBool==1)={'Active'};
                P(thermoModelMetBool==0)={'Inactive'};

                T = cell(nMet,1);
                if 1
                    T(metWeights < 0)={'Active'};
                    T(metWeights >= 0)={'Inactive'};
                    T=T(metWeights~=0);
                    P=P(metWeights~=0);
                else
                    T(presentAbsentMet==1)={'Active'};
                    T(presentAbsentMet==-1)={'Inactive'};
                    T=T(presentAbsentMet~=0);
                    P=P(presentAbsentMet~=0);
                end
                %             usage:
                %             % Find the confusionmat matrix first
                %             [C,order] = confusionmat(TestResults,PredictionResults)
                %             % Then, plot
                %             plotConfMat(C) plots the confmat with integers 1 to n as class labels
                %             plotConfMat(C, order) plots the confmat with the specified labels

                %    CM = CONFUSIONMAT(G,GHAT) returns the confusion matrix CM determined

                %    by the known group labels G and the predicted group labels GHAT.
                if any(metWeights == 0) && 0
                    labels = {'Active';'Inactive';'Unspecified'};
                else
                    labels = {'Active';'Inactive'};
                end
                C = confusionmat(P,T,'order',labels);
                plotConfMat(C,labels);
                accuracy = sum(diag(C),1)/sum(sum(C,1));
                title(['Metabolite confusion matrix, accuracy = ' num2str(accuracy)])
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function dspace = dspacing(h)
% space between crystal planes

dspace = 1./norm(h);

author: mtex-toolbox

repo: mtex

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


text:
function [var_prc,pci,Nneeded,pcitarget] = var_prctile(p,n_in_sample,varargin)
% :Usage:
% ::
%
%     [var_prc,pci,Nneeded,pcitarget] = var_prctile(p,n_in_sample,['nboot',nboot],['x',data])
%
%     [var_prc,pci,Nneeded,pcitarget] = var_prctile(p,50,'nboot',5000)
%
% :Inputs:
%
%   **p:**
%        is desired prctile of data (threshold)
%
%   **nboot:**
%        is number of bootstrap samples (if bootstrapping)
%
%   **n_in_sample:**
%        is number of obs. in original sample
%
%   **x:**
%        is data sample of distribution of interest (empirical pdf based on
%        this)
%
%        Note: using empirical PDF/CDF depends a great deal on choice of h (see
%        code).
%
% :Examples:
% ::
%
%    Nneeded = [];
%    for p = [.05:-.001:.001]
%        [var_prc,pci,Nneeded(end+1),pcitarget] = var_prctile(x,p);
%    end
%    figure;
%    plot([.05:-.001:.001],Nneeded)
%
% ..
%    tor wager, jan 2007
% ..

nboot = Inf;
norm_model = 1;
x = [];
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case 'nboot', nboot = varargin{i+1};                    % # bootstrap samples
            case {'x','data'}, x = varargin{i+1}; norm_model = 0;   % use empirical PDF

            otherwise
                error('Unknown string option');
        end
    end
end

% ------------------------------------------------------------
% * get PDF
% ------------------------------------------------------------

if norm_model
    % get normal PDF at p-th percentile
    % -------------------------------------------
    prc = norminv(p);
    pdfval = normpdf(prc);
else
    % empirical estimate of PDF from x (data)
    % -------------------------------------------
    % prc is x-score at prctile of interest
    prc = prctile(x,100*p);

    % now we need pdf at that prctile.
    % pick some unit h, and differentiate around prc

    h = max(.01,1000 ./ length(x));  % enough so we have reasonable idea
    h = min(h,p);
    h = max(h,1000 ./ length(x));

    pdfval = 0;

    while pdfval == 0
        pdfval = emppdf(x,prc,h);

        h = h + .01;
    end
end

% ------------------------------------------------------------
% * get variance
% ------------------------------------------------------------

%for normal: fit = ( 1./ (normpdf(norminv(p)).^2) ) .* (p *(1-p)./N);
% from Brown and Wolfe, asymptotic
%var_prc = ( 1./ (pdfval.^2) ) .* (p *(1-p)./N);

% from Martin's bootstrap book, Ch 19, Eq. 19.7, p. 275
var_prc = ( 1./ (pdfval.^2) ) .* (p *(1-p)) .* (1./nboot + 1./n_in_sample);

if nargout < 2, return, end

% ------------------------------------------------------------
% * get confidence interval
% ------------------------------------------------------------

halfci_prc = 1.96 .* sqrt(var_prc);

% lower and upper bounds on p-value derived from distribution of x
pci = get_pval_ci(x,prc,halfci_prc);

if nargout < 3, return, end

% now get nboot needed to achieve target
% ---------------------------------------

% expression for tolerance: if upper bound of p-value is within 10% of p-value
maxN = 50000;               % maximum number of iterations
Nneeded = 500;              % starting estimate for nboot needed (should be fairly low)
Nstep = 500;                % increase Nneeded in units of Nstep until satisfied
p_upper_bound = p + .1 * p; % upper bound for p desired; now 10% larger than p
pcitarget = Inf;            % target conf. interval for p-values with Nneeded boot samples

% computations we don't have to repeat.  Divide by sqrt(N) to get halfci of prc
halfci_squared_determiner = 1.96 .* sqrt( ( 1./ (pdfval.^2) ) .* p *(1-p) );

while ( max(pcitarget) > p_upper_bound ) && ( Nneeded < maxN )

    Nneeded = Nneeded + Nstep;
    halfci_prc = halfci_squared_determiner ./ sqrt(Nneeded);
    pcitarget = get_pval_ci(x,prc,halfci_prc);

end

end

function pci = get_pval_ci(x,prc,halfci_prc)
% get confidence interval for p-value, based on data x, x-axis-value prc,
% and confidence half-interval for prc

if isempty(x)
    % normal CDF
    pci = [max(0,normcdf(prc - halfci_prc)) min(1,normcdf(prc + halfci_prc))];
else
    % empirical CDF
    pci = [max(0,empcdf(x,prc - halfci_prc)) min(1,empcdf(x,prc + halfci_prc))];
end

end

function pdfval = emppdf(x,prc,h)
% empirical PDF of data numerically differentiated in a window h
pdfval = ( sum(x <= prc + h) - sum(x <= prc - h) ) ./ (length(x) * 2 * h);

end

function p = empcdf(x,prc)
% empirical CDF of data x at prctile prc

p = sum(x <= prc) ./ length(x);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function h = drawCoorSys(origin, R, length, width)
    % R - orintation of frame
    x_axis = R(:,1);
    y_axis = R(:,2);
    z_axis = R(:,3);

    % draw three originamera axis
    h = hggroup;
    line([origin(1) origin(1) + length*x_axis(1)], [origin(2) origin(2) + length*x_axis(2)]...
        , [origin(3) origin(3) + length*x_axis(3)], 'Color', 'r', 'LineWidth', width, 'Parent', h);
    line([origin(1) origin(1) + length*y_axis(1)], [origin(2) origin(2) + length*y_axis(2)]...
        , [origin(3) origin(3) + length*y_axis(3)], 'Color', 'g', 'LineWidth', width, 'Parent', h);
    line([origin(1) origin(1) + length*z_axis(1)], [origin(2) origin(2) + length*z_axis(2)]...
        , [origin(3) origin(3) + length*z_axis(3)], 'Color', 'b', 'LineWidth', width, 'Parent', h);
end

author: zhixuany

repo: HUMBI

save_path: github-repos/MATLAB/zhixuany-HUMBI


text:
function [rs, eigvals, eigvecs] = pca(s, mode, maxpercent)

%tstoolbox/@signal/pca
%   Syntax:
%     * [rs, eigvals, eigvecs] = pca(s) => mode='normalized' , maxpercent
%       = 95
%     * [rs, eigvals, eigvecs] = pca(s, mode) => maxpercent = 95
%     * [rs, eigvals, eigvecs] = pca(s, mode, maxpercent)
%
%   Input arguments:
%     * each row of data is one 'observation', e.g. the sample values of
%       all channels in a multichannel measurement at one point in time
%     * mode can be one of the following : 'normalized' (default), 'mean',
%       'raw'
%          + in mode 'normalized' each column of data is centered by
%            removing its mean and then normalized by dividing through its
%            standard deviation before the covariance matrix is calculated
%          + in mode 'mean' only the mean of every column of data is
%            removed
%          + in mode 'raw' no preprocessing is applied to data
%     * maxpercent gives the limit of the accumulated percentage of the
%       resulting eigenvalues, default is 95 %
%
%   Principal component analysis of column orientated data set.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

%
%   principal component analysis of column orientated data set <data>
%
%   input arguments :
%
%   - each row of data is one 'observation', e.g. the sample values of
%     all channels in a multichannel measurement at one point in time
%
%   - mode can be one of the following : 'normalized' (default), 'mean', 'raw'
%     - in mode 'normalized' each column of data is centered by removing its mean
%       and then normalized by dividing through its standard deviation before
%       the covariance matrix is calculated
%     - in mode 'mean' only the mean of every column of data is removed
%     - in mode 'raw' no preprocessing is applied to data
%
%   - maxpercent gives the limit of the accumulated percentage of the resulting
%     eigenvalues, default is 95 %
%
% [rs, eigvals, eigvecs] = pca(s)       => mode='normalized' , maxpercent = 95
% [rs, eigvals, eigvecs] = pca(s, mode)                     => maxpercent = 95
% [rs, eigvals, eigvecs] = pca(s, mode, maxpercent)
%
% C.Merkwirth,U.Parlitz,W.Lauterborn  DPI Goettingen 1998

narginchk(1,3);

if nargin < 2
	mode = 'normalized';
end
if nargin < 3
	maxpercent = 95;
end

if ndim(s)~=2
	error('pca needs a signal with two dimensions as input')
end

dat = data(s);
n = dlens(s,1);
m = dlens(s,2);

htext = {''};
htext{end+1} = 'Applied Karhunen-Loeve Transform (pca)';
htext{end+1} = ['Tried to capture ' num2str(maxpercent) ' percent of total variance'];

mpercent = maxpercent/100;
mode = lower(mode); 		% no problems with uppercase letters

if strncmp(mode, 'r',1)
	mode = 'raw';
	htext{end+1} = 'No data preprocessing';
elseif strncmp(mode, 'm',1)
	mode = 'mean';
	htext{end+1} = 'Removed mean from data set';
	mn = mean(dat);
	dat = dat - repmat(mn, n, 1);
else
	mode = 'normalized';
	htext{end+1} = 'Removed mean and normalized data';
	mn = mean(dat);
	dv = std(dat);
	dat = dat - repmat(mn, n, 1);
	dat = dat ./ repmat(dv, n, 1);
end

if n>m
	htext{end+1} = 'Used direct method to compute covariance matrix';
	K =   dat' * dat;       % oder K = corrcoef(dat)
	[Q,D] = eig(K);
 else
    htext{end+1} = 'using indirect method to compute covariance matrix';
 	C = dat * dat';
 	[Q,D] = eig(C);
end

[evalues,index] = sort(diag(D));
evalues = flipud(evalues);
index = flipud(index);

total = sum(evalues);
rlvm = min(find(cumsum(evalues) >= (total*mpercent)));

if isempty(rlvm)			% in case percentage was choosen over 100 %,
	rlvm = length(evalues); % return all eigenvalues
end

frvals = evalues(1:rlvm);

if n>m
	frvecs = Q(:, index(1:rlvm));
	trnsfrmd=dat*frvecs;
else
 	scalefac = 1./ sqrt(evalues(1:rlvm));
 	for i = 1:rlvm
 		P(:,i) = Q(:,index(i)) * scalefac(i);
 	end
 	frvecs = dat' * P; % '
	trnsfrmd=C*P;
end

a = achse(unit, 1,1);
a = setname(a, 'Mode');

rs = signal(core(trnsfrmd), s);
rs = addhistory(rs, htext);
rs = addcommandlines(rs, 's = pca(s', mode, maxpercent);
rs = setaxis(rs, 2, a);
%rs = settype(rs, 'Transformed data set');
rs = setplothint(rs, 'multigraph');

eigvals = signal(core(frvals),s);
eigvals = addhistory(eigvals, htext);
eigvals = addcommandlines(eigvals, '[dummy, s] = pca(s', mode, maxpercent);
eigvals = setaxis(eigvals, 1, a);
%eigvals = settype(eigvals, 'Eigenvalues');
eigvals = setplothint(eigvals, 'bar');

eigvecs = signal(core(frvecs),s);
eigvecs = addhistory(eigvecs, htext);
eigvecs = addcommandlines(eigvecs, '[dummy, dummy2, s] = pca(s', mode, maxpercent);
eigvecs = setaxis(eigvecs, 2, a);
eigvecs = setaxis(eigvecs, 1, getaxis(s,2));
%eigvecs = settype(eigvecs, 'Eigenvectors');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [retRes,retAlgo] = training(retAlgo,dat,lossType)

%  Does cross validation and returns the results as a data object
%  plus the updated group (after training models).

  rndState = rand( 'state');
  if retAlgo.algorithm.verbosity
    disp(['training ' get_name(retAlgo) '.... '])
  end
  if nargin==1
    dat=data; % <-- in case data is not specified
  end;

  if nargin<3
    lossType=[];  % <-- in case loss is not specified
  end

  untrained = retAlgo.child{1};
  if isa(untrained,'cell')
    untrained = group(untrained); % <-- convert retAlgo.child{1} to group
  end;

  [num,vecDim,outDim]=get_dim(dat);
  orig_name=get_name(dat);
  retRes = {};

  retAlgo.trialbytrial=[];
  for loop=1:retAlgo.repeats % <-- do cross validation vecDim times

    rand('seed',loop);
    if retAlgo.balanced     % <-- balance the number of positives/negatives labels in folds
      labels=get_y(dat);
      fin1=find(labels==1);
      fin2=find(labels==-1);

      pe1=randperm(sum(labels==1));
      pe1=fin1(pe1);

      pe2=randperm(sum(labels==-1));
      pe2=fin2(pe2);
      perm=[pe1 ;pe2];
    else
      perm=randperm(num);
    end

    trialbytrial = [];
    for i=1:retAlgo.folds
      if retAlgo.train_on_fold==0
        tst=perm(i:retAlgo.folds:num);
        trn=setdiff(1:num,tst);
      else  % <--- take fraction as training data
        trn=perm(i:retAlgo.folds:num);
        tst=setdiff(1:num,trn);
      end

      if retAlgo.store_all
        indx=i+(loop-1)*(retAlgo.folds);
      else
        indx=1;
      end;

      str=[' fold=' num2str(i+(loop-1)*(retAlgo.folds))];
      if retAlgo.output_train_error==0  % <--- testing on left out fold
        dat=set_name(dat,[ orig_name ' -> cv' str]);
      else
        dat=set_name(dat,[ orig_name ' -> cv' str ' output_train_error=1']);
      end

      if isfield(struct(untrained), 'test_on_the_fly')
        if isstruct(untrained.test_on_the_fly)
          %% untrained.test_on_the_fly.data = set_name(get(dat, tst), [orig_name ' -> cv' str]);
          % line above fails because subsasgn has suddenly stopped working recursively
          % work around it as follows:
          sss = untrained.test_on_the_fly;
          sss.data = set_name(get(dat, tst), [orig_name ' -> cv' str]);
          untrained.test_on_the_fly = sss;
          % some algorithms train another algorithm multiple times in the course of
          % their own training: for example, RFE trains a classifier once per
          % feature elimination step. It is often handy to evaluate cv-error at each
          % stage: previously this involved either re-training the individual classifiers
          % afterwards (time-consuming), or storing a whole series of trained
          % classifiers (memory-consuming). This solution gives the wrapper
          % a sneak preview of the test data for this fold, so that cv error can be
          % evaluated as it goes. Of course, the sneak preview should NOT be used
          % for training!
        end
      end
      if retAlgo.algorithm.verbosity
        str = sprintf('training cv fold %d of %d', i, retAlgo.folds);
        if retAlgo.repeats > 1, str = sprintf('%s (repeat %d of %d)', str, loop, retAlgo.repeats); end
        disp(str)
      end

      if retAlgo.output_train_error
        [res, retAlgo.child{indx}] = train(untrained, get(dat, trn), lossType);
      else
        [res, retAlgo.child{indx}] = traintest(untrained, dat, trn, tst, lossType);
      end
      if retAlgo.store_trialbytrial
        trialbytrial = [trialbytrial; tst(:) repmat(i, length(tst), 1) res.Y res.X];
      end
      %	[res retAlgo.child{indx}]=train(untrained,get(dat,trn));
      %    if retAlgo.output_train_error==0  % <--- test on left out fold
      %      [res]=test(retAlgo.child{indx},get(dat,tst));
      %    end
      %    if ~isempty(lossType)        % <--- calculate loss on results
      %      res=loss(res,lossType,[],1);
      %    end

      if ~isa(res, 'cell'), res = {res}; end
      retRes = [retRes res(:)'];

    end
    if retAlgo.store_trialbytrial
      retAlgo.trialbytrial = cat(3, retAlgo.trialbytrial, sortrows(trialbytrial));
    end
  end

  retRes=group(retRes);

  rand( 'state', rndState);  % restore the random generator


text:
function f = uplus(f)
%UPLUS   Unary plus for a SEPARABLEAPPROX.

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

% The algorithm below is advocated by Vladimir Rohklin:

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function u = acsc(a)
%ACSC         Slope inverse cosecant acsc(a)
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  INTLAB_SLOPE = getappdata(0,'INTLAB_SLOPE');

  u = a;

  u.r = acsc(a.r);
  indexc = 1:INTLAB_SLOPE.NUMVAR;
  indexr = 2:INTLAB_SLOPE.NUMVAR+1;
  Xxs = hull(a.r(:,indexc),a.r(:,indexr));
  Index = 1:size(a.r.inf,1);

  index = all( a.r.sup<=0 , 2);
  if any(index)
    aindex.r = a.r(index,:);
    aindex.s = a.s(index,:);
    u.s(index,:) = ...
      slopeconvexconcave('acsc','-1./(abs(%).*sqrt(sqr(%)-1))',aindex,0);
    Index(index) = 0;
  end

  index = all( a.r.inf>=0 , 2);
  if any(index)
    aindex.r = a.r(index,:);
    aindex.s = a.s(index,:);
    u.s(index,:) = ...
      slopeconvexconcave('acsc','-1./(abs(%).*sqrt(sqr(%)-1))',aindex,1);
    Index(index) = 0;
  end

  if any(Index)
    Index( Index==0 ) = [];
    Xxs = Xxs(Index);
    u.s(Index,:) = - a.s(Index,:) ./ ( abs(Xxs) .* sqrt( sqr(Xxs)-1 ) );
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


function in = iresponse_proxy(in)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'iresponse_proxy' computes the impulse response functions to shock
% identified via instrumental variables
% Codes are written based on the Mertens Ravn (2013, AER) source codes,
% modified to
%1) handle different length of VAR and factors
%2) Multiple instruments to explain the same shock

% Filippo Ferroni, 6/1/2015
% Revised, 2/15/2017
% Revised, 3/21/2018
% Revised, 9/11/2019
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y         = in.vars(in.p+1:end,:);
[T,n]     = size(Y);
[T_m,n_m] = size(in.proxies);

% number of proxies
k  = 1;
%Assuming proxies start at least p periods later
instrument  = in.proxies(1:end,:);

% Identification
%%%%%%%%%%%%%%%%%

% covariance of the reduced form shocks
in.Sigma_m = in.Sigma;

% Instrument on VAR residuals
Phib = [ones(T_m,1) instrument]\in.res(T-T_m-in.T_m_end+1:T-in.T_m_end,:);

% Fitted values of the identified shocks
% here is where thes prior on beta should enter , instead of Phib(:,1)
% m = beta e1 + v
uhat1           =   [ones(T_m,1) instrument]*Phib(:,1);

% regress the fitted values on the other reduced form shocks
b21ib11_TSLS    =   [ones(T_m,1) uhat1]\in.res(T-T_m-in.T_m_end+1:T-in.T_m_end,2:end);
b21ib11_TSLS    =   b21ib11_TSLS(2:end,:)';
b21ib11         =   b21ib11_TSLS;

% Identification of b11 and b12 from the covariance matrix of the VAR
Sig11   = in.Sigma_m(1:k,1:k);
Sig21   = in.Sigma_m(k+1:n,1:k);
Sig22   = in.Sigma_m(k+1:n,k+1:n);
ZZp     = b21ib11*Sig11*b21ib11'-(Sig21*b21ib11'+b21ib11*Sig21')+Sig22;
b12b12p = (Sig21- b21ib11*Sig11)'*(ZZp\(Sig21- b21ib11*Sig11));
b11b11p = Sig11-b12b12p;
b11     = sqrt(b11b11p);
in.b1   = [b11; b21ib11*b11];
in.Phib = Phib;

% Impulse Responses
%%%%%%%%%%%%%%%%%%%%
% initial shock: eps(1,1)=1
irs(in.p+1,:) = in.b1(:,1);

for jj = 2:in.irhor%+max(max(VAR.term_spreads_matur),max(VAR.real_rates_init+VAR.real_rates_matur-1))
    lvars = (irs(in.p+jj-1:-1:jj,:))';
    irs(in.p+jj,:) = lvars(:)'*in.Phi(1:in.p * n,:);
end

in.irs   = irs(in.p+1:in.p+in.irhor,:);
in.uhat1 = uhat1;

if in.compute_F_stat == 1
    % F-test
    %%%%%%%%%%%%%%%%%%%%%%%%%
    XX_m      = [ones(T_m,1) instrument];
    Res_m     = in.res( T-T_m-in.T_m_end+1 : T-in.T_m_end , :)- XX_m * Phib;
    Res_const = in.res( T-T_m-in.T_m_end+1 : T-in.T_m_end,1) - ...
        ones(T_m,1)*(ones(T_m,1)\in.res(T-T_m-in.T_m_end+1 : T-in.T_m_end,1));

    SST_m      = Res_const'*Res_const;
    SSE_m      = Res_m(:,1)'*Res_m(:,1);
    in.F_m     = ((SST_m-SSE_m)/n_m)/ ...
        (SSE_m/(length(in.res(T-T_m-in.T_m_end+1 : T-in.T_m_end,1))-(n_m+1)));
    in.R2_m    = (1-SSE_m/(SST_m));
    in.R2adj_m = in.R2_m-...
        (n_m / ((length(in.res(T-T_m-in.T_m_end+1 : T-in.T_m_end,1)) -(n_m+1))))*(1-in.R2_m);

    % Calculate robust standard errors
    SS_m = zeros(n_m+1,n_m+1);
    for ii = 1 : T_m
        SS_m = SS_m+1/T_m * XX_m(ii,:)'*XX_m(ii,:)*Res_m(ii,1)^2;
    end
    Avarb_m    = inv(1/T_m * XX_m'*XX_m)*SS_m *inv(1/T_m * XX_m'*XX_m);
    RR_m       = [zeros(n_m,1) eye(n_m)];
    WW_m       = T_m*(RR_m*Phib(:,1))' * inv(RR_m*Avarb_m*RR_m') * (RR_m*Phib(:,1));
    in.F_m_rob = WW_m/n_m;
end


text:
function handle = colorbar_axis(s,loc,clabel,rlab1,rlab2, fontsize)
% COLORBAR - Display color bar (color scale).
%
%   This function differs from colorbar(loc,clabel,rlab1,rlab2) in that
%   only a portion of the colorbar defined by the plot axis "paxis" is shown
%
%       s.dBlims is the [2] vector containing the plotting axis limits
%       paxis   is the [2] vector containing the plotting axis limits
%               that are used to limit the portion of the colorbar shown
%
%	COLORBAR('vert') appends a vertical color scale to
%	the current axis. COLORBAR('horiz') appends a
%	horizontal color scale.
%
%	COLORBAR(H) places the colorbar in the axes H. The
%	colorbar will be horizontal if the axes width > height.
%
%	COLORBAR without arguments either adds a new vertical
%	color scale or updates an existing colorbar.
%
%	H = COLORBAR(...) returns a handle to the colorbar axis.
%
%       clabel  is the string containing the colorbar axis label
%       rlab1   is the string label for the lower limit plotted
%       rlab2   is the string label for the upper limit plotted
%
%

% AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

%return

%set (gca, 'FontSize', 8)
if ~exist('fontsize','var')
    fontsize = 8;
end;
if nargin<1,
  echo ' Error:  vector of plot-axis limits required'
  echo '         else,  use: '
  echo '    colorbar_axis(spectralobject,loc,clabel,rlab1,rlab2)'
  return
end

if nargin>1,
  lower=s.dBlims(1);
  upper=s.dBlims(2);
  paxis = s.dBlims;
end

if nargin<2,
    loc = 'vert';
end

ax = [];
if nargin==2,
  if ~ischar(loc),
    ax = loc;
    if ~strcmp(get(ax,'type'),'axes'),
      error('Requires axes handle.');
    end
    rect = get(ax,'position');
    if rect(3) > rect(4), loc = 'horiz'; else loc = 'vert'; end
  end
end

% Determine color limits by context.  If any axes child is an image
% use scale based on size of colormap, otherwise use current CAXIS.

ch = get(gca,'children');
hasimage = 0; t = [];
for i=1:length(ch),
  if strcmp(get(ch(i),'type'),'image'),
    hasimage = 1;
    t = get(ch(i),'UserData'); % Info stored by imshow or imagesc
  elseif strcmp(get(ch(i),'Type'),'surface'), % Texturemapped surf?
    if strcmp(get(ch(i),'FaceColor'),'texturemap')
      hasimage = 2;
    t = get(ch(i),'UserData'); % Info stored by imshow or imagesc
    end
  end
end
if hasimage,
  if isempty(t),
    t = [0.5 size(colormap,1)+0.5];
  end
else
  t = paxis;
  cmin=t(1);
  cmax=t(2);
end

h = gca;

if nargin==1,
  % Search for existing colorbar
  ch = get(gcf,'children'); ax = [];
  for i=1:length(ch),
    d = get(ch(i),'userdata');
    if numel(d)==1,
        if d==h,
            ax = ch(i);
            pos = get(ch(i),'Position');
            if pos(3)<pos(4), loc = 'vert';
            else loc = 'horiz';
            end
            break;
        end
    end
  end
end

if strcmp(get(gcf,'NextPlot'),'replace'),
  set(gcf,'NextPlot','add')
end
set(gca,'FontSize',fontsize);
fsize=get(gca,'FontSize');          %get fontsize and shift if>16
%vshift=0;
%if fsize>16, vshift=.25; end

%------------------------------------------------------------------------

if loc(1)=='v',        % Append VERTICAL scale to right of current plot
  stripe = 0.075; edge = 0.02;

  if isempty(ax),
    pos = get(h,'Position');
    [az,el] = view;
    if all([az,el]==[0 90]), space = 0.05; else space = .1; end
    set(h,'Position',[pos(1) pos(2) pos(3)*(1-stripe-edge-space) pos(4)])
    rect = [pos(1)+(1-stripe-edge)*pos(3) pos(2) stripe*pos(3)*.5 pos(4)];

    % Create axes for stripe
    ax = axes('Position', rect);
  else
    axes(ax);
  end

  % Create color stripe
  n = size(colormap,1);
  image([0 1],t,[1:n]'); set(ax,'Ydir','normal')

  if nargin>2,
    set(ax,'ylabel',text(0,0,clabel, 'FontSize', fontsize));
  end

  xpos = get(ax,'Xlim');
  ypos = get(ax,'Ylim');
  if nargin>3,
    xshift=.5*(xpos(2)-xpos(1));
    yshift=.05*(ypos(2)-ypos(1));
    if nargin==4,        %put color range  max label on colorbar
      text(xpos(2)-xshift,ypos(2)+yshift,0.,rlab1);
    end
    if nargin==5,        %put color range min, max labels on colorbar
      text(xpos(1)-xshift,ypos(1)-.75*yshift,0.,rlab1);
      text(xpos(1)-xshift,ypos(2)+yshift,0.,rlab2);
    end
  end
  %d=date;
  yshift=.12*(ypos(2)-ypos(1));
%  text(xpos(1)+2.0,ypos(1)-yshift,0.,d)

  % Create color axis
  ylim = get(ax,'ylim');   % Note: axis ticlabel range will be truncated
  units = get(ax,'Units'); set(ax,'Units','pixels');
  pos = get(ax,'Position');
  set(ax,'Units',units);
  yspace = get(ax,'FontSize')*(ylim(2)-ylim(1))/pos(4)/2;
  xspace = .5*get(ax,'FontSize')/pos(3);
  yticks = get(ax,'ytick');
  ylabels = get(ax,'yticklabel');
  labels = []; width = [];

set (gca, 'FontSize', fontsize)  %This sets the font size
  for i=1:length(yticks),
    labels = [labels;text(1+0*xspace,yticks(i),deblank(ylabels(i,:)), ...
         'HorizontalAlignment','right', ...
         'VerticalAlignment','middle', ...
         'FontName',get(ax,'FontName'), ...
         'FontSize',get(ax,'FontSize'), ...
         'FontAngle',get(ax,'FontAngle'), ...
         'FontWeight',get(ax,'FontWeight'))];
    width = [width;get(labels(i),'Extent')];
  end

  % Shift labels over so that they line up
  [dum,k] = max(width(:,3)); width = width(k,3);
  for i=1:length(labels),
    pos = get(labels(i),'Position');
    set(labels(i),'Position',[pos(1)+width pos(2:3)])
  end

  % If we need an exponent then draw one
  [ymax,k] = max(abs(yticks));
  if abs(abs(str2num(ylabels(k,:)))-ymax)>sqrt(eps),
    ex = log10(max(abs(yticks)));
    ex = sign(ex)*ceil(abs(ex));
    l = text(0,ylim(2)+2*yspace,'x 10', ...
         'FontName',get(ax,'FontName'), ...
         'FontSize',get(ax,'FontSize'), ...
         'FontAngle',get(ax,'FontAngle'), ...
         'FontWeight',get(ax,'FontWeight'));
    width = get(l,'Extent');
    text(width(3)-xspace,ylim(2)+3.2*yspace,num2str(ex), ...
         'FontName',get(ax,'ExpFontName'), ...
         'FontSize',get(ax,'ExpFontSize'), ...
         'FontAngle',get(ax,'ExpFontAngle'), ...
         'FontWeight',get(ax,'ExpFontWeight'));
  end

  set(ax,'yticklabelmode','manual','yticklabel','')
  set(ax,'xticklabelmode','manual','xticklabel','')

%set(gca,'ytick',[])
%------------------------------------------------------------------------
%------------------------------------------------------------------------
else            % Append HORIZONTAL scale to bottom of current plot

  if isempty(ax),
    pos = get(h,'Position');           %[left,bottom,width,height]
    stripe = 0.05; space = 0.1;       %stripe = 0.075
    ori=get(gcf,'paperorientation');
    if fsize<=16,
        sfact=1;
    else
        sfact=2;
    end
%    if ori=='landscape', stripe = 0.05; end
%    if fsize<=26,
%      set(h,'Position',...
%        [pos(1) pos(2)+(stripe+space)*pos(4) pos(3) (1-stripe-space)*pos(4)])
%    else
      set(h,'Position',...
     [pos(1) pos(2)+(sfact*stripe+space)*pos(4) pos(3) (1-stripe-space)*pos(4)])
%    end
    rect = [pos(1) pos(2) pos(3) stripe*pos(4)*.8];

    % Create axes for stripe
    ax = axes('Position', rect);
  else
    axes(ax);
  end

  t = paxis;
  cmin=t(1);
  cmax=t(2);

  % Create color stripe
  n = size(colormap,1);
  if paxis(1)>cmin,
    diff=paxis(1)-cmin;
    cstep=(cmax-cmin)/n;
    n1=round(diff/cstep);
  else
    n1=1;
  end
  if paxis(2)<cmax,
    diff=cmax-paxis(2);
    cstep=(cmax-cmin)/n;
    n2=n-round(diff/cstep);
  else
    n2=n;
  end

%  image(t,[0 1],[1:n]); set(ax,'Ydir','normal')
  image(t,[0 1],[n1:n2]); set(ax,'Ydir','normal')
  set(ax,'yticklabelmode','manual')
  set(ax,'yticklabel','')

set (gca, 'FontSize', fontsize)  %This sets the font size
  if nargin>2,
   set(ax,'xlabel',text(0,0,clabel, 'FontSize', fontsize));
  end

  xpos = get(ax,'Xlim');
  ypos = get(ax,'Ylim');
  if nargin>3,
    xshift=.025*(xpos(2)-xpos(1));
    if nargin==4,        %put color range  max label on colorbar
      text(xpos(2)-xshift,ypos(1)-1.0,0.,rlab1);
    end
    if nargin==5,        %put color range min, max labels on colorbar
      text(xpos(1)-xshift,ypos(1)-1.0,0.,rlab1);
      text(xpos(2)-xshift,ypos(1)-1.0,0.,rlab2);
    end
  end
  d=date;
  xshift=.15*(xpos(2)-xpos(1));
%  text(xpos(1)-xshift,ypos(1)-3.5,0.,d)

%set(gca,'xtick',[])
end
%------------------------------------------------------------------------
%------------------------------------------------------------------------

set(ax,'userdata',h)
set(gcf,'CurrentAxes',h)
set(gcf,'Nextplot','Replace')

if nargout>0, handle = ax; end

author: geoscience-community-codes

repo: GISMO


text:
function [beta, p] = logist2Fit(y, x, addOne, w)
% LOGIST2FIT 2 class logsitic classification
% function beta = logist2Fit(y,x, addOne)
%
% y(i) = 0/1
% x(:,i) = i'th input - we optionally append 1s to last dimension
% w(i) = optional weight
%
% beta(j)- regression coefficient

if nargin < 3, addOne = 1; end
if nargin < 4, w = 1; end

Ncases = size(x,2);
if Ncases ~= length(y)
  error(sprintf('size of data = %dx%d, size of labels=%d', size(x,1), size(x,2), length(y)))
end
if addOne
  x = [x; ones(1,Ncases)];
end
[beta, p] = logist2(y(:), x', w(:));
beta = beta(:);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function z = mldivide( x, y )

%   Disciplined convex/geomtric programming information for MLDIVIDE:
%      The MLDIVIDE operation X\Y can be employed with X as a scalar. In
%      that case, it is equivalent to the LDIVIDE operation X.\Y, and
%      must obey the same rules as outlined in the help for CVX/LDIVIDE.
%
%      When X is a matrix, the MRDIVIDE operation X\Y is equivalent to
%      inv(X)*Y for both DCP and DGP purposes. The inv() operation is
%      not supported for non-constant expressions, so X must be both
%      constant and nonsingular. The resulting matrix multiplication
%      must obey the same rules as outlined in the help for CVX/MTIMES.

sz = size( x );
if all( sz == 1 ),
    z = ldivide( x, y, '\' );
elseif length( sz ) > 2,
    cvx_throw( 'Inputs must be 2-D, or at least one input must be scalar.' );
elseif sz( 1 ) ~= sz( 2 ),
    cvx_throw( 'Non-square matrix divisors are not supported in CVX.' );
elseif ~cvx_isconstant( x ),
    cvx_throw( 'Matrix divisors must be constant.' );
else
    z = mtimes( cvx_constant( x ), y, '\' );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
function b = transpose(a)
% function B=transpose(A)
%
% DESCRIPTION
%   Returns the non-conjugate transpose of a polynomial.
%
% INPUTS
%   A: polynomial
%
% OUTPUTS
%   B: transpose of A.
%
% SYNTAX
%   B= transpose(A);
%   B= A.';
%
% See also ctranspose

% 6/7/2002: PJS  Initial Coding

b=a;

% Get Dimensions
[nrb,ncb]=size(b);
bcoef = b.coefficient;

idx = [];
tempidx = 0:nrb:(ncb-1)*nrb;
for i1 = 1:nrb
    idx = [idx i1+tempidx];
end
bcoef = bcoef(:,idx);
b.coefficient = bcoef;
b.matdim = [ncb nrb];


text:
function X = vgg_intersect_quadrics(A, B, C);
%
% function X = vgg_intersect_quadrics(A, B, C);
%
% Purpose:
%   Compute intersections of three quadrics in projective 3-space.
%
% Input:
%   A, B, C are symmetric 4x4 matrices.
%
% Output:
%   X is 4-by-8 complex. Each column is a root.
%
% The function works by converting the problem into an 8x8 generalized
% eigen-problem.
%
% fsm@robots.ox.ac.uk
%

if any(size(A) ~= [4 4])
  error('bad A')
end
if any(size(B) ~= [4 4])
  error('bad B')
end
if any(size(C) ~= [4 4])
  error('bad C')
end

% symmetrize and normalize.
A = A + A.'; A = A / norm(A, 'fro'); a = [A(1, 1) 2*A(1, 2) 2*A(1, 3) 2*A(1, 4) A(2, 2) 2*A(2, 3) 2*A(2, 4) A(3, 3) 2*A(3, 4) A(4, 4)].';
B = B + B.'; B = B / norm(B, 'fro'); b = [B(1, 1) 2*B(1, 2) 2*B(1, 3) 2*B(1, 4) B(2, 2) 2*B(2, 3) 2*B(2, 4) B(3, 3) 2*B(3, 4) B(4, 4)].';
C = C + C.'; C = C / norm(C, 'fro'); c = [C(1, 1) 2*C(1, 2) 2*C(1, 3) 2*C(1, 4) C(2, 2) 2*C(2, 3) 2*C(2, 4) C(3, 3) 2*C(3, 4) C(4, 4)].';

% make multiplication tables (for degrees 1+2=3, 2+2=4 and 3+1=4)
if 1
  ijk = [
    1
    22
    43
    64
    82
    105
    126
    147
    163
    186
    208
    229
    244
    267
    289
    310
    325
    351
    372
    393
    406
    432
    454
    475
    487
    513
    535
    556
    568
    594
    617
    638
    649
    675
    698
    719
    730
    756
    779
    800
    ];
  X12 = zeros(20, 4, 10); X12(ijk) = 1;

  ijk = [
    1
    37
    73
    109
    145
    181
    217
    253
    289
    325
    352
    390
    426
    462
    501
    537
    573
    609
    645
    681
    703
    741
    778
    814
    852
    889
    925
    962
    998
    1034
    1054
    1092
    1129
    1165
    1203
    1240
    1276
    1313
    1349
    1385
    1405
    1446
    1482
    1518
    1561
    1597
    1633
    1669
    1705
    1741
    1756
    1797
    1834
    1870
    1912
    1949
    1985
    2022
    2058
    2094
    2107
    2148
    2185
    2221
    2263
    2300
    2336
    2373
    2409
    2445
    2458
    2499
    2537
    2573
    2614
    2652
    2688
    2726
    2762
    2798
    2809
    2850
    2888
    2924
    2965
    3003
    3039
    3077
    3113
    3149
    3160
    3201
    3239
    3275
    3316
    3354
    3390
    3428
    3464
    3500
    ];
  X22 = zeros(35, 10, 10); X22(ijk) = 1;

  ijk = [
    1
    37
    73
    109
    145
    181
    217
    253
    289
    325
    361
    397
    433
    469
    505
    541
    577
    613
    649
    685
    702
    740
    776
    812
    851
    887
    923
    959
    995
    1031
    1071
    1107
    1143
    1179
    1215
    1251
    1287
    1323
    1359
    1395
    1403
    1441
    1478
    1514
    1552
    1589
    1625
    1662
    1698
    1734
    1772
    1809
    1845
    1882
    1918
    1954
    1991
    2027
    2063
    2099
    2104
    2142
    2179
    2215
    2253
    2290
    2326
    2363
    2399
    2435
    2473
    2510
    2546
    2583
    2619
    2655
    2692
    2728
    2764
    2800
    ];
  X31 = zeros(35, 20, 4); X31(ijk) = 1;

  clear ijk
end

% compose ideal in degree 2.
I2 = [a b c];

% compute ideal in degrees 3 and 4.
I3 = reshape(reshape(X12, [20* 4 10]) * I2, [20 3* 4]);
I4 = reshape(reshape(X22, [35*10 10]) * I2, [35 3*10]);

% the columns of C3 and C4 span complements of I3 and I4, respectively.
[U,S,V] = svd(I3.'); C3 = V(:, 13:20); %W3 = diag(S).'
[U,S,V] = svd(I4.'); C4 = V(:, 28:35); %W4 = diag(S).'

% make the four multiplication operators.
M = cell(1, 4);
for k=1:4
  M{k} = C4.' * X31(:, :, k) * C3;
end

% FIXME: here we only use M{1} and M{2} to compute the Schur
% decomposition although we use all four multiplication operators to
% extract the root coordinates.
[AA, BB, Q, Z] = qz(M{1}, M{2});
for k=1:4
  M{k} = Q*M{k}*Z;
end
%M{:}

X = zeros(4, 8);
X(1, :) = diag(M{1}).';
X(2, :) = diag(M{2}).';
X(3, :) = diag(M{3}).';
X(4, :) = diag(M{4}).';

% normalize
X = X ./ repmat(max(abs(X), [], 1), [4 1]);

return;

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function [x,err,it] = perform_conjugate_gradient(A,y,options)

% perform_conjugate_gradient - perform (bi)-conjugate gradient
%
%     [x,err,k] = perform_conjugate_gradient(A,y,options);
%
%   Solves for A*x=y.
%   Works both for vector x,y and matrix x,y (parralel soving of many
%   linear equations).
%
%   Important: the algorithm assumes that the matrix is symmetric definite
%   positive. If it is not the case, then you must set options.is_sdp=0,
%   and the algorithm will use a bi-conjugate gradient descent (x2 slower).
%
%   A can be a matrix or a callback function y=A(x,options).
%   In this case, you have to set options.ncols as the number of columns of
%   A (if it is a callback).
%   In the non-symmetric case, then options.transpose=1 for the computation
%   of A*x and options.transpose=-1 for the computation of A'*x.
%
%   err monitors the decreasing of |A*x-y|.
%   k is the total number of iterations.
%
%   You can set:
%       options.x is an initial guess
%       options.epsilon is maximum error
%       options.niter_max is maximum number of iterations
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
niter = getoptions(options, 'niter_max', 100);
epsilon = getoptions(options, 'epsilon', 1e-5);
is_sdp = getoptions(options, 'is_sdp', 1);

use_callback = 0;
if not(isnumeric(A))
    use_callback = 1;
end
if isfield(options, 'x')
    x = options.x;
else
    if use_callback==0
        x = zeros(size(A,2),1);
    else
        if isfield(options, 'ncols')
            x = zeros(options.ncols, 1);
        else
            error('You have to specify options.ncols');
        end
    end
end

if norm(y(:), 'fro') == 0
    normb = epsilon;
else
    normb = epsilon; % * sum(y(:).^2);
end

if use_callback==0
    r  = y - A*x;
else
    options.transpose = 1;
    r  = y - feval(A,x,options);
end
p = r;
if is_sdp==0
    rr = r;
    pp = p;
end
r0 = sum(r(:).^2);

err = [sum(r0)];
for it=1:niter

    % auxiliary vector
    if use_callback==0
        w  = A*p;
    else
        options.transpose = 1;
        % options.transpose = -1;
        w  = feval(A,p,options);
    end

    if is_sdp==1
        d = sum(p(:) .* w(:));
    else
        d = sum(pp(:) .* w(:));
    end
    I = find(abs(d)<eps); d(I) = 1;
    alpha = repmat( r0 ./ d, [size(x,1) 1] );              % found optimal alpha in line search
    x = x + alpha.*p;                       % new guess
    r = r - alpha.*w;                       % the residual is in fact r=b-A*x
    if is_sdp==0
        if use_callback==0
            w  = A'*pp;
        else
            options.transpose = -1;
            w  = feval(A,pp,options);
        end
        rr = rr-alpha.*w;
    end

    rprev = r0;                             % save norm of the old residual
    if is_sdp==0
        r0 = sum(r.*rr);
    else
        r0 = sum(r.^2);                         % compute norm of new residual
    end

    err(end+1) = sqrt( sum(r.^2) );
    if err(end)<normb
        return;
    end

    % search direction
    beta = r0./rprev;
    p = r + repmat(beta, [size(x,1) 1]).*p;
    if is_sdp==0
        pp = rr + repmat(beta, [size(x,1) 1]).*pp;
    end

end

author: gpeyre

repo: matlab-toolboxes


text:
classdef OptimizerAugmentedLagrangian < Optimizer

    properties (GetAccess = public, SetAccess = protected)
        type = 'Augmented Lagrangian';
    end

    properties (Access = private)
        tau
        lineSearchTrials
        lineSearch
        costOld
        upperBound
        lowerBound
        tol = 1e-8
        nX
        nConstr
        hasConverged
        acceptableStep
        oldDesignVariable
        oldCost
        incrementalScheme
        hasFinished
        mOld
        meritNew
        penalty
        meritGradient

        globalCost
        globalConstraint
        globalCostGradient
        globalMerit
        globalLineSearch
        globalDual
        globalDesignVar
    end

    methods (Access = public)

        function obj = OptimizerAugmentedLagrangian(cParams)
            obj.initOptimizer(cParams);
            obj.init(cParams);
            obj.outputFunction.monitoring.create(cParams);
            obj.createPrimalUpdater(cParams);
            obj.createDualUpdater(cParams);
            obj.prepareFirstIter();
        end

        function obj = solveProblem(obj)
            obj.hasConverged = false;
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();
%             obj.saveVariablesForAnalysis();
            obj.hasFinished = 0;
            obj.printOptimizerVariable();
            while ~obj.hasFinished
%             while ~obj.hasConverged
                obj.update();
                obj.updateIterInfo();
                obj.updateMonitoring();
                obj.checkConvergence();
                obj.printOptimizerVariable();
%                 obj.saveVariablesForAnalysis();
            end
        end

    end

    methods(Access = private)

        function init(obj,cParams)
            obj.upperBound             = cParams.uncOptimizerSettings.ub;
            obj.lowerBound             = cParams.uncOptimizerSettings.lb;
            obj.cost                   = cParams.cost;
            obj.constraint             = cParams.constraint;
            obj.nConstr                = cParams.constraint.nSF;
            obj.designVariable         = cParams.designVar;
            obj.dualVariable           = cParams.dualVariable;
            obj.incrementalScheme      = cParams.incrementalScheme;
            obj.nX                     = length(obj.designVariable.value);
            obj.maxIter                = cParams.maxIter;
            obj.hasConverged           = false;
            obj.nIter                  = 0;
        end

        function prepareFirstIter(obj)
            obj.cost.computeFunctionAndGradient();
            obj.costOld = obj.cost.value;
            obj.designVariable.updateOld();
            obj.dualVariable.value = zeros(obj.nConstr,1);
            obj.penalty            = 10;
        end

        function obj = update(obj)
            x0 = obj.designVariable.value;
            obj.designVariable.update(x0);
            obj.saveOldValues(x0);
            obj.mOld = obj.computeMeritFunction(x0);
            obj.calculateInitialStep();
            obj.acceptableStep   = false;
            obj.lineSearchTrials = 0;
            obj.computeMeritGradient();
            while ~obj.acceptableStep
                x = obj.updatePrimal();
                obj.checkStep(x,x0);
            end
            obj.updateOldValues(x);
        end

        function displayIter(obj,x)
            m = obj.designVariable.mesh;
            bm = m.createBoundaryMesh();
            s.backgroundMesh = m;
            s.boundaryMesh   = bm;
            um = UnfittedMesh(s);
            um.compute(x);
            figure()
            um.plot();
        end

        function obj = calculateInitialStep(obj)
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();
            x       = obj.designVariable.value;
            l       = obj.dualVariable.value;
            DJ      = obj.cost.gradient;
            Dg      = obj.constraint.gradient;
            g       = obj.constraint.value;
            p       = obj.penalty;
            DmF     = DJ + Dg*(l + p*g);
            if obj.nIter == 0
                factor = 1;
                obj.primalUpdater.computeFirstStepLength(DmF,x,factor);
            else
                factor = 1.05;
                obj.primalUpdater.increaseStepLength(factor);
            end
        end

        function x = updatePrimal(obj)
            x   = obj.designVariable.value;
            g   = obj.meritGradient;
            x   = obj.primalUpdater.update(g,x);
        end

        function computeMeritGradient(obj)
            Dh    = obj.constraint.gradient;
            DJ    = obj.cost.gradient;
            l     = obj.dualVariable.value;
            p     = obj.penalty;
            gPlus = obj.defineConstraintValue();
            g     = (DJ + Dh*(l + p*gPlus));
            obj.meritGradient = g;
        end

        function mF = computeMeritFunction(obj,x)
            obj.designVariable.update(x)
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();
            J      = obj.cost.value;
            gPlus  = obj.defineConstraintValue();
            l      = obj.dualVariable.value;
            rho    = obj.penalty;
            mF     = J + l'*gPlus + 0.5*rho*(gPlus'*gPlus);
        end

        function c = defineConstraintValue(obj)
            c   = obj.constraint.value;
            l   = obj.dualVariable.value;
            rho = obj.penalty;
            for i = 1:obj.nConstr
                switch obj.constraintCase{i}
                    case 'EQUALITY'

                    case 'INEQUALITY'
                        c(i) = max(c(i),-l/rho);
                end
            end
        end

        function checkStep(obj,x,x0)
            mNew = obj.computeMeritFunction(x);
            if mNew < obj.mOld
                obj.acceptableStep = true;
                obj.dualUpdater.updatePenalty(obj.penalty);
                obj.dualUpdater.update();
                obj.meritNew = mNew;
            elseif obj.primalUpdater.isTooSmall()
%                 error('Convergence could not be achieved (step length too small)')
                warning('Convergence could not be achieved (step length too small)')
                obj.acceptableStep = true;
                obj.meritNew = mNew; % Provisional value
            else
                obj.primalUpdater.decreaseStepLength();
                obj.designVariable.update(x0);
                obj.lineSearchTrials = obj.lineSearchTrials + 1;
            end
        end

        function obj = saveOldValues(obj,x)
            obj.designVariable.update(x);
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();
            obj.oldCost            = obj.cost.value;
            obj.oldDesignVariable  = x;
        end

        function obj = updateOldValues(obj,x)
            obj.designVariable.update(x);
            obj.cost.computeFunctionAndGradient();
            obj.constraint.computeFunctionAndGradient();
        end

        function obj = checkConvergence(obj)
           if abs(obj.meritNew - obj.mOld) < obj.tol && obj.checkConstraint()
               obj.hasConverged = true;
           else

           end

        end

        function obj = updateMonitoring(obj)
            s.nIter            = obj.nIter;
            s.tau              = obj.primalUpdater.tau;
            s.lineSearch       = obj.lineSearch;
            s.lineSearchTrials = obj.lineSearchTrials;
            s.oldCost          = obj.oldCost;
            s.hasFinished      = obj.hasFinished;
            s.meritNew         = obj.meritNew;
            obj.outputFunction.monitoring.compute(s);
        end

        function updateIterInfo(obj)
            obj.increaseIter();
            obj.updateStatus();
        end

        function increaseIter(obj)
            obj.nIter = obj.nIter + 1;
        end

        function updateStatus(obj)
            obj.hasFinished = obj.hasConverged || obj.hasExceededStepIterations();
        end

        function itHas = hasExceededStepIterations(obj)
            iStep = obj.incrementalScheme.iStep;
            nStep = obj.incrementalScheme.nSteps;
            itHas = obj.nIter >= obj.maxIter*(iStep/nStep);
        end

        function saveVariablesForAnalysis(obj)
            i                           = obj.nIter + 1;
            obj.globalCost(i)           = obj.cost.value;
            obj.globalConstraint(:,i)   = obj.constraint.value;
            obj.globalCostGradient(i)   = norm(obj.cost.gradient);
            obj.globalMerit(i)          = obj.meritNew;
            obj.globalLineSearch(i)     = obj.primalUpdater.tau;
            obj.globalDual(:,i)         = obj.dualVariable.value;
            obj.globalDesignVar(:,i)    = obj.designVariable.value;
            if obj.hasConverged
                c = obj.globalCost;
                h = obj.globalConstraint;
                g = obj.globalCostGradient;
                m = obj.globalMerit;
                t = obj.globalLineSearch;
                d = obj.globalDual;
                v = obj.globalDesignVar;
                save('name.mat',"c","g","h","d","v");
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function obj = setInitialValue(obj, x)
    % Sets a typical initial value for the NLP variables. This value will
    % be used to assign the initial value for the NLP solver.
    %
    % @note If the initial value is not given, then use the middle value of
    % the lower/upper boundaries.
    %
    % Parameters:
    % x: an array of initial value of the NLP variables @type colvec

    % determine the typical value
    if nargin > 1

        if isscalar(x)
            x = x*ones(obj.Dimension,1);
        else
            if isrow(x)
                x = transpose(x);
            end
            validateattributes(x,{'double'},...
                {'vector','numel',obj.Dimension,'real'},...
                'NlpVariable.setInitialValue','InitialValue');
        end
        obj.InitialValue = x;
    else
        % preallocate

        lb_tmp = obj.LowerBound;
        ub_tmp = obj.UpperBound;

        % replace infinity with very high numbers
        lb_tmp(lb_tmp==-inf) = -1e5;
        ub_tmp(ub_tmp==inf)  = 1e5;

        obj.InitialValue = (ub_tmp + lb_tmp)/2;
    end

end

author: ayonga

repo: frost-dev

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


text:
function geometry_test0757 ( )

%*****************************************************************************80
%
%% TEST0757 tests POLYGON_ANGLES_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;
  n = 6;

  v = [ ...
    0.0, 0.0; ...
    1.0, 0.0; ...
    2.0, 1.0; ...
    3.0, 0.0; ...
    3.0, 2.0; ...
    1.0, 2.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0757\n' );
  fprintf ( 1, '  For a polygon in 2D:\n' );
  fprintf ( 1, '  POLYGON_ANGLES_2D computes the angles.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of polygonal vertices = %d\n', n );

  r8mat_transpose_print ( dim_num, n, v, '  The polygon vertices:' );

  angle = polygon_angles_2d ( n, v );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Polygonal angles in degrees:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %6d  %14f\n', i, radians_to_degrees ( angle(i) ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fgplvmClassVisualise(call)

% FGPLVMCLASSVISUALISE Callback function for visualising data in 2-D.

% FGPLVM

global visualiseInfo

switch call
 case 'click'
  visualiseInfo.clicked = ~visualiseInfo.clicked;
 case 'move'
  if visualiseInfo.clicked
    [x, y]  = localCheckPointPosition(visualiseInfo);
    if ~isempty(x)
      % This should be changed to a model specific visualisation.
      set(visualiseInfo.latentHandle, 'xdata', x, 'ydata', y);
      fhandle = str2func([visualiseInfo.model.type 'PosteriorMeanVar']);
      if str2num(version('-release'))>13
        [mu, varsigma] = fhandle(visualiseInfo.model, [x y]);
      else
        [mu, varsigma] = feval(fhandle, visualiseInfo.model, [x y]);
      end
      if isfield(visualiseInfo.model, 'noise')
        Y = noiseOut(visualiseInfo.model.noise, mu, varsigma);
      else
        Y = mu;
      end
%      if str2num(version('-release'))>13
%        visualiseInfo.visualiseModify(visualiseInfo.visHandle, Y, ...
%                                      visualiseInfo.varargin{:});
%      else
        feval(visualiseInfo.visualiseModify, visualiseInfo.visHandle, ...
              Y, visualiseInfo.varargin{:});
%      end
      visualiseInfo.latentPos=[x, y];
    end
  end
end

function point = localGetNormCursorPoint(figHandle)

point = get(figHandle, 'currentPoint');
figPos = get(figHandle, 'Position');
% Normalise the point of the curstor
point(1) = point(1)/figPos(3);
point(2) = point(2)/figPos(4);

function [x, y] = localGetNormAxesPoint(point, axesHandle)

position = get(axesHandle, 'Position');
x = (point(1) - position(1))/position(3);
y = (point(2) - position(2))/position(4);
lim = get(axesHandle, 'XLim');
x = x*(lim(2) - lim(1));
x = x + lim(1);
lim = get(axesHandle, 'YLim');
y = y*(lim(2) - lim(1));
y = y + lim(1);

function [x, y] = localCheckPointPosition(visualiseInfo)

% Get the point of the cursor
point = localGetNormCursorPoint(gcf);

% get the position of the axes
position = get(visualiseInfo.plotAxes, 'Position');

% Check if the pointer is in the axes
if point(1) > position(1) ...
      & point(1) < position(1) + position(3) ...
      & point(2) > position(2) ...
      & point(2) < position(2) + position(4);

  % Rescale the point according to the axes
  [x y] = localGetNormAxesPoint(point, visualiseInfo.plotAxes);

  % Find the nearest point
else
  % Return nothing
  x = [];
  y = [];
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function out = isempty(F)
%ISEMPTY   Empty boolean check for a SPHEREFUNV object.
%   ISEMPTY(F) returns 1 if every component of F is an empty SPHEREFUN, and
%   return 0 otherwise.

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

if ( isempty( F.components) )
    out = 1;
    return
end

% Take isempty of each component:
out = cellfun(@isempty, F.components, 'UniformOutput', false);
out = all(cell2mat(out));

end


text:
function [dx, dy] = evalshift(p, pval, L1, L2, L3)

% function [dx, dy] = evalshift(p, pval, L1, L2, L3)
% Purpose: compute two-dimensional Warp & Blend transform

% 1) compute Gauss-Lobatto-Legendre node distribution
gaussX = -JacobiGL(0,0,p);

% 2) compute blending function at each node for each edge
blend1 = L2.*L3; blend2 = L1.*L3; blend3 = L1.*L2;

% 3) amount of warp for each node, for each edge
warpfactor1 = 4*evalwarp(p, gaussX, L3-L2);
warpfactor2 = 4*evalwarp(p, gaussX, L1-L3);
warpfactor3 = 4*evalwarp(p, gaussX, L2-L1);

% 4) combine blend & warp
warp1 = blend1.*warpfactor1.*(1 + (pval*L1).^2);
warp2 = blend2.*warpfactor2.*(1 + (pval*L2).^2);
warp3 = blend3.*warpfactor3.*(1 + (pval*L3).^2);

% 5) evaluate shift in equilateral triangle
dx = 1*warp1 + cos(2*pi/3)*warp2 + cos(4*pi/3)*warp3;
dy = 0*warp1 + sin(2*pi/3)*warp2 + sin(4*pi/3)*warp3;
return;

repo: Aero-matlab


text:
classdef PrimalUpdaterFactory < handle

    methods (Access = public, Static)

        function p = create(cParams)
            switch cParams.optimizerNames.primal
                case 'SLERP'
                    p = SLERP(cParams);
                case 'PROJECTED GRADIENT'
                    p = ProjectedGradient(cParams);
                case 'HAMILTON-JACOBI'
                    p = HamiltonJacobi(cParams);
            end
        end

    end

end


text:
function inside = circle_sector_contains_point_2d ( r, center, theta1, theta2, p )

%*****************************************************************************80
%
%% CIRCLE_SECTOR_CONTAINS_POINT_2D : is a point inside a circular sector?
%
%  Discussion:
%
%    A circular sector is formed by a circular arc, and the two straight line
%    segments that join its ends to the center of the circle.
%
%    A circular sector is defined by
%
%      ( X - CENTER(1) )**2 + ( Y - CENTER(2) )**2 = R**2
%
%    and
%
%      Theta1 <= Theta <= Theta2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the circle.
%
%    Input, real CENTER(2,1), the center of the circle.
%
%    Input, real THETA1, THETA2, the angles defining the arc,
%    in radians.  Normally, THETA1 < THETA2.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, logical INSIDE, is TRUE if the point is inside or on the
%    circular sector, FALSE otherwise.
%
  inside = 0;
%
%  Is the point inside the (full) circle?
%
  if ( ( p(1,1) - center(1,1) ) * ( p(1,1) - center(1,1) ) ...
     + ( p(2,1) - center(2,1) ) * ( p(2,1) - center(2,1) ) <= r * r )
%
%  Is the point's angle within the arc's range?
%  Try to force the angles to lie between 0 and 2 * PI.
%
    theta = r8_atan ( p(2,1) - center(2,1), p(1,1) - center(1,1) );

    if ( r8_modp ( theta  - theta1,  2.0 * pi ) <= ...
         r8_modp ( theta2 - theta1,  2.0 * pi ) )

      inside = 1;

    end

  end

  return
end


text:
function varargout = createIcosahedron()
%CREATEICOSAHEDRON Create a 3D mesh representing an Icosahedron.
%
%   MESH = createIcosahedron;
%   [V, E, F] = createIcosahedron;
%   Create a solid with 12 vertices, and 20 triangular faces. Faces are
%   oriented outwards of the mesh.
%
%   [V, F] = createIcosahedron;
%   Returns only the vertices and the face vertex indices.
%
%   MESH = createIcosahedron;
%   Returns the data as a mesh structure, with fields 'vertices', 'edges'
%   and 'faces'.
%
%   Example
%     [n, e, f] = createIcosahedron;
%     drawMesh(n, f);
%
%   See also
%   meshes3d, drawMesh
%   createCube, createOctahedron, createDodecahedron, createTetrahedron
%

%   ---------
%   author: David Legland
%   mail: david.legland@inra.fr
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/03/2005.
%

%   HISTORY
%   2007-01-04 remove unused variables
%   2010-12-06 format output, orient normals outwards

%% Initialisations

theta = 2*pi/5;
l = 1/sin(theta/2)/2;
z1 = sqrt(1-l*l);

t1 = (0:2*pi/5:2*pi*(1-1/5))';
x1 = l*cos(t1);
y1 = l*sin(t1);

t2 = t1 + 2*pi/10;
x2 = l*cos(t2);
y2 = l*sin(t2);

h = sqrt(l*l-.5*.5);
z2 = sqrt(3/4 - (l-h)*(l-h));

%% Create mesh data

nodes = [0 0 0;...
    [x1 y1 repmat(z1, [5 1])]; ...
    [x2 y2 repmat(z1+z2, [5 1])]; ...
    0 0 2*z1+z2];

edges = [...
    1 2;1 3;1 4;1 5;1 6; ...
    2 3;3 4;4 5;5 6;6 2; ...
    2 7;7 3;3 8;8 4;4 9;9 5;5 10;10 6;6 11;11 2; ...
    7 8;8 9;9 10;10 11;11 7; ...
    7 12;8 12;9 12;10 12;11 12];

% faces are ordered to have normals pointing outside of the mesh
faces = [...
    1 3  2 ; 1 4  3 ; 1  5  4 ;  1  6  5 ;  1 2  6;...
    2 3  7 ; 3 4  8 ; 4  5  9 ;  5  6 10 ;  6 2 11;...
    7 3  8 ; 8 4  9 ; 9  5 10 ; 10  6 11 ; 11 2  7;...
    7 8 12 ; 8 9 12 ; 9 10 12 ; 10 11 12 ; 11 7 12];

% format output
varargout = formatMeshOutput(nargout, nodes, edges, faces);

author: rpng

repo: lips


text:
function y = sfb2D_A(lo, hi, sf, d)

% 2D Synthesis Filter Bank
% (along single dimension only)
%
% y = sfb2D_A(lo, hi, sf, d);
% sf - synthesis filters
% d  - dimension of filtering
% see afb2D_A

lpf = sf(:, 1);     % lowpass filter
hpf = sf(:, 2);     % highpass filter

if d == 2
   lo = lo';
   hi = hi';
end

N = 2*size(lo,1);
L = length(sf);
y = upfirdn(lo, lpf, 2, 1) + upfirdn(hi, hpf, 2, 1);
y(1:L-2, :) = y(1:L-2, :) + y(N+[1:L-2], :);
y = y(1:N, :);
y = cshift2D(y, 1-L/2);

if d == 2
   y = y';
end

author: thomaskuestner


text:
function test_bug1408

% MEM 4gb
% WALLTIME 00:10:00
% DEPENDENCY ft_preprocessing preproc ft_preproc_bandpassfilter ft_preproc_bandstopfilter ft_preproc_dftfilter ft_preproc_highpassfilter ft_preproc_lowpassfilter

nchans   = 200;
nsamples = 1e6;

dat = rand(nchans,nsamples); datmean = rand(nchans,1);
tic
dat = dat - datmean(:,ones(1,nsamples));
t(1) = toc;

dat = rand(nchans,nsamples); datmean = rand(nchans,1);
tic
dat = dat - repmat(datmean,[1 nsamples]);
t(2) = toc;

dat = rand(nchans,nsamples); datmean = rand(nchans,1);
tic
for ichan = 1:nchans
  dat(ichan,:) = dat(ichan,:) - datmean(ichan);
end
t(3) = toc;

dat = rand(nchans,nsamples); datmean = rand(nchans,1);
tic
dat = dat';
for ichan = 1:nchans
  dat(:,ichan) = dat(:,ichan) - datmean(ichan);
end
dat = dat';
t(4) = toc;

dat = rand(nchans,nsamples); datmean = rand(nchans,1);
tic
dat = bsxfun(@minus, dat, datmean);
t(5) = toc;

dat = rand(nchans,nsamples); datmean = rand(nchans,1);
tic
for isample = 1:nsamples
  dat(:,isample) = dat(:, isample) - datmean;
end
t(6) = toc;

[minval, minindx] = min(t);
if minindx~=6
  warning('unexpected winner of the speed test');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data          = [];
data.label    = {'1', '2', '3'};
data.trial{1} = rand(3,1000);
data.time{1}  = (0:999)/1000;
data.fsample  = 1000;

cfg = [];
cfg.demean = 'yes';
dataout = ft_preprocessing(cfg, data);

cfg = [];
cfg.detrend = 'yes';
dataout = ft_preprocessing(cfg, data);

cfg = [];
cfg.polyremoval = 'yes';
cfg.polyorder = 3;
dataout = ft_preprocessing(cfg, data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dat = rand(10,300)+ 30;
Fs  = 1000;
Fl  = 50;
Flp = 10;
Fhp = 70;
Fbp = [40 70];

filt = ft_preproc_bandpassfilter(dat,Fs,Fbp); assert(all(mean(filt,2)<1)); % the DC component should disappear
filt = ft_preproc_bandstopfilter(dat,Fs,Fbp); assert(all(mean(filt,2)>0));
filt = ft_preproc_dftfilter(dat,Fs,Fl);       assert(all(mean(filt,2)>0));
filt = ft_preproc_highpassfilter(dat,Fs,Fhp); assert(all(mean(filt,2)<1)); % the DC component should disappear
filt = ft_preproc_lowpassfilter(dat,Fs,Flp);  assert(all(mean(filt,2)>0));


text:
function [result] = VBA_susceptibility(posterior,out,ind)
% scores the susceptibility of u->y relationships w.r.t. model parameters
% [result] = VBA_susceptibility(posterior,out,ind)
% IN:
%   - out/posterior: the output structures of the VBA inversion
%   - ind: a structure containing the following fields:
%       .u  : indices of inputs of interests
%       .phi: indices of observation parameters, whose importance in the
%       u->y relationship one wants to score
%       .theta: idem, for evolution parameters.
%       .y : indices of observations of interest
% OUT:
%   - susceptibility: nuXnphi and nuxntheta relative susceptibility
%   matrices (w.r.t. observation/evolution parameters), where nu and
%   nphi/ntehta are the number of considered inputs and parameters, respectively.
%   - susceptibility: nuXntphi qnd nuxntheta relative susceptibility matrices
%   (w.r.t. observation/evolution  parameters).

% ========================================================================================
% check parameters
% ========================================================================================

dim = out.dim;

% ________________________________________________________________________________________
% check if models has inputs

if dim.u == 0
    disp('*** Susceptibility analysis: the system does not receive any input!')
    result.susceptibility = struct('phi',[],'theta',[]);
    result.specificity    = struct('phi',[],'theta',[]);
    return
end

% ________________________________________________________________________________________
% set default parameters of interest

if ~exist('ind'), ind = struct; end

% check if is DCM
try, isDCM=out.options.inF.fullDCM; catch, isDCM=0; end

% complete option strucutre
if isDCM
    % keep only DCM connections
   idxConnect = 1:(out.options.inF.indself-1);
   % idxConnect = setdiff(idxConnect,out.options.inF.indC); % trivial effects (or maybe not)
   if isfield(out.options.inF,'indhself')
      % behaviour DCM
      idxResp = [out.options.sources(2:end).out] ;
   else
      % normal DCM
      idxResp = 1:dim.p;
   end
   ind = VBA_check_struct(ind,        ...
        'u'     , 1:dim.u       , ...
        'phi'   , []            , ...
        'theta' , idxConnect    , ...
        'y'     , idxResp         ...
    );
else
    % all params
    ind = VBA_check_struct(ind,       ...
        'u'     , 1:dim.u       , ...
        'phi'   , 1:dim.n_phi   , ...
        'theta' , 1:dim.n_theta , ...
        'y'     , 1:dim.p         ...
    );
end

% get rid of fixed parameters
fixedTheta = find(diag(out.options.priors.SigmaTheta)==0);
ind.theta = setdiff(ind.theta,fixedTheta);
fixedPhi = find(diag(out.options.priors.SigmaPhi)==0);
ind.phi = setdiff(ind.phi,fixedPhi);

% get rid of session param
if isfield(out.options, 'multisession')
    ind.u = setdiff(ind.u,dim.u);
end
% ========================================================================================
% shortcuts
% ========================================================================================

%options = out.options;
%f_fname = options.f_fname;
%g_fname = options.g_fname;
%y = out.y;

% ========================================================================================
% initialization
% ========================================================================================

% ========================================================================================
% Mutilating the models and deriving the output distortion
% ========================================================================================

% ________________________________________________________________________________________
% 1 - switching off inputs
vu = zeros(length(ind.u),1);
parfor iu=1:length(ind.u)
    mutilated_input_idx = ind.u(iu);
    vu(iu) = mutilation_score(posterior,out,ind,'u',mutilated_input_idx);
end

% dealing with parameters
for paramType = {'phi','theta'}
    paramType = paramType{1} ;

    if ~isempty(ind.(paramType)) > 0

        % compute mutilation score
        vparam  = zeros(1,numel(ind.(paramType)));
        % ________________________________________________________________________________
        % 2 - switching off system parameters
        parfor ip=1:length(ind.(paramType))
            mutilated_param_idx = ind.(paramType)(ip);
            vparam(ip) = mutilation_score(posterior,out,ind,paramType,mutilated_param_idx);
        end

        % _______________________________________________________________________________
        % 3 - bilateral mutilations
        vuparam = zeros(numel(ind.u),numel(ind.(paramType)));
        for iu=1:length(ind.u)
           parfor ip=1:length(ind.(paramType))
               mutilated_param_idx = ind.(paramType)(ip);
               mutilated_input_idx = ind.u(iu);
               vuparam_iu(ip) = mutilation_score(posterior,out,ind, ...
               paramType,mutilated_param_idx,'u',mutilated_input_idx);
           end
           vuparam(iu,:) = vuparam_iu ;
        end

        % ________________________________________________________________________________
        % compute susceptibility and specificity
        dvuparam = vuparam;
        dvu      = repmat(vu,1,length(ind.(paramType)));
        dvparam  = repmat(vparam,length(ind.u),1);

        % relative susceptibility
        susceptibility.raw.(paramType)  = 1 + (dvparam - dvuparam)      ;
        susceptibility.norm.(paramType) = 1 + (dvparam - dvuparam)./dvu ;
       % relative specificity
        specificity.raw.(paramType)  = 1 + (dvu - dvuparam)      ;
        specificity.norm.(paramType) = 1 + (dvu - dvuparam)./dvparam ;

        % interaction scores
        inter =  ((dvparam+dvu)-dvuparam) ;
        interaction.(paramType) = inter;
        interaction_norm.(paramType)  =  inter./min(dvu,dvparam);
        interaction_normU.(paramType) =  inter./dvu;
        interaction_normP.(paramType) =  inter./dvparam ;
        ii = max(inter,0) ;
        interaction_normC.(paramType) =  ii./ repmat(sum(ii),size(ii,1),1) ;

    else
        susceptibility.(paramType) = [];
        specificity.(paramType) = [];
        interaction.(paramType) = [];
        interaction_norm.(paramType) = [];
        interaction_normP.(paramType) = [];
        interaction_normU.(paramType) = [];
        interaction_normC.(paramType) = [];
    end

end

result.ind = ind ;
result.susceptibility = susceptibility;
result.specificity    = specificity;
result.interaction.raw     = interaction;
result.interaction.norm    = interaction_norm;
result.interaction.normP   = interaction_normP;
result.interaction.normU   = interaction_normU;
result.interaction.normC   = interaction_normC;

end

% ========================================================================================
% Additional functions
% ========================================================================================

% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% compute the loss of explained variance of a mutilated model
function v = mutilation_score(posterior,out,ind,varargin)

    % default is no mutilation
    mutilation = VBA_check_struct(struct,'u',0,'phi',0,'theta',0);

    % set mutilation from args
    for iArg = 1:2:numel(varargin)
        mutilation.(varargin{iArg}) = varargin{iArg+1} ;
    end

    % compute mutilated model trajectory
    [muy] = laplace_mutilated(posterior,out,mutilation);

    % compute score
    y = out.y ;
    vfull = score(y,out.suffStat.gx,out.options,ind) ;
    v = vfull - score(y,muy,out.options,ind) ;
end

% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% compute the loss of explained variance of a trajectory
function v = score(y,muy,options,ind)

    % select observation of interest
    y      =   y(ind.y,:);
    muy    = muy(ind.y,:);
    isYout = options.isYout(ind.y,:);

    % get rid of unused data
    y(isYout==1)   = [] ;
    muy(isYout==1) = [] ;

    % compute explained variance % ### TODO: multiline data
    dy = VBA_vec(y) - VBA_vec(muy);
    v =  1 - mean(dy.^2)./var(y) ;
end

% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% compute the trajectory moments for a mutilated model
function [muy] = laplace_mutilated(posterior,out,mutilation)

    mutilation = VBA_check_struct(mutilation,'u',0,'phi',0,'theta',0);

    mutilation_ratio = 0;
    % define input
    u = out.u;
    if all(mutilation.u > 0)
        u(mutilation.u,:) = mutilation_ratio*u(mutilation.u,:);
    end

    % define parameters
    opt = out.options;
    opt.priors = posterior;

    if mutilation.theta > 0
        opt.priors.muTheta(mutilation.theta)      = mutilation_ratio*opt.priors.muTheta(mutilation.theta);
        opt.priors.SigmaTheta(mutilation.theta,:) = 0;
        opt.priors.SigmaTheta(:,mutilation.theta) = 0;
    end

    if mutilation.phi > 0
        opt.priors.muPhi(mutilation.phi)      = mutilation_ratio*opt.priors.muPhi(mutilation.phi);
        opt.priors.SigmaPhi(mutilation.phi,:) = 0;
        opt.priors.SigmaPhi(:,mutilation.phi) = 0;
    end

   % compute trajectory moments
   [muy] = VBA_getLaplace(u,out.options.f_fname,out.options.g_fname,out.dim,opt,0,'skip');
   muy = reshape(muy,size(out.y));
end

author: MBB-team

repo: VBA-toolbox


text:
%demoL2iPotts_DecLap
% Reconstruction of a blurred sparse signal from incomplete measurements under
% Laplacian noise using the inverse L2-Potts functional

% create signal
groundTruth = loadSparse('sig1');
n = numel(groundTruth);

% create Gaussian convolution matrix
K = convkernel('gaussian', 51, 5);
Afull = spconvmatrix(K, n);

% select random measurements
idx = sort(randidx(numel(groundTruth), 0.5));
A = Afull(idx, :);

% create blurred and noisy signal (Gaussian noise)
fBlurry = A * groundTruth;
sigma = 0.05;
fNoisy = fBlurry + sigma* randn(size(fBlurry));

% reconstruction
gamma = 0.025;
[u, nSpikes] = minL2iSpars(fNoisy, gamma, A);

% show result
showSparse(fNoisy, u, groundTruth)

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
classdef TestLDA
    %TestLDA

    methods (Static)
        function test_1
            % data
            X = [randn(50,5)-1; randn(50,5)+1];
            labels = [ones(50,1)*1; ones(50,1)*2];
            [N,d] = size(X);

            % LDA
            lda = cv.LDA();
            lda.compute(X, labels);
            evals = lda.eigenvalues;
            evecs = lda.eigenvectors;
            ncomponents = numel(evals);

            assert(isvector(evals) && ismatrix(evecs));
            assert(isequal(size(evecs), [d ncomponents]))

            % project/reconstruct
            if ncomponents == 1
                % avoid edge case when only 1 component was chosen
                return
            end
            A = randn(40,d);
            Z = lda.project(A);
            AA = lda.reconstruct(Z);
            assert(isequal(size(Z), [size(A,1) ncomponents]))
            assert(isequal(size(AA), size(A)));
        end

        function test_2
            X = randn(50,4);
            labels = randi([1 3], [50 1]);
            lda = cv.LDA('NumComponents',2);
            lda.compute(X, labels)
            assert(numel(lda.eigenvalues) == 2);
            assert(size(lda.eigenvectors,2) == 2);
        end

        function test_3
            X = randn(50,4);
            k = 3;
            mn = mean(X);
            [V,~] = eig(cov(bsxfun(@minus,X,mn)));

            P = cv.LDA.subspaceProject(V(:,1:k), mn, X);
            validateattributes(P, {'double'}, {'size',[50 k]});

            R = cv.LDA.subspaceReconstruct(V(:,1:k), mn, P);
            validateattributes(R, {'double'}, {'size',[50 4]});
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [lossValue_images, unaryDerivative, pairwiseDerivative, predictions] = vl_logisticScoreLoss_pairwiseCompactModel( unaryPotentials, pairwisePotentials, labels, dzdy, computeMaxMarginals )
%vl_logisticScoreLoss_pairwiseCompactModel implements the logistic loss on top of the structured scores
%
% The joint score:
%   S(y, theta) = \sum_i y_i * theta^U_i + \sum_ij y_i * y_j * \theta^P_{i,j,k_ij}
% where y_i \in \{0,1\} are the variables and theta^U, \theta_P - potentials
% i indexes the nodes, ij - the edges; k_ij - the cluster index of edge ij
%
% The individual scores are based on the max-marginals of the joint score:
%   s_i = \max_{ y: y_i = 1 } S(y, theta) - \max_{ y: y_i = 0 } S(y, theta)
%
% The final loss is computed as the sum of logistics on top of the individual scores:
%   loss = \sum_{i is positive} v( s_i ) + \sum_{i is negative} v( -s_i )
%
% If the number of nodes is <= 20 than the max-marginals are computed exacly using the exhaustive search, otherwise the approximations are used, see computeMinMarginalsPairwiseBinary.m
%
% Usage:
%   [lossValue_images, unaryDerivative, pairwiseDerivative, predictions] = vl_logisticScoreLoss_pairwiseCompactModel( unaryPotentials, pairwisePotentials, labels, dzdy, computeMaxMarginals )
%
% Input:
%   unaryPotentials - unary potentials parameterized by \theta^U, double[1 x 1 x 1 x numNodes], where numNodes is the number of all nodes in the batch
%   pairwisePotentials - pairwise potentials parameterizws by \theta^P, double[1 x 1 x numClusters x numNodes]
%   labels - cell array providing information about the batch, cell[numImages x 1], each cell has
%       candidateBatchIds - indices of the nodes of the current image, double[ numNodesInImage x 1 ]
%       classGroundTruth - class labels of the nodes of the current image, double[ numNodesInImage x 1 ] with labels 0 or 1
%       instanceGroundTruth - indices of the gound-truth objects of the current image, double[ numNodesImage x 1 ] of 0 (background), 1, ..., numInstances
%       clusteredEdges - information about clusters of the edges, structure with fields:
%           bbIds - indices of the nodes connected with an edge, double[numEdges x 2] of 1,...,numNodesInImage
%           clusterId - IDs of clusters for the edges, double[numEdges x 1]
%   dzdy - if empty, only forward pass is computed, otherwise the gradient is multiplied by dzdy
%   computeMaxMarginals - flag showing whether to compute output max-marginals (default: false)
%
% Output:
%   lossValue_images - values of the loss for all images in a batch, double[numImages x 1]
%   unaryDerivative - derivatives w.r.t. the \theta^U, same size as unaryPotentials
%   pairwiseDerivative - derivatives w.r.t. the \theta^P, same size as pairwisePotentials
%   predictions - if computeMaxMarginals == false than labels for the nodes, double[numNodesImage x 1]
%           if computeMaxMarginals == true than cell[numImages x 1], each cell has
%              bestLabeling - labels for the nodes, double[numNodesImage x 1]
%              maxMarginals - max-marginals

if ~exist('computeMaxMarginals', 'var') || isempty(computeMaxMarginals)
    computeMaxMarginals = false;
end

%% preparation
numImages = length(labels);
numNodes = zeros(numImages, 1);
numEdges = zeros(numImages, 1);
numLabels = 2;
for iImage = 1 : numImages
    numNodes( iImage ) = length( labels{iImage}.candidateBatchIds );
    numEdges( iImage ) = size(labels{ iImage }.clusteredEdges.bbIds, 1);
end

if size(unaryPotentials, 1) ~= 1 ||  size(unaryPotentials, 2) ~= 1 || size(unaryPotentials, 3) ~= 1 || size(unaryPotentials, 4) ~= sum(numNodes)
    error('Unary potentials are of incorrect size');
end
numClusters = size(pairwisePotentials, 3);
if size(pairwisePotentials, 1) ~= 1 ||  size(pairwisePotentials, 2) ~= 1 || size(pairwisePotentials, 3) ~= numClusters * 1 || size(pairwisePotentials, 4) ~= sum(numEdges)
    error('Pairwise potentials are of incorrect size');
end

if ~isempty(dzdy)
    % compute derivatives
    unaryDerivative = zeros( size(unaryPotentials), 'like', unaryPotentials );
    pairwiseDerivative = zeros( size(pairwisePotentials), 'like', pairwisePotentials );
else
    unaryDerivative = [];
    pairwiseDerivative = [];
end

%% start computations
nodeStartIndex = 0;
edgeStartIndex = 0;
lossValue_images = nan(numImages, 1);
predictions = cell(numImages, 1);
for iImage = 1 : numImages
    %% extract potentials
    nodeIds = nodeStartIndex + (1 : numNodes( iImage ));
    edgeIds = edgeStartIndex + (1 : numEdges( iImage ));
    nodeStartIndex = nodeStartIndex + numNodes( iImage );
    edgeStartIndex = edgeStartIndex + numEdges( iImage );

    % extract unaries
    unaries = reshape( unaryPotentials(:,:,:,nodeIds), 1, numNodes( iImage ) );

    % extract pairwise potentials by mapping only one cluster to the edge
    edgeEnds = labels{ iImage }.clusteredEdges.bbIds;
    edgeClusters = labels{ iImage }.clusteredEdges.clusterId;

    pairwise = pairwisePotentials(:,:,:,edgeIds);

    potentialsId = edgeClusters(:) + numClusters * (0 : numEdges( iImage ) - 1)';
    pairwise = reshape( pairwise( potentialsId ), 1, numEdges( iImage ) );

    % prepare terms for energy minimizationsc
    unaryTerms = [ -double(unaries(:)), zeros(numel(unaries), 1)];
    pairwiseTerms = [ edgeEnds, -double(pairwise(:)), zeros(numel(pairwise), 3) ];

    % labels{iImage}.classGroundTruth - 0 for bkg, 1 for obj
    % labels{iImage}.instanceGroundTruth - 0 for bkg, i for detection object #i, exactly one detection for object i should be present
    % labelsBinary - 1 - head, 2 - bkg

    % compute all the max-marginals
    if numNodes( iImage ) <= 20
        [minMarginals, minMarginals_args] = computeMinMarginalsBinaryMex( unaryTerms, pairwiseTerms );
    else
        [minMarginals, ~, minMarginals_args] = computeMinMarginalsPairwiseBinary( unaryTerms, pairwiseTerms );
    end
    maxMarginals = -minMarginals;

    scores = maxMarginals(:, 1) - maxMarginals(:, 2);
    [~, bestLabeling] = max( maxMarginals, [], 2 );

    % compute the logistic loss
    labelsBinary_GT = 2 * labels{iImage}.classGroundTruth - 1; % convert GT from (0 for bkg, 1 for obj) to (1 - head, -1 - bkg)
    expGtScores = exp( - scores .* labelsBinary_GT );

    numBkg = sum( labels{iImage}.classGroundTruth == 0 );
    numObj = sum( labels{iImage}.classGroundTruth == 1 );
    weightObj = 1 / numObj;
    weightBkg = 1 / numBkg;
    weights = weightBkg * ones( numNodes( iImage ), 1 );
    weights( labels{iImage}.classGroundTruth(:) == 1 ) = weightObj;
    weights = weights(:) / sum(weights(:)) * numel(weights);

    lossValue_images(iImage) = sum( log( 1 + expGtScores(:) ) .* weights(:) );

    if ~computeMaxMarginals
        predictions{iImage} = bestLabeling;
    else
        predictions{iImage}.maxMarginals = maxMarginals;
        predictions{iImage}.bestLabeling = bestLabeling;
    end

    if ~isempty(dzdy)
        % compute derivatives
        curUnaryDerivative = zeros(1, 1, 1, numNodes( iImage ), 'like', unaryDerivative);
        curPairwiseDerivative = zeros(1, 1, numClusters, numEdges( iImage ), 'like', pairwiseDerivative);

        % derivative of the loss w.r.t. the scores
        scoreDerivatives = (expGtScores(:) ./ (1 + expGtScores(:))) .* (-labelsBinary_GT(:)) .* weights(:);

        % derivative of the loss w.r.t. the max-marginals
        mmDerivates = [ ones( numNodes( iImage ), 1 ), -ones( numNodes( iImage ), 1 ) ];
        mmDerivates = bsxfun(@times, mmDerivates, scoreDerivatives);

        % derivative of the loss w.r.t. the potentials
        for iNode = 1 : numNodes( iImage )
            for iLabel = 1 : 2
                curArg = permute( minMarginals_args(iNode, iLabel, :), [3, 1, 2]);

                % unaries
                curMask = curArg(:) == 0; % the node belongs to the head class: 0 in this notaion
                curUnaryDerivative( curMask ) = curUnaryDerivative( curMask ) + 1 * mmDerivates(iNode, iLabel);

                % pairwise
                curMask = curArg( edgeEnds(:,1) ) == 0 & curArg( edgeEnds(:,2) ) == 0; % both incident nodes belong to the head class: 0 in this notaion
                curIds = edgeClusters(curMask) + numClusters * (find(curMask) - 1);
                curPairwiseDerivative( curIds ) = curPairwiseDerivative( curIds ) + 1 * mmDerivates(iNode, iLabel);

            end
        end

        unaryDerivative(:,:,:,nodeIds) = curUnaryDerivative;
        pairwiseDerivative(:,:,:,edgeIds) = curPairwiseDerivative;
    end
end

if ~isempty(dzdy)
    unaryDerivative = unaryDerivative * dzdy;
    pairwiseDerivative = pairwiseDerivative * dzdy;
end

end

author: aosokin


text:
function [SD] = Demodulator(RxIn, PN, MF, Walsh);
%
% DEMODULATOR			This function performs demodulation of the forward
%                                       Channel packet, based on RAKE Receiver
% 						Block Diagram
%						Input Signal -> [Matched Filter] -> [Sampler] -> [RAKE Receiver] -> [Walsh] -> [DeSpreading]
%
% 						Inputs: RxIn - input signal (I/Q) analoge
%                                    PN - PN sequence (used for De-spreading)
%                                    MF - matched filter taps
%                                    Walsh - Used row of Walsh matrix for  recovering
%
% 						Outputs: SD - Soft Decisions of RAKE receiver
%

global R
N = length(RxIn)/R;

%--------------- Matched Rx Filter (Analog) --------------------
L = length(MF);
L_2 = floor(L/2);
rr = conv(flipud(conj(MF)), RxIn);
rr = rr(L_2+1: end - L_2);

%----------- Rx Symbols Sampling ------------
% R = 1.2288 Mcps
Rx = sign(real(rr(1:R:end))) + j*sign(imag(rr(1:R:end)));

%----------- RAKE Receiver ------------
Rx = reshape(Rx, 64, N/64); 				% -------- column oriented

%-------- Walsh recovering ---------
Walsh = ones(N/64, 1)*sign(Walsh'-1/2);	%--- row oriented Walsh
PN = reshape(PN, 64, N/64)'; 					%--- conjugated row oriented PN sequence
PN = PN.*Walsh;									%--- Walsh Orthogonalization

%---------- Despreading (Correlate and Sum)
% Input Rate = 1.2288 Mpbs, Output Rate = 19.2 KBps
SD= PN*Rx;
SD= real(diag(SD));  % Find Soft Decisions (on main diagonal)

author: Sable

repo: mcbench-benchmarks