text:
function b = isPointInPolygon(point, poly)
%ISPOINTINPOLYGON Test if a point is located inside a polygon
%
%   B = isPointInPolygon(POINT, POLYGON)
%   Returns true if the point is located within the given polygon.
%
%   This function is simply a wrapper for the function inpolygon, to avoid
%   decomposition of point and polygon coordinates.
%
%   Example
%     pt1 = [30 20];
%     pt2 = [30 5];
%     poly = [10 10;50 10;50 50;10 50];
%     isPointInPolygon([pt1;pt2], poly)
%     ans =
%          1
%          0
%
%     poly = [0 0; 10 0;10 10;0 10;NaN NaN;3 3;3 7;7 7;7 3];
%     pts = [5 1;5 4];
%     isPointInPolygon(pts, poly);
%     ans =
%          1
%          0
%
%
%   See also
%   points2d, polygons2d, inpolygon, isPointInTriangle

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

%   HISTORY
%   2013-04-24 add support for multiply connected polygons

% In case of a multiple polygon, decompose into a set of contours, and
% performs test for each contour
if iscell(poly) || any(isnan(poly(:)))
    % transform as a cell array of simple polygons
    polygons = splitPolygons(poly);
    N = length(polygons);
    Np = size(point, 1);

    % compute orientation of polygon, and format to have Np*N matrix
    areas = zeros(N, 1);
    for i = 1:N
        areas(i) = polygonArea(polygons{i});
    end
    ccw = areas > 0;
    ccw = repmat(ccw', Np, 1);

    % test if point inside each polygon
    in = false(size(point, 1), N);
    for i = 1:N
        poly = polygons{i};
        in(:, i) = inpolygon(point(:,1), point(:,2), poly(:,1), poly(:,2));
    end

    % count polygons containing point, weighted by polygon orientation
    b = sum(in.*(ccw==1) - in.*(ccw==0), 2) > 0;

else
    % standard test for simple polygons
    b = inpolygon(point(:,1), point(:,2), poly(:,1), poly(:,2));
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ a, det, inert ] = chidi ( a, lda, n, ipvt, job )

%*****************************************************************************80
%
%% CHIDI computes the determinant and inverse of a matrix factored by CHIFA.
%
%  Discussion:
%
%    CHIDI computes the determinant, inertia (number of positive, zero,
%    and negative eigenvalues) and inverse of a complex hermitian matrix
%    using the factors from CHIFA.
%
%    A division by zero may occur if the inverse is requested
%    and CHICO has set RCOND == 0.0 or CHIFA has set INFO /= 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex A(LDA,N); the factored matrix from CHIFA.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer IPVT(N), the pivot vector from CHIFA.
%
%    Input, integer JOB, has the decimal expansion ABC where:
%    if C /= 0, the inverse is computed,
%    if B /= 0, the determinant is computed,
%    if A /= 0, the inertia is computed.
%    For example, JOB = 111 gives all three.
%
%    Output, complex A(LDA,N); if the inverse was requested, A contains
%    the inverse matrix.  The strict lower triangle of A is never
%    referenced.
%
%    Output, real DET(2), the determinant of the original matrix.
%    Determinant = DET(1) * 10.0**DET(2) with 1.0 <= abs ( DET(1) ) < 10.0
%    or DET(1) = 0.0.
%
%    Output, integer INERT(3), the inertia of the original matrix.
%    INERT(1) = number of positive eigenvalues.
%    INERT(2) = number of negative eigenvalues.
%    INERT(3) = number of zero eigenvalues.
%
  inert = [];
  det = [];

  noinv = floor ( mod ( job,   10 )       ) == 0;
  nodet = floor ( mod ( job,  100 ) /  10 ) == 0;
  noert = floor ( mod ( job, 1000 ) / 100 ) == 0;

  if ( ~nodet || ~noert )

    if ( ~noert )
      inert(1:3) = 0;
    end

    if ( ~nodet )
      det(1) = 1.0;
      det(2) = 0.0;
    end

    t = 0.0;

    for k = 1 : n

      d = real ( a(k,k) );
%
%  Check if 1 by 1.
%
      if ( ipvt(k) <= 0 )
%
%  2 by 2 block
%  Use DET = ( D / T * C - T ) * T, T = abs ( S )
%  to avoid underflow/overflow troubles.
%  Take two passes through scaling.  Use T for flag.
%
        if ( t == 0.0 )
          t = abs ( a(k,k+1) );
          d = ( d / t ) * real ( a(k+1,k+1) ) - t;
        else
          d = t;
          t = 0.0;
        end

      end

      if ( ~noert )
        if ( 0.0 < d )
          inert(1) = inert(1) + 1;
        elseif ( d < 0.0 )
          inert(2) = inert(2) + 1;
        elseif ( d == 0.0 )
          inert(3) = inert(3) + 1;
        end
      end

      if ( ~nodet )

        det(1) = det(1) * d;

        if ( det(1) ~= 0.0 )

          while ( abs ( det(1) ) < 1.0 )
            det(1) = det(1) * 10.0;
            det(2) = det(2) - 1.0;
          end

          while ( 10.0 <= abs ( det(1) ) )
            det(1) = det(1) / 10.0;
            det(2) = det(2) + 1.0;
          end

        end

      end

    end

  end
%
%  Compute inverse(A).
%
  if ( ~noinv )

    k = 1;

    while ( k <= n )

      km1 = k - 1;

      if ( 0 <= ipvt(k) )
%
%  1 by 1
%
        a(k,k) = 1.0 / real ( a(k,k) );

        if ( 1 <= km1 )

          work(1:km1) = a(1:km1,k);

          for j = 1 : km1
            a(j,k) = work(1:j) * conj ( a(1:j,j) );
            a(1:j-1,k) = a(1:j-1,k) + work(j) * a(1:j-1,j);
          end

          a(k,k) = a(k,k) + real ( conj ( work(1:km1) ) * a(1:km1,k) );

        end

        kstep = 1;

      else
%
%  2 by 2
%
        t = abs ( a(k,k+1) );
        ak = real ( a(k,k) ) / t;
        akp1 = real ( a(k+1,k+1) ) / t;
        akkp1 = a(k,k+1) / t;
        d = t * ( ak * akp1 - 1.0 );
        a(k,k) = akp1 / d;
        a(k+1,k+1) = ak / d;
        a(k,k+1) = -akkp1 / d;

        if ( 1 <= km1 )

          work(1:km1) = a(1:km1,k+1);

          for j = 1 : km1
            a(j,k+1) = work(1:j) * conj ( a(1:j,j) );
            a(1:j-1,k+1) = a(1:j-1,k+1) + work(j) * a(1:j-1,j);
          end

          a(k+1,k+1) = a(k+1,k+1) + ...
            real ( conj ( work(1:km1) ) * a(1:km1,k+1) );

          a(k,k+1) = a(k,k+1) ...
            + transpose ( conj ( a(1:km1,k) ) ) * a(1:km1,k+1);

          work(1:km1) = a(1:km1,k);

          for j = 1 : km1
            a(j,k) = work(1:j) * conj ( a(1:j,j) );
            a(1:j-1,k) = a(1:j-1,k) + work(j) * a(1:j-1,j);
          end

          a(k,k) = a(k,k) + real ( conj ( work(1:km1) ) * a(1:km1,k) );

        end

        kstep = 2;

      end
%
%  Swap
%
      ks = abs ( ipvt(k) );

      if ( ks ~= k )

        temp       = a(1:ks,ks);
        a(1:ks,ks) = a(1:ks,k);
        a(1:ks,k)  = temp;

        for j = k : -1 : ks
          temp    = conj ( a(j,k) );
          a(j,k)  = conj ( a(ks,j) );
          a(ks,j) = temp;
        end

        if ( kstep ~= 1 )
          temp      = a(ks,k+1);
          a(ks,k+1) = a(k,k+1);
          a(k,k+1)  = temp;
        end

      end

      k = k + kstep;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear all
close all
path(path,'..\..\..\FUZZCLUST')
%data set
data.X = nDexample(5,400,2,1)

%parameters
param.c=3;
param.m=2;
param.e=1e-6;
param.ro=ones(1,param.c);
param.val=1;
%normalization
data=clust_normalize(data,'range');

result = GKclust(data,param);
result = validity(result,data,param);

plot(data.X(:,1),data.X(:,2),'b.',result.cluster.v(:,1),result.cluster.v(:,2),'ro');
hold on
%draw contour-map
new.X=data.X;
eval=clusteval(new,result,param);
result.validity

author: Sable

repo: mcbench-benchmarks


text:
function mesh = genMesh3DType2(domain, nx, ny, nz)

%% Usage: mesh structure of a uniform cubic partition
%
% INPUTS:
% domain --- cubic domain = [xmin, xmax, ymin, ymax, zmin, zmax].
% nx --- the number of uniform partition in x direction.
% ny --- the number of uniform partition in y direction.
% nz --- the number of uniform partition in z direction.
%
% option.meshinfo --- basic (default): generate only p and t.
%                     all: enriched mesh information
%
% OUTPUTS:
% mesh --- a struct data contains mesh information.
%
% Last Modified: 08/07/2020 by Xu Zhang
%
%                A8-------------------A7        The Cube is divided into
%                /|                   /|        six congruent tetrahedrons
%               / |                  / |
%              /  |                 /  |
%             /   |                /   |        (1) A1-A2-A3-A7
%            /    |               /    |        (2) A1-A6-A2-A7
%          A5-----+-------------A6     |        (3) A1-A5-A6-A7
%           |     |             |      |        (4) A1-A8-A5-A7
%           |     |             |      |        (5) A1-A4-A8-A7
%           |     |             |      |        (6) A1-A3-A4-A7
%           |     A4------------+------A3
%           |     /             |      /
%           |    /              |     /
%           |   /               |    /
%           |  /                |   /
%           | /                 |  /
%           |/                  | /
%          A1-------------------A2
%
%% 1. Generate basic mesh info: p t
[p,T] = genMesh3DRectPT(domain, nx, ny, nz);
c1 = zeros(5,4);
c1(1,:) = [1,2,4,5];
c1(2,:) = [2,7,5,6];
c1(3,:) = [3,2,4,7];
c1(4,:) = [5,7,4,8];
c1(5,:) = [2,4,5,7];
c2 = zeros(5,4);
c2(1,:) = [1,6,8,5];
c2(2,:) = [1,3,6,2];
c2(3,:) = [8,6,3,7];
c2(4,:) = [1,8,3,4];
c2(5,:) = [1,3,8,6];

CBxy = zeros(nx,ny);
CBxy(1:2:nx,2:2:ny)=1;
CBxy(2:2:nx,1:2:ny)=1;
CB = repmat(CBxy,1,1,nz);
CB(:,:,1:2:nz) = repmat(1-CBxy,1,1,length(1:2:nz));
CB = reshape(CB,[],1);
tid1 = find(CB==1); tid2 = find(CB==0);

t = T(:,c1(1,:));
t(tid2,:) = T(tid2,c2(1,:));
for i = 2:5
    tmp = T(:,c1(i,:));
    tmp(tid2,:) = T(tid2,c2(i,:));
    t = [t;tmp];
end

mesh = struct('p',p,'t',t,'T',T);

repo: ifem


text: %% ODFAnalysis_index.m

author: mtex-toolbox

repo: mtex

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


text:
function hmm = updateOmega(hmm,Gamma,residuals,XX,XXGXX,XW,Tfactor,rangeK)

is_gaussian = hmm.train.order == 0; % true if Gaussian observation model is being used
K = hmm.K; ndim = hmm.train.ndim;
if nargin < 8 || isempty(rangeK), rangeK = 1:K; end
if nargin < 7, Tfactor = 1; end
setstateoptions
T = size(residuals,1);
if isfield(hmm.train,'B'), Q = size(hmm.train.B,2);
else Q = ndim; end

if (strcmp(hmm.train.covtype,'uniquediag') || strcmp(hmm.train.covtype,'shareddiag')) ...
        && hmm.train.uniqueAR
    % all are AR and there's a single covariance matrix
    hmm.Omega.Gam_rate = hmm.prior.Omega.Gam_rate;
    for k = rangeK
        if ~isempty(XW)
            XWk = XW(:,:,k);
        else
            XWk = zeros(size(residuals));
        end
        e = (residuals - XWk).^2;
        swx2 = zeros(T,ndim);
        for n=1:ndim
            ind = n:ndim:size(XX,2);
            tmp = XX(:,ind) * hmm.state(k).W.S_W;
            swx2(:,n) = sum(tmp .* XX(:,ind),2);
        end
        hmm.Omega.Gam_rate = hmm.Omega.Gam_rate + ...
            0.5 * Tfactor * sum( bsxfun(@times, e + swx2, Gamma(:,k)) );
    end
    hmm.Omega.Gam_shape = hmm.prior.Omega.Gam_shape + 0.5 * Tfactor * T;

elseif strcmp(hmm.train.covtype,'uniquediag') || strcmp(hmm.train.covtype,'shareddiag')
    hmm.Omega.Gam_rate(regressed) = hmm.prior.Omega.Gam_rate(regressed);
    for k = rangeK
        if ~isempty(XW)
            XWk = XW(:,:,k);
        else
            XWk = zeros(size(residuals));
        end
        e = (residuals(:,regressed) - XWk(:,regressed)).^2;
        swx2 = zeros(T,ndim);
        if ~isempty(hmm.state(k).W.Mu_W)
            for n = 1:ndim
                if ~regressed(n), continue; end
                if ndim==1
                    tmp = XX(:,Sind(:,n)) * hmm.state(k).W.S_W;
                else
                    tmp = XX(:,Sind(:,n)) * permute(hmm.state(k).W.S_W(n,Sind(:,n),Sind(:,n)),[2 3 1]);
                end
                swx2(:,n) = sum(tmp .* XX(:,Sind(:,n)),2);
            end
        end
        hmm.Omega.Gam_rate(regressed) = hmm.Omega.Gam_rate(regressed) + ...
            0.5 * Tfactor * sum( bsxfun(@times, e + swx2(:,regressed), Gamma(:,k)));
    end
    hmm.Omega.Gam_shape = hmm.prior.Omega.Gam_shape + 0.5 * Tfactor * T;

elseif strcmp(hmm.train.covtype,'uniquefull') || strcmp(hmm.train.covtype,'sharedfull')
    hmm.Omega.Gam_rate(regressed,regressed) = hmm.prior.Omega.Gam_rate(regressed,regressed);
    for k = rangeK
        if ~isempty(XW)
            XWk = XW(:,:,k);
        else
            XWk = zeros(size(residuals));
        end
        e = (residuals(:,regressed) - XWk(:,regressed));
        e = (bsxfun(@times,e,Gamma(:,k)))' * e;
        if all(S(:)==1)
            swx2 =  zeros(ndim,ndim);
            if ~isempty(hmm.state(k).W.Mu_W)
                for n1=find(regressed)
                    for n2=find(regressed)
                        if n2<n1, continue, end
                        if hmm.train.pcapred>0
                            index1 = (0:hmm.train.pcapred+(~hmm.train.zeromean)-1) * ndim + n1;
                            index2 = (0:hmm.train.pcapred+(~hmm.train.zeromean)-1) * ndim + n2;
                        else
                            index1 = (0:length(orders)*Q+(~hmm.train.zeromean)-1) * ndim + n1;
                            index2 = (0:length(orders)*Q+(~hmm.train.zeromean)-1) * ndim + n2;
                            index1 = index1(Sind(:,n1)); index2 = index2(Sind(:,n2));
                        end
                        swx2(n1,n2) = sum(sum(hmm.state(k).W.S_W(index1,index2) .* XXGXX{k}(Sind(:,n1),Sind(:,n2))));
                        swx2(n2,n1) = swx2(n1,n2);
                    end
                end
            end
            hmm.Omega.Gam_rate(regressed,regressed) = hmm.Omega.Gam_rate(regressed,regressed) ...
                + Tfactor * (e + swx2(regressed,regressed));
        else % multivariate regression model
            index_iv = sum(S,2)>0; % the independent variables
            xdim = sum(index_iv);
            ydim = sum(regressed);
            S_W = hmm.state(k).W.S_W(logical(S(:)),logical(S(:)));
            L = chol(S_W)';
            XGX = (bsxfun(@times,residuals(:,index_iv),Gamma(:,k)))' * residuals(:,index_iv);
            tracesum = zeros(ydim,ydim);
            for iL = 1:xdim*ydim
                vecinvL = reshape(L(:,iL),xdim,ydim);
                tracesum = tracesum + vecinvL'*XGX*vecinvL;
            end
            hmm.Omega.Gam_rate(regressed,regressed) = hmm.Omega.Gam_rate(regressed,regressed) + tracesum + e;
        end
    end
    hmm.Omega.Gam_rate(regressed,regressed) = (hmm.Omega.Gam_rate(regressed,regressed) + ...
        hmm.Omega.Gam_rate(regressed,regressed)') / 2;
    hmm.Omega.Gam_irate(regressed,regressed) = inv(hmm.Omega.Gam_rate(regressed,regressed));
    hmm.Omega.Gam_shape = hmm.prior.Omega.Gam_shape + Tfactor * T;

else % state dependent

    Gammasum = sum(Gamma);
    for k = rangeK
        if ~hmm.train.active, continue; end
        if ~isempty(XW)
            XWk = XW(:,:,k);
        else
            XWk = zeros(size(residuals));
        end
        if train.uniqueAR
            e = (residuals - XWk).^2;
            swx2 = zeros(T,ndim);
            for n = 1:ndim
                ind = n:ndim:size(XX,2);
                swx2(:,n) = sum(XX(:,ind) * hmm.state(k).W.S_W .* XX(:,ind),2);
            end
            hmm.state(k).Omega.Gam_rate = hmm.state(k).prior.Omega.Gam_rate + ...
                0.5 * Tfactor * sum( bsxfun(@times, e + swx2, Gamma(:,k)) );
            hmm.state(k).Omega.Gam_shape = hmm.state(k).prior.Omega.Gam_shape ...
                + 0.5 * Tfactor * Gammasum(k);

        elseif strcmp(train.covtype,'diag')
            e = (residuals(:,regressed) - XWk(:,regressed)).^2;
            swx2 = zeros(T,ndim);
            if ~isempty(hmm.state(k).W.Mu_W)
                for n=1:ndim
                    if ~regressed(n), continue; end
                    if ndim==1
                        swx2(:,n) = sum(XX(:,Sind(:,n)) * hmm.state(k).W.S_W(Sind(:,n),Sind(:,n)) ...
                            .* XX(:,Sind(:,n)),2);
                    else
                        swx2(:,n) = sum(XX(:,Sind(:,n)) * permute(hmm.state(k).W.S_W(n,Sind(:,n),Sind(:,n)),[2 3 1]) ...
                            .* XX(:,Sind(:,n)),2);
                    end
                end
            end
            hmm.state(k).Omega.Gam_rate(regressed) = hmm.state(k).prior.Omega.Gam_rate(regressed) + ...
                0.5 * Tfactor * sum( bsxfun(@times, e + swx2(:,regressed), Gamma(:,k)) );
            hmm.state(k).Omega.Gam_shape = hmm.state(k).prior.Omega.Gam_shape + 0.5 * Tfactor * Gammasum(k);

        else % full
            if is_gaussian
                if hmm.train.zeromean % If zeromean == 1 for Gaussian model, then XWk=0 and we don't need to subtract at all
                    e = residuals(:,regressed);
                else % If zeromean == 0 for Gaussian model, then XWk has all the same rows, and bsxfun() is a fair bit faster
                    e = bsxfun(@minus,residuals(:,regressed),XWk(1,regressed));
                end
            else
                e = (residuals(:,regressed) - XWk(:,regressed));
            end
            e = (bsxfun(@times,e,Gamma(:,k)))' * e;
            if all(S(:)==1)
                swx2 =  zeros(ndim,ndim);
                if ~isempty(hmm.state(k).W.Mu_W)
                    orders = formorders(train.order,train.orderoffset,train.timelag,train.exptimelag);
                    if isempty(orders)
                        swx2 = hmm.state(k).W.S_W * XXGXX{k};
                    else
                        for n1=find(regressed)
                            for n2=find(regressed)
                                if n2<n1, continue, end
                                if hmm.train.pcapred>0
                                    index1 = (0:hmm.train.pcapred+(~train.zeromean)-1) * ndim + n1;
                                    index2 = (0:hmm.train.pcapred+(~train.zeromean)-1) * ndim + n2;
                                else
                                    index1 = (0:length(orders)*Q+(~train.zeromean)-1) * ndim + n1;
                                    index2 = (0:length(orders)*Q+(~train.zeromean)-1) * ndim + n2;
                                    index1 = index1(Sind(:,n1)); index2 = index2(Sind(:,n2));
                                end
                                swx2(n1,n2) = sum(sum(hmm.state(k).W.S_W(index1,index2) .* XXGXX{k}(Sind(:,n1),Sind(:,n2))));
                                swx2(n2,n1) = swx2(n1,n2);
                            end
                        end
                    end
                end
                hmm.state(k).Omega.Gam_rate(regressed,regressed) = hmm.state(k).prior.Omega.Gam_rate(regressed,regressed) + ...
                    Tfactor * (e + swx2(regressed,regressed));
            else % multivariate regression model
                index_iv = sum(S,2)>0; % the independent variables
                xdim = sum(index_iv);
                ydim = sum(regressed);
                S_W = hmm.state(k).W.S_W(logical(S(:)),logical(S(:)));
                L = chol(S_W)';
                XGX = (bsxfun(@times,residuals(:,index_iv),Gamma(:,k)))' * residuals(:,index_iv);
                tracesum = zeros(ydim,ydim);
                for iL = 1:xdim*ydim
                    vecinvL = reshape(L(:,iL),xdim,ydim);
                    tracesum = tracesum + vecinvL'*XGX*vecinvL;
                end

                hmm.state(k).Omega.Gam_rate(regressed,regressed) = tracesum + e + ...
                    hmm.state(k).prior.Omega.Gam_rate(regressed,regressed);
            end

            hmm.state(k).Omega.Gam_shape = hmm.state(k).prior.Omega.Gam_shape + Tfactor * Gammasum(k);
            hmm.state(k).Omega.Gam_irate(regressed,regressed) = inv(hmm.state(k).Omega.Gam_rate(regressed,regressed));

            %if train.FC
            %    C0 = hmm.state(k).Omega.Gam_rate(regressed,regressed) / hmm.state(k).Omega.Gam_shape;
            %    C = hmm.train.A' * corrcov(hmm.train.A * C0 * hmm.train.A',1) * hmm.train.A;
            %    %iC = hmm.train.A' * pinv(corrcov(hmm.train.A * C0 * hmm.train.A'),1e-10) * hmm.train.A;
            %    iC = hmm.train.A' * ( (corrcov(hmm.train.A * C0 * hmm.train.A',1)+1e-4*eye(size(hmm.train.A,1))) ...
            %           \ hmm.train.A);
            %    hmm.state(k).Omega.Gam_rate(regressed,regressed) = C * hmm.state(k).Omega.Gam_shape;
            %    hmm.state(k).Omega.Gam_irate(regressed,regressed) = iC / hmm.state(k).Omega.Gam_shape;
            %end

            % ensuring symmetry
            hmm.state(k).Omega.Gam_rate(regressed,regressed) = (hmm.state(k).Omega.Gam_rate(regressed,regressed) + ...
                hmm.state(k).Omega.Gam_rate(regressed,regressed)') / 2;
            hmm.state(k).Omega.Gam_irate(regressed,regressed) = (hmm.state(k).Omega.Gam_irate(regressed,regressed) + ...
                hmm.state(k).Omega.Gam_irate(regressed,regressed)') / 2;

        end
    end

end

end

author: OHBA-analysis

repo: HMM-MAR


text:
addpath('/home/liisa/080812/');
addpath('/home/rheintz/matlab/');

global para;

s=19;

% beadsAt=[180 504; 207 489; 253 236; 399 276; 262 54; 364 49; 393 168; 410 156; 300 139; 190 138; 376 282; 182 459; 333 449; 367 467; 415 453; 204 550; 175 579];

myres=para.res_resIm
fprintf('Please select bead coordinates by left clicking. Finish with right click\n');
beadsAt=dipgetcoords(100);
fprintf('%d beads selected\n',size(beadsAt,1));
beadsAt(end,:)=[];

[sumBdsHighRes,mybdsHRes,sumBdsSum,mybdsSum]=MeanFromCoord(myres,beadsAt,s,para.res_sumIm)

%myresSum=para.res_sumIm;
%[sumBdsSum,mybdsSum]=MeanFromCoord(myresSum,beadsAt,s)

%im = readtimeseries('beads1_000');
%myres=squeeze(sum(im(:,:,3:5)));
%[sumBdsRaw,mybdsRaw]=MeanFromCoord(myres,beadsAt,s)

% [params,res,fitted]=FitDataNDFast([10 10 0; 1000 s/2 s/2],sumBdsSum,300,'idiv')

[fittedSum,paramsSum,idiv,myfunct] = FitDataND('c(1)*exp(-((x{1}-c(2)).^2+(x{2}-c(3)).^2)/c(4))+c(5)',[max(sumBdsSum) 0 0 20 min(sumBdsSum)],sumBdsSum,1000);
sqrt(paramsSum(4))*para.res_resampledPixelsize(1) * 2*sqrt(log(2))  % to get FWHM
sumBdsSum-fittedSum

[fittedHRes,paramsHRes,idiv,myfunct] = FitDataND('c(1)*exp(-((x{1}-c(2)).^2+(x{2}-c(3)).^2)/c(4))+c(5)',[max(sumBdsHighRes) 0 0 20 min(sumBdsHighRes)],sumBdsHighRes,1000);
sqrt(paramsHRes(4))*para.res_resampledPixelsize(1)* 2*sqrt(log(2))  % to get FWHM
sumBdsHighRes-fittedHRes

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [V_RF,V_U] = SIPEVD_method(H, Vn)

global Nrf;
[Nt, Ns, Nk] = size(H);
X = zeros(Nt,Nt,Nk);
V_U = zeros(Nrf,Ns,Nk);

for k = 1:Nk
    %X(:,:,k) = (1/Vn(k)/Nt * H(:,:,k)*H(:,:,k)' + eye(Nt))^(-1);
   % X(:,:,k) = eye(Nt) - 1/Vn(k)/Nt * H(:,:,k)*H(:,:,k)';
   X(:,:,k) = 1/Vn(k)/Nt * H(:,:,k) * inv(eye(Ns)+1/Vn(k)/Nt*H(:,:,k)'*H(:,:,k))*H(:,:,k)';
end

X = sum(X,3);
[V,D] = eig(X);
% get the Nrf largest eigenvector
[~,IX] = sort(diag(D),'descend');
V = V(:,IX);
V_RF =  exp(1i*angle(V(:,1:Nrf)));

for k = 1:Nk
    V_U(:,:,k) = inv(V_RF'*H(:,:,k) * H(:,:,k)'* V_RF+  Vn(k) *(V_RF)'*V_RF)*V_RF'*H(:,:,k);
end

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
function D = dmatrix(P)

% Distance Matrix (Eucledean distance)
%
% A simple function to calculate the distances of all pairs of points in a
% given group, as rows of a matrix P. The advantage of this method is that it
% exhibits a small computational cost.
%
% Created by Dr. Zacharias Voulgaris, 2009

N = size(P,1);
d(N,N) = 0;

for i = 1:(N-1)
    X = ones((N-i),1)*P(i,:);
    d((i+1):N,i) = sqrt( sum((X - P((i+1):N,:)).^2,2) );
end

D =d + d';


text:
function symb_pvec = sdisplay2(p,names)

% TODO: sdpvar a b
%       x = (a+b)^0.3        -- writes "mpower_internal"
%
% TODO: sdpvar h x k
%       sdisplay2(h*x - k)   -- misses the minus in front of "k"

if nargin < 2
    LinearVariables = 1:yalmip('nvars');
    x = recover(LinearVariables);
    names = cell(length(x),1);

    W = evalin('caller','whos');
    for i = 1:size(W,1)
        if strcmp(W(i).class,'sdpvar') | strcmp(W(i).class,'ncvar')
            % Get the SDPVAR variable
            thevars = evalin('caller',W(i).name);

            % Distinguish 4 cases
            % 1: Sclalar varible x
            % 2: Vector variable x(i)
            % 3: Matrix variable x(i,j)
            % 4: Variable not really defined
            if is(thevars,'scalar') & is(thevars,'linear') & length(getvariables(thevars))==1 & isequal(getbase(thevars),[0 1])
                index_in_p = find(ismember(LinearVariables,getvariables(thevars)));

                if ~isempty(index_in_p)
                    already = ~isempty(names{index_in_p});
                    if already
                        already = ~strfind(names{index_in_p},'internal');
                        if isempty(already)
                            already = 0;
                        end
                    end
                else
                    already = 0;
                end

                if ~isempty(index_in_p) & ~already
                    % Case 1
                    names{index_in_p}=W(i).name;
                end
            elseif is(thevars,'lpcone')

                if size(thevars,1)==size(thevars,2)
                    % Case 2
                    vars = getvariables(thevars);
                    indicies = find(ismember(vars,LinearVariables));
                    for ii = indicies
                        index_in_p = find(ismember(LinearVariables,vars(ii)));

                        if ~isempty(index_in_p)
                            already = ~isempty(names{index_in_p});
                            if already
                                already = ~strfind(names{index_in_p},'internal');
                                if isempty(already)
                                    already = 0;
                                end
                            end
                        else
                            already = 0;
                        end

                        if ~isempty(index_in_p) & ~already
                            B = reshape(getbasematrix(thevars,vars(ii)),size(thevars,1),size(thevars,2));
                            [ix,jx,kx] = find(B);
                            ix=ix(1);
                            jx=jx(1);
                            names{index_in_p}=[W(i).name '(' num2str(ix) ',' num2str(jx) ')'];
                        end
                    end

                else
                    % Case 3
                    vars = getvariables(thevars);
                    indicies = find(ismember(vars,LinearVariables));
                    for ii = indicies
                        index_in_p = find(ismember(LinearVariables,vars(ii)));

                        if ~isempty(index_in_p)
                            already = ~isempty(names{index_in_p});
                            if already
                                already = ~strfind(names{index_in_p},'internal');
                                if isempty(already)
                                    already = 0;
                                end
                            end
                        else
                            already = 0;
                        end

                        if ~isempty(index_in_p) & ~already
                            names{index_in_p}=[W(i).name '(' num2str(ii) ')'];
                        end
                    end
                end

            elseif is(thevars,'sdpcone')
                % Case 3
                vars = getvariables(thevars);
                indicies = find(ismember(vars,LinearVariables));
                for ii = indicies
                    index_in_p = find(ismember(LinearVariables,vars(ii)));
                    if ~isempty(index_in_p)
                        already = ~isempty(names{index_in_p});
                        if already
                            already = ~strfind(names{index_in_p},'internal');
                        end
                    else
                        already = 0;
                    end

                    if ~isempty(index_in_p) & ~already
                        B = reshape(getbasematrix(thevars,vars(ii)),size(thevars,1),size(thevars,2));
                        [ix,jx,kx] = find(B);
                        ix=ix(1);
                        jx=jx(1);
                        names{index_in_p}=[W(i).name '(' num2str(ix) ',' num2str(jx) ')'];
                    end
                end

            else
                % Case 4
                vars = getvariables(thevars);
                indicies = find(ismember(vars,LinearVariables));

                for i = indicies
                    index_in_p = find(ismember(LinearVariables,vars(i)));
                    if ~isempty(index_in_p) & isempty(names{index_in_p})
                        names{index_in_p}=['internal(' num2str(vars(i)) ')'];
                    end
                end

            end
        end
    end
end

[mt,vt] = yalmip('monomtable');
ev = yalmip('extvariables');

for i = 1:size(p, 1)
    for j = 1:size(p, 2)
        symb_pvec{i, j} = symbolicdisplay(p(i, j), names, vt, ev, mt);
    end
end

%------------------------------------------------------------------------
function expression = symbolicdisplay(p,names,vt,ev,mt)

sp = size(p);
if any(sp > 1)
    out = '[';
else
    out = '';
end
p_orig = p;
for i1 = 1:sp(1)
    for i2 = 1:sp(2)
        p = p_orig(i1, i2);
        basis = getbase(p);
        if basis(1)~=0
            expression = [num2str(basis(1)) '+'];
        else
            expression = [''];
        end

        [dummy, variables, coeffs] = find(basis(2:end));
        variables = getvariables(p);
        for i = 1:length(coeffs)
            if coeffs(i)==1
                expression = [expression symbolicmonomial(variables(i), ...
                    names,vt,ev,mt) '+'];
            else
                expression = [expression num2str(coeffs(i)) '*' ...
                    symbolicmonomial(variables(i),names,vt,ev,mt) '+'];
            end
        end
        expression(end) = [];
        out = [out expression ','];
    end
    out(end) = ';';
end
out(end) = [];
if any(sp > 1)
    out = [out ']'];
end
expression = out;

%------------------------------------------------------------------------
function s = symbolicmonomial(variable,names,vt,ev,mt)

terms = find(mt(variable,:));
if ismember(variable,ev)
    q = yalmip('extstruct',variable);
    s = [q.fcn '(' symbolicdisplay(q.arg{1},names,vt,ev,mt)];
    for i = 2:length(q.arg)-1
        s = [s ',' symbolicdisplay(q.arg{i}, names, vt, ev, mt)];
    end
    s = [s ')'];

elseif ~vt(variable)
    % Linear expression
    s = names{variable};
else
    % Fancy display of a monomial
    s = [''];
    for i = 1:length(terms)
        if mt(variable,terms(i)) == 1
            exponent = '';
        else
            exponent = ['^' num2str(mt(variable,terms(i)))];
        end
        s = [s symbolicmonomial(terms(i),names,vt,ev,mt) exponent '*'];
    end
    s(end)=[];
end
% s = strrep(s,'^1+','+');
% s = strrep(s,'^1*','*');

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Testing Data points L = 495

% TO RUN
% ---------
% mydata = mgts(3000);
% n = 0;
% RMSE = ttsf(mydata)

function RMSE = ttsf(mydata)               % Test Time Series Forecasting
L=400;                                                           % No. of Test Points
for n = 0:L
    fcs(n+505) = st1(mydata, n);
end

% Plot My Data
t=1001:2000;
subplot(2, 1, 1);
plot(t, mydata)
h = legend('My Data', 1);

% Plot Forecasted & Test Data Simultaneously
t = 505:505+L;
tsd = mydata(505:505+L);
fcd = fcs(505:505+L);
subplot(2, 1, 2);
plot(t, fcd, '-k', t, tsd, ':b')
h = legend('Forecast','Test Data',2);

ErrS = 0;
for t=1:L
    ErrS = ErrS + (tsd(t) - fcd(t))^2;
end

RMSE = sqrt(ErrS/L);


text:
function [ yval, ypval, yppval ] = spline_cubic_val ( n, t, y, ypp, tval )

%*****************************************************************************80
%
%% SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point.
%
%  Discussion:
%
%    SPLINE_CUBIC_SET must have already been called to define the
%    values of YPP.
%
%    For any point T in the interval T(IVAL), T(IVAL+1), the form of
%    the spline is
%
%      SPL(T) = A
%             + B * ( T - T(IVAL) )
%             + C * ( T - T(IVAL) )**2
%             + D * ( T - T(IVAL) )**3
%
%    Here:
%      A = Y(IVAL)
%      B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
%        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
%      C = YPP(IVAL) / 2
%      D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl de Boor,
%    A Practical Guide to Splines,
%    Springer Verlag, 1978.
%
%  Parameters:
%
%    Input, integer N, the number of data values.
%
%    Input, real T(N), the knot values.
%
%    Input, real Y(N), the data values at the knots.
%
%    Input, real YPP(N), the second derivatives of the spline at the knots.
%
%    Input, real TVAL, a point, typically between T(1) and T(N), at
%    which the spline is to be evalulated.  If TVAL lies outside
%    this range, extrapolation is used.
%
%    Output, real YVAL, YPVAL, YPPVAL, the value of the spline, and
%    its first two derivatives at TVAL.
%

%
%  Determine the interval [T(LEFT), T(RIGHT)] that contains TVAL.
%  Values below T(1) or above T(N) use extrapolation.
%
  [ left, right ] = r8vec_bracket ( n, t, tval );
%
%  Evaluate the polynomial.
%
  dt = tval - t(left);
  h = t(right) - t(left);

  yval = y(left) ...
       + dt * ( ( y(right) - y(left) ) / h ...
              - ( ypp(right) / 6.0 + ypp(left) / 3.0 ) * h ...
       + dt * ( 0.5 * ypp(left) ...
       + dt * ( ( ypp(right) - ypp(left) ) / ( 6.0 * h ) ) ) );

  ypval = ( y(right) - y(left) ) / h ...
       - ( ypp(right) / 6.0 + ypp(left) / 3.0 ) * h ...
       + dt * ( ypp(left) ...
       + dt * ( 0.5 * ( ypp(right) - ypp(left) ) / h ) );

  yppval = ypp(left) + dt * ( ypp(right) - ypp(left) ) / h;

  return
end


text:
function [ t_start, y_start ] = p23_start ( neqn )

%*****************************************************************************80
%
%% P23_START returns the starting point for problem p23.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_START, Y_START(NEQN,1), the initial data.
%
  y_start = zeros ( neqn, 1 );

  t_start = 0.0;
  y_start = [ 0.0; 0.0 ];

  return
end


text:
clear all;
close all;
clc;
rng default;

% Ambient space dimension
M = 50;
% Number of subspaces
K = 10;
% common dimension for each subspace
D = 20;

% Construct bases for random subspaces
bases = spx.data.synthetic.subspaces.random_subspaces(M, K, D);
% samples angles between subspaces
angles_matrix = spx.la.spaces.smallest_angles_deg(bases)
angles = spx.matrix.off_diag_upper_tri_elements(angles_matrix)';
min(angles)

% Number of points on each subspace
Sk = 4 * D;

cluster_sizes = Sk * ones(1, K);
% total number of points
S = sum(cluster_sizes);

points_result = spx.data.synthetic.subspaces.uniform_points_on_subspaces(bases, cluster_sizes);
X0 = points_result.X;

% noise level
sigma = 0.5;
% Generate noise
Noise = sigma * spx.data.synthetic.uniform(M, S);
% Add noise to signal
X = X0 + Noise;
% Normalize noisy signals.
X = spx.norm.normalize_l2(X);
% labels assigned to the data points
true_labels = spx.cluster.labels_from_cluster_sizes(cluster_sizes);
cvx_solver sdpt3
cvx_quiet(true);

% storage for coefficients
Z = zeros(S, S);
start_time = tic;
fprintf('Processing %d signals\n', S);
for s=1:S
    fprintf('.');
    if (mod(s, 50) == 0)
        fprintf('\n');
    end
    x = X(:, s);
    cvx_begin
    % storage for  l1 solver
    variable z(S, 1);
    minimize norm(z, 1)
    subject to
    x == X*z;
    z(s) == 0;
    cvx_end
    Z(:, s)  = z;
end
elapsed_time  = toc(start_time);
fprintf('\n Time spent: %.2f seconds\n', elapsed_time);
W = abs(Z) + abs(Z).';
clustering_result = spx.cluster.spectral.simple.normalized_symmetric(W);
cluster_labels = clustering_result.labels;
% Time to compare the clustering
comparsion_result = spx.cluster.clustering_error_hungarian_mapping(cluster_labels, true_labels, K);
clustering_error_perc = comparsion_result.error_perc;
clustering_acc_perc = 100 - comparsion_result.error_perc;

spr_stats = spx.cluster.subspace.subspace_preservation_stats(Z, cluster_sizes);
spr_error = spr_stats.spr_error;
spr_flag = spr_stats.spr_flag;
spr_perc = spr_stats.spr_perc;
elapsed_time  = toc(start_time);
fprintf('\nclustering error: %0.2f %%, clustering accuracy: %0.2f %% \n'...
    , clustering_error_perc, clustering_acc_perc);
fprintf('mean spr error: %0.2f, preserving : %0.2f %%\n', spr_stats.spr_error, spr_stats.spr_perc);
fprintf('elapsed time: %0.2f sec', elapsed_time);
fprintf('\n\n');

author: indigits

repo: sparse-plex


text:
function res = type(this, value)
% Method for and getting/setting EEG file type
% FORMAT res = type(this, value)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: type.m 1267 2008-03-28 12:12:14Z vladimir $

if nargin == 1
    res = this.type;
else
    switch value
        case 'continuous'
            if ntrials(this)>1
                error('Continuous file can only have one trial');
            end
        case 'single'
        case 'evoked' % Add additional checks here
        case 'grandmean'
        otherwise
            error('Unrecognized type');
    end

    this.type = value;
    res = this;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [passed_clusters] = plot(this,p_skip,cl,cluster_draw_thresh,only_cluster)

% plotDatabase(this,p_skip,cl,cluster_draw_thresh,only_cluster)
%
% this is pathway database
% 'cl' is vector with clustered path labels

if ~exist('p_skip','var')
    p_skip=1;
end
if ~exist('only_cluster','var')
    only_cluster = 0;
end
if ~exist('cluster_draw_thresh','var')
    cluster_draw_thresh = 0;
end

if exist('cl','var')
    cluster_sizes = hist(cl,[1:max(cl)]);
    npaths = sum(cluster_sizes);
    if (cluster_draw_thresh == 0)
        cluster_draw_thresh = floor(0.05*npaths);
    end
    n = sum(cluster_sizes > cluster_draw_thresh);
    passed_clusters = find(cluster_sizes > cluster_draw_thresh);
else
    passed_clusters = 0;
end

hold off;

if exist('cl','var')
    % Iterate over all clusters
    for c = 1:length(passed_clusters)
        ci = find( cl == passed_clusters(c) );
        % Iterate over all paths in clusters
        for i = 1:length(ci)
            p = ci(i);
            npts = length(this.pathways(p).xpos(1:p_skip:end));
            % Maybe I am supposed to only draw one particular cluster out of
            % the ones that passed
            if ( only_cluster == 0 || cl(p) == only_cluster)
                label = ones(1,npts)*c;
                tubeplot(this.pathways(p).xpos(1:p_skip:end),this.pathways(p).ypos(1:p_skip:end),this.pathways(p).zpos(1:p_skip:end),0.5,label);
                hold on;
            end
        end
    end
    colormap( jet );
    caxis( [1 length(passed_clusters)+1] );
else
    % Setting up rendering parameters
    for p = 1:length(this.pathways)
%         if p > 500
%             p
%         end
        if(rand > 0 && length(this.pathways(p).xpos) >= 10)
            %tubeplot(pd(p).xpos,pd(p).ypos,pd(p).zpos,0.5,pd(p).xpos,10);
            %tubeplot(pd(p).xpos,pd(p).ypos,pd(p).zpos,0.5,ones(length(pd(p
            %).xpos),1));
            tubeplot(this.pathways(p).xpos(1:p_skip:end),this.pathways(p).ypos(1:p_skip:end),this.pathways(p).zpos(1:p_skip:end),0.5);
            hold on;
        end
    end
    colormap(jet);
end

axis([0 73*2 0 98*2 0 66*2])
daspect([1 1 1])

lighting phong;
hidden off;
%shading interp;

%plot3(45,0,50,'ro','LineWidth',10);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% LFBuild4DFreqHypercone - construct a 4D hypercone passband filter in the frequency domain
%
% Usage:
%
%     [H, FiltOptions] = LFBuild4DFreqHypercone( LFSize, Slope, BW, FiltOptions )
%     H = LFBuild4DFreqHypercone( LFSize, Slope, BW )
%
% This file constructs a real-valued magnitude response in 4D, for which the passband is a hypercone.
%
% Once constructed the filter must be applied to a light field, e.g. using LFFilt4DFFT. The
% LFDemoBasicFilt* files demonstrate how to contruct and apply frequency-domain filters.
%
% A more technical discussion, including the use of filters for denoising and volumetric focus, and
% the inclusion of aliases components, is included in:
%
% [2] D.G. Dansereau, O. Pizarro, and S. B. Williams, "Linear Volumetric Focus for Light Field
% Cameras," to appear in ACM Transactions on Graphics (TOG), vol. 34, no. 2, 2015.
%
% Inputs:
%
%     LFSize : Size of the frequency-domain filter. This should match or exceed the size of the
%     light field to be filtered. If it's larger than the input light field, the input is
%     zero-padded to match the filter's size by LFFilt4DFFT.
%
%     Slope : The slope of the planar passband. If different slopes are desired in s,t and u,v,
%     the optional aspect parameter should be used.
%
%     BW : 3-db Bandwidth of the planar passband.
%
%     [optional] FiltOptions : struct controlling filter construction
%               SlopeMethod : 'Skew' or 'Rotate' default 'skew'
%                 Precision : 'single' or 'double', default 'single'
%                   Rolloff : 'Gaussian' or 'Butter', default 'Gaussian'
%                     Order : controls the order of the filter when Rolloff is 'Butter', default 3
%                  Aspect4D : aspect ratio of the light field, default [1 1 1 1]
%                    Window : Default false. By default the edges of the passband are sharp; this adds
%                             a smooth rolloff at the edges when used in conjunction with Extent4D or
%                             IncludeAliased.
%                  Extent4D : controls where the edge of the passband occurs, the default [1 1 1 1]
%                             is the edge of the Nyquist box. When less than 1, enabling windowing
%                             introduces a rolloff after the edge of the passband. Can be greater
%                             than 1 when using IncludeAliased.
%            IncludeAliased : default false; allows the passband to wrap around off the edge of the
%                             Nyquist box; used in conjunction with Window and/or Extent4D. This can
%                             increase processing time dramatically, e.g. Extent4D = [2,2,2,2]
%                             requires a 2^4 = 16-fold increase in time to construct the filter.
%                             Useful when passband content is aliased, see [2].
%
% Outputs:
%
%                 H : real-valued frequency magnitude response
%       FiltOptions : The filter options including defaults, with an added PassbandInfo field
%                     detailing the function and time of construction of the filter
%
% User guide: LFToolbox.pdf
% See also:  LFDemoBasicFiltGantry, LFDemoBasicFiltIllum, LFDemoBasicFiltLytroF01,
% LFBuild2DFreqFan, LFBuild2DFreqLine, LFBuild4DFreqDualFan, LFBuild4DFreqHypercone,
% LFBuild4DFreqHyperfan, LFBuild4DFreqPlane, LFFilt2DFFT, LFFilt4DFFT, LFFiltShiftSum

% Copyright (c) 2013-2020 Donald G. Dansereau

function [H, FiltOptions] = LFBuild4DFreqHypercone( LFSize, BW, FiltOptions )

FiltOptions = LFDefaultField('FiltOptions', 'HyperconeMethod', 'Rotated'); % 'Direct', 'Rotated'

DistFunc = @(P, FiltOptions) DistFunc_4DCone( P, FiltOptions );
[H, FiltOptions] = LFHelperBuild4DFreq( LFSize, BW, FiltOptions, DistFunc );

TimeStamp = datestr(now,'ddmmmyyyy_HHMMSS');
FiltOptions.PassbandInfo = struct('mfilename', mfilename, 'time', TimeStamp, 'VersionStr', LFToolboxVersion);

end

%-----------------------------------------------------------------------------------------------------------------------
function Dist = DistFunc_4DCone( P, FiltOptions )
switch( lower(FiltOptions.HyperconeMethod ))
	case 'direct'
		Dist = (P(1,:).*P(4,:) - P(2,:).*P(3,:)).^2 * 4;

	case 'rotated'
		R = 1/sqrt(2) .* [1,0,0,1; 0,1,1,0; 0,1,-1,0; 1,0,0,-1];
		P = R*P;
		Dist = (sqrt(P(1,:).^2 + P(3,:).^2) - sqrt(P(2,:).^2 + P(4,:).^2)).^2 /2;

	otherwise
		error('Unrecognized hypercone construction method');
end
end

repo: LFToolbox


text:
function [seminnrank,U,V] = seminonnegativerank(M);
% [seminnrank,U,V] = seminonnegativerank(M);
%
% Computes the semi-nonnegative rank and a corresponding factorization of
% matrix M: M = UV, where U has r columns and V >= 0 has r rows and
% r = semi-nonnegative rank of M.
%
% See Corollary 3 in
% N. Gillis, A. Kumar, Exact and Heuristic Algorithms for Semi-Nonnegative
% Matrix Factorization, arXiv, 2014
%
% ****** Input ******
%   M     : m-by-n matrix
%
% ****** Output ******
%   seminnrank  : the (numerical) semi-nonnegative rank of M
%   (U,V)       : M = UV, V >= 0 and U (resp. V) has 'seminnrank' columns
%                  (resp. rows)

r = rank(M);
if issparse(M) ~= 1
    [U,S,V] = svd(M);
    V = V(:,1:r)';
else
    [U,S,V] = svds(M,r);
    V = V';
end

[y,eflag] = linsys_semiNMF(V,0);
if eflag == 1
    seminnrank = r;
    if nargout >= 2
        V = LPinitSemiNMF(M,r);
        U = M/V;
    end
else
    seminnrank = r + 1;
    if nargout >= 2
        [U,V] = SVDinitSemiNMF(M,r+1);
    end
end

author: jwyang


text:
function learningConvergence

k = 10;
iter = 1:200;

lblRate  = LowerBoundLipshitzFirstOrder(iter);
lbscRate = lowerBoundStronglyConvexFirstOrder(iter,k);
sRate = subgradientConvergence(iter);
nRate = stronglyConvexNesterovConvergence(iter,k);
qRate = quadraticConvergence(iter,1/k);
lgRate = LipschitzGradientConvergence(iter);
scRate = stronglyConvexGradientConvergence(iter,k);

colors =    {
            [0.9290, 0.6940, 0.1250];
            [0, 0.4470, 0.7410];
          	[0.8500, 0.3250, 0.0980];
            [0, 0.4470, 0.7410];
            [0.8500, 0.3250, 0.0980];
          	[0.4940, 0.1840, 0.5560];
          	[0.4660, 0.6740, 0.1880];
          	[0.3010, 0.7450, 0.9330]};
style = {'-';'--';'--';'-';'-';'-';'-';};

rates = [sRate;lblRate;lbscRate;lgRate;scRate;nRate;qRate];

figure(1)
subplot(1,2,1);
hold on
for i = 1:size(rates,1)
    plot(iter,rates(i,:),'Color',colors{i},'LineStyle',style{i});
end
legend('Subgradient','LowerBound Lipshitz','LowerBound Strongly Convex', 'Lipschitz Diff Gradient','StronglyConvex Diff Gradient','Strongly Convex Diff Nesterov','Quadratic Diff')

figure(1)
subplot(1,2,2);
hold on
for i = 1:size(rates,1)
h = plot(iter,rates(i,:),'Color',colors{i},'LineStyle',style{i});
set(gca,'yscale','log');
end

legend('Subgradient','LowerBound Lipshitz','LowerBound Strongly Convex', 'Lipschitz Diff Gradient','StronglyConvex Diff Gradient','Strongly Convex Diff Nesterov','Quadratic Diff')

end

function rate = subgradientConvergence(x)

rate = 1./sqrt(x);

end

function rate = LowerBoundLipshitzFirstOrder(x)

rate = 1./(x + 1).^2;

end

function rate = LipschitzGradientConvergence(x)

rate = 1./(x + 1);

end

function rate = stronglyConvexGradientConvergence(x,k)

rate = (((k) - 1)/((k)+1)).^(x);
end

function rate = quadraticConvergence(x,k)

rate = k.^(2.^x);

end

function rate = lowerBoundStronglyConvexFirstOrder(x,k)

rate = ((sqrt(k) - 1)/(sqrt(k) + 1)).^(2*x);

end

function rate = stronglyConvexNesterovConvergence(x,k)

rate = ((sqrt(k) - 1)/(sqrt(k)+1)).^(x);

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [H,Rest,activeUEs] = functionExampleSetup_Quadriga(L,Kdrop,B,noiseVariancedBm,Kmax,f,M,polarizations)
%This function generates the channel realizations between UEs at random
%locations and the BSs in the running example, defined in Section 4.1.3.
%BSs with cylindrical arrays and channel are generated using QuaDRiGa from
%the Fraunhofer Heinrich Hertz Institute. The Urban Microcell NLOS scenario
%is used for channel modeling. QuaDRiGa needs to be installed separately
%(http://www.quadriga-channel-model.de) and is delivered with a separate
%license. This function has been tested using QuaDRiGa version 1.4.8-571.
%
%INPUT:
%L                = Number of BSs and cells
%Kdrop            = Number of UEs to be dropped in the square around a BS
%B                = Bandwidth in Hz
%noiseVariancedBm = Noise variance in dBm
%Kmax             = Maximum number of UEs served by a BS
%f                = Pilot reuse factor, giving pilot length Kmax*f
%M                = Number of BS antennas
%polarizations    = Select number of antenna polarizations (1 or 2)
%
%OUTPUT:
%H         = M x 400 x K x L x L matrix with the channel realizations over
%            400 subcarriers at one time instance
%Rest      = M x M x K x L x L matrix with estimates of the spatial
%            correlation matrices for all UEs in the network.
%            Rest(:,:,k,j,l) is the correlation matrix for the channel
%            between UE k in cell j and the BS in cell l.
%activeUEs = Kmax x L with zeros and ones. activeUEs(k,l)==1 means that
%            pilot k is used by a UE in cell l
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%% Create new Quadriga layout

%Set irrelevant parameters
s = simulation_parameters;
s.sample_density = 1;
s.use_absolute_delays = 1;

%Set center frequency
center_frequency = 2e9;
s.center_frequency = center_frequency;

%Number of subcarriers
nbrOfSubcarriers = 400;

%Create new layout from general parameters
lay = layout(s);

%Generate BSs
lay.no_tx = L;

%Set BS heights
lay.tx_position(3,:) = 25;

%Set the length in meters of the total square area
squareLength = 1000;

%Number of BSs per dimension
nbrBSsPerDim = sqrt(L);

%% Deploy BSs

%Minimum distance between BSs and UEs
minDistance = 35;

%Distance between BSs in vertical/horizontal direction
interSiteDistance = squareLength/nbrBSsPerDim;

%Deploy BSs on the grid
locationsGridHorizontal = repmat(interSiteDistance/2:interSiteDistance:squareLength-interSiteDistance/2,[nbrBSsPerDim 1]);
locationsGridVertical = locationsGridHorizontal';
BSpositions = locationsGridHorizontal(:) + 1i*locationsGridVertical(:);

for j = 1:length(BSpositions)

    lay.tx_position(1:2,j) = [real(BSpositions(j)); imag(BSpositions(j))];

end

%% Create a circular antenna array for each BS
M_V = 5; %Number of vertical antennas
M_H = M/M_V; %Number of antennas on each horizontal circle

%Define the antenna spacing (in number of wavelengths)
antennaSpacing = 1/2; %Half wavelength distance

if polarizations == 1
    PolarizationIndicator = 1; %Single polarization antennas switching between vertical and horizontal
elseif polarizations == 2
    PolarizationIndicator = 3; %Dual +/-45deg polarized antennas
end

%Compute height of array
arrayHeight = (M_V-1)*antennaSpacing*3e8/center_frequency;

if (PolarizationIndicator==1) %Single polarized elements

    circumference = M_H*antennaSpacing*3e8/center_frequency;
    radius = circumference/(2*pi);
    delta_angle = 2*pi/M_H;

    %Go through all BSs
    for b = 1:L

        %Create rectangular array of size M_V x M_H
        lay.tx_array(b).generate('3gpp-3d', 1, M_V, M_H, center_frequency, PolarizationIndicator, 0, antennaSpacing);

        %Place antennas on a circle and rotate radiation patters
        for i = 1:M_V
            for j = 1:M_H
                indices = (i-1)*M_H + j;
                angle = (j-1)*delta_angle;
                lay.tx_array(b).element_position(1, indices) = radius*cos(angle);
                lay.tx_array(b).element_position(2, indices) = radius*sin(angle);
                lay.tx_array(b).element_position(3, indices) = (i-1)*antennaSpacing*3e8/center_frequency - arrayHeight/2;
                lay.tx_array(b).rotate_pattern(rad2deg(angle), 'z', indices, 0);

                if mod(indices,2) == 0 %Switch between vertical and horizontal polarization
                    lay.tx_array(b).rotate_pattern(90, 'y', indices, 2);
                end
            end
        end
    end

elseif (PolarizationIndicator==3) %Dual polarized elements

    circumference = M_H/2*antennaSpacing*3e8/center_frequency;
    radius = circumference/(2*pi);
    delta_angle = 2*pi/(M_H/2);

    %Go through all BSs
    for b = 1:L

        %Create rectangular array of size M_V x M_H/2
        lay.tx_array(b).generate('3gpp-3d', 1, M_V, M_H/2, center_frequency, PolarizationIndicator, 0, antennaSpacing);

        %Place antennas on a circle and rotate radiation patters (while keeping co-located antennas together)
        for i = 1:M_V
            for j = 1:M_H/2
                indices = (i-1)*M_H + 2*j-1 : (i-1)*M_H + 2*j;
                angle = (j-1)*delta_angle;
                lay.tx_array(b).element_position(1, indices) = radius*cos(angle);
                lay.tx_array(b).element_position(2, indices) = radius*sin(angle);
                lay.tx_array(b).element_position(3, indices) = (i-1)*antennaSpacing*3e8/center_frequency - arrayHeight/2;
                lay.tx_array(b).rotate_pattern(rad2deg(angle), 'z', indices, 0);
            end
        end
    end
end

%Compute all nine alternatives of the BS locations when using wrap around
wrapHorizontal = repmat([-squareLength 0 squareLength],[3 1]);
wrapVertical = wrapHorizontal';
wrapLocations = wrapHorizontal(:)' + 1i*wrapVertical(:)';

%Compute the exact dimension of the square where the users are located
maxDistance = interSiteDistance;

%Prepare to put out UEs in the cells
UEpositions = zeros(Kdrop,L);
UEpositionsWrapped = zeros(Kdrop,L,length(wrapLocations));
perBS = zeros(L,1);

%Go through all the cells
for l = 1:L

    %Put out K UEs in the cell, uniformly at random. The procedure is
    %iterative since UEs that do not satisfy the minimum distance are
    %replaced with new UEs
    while perBS(l)<Kdrop

        %Put out users
        UEremaining = Kdrop-perBS(l);
        posX = rand(UEremaining,1)*maxDistance - maxDistance/2;
        posY = rand(UEremaining,1)*maxDistance - maxDistance/2;
        posXY = posX + 1i*posY;

        %Keep those that satisfy the minimum distance
        posXY = posXY(abs(posXY)>=minDistance);

        %Store new UEs
        UEpositions(perBS(l)+1:perBS(l)+length(posXY),l) = posXY + BSpositions(l);
        perBS(l) = perBS(l)+length(posXY);

    end

    %Create alternative UE positions using wrap around
    for k = 1:Kdrop

        UEpositionsWrapped(k,l,:) = UEpositions(k,l) + wrapLocations;

    end

end

%% Configure UEs

Ktotal = Kdrop*L*length(wrapLocations); %Total number of UEs

%Define UE heights
UE_heights = 1.5*ones(Kdrop,L);
UE_heightsWrapped = repmat(UE_heights,[1 1 length(wrapLocations)]);

%Generate UEs
lay.no_rx = Ktotal;

%Define UE antennas
lay.rx_array.generate('omni');

%% Simulate channels

% Randomly distribute UEs
lay.rx_position =  [real(UEpositionsWrapped(:))'; imag(UEpositionsWrapped(:))'; UE_heightsWrapped(:)'];

%Define tracks for each UE, assuming a fixed UE location
for k=1:Ktotal
    lay.track(k).generate('linear',0,0) %Define a linear track consisting of only one position
    lay.track(k).scenario = '3GPP_3D_UMi_NLOS'; %Select the Urban Microcell NLOS scenario
end

%Generate pilot patterns
if f == 1

    pilotPattern = ones(L,1);

elseif f == 2 %Only works for 16 BSs

    pilotPattern = kron(ones(2,1),[1; 2; 1; 2; 2; 1; 2; 1]);

end

%Randomize pilot allocation in each cell
randOrder = zeros(Kmax*f,L);

for j = 1:L

    randOrder(1+(pilotPattern(j)-1)*Kmax:pilotPattern(j)*Kmax,j) = randperm(Kmax)+(pilotPattern(j)-1)*Kmax;

end

%Compute variance and standard deviation of the noise
noiseVar = 10^(noiseVariancedBm/10);
noiseStd = sqrt(noiseVar);

%Prepare to store channel realizations
H = zeros(M,nbrOfSubcarriers,Kmax,L,L);
Rest = zeros(M,M,Kmax,L,L);
perBS = zeros(L,1);
activeUEs = zeros(Kmax,L);

%% Go through all cells
for j = 1:L

    %Output simulation progress
    disp([num2str(j) ' cells generated out of ' num2str(L)]);

    %Go through all UEs
    for k = 1:Kdrop

        Huser = zeros(M,nbrOfSubcarriers,1,1,L);
        Ruser = zeros(M,M,1,1,L);

        %Extract the channels to all BSs
        for l = 1:L

            [~,minr] = min(abs(UEpositionsWrapped(k,j,:)-BSpositions(l)));

            userind = k+(j-1)*Kdrop+(minr-1)*Kdrop*L;

            [ h_channel, ~ ] = lay.get_channels_seg(l, userind);
            Hextract = h_channel.fr(B, nbrOfSubcarriers);

            Huser(:,:,1,1,l) = reshape(Hextract,[M nbrOfSubcarriers])/noiseStd;
            Ruser(:,:,1,1,l) = diag(mean(abs(Huser(:,:,1,1,l)).^2,2)/noiseVar);

        end

        %Determine which BS should serve the UE
        [~,bestBS] = max(mean(sum(abs(Huser(:,:,1,1,:)).^2,1),2));

        %Check if the selected BS has pilots available
        if perBS(bestBS)<Kmax

            %Add the UE to the cell of the selected BS
            perBS(bestBS) = perBS(bestBS) + 1;
            H(:,:,randOrder(perBS(bestBS)+(pilotPattern(bestBS)-1)*Kmax,bestBS),bestBS,:) = Huser;
            Rest(:,:,randOrder(perBS(bestBS)+(pilotPattern(bestBS)-1)*Kmax,bestBS),bestBS,:) = Ruser;
            activeUEs(randOrder(perBS(bestBS)+(pilotPattern(bestBS)-1)*Kmax,bestBS),bestBS) = 1;

        end

    end

end

author: emilbjornson

repo: massivemimobook


text:
function [z,p,sig,pt] = fisherp(p,varargin)
% :Usage:
% ::
%
%     function [z,p,sig,pt] = fisherp(p,[alph])
%
% :Inputs:
%
%   **p:**
%        values in 4-D array
%
%        1st 3 dims are within images, dim4 = image
%
% :Optional: alpha value for thresholding
%
% :Outputs:
%
%   **z:**
%        Fisher's combined test statistic, compare to normal
%
%   **p:**
%        p-values for combined test
%
%   **sig:**
%        signficance 1 / 0 binary mask, p < .05 (or alph) FDR-corr
%
%   **pt:**
%        p-value threshold for FDR corrected significance at alph
%
% :Described in:
% Lazar, N. A., Luna, B., Sweeney, J. A., & Eddy, W. F. (2002).
% Combining brains: a survey of methods for statistical pooling
% of information. Neuroimage, 16(2), 538-550.
%
% Stouffer, S. A., Suchman, E. A., DeVinney, L. C., Star, S. A., and
% Williams, R. M. 1949. The American Soldier: Vol. I. Adjustment
% During Army Life. Princeton University Press, Princeton.
%
% Threshold is determined with False Discovery Rate (Benjamini & Hochberg, 1995)
%
% ..
%    tor wager
% ..

if length(varargin) > 0, alph = varargin{1};,else,alph = 0.05;,end

lastdim = length(size(p));
k = size(p,lastdim);

z = -2*sum(log(p),lastdim);
p = 1 - chi2cdf(z(:),2*k);  % distributed as chi-square with 2k df
p = reshape(p,size(z));

% eliminate voxels outside of brain, all 0, 1, or all NaN values
p(all(p == 0,lastdim)) = NaN;
p(all(p == 1,lastdim)) = NaN;

pp = p(:); pp(isnan(pp)) = [];
pt = FDR(pp,alph);
if isempty(pt), pt = -Inf;, end

%z = sum(norminv(1 - p),lastdim) ./ sqrt(k);
%p = normcdf(1-z);

%if alph, sig = p <= alph;,end
sig = p <= pt;

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ n_data, n, c ] = tau_values ( n_data )

%*****************************************************************************80
%
%% TAU_VALUES returns some values of the Tau function.
%
%  Discussion:
%
%    TAU(N) is the number of divisors of N, including 1 and N.
%
%    In Mathematica, the function can be evaluated by:
%
%      DivisorSigma[1,n]
%
%  First values:
%
%     N   TAU(N)
%
%     1    1
%     2    2
%     3    2
%     4    3
%     5    2
%     6    4
%     7    2
%     8    4
%     9    3
%    10    4
%    11    2
%    12    6
%    13    2
%    14    4
%    15    4
%    16    5
%    17    2
%    18    6
%    19    2
%    20    6
%
%  Formula:
%
%    If the prime factorization of N is
%
%      N = P1^E1 * P2^E2 * ... * PM^EM,
%
%    then
%
%      TAU(N) = ( E1 + 1 ) * ( E2 + 1 ) * ... * ( EM + 1 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the argument of the Tau function.
%
%    Output, integer C, the value of the Tau function.
%
  n_max = 20;

  c_vec = [ ...
    1,  2,  2,  3,  2,  4,  2,  4,  3,  4, ...
    2, 12, 12,  4, 18, 24,  2,  8, 14, 28 ];

  n_vec = [ ...
      1,   2,   3,   4,   5,   6,   7,   8,   9,  10, ...
     23,  72, 126, 226, 300, 480, 521, 610, 832, 960 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    c = 0;
  else
    n = n_vec(n_data);
    c = c_vec(n_data);
  end

  return
end


text:
function errorMsg = import_anatomy_bs(iSubject, BsDir, nVertices, isInteractive, sFid, isVolumeAtlas, isKeepMri)
% IMPORT_ANATOMY_BS: Import a full BrainSuite folder as the subject's anatomy.
%
% USAGE:  errorMsg = import_anatomy_bs(iSubject, BsDir=[ask], nVertices=[ask], isInteractive=1, sFid=[], isVolumeAtlas=1, isKeepMri=0)
%
% INPUT:
%    - iSubject      : Indice of the subject where to import the MRI
%                      If iSubject=0 : import MRI in default subject
%    - BsDir         : Full filename of the BrainSuite folder to import
%    - nVertices     : Number of vertices in the file cortex surface
%    - isInteractive : If 0, no input or user interaction
%    - sFid          : Structure with the fiducials coordinates
%                      Or full MRI structure with fiducials defined in the SCS structure, to be registered with the FS MRI
%    - isVolumeAtlas : If 1, imports the svreg atlas as a set of surfaces
%    - isKeepMri     : 0=Delete all existing anatomy files
%                      1=Keep existing MRI volumes (when running segmentation from Brainstorm)
%                      2=Keep existing MRI and surfaces
% OUTPUT:
%    - errorMsg : String: error message if an error occurs

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author   : Francois Tadel, 2012-2022
% Modified : Andrew Krause, 2013

%% ===== PARSE INPUTS =====
% Keep MRI
if (nargin < 7) || isempty(isKeepMri)
    isKeepMri = 0;
end
% Import volume atlas
if (nargin < 6) || isempty(isVolumeAtlas)
    isVolumeAtlas = 1;
end
% Fiducials
if (nargin < 5) || isempty(sFid)
    sFid = [];
end
% Interactive / silent
if (nargin < 4) || isempty(isInteractive)
    isInteractive = 1;
end
% Ask number of vertices for the cortex surface
if (nargin < 3) || isempty(nVertices)
    nVertices = [];
end
% Initialize returned variables
errorMsg = [];
% Ask folder to the user
if (nargin < 2) || isempty(BsDir)
    % Get default import directory and formats
    LastUsedDirs = bst_get('LastUsedDirs');
    % Open file selection dialog
    BsDir = java_getfile( 'open', ...
        'Import BrainSuite folder...', ...     % Window title
        bst_fileparts(LastUsedDirs.ImportAnat, 1), ...           % Last used directory
        'single', 'dirs', ...                  % Selection mode
        {{'.folder'}, 'BrainSuite folder', 'BsDir'}, 0);
    % If no folder was selected: exit
    if isempty(BsDir)
        return
    end
    % Save default import directory
    LastUsedDirs.ImportAnat = BsDir;
    bst_set('LastUsedDirs', LastUsedDirs);
end
% Unload everything
bst_memory('UnloadAll', 'Forced');

%% ===== DELETE PREVIOUS ANATOMY =====
% Get subject definition
sSubject = bst_get('Subject', iSubject);
% Check for existing anatomy
if (~isempty(sSubject.Anatomy) && (isKeepMri == 0)) || (~isempty(sSubject.Surface) && (isKeepMri < 2))
    % Ask user whether the previous anatomy should be removed
    if isInteractive
        isDel = java_dialog('confirm', ['Warning: There is already an anatomy defined for this subject.' 10 10 ...
            'Are you sure you want to delete the previous MRI and surfaces ?' 10 10], 'Import BrainSuite folder');
    else
        isDel = 1;
    end
    % If user canceled process
    if ~isDel
        bst_progress('stop');
        return;
    end
    % Delete anatomy
    sSubject = db_delete_anatomy(iSubject, isKeepMri);
end

%% ===== ASK NB VERTICES =====
if isempty(nVertices)
    nVertices = java_dialog('input', 'Number of vertices on the cortex surface:', 'Import BrainSuite folder', [], '15000');
    if isempty(nVertices)
        return
    end
    nVertices = str2double(nVertices);
end
% Number for each hemisphere
nVertHemi = round(nVertices / 2);

%% ===== PARSE BRAINSUITE FOLDER =====
% Find MRI
FilePrefix = get_fileprefix(BsDir);
if isempty(FilePrefix)
    errorMsg = [errorMsg 'Could not determine file prefix from BFC file' 10 10 ...
        'Make sure SVREG was executed in BrainSuite before importing the data in Brainstorm.'];
    if isInteractive
        bst_error(['Could not import BrainSuite folder: ' 10 10 errorMsg], 'Import BrainSuite folder', 0);
    end
    return;
end
T1File = {file_find(BsDir, [FilePrefix '.nii.gz']), ...
           file_find(BsDir, [FilePrefix '.nii']), ...
           file_find(BsDir, [FilePrefix '.img.gz']),...
           file_find(BsDir, [FilePrefix '.img'])};
T1File = [T1File{find(~cellfun(@isempty, T1File))}];
if isempty(T1File)
    errorMsg = [errorMsg 'MRI file was not found: ' FilePrefix '.*' 10];
end

% Find volume segmentation file
BsDirMultiParc = fullfile(BsDir,'multiparc');

SvregFile = file_find(BsDir, [FilePrefix '.svreg.label.nii.gz']);
OtherSvregFiles = file_find(BsDirMultiParc, [FilePrefix '.svreg.*.label.nii.gz'], 2, 0);

% Find surfaces
HeadFile        = file_find(BsDir, [FilePrefix '.scalp.dfs']);
InnerSkullFile  = file_find(BsDir, [FilePrefix '.inner_skull.dfs']);
OuterSkullFile  = file_find(BsDir, [FilePrefix '.outer_skull.dfs']);
TessLhFile      = file_find(BsDir, [FilePrefix '.left.pial.cortex.svreg.dfs']);
TessRhFile      = file_find(BsDir, [FilePrefix '.right.pial.cortex.svreg.dfs']);
TessLwFile      = file_find(BsDir, [FilePrefix '.left.inner.cortex.svreg.dfs']);
TessRwFile      = file_find(BsDir, [FilePrefix '.right.inner.cortex.svreg.dfs']);
TessLsphFile    = file_find(BsDir, [FilePrefix '.left.mid.cortex.svreg.dfs']);
TessRsphFile    = file_find(BsDir, [FilePrefix '.right.mid.cortex.svreg.dfs']);
TessLAtlsphFile = file_find(BsDir, 'atlas.left.mid.cortex.svreg.dfs');
TessRAtlsphFile = file_find(BsDir, 'atlas.right.mid.cortex.svreg.dfs');

% Find labels
AnnotLhFiles = file_find(BsDirMultiParc, [FilePrefix '.left.mid.cortex.svreg.*.dfs'], 2, 0);
AnnotRhFiles = file_find(BsDirMultiParc, [FilePrefix '.right.mid.cortex.svreg.*.dfs'], 2, 0);

if isempty(AnnotLhFiles)
    AnnotLhFiles = TessLhFile;
end

if isempty(AnnotRhFiles)
    AnnotRhFiles = TessRhFile;
end

if isempty(HeadFile)
    %errorMsg = [errorMsg 'Scalp file was not found: ' FilePrefix '.scalp.dfs' 10];
    disp(['BST> Warning: Scalp file was not found: ' FilePrefix '.scalp.dfs']);
end
if isempty(InnerSkullFile) && isempty(OuterSkullFile)
    %errorMsg = [errorMsg 'Inner or Outer Skull File not found' 10];
    disp('BST> Warning: Inner or Outer Skull File not found.');
end
if isempty(TessLhFile)
    errorMsg = [errorMsg 'Surface file was not found: ' FilePrefix '.left.pial.cortex.svreg.dfs' 10];
end
if isempty(TessRhFile)
    errorMsg = [errorMsg 'Surface file was not found: ' FilePrefix '.right.pial.cortex.svreg.dfs' 10];
end

% Report errors
if ~isempty(errorMsg)
    if isInteractive
        bst_error(['Could not import BrainSuite folder: ' 10 10 errorMsg], 'Import BrainSuite folder', 0);
    end
    return;
end

%% ===== IMPORT MRI =====
if isKeepMri && ~isempty(sSubject.Anatomy)
    BstT1File = file_fullpath(sSubject.Anatomy(sSubject.iAnatomy).FileName);
else
    bst_progress('text', 'Loading MRI...');
    % Read MRI
    [BstT1File, sMri] = import_mri(iSubject, T1File);
    if isempty(BstT1File)
        errorMsg = 'Could not import BrainSuite folder: MRI was not imported properly';
        if isInteractive
            bst_error(errorMsg, 'Import BrainSuite folder', 0);
        end
        return;
    end
    % Enforce it as the permanent default MRI
    sSubject = db_surface_default(iSubject, 'Anatomy', 1, 0);
end

%% ===== DEFINE FIDUCIALS / MNI NORMALIZATION =====
% Set fiducials and/or compute linear MNI normalization
[isComputeMni, errCall] = process_import_anatomy('SetFiducials', iSubject, BsDir, BstT1File, sFid, isKeepMri, isInteractive);
% Error handling
if ~isempty(errCall)
    errorMsg = [errorMsg, errCall];
    if isempty(isComputeMni)
        if isInteractive
            bst_error(errorMsg, 'Import BrainSuite folder', 0);
        end
        return;
    end
end

%% ===== IMPORT SURFACES =====
% Left pial
if ~isempty(TessLhFile)
    % Import file
    [iLh, BstTessLhFile, nVertOrigL] = import_surfaces(iSubject, TessLhFile, 'DFS', 0);
    BstTessLhFile = BstTessLhFile{1};
    % Load atlas
    if ~isempty(AnnotLhFiles)
        bst_progress('start', 'Import BrainSuite folder', 'Loading atlas: left pial...');
        [sAllAtlas, err] = import_label(BstTessLhFile, AnnotLhFiles, 1);
        errorMsg = [errorMsg err];
    end
    % Load registration square
    if ~isempty(TessLsphFile)
        bst_progress('start', 'Import BrainSuite folder', 'Loading registered square: left pial...');
        [TessMat, err] = tess_addsquare(BstTessLhFile, TessLsphFile, TessLAtlsphFile);
        errorMsg = [errorMsg err];
    end
    % Downsample
    bst_progress('start', 'Import BrainSuite folder', 'Downsampling: left pial...');
    [BstTessLhLowFile, iLhLow, xLhLow] = tess_downsize(BstTessLhFile, nVertHemi, 'reducepatch');
end
% Right pial
if ~isempty(TessRhFile)
    % Import file
    [iRh, BstTessRhFile, nVertOrigR] = import_surfaces(iSubject, TessRhFile, 'DFS', 0);
    BstTessRhFile = BstTessRhFile{1};
    % Load atlas
    if ~isempty(AnnotRhFiles)
        bst_progress('start', 'Import BrainSuite folder', 'Loading atlas: right pial...');
        [sAllAtlas, err] = import_label(BstTessRhFile, AnnotRhFiles, 1);
        errorMsg = [errorMsg err];
    end
    % Load registration square
    if ~isempty(TessRsphFile)
        bst_progress('start', 'Import BrainSuite folder', 'Loading registered square: right pial...');
        [TessMat, err] = tess_addsquare(BstTessRhFile, TessRsphFile, TessRAtlsphFile);
        errorMsg = [errorMsg err];
    end
    % Downsample
    bst_progress('start', 'Import BrainSuite folder', 'Downsampling: right pial...');
    [BstTessRhLowFile, iRhLow, xRhLow] = tess_downsize(BstTessRhFile, nVertHemi, 'reducepatch');
end
% Left white matter
if ~isempty(TessLwFile)
    % Import file
    [iLw, BstTessLwFile] = import_surfaces(iSubject, TessLwFile, 'DFS', 0);
    BstTessLwFile = BstTessLwFile{1};
    % Load atlas
    if ~isempty(AnnotLhFiles)
        bst_progress('start', 'Import BrainSuite folder', 'Loading atlas: left white...');
        [sAllAtlas, err] = import_label(BstTessLwFile, AnnotLhFiles, 1);
        errorMsg = [errorMsg err];
    end
    % Load registration square
    if ~isempty(TessLsphFile)
        bst_progress('start', 'Import BrainSuite folder', 'Loading registered square: left white...');
        [TessMat, err] = tess_addsquare(BstTessLwFile, TessLsphFile, TessLAtlsphFile);
        errorMsg = [errorMsg err];
    end
    % Downsample
    bst_progress('start', 'Import BrainSuite folder', 'Downsampling: left white...');
    [BstTessLwLowFile, iLwLow] = tess_downsize(BstTessLwFile, nVertHemi, 'reducepatch');
end
% Right white matter
if ~isempty(TessRwFile)
    % Import file
    [iRw, BstTessRwFile] = import_surfaces(iSubject, TessRwFile, 'DFS', 0);
    BstTessRwFile = BstTessRwFile{1};
     % Load atlas
    if ~isempty(AnnotRhFiles)
        bst_progress('start', 'Import BrainSuite folder', 'Loading atlas: right inner...');
        [sAllAtlas, err] = import_label(BstTessRwFile, AnnotRhFiles, 1);
        errorMsg = [errorMsg err];
    end
    % Load registration square
    if ~isempty(TessRsphFile)
        bst_progress('start', 'Import BrainSuite folder', 'Loading registered square: right white...');
        [TessMat, err] = tess_addsquare(BstTessRwFile, TessRsphFile, TessRAtlsphFile);
        errorMsg = [errorMsg err];
    end
    % Downsample
    bst_progress('start', 'Import BrainSuite folder', 'Downsampling: right white...');
    [BstTessRwLowFile, iRwLow] = tess_downsize(BstTessRwFile, nVertHemi, 'reducepatch');
end
% Process error messages
if ~isempty(errorMsg)
    if isInteractive
        bst_error(errorMsg, 'Import BrainSuite folder', 0);
    end
    return;
end

%% ===== MERGE SURFACES =====
rmFiles = {};
rmInd   = [];
% Merge hemispheres: pial
if ~isempty(TessLhFile) && ~isempty(TessRhFile)
    % Hi-resolution surface
    CortexHiFile  = tess_concatenate({BstTessLhFile,    BstTessRhFile},    sprintf('cortex_%dV', nVertOrigL + nVertOrigR), 'Cortex');
    CortexLowFile = tess_concatenate({BstTessLhLowFile, BstTessRhLowFile}, sprintf('cortex_%dV', nVertices), 'Cortex');
    % Delete separate hemispheres
    rmFiles = cat(2, rmFiles, {BstTessLhFile, BstTessRhFile, BstTessLhLowFile, BstTessRhLowFile});
    rmInd   = [rmInd, iLh, iRh, iLhLow, iRhLow];
end
% Merge hemispheres: white
if ~isempty(TessLwFile) && ~isempty(TessRwFile)
    % Hi-resolution surface
    WhiteHiFile  = tess_concatenate({BstTessLwFile,    BstTessRwFile},    sprintf('white_%dV', nVertOrigL + nVertOrigR), 'Cortex');
    WhiteLowFile = tess_concatenate({BstTessLwLowFile, BstTessRwLowFile}, sprintf('white_%dV', nVertices), 'Cortex');
    % Delete separate hemispheres
    rmFiles = cat(2, rmFiles, {BstTessLwFile, BstTessRwFile, BstTessLwLowFile, BstTessRwLowFile});
    rmInd   = [rmInd, iLw, iRw, iLwLow, iRwLow];
end

%% ===== HEAD AND SKULL SURFACES =====
bst_progress('start', 'Import BrainSuite folder', 'Importing scalp and skull surfaces...');
% Head
if ~isempty(HeadFile)
    % Import file
    bst_progress('start', 'Import BrainSuite folder', 'Imported scalp surface...');
    [iHead, BstHeadHiFile] = import_surfaces(iSubject, HeadFile, 'DFS', 0);
    BstHeadHiFile = BstHeadHiFile{1};
    % Downsample
    bst_progress('start', 'Import BrainSuite folder', 'Downsampling: scalp...');
    BstHeadFile = tess_downsize( BstHeadHiFile, 1082, 'reducepatch' );
    % Load MRI
    bst_progress('start', 'Import BrainSuite folder', 'Filling holes in the head surface...');
    sMri = in_mri_bst(BstT1File);
    % Load head surface
    sHead = in_tess_bst(BstHeadFile);
    % Remove holes
    [sHead.Vertices, sHead.Faces] = tess_fillholes(sMri, sHead.Vertices, sHead.Faces, 2, 2);
    % Save back to file
    sHeadNew.Vertices = sHead.Vertices;
    sHeadNew.Faces = sHead.Faces;
    sHeadNew.Comment = sHead.Comment;
    bst_save(file_fullpath(BstHeadFile), sHeadNew, 'v7');
    % Delete initial file
    rmFiles = cat(2, rmFiles, BstHeadHiFile);
    rmInd   = [rmInd, iHead];
% Or generate one from Brainstorm
else
    % Generate head surface
    BstHeadFile = tess_isohead(iSubject, 10000, 0, 2);
end

% Inner Skull
if ~isempty(InnerSkullFile)
    % Import file
    [iIs, BstInnerSkullHiFile] = import_surfaces(iSubject, InnerSkullFile, 'DFS', 0);
    BstInnerSkullHiFile = BstInnerSkullHiFile{1};
    % Downsample
    bst_progress('start', 'Import BrainSuite folder', 'Downsampling: inner skull...');
    BstInnerSkullFile = tess_downsize(BstInnerSkullHiFile, 1000, 'reducepatch');
    % Delete initial file
    rmFiles = cat(2, rmFiles, BstInnerSkullHiFile);
    rmInd   = [rmInd, iIs];
end
if ~isempty(OuterSkullFile)
    % Import file
    [iOs, BstOuterSkullHiFile] = import_surfaces(iSubject, OuterSkullFile, 'DFS', 0);
    BstOuterSkullHiFile = BstOuterSkullHiFile{1};
    % Downsample
    bst_progress('start', 'Import BrainSuite folder', 'Downsampling: outer skull...');
    BstOuterSkullFile = tess_downsize(BstOuterSkullHiFile, 1000, 'reducepatch');
    % Delete initial file
    rmFiles = cat(2, rmFiles, BstOuterSkullHiFile);
    rmInd   = [rmInd, iOs];
end

% Delete intermediary files
if ~isempty(rmFiles)
    % Delete files
    file_delete(file_fullpath(rmFiles), 1);
    % Update subject definition
    sSubject = bst_get('Subject', iSubject);
    sSubject.Surface(rmInd) = [];
    bst_set('Subject', iSubject, sSubject);
    % Refresh tree
    panel_protocols('UpdateNode', 'Subject', iSubject);
    panel_protocols('SelectNode', [], 'subject', iSubject, -1 );
end

%% ===== IMPORT SVREG ATLAS =====
if isVolumeAtlas && ~isempty(SvregFile)
    % Import atlas as volume
    [BstSvregFile, sMriSvreg] = import_mri(iSubject, SvregFile, 'ALL-ATLAS', 0, 1, 'svreg');
    % Import other label volumes
    for iFile = 1:length(OtherSvregFiles)
        st=strfind(OtherSvregFiles{iFile},'.svreg.');
        ed=strfind(OtherSvregFiles{iFile},'.label.nii.gz');
        AtlasName = OtherSvregFiles{iFile}(st+7:ed-1);
        import_mri(iSubject, OtherSvregFiles{iFile}, 'ALL-ATLAS', 0, 1, AtlasName);
    end
    % Import atlas
    SelLabels = {...
        'Accumbens L', 'Hippocampus L', 'Pallidum L', 'Putamen L', 'Thalamus L', ...
        'Accumbens R', 'Hippocampus R', 'Pallidum R', 'Putamen R', 'Thalamus R', ...
        'Brainstem', 'Cerebellum'};
    [iSvreg, BstSvregFile] = import_surfaces(iSubject, SvregFile, 'MRI-MASK', 0, [], SelLabels, 'subcortical');
    % Extract cerebellum only
    try
        BstCerebFile = tess_extract_struct(BstSvregFile{1}, {'Cerebellum'}, 'svreg | cerebellum');
    catch
        BstCerebFile = [];
    end
    % If the cerebellum surface can be reconstructed
    if ~isempty(BstCerebFile)
        % Downsample cerebllum
        [BstCerebLowFile, iCerLow, xCerLow] = tess_downsize(BstCerebFile, 2000, 'reducepatch');
        % Merge with low-resolution pial
        BstMixedLowFile = tess_concatenate({CortexLowFile, BstCerebLowFile}, sprintf('cortex_cereb_%dV', length(xLhLow) + length(xRhLow) + length(xCerLow)), 'Cortex');
        % Rename mixed file
        oldBstMixedLowFile = file_fullpath(BstMixedLowFile);
        BstMixedLowFile    = bst_fullfile(bst_fileparts(oldBstMixedLowFile), 'tess_cortex_pialcereb_low.mat');
        file_move(oldBstMixedLowFile, BstMixedLowFile);
        % Delete intermediate files
        file_delete({file_fullpath(BstCerebFile), file_fullpath(BstCerebLowFile)}, 1);
        db_reload_subjects(iSubject);
    end
end

%% ===== UPDATE GUI =====
% Set default cortex
if ~isempty(TessLhFile) && ~isempty(TessRhFile)
    [sSubject, iSubject, iSurface] = bst_get('SurfaceFile', CortexLowFile);
    db_surface_default(iSubject, 'Cortex', iSurface);
end
% Set default scalp
db_surface_default(iSubject, 'Scalp');
% Set default skulls
db_surface_default(iSubject, 'OuterSkull');
db_surface_default(iSubject, 'InnerSkull');
% Update subject node
panel_protocols('UpdateNode', 'Subject', iSubject);
% Save database
db_save();
% Unload everything
bst_memory('UnloadAll', 'Forced');
% Give a graphical output for user validation
if isInteractive
    % Display the downsampled cortex and the head
    hFig = view_surface(BstHeadFile);
    view_surface(CortexLowFile);
    % Set orientation
    figure_3d('SetStandardView', hFig, 'left');
end
% Close progress bar
bst_progress('stop');

end

%% ======================================================================================
%  ===== HELPER FUNCTIONS ===============================================================
%  ======================================================================================
%% ===== GET PREFIX OF FILENAMES =====
function FilePrefix = get_fileprefix(BsDir)
    % Determine file prefix based on left cortex file
    TestFile = file_find(BsDir, '*.left.pial.cortex.svreg.dfs');
    if ~isempty(TestFile)
        [tmp, FilePrefix] = bst_fileparts(TestFile(1:end-27));
    else
        FilePrefix = [];
    end
end

author: brainstorm-tools


text:
function imdb = setupScene67(datasetDir, varargin)
% SETUPSCENE67    Setup MIT Scene 67 dataset
%    This is similar to SETUPGENERIC(), with modifications to setup
%    MIT Scene 67 according to the standard evaluation protocols. The
%    function supports only the LITE option.
%
%    See: SETUPGENERIC().

% Author: Andrea Vedaldi

% Copyright (C) 2013 Andrea Vedaldi
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

opts.lite = false ;
opts = vl_argparse(opts, varargin) ;

% Download and unpack
vl_xmkdir(datasetDir) ;
if exist(fullfile(datasetDir, 'Images', 'airport_inside'))
  % ok
else
  url = 'http://groups.csail.mit.edu/vision/LabelMe/NewImages/indoorCVPR_09.tar' ;
  trainImagesUrl = 'http://web.mit.edu/torralba/www/TrainImages.txt' ;
  testImagesUrl = 'http://web.mit.edu/torralba/www/TestImages.txt' ;
  fprintf('Downloading MIT Scene 67 data to ''%s''. This will take a while.', datasetDir) ;
  urlwrite(trainImagesUrl, fullfile(datasetDir, 'TrainImages.txt')) ;
  urlwrite(testImagesUrl, fullfile(datasetDir, 'TestImages.txt')) ;
  untar(url, datasetDir) ;
end

% Construct image database imdb structure
imdb.meta.sets = {'train', 'val', 'test'} ;

trainNames = textread(fullfile(datasetDir, 'TrainImages.txt'),'%s','delimiter','\n') ;
testNames = textread(fullfile(datasetDir, 'TestImages.txt'),'%s','delimiter','\n') ;
names = [trainNames; testNames]' ;
sets = [ones(1,numel(trainNames)), 3*ones(1,numel(testNames))] ;
imdb.images.id = 1:numel(names) ;
[imdb.images.name,perm] = sort(names) ;
imdb.images.set = sets(perm) ;

a = regexp(imdb.images.name, '^([^/]+)*/.*$', 'tokens') ;
for i = 1:numel(names)
  a{i} = a{i}{1}{1} ;
end
[imdb.meta.classes, ~, imdb.images.class] = unique(a) ;
imdb.images.class = imdb.images.class(:)' ;
imdb.imageDir = fullfile(datasetDir, 'Images') ;

if opts.lite
  ok = {} ;
  for c = 1:3
    ok{end+1} = vl_colsubset(find(imdb.images.class == c & imdb.images.set == 1), 5) ;
    ok{end+1} = vl_colsubset(find(imdb.images.class == c & imdb.images.set == 3), 5) ;
  end
  ok = cat(2, ok{:}) ;
  imdb.meta.classes = imdb.meta.classes(1:3) ;
  imdb.images.id = imdb.images.id(ok) ;
  imdb.images.name = imdb.images.name(ok) ;
  imdb.images.set = imdb.images.set(ok) ;
  imdb.images.class = imdb.images.class(ok) ;
end

author: jianxiongxiao

repo: ProfXkit


text:
function y = inter_xls(Y,x,ta,sc,type,ip,d,flax1,flax2,file_name)
% PURPOSE: Interface via Excel Link for univariate temporal disaggregation
% -----------------------------------------------------------------------
% SYNTAX: y = inter_xls(Y,x,ta,sc,type,ip,d,flax1,flax2,file_name);
% -----------------------------------------------------------------------
% INPUT
%           SELECTION OF THE METHOD
%
% Common parameters:
%        ta: type of disaggregation
%            ta=1 ---> sum (flow)
%            ta=2 ---> average (index)
%            ta=3 ---> last element (stock) ---> interpolation
%            ta=4 ---> first element (stock) ---> interpolation
%        sc: number of high frequency data points for each low frequency data point
%            sc= 4 ---> annual to quarterly
%            sc=12 ---> annual to monthly
%            sc= 3 ---> quarterly to monthly
%
% Specific parameters:
%
% ==> Boot-Feibes-Lisman, Denton (additive variant, standard solution):
%        d: objective function to be minimized: volatility of ...
%            d=0 ---> levels (only BFL)
%            d=1 ---> first differences
%            d=2 ---> second differences
%
% ==> Chow-Lin, Fernandez, Litterman, Santos-Cardoso:
%           opC = -1 ---> pretesting intercept signficiance
%
% ==> Chow-Lin, Litterman, Santos-Cardoso:
%        type: estimation method:
%            type=0 ---> weighted least squares
%            type=1 ---> maximum likelihood
%        innovational parameter rl = []. Default: range [.05 .99], 100 points grid.
%
% ==> Chow-Lin, Litterman, Santos-Cardoso when innovational parameter is supplied:
%         ip, such as -1 < ip < 1
%
%           INPUT DATA:
%
% Common:
%       Y: Nx1 --> Low-frequency time series (to be temporally disaggregated)
% Specific:
%       x: nx1 --> Denton, n=sc*N (extrapolation is not feasible)
%       x: nxp --> Fernandez, Chow-Lin, Litterman, Santos-Cardoso
%          p >= 1, n >= sc*N (extrapolation is feasible)
%
% -----------------------------------------------------------------------
% OUTPUT: y: nxi
%
%       i=1 brief --> only temporally disaggregated series (all procedures)
%       i=5 normal --> temporally disaggregated series, standard errors of estimates,
%                  one-sigma upper and lower limits and residuals.
%                  Available for Fernandez, Chow-Lin, Litterman, Santos-Cardoso
%       i=5 detailed --> normal + ASCII file with model results (all
%                   procedures). A name for the output file should be supplied.
%
% -----------------------------------------------------------------------
% LIBRARY: bfl, denton_uni, fernandez, chowlin, litterman, ssc,
% td_uni_print, td_print

% written by:
% Ana Abad(*) & Enrique M. Quilis(**)
%   (*) National Statistical Institute
%   (**) Macroeconomic Research Department
%        Ministry of Economy and Competitiveness
%        <enrique.quilis@mineco.es>

% Version 2.3 (January, 2013)

% -----------------------------------------------------------------------
% SELECTION OF THE METHOD

% Intercept pretesting
opC = -1;

% Range for grid search: default case: search is performed
% on the range [.05 .99], 100 points grid. Default.
rl = [];

switch flax1
    case 1
        % Boot-Feibes-Lisman
        res=bfl(Y,ta,d,sc);
    case 2
        % Denton
        op1 = 1; %Additive variant
        if (d == 0); d=1; end; %New version requires d=1 or d=2.
        res=denton_uni(Y,x,ta,d,sc,op1);
    case 3
        % Fernandez
         res=fernandez(Y,x,ta,sc,opC);
     case 4
         % Chow-Lin
         res=chowlin(Y,x,ta,sc,type,opC,rl);
     case 5
         % Litterman
         res=litterman(Y,x,ta,sc,type,opC,rl);
     case 6
         % Santos-Cardoso
         res=ssc(Y,x,ta,sc,type,opC,rl);
     case 7
         % Chow-Lin, fixed innovational parameter (ip)
         res=chowlin(Y,x,ta,sc,type,opC,ip);
     case 8
         % Litterman, fixed innovational parameter (ip)
         res=litterman(Y,x,ta,sc,type,opC,ip);
     case 9
         % Santos-Cardoso, fixed innovational parameter (ip)
         res=ssc(Y,x,ta,sc,type,opC,ip);
  end

% -----------------------------------------------------------------------
% SELECTION OF OUTPUT

 switch flax2
     case 1
         % Brief output
         y = res.y;
     case 2
         % Normal output
         switch res.meth
             case {'Boot-Feibes-Lisman'}
                 y = [res.y zeros(res.sc*res.N,4)];
             case {'Denton'}
                  y = [res.y zeros(res.sc*res.N,3) res.u];
             case {'Fernandez','Chow-Lin','Litterman','Santos Silva-Cardoso'}
                 y = [res.y res.y_dt res.y_lo res.y_up res.u];
         end
     case 3
         % Detailed output
         switch res.meth
             case {'Boot-Feibes-Lisman'}
                 y = [res.y zeros(res.sc*res.N,4)];
                 tduni_print(res,file_name);
             case {'Denton'}
                  y = [res.y zeros(res.sc*res.N,3) res.u];
                  tduni_print(res,file_name);
             case {'Fernandez','Chow-Lin','Litterman','Santos Silva-Cardoso'}
                 y = [res.y res.y_dt res.y_lo res.y_up res.u];
                 td_print(res,file_name);
         end
 end


text:
function [km3] = gal2km3(gal)
% Convert volume from US liquid gallons to cubic kilometers.
% Chad Greene 2012
km3 = gal*3.785411784e-12;


text:
function bernstein_polynomial_test02 ( )

%*****************************************************************************80
%
%% BERNSTEIN_POLYNOMIAL_TEST02 tests BERNSTEIN_POLY_AB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 July 2011
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BERNSTEIN_POLYNOMIAL_TEST02\n' );
  fprintf ( 1, '  BERNSTEIN_POLY_AB evaluates Bernstein polynomials over an\n' );
  fprintf ( 1, '  arbitrary interval [A,B].\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here, we demonstrate that \n' );
  fprintf ( 1, '    BPAB(N,K,A1,B1)(X1) = BPAB(N,K,A2,B2)(X2)\n' );
  fprintf ( 1, '  provided only that\n' );
  fprintf ( 1, '    (X1-A1)/(B1-A1) = (X2-A2)/(B2-A2).\n' );

  x = 0.3;
  a = 0.0;
  b = 1.0;
  bern = bernstein_poly_ab ( n, a, b, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     N     K     A        B        X       BPAB(N,K,A,B)(X)\n' );
  fprintf ( 1, '\n' );
  for k = 0 : n
    fprintf ( 1, '  %4d  %4d  %7.4f  %7.4f  %7.4f  %14.6g\n', n, k, a, b, x, bern(k+1) );
  end

  x = 1.3;
  a = 1.0;
  b = 2.0;
  bern = bernstein_poly_ab ( n, a, b, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     N     K     A        B        X       BPAB(N,K,A,B)(X)\n' );
  fprintf ( 1, '\n' );
  for k = 0 : n
    fprintf ( 1, '  %4d  %4d  %7.4f  %7.4f  %7.4f  %14.6g\n', n, k, a, b, x, bern(k+1) );
  end

  x = 2.6;
  a = 2.0;
  b = 4.0;
  bern = bernstein_poly_ab ( n, a, b, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     N     K     A        B        X       BPAB(N,K,A,B)(X)\n' );
  fprintf ( 1, '\n' );
  for k = 0 : n
    fprintf ( 1, '  %4d  %4d  %7.4f  %7.4f  %7.4f  %14.6g\n', n, k, a, b, x, bern(k+1) );
  end

  return
end


text:
function jed = cws_to_jed_gps ( c, w, s )

%*****************************************************************************80
%
%% CWS_TO_JED_GPS converts a GPS CWS date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer C, integer W, real S,
%    the GPS cycle/week/second date.
%
%    Output, real JED, the corresponding Julian Ephemeris Date.
%
  jed_epoch = epoch_to_jed_gps ( );

  d = ( 7 * ( 1024 * c + w ) ) + s / ( 24.0 * 60.0 * 60.0 );

  jed = jed_epoch + d;

  return
end


text:
function yaw=yaw_angle(sat,type,opt,beta,mu) %#ok

yaw=yaw_nominal(beta,mu);

return

repo: GINav


text:
function res = transformLine3d(line, trans)
%TRANSFORMLINE3D Transform a 3D line with a 3D affine transform.
%
%   LINE2 = transformLine3d(LINE1, TRANS)
%
%   Example
%   P1 = [10 20 30];
%   P2 = [30 40 50];
%   L = createLine3d(P1, P2);
%   T = createRotationOx(P1, pi/6);
%   L2 = transformLine3d(L, T);
%   figure; hold on;
%   axis([0 100 0 100 0 100]); view(3);
%   drawPoint3d([P1;P2]);
%   drawLine3d(L, 'b');
%   drawLine3d(L2, 'm');
%
%   See also:
%   lines3d, transforms3d, transformPoint3d, transformVector3d
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2008-11-25,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

res = [...
    transformPoint3d(line(:, 1:3), trans) ...   % transform origin point
    transformVector3d(line(:,4:6), trans)];     % transform direction vect.

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function bin_pos = createIndexVector(keypts, match_binsize, x_bin_num, y_bin_num)
%CREATEINDEXVECTOR Find the bin position to which each keypoint belongs to
%
% INPUT:
%   - keypts(1, N): input feature keypoints
%   - match_binsize: matching bin width/height
%   - x_bin_num: number of bins along x-direction
%   - y_bin_num: number of bins along y-direction
%
% OUTPUT:
%   - bin_pos(x_bin_num, y_bin_num, 4): cell containing indices of keypoint

% allocate memory
num = size(keypts, 2);
% number of classes
classes_num = 4;

% create a array for bin positions/class
bin_pos = cell(x_bin_num, y_bin_num, classes_num);

% iterate over all keypoints
for i = 1:num
    % coordinate of keypoint
    x = keypts(i).location(1);
    y = keypts(i).location(2);
    c = keypts(i).class;
    % find bin position
    bin_x = min(ceil(x/match_binsize), x_bin_num);
    bin_y = min(ceil(y/match_binsize), y_bin_num);
    % add keypoint index to corresponding bin_pos
    bin_pos{bin_x, bin_y, c} = horzcat(bin_pos{bin_x, bin_y, c}, i);
end

end

repo: Stereo-Odometry-SOFT


function [status] = execute_hash_stage(status, ipp)
    % test stage: check recorded hash checksum
    calculated = '';
    expected = '';
    try
        expected = getReferenceHash(status, ipp);
        calculated = calculateMD5Hash(status.tikzStage.texFile);

        % do the actual check
        if ~strcmpi(expected, calculated)
            % throw an error to signal the testing framework
            error('testMatlab2tikz:HashMismatch', ...
                'The hash "%s" does not match the reference hash "%s"', ...
                calculated, expected);
        end
    catch %#ok
        e = lasterror('reset'); %#ok
        [status.hashStage] = errorHandler(e);
    end
    status.hashStage.expected = expected;
    status.hashStage.found    = calculated;
end
% ==============================================================================
function hash = getReferenceHash(status, ipp)
    % retrieves a reference hash from a hash table
    % WARNING: do not make `hashTable` persistent, since this is slower

    hashTable = loadHashTable(ipp.Results.testsuite);
    if isfield(hashTable.contents, status.function)
        hash = hashTable.contents.(status.function);
    else
        hash = '';
    end
end


text:
% ######################################################################################################################################################################################
% We are here talking about spatio-temporal detections, i.e. a set of ground-truth bounding boxes that
%  I will denote by g_t, with t between t_g^b and t_g^e (beginning and end time of the ground-truth)
% versus a detection which is also a set of bounding boxes, denoted by d_t, with t between t_d^e et t_d^e.
%
% a) temporal iou =  T_i / T_u
%  this is the intersection over union between the timing of the the tubes,
% ie mathematically T_i / T_u with
% the intersection T_i = max(0,   max(t_g^b,t_d^b)-min(t_d^e,t_g^e) )
% and the union T_u = min(t_g^b,t_d^b)-max(t_d^e,t_g^e)
%
% b) for each t between max(tgb,tdb)-min(tde,tge), we compute the IoU between g_t and d_t, and average them
%
% Multiplying (a) and (b) is the same as computed the average of the spatial iou over all frames in T_u of the two tubes, with a spatial iou of 0 for frames where only one box exists.
% c) as this is standard in detection problem, if there are multiple detections for the same groundtruth detection, the first one is counted as positive and the other ones as negatives
% ######################################################################################################################################################################################
%{
gt_fnr = 1xn doube
gt_bb = nx4 doubld - [x y w h]
dt_fnr = 1xm double
dt_bb = mx4 double - [x y w h]
%}
% -------------------------------------------------------------------------
function st_iou = compute_spatio_temporal_iou(gt_fnr, gt_bb, dt_fnr, dt_bb)
% -------------------------------------------------------------------------

% time gt begin
tgb = gt_fnr(1);
% time gt end
tge = gt_fnr(end);
%time dt begin
tdb = dt_fnr(1);
tde = dt_fnr(end);
% temporal intersection
T_i = double(max(0, min(tge,tde)-max(tgb,tdb)));

if T_i>0
    T_i = T_i +1;
    % temporal union
    T_u = double(max(tge,tde) - min(tgb,tdb)+1);
    %temporal IoU
    T_iou = T_i/T_u;
    % intersect frame numbers
    int_fnr = max(tgb,tdb):min(tge,tde);

    % find the ind of the intersected frames in the detected frames
    [~,int_find_dt] = ismember(int_fnr, dt_fnr);
    [~,int_find_gt] = ismember(int_fnr, gt_fnr);

    assert(length(int_find_dt)==length(int_find_gt));

    iou = zeros(length(int_find_dt),1);
    for i=1:length(int_find_dt)
        if int_find_gt(i)<1
%             fprintf('error ')
            pf = pf;
        else
            pf = i;
        end

        gt_bound = gt_bb(int_find_gt(pf),:);
        dt_bound = dt_bb(int_find_dt(pf),:)+1;

        % gt_bound = [gt_bound(:,1:2) gt_bound(:,3:4)-gt_bound(:,1:2)];
        % dt_bound = [dt_bound(:,1:2) dt_bound(:,3:4)-dt_bound(:,1:2)];
        iou(i) = inters_union(double(gt_bound),double(dt_bound));
    end
    % finalspatio-temporal IoU threshold
    st_iou = T_iou*mean(iou);
else
    st_iou =0;
end
% % iou_thresh = 0.2,...,0.6 % 'Learing to track paper' takes 0.2 for UCF101 and 0.5 for JHMDB
% if delta >= iou_thresh
%     % consider this tube as valid detection
% end

end

% -------------------------------------------------------------------------
function iou = inters_union(bounds1,bounds2)
% -------------------------------------------------------------------------

inters = rectint(bounds1,bounds2);
ar1 = bounds1(:,3).*bounds1(:,4);
ar2 = bounds2(:,3).*bounds2(:,4);
union = bsxfun(@plus,ar1,ar2')-inters;

iou = inters./(union+eps);

end

author: gurkirt

repo: realtime-action-detection