text:
% --------------------------------------------
% demo_rand: random low-rank and sparse
% --------------------------------------------
% This demo tests the matlab code lmafit_sms.m
% for solving sparse matrix separation problem.
% July 25, 2010. Yin Zhang

function demo_rand(Solver,m,n,r,impulse,tol)

if nargin < 1; Solver = 1; end
if nargin < 2; m = 400; end
if nargin < 3; n = m; end
if nargin < 4; r = round(.1*min(m,n)); end
if nargin < 5; impulse = 0.1; end
if nargin < 6; tol = 1e-4; end

if ~exist('PROPACK','dir');
   path(path,genpath(pwd));
   warning off all;
end

sigma = 0e-3; % noise STD
tol = max(tol, .1*sigma);
tol = min(1e-1,tol);
maxit = 200;

% lmafit parameter setup
opts.tol = tol;
opts.maxit = maxit;
opts.est_rank = 1;
opts.print = 0;
% initial rank estimate
k = round(0.25*min(m,n));

% random low-rank and sparse
XY0 = randn(m,r)*randn(r,n);            % random low-rank
Z0 = rand(m,n) < impulse;               % random sparse pattern
Z0 = Z0.*randn(m,n);                    % Gaussian values
%disp(max(abs(XY0(:)))/max(abs(Z0(:))))

% data matrix D (with possible noise)
D = XY0 + Z0;
D = D + sigma*randn(m,n);

beta = []; %beta = 1.2.^(0:20);
t = 1;

fprintf('\n[m n r k] = [%i, %i, %i, %i]\n',m,n,r,k);
fprintf('noise sigma: %6.2e\n',sigma);
fprintf('---------------------------------\n')

% call solvers
for sol = Solver
    tic
    switch sol
        case 1; fprintf('   ... LMafit_sms ...\n')
            [X,Y,Z,out] = lmafit_sms_v1(D,k,opts,beta);
            XY = X*Y; iter = out.iter; rk = size(X,2);
        case 2; fprintf('   ... IALM_rpca  ...\n')
            lambda = t/sqrt(m);
            fprintf('lambda = %.2f/sqrt(m)\n',t);
            [XY,Z,iter] = inexact_alm_rpca(D,lambda,tol/10,maxit);
            rk = rank(XY);
    end
    toc
    % error calculation function
    ferr = @(U,V) norm(U-V,'fro')/norm(V,'fro');
    fprintf('XY RelErr = %8.3e\n',ferr(XY,XY0));
    fprintf('Z  RelErr = %8.3e\n',ferr( Z,Z0 ));
    fprintf('A  RelErr = %8.3e\n',ferr(XY+Z,D));
    fprintf('iter %i: final rank = %i\n',iter,rk);
    fprintf('---------------------------------\n')
end
fprintf('\n\n')

author: andrewssobral

repo: lrslibrary


text:
%TRINTERP Interpolate homogeneous transformations
%
% T = TRINTERP(T0, T1, S) is a homogeneous transform (4x4) interpolated
% between T0 when S=0 and T1 when S=1.  T0 and T1 are both homogeneous
% transforms (4x4).  Rotation is interpolated using quaternion spherical
% linear interpolation (slerp).  If S (Nx1) then T (4x4xN) is a sequence of
% homogeneous transforms corresponding to the interpolation values in S.
%
% T = TRINTERP(T1, S) as above but interpolated between the identity matrix
% when S=0 to T1 when S=1.
%
% See also CTRAJ, QUATERNION.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function T = trinterp(A, B, C)

    if nargin == 3
        %	TR = TRINTERP(T0, T1, r)
        T0 = A; T1 = B; r = C;

        if length(r) > 1
            T = [];
            for rr=r(:)'
                TT = trinterp(T0, T1, rr);
                T = cat(3, T, TT);
            end
            return;
        end

        q0 = Quaternion(T0);
        q1 = Quaternion(T1);

        p0 = transl(T0);
        p1 = transl(T1);

        qr = q0.interp(q1, r);
        pr = p0*(1-r) + r*p1;
    elseif nargin == 2
    %	TR = TRINTERP(T, r)
        T0 = A; r = B;

        if length(r) > 1
            T = [];
            for rr=r(:)'
                TT = trinterp(T0, rr);
                T = cat(3, T, TT);
            end
            return;
        end

        q0 = Quaternion(T0);
        p0 = transl(T0);

        qr = q0.scale(r);
        pr = r*p0;
    else
        error('must be 2 or 3 arguments');
    end
    T = rt2tr(qr.R, pr);

author: Allopart

repo: rbpf-gmapping


text:
function jed = ymdf_to_jed_julian ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_JED_JULIAN converts a Julian YMDF date to a JED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 December 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Edward Richards,
%    Algorithm E,
%    Mapping Time, The Calendar and Its History,
%    Oxford, 1999, pages 323-324.
%
%  Parameters:
%
%    Input, integer Y, integer M, integer D, real F, the YMDF date.
%
%    Output, real JED, the Julian Ephemeris Date.
%

%
%  Check the date.
%
  [ y, m, d, f, ierror ] = ymdf_check_julian ( y, m, d, f );

  if ( ierror ~= 0 )
    jed = -1.0;
    return
  end
%
%  Account for the missing year 0 by moving negative years up one.
%
  y2 = y_common_to_astronomical ( y );
%
%  Convert the calendar date to a computational date.
%
  y_prime = y2 + 4716 - floor ( ( 14 - m ) / 12 );
  m_prime = mod ( m + 9, 12 );
  d_prime = d - 1;
%
%  Convert the computational date to a JED.
%
  j1 = floor ( ( 1461 * y_prime ) / 4 );

  j2 = floor ( ( 153 * m_prime + 2 ) / 5 );

  jed = j1 + j2 + d_prime - 1401 - 0.5;
  jed = jed + f;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [estimate] = ft_inverse_eloreta(sourcemodel, sens, headmodel, dat, C, varargin)

% FT_INVERSE_ELORETA estimates the source activity using eLORETA
%
% Use as
%   [estimate] = ft_inverse_eloreta(sourcemodel, sens, headmodel, dat, cov, ...)
% where
%   sourcemodel is the input source model, see FT_PREPARE_SOURCEMODEL
%   sens        is the gradiometer or electrode definition, see FT_DATATYPE_SENS
%   headmodel   is the volume conductor definition, see FT_PREPARE_HEADMODEL
%   dat         is the data matrix with the ERP or ERF
%   cov         is the data covariance or cross-spectral density matrix
% and
%   estimate    contains the estimated source parameters
%
% Additional input arguments should be specified as key-value pairs and can include
%   'keepfilter'       = remember the spatial filter,    can be 'yes' or 'no'
%   'keepleadfield'    = remember the forward computation,  can be 'yes' or 'no'
%   'keepmom'          = remember the dipole moment,        can be 'yes' or 'no'
%   'lambda'           = scalar, regularisation parameter (default = 0.05)
%
% These options influence the forward computation of the leadfield
%   'reducerank'      = 'no' or number  (default = 3 for EEG, 2 for MEG)
%   'backproject'     = 'yes' or 'no', in the case of a rank reduction this parameter determines whether the result will be backprojected onto the original subspace (default = 'yes')
%   'normalize'       = 'no', 'yes' or 'column' (default = 'no')
%   'normalizeparam'  = parameter for depth normalization (default = 0.5)
%   'weight'          = number or Nx1 vector, weight for each dipole position to compensate for the size of the corresponding patch (default = 1)
%
% If the dipole definition only specifies the dipole location, a rotating dipole
% (regional source) is assumed on each location. If a dipole moment is specified, its
% orientation will be used and only the strength will be fitted to the data.
%
% This implements:
% - R.D. Pascual-Marqui; Discrete, 3D distributed, linear imaging methods of electric
%   neuronal activity. Part 1: exact, zero error localization. arXiv:0710.3341
%   2007-October-17, http://arxiv.org/pdf/0710.3341
%
% See also FT_SOURCEANALYSIS, FT_PREPARE_HEADMODEL, FT_PREPARE_SOURCEMODEL

% Copyright (C) 2013, Marlene Boenstrup, Jan-Mathijs Schoffelen and Guido Nolte
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if mod(nargin-5,2)
  % the first 5 arguments are fixed, the other arguments should come in pairs
  ft_error('invalid number of optional arguments');
end

% get the optional input arguments, or use defaults
keepfilter      = ft_getopt(varargin, 'keepfilter', 'no');
keepmom         = ft_getopt(varargin, 'keepmom', 'yes');
keepleadfield   = ft_getopt(varargin, 'keepleadfield', 'no');
lambda          = ft_getopt(varargin, 'lambda', 0.05);

% construct the low-level options for the leadfield computation as key-value pairs, these are passed to FT_COMPUTE_LEADFIELD
leadfieldopt = {};
leadfieldopt = ft_setopt(leadfieldopt, 'reducerank',     ft_getopt(varargin, 'reducerank'));
leadfieldopt = ft_setopt(leadfieldopt, 'backproject',    ft_getopt(varargin, 'backproject'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalize',      ft_getopt(varargin, 'normalize'));
leadfieldopt = ft_setopt(leadfieldopt, 'normalizeparam', ft_getopt(varargin, 'normalizeparam'));
leadfieldopt = ft_setopt(leadfieldopt, 'weight',         ft_getopt(varargin, 'weight'));

% convert the yes/no arguments to the corresponding logical values
keepfilter     = istrue(keepfilter);
keepmom        = istrue(keepmom);
keepleadfield  = istrue(keepleadfield);

% flags to avoid calling isfield repeatedly in the loop over grid positions (saves a lot of time)
hasmom        = isfield(sourcemodel, 'mom');
hasleadfield  = isfield(sourcemodel, 'leadfield');
hasfilter     = isfield(sourcemodel, 'filter');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the dipole positions that are inside/outside the brain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isfield(sourcemodel, 'inside')
  if hasfilter
    sourcemodel.inside = ~cellfun(@isempty, sourcemodel.filter);
  elseif hasleadfield
    sourcemodel.inside = ~cellfun(@isempty, sourcemodel.leadfield);
  else
    sourcemodel.inside = ft_inside_headmodel(sourcemodel.pos, headmodel);
  end
end

% convert to logical representation
sourcemodel = fixinside(sourcemodel);

% keep the original details on inside and outside positions
originside = sourcemodel.inside;
origpos    = sourcemodel.pos;

% select only the dipole positions inside the brain for scanning
sourcemodel.pos    = sourcemodel.pos(originside,:);
sourcemodel.inside = true(size(sourcemodel.pos,1),1);

if hasmom
  sourcemodel.mom = sourcemodel.mom(:,originside);
end

if hasfilter
  ft_info('using precomputed filters\n');
  sourcemodel.filter = sourcemodel.filter(originside);
elseif hasleadfield
  ft_info('using precomputed leadfields\n');
  sourcemodel.leadfield = sourcemodel.leadfield(originside);
else
  ft_info('computing forward model on the fly\n');
  if hasmom
    for i=size(sourcemodel.pos,1)
      % compute the leadfield for a fixed dipole orientation
      sourcemodel.leadfield{i} = ft_compute_leadfield(sourcemodel.pos(i,:), sens, headmodel, leadfieldopt{:}) * sourcemodel.mom(:,i);
    end
  else
    for i=1:size(sourcemodel.pos,1)
      % compute the leadfield
      sourcemodel.leadfield{i} = ft_compute_leadfield(sourcemodel.pos(i,:), sens, headmodel, leadfieldopt{:});
    end
  end
end

% use existing filters, or compute them
if ~hasfilter
  % deal with reduced rank
  % check the rank of the leadfields, and project onto the lower dimensional
  % subspace if the number of columns per leadfield > rank: this to avoid
  % numerical issues in the filter computation
  rank_lf = zeros(1,size(sourcemodel.pos,1));
  for i=1:size(sourcemodel.pos,1)
    rank_lf(i) = rank(sourcemodel.leadfield{i});
  end
  if ~all(rank_lf==rank_lf(1))
    ft_error('the forward solutions have a different rank for each location, which is not supported');
  end
  if rank_lf(1)<size(sourcemodel.leadfield{1})
    ft_notice('the forward solutions have a rank of %d, but %d orientations\n',rank_lf(1),size(sourcemodel.leadfield{1},2));
    ft_notice('projecting the forward solutions on the lower dimensional subspace\n');
    for i=1:size(sourcemodel.pos,1)
      [u,s,v{i}] = svd(sourcemodel.leadfield{i}, 'econ');
      sourcemodel.leadfield{i} = sourcemodel.leadfield{i}*v{i}(:,1:rank_lf(i));
    end
  end

  % convert the leadfield into Nchan*Ndip*Nori
  [Nchan, Nori] = size(sourcemodel.leadfield{1});
  Ndip          = numel(sourcemodel.leadfield);
  leadfield     = permute(reshape(cat(2,sourcemodel.leadfield{:}),Nchan,Nori,Ndip),[1 3 2]);

  filt = mkfilt_eloreta(leadfield, lambda);
  for i=1:size(sourcemodel.pos,1)
    sourcemodel.filter{i,1} = squeeze(filt(:,i,:))';
  end
end

% get the power
siz_C  = [size(C) 1 1]; % C can have both a freq and time dimension
sourcemodel.pow = zeros([size(sourcemodel.pos,1),siz_C(3:4)]);
sourcemodel.ori = cell(size(sourcemodel.pos,1),1);
for i=1:size(sourcemodel.pos,1)
  sourcemodel.ori{i} = zeros([size(sourcemodel.filter{i},1) siz_C(3:4)]);
  for j=1:siz_C(3)
    for k=1:siz_C(4)
      csd               = sourcemodel.filter{i}*C(:,:,j,k)*sourcemodel.filter{i}';
      [u,s,v]           = svd(real(csd));
      sourcemodel.pow(i,j,k)    = s(1);
      sourcemodel.ori{i}(:,j,k) = u(:,1);
    end
  end
end

% get the dipole moment
if keepmom && ~isempty(dat)
  siz = [size(dat) 1];
  % remove the dipole moment from the input
  if hasmom
    sourcemodel = rmfield(sourcemodel, 'mom');
  end
  for i=1:size(sourcemodel.pos,1)
    sourcemodel.mom{i} = reshape(sourcemodel.filter{i}*dat(:,:), [size(sourcemodel.filter{i},1) siz(2:end)]);
  end
end

% reassign the estimated values over the inside and outside grid positions
estimate.inside  = originside;
estimate.pos     = origpos;
if isfield(sourcemodel, 'pow')
  estimate.pow( originside,:,:) = sourcemodel.pow;
  estimate.pow(~originside,:,:) = nan;
end
if isfield(sourcemodel, 'ori') % here ori is cell
  estimate.ori( originside) = sourcemodel.ori;
  estimate.ori(~originside) = {[]};
end
if isfield(sourcemodel, 'leadfield') && keepleadfield
  estimate.leadfield( originside) = sourcemodel.leadfield;
  estimate.leadfield(~originside) = {[]};
end
if isfield(sourcemodel, 'filter') && keepfilter
  estimate.filter( originside) = sourcemodel.filter;
  estimate.filter(~originside) = {[]};
end
if isfield(sourcemodel, 'mom') && keepmom
  estimate.mom( originside) = sourcemodel.mom;
  estimate.mom(~originside) = {[]};
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function perm=field3d(Nx,Ny,Nz,k_avg,V_dp,clx,cly,clz)
% Written by Ali A. Eftekhari
% See the license file
% Does not work! Needs lots of improvements
Lx=1.0; % domain length
x=linspace(-Lx/2.0,Lx/2.0,Nx);
y=linspace(-Lx/2.0,Lx/2.0,Ny);
z=linspace(-Lx/2.0,Lx/2.0,Nz);
[X,Y,Z]=ndgrid(x,y,z);
s=-log(1-V_dp); % standard deviation
mu=log(k_avg)-s*s/2.0; % mean for a log-random field
Z=s*randn(Nx,Ny,Nz); % normal distribution
F = exp(-(X.^2/(clx*clx/2.0)+Y.^2/(cly*cly/2.0)+Z.^2/(clz*clz/2.0)));
% Gaussian filter
f =2.0/sqrt(pi)*Lx/(Nx*Ny*Nz)^(1/3)/(clx*cly*clz)^(1/3).*ifftn(fft2(Z).*fftn(F)); % another filter
perm=exp(mu+real(f)); % perm field

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
%% Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion
%% This MATLAB code was written based on the C++ code at link:
%% http://mathlib.zfx.info/html_eng/SHRotate_8h-source.html
%% which I believe was again based on the implementation from Don Williamson.
%% http://www.donw.co.uk/home/downloads/SHRotation.zip
%% The original algorithm was given in the following references:
%% [1] Joseph Ivanic and Klaus Ruedenberg
%%     J. Phys. Chem. 1996, 100, 6342-5347
%%
%% [2] Additions and Corrections (to previous paper)
%%     Joseph Ivanic and Klaus Ruedenberg
%%     J. Phys. Chem. A, 1998, Vol. 102, No. 45, 9099
%% Todo:
%% http://www.cgg.cvut.cz/~xkrivanj/papers/2005-rapport1728/rr1728-irisa-kr
%% ivanek-shrotation.pdf

function [R] = SHRotate(r, degree)

%%=============================================================
%% Project:   Spherical Harmonics
%% Module:    $RCSfile: SHRotate.m,v $
%% Language:  MATLAB
%% Author:    $Author: bjian $
%% Date:      $Date: 2007/12/27 06:23:35 $
%% Version:   $Revision: 1.4 $
%%=============================================================

% The first 1x1 sub-matrix is always 1
R{1} = 1;

% The l=1 band is rotated by a permutation of the original matrix
R1(-1+2,-1+2) = r(2,2);
R1(-1+2, 0+2) = r(3,2);
R1(-1+2, 1+2) = r(1,2);

R1( 0+2,-1+2) = r(2,3);
R1( 0+2, 0+2) = r(3,3);
R1( 0+2, 1+2) = r(1,3);

R1( 1+2,-1+2) = r(2,1);
R1( 1+2, 0+2) = r(3,1);
R1( 1+2, 1+2) = r(1,1);

R{1+1} = R1;

% Calculate each block of the rotation matrix for each subsequent band
for band=2:degree
    for m=-band:band
        for n=-band:band
            R{band+1}(m+band+1,n+band+1) = M(band,m,n,R);
        end
    end
end

function [ret] = M(l, m, n, R)

d = (m==0);
if (abs(n)==l)
    denom = (2*l)*(2*l-1);
else
    denom = (l*l-n*n);
end

u = sqrt((l*l-m*m)/denom);

v = sqrt((1+d)*(l+abs(m)-1)*(l+abs(m))/denom)*(1-2*d)*0.5;
w = sqrt((l-abs(m)-1)*(l-abs(m))/denom)*(1-d)*(-0.5);

if (u~=0)
    u = u*U(l,m,n,R);
end

if (v~=0)
    v = v*V(l,m,n,R);
end

if (w~=0)
    w = w*W(l,m,n,R);
end

ret = u+v+w;

function [ret] = U(l,m,n,R)

ret = P(0,l,m,n,R);

function [ret] = P(i,l,a,b,R)

ri1 = R{2}(i+2,1+2);
rim1 = R{2}(i+2,-1+2);
ri0 = R{2}(i+2,0+2);

if (b==-l)
    ret = ri1*R{l}(a+l,1) + rim1*R{l}(a+l, 2*l-1);
else
    if (b==l)
        ret = ri1*R{l}(a+l,2*l-1) - rim1*R{l}(a+l, 1);
    else
        ret = ri0*R{l}(a+l,b+l);
    end
end

function [ret] = V(l,m,n,R)
if (m==0)
    p0 = P(1,l,1,n,R);
    p1 = P(-1,l,-1,n,R);
    ret = p0+p1;
else
    if (m>0)
        d = (m==1);
        p0 = P(1,l,m-1,n,R);
        p1 = P(-1,l,-m+1,n,R);
        ret = p0*sqrt(1+d) - p1*(1-d);
    else
        d = (m==-1);
        p0 = P(1,l,m+1,n,R);
        p1 = P(-1,l,-m-1,n,R);
        ret = p0*(1-d) + p1*sqrt(1+d);
    end
end

function [ret] = W(l,m,n,R)
if (m==0)
    error('never gets called')
else
    if (m>0)
        p0 = P(1,l,m+1,n,R);
        p1 = P(-1,l,-m-1,n,R);
        ret = p0 + p1;
    else
        p0 = P(1,l,m-1,n,R);
        p1 = P(-1,l,-m+1,n,R);
        ret = p0 - p1;
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function [winscores]=mvg_windowbec_fast(img, windows)

if size(img,3)==3
  if isinteger(img)
    grayimage=double(rgb2gray(img))/255.0;
    img=double(img)/255.0;
  else
    grayimage=0.2989*img(:,:,1)+0.5879*img(:,:,2)+0.1140*img(:,:,3);
  end
else
  grayimage=img;
end

gmagth=0.02;%0.025;
nsubbox=3;
nbin=4; % 4, 6 or 8
binids=1:nbin;
binbounds=[1:2:2*nbin]/(2*nbin)*180;
binbounds=binbounds/180*pi;
bincenters=[0:2:(2*nbin-1)]/(2*nbin)*180;
bincenters=bincenters/180*pi;
binvects=[cos(bincenters);sin(bincenters)];

nr=size(grayimage,1);
nc=size(grayimage,2);
%edgeim=edge(grayimage,'canny');
[X,Y]=meshgrid(1:nc,1:nr);

[gv,Ngv,gvx,gvy]=gaussianderiv1D_(1);
Gx=zeros(nr,nc);
Gy=zeros(nr,nc);
if size(img,3)>1
  for cdim=1:size(img,3)
    Gxc=conv2(gv,gvx,img(:,:,cdim),'same');
    Gindicator=logical(abs(Gxc)>abs(Gx));
    Gx=Gindicator.*Gxc+~Gindicator.*Gx;
    Gyc=conv2(gvy,gv,img(:,:,cdim),'same');
    Gindicator=logical(abs(Gyc)>abs(Gy));
    Gy=Gindicator.*Gyc+~Gindicator.*Gy;
  end
else
  Gx=conv2(gv,gvx,grayimage,'same');%.*edgeim;
  Gy=conv2(gvy,gv,grayimage,'same');%.*edgeim;
end

Gx(1:Ngv,:)=0;Gx((nr-Ngv+1):nr,:)=0;
Gx(:,1:Ngv)=0;Gx(:,(nc-Ngv+1):nc)=0;
Gy(1:Ngv,:)=0;Gy((nr-Ngv+1):nr,:)=0;
Gy(:,1:Ngv)=0;Gy(:,(nc-Ngv+1):nc)=0;

Gmag=sqrt(Gx.^2+Gy.^2);
Gxu=Gx./(Gmag+eps);
Gyu=Gy./(Gmag+eps);

[cedge,cannyth]=edge(grayimage,'canny');
[cedge,cannyth]=edge(grayimage,'canny',0.5*cannyth);

cedgeint=integralimage(double(cedge));

Gvalid=conv2(gv,gv,double(cedge).*(Gmag>gmagth),'same');

Gxvalid=Gx.*Gvalid;
Gyvalid=Gy.*Gvalid;

Gs(:,:,1)=Gxvalid;
Gs(:,:,2)=Gyvalid;
%keyboard
%figure;sc(Gs,'flow');
%keyboard

oribins=zeros(nr,nc,nbin);
tmp=zeros(nr,nc);
gvs=gaussianderiv1D_(1.5);
oribinintim=zeros(nr,nc,nbin);
for i=1:nbin
  tmp=Gmag.*Gvalid.*reshape(abs([Gxu(:) Gyu(:)]*binvects(:,i)),nr,nc);
  tmp=conv2(gvs,gvs,tmp,'same');
  oribins(:,:,i)=tmp;
  %oribinintim(:,:,i)=integralimage(tmp);
  %figure;sc(tmp);
end
oribinsum=sum(oribins,3);
oribinsummax=max(oribinsum(:));
oribins=oribins/oribinsummax;
oribinsum=oribinsum/oribinsummax;
for i=1:nbin
  oribinintim(:,:,i)=integralimage(oribins(:,:,i));
  %figure;sc(tmp);
end

wincx=0.5*(windows(:,1)+windows(:,3));
wincy=0.5*(windows(:,2)+windows(:,4));
nw=size(windows,1);

%for i=1:4
  winhists=zeros(2*nsubbox,nbin,2*nsubbox);
  winhistsums=zeros(2*nsubbox,2*nsubbox);
%end

winscores=zeros(nw,1);
winareas=zeros(nw,1);
nedgepixels=zeros(nw,1);

onesv=ones(2*nsubbox+1,1);
onesh=ones(1,2*nsubbox+1);

subboxind=1:(2*nsubbox);
intersection=0;
suma=0;sumb=0;

BoxWeights=[ones(1,6); 1 0.5*ones(1,4) 1; 1 0.5 0 0 0.5 1];
BoxWeights=[BoxWeights;flipud(BoxWeights)];
BoxIndices=[2 3*ones(1,4) 4; 1 2 3 3 4 1; 1 1 2 4 1 1;1 1 4 2 1 1; ...
	    1 4 3 3 2 1; 4 3*ones(1,4) 2];

winscores=bescores(windows,oribinintim,nbin,nsubbox,BoxWeights,BoxIndices);
%fprintf('BE scores computed in C \n');
if 0
for i=1:nw
  %fprintf('%d / %d \n',i,nw);
  xa=windows(i,1);ya=windows(i,2);
  xb=windows(i,3);yb=windows(i,4);
  xaa=max(1,xa-1);yaa=max(1,ya-1);

  wi=xb-xa+1;
  hi=yb-ya+1;
  winareas(i,1)=wi*hi;
  if wi<(2*nsubbox) | hi<(2*nsubbox)
    continue;
  end

  xs=[-0.5:1/(2*nsubbox):0.5]*wi;
  ys=[-0.5:1/(2*nsubbox):0.5]*hi;
  xsr=ceil(xs+wincx(i));
  ysr=ceil(ys+wincy(i))';

  Xstart=max(1,onesv*xsr-1);
  Ystart=max(1,ysr*onesh-1);

  for j=subboxind
    for k=subboxind
      for l=binids
	winhists(j,l,k)=oribinintim(Ystart(j+1,k+1),Xstart(j+1,k+1),l)+oribinintim(Ystart(j,k),Xstart(j,k),l)-oribinintim(Ystart(j+1,k),Xstart(j+1,k),l)-oribinintim(Ystart(j,k+1),Xstart(j,k+1),l);
      end
      winhistsum(j,k)=sum(winhists(j,:,k));
    end
  end

  for j=1:(2*nsubbox)
    for k=1:(2*nsubbox)
      winscores(i,1)=winscores(i,1)+BoxWeights(j,k)*winhists(j,BoxIndices(j,k),k);
    end
  end

end
end
winscoremax=max(winscores);
winscores=(winscores./winscoremax);
%winscoremaxtmp=max(winscorestmp);
%winscorestmp=(winscorestmp./winscoremaxtmp);

%keyboard
%%%%%%%%%%%%%%%%%%%%%%%%
% Additional functions %
%%%%%%%%%%%%%%%%%%%%%%%%

function [g,N,gx,gy]=gaussianderiv1D_(sigma,N)

if nargin<2
  N=4*sigma+1;
end

% make sure that N is odd
N=2*floor(N/2)+1;
sigma2=sigma^2;

t=1:N;
mu=(N-1)/2+1;
g=1/(sqrt(2*pi)*sigma)*exp(-0.5*(t-mu).^2/sigma2);
gx=-(t-mu)/sigma2.*g;
gy=gx';

author: Cloud-CV

repo: object-proposals


text:
%TEST_TIMING_DGT_FAC  Test timing factorization DGTs
%
%   This script test the timing DGTs by comparing the results to the
%   comp_dgt_fac in the main toolbox. Therefore, the correctness of
%   comp_dgt_fac must be verified first.

routinemax=7;

Lr=[24,16,144,108,144,24,135,35,77,20];
ar=[ 4, 4,  9,  9, 12, 6,  9, 5, 7, 1];
Mr=[ 6, 8, 16, 12, 24, 8,  9, 7,11,20];

test_failed=0;

disp('--- Used subroutines ---');

for ii=1:routinemax
  which(['mex_dgt_fac_',num2str(ii)])
end;

for ii=1:length(Lr);

  L=Lr(ii);

  M=Mr(ii);
  a=ar(ii);

  b=L/M;
  N=L/a;
  c=gcd(a,M);
  d=gcd(b,N);
  p=a/c;
  q=M/c;

  for W=1:3

    for R=1:3

      for rtype=1:2
        if rtype==1
          rname='REAL ';
          f=rand(L,W);
          g=rand(L,R);
        else
          rname='CMPLX';
          f=crand(L,W);
          g=crand(L,R);
        end;

        gf=comp_wfac(g,a,M);
        cc=comp_dgt_fac(f,gf,a,M);

        for rout=1:routinemax
          cc2=feval(['mex_dgt_fac_',num2str(rout)],f,gf,a,M,0);

          res=norm(cc(:)-cc2(:));

          failed='';
          if res>10e-10
            failed='FAILED';
            test_failed=test_failed+1;
          end;

          s=sprintf('DGT  %s %i L:%3i W:%2i R:%2i a:%3i b:%3i c:%3i d:%3i p:%3i q:%3i %0.5g %s',rname,rout,L,W,R,a,b,c,d,p,q,res,failed);
          disp(s)

        end;

      end;

    end;

  end;

end;

test_failed

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function e = gsp_assert_test(X1,X2,tol, name)

if X1==0
    if norm(X2(:),'fro')<tol
        fprintf(['Test ',name, ': OK\n']);
        e = 0;
    else
        warning(['Test ',name, ': Error, badness ',...
            num2str(norm(X2(:),'fro')),...
            ', tol ',num2str(tol)]);
        e = 1;
    end
else
    if norm(X1(:)-X2(:),'fro')/norm(X2(:),'fro')<tol
        fprintf(['Test ',name, ': OK\n']);
        e = 0;
    else
        warning(['Test ',name, ': Error, badness ',...
            num2str(norm(X1(:)-X2(:),'fro')/norm(X2(:),'fro')), ...
            ', tol ',num2str(tol)]);
        e = 1;
    end
end

end

repo: gspbox


text:
clear; clc;
dataset='SearchSnippets'; %SearchSnippets, StackOverflow, Biomedical
method = 'SVM'; % SVM(Support Vector Machine)
isKernel = 0; %1: SVM Gaussian Kernel; 0: SVM Linear Kernel.
Weighting = 'TF'; %TF, TFIDF or AE(Average Embedding)
dataStr=['./../../dataset/',dataset,'-lite.mat'];
load(dataStr);
rand('state',0)
randn('state',0)
%%
disp('Step 1 get the train and test data ...')
% TFIDF
if (strcmp(Weighting,'TF'))
    testFea = fea(testIdx, :);
    trainFea = fea(trainIdx, :);
    testGnd = gnd(testIdx, :);
    trainGnd = gnd(trainIdx, :);
elseif (strcmp(Weighting,'TFIDF'))
    testFea = fea(testIdx, :);
    trainFea = fea(trainIdx, :);
    testGnd = gnd(testIdx, :);
    trainGnd = gnd(trainIdx, :);
    [trainFea, testFea] = tf_idf_tt(trainFea, testFea);
elseif (strcmp(Weighting,'AE'))
    dataStr=['./../../dataset/',dataset,'-STC2.mat'];
    load(dataStr);
    parameters.wordDim = 48;
    parameters.vocSize = size_vocab;
    disp('AE: Generate verage embedding vectors ...')
    % Step a. Generate word vector sets
    CR_E = randi([-25,25],parameters.wordDim,parameters.vocSize)/100;
    disp(strcat('Number of weights E:',num2str(size(CR_E))));
    vocab_emb_length = length(vocab_emb_Word2vec_48(1,:));
    if vocab_emb_length > size_vocab
        error(['Error, and the size fo vocab_emb is:',vocab_emb_length])
    end
    CR_E(1:parameters.wordDim,vocab_emb_Word2vec_48_index) = vocab_emb_Word2vec_48(1:parameters.wordDim,1:vocab_emb_length);
    % Step b. Average Embedding
    textSize = length(fea_All(:,1));
    fea =[];
    for i=1:textSize
        tmp_fea_vector_weight = repmat(fea_All(i,find(fea_All(i,:)>0)),parameters.wordDim,1);
        tmp_fea_vector_matrix = CR_E(:,find(fea_All(i,:)>0)) .* tmp_fea_vector_weight;
        tmp_fea_vector = sum(tmp_fea_vector_matrix,2);
        fea(i,:) = tmp_fea_vector';
        if mod(i,2000) == 0
            disp(['has averaged embedding number:',num2str(i)]);
        end
    end
    testFea = fea(testIdx, :);
    trainFea = fea(trainIdx, :);
    testGnd = gnd(testIdx, :);
    trainGnd = gnd(trainIdx, :);
    %
end
testFea = normalize(testFea);
trainFea = normalize(trainFea);
%%
disp('step 2 train model ...')
if (strcmp(method,'SVM'))
    trainFea = sparse(trainFea);
    testFea = sparse(testFea);
    if ~isKernel
        disp('start train linear SVM model ...')
        model = train(trainGnd, trainFea, '-q');
        disp('start predict test data via linear SVM ...')
        disp('step 3 predict test data ...')
        [predict_label, accuracy, predict_scores] = predict(testGnd, testFea, model, '-b 1');
    else
        disp('start train kernel SVM model ...')
        model = svmtrain(trainGnd, trainFea, '-t 0');
        disp('start predict test data via kernel SVM model ...')
        disp('step 3 predict test data ...')
        [predict_label, accuracy, predict_scores] = svmpredict(testGnd, testFea, model);
    end
end
AC = length(find(predict_label == testGnd))/length(testGnd)*100;
disp(['Accuracy is ',num2str(AC)])

author: jacoxu


text:
function [X_poly] = polyFeatures(X, p)
%POLYFEATURES Maps X (1D vector) into the p-th power
%   [X_poly] = POLYFEATURES(X, p) takes a data matrix X (size m x 1) and
%   maps each example into its polynomial features where
%   X_poly(i, :) = [X(i) X(i).^2 X(i).^3 ...  X(i).^p];
%

% You need to return the following variables correctly.
X_poly = zeros(numel(X), p);

% ====================== YOUR CODE HERE ======================
% Instructions: Given a vector X, return a matrix X_poly where the p-th
%               column of X contains the values of X to the p-th power.
%
%

X_poly(:,1) = X;
for i=2:p
  X_poly(:,i) = X.*X_poly(:,i-1);
end

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

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function data_num = p03_data_num ( )

%*****************************************************************************80
%
%% P03_DATA_NUM returns the number of data points for problem p03.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DATA_NUM, the number of data points.
%
  data_num = 11;

  return
end


text:
%%***********************************************************************
%% gencorrmat2: generate correlation matrix
%%
%%  B = gencorrmat(n,options);
%%
%%  options = 'Randcorr'
%%          = 'Randcorr2'
%%          = 'Diag'
%%          = 'AR1'
%%          = 'CompSym'
%%          = 'Rand'
%%***********************************************************************
%% SDPNAL+
%% Copyright (c) 2014 by
%% Liuqin Yang, Defeng Sun, and Kim-Chuan Toh
%%***********************************************************************

  function B = gencorrmat(n,options,alpha)

   rng('default')
   if nargin == 2; alpha = 0; end

   if strcmp(options,'Randcorr') || (options == 1)
      beta = 10^(-4/(n-1));
      xx = beta.^[0:n-1];
      xx = n*xx/sum(xx);
      B = gallery('randcorr',xx);
   elseif strcmp(options,'Randcorr2') || (options == 2)
      n2 = n/2;
      beta = 10^(-4/(n-1));
      xx = [beta.^[0:n2-1], zeros(1,n2)];
      xx = n*xx/sum(xx);
      B = gallery('randcorr',xx);
   elseif strcmp(options,'Diag') || (options == 3)
      n2 = n/2;
      tmp = [ones(n2), zeros(n2); zeros(n2), eye(n2)];
      B = diag(1e4*(2*rand(n,1)-1));
   elseif strcmp(options,'AR1') || (options == 4)
      tmp = -diag(ones(n,1)) - 0.8*diag(ones(n-1,1),-1);
      Dinv = (1/0.01)*speye(n,n);
      Binv = tmp'*Dinv*tmp;
      B = inv(Binv); B = 0.5*(B+B');
      diagB = diag(B);
      B = diag(1./sqrt(diagB))*B*diag(1./sqrt(diagB));
   elseif strcmp(options,'CompSym') || (options == 5)
      sigma = 1; rho = 0.5;
      dd  = sigma^2*(1-[0:n-1]*(rho^2)./(1+[0:n-1]*rho));
      tmp = zeros(n,n);
      for j=2:n
         tmp(1:j-1,j) = rho./(1+[1:j-1]'*rho);
      end
      tmp = -eye(n) -tmp;
      Dinv = diag(1./dd);
      Binv = tmp'*Dinv*tmp;
      B = inv(Binv); B = 0.5*(B+B');
      diagB = diag(B);
      B = diag(1./sqrt(diagB))*B*diag(1./sqrt(diagB));
   elseif strcmp(options,'Rand') || (options == 6)
      B = 2*rand(n)-1;
      B = 0.5*(B+B');
   end
%%
%% add noise
%%
    E = 2*(rand(n)-0.5); E = triu(E,1)+ triu(E,1)';
    B = (1-alpha)*B + alpha*E;
    B = min(B,1);
    B = max(-1,B);
    for i=1:n; B(i,i) = 1; end
%%****************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
%COFICOSTFUNC Collaborative filtering cost function
%   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
%   num_features, lambda) returns the cost and gradient for the
%   collaborative filtering problem.
%

% Unfold the U and W matrices from params
X = reshape(params(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(params(num_movies*num_features+1:end), ...
                num_users, num_features);

% You need to return the following values correctly
J = 0;
X_grad = zeros(size(X));
Theta_grad = zeros(size(Theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost function and gradient for collaborative
%               filtering. Concretely, you should first implement the cost
%               function (without regularization) and make sure it is
%               matches our costs. After that, you should implement the
%               gradient and use the checkCostFunction routine to check
%               that the gradient is correct. Finally, you should implement
%               regularization.
%
% Notes: X - num_movies  x num_features matrix of movie features
%        Theta - num_users  x num_features matrix of user features
%        Y - num_movies x num_users matrix of user ratings of movies
%        R - num_movies x num_users matrix, where R(i, j) = 1 if the
%            i-th movie was rated by the j-th user
%
% You should set the following variables correctly:
%
%        X_grad - num_movies x num_features matrix, containing the
%                 partial derivatives w.r.t. to each element of X
%        Theta_grad - num_users x num_features matrix, containing the
%                     partial derivatives w.r.t. to each element of Theta
%

J = (1 / 2) * sum(sum((X * Theta' .* R - Y) .^ 2)) + ...
    (lambda / 2) * sum(sum(Theta .^ 2)) + (lambda / 2) * sum(sum(X .^ 2));

X_grad = (X * Theta' .* R - Y) * Theta + lambda * X;

Theta_grad = (X * Theta' .* R - Y)' * X + lambda * Theta;

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

grad = [X_grad(:); Theta_grad(:)];

end

repo: machine-learning


text:
function [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters)
%GRADIENTDESCENT Performs gradient descent to learn theta
%   theta = GRADIENTDESENT(X, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCost) and gradient here.
    %
    theta_temp = theta;
    for j = 1:size(X, 2)
        theta_temp(j) = theta(j)-alpha*(1/m)*(X*theta - y)' * X(:, j);
    end
    theta = theta_temp;

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

    % Save the cost J in every iteration
    J_history(iter) = computeCost(X, y, theta);

end

end

author: zlotus


text:
classdef Crop < dagnn.ElementWise
%CROP DagNN cropping layer.
%    This is a pecurial layer from FCN. It crops inputs{1} to
%    match the size of inputs{2} (starting with a base crop amount).
%    A future version

  properties
    crop = [0 0]
  end

  properties (Transient)
    inputSizes = {}
  end

  methods
    function crop = getAdaptedCrops(obj)
      cropv = obj.inputSizes{1}(1) - obj.inputSizes{2}(1) ;
      cropu = obj.inputSizes{1}(2) - obj.inputSizes{2}(2) ;
      cropv1 = max(0, cropv - obj.crop(1)) ;
      cropu1 = max(0, cropu - obj.crop(2)) ;
      crop = [cropv - cropv1, cropv1, cropu - cropu1, cropu1] ;
    end

    function outputs = forward(obj, inputs, params)
      obj.inputSizes = cellfun(@size, inputs, 'UniformOutput', false) ;
      adjCrop = obj.getAdaptedCrops() ;
      outputs{1} = vl_nncrop(inputs{1}, adjCrop) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      adjCrop = obj.getAdaptedCrops() ;
      derInputs{1} = vl_nncrop(inputs{1}, adjCrop, derOutputs{1}, obj.inputSize{1}) ;
      derInputs{2} = [] ;
      derParams = {} ;
    end

    function reset(obj)
      obj.inputSizes = {} ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      obj.inputSizes = inputSizes ;
      crop = obj.getAdaptedCrops() ;
      outputSizes{1} = inputSizes{1} - [crop(1)+crop(2), crop(3)+crop(4), 0, 0]
    end

    function rfs = getReceptiveFields(obj)
      rfs(1,1).size = [1 1] ;
      rfs(1,1).stride = [1 1] ;
      rfs(1,1).offset = 1 + obj.crop ;
      rfs(2,1).size = [] ;
      rfs(2,1).stride = [] ;
      rfs(2,1).offset = [] ;
    end

    function obj = Crop(varargin)
      obj.load(varargin) ;
    end
  end
end

author: thomas-koehler

repo: SupER


text:
function T_plotting = CalculatePlottingMatrix(T_new,x_intervals,y_intervals)

a = 1;

for x_index = 2:1:(x_intervals-1)
    b = 1;
    for y_index = 2:1:(y_intervals-1)
         T_plotting(a,b) = T_new(x_index,y_index);
         b = b+1;
    end

    a = a+1;
end

end


text:
% op_takecoils.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% out=op_takecoils(in,index);
%
% DESCRIPTION:
% Extract the rf coil channels with indices corresponding to the 'index' input
% array.
%
% INPUTS:
% in     = input data in matlab structure format.
% index  = vector indicating the indices of the rf channels you would like to extract.
%
% OUTPUTS:
% out    = Output dataset consisting of averages extracted from the input.

function out=op_takecoils(in,index);

if in.dims.coils==0
    %Can't take average because there are none:
    error('ERROR:  There aren''t multiple coils in this dataset!  Aborting!');
elseif in.dims.coils==1
    %SHOULD NEVER HAPPEN (Time dimension should always be dim=1)
    error('ERROR:  dims.coils==1.  This should never happen!  Aborting!');
elseif in.dims.coils==2
    fids=in.fids(:,index,:,:,:);
elseif in.dims.coils==3;
    fids=in.fids(:,:,index,:,:);
elseif in.dims.coils==4;
    fids=in.fids(:,:,:,index,:);
elseif in.dims.coils==5
    fids=in.fids(:,:,:,:,index);
end

fids=squeeze(fids);

%re-calculate Specs using fft
specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

%change the dims variables
dims.t=in.dims.t;
if length(index)==1
    dims.coils=0;
    if in.dims.averages > in.dims.coils
        dims.averages=in.dims.averages-1;
    else
        dims.averages=in.dims.averages;
    end
    if in.dims.subSpecs > in.dims.coils
        dims.subSpecs=in.dims.subSpecs-1;
    else
        dims.subSpecs=in.dims.subSpecs;
    end
    if in.dims.extras > in.dims.coils
        dims.extras=in.dims.extras-1
    else
        dims.extras=in.dims.extras;
    end
elseif length(index)>1
    dims.coils=in.dims.coils;
    dims.averages=in.dims.averages;
    dims.subSpecs=in.dims.subSpecs;
    dims.extras=in.dims.extras;
end

%re-calculate the sz variable
sz=size(fids);

%FILLING IN DATA STRUCTURE
out=in;
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.dims=dims;

%FILLING IN THE FLAGS
out.flags=in.flags;
out.flags.writtentostruct=1;
if length(index)<2
    out.flags.addedrcvrs=1;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [in3] = cm32in3(cm3)
% Convert volume from cubic centimeters to cubic inches.
% Chad Greene 2012
in3 = cm3*0.061023744095;


text:
function [sigma_draw, Ft_draw] = VARdrawpost(VAR)
% =======================================================================
% Draw from the posterior distribution of a VAR model
% =======================================================================
% [sigma_draw, Ft_draw] = VARdrawpost(VAR)
% -----------------------------------------------------------------------
% INPUT
%   - VAR: structure, result of VARmodel function
% -----------------------------------------------------------------------
% OUPUT
%   - sigma_draw: draw from the posterior of VCV matrix of residuals of VAR
%   - Ft_draw: draw from the posterior of Ft of VAR
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

%% Get relevant parameters from VAR structure
%===============================================
nobs = VAR.nobs;
k    = [];
nvar = VAR.nvar;
X    = VAR.X;

%% OLS estimates
%============================================
Ft_hat = VAR.Ft;
sigma_hat = VAR.sigma;
inv_sigma_hat = inv(sigma_hat);

%% Draw the VCV matrix (sigma)
%============================================
inv_sigma_draw = wishrnd(inv_sigma_hat/nobs,nobs);
sigma_draw     = inv(inv_sigma_draw);

%% Draw the coeffiecient matrix (Ft)
%============================================
aux1 = inv(X'*X);
aux2 = kron(sigma_draw, aux1);
Fthat_vec = Ft_hat(:);
Ftdraw = mvnrnd(Fthat_vec,aux2);
Ft_draw = reshape(Ftdraw,k,nvar);

author: ambropo

repo: VAR-Toolbox

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


text:
function test_suite=test_montecarlo_phase_stat
% tests for test_phase_stat
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function r=randint()
    r=ceil(rand()*10+10);

function test_phase_stat_basics
    ds=generate_random_phase_dataset(40+randint(),'small');
    nsamples=size(ds.samples,1);

    methods={'param','nonparam_nan','nonparam',''};
    outputs={'pbi','pos','pop'};

    for k=1:numel(outputs)
        method=methods{k};
        for j=1:numel(outputs)
            output=outputs{j};

            opt=struct();
            opt.niter=randint();
            opt.progress=false;

            opt.permuter_func=@(iter) deterministic_permute(nsamples,...
                                                        opt.niter,iter);
            opt.output=output;
            opt.seed=randint();

            is_parametric=false;
            extreme_tail_is_nan=true;
            switch method
                case 'param'
                    opt.zscore='parametric';
                    is_parametric=true;

                case 'nonparam_nan'
                    opt.zscore='non_parametric';

                case 'nonparam'
                    opt.extreme_tail_set_nan=false;
                    extreme_tail_is_nan=false;

                case ''
                    % deafults,ok

                otherwise
                    assert(false);
            end

            expected_samples=compute_expected_samples(ds,output,...
                                        opt.niter,opt.permuter_func,...
                                        is_parametric,...
                                        extreme_tail_is_nan);
            result=cosmo_montecarlo_phase_stat(ds,opt);

            assertElementsAlmostEqual(expected_samples,result.samples,...
                                            'absolute',1e-5);
        end
    end

function test_random_data_nonparam_uniformity
% when getting phase stats for random data, z-scores must follow some sort
% of z-like distribution
    ds=generate_random_phase_dataset(40+randint(),'big');
    nsamples=size(ds.samples,1);

    methods={'nonparam_nan','nonparam',''};
    outputs={'pbi','pos','pop'};
    for k=1:numel(outputs)
        method=methods{k};
        for j=1:numel(outputs)
            output=outputs{j};

            opt=struct();
            opt.niter=50+randint();
            opt.output=output;
            opt.seed=[];
            opt.permuter_func=@(unused)nondeterministic_permute(nsamples);
            opt.progress=false;

            if ~isempty(method)
                opt.zscore='non_parametric';
                opt.extreme_tail_set_nan=~strcmp(method,'nonparam');
            end

            stat_ds=cosmo_montecarlo_phase_stat(ds,opt);
            samples=stat_ds.samples;
            nan_msk=isnan(samples);
            assert(mean(nan_msk)<.2); % not too many nans

            z_sorted=sort(samples(~nan_msk));
            n_z=numel(z_sorted);

            p_uniform=(.5:n_z)/n_z;
            z_uniform=cosmo_norminv(p_uniform);

            r2=var(z_sorted);
            r2_resid=var(z_sorted-z_uniform);
            F=r2/r2_resid;
            assert(F>10);
        end
    end

function ds=generate_random_phase_dataset(nsamples_per_class,size_str)
    ds=cosmo_synthetic_dataset('ntargets',2,...
                                    'nchunks',nsamples_per_class,...
                                    'size',size_str,...
                                    'seed',0);
    sz=size(ds.samples);
    ds.samples=randn(sz)+1i*randn(sz);
    ds.sa.chunks(:)=1:sz(1);

function samples=compute_expected_samples(ds,output,...
                                        niter,permuter_func,...
                                        is_parametric,...
                                        extreme_tail_is_nan)

    stat_orig=cosmo_phase_stat(ds,'output',output);
    [nsamples,nfeatures]=size(ds.samples);

    stat_null_cell=cell(niter,1);
    for iter=1:niter
        rp=permuter_func(iter);
        ds_null=ds;
        ds_null.sa.targets=ds.sa.targets(rp);
        stat=cosmo_phase_stat(ds_null,'output',output);
        stat_null_cell{iter}=stat;
    end

    stat_null=cosmo_stack(stat_null_cell);

    if is_parametric
        mu=mean(stat_null.samples,1);
        sd=std(stat_null.samples,[],1);

        samples=(stat_orig.samples-mu)./sd;
    else
        count_gt=sum(bsxfun(@gt,stat_orig.samples,stat_null.samples),1);
        count_lt=sum(bsxfun(@lt,stat_orig.samples,stat_null.samples),1);

        msk_gt=count_gt>niter/2;
        msk_lt=count_lt>niter/2;

        p=zeros(1,nfeatures)+.5;
        p(msk_gt)=count_gt(msk_gt)/(1+niter);
        p(msk_lt)=1-count_lt(msk_lt)/(1+niter);

        min_p=1/(1+niter)+1e-10;

        assert(all(p>=min_p-2e-10));
        assert(all((1-p)>=(min_p-2e-10)));

        if extreme_tail_is_nan
            p(count_gt==niter | count_lt==niter)=NaN;
        end

        samples=cosmo_norminv(p);
    end

function func=get_determistic_permute_func(ntargets,niter)
    func=@(iter) deterministic_permute(ntargets,niter,iter);

function targets_idxs=deterministic_permute(ntargets,niter,iter)
    persistent cached_rand_vec;
    persistent cached_args;

    args={ntargets,niter};

    if ~isequal(args,cached_args)
        cached_rand_vec=cosmo_rand(ntargets,1,'seed',ntargets*niter);

        cached_args=args;
    end

    rand_vals=cached_rand_vec+iter/niter;
    msk=rand_vals>1;
    rand_vals(msk)=rand_vals(msk)-1;

    [unused,targets_idxs]=sort(rand_vals,1);

function target_idxs=nondeterministic_permute(ntargets)
    rand_vals=randn(ntargets,1);
    [unused,target_idxs]=sort(rand_vals);

function test_monte_carlo_phase_stat_seed
    ds=generate_random_phase_dataset(20,'tiny');

    opt=struct();
    opt.niter=10+randint();
    opt.output='pbi';
    opt.progress=false;

    % different results with empty seeed
    opt.seed=[];
    r1=cosmo_montecarlo_phase_stat(ds,opt);
    attempt=10;
    while attempt>0
        attempt=attempt-1;
        assert(attempt>0,'results are always the same');
        r2=cosmo_montecarlo_phase_stat(ds,opt);
        if ~isequal(r1.samples,r2.samples)
            break;
        end
    end

    % fixed seed, same result
    opt.seed=randint();
    r1=cosmo_montecarlo_phase_stat(ds,opt);
    r2=cosmo_montecarlo_phase_stat(ds,opt);
    assertElementsAlmostEqual(r1.samples,r2.samples);

    % different seed, different result
    attempt=10;
    while attempt>0
        opt.seed=opt.seed+1;
        attempt=attempt-1;
        assert(attempt>0,'results are always the same');
        r2=cosmo_montecarlo_phase_stat(ds,opt);
        if ~isequal(r1.samples,r2.samples)
            break;
        end
    end

function test_montecarlo_phase_stat_exceptions()
    func=@cosmo_montecarlo_phase_stat;
    aet=@(x,varargin)assertExceptionThrown(@()...
                func(x,varargin{:}),'');
    extra_args=cosmo_structjoin({'progress',false,...
                                'niter',3,...
                                'output','pbi'});
    aet_arg=@(x,varargin)aet(x,extra_args,varargin{:});

    % valid
    ds=generate_random_phase_dataset(5,'tiny');
    func(ds,extra_args); % ok

    % unbalanced targets
    bad_ds=ds;
    i=find(ds.sa.targets==2,1,'first');
    bad_ds.sa.targets(i)=1;
    aet_arg(bad_ds);

    % invalid output
    aet_arg(ds,'output','foo');

    % valid zscore
    func(ds,extra_args,'zscore','parametric');
    func(ds,extra_args,'zscore','non_parametric');

    % invalid zscore
    aet_arg(ds,'zscore','nonparametric');
    aet_arg(ds,'zscore','foo');

    % invalid niter
    aet_arg(ds,'niter',.3);
    aet_arg(ds,'niter',-3);
    aet_arg(ds,'niter',[2 2]);
    aet_arg(ds,'niter','f');

    % valid output
    func(ds,extra_args,'output','pbi');
    func(ds,extra_args,'output','pos');
    func(ds,extra_args,'output','pop');

    % invalid output
    aet_arg(ds,'output','foo');
    aet_arg(ds,'output',2);

    % missing fields
    aet(ds,rmfield(extra_args,'niter'));
    aet(ds,rmfield(extra_args,'output'));

    % valid extreme_tail_set_nan
    func(ds,extra_args,'extreme_tail_set_nan',true);
    func(ds,extra_args,'extreme_tail_set_nan',false);

    % invalid func(ds,extra_args,'extreme_tail_set_nan',true);
    aet(ds,extra_args,'extreme_tail_set_nan',2);
    aet(ds,extra_args,'extreme_tail_set_nan','foo');

function test_unit_length_exception
    ds=cosmo_synthetic_dataset('nchunks',10);
    sample_size=size(ds.samples);
    ds.sa.chunks(:)=1:sample_size(1);

    rand_func=single(randn(sample_size));

    ds.samples=rand_func() + 1i*rand_func();

    opt=struct();
    opt.output='pos';
    opt.niter=100;
    opt.progress=false;

    % should not raise an exception
    cosmo_montecarlo_phase_stat(ds,opt);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function edl_writescript(filename, guid, workdir, ctrlpath, props)
%EDL_WRITESCRIPT Writes an EDL script
%
% $ Syntax $
%   - edl_writescript(filename, guid, workdir, ctrlpath, props)
%
% $ Arguments $
%   - filename:     the filename of the destination script
%   - guid:         the guid string assigned to the script
%   - workdir:      the root working directory of the experiments
%   - ctrlpath:     the path of control file (r.t. script's parent)
%   - props:        the properties of the experiment parameters
%                   (no need of internal_index)
%
% $ Description $
%   - edl_writescript(filename, guid, workdir, ctrlpath, props) writes the
%     script according to the information provided. It will also add the
%     internal index to each entry.
%
% $ Remarks $
%   - The function will also creates the initial control file.
%
% $ History $
%   - Created by Dahua Lin, on Aug 14, 2006
%

%% parse and verify input

if nargin < 5
    raise_lackinput('edl_writescript', 5);
end

%% Prepare elements

doctag = 'ExpScript';
nodetag = 'Entry';

attribs.guid = guid;
attribs.workdir = workdir;
attribs.ctrlpath = ctrlpath;

n = length(props);
for i = 1 : n
    props(i).internal_index = i;
end

%% Write

edl_writeprops(doctag, attribs, nodetag, props, filename);
cpath = sladdpath(ctrlpath, slfilepart(filename, 'parent'));
edl_initctrlfile(cpath, guid, n);

author: lmthang


text:
function [OP,MM] = CurvedPoissonIPDG2D()

% function [OP,MM] = CurvedPoissonIPDG2D()
% Purpose: Set up a discrete Poisson matrix and mass matrix
%          using Interior Penalty Discontinuous Galerkin (IPDG).

Globals2D;

NGauss = gauss.NGauss;

% build DG derivative matrices
MM  = zeros(K*Np*Np, 3);  OP = zeros(K*Np*Np*(1+Nfaces), 3);

% global node numbering
entries = (1:Np*Np)'; entriesMM = (1:Np*Np)';
for k1=1:K
  if(~mod(k1,200)) k1, end;

  % Location of k1'th diagonal block entries in OP matrix
  rows1 = ((k1-1)*Np+1:k1*Np)'*ones(1,Np); cols1 = rows1';

  % Extract local mass matrix and cubature weights
  locmm = cub.mm(:,:,k1);
  cw  = spdiags(cub.W (:,k1), 0, cub.Ncub, cub.Ncub);

  % Evaluate derivatives of Lagrange basis functions at cubature nodes
  [cDx, cDy]  = PhysDmatrices2D(x(:,k1), y(:,k1), cub.V);

  % Evaluate local stiffness matrix
  OP11 = cDx'*cw*cDx + cDy'*cw*cDy;

  % Build element-to-element parts of stiffness matrix for element k1
  for f1=1:Nfaces

    % Find neighbor
    k2 = EToE(k1,f1); f2 = EToF(k1,f1);

    idsM = (f1-1)*NGauss+1:f1*NGauss;

    % Extract Lagrange basis function -> Gauss node interpolation matrix
    gVM = gauss.finterp(:,:,f1);
    gVP = gauss.finterp(:,:,f2);    gVP = gVP(NGauss:-1:1,:);

    % Evaluate spatial derivatives of  Lagrange basis function at Gauss nodes
    [gDxM, gDyM] = PhysDmatrices2D(x(:,k1), y(:,k1),gVM);
    [gDxP, gDyP] = PhysDmatrices2D(x(:,k2), y(:,k2),gVP);

    % Evaluate normals at Gauss nodes on face
    gnx = spdiags(gauss.nx(idsM, k1), 0, NGauss, NGauss);
    gny = spdiags(gauss.ny(idsM, k1), 0, NGauss, NGauss);
    gw  = spdiags(gauss.W(idsM, k1),  0, NGauss, NGauss);

    % Compute normal derivatives of Lagrange basis functions at Gauss nodes
    gDnM = gnx*gDxM + gny*gDyM;
    gDnP = gnx*gDxP + gny*gDyP;

    % Locate global numbers of Lagrange nodes in neighboring element
    cols2 = ones(Np,1)*((k2-1)*Np+1:k2*Np);

    % Find minimum height of two elements sharing this face
    hinv = max(Fscale( 1 + (f1-1)*Nfp, k1), Fscale( 1 + (f2-1)*Nfp, k2));

    % Set penalty scaling
    gtau = 20*(N+1)*(N+1)*hinv;

    % Determine type of face
    switch(BCType(k1,f1))
      case {Dirichlet}
        % Dirichlet boundary face variational terms
	OP11 = OP11 + ( gVM'*gw*gtau*gVM - gVM'*gw*gDnM - gDnM'*gw*gVM);

      case {Neuman}
        % Do nothing
      otherwise
	% Interior face variational terms for stiffness matrix
	OP11 = OP11 + 0.5*( gVM'*gw*gtau*gVM - gVM'*gw*gDnM - gDnM'*gw*gVM );

	OP12 =      - 0.5*( gVM'*gw*gtau*gVP + gVM'*gw*gDnP - gDnM'*gw*gVP );

        % Store self-neighbor interaction term in global stiffness matrix
        OP(entries(:), :) = [rows1(:), cols2(:), OP12(:)];
	entries = entries + Np*Np;
    end
  end

  % Store k1'th self-self interaction term in global stiffness matrix
  OP(entries(:), :)   = [rows1(:), cols1(:), OP11(:)];
  MM(entriesMM(:), :) = [rows1(:), cols1(:), locmm(:)];
  entries = entries + Np*Np; entriesMM = entriesMM + Np*Np;
end

% Convert OP and MM from coordinate storage format to Matlab's intrinsic sparse matrix format
OP = OP(1:max(entries)-Np*Np,:);
OP = myspconvert(OP, Np*K, Np*K, 1e-15); MM = myspconvert(MM, Np*K, Np*K, 1e-15);
return

repo: Aero-matlab


text:
function [bbci, data]= bbci_calibrate_NIRS(bbci, data)
%BBCI_CALIBRATE_NIRS - Calibrate for ERP-design with NIRS
%
%This function is called by bbci_calibrate
%(if BBCI.calibate.analyze_fcn is set to @bbci_calibrate_NIRS_ERP).
%Via BBCI.calibrate.settings, the details can be specified, se below.
%
%Synopsis:
% [BBCI, DATA]= bbci_calibrate_NIRS_ERP(BBCI, DATA)
%
%Arguments:
%  BBCI -  the field 'calibrate.settings' holds parameters specific to
%          calibrate CSP-based BCI processing.
%  DATA -  holds the calibration data
%
%Output:
%  BBCI - Updated BBCI structure in which all necessary fields for
%     online operation are set, see bbci_apply_structures.
%  DATA - As input but added some information of the analysis that might
%     be reused in a second run
%
%BBCI.calibrate.settings may include the following parameters:
%  visu_clab:   Clab used for visualization
%  lp:          butterworth lowpass filter frequency
%  derivative   use derivative (slope) of signal as feature
%  classes: [1x2 CELL of CHAR] Names of the two classes, which are to be
%     discriminated. For classes = 'auto', all pairs of available classes
%     are investigated, and the one with best xvalidation performance is
%     chosen. Default is 'auto'.   **** TO BE IMPLEMENTED ****
% Display options:
%  zscore:      (for display purposes) channels are normalized to zero mean
%               and unit variance
%  grid:        provide a string specifying the grid layout for grid_plot
%
% ------------
%  ival: [1x2 DOUBLE] interval on which CSP is performed, 'auto' means
%     automatic selection. Default is 'auto'.
%  band: [1x2 DOUBLE] frequency band on which CSP is performed, 'auto' means
%     automatic selection. Default is 'auto'.
%  clab: [CELL] Labels of the channels that are used for classification,
%     default {'*'}.
%  nPatters: [INT>0] number of CSP patterns which are considered from each
%     side of the eigenvalue spectrum. Note, that not neccessarily all of
%     these are used for classification, see settings.pattern.
%     Default is 3.
%  patterns: [1xN DOUBLE or 'auto'] vector specifying the indices of the
%     CSP filters that should be used for classification; 'auto' means
%     automatic selection. Default is 'auto'.
%  model: [CHAR or CELL] Classification model.
%     Default {'RLDAshrink', 'gamma',0, store_means',1, 'scaling',1}.
%  reject_artifacts:
%  reject_channels:
%  reject_artifacts_opts: cell array of options which are passed to
%     reject_varEventsAndChannels.
%  check_ival: interval which is checked for artifacts
%  do_laplace: do Laplace spatial filtering for automatic selection
%     of ival/band. If opt.do_laplace is set to 0, the default value of
%     will also be set to opt.visu_laplace 0 (but can be overwritten).
%  visu_laplace: do Laplace filtering for grid plots of Spectra/ERDs.
%     If visu_laplace is a string, it is passed to proc_laplace. This
%     can be used to use alternative geometries, like 'vertical'.
%  visu_band: Frequency range for which spectrum is shown.
%  visu_ival: Time interval for which ERD/ERS curves are shown.
%  grd: grid to be used in the grid plots of spectra and ERDs.
%
%Only the figures of the chosen class combination are visible, while the
%others are hidden. In order to make them visible, type
%>> set(cat(2, data.all_results.figure_handles), 'Visible','on')
%
%You might like to modify bbci.feature.ival after running this function.

calibrate = bbci.calibrate;
opt= calibrate.settings;
[opt, isdefault]= ...
    set_defaults(opt, ...
                 'classes', 'auto', ...
                 'clab', {'*'}, ...
                 'classifier','RLDAshrink',...
                 'ref_ival',[-2000 0], ...
                 'ival',[-2000,10000], ...
                 'nIvals', 5, ...
                 'signal','both', ...
                 'doLowpass', 1, ...
                 'derivative', 0, ...
                 'baseline',1, ...
                 'zscore',1, ...
                 'grid',[] , ...
                 'spectra', 1, ...
                 'plot',1 ...
                 );

opt_grid = {'scaleGroup' 'all'};
opt_tit = {'Fontsize' 10 'Fontweight' 'bold'};

%% Build dummy grid for grid_plot
clab = strhead(data.cnt.clab(1:end/2));
clab = {'scale',clab{:},'legend'};
ii=1;
grd = '';
for r=1:ceil(sqrt(numel(clab)-2)) % nr of rows
  for c=1:ceil(sqrt(numel(clab)-2))
    if ii<=numel(clab)
      grd = [grd clab{ii} ','];
    else
      grd = [grd '_,'];
    end
    ii=ii+1;
  end
  grd= [grd(1:end-1) '\n'];
end
grd = sprintf(grd);
mnt= mnt_setGrid(data.mnt, grd);

%% *** Preprocessing ***
if ~isempty(opt.clab)
  data.cnt= proc_selectChannels(data.cnt,opt.clab);
end

% LP filter data [IS IT POSSIBLE TO SET BBCI.SIGNAL AND DO THIS
% AUTOMATICALLY?]
if opt.doLowpass
  fprintf('Lowpass filtering signal\n')
  [filt_b,filt_a]=butter(3, opt.lp*2/data.cnt.fs,'low');
  data.cnt = proc_filtfilt(data.cnt,filt_b,filt_a);
end

% oxy = dat;
% oxy.x = oxy.x(:,1:end/2);
% oxy.clab = oxy.clab(:,1:end/2);
% deoxy = dat;
% deoxy.x = dat.x(:,end/2+1:end);
% deoxy.clab = deoxy.clab(:,end/2+1:end);

if opt.spectra
  % Spectra
  spec_band = [0.01 3];     % Frequency band for spectrum plot
  spec_win = kaiser(250,2);   % window for proc_spectrumspecoxy = proc_spectrum(oxy, spec_band,spec_win);
  spec = proc_spectrum(data.cnt, spec_band,spec_win);
end

% derivative feat?
if opt.derivative
   fprintf('Calculating derivative of NIRS signal\n')

   data.cnt.x=diff(data.cnt.x);
%    deoxy.x=diff(deoxy.x); %*data.cnt.fs;
%    oxy.x=diff(oxy.x);
end

%% Epoch data
data.cnt = cntToEpo(data.cnt,data.mrk,opt.ival);

% baseline
if opt.baseline
  fprintf('Performing baseline correction\n')
  data.cnt = proc_baseline(data.cnt,opt.ref_ival);
end

% oxy.clab = untex(oxy.clab);
% deoxy.clab = untex(deoxy.clab);
%
%% *** Neurophysiology ***
% Split in oxy/deoxy or WL1/WL2
% oxy = data.cnt;
% oxy.x = oxy.x(:,1:end/2,:);
% oxy.clab = oxy.clab(:,1:end/2);
% deoxy = data.cnt;
% deoxy.x = data.cnt.x(:,end/2+1:end,:);
% deoxy.clab = deoxy.clab(:,end/2+1:end);

opt_tit = {'Fontsize' 10 'Fontweight' 'bold'};

%% make plots
if opt.plot
    %% Power spectrum of CNT data
    if opt.spectra
      fig_set(2),clf
      subplot(2,1,1)
      plot(spec.t,spec.x(:,1:end/2)), title(strtail(spec.clab{1}),opt_tit{:})
      xlabel(sprintf('Frequency [%s]',spec.xUnit)),ylabel(sprintf('Power [%s]',spec.yUnit))
      subplot(2,1,2)
      plot(spec.t,spec.x(:,1:end/2)), title(strtail(spec.clab{1}),opt_tit{:})
      xlabel(sprintf('Frequency [%s]',spec.xUnit)),ylabel(sprintf('Power [%s]',spec.yUnit))
      annotation('textbox',[.01 .9 .1 .1 ],'String','Power spectra of CNT data ');
    end

    %% z-score normalization
    if opt.zscore
      fprintf('Calculating z-scores\n')
      m= mean(data.cnt.x,1);   % mean
      s= std(data.cnt.x,1);    % std
      len_x = length(data.cnt.x);
      data.cnt.x= (data.cnt.x - repmat(m,[len_x 1])) ./ repmat(s,[len_x 1]);
      clear m s
    end

    %% Grid plot -- signal
    fig_set(1),clf
    H = grid_plot(proc_selectChannels(data.cnt,1:numel(data.cnt.clab)/2), mnt, defopt_scalp_power,opt_grid{:});
    set(gcf,'Name',['Epoched signal ' strtail(data.cnt.clab{1})])
    % printFigure(['grid_oxy' saveapp], [36 20], opt_fig);

    fig_set(2),clf
    H = grid_plot(proc_selectChannels(data.cnt,(numel(data.cnt.clab)/2)+1:numel(data.cnt.clab)), mnt, defopt_scalp_power,opt_grid{:});
    set(gcf,'Name',['Epoched signal ' strtail(data.cnt.clab{end})])
    % printFigure(['grid_deoxy' saveapp], [36 20], opt_fig);

    %% R-Square plots
    dat_r = proc_r_square_signed(data.cnt,'multiclass_policy','pairwise');
    npairs = size(dat_r.y,1);
    nclab = numel(data.cnt.clab)/2;
    for ii=1:npairs
    for od=[1 2]   %oxy-deoxy
      figure
      imagesc(dat_r.x(:,(od-1)*nclab+1:od*nclab,ii)')
      title([dat_r.className{ii} ' ' strtail(dat_r.clab{od*nclab})],opt_tit{:})
      ylabel('Channel',opt_tit{:}),xlabel('Time',opt_tit{:})
      xt=  str2num(get(gca,'XTickLabel'));
    %   yt=  str2num(get(gca,'YTickLabel'));
      set(gca,'XTickLabel',dat_r.t(xt))
    %   set(gca,'YTickLabel',strhead(oxy_r.clab))
      colorbar
    end
    end
    clear dat_r

    %% Wavelets
    spec = proc_selectChannels(data.cnt,1:nclab);
    spec = proc_specgram(spec, [0.01 3]);

    % data.cnt = proc_selectChannels(data.cnt,1:nclab);
    % wave_band = [.01:.4:3];
    % for w=wave_band
    %   if w==wave_band(1)
    %     wave = proc_wavelets(data.cnt,w,'vectorize',1);
    %   else
    %     w1 = proc_wavelets(data.cnt,w,'vectorize',1);
    %   end
    % end
end
%% *** Classification ***
% BBCI.FEATURE selection
bbci.feature.fcn = {@proc_jumpingMeans};

% class combination
if isequal(opt.classes, 'auto'),
  class_combination= nchoosek(1:size(data.mrk.y,1), 2);
else
  class_combination= find(ismember(data.mrk_all.className, opt.classes,'legacy'));
  if length(class_combination) < length(opt.classes),
    error('Not all specified classes were found.');
  end
  if length(class_combination) ~= 2,
    error('This calibration is only for binary classification.');
  end
end

ival_cfy={};

for ci= 1:size(class_combination,1)

    fv=proc_selectClasses(data.cnt,class_combination(ci,:));
    fprintf('\nComparing classes %s vs %s.',fv.className{1},fv.className{2})

    % find the feature/time intervals
    fvr = proc_r_square_signed(fv);
    [ival_cfy{ci},nfo]= select_time_intervals(fvr,'nIvals',opt.nIvals,'clab_pick_peak', fv.clab,...
      'visualize',0,'intersample_timing',1);

    bbci.feature.param = { {ival_cfy{ci}} };
    % Get feature vector ans train classifier
    fv = bbci_calibrate_evalFeature(fv, bbci.feature);

    bbci.classifier.C= trainClassifier(fv, opt.classifier);

%     % Xvalidation
%     [loss, loss_std]= xvalidation(fv, opt.classifier, ...
%        'sample_fcn', {'divisions', [opt.nShuffles opt.nFolds]},'save_classifier',0);
%
%     mean_loss(ci)= loss;
%     std_loss(ci)= loss_std;

    cfy_fields= {'signal', 'feature', 'classifier', 'quit_condition'};
    bbci_all(ci)= struct_copyFields(bbci, cfy_fields);

    clear fv*
end

%% Choose best binary combination of classes
% Copy parameters into BBCI struct

nComb= size(class_combination,1);
if nComb > 1,
%   [dummy, bi]= min(mean_loss + 0.1*std_loss);
bi=2;
  bbci= struct_copyFields(bbci, bbci_all(bi), cfy_fields);
%     bbci.feature.param = { {ival_cfy{ci}} };
end

comb = class_combination(bi,:);
fprintf('Chosen (best) combination: classes %s vs %s.\n',data.cnt.className{comb(1)},data.cnt.className{comb(2)})
bbci.classifier.classes=data.cnt.className(comb);
%% Set the BBCI struct

% BBCI.SIGNAL preprocessing (really just for online data?)
% bbci.signal.clab = opt.clab;
bbci.signal.clab = data.cnt.clab;

% data.signal.clab = opt.clab;
bbci.signal.proc = {};
% if ~isempty(opt.clab)
%   bbci.signal.proc = {bbci.signal.proc{:} {@proc_selectChannels,opt.clab}};
% end
if opt.doLowpass
  bbci.signal.proc = {bbci.signal.proc{:} {@proc_filtfilt,filt_b,filt_a}};
end
if opt.derivative
  bbci.signal.proc = {bbci.signal.proc{:} {@diff}};
end

% SOMEWHERE: proc_baseline ?

% BBCI.FEATURE selection
bbci.feature.ival = opt.ival;

author: bbci


text:
function b = isPointOnRay(point, ray, varargin)
%ISPOINTONRAY Test if a point belongs to a ray.
%
%   B = isPointOnRay(PT, RAY);
%   Returns 1 if point PT belongs to the ray RAY.
%   PT is given by [x y] and RAY by [x0 y0 dx dy].
%
%   If PT is a N-by-2 array, and RAY is a M-by-4 array, then the result is
%   a N-by-M array containing the result of each pair-wise test.
%
%   B = isPointOnRay(PT, RAY, TOL);
%   Specifies the tolerance to use for testing if point is on the ray.
%
%   Example
%     ray = [10 20 3 4];
%     % test for a point on the ray
%     p1 = [16 28];
%     isPointOnRay(p1, ray)
%     ans =
%       logical
%        0
%     % test for a point on the supporting line but "before" the origin
%     p2 = [7 16];
%     isPointOnRay(p1, ray)
%     ans =
%       logical
%        0
%
%   See also
%   rays2d, points2d, isPointOnLine, isPointOnEdge
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-10-31
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

% extract computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% number of rays and points
Nr = size(ray, 1);
Np = size(point, 1);

% if several rays or several points, adapt sizes of arrays
x0 = repmat(ray(:,1)', Np, 1);
y0 = repmat(ray(:,2)', Np, 1);
dx = repmat(ray(:,3)', Np, 1);
dy = repmat(ray(:,4)', Np, 1);
xp = repmat(point(:,1), 1, Nr);
yp = repmat(point(:,2), 1, Nr);

% test if points belongs to the supporting line
b1 = abs((xp-x0).*dy - (yp-y0).*dx) ./ (dx.*dx + dy.*dy) < tol;

% check if points lie the good direction on the rays
ind     = abs(dx) > abs(dy);
t       = zeros(size(b1));
t(ind)  = (xp(ind) - x0(ind)) ./ dx(ind);
t(~ind) = (yp(~ind) - y0(~ind)) ./ dy(~ind);

% combine the two tests
b = b1 & (t >= 0);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function res = eq(a,b)
%EQ           Implements  a == b  elementwise for intervals a and b
%
% a and b must be either both real or both complex
%

% written  10/16/98     S.M. Rump
% modified 11/30/98     S.M. Rump
% modified 06/24/98     S.M. Rump  multi-dimensional arrays
% modified 09/01/00     S.M. Rump  result array
%                                  comparison real/complex
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if ~isa(a,'intval')
    if b.complex
      res = ( b.mid==a ) & ( b.rad==0 );
    else
      res = ( b.inf==a ) & ( b.sup==a );
    end
    return
  end

  if ~isa(b,'intval')
    if a.complex
      res = ( a.mid==b ) & ( a.rad==0 );
    else
      res = ( a.inf==b ) & ( a.sup==b );
    end
    return
  end

  if ( a.complex ~= b.complex )
    if a.complex
      error('intval comparison == of complex and real')
    else
      error('intval comparison == of real and complex')
    end
  end

  if a.complex
    res = (a.mid==b.mid) & (a.rad==b.rad);
  else
    res = (a.inf==b.inf) & (a.sup==b.sup);
  end

repo: OpenMAS


text:
 % Left Stance Domain
 %
 % Contact: Left Toe
function domain = LeftStance(model, load_path)
    % construct the right stance domain of Cassie
    %
    % Parameters:
    % model: the right body model of Cassie robot

    %% first make a copy of the robot model
    %| @note Do not directly assign with the model variable, since it is a
    %handle object.
    domain = copy(model);
    % set the name of the new copy
    domain.setName('LeftStance');

    % Extract state variables
    q = domain.States.x;

    if nargin < 2
        load_path = [];
    end

    %% Add contact
    % left foot point contact
    [left_foot, fric_coef] = sys.frames.LeftFoot(model);

    p_left_foot = getCartesianPosition(domain, left_foot);
    r_left_foot  = getRelativeEulerAngles(domain, left_foot );

    constr = [p_left_foot(1);p_left_foot(2);p_left_foot(3);r_left_foot(3)];
    hol = HolonomicConstraint(domain, constr, 'LeftFoot',...
        'ConstrLabel',{{'LeftFootX','LeftFootY','LeftFootZ','LeftFootYaw'}},...
        'DerivativeOrder',2);
    domain = addHolonomicConstraint(domain,hol);
    %     domain = addContact(domain,left_foot,fric_coef, geom, load_path);
    if ~isempty(fric_coef)

        f = domain.Inputs.ConstraintWrench.fLeftFoot;
        % get the friction cone constraint
        %         [friction_cone, fc_label, auxdata] = getFrictionCone(left_foot, f, fric_coef);
        mu = SymVariable('mu');
        gamma = SymVariable('gamma');
        % x, y, z, roll, yaw
        constr = [f(3) - 300; % fz >= 0
            f(1) + (mu/sqrt(2))*f(3);  % -mu/sqrt(2) * fz < fx
            -f(1) + (mu/sqrt(2))*f(3); % fx < mu/sqrt(2) * fz
            f(2) + (mu/sqrt(2))*f(3);  % -mu/sqrt(2) * fz < fu
            -f(2) + (mu/sqrt(2))*f(3); % fy < mu/sqrt(2) * fz
            f(4) + gamma * f(3);       % -gamma * fz < wy
            -f(4) + gamma * f(3)];     % wy < gamma * fz

        % create a symbolic function object
        friction_cone = SymFunction(['u_friction_cone_', left_foot.Name],...
            constr,{f},{[mu;gamma]});

        % create the label text
        fc_label = {'normal_force';
            'friction_x_pos';
            'friction_x_neg';
            'friction_y_pos';
            'friction_y_neg';
            'tor_firction_neg';
            'tor_firction_pos';
            };

        auxdata = [fric_coef.mu; fric_coef.gamma];
        % create an unilateral constraint object
        fc_cstr = UnilateralConstraint(domain, friction_cone,...
            ['fc' left_foot.Name], 'fLeftFoot', ...
            'ConstrLabel',{fc_label(:)'},...
            'AuxData',auxdata);
        % add as a set of unilateral constraitns
        domain = addUnilateralConstraint(domain, fc_cstr);
    end

    %% Add event
    % height of non-stance foot (left toe)
    [right_foot_frame] = sys.frames.RightFoot(model);
    p_swingFoot = getCartesianPosition(domain, right_foot_frame);
    h_nsf = UnilateralConstraint(domain,p_swingFoot(3),'rightFootHeight','x');
    domain = addEvent(domain, h_nsf);

    % phase variable: time
    t = SymVariable('t');
    p = SymVariable('p',[2,1]);
    tau = (t-p(1))/(p(2)-p(1));

    % relative degree two outputs:
    ya_2 = [ q('qHLeft');
        (q('qALeft')+q('qBLeft'))./2;
        -q('qALeft')+q('qBLeft');
        q('qHRight');
        (q('qARight')+q('qBRight'))./2;
        -q('qARight')+q('qBRight');];

    y2 = VirtualConstraint(domain,ya_2,'output','DesiredType','Bezier','PolyDegree',5,...
        'RelativeDegree',2,'PhaseType','TimeBased',...
        'PhaseVariable',tau,'PhaseParams',p,'Holonomic',true, 'LoadPath', load_path);

    domain = addVirtualConstraint(domain,y2);

    domain.PreProcess = @sim.LeftStancePreProcess;
end

author: ayonga

repo: frost-dev

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


text:
function [ y,dzdy ] = pad_data_1d(  I,P,dzdy)
%PAD_DATA Summary of this function goes here
%   Detailed explanation goes here

    if(length(P)==1)
        P=ones(1,2)*P;
    end

    if isempty(dzdy)
        %%forward
        original_size=size(I);
        new_size=original_size;
        new_size(1)=new_size(1)+P(1)+P(2);
        y=zeros(new_size,'like',I);
        y(1+P(1):P(1)+original_size(1),:,:)=I;

    else
        %backward
        new_size=size(dzdy);
        original_size=new_size;
        original_size(1)=original_size(1)-P(1)-P(2);
        y=I(1+P(1):P(1)+original_size(1),:,:);
        dzdy=dzdy(1+P(1):P(1)+original_size(1),:,:);

    end

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function planC = uniformizeScanSupInf(planC, tMin, tMax, optS, hBar)
%"uniformizeScanSupInf"
%    Creates the superior and inferior scan arrays so that they
%   are uniform, consistent with the rest of the scan array.
%
%Latest modifications:
% 16 Aug 02, V H Clark, first version.
% 09 Apr 03, JOD, added hBar to input parameter list.
% 18 Feb 05, JRA, Added support for multiple scans.
%
%Usage:
%   function planC = uniformizeScanSupInf(planC, tMin, tMax, optS, hBar)

indexS = planC{end};

for scanNum=1:length(planC{indexS.scan})
    scanStruct = planC{indexS.scan}(scanNum);

	uniformScanInfo = planC{indexS.scan}(scanNum).uniformScanInfo;
	sliceNumSup = uniformScanInfo.sliceNumSup; %superior slice number of original CT scan still being used
	sliceNumInf = uniformScanInfo.sliceNumInf; %inferior slice number of original CT scan still being used
	uniformSliceThickness = uniformScanInfo.sliceThickness;
	scanArray = planC{indexS.scan}(scanNum).scanArray;
	scanInfo = planC{indexS.scan}(scanNum).scanInfo;

	[scanArraySup, scanArrayInf, uniformScanFirstZValue] = uniformizeScanEnds(scanStruct, sliceNumSup, sliceNumInf, uniformSliceThickness, tMin, tMax, optS, hBar);

	uniformScanInfo.firstZValue = uniformScanFirstZValue;
	uniformScanInfo.supInfScansCreated = 1;

	planC{indexS.scan}(scanNum).scanArraySuperior = scanArraySup;
	planC{indexS.scan}(scanNum).scanArrayInferior = scanArrayInf;

	planC{indexS.scan}(scanNum).uniformScanInfo = uniformScanInfo;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function e = vgg_rms_rrror(M)
% e = vgg_rms_rrror(M)
%
% Get RMS diff from zero of matrix or vector M

e = sqrt(sum(sum(M.*M)) / prod(size(M)));

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
%*****************************************************************************
% DSDP5:  Dual Scaling Algorithm for Positive Semidefinite Programming
% Copyright (c) 2002 by
% S. J. Benson, Y. Ye
% Last modified: 20 January 2004
%*****************************************************************************
%
% Converts data from DSDP2 format to DSDP5 format and calls solver.
%
%
% > DSDP(C,A,b) attempts to solve the positive semidefinite program
%      MINIMIZE trace(C*X) SUCH THAT A(:,i)'*X*A(:,i) = b_i, i=1,...,m and X >= 0
%      using a dual scaling algorithm.  The first argument is an n x n symmetric
%      matrix.  The second argument is an n x m matrix, and the third
%      argument b is a dense column vector of length m.
%
% > DSDP(C,A,b,y0) specifies an initial dual vector y0.
%
% > DSDP(C,A,b,CC,AA) adds additional linear constraints to the dual problem.
%
% > DSDP(C,A,b,CC,AA,y0) specifies an initial dual vector y0.
%
% > [Y] = DSDP() returns the dual solution vector.
%
% > [Y,X] = DSDP() returns the primal and dual solution
%
% > [Y,X,V] = DSDP() returns the rank reduction vector.
%
%
%*****************************************************************************

function [Y,DY,K,V] = dsdp2(C,A,b,CC,AA,y0);

  m=length(b);

  AC=cell(2,3);
  n=size(C,1);
  nn=n*(n+1)/2;
  m=length(b);
  AAC=sparse(nn,m+1);
  for i=1:m, ai=A(:,i); AAC=[AAC dvec(ai*ai')]; end;
  AAC=[AAC dvec(C)];
  AC{1,1}='SDP';
  AC{1,2}=size(C,1);
  AC{2,1}='LP';
  AC{2,2}=0;

  if nrhs>4
     [n1,n2]=size(C{j});
     if (n1==1 | n2==1)
       AAC=sparse([]);
       for i=1:m, AAC=[AAC sparse(AA(:,i)]; end;
       AAC=[AAC sparse(CC)'];
       AC{2,1}='LP';
       AC{2,2}=length(CC);
       AC{2,3}=AAC;
       end;
  end;

  [STAT,y,X]=dsdp(b,AC,OPTIONS,y0);

  XX=cell(p,1);
  for j=1:p,
     [n1,n2]=size(C{j})
     if (n1==1 | n2==1)
       XX{j}=X{j};
     else
       XX{j}=dmat(X{j});
     end;
  end;

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [y, details] = newton(functionf, dfunctionf, yest, yold, ...
                               tolabs, tolrel, maxit, printConsole)
%newton Newton iteration with Armijo step control for the nonlinear function
% and Jacobi. Starting from an estimate yest the function computes with a newton
% iteration the solution matrix y of f(y)=0 up to
% tol = tolrel*||f(y)||+tolabs.
%
% Call:   [matrix scalar scalar] =
%                newton(matrix, matrix, scalar, scalar, scalar, structure)
% Input:  matrix yest, size(yest) = (M,N) - estimate of the solution
% Input:  matrix yold, size(yold) = size(y) - parameter of fideq and dfideq
%         (state at previous time level)
% Input:  scalar tolabs>0 - absolute tolerance
% Input:  scalar tolrel>0 - relative tolerance
% Input:  scalar maxit>=1 - maximal number of iterations
% Input:  printConsole - > 0: print iteration information on console output
% Output: matrix y, size(y) = (M,N) - solution
% Output: scalar res, residuum ||f(y)||
% Output: scalar iter, number of iterations
%
% Author: Raimund Wegener, Stefan Schießl
% Date: 12.10.2011

if (nargin < 8)
    printConsole = 1;
end

%%
% initialization
sigma = 1e-2; alpha = 0.5; lambda = 1;

y = yest; f = functionf(y,yold); psi = 0.5*f(:)'*f(:); res = sqrt(2*psi);
tol= tolrel * res + tolabs;

%%
%iteration
[M,N] = size(y);
iter = 0;
% Flag for the message text
armijoMessage = '';
while (res > tol) && (iter < maxit)
    [i,j,S] = dfunctionf(y);
    dfdy = sparse(i,j,S,M*N,M*N);
    % '\' only gives a warning for singular matrices, generate an error
    % for that.
    lastwarn('');
    % bandden:  backslash uses band solver if band density is > bandden.
    % if bandden = 1.0, never use band solver, if bandden = 0.0, always use
    % band solver.
    spparms('bandden',0);
    del = -dfdy\f(:);
    assert(~strcmp(lastwarn, 'Matrix is singular to working precision.'), ...
        'Step:NotSolvable', 'Matrix is singular to working precision.');
    assert(~strcmp(lastwarn, 'Matrix is close to singular or badly scaled.'), ...
        'Step:NotSolvable', 'Matrix is close to singular or badly scaled.');

    y_ = y + del;
    f_ = functionf(y_,yold);
    psi_ = 0.5*f_(:)'*f_(:);
    if (psi_ > (1-2*sigma)*psi)
        lambda = min(lambda/alpha,1);
        y_ = y + lambda*del; f_ = functionf(y_,yold); psi_ = 0.5*f_(:)'*f_(:);
        while psi_ > (1-2*sigma*lambda)*psi
            lambda = lambda*alpha;
            % step control with Armijo rule
            y_ = y + lambda*del; f_ = functionf(y_,yold); psi_ = 0.5*f_(:)'*f_(:);

        end
        armijoMessage = ['newton/armijo: iter=' num2str(iter) ',\t lambda=' num2str(lambda) '\n'];
        if (printConsole > 1)
            fprintf(armijoMessage)
        end
    end

    y = y_; f = f_; psi = psi_; res = sqrt(2*psi);
    iter = iter + 1;
end
message = ['newton: res=' num2str(res) ', iter=' num2str(iter) '\n' armijoMessage];
if (printConsole > 0)
    fprintf(message)
end

if (isnan(res))
   error('Newton:NoConvergence', 'No solution found.');
end

if (res > tol)
   error('Newton:AccuracyNotReached', ['newton/armijo: Results are inaccurate (Res: ' num2str(res) ').']);
end

details.res = res;
details.iter = iter;
details.message = message;

end

author: Sable

repo: mcbench-benchmarks


text:
function seed = latin_edge_test01 ( seed )

%*****************************************************************************80
%
%% LATIN_EDGE_TEST01 tests LATIN_EDGE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, an initial seed for the random number generator.
%
%    Output, integer SEED, the updated random number seed.
%
  dim_num = 2;
  point_num = 11;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LATIN_EDGE_TEST01\n' );
  fprintf ( 1, '  LATIN_EDGE chooses a Latin cell arrangement,\n' );
  fprintf ( 1, '  which includes the edge points.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension = %d\n', dim_num );
  fprintf ( 1, '  Number of points =  %d\n', point_num );
  fprintf ( 1, '  Using seed = %d\n', seed );

  [ x, seed ] = latin_edge ( dim_num, point_num, seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The Latin Edge Square points:\n' );
  fprintf ( 1, '\n' );

  for j = 1: point_num
    for i = 1: dim_num
      fprintf ( 1, '%10f  ', x(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function labelImg = segToLabel(u, conn)
%SEGTOLABEL Converts the output image of a segmentation function (e.g. minL2Potts2DADMM) to a
%label image with integer labels (each connected component gets a unique integer)
% u: image to label (typically the output of minL2Potts2DADMM))
% conn: connectivity as in bwconncomp (default: 8)

if ~exist('conn', 'var')
    conn = 8;
end

[m,n,c] = size(u);
vals = unique(reshape(u, m*n, c), 'rows');
labelImg = zeros(m, n);
labelCounter = 0;
for i=1:size(vals,1)
    bw_mult = (vals(i,:) == reshape(u, m*n, c));
    bw = bw_mult(:,1);
    CC = bwconncomp(reshape(double(bw), m, n), conn);
    for j = 1:CC.NumObjects
        labelImg(CC.PixelIdxList{j}) = labelCounter;
        labelCounter = labelCounter + 1;
    end
end

end

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
function output = calllpsolve(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
integer_variables = interfacedata.integer_variables;
binary_variables = interfacedata.binary_variables;
semicont_variables = interfacedata.semicont_variables;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

n = length(c);
% Bounded variables converted to constraints
if ~isempty(ub)
    LB = lb;
    UB = ub;
else
    if isempty(integer_variables) & isempty(binary_variables)
        LB = -ones(n,1)*inf;;
        UB = ones(n,1)*inf;;
    else
        %LP_SOLVE FAILS IF BOUNDS NOT EXPLICIT
        [LB,UB,used_rows] = find_lp_bounds(F_struc,K);
        LB(isinf(LB)) = -1e12;
        UB(isinf(UB)) = 1e12;
        F_struc(K.f+used_rows,:)=[];
        K.l = K.l - length(used_rows);
        LB(binary_variables) = max(LB(binary_variables),0);
        UB(binary_variables) = min(UB(binary_variables),0);
    end
end

if options.showprogress;showprogress('Calling LPSOLVE',options.showprogress);end

f = - full(c);            % Must be full
A = - F_struc(:,2:end);
b = full(F_struc(:,1));   % Must be full
e = -ones(size(A,1),1);
e(1:K.f) = 0;

xint = uniquestripped([integer_variables binary_variables]);

% Call mex-interface
solvertime = clock;

if K.f>0
    Aeq = A(1:K.f,:);
    beq = b(1:K.f);
    Alin = A(K.f+1:end,:);
    blin = b(K.f+1:end);
    [LB,UB,Alin,blin] = remove_bounds_from_Ab(Alin,blin,LB,UB);
    A = [Aeq;Alin];
    b = [beq;blin];
    e = -ones(size(A,1),1);
    e(1:K.f) = 0;
    options.saveduals = 0;
else
    [LB,UB,A,b] = remove_bounds_from_Ab(A,b,LB,UB);
    e = -ones(size(A,1),1);
    options.saveduals = 0;
end
if ~isempty(semicont_variables)
    redundant = find(LB<=0 & UB>=0);
    semicont_variables = setdiff(semicont_variables,redundant);
end

% LPSOLVE assumes semi-continuous variables only can take negative values so
% we negate semi-continuous violating this
NegativeSemiVar = [];
if ~isempty(semicont_variables)
    NegativeSemiVar = find(UB(semicont_variables) < 0);
    if ~isempty(NegativeSemiVar)
        temp = UB(semicont_variables(NegativeSemiVar));
        UB(semicont_variables(NegativeSemiVar)) = -LB(semicont_variables(NegativeSemiVar));
        LB(semicont_variables(NegativeSemiVar)) = -temp;
        A(:,semicont_variables(NegativeSemiVar)) = -A(:,semicont_variables(NegativeSemiVar));
        f(semicont_variables(NegativeSemiVar)) = -f(semicont_variables(NegativeSemiVar));
    end
end
if options.savedebug
    save lpsolvedebug f A b e UB LB xint
end

lp = create_lp_solve_model(A,b,f,xint,LB,UB,e,options);
if ~isempty(K.sos)
    for i = 1:length(K.sos.type)
       lp_solve('add_SOS', lp, 'Dummy', str2num(K.sos.type(i)), 1, K.sos.variables{i}, K.sos.weight{i});
    end
end
if ~isempty(semicont_variables)
    lp_solve('set_semicont', lp, semicont_variables)
end

try
    if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
    solvertime = tic;
    result=lp_solve('solve', lp);
    solvertime = toc(solvertime);
    if result == 0 | result == 1 | result == 11 | result == 12
        [obj, x, duals] = lp_solve('get_solution', lp);
    else
        obj = [];
        x = zeros(length(c),1);
        duals = [];
    end
    lp_solve('delete_lp', lp);
catch
    obj = [];
    x = zeros(length(c),1);
    duals = [];
    result = -1;
    lp_solve('delete_lp', lp);
end

% LPSOLVE assumes semi-continuous variables only can take negative values so
% we negate semi-continuous violating this
 if length(x) == length(c)
     if ~isempty(NegativeSemiVar)
        x(NegativeSemiVar) = -x(NegativeSemiVar);
     end
 end

if options.saveduals & isempty(integer_variables)
    D_struc = duals;
else
    D_struc = [];
end

switch result
    case {0,1}
        problem = 0; % OPTIMAL
    case 2
        problem = 1; % INFEASIBLE
    case 3
        problem = 2; % UNBOUNDED
    case {7,12,13}
        problem = 3; % RUN OUT OF TIME OR SIMILIAR
    case 5
        problem = 4;
    case {-2,10,11}
        problem = 11;
    otherwise
        problem = -1;
end

% Save all data sent to solver?
if options.savesolverinput
	solverinput.A = A;
	solverinput.f = f;
	solverinput.b = b;
	solverinput.LB = LB;
	solverinput.UB = UB;
    solverinput.xint = xint;
else
	solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
	solveroutput.x = x;
    solveroutput.obj = obj;
    solveroutput.duals = duals;
    solveroutput.result = result;
else
	solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% file: get_prop_diff_ci1.m
% one-sided confidence intervals for the difference of proportions

% 030113 tdr created

function ci = get_prop_diff_ci1(x1,n1,x2,n2,alpha,method,verbose)

p1_hat = x1/n1; p2_hat = x2/n2; delta_p_hat = p1_hat - p2_hat;
a = find_lower_limit(x1,n1,x2,n2, alpha);
b = find_upper_limit(x1,n1,x2,n2, alpha);
ci = [delta_p_hat a b];

% ------------------------------------------------------------
% find lower limit

function limit = find_lower_limit(x1,n1,x2,n2,alpha)
% alpha is 1-Pr{delta_p >= a}
p1_hat = x1/n1; p2_hat = x2/n2; delta_p_hat = p1_hat - p2_hat;
too_small = -1; too_big = 1;

tolerance = 1e-6;
max_count = 50;

count = 0;
guess = (too_big + too_small)/2;
alpha_guess = 1- prop_diff(x1,n1,x2,n2,guess);
while (abs(alpha_guess - alpha) > tolerance & (count < max_count))
    if (alpha_guess > alpha)
        too_big = guess;
        guess = (too_small + guess)/2;
    else
        too_small=guess;
        guess = (too_big + guess)/2;
    end;
    count = count+1;
    alpha_guess = 1 - prop_diff(x1,n1,x2,n2,guess);
end;

limit = guess;

% ------------------------------------------------------------
% find upper limit

function limit = find_upper_limit(x1,n1,x2,n2,alpha)
% alpha is Pr{delta_p >= b}
p1_hat = x1/n1; p2_hat = x2/n2; delta_p_hat = p1_hat - p2_hat;
too_small = -1; too_big = 1;

tolerance = 1e-6;
max_count = 50;
count = 0;
guess = (too_big + too_small)/2;
alpha_guess = prop_diff(x1,n1,x2,n2,guess);
while (abs(alpha_guess - alpha) > tolerance & (count < max_count))
    if (alpha_guess < alpha)
        too_big = guess;
        guess = (too_small + guess)/2;
    else
        too_small=guess;
        guess = (too_big + guess)/2;
    end;
    count = count+1;
    alpha_guess = prop_diff(x1,n1,x2,n2,guess);
end;

limit = guess;


text:
function A = warmUpExercise()
%WARMUPEXERCISE Example function in octave
%   A = WARMUPEXERCISE() is an example function that returns the 5x5 identity matrix

A = [];
% ============= YOUR CODE HERE ==============
% Instructions: Return the 5x5 identity matrix
%               In octave, we return values by defining which variables
%               represent the return values (at the top of the file)
%               and then set them accordingly.

A = eye(5);

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

end

repo: machine-learning-coursera


function [ DP_cost, DP_time, av_cost, av_time ] = DP_cut_k( cutCostsFile, cut_num, search_k)
%DP_CUT Summary of this function goes here
%   Detailed explanation goes here

    file_t = fopen(cutCostsFile);
    cut_costs = fscanf(file_t, '%d');
    fclose(file_t);

    % get the initial averaging-cut

    cut_length = length(cut_costs) / cut_num;

    for i = 1:cut_num-1
       initial_cuts(i,:) = floor(i*cut_length);
    end

    % start Dynamic Programming
    % find the best/min value in the domain near cut_points as initial
    % values

    tic;
    best_values = [];
    best_cuts = [];
    search_domain = floor(cut_length / search_k);

    for i = 1:length(initial_cuts)

        domain = cut_costs(initial_cuts(i,:)-search_domain:initial_cuts(i,:)+search_domain,:);

        [min_value, min_cut] = min(domain);

        best_values(i,:) = min_value;

        best_cuts(i,:) = initial_cuts(i,:) + (min_cut - search_domain - 1);

    end

    DP_cost = sum(best_values);
    DP_time = toc;

    tic
    % get average cost sum
    av_sum = 0;
    for i = 1:length(initial_cuts)
       av_sum = av_sum + cut_costs(initial_cuts(i));
    end
    av_cost = av_sum;
    av_time = toc;

end


text:
% Author: Guosheng Lin (guosheng.lin@gmail.com)

% This is an example of resume training with lower learning rate.

function demo_refinenet_train_reduce_learning_rate()

rng('shuffle');

addpath('./my_utils');
dir_matConvNet='../libs/matconvnet/matlab';
run(fullfile(dir_matConvNet, 'vl_setupnn.m'));

run_config=[];
ds_config=[];

run_config.use_gpu=true;
% run_config.use_gpu=false;
run_config.gpu_idx=1;

% use a random model name:
% model_name=['model_' datestr(now, 'YYYYmmDDHHMMSS')];

% or specify a model_name for training, if the cached snapshot files existed, it will resume training
% model_name='model_20161219094311_example';

% update the model name when reducing learning rate for training, e.g.,
model_name='model_20161219094311_example_epoch70_low_learn_rate';

ds_name='voc2012_trainval';
gen_ds_info_fn=@my_gen_ds_info_voc;

% control the size of input images
run_config.input_img_short_edge_min=450;
run_config.input_img_short_edge_max=1100;

run_config.input_img_scale=1;

ds_config.use_dummy_gt=false;
run_config.use_dummy_gt=ds_config.use_dummy_gt;
ds_config.use_custom_data=false;

ds_config.ds_name=ds_name;
ds_config.gen_ds_info_fn=gen_ds_info_fn;
ds_config.ds_info_cache_dir=fullfile('../datasets', ds_name);

run_config.root_cache_dir=fullfile('../cache_data', ds_name, model_name);
mkdir_notexist(run_config.root_cache_dir);

run_config.model_name=model_name;

diary_dir=run_config.root_cache_dir;
mkdir_notexist(diary_dir);
diary(fullfile(diary_dir, 'output.txt'));
diary on

run_dir_name=fileparts(mfilename('fullpath'));
[~, run_dir_name]=fileparts(run_dir_name);
run_config.run_dir_name=run_dir_name;
run_config.run_file_name=mfilename();

fprintf('\n\n\n=====================================================================\n');
disp('gen ds_info');

disp('ds_config:');
disp(ds_config);

ds_info=gen_dataset_info(ds_config);

my_diary_flush();

fprintf('\n\n\n=====================================================================\n');
disp('run network');

run_config.run_evaonly=false;

% settings for training:

% run_config.trained_model_path=[];
% run_config.learning_rate=5e-4;

% init from a cached model
% e.g., can be used for resuming network training with a lower learning rate by
% initilizing from a cached model.
% the following is an example for loading a cached model and using a lower learning
% to continue the training

run_config.trained_model_path='../cache_data/voc2012_trainval/model_20161219094311_example/model_cache/epoch_70';
run_config.learning_rate=5e-5;

run_config.epoch_num=600;

% turn on this option to cache all data into memory, if it's possible
% run_config.cache_data_mem=true;
run_config.cache_data_mem=false;

% random crop training:
run_config.crop_box_size=400;

% for cityscape, using a larger crop:
% run_config.crop_box_size=600;

% evaluate step: do evaluation every 5 epochs, can be set to 10:
run_config.eva_run_step=5;
run_config.snapshot_step=1;

% choose ImageNet pre-trained resnet:
run_config.init_resnet_layer_num=50;
% run_config.init_resnet_layer_num=101;
% run_config.init_resnet_layer_num=152;

% generate network
run_config.gen_network_fn=@gen_network_main;

run_config.gen_net_opts_fn=@gen_net_opts_model_type1;

% uncomment the following for debug:
% select a subset for both training and evaluation.
% ds_info.train_idxes=ds_info.train_idxes(1:10);
% ds_info.test_idxes=ds_info.train_idxes;
% run_config.snapshot_step=10;

train_opts=run_config.gen_net_opts_fn(run_config, ds_info.class_info);

disp('train_opts:');
disp(train_opts);

disp('train_opts.eva_param:');
disp(train_opts.eva_param);

my_diary_flush();

imdb=my_gen_imdb(train_opts, ds_info);

data_norm_info=[];
data_norm_info.image_mean=128;

imdb.ref.data_norm_info=data_norm_info;

if run_config.use_gpu
	gpu_num=gpuDeviceCount;
	if gpu_num>=1
		gpuDevice(run_config.gpu_idx);
    else
        error('no gpu found!');
	end
end

[net_config, net_exp_info]=prepare_running_model(train_opts);

% net_config can be changed here before running the network

my_net_tool(train_opts, imdb, net_config, net_exp_info);

fprintf('\n\n--------------------------------------------------\n\n');
disp('results are saved in:');
disp(run_config.root_cache_dir);

my_diary_flush();
diary off

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function xf = shift_sample(xf, shift, kx, ky)

% Shift a sample in the Fourier domain. The shift should be normalized to
% the range [-pi, pi].

shift_exp_y =exp((1i * shift(1)) * ky);
shift_exp_x = exp((1i * shift(2)) * kx);
xf =bsxfun(@times, bsxfun(@times, xf, shift_exp_y), shift_exp_x);

author: Daikenan

repo: ASRCF

save_path: github-repos/MATLAB/Daikenan-ASRCF