text:
function [f,x,u]=ksdensityw(y,w,varargin)
%KSDENSITY Compute density estimate
%   [F,XI]=KSDENSITY(X) computes a probability density estimate of the sample
%   in the vector X.  F is the vector of density values evaluated at the
%   points in XI.  The estimate is based on a normal kernel function, using a
%   window parameter (bandwidth) that is a function of the number of points
%   in X.  The density is evaluated at 100 equally-spaced points covering
%   the range of the data in X.
%
%   F=KSDENSITY(X,XI) specifies the vector XI of values where the density
%   estimate is to be evaluated.
%
%   [F,XI,U]=KSDENSITY(...) also returns the bandwidth of the kernel smoothing
%   window.
%
%   DWH: w are the weights for y (w = 1/n for ksdenstiy)
%
%   [...]=KSDENSITY(...,'PARAM1',val1,'PARAM2',val2,...) specifies parameter
%   name/value pairs to control the density estimation.  Valid parameters
%   are the following:
%
%      Parameter    Value
%      'kernel'     The type of kernel smoother to use, chosen from among
%                   'normal' (default), 'box', 'triangle', and
%                   'epanechinikov'.
%      'npoints'    The number of equally-spaced points in XI.
%      'width'      The bandwidth of the kernel smoothing window.  The default
%                   is optimal for estimating normal densities, but you
%                   may want to choose a smaller value to reveal features
%                   such as multiple modes.
%
%   In place of the kernel functions listed above, you can specify another
%   function by using @ (such as @normpdf) or quotes (such as 'normpdf').
%   The function must take a single argument that is an array of distances
%   between data values and places where the density is evaluated, and
%   return an array of the same size containing corresponding values of
%   the kernel function.
%
%   Example:
%      x = [randn(30,1); 5+randn(30,1)];
%      [f,xi] = ksdensity(x);
%      plot(xi,f);
%   This example generates a mixture of two normal distributions, and
%   plots the estimated density.
%
%   See also HIST, @.

% Reference:
%   A.W. Bowman and A. Azzalini (1997), "Applied Smoothing
%      Techniques for Data Analysis," Oxford University Press.

%   Copyright 1993-2002 The MathWorks, Inc.
%   $Revision: 1.7 $  $Date: 2002/03/21 20:36:29 $

% Get y vector and its dimensions
if (prod(size(y)) > length(y)), error('X must be a vector'); end
y = y(:);
y(isnan(y)) = [];
n = length(y);
ymin = min(y);
ymax = max(y);

% Maybe x was specified, or maybe not
if ~isempty(varargin)
   if ~ischar(varargin{1})
      x = varargin{1};
      varargin(1) = [];
   end
end

% Process additional name/value pair arguments
okargs = {'width' 'npoints' 'kernel'};
defaults = {[] [] 'normal'};
[emsg,u,m,kernel] = statgetargs(okargs, defaults, varargin{:});
error(emsg);

% Default window parameter is optimal for normal distribution
if (isempty(u)),
   med = median(y);
   sig = median(abs(y-med)) / 0.6745;
   if sig<=0, sig = ymax-ymin; end
   if sig>0
      u = sig * (4/(3*n))^(1/5);
   else
      u = 1;
   end
end

% Check other arguments or get defaults.
if ~exist('x','var')
   if isempty(m), m=100; end
   x = linspace(ymin-2*u, ymax+2*u, m);
elseif (prod(size(x)) > length(x))
   error('XI must be a vector');
end
xsize = size(x);
x = x(:);
m = length(x);

okkernels = {'normal' 'epanechinikov' 'box' 'triangle'};
if isempty(kernel)
   kernel = okkernels{1};
elseif ~(isa(kernel,'function_handle') | isa(kernel,'inline'))
   if ~ischar(kernel)
      error('Smoothing kernel must be a function.');
   end
   knum = strmatch(lower(kernel), okkernels);
   if (length(knum) == 1)
      kernel = okkernels{knum};
   end
end

blocksize = 1e6;
if n*m<=blocksize
   % Compute kernel density estimate in one operation
   z = (repmat(x',n,1)-repmat(y,1,m))/u;
   w2 = repmat(w, 1, m);
   f = sum(feval(kernel, z).*w2,1);
else
   % For large vectors, process blocks of elements as a group
   M = max(1,ceil(blocksize/n));
   mrem = rem(m,M);
   if mrem==0, mrem = min(M,m); end
   x = x';

   f = zeros(1,m);
   ii = 1:mrem;
   z = (repmat(x(ii),n,1)-repmat(y,1,mrem))/u;
   w2 = repmat(w, 1, mrem);
   f(ii) = sum(feval(kernel, z).*w2,1);
   z = zeros(n,M);
   w2 = zeros(n,M);
   for j=mrem+1:M:m
      ii = j:j+M-1;
      z(:) = (repmat(x(ii),n,1)-repmat(y,1,M))/u;
      w2(:) = repmat(w, 1, M);
      f(ii) = sum(feval(kernel, z).*w2,1);
   end
end

f = reshape(f./u, xsize);

% -----------------------------
% The following are functions that define smoothing kernels k(z).
% Each function takes a single input Z and returns the value of
% the smoothing kernel.  These sample kernels are designed to
% produce outputs that are somewhat comparable (differences due
% to shape rather than scale), so they are all probability
% density functions with unit variance.
%
% The density estimate has the form
%    f(x;k,h) = mean over i=1:n of k((x-y(i))/h) / h

function f = normal(z)
%NORMAL Normal density kernel.
%f = normpdf(z);
f = exp(-0.5 * z .^2) ./ sqrt(2*pi);

function f = epanechinikov(z)
%EPANECHINIKOV Epanechinikov's asymptotically optimal kernel.
a = sqrt(5);
z = max(-a, min(z,a));
f = .75 * (1 - .2*z.^2) / a;

function f = box(z)
%BOX    Box-shaped kernel
a = sqrt(3);
f = (abs(z)<=a) ./ (2 * a);

function f = triangle(z)
%TRIANGLE Triangular kernel.
a = sqrt(6);
z = abs(z);
f = (z<=a) .* (1 - z/a) / a;

author: zouchuhang

repo: LayoutNet


text:
function [gx] = g_nl0(x,P,u,in)

try
    gx = P(1)*exp(P(2)*in.x);
catch
    gx = x(1)*exp(x(2)*P);
end

author: MBB-team

repo: VBA-toolbox


text:
function [str,iter]=display_simple(comp,tab,iter)

  global display_tree_array_indexing;
  if isempty(display_tree_array_indexing)
      display_tree_array_indexing=2;
  end;
  global display_tree_index_at_front;
  if isempty(display_tree_index_at_front)
      display_tree_index_at_front=1;
  end;

  if nargin==3
   if isempty(iter)
    iter{1}=[];
   else
      i=length(iter);
      iter{i+1}=iter{i};  %% <--- new member of tree
   end

   if display_tree_array_indexing==0
     if isempty(iter{length(iter)})
      index=[];
     else
      index=[num2str(length(iter)-1) ':'];
     end
     if display_tree_index_at_front
      str=[ index get_name(comp) '\n'];
     else
      str=[ get_name(comp) ' (' index(length(index))  ')\n'];
     end;
   else
     index=[];
     if ~isempty(iter{length(iter)})
      index= num2str(iter{length(iter)});
      a=find(index==32);
      a=a(1:2:length(a));
      a=setdiff([1:length(index)],a);
      index=index(a);

      if display_tree_array_indexing==2
         index= num2str(iter{length(iter)});
	     fTmp=find(index==' ');
         fTmp=max(setdiff(fTmp,length(index)));
	        if ~isempty(fTmp)
	            index=index(fTmp+1:length(index));
        	end
      end
      index=['[' index ']:'];
     end;

     if display_tree_index_at_front
      str=[ index  get_name(comp) '\n'];
     else
      str=[ get_name(comp)  ' ' index(1:length(index)-1)  '\n'];
     end;
   end

  else
   disp([get_name(comp)]);
  end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [u,eqn,info,node,elem,Neumann] = fracLap1d(square,h,pde,option)

global s
alpha = 1-2*s;
if s == 0.5
    gamma = 1;
else
    gamma = 3/(2*s)+0.1;
end

tic;
%% Mesh
if option.gradmesh
   [node,elem,y] = squaregradmeshquad(square,h,gamma);
else
   [node,elem] = squarequadmesh(square,h);
    y0 = square(3); y1 = square(4);
    y = (y0:h:y1)';
end
N = size(node,1); NT = size(elem,1); Ndof = N;
% x0 = square(1); x1 = square(2);

%% Compute Local matrix
hy = diff(y);
NTy = length(hy);
NTx = NT/NTy;
% stiffness matrix in y direction
Ay = zeros(NTy,3);
intyalpha = diff(y.^(alpha+1)/(alpha+1));
Ay(:,1) = intyalpha./hy.^2;  % Ay(1,1)
Ay(:,2) = -Ay(:,1);          % Ay(1,2) and Ay(2,1)
Ay(:,3) = Ay(:,1);           % Ay(2,2)
% mass matrix in y direction
My = zeros(NTy,3);
a3 = diff(y.^(alpha+3)/(alpha+3));
a2 = diff(y.^(alpha+2)/(alpha+2));
a1 = diff(y.^(alpha+1)/(alpha+1));
yleft = y(1:end-1);
yright = y(2:end);
My(:,1) = (a3 - 2*yleft.*a2 + yleft.^2.*a1)./hy.^2;  % My(1,1)
My(:,2) = (-a3 + (yleft+yright).*a2 - yleft.*yright.*a1)./hy.^2;    % My(1,2) and My(2,1)
My(:,3) = (a3 - 2*yright.*a2 + yright.^2.*a1)./hy.^2;  % My(2,2)
% stiffness matrix in x direction
hx = h;
Ax = 1/hx*[1 -1; -1 1];
Mx = hx/6*[2 1; 1 2];

%% Assemble stiffness matrix
% index map
ixiy = [1 1; 2 1; 2 2; 1 2];
ij2k(1,1) = 1; ij2k(1,2) = 2; ij2k(2,1) = 2; ij2k(2,2) = 3;
% generate sparse pattern
ii = zeros(10*NT,1); jj = zeros(10*NT,1);
index = 0;
for i = 1:4
    for j = i:4
        ii(index+1:index+NT) = double(elem(:,i));
        jj(index+1:index+NT) = double(elem(:,j));
        index = index + NT;
    end
end
% compute non-zeros
sA = zeros(10*NT,1);
index = 0;
for i = 1:4
    for j = i:4
        ix = ixiy(i,1); iy = ixiy(i,2);
        jx = ixiy(j,1); jy = ixiy(j,2);
        Aij = Ax(ix,jx)*My(:,ij2k(iy,jy)) + Mx(ix,jx)*Ay(:,ij2k(iy,jy));
        sA(index+1:index+NT,1) = repmat(Aij,NTx,1);
        index = index + NT;
    end
end

% assemble the matrix
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Ndof,Ndof);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Ndof,Ndof);
A = A + AU + AU';
clear Aij ii jj AU

%% Right hand side
% f = 0
b = zeros(Ndof,1);

%% Set up boundary conditions
% find boundary nodes and Neumann edges
nx = NTx + 1;
ny = NTy + 1;
fixedNode = [1:ny (2:nx-1)*ny (nx-1)*ny+(1:ny)];
bottom = transpose(1:ny:(nx-1)*ny+1);
Neumann = [bottom(1:end-1) bottom(2:end)];
isBdNode = false(N,1);
isBdNode(fixedNode) = true;
freeNode = find(~isBdNode);

%
% bdFlag = setboundary(node,elem,'Dirichlet','abs(y)>eps','Neumann','abs(y)<=eps');
% [fixedNode,Neumann,isBdNode] = findboundary(elem,bdFlag);

% Modify the matrix to include the Dirichlet boundary condition
bdidx = zeros(Ndof,1);
bdidx(fixedNode) = 1;
Tbd = spdiags(bdidx,0,Ndof,Ndof);
T = spdiags(1-bdidx,0,Ndof,Ndof);
A = T*A*T + Tbd;

% Neumann boundary condition
el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
if ~isfield(option,'gNquadorder')
    option.gNquadorder = 4;   % default order exact for linear gN
end
[lambdagN,weightgN] = quadpts1(option.gNquadorder);
phigN = lambdagN;                 % linear bases
nQuadgN = size(lambdagN,1);
ge = zeros(size(Neumann,1),2);
for pp = 1:nQuadgN
    % quadrature points in the x-y coordinate
    ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
         + lambdagN(pp,2)*node(Neumann(:,2),:);
    gNp = pde.g_N(ppxy);
    for igN = 1:2
        ge(:,igN) = ge(:,igN) + weightgN(pp)*phigN(pp,igN)*gNp;
    end
end
ge = ge.*repmat(el,1,2);
b = b + accumarray(Neumann(:), ge(:),[Ndof,1]);

% Modify right handside for Dirichlet boundary condition
% Neumann edges are considered as open set. So the corner points should be
% set as Dirichlet boundary condition!
b(fixedNode) = 0;

eqn = struct('A',A,'b',b,'freeNode',freeNode);

%% Record assembling time
assembleTime = toc;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve
u = zeros(Ndof,1);
switch option.solver
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
        return
    case 'direct'
        u(freeNode) = A(freeNode,freeNode)\b(freeNode);
        residual = norm(b - A*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'amg'
        option.solver = 'CG';
        [u(freeNode),info] = amg(A(freeNode,freeNode),b(freeNode),option);
end
info.assembleTime = assembleTime;

%% Compute error using boundary integral
[lambdagN,weightgN] = quadpts1(option.gNquadorder);
phigN = lambdagN;                 % linear bases
nQuadgN = size(lambdagN,1);
err = zeros(size(Neumann,1),1);
for pp = 1:nQuadgN
    % quadrature points in the x-y coordinate
    ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
         + lambdagN(pp,2)*node(Neumann(:,2),:);
    uhp = u(Neumann(:,1))*phigN(pp,1) + u(Neumann(:,2))*phigN(pp,2);
    up = pde.exactu(ppxy);
    gNp = pde.g_N(ppxy);
    err = err + weightgN(pp)*gNp.*(up - 2*uhp);
end
err = sum(err.*el) + u'*A*u;
errH1 = sqrt(abs(err));
info.errH1 = errH1;

repo: ifem


text:
function varargout = flipdim(varargin)
%FLIPDIM   Flip/reverse a CHEBFUN2 in a chosen direction.
%   G = FLIPDIM(F, DIM) returns a CHEBFUN2 G with the same domain as F but
%   reversed in a direction, i.e., G(x,y)=F(x, c+d-y). If DIM = 2 (default) then
%   G(x,y) = F(x, c+d-y).  Otherwise DIM = 1 and G(x,y) = F(a+b-x, y). The
%   domain of F is [a, b, c, d].
%
% See also FLIPLR, FLIPUD.

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

[varargout{1:nargout}] = flipdim@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function test_bug1443

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_rejectcomponent ft_componentanalysis

load(dccnpath('/home/common/matlab/fieldtrip/data/test/latest/raw/meg/preproc_ctf151.mat'));

cfg = [];
cfg.method = 'fastica';
cfg.numcomponent = 14; % to make it go fast
cfg.randomseed = 13; % so we get the same output each time
comp = ft_componentanalysis(cfg,data);

assert(ft_senstype(comp,'ctf151'))

cfg = [];
cfg.component = 2; % chosen randomly
rej1 = ft_rejectcomponent(cfg, comp, data);
rej2 = ft_rejectcomponent(cfg, comp);

norm(rej2.grad.tra-rej1.grad.tra)/norm(rej2.grad.tra);
figure; imagesc(rej2.grad.tra - rej1.grad.tra); caxis([-1 1])

load standard_sourcemodel3d10mm
load standard_singleshell

cfg = [];
cfg.sourcemodel = sourcemodel;
cfg.headmodel = vol;

cfg.grad = rej1.grad;
grid1 = ft_prepare_leadfield(cfg, rej1);

cfg.grad = rej2.grad;
grid2 = ft_prepare_leadfield(cfg, rej2);

assert(~isequaln(grid1.leadfield,grid2.leadfield))

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function K = kernelD(xp0,yp0,len)

D  = size(xp0,1);
N  = size(xp0,2);
M  = size(yp0,2);

% split M into chunks if on GPU to reduce memory usage
if isa(xp0,'gpuArray')
    K=gpuArray.zeros(N,M);
    cs  = 60;
elseif N > 10000
    K = zeros(N,M);
    cs = 10000;
else
    K= zeros(N,M);
    cs  = M;
end

for i = 1:ceil(M/cs)
    ii = [((i-1)*cs+1):min(M,i*cs)];
    mM = length(ii);
    xp = repmat(xp0,1,1,mM);
    yp = reshape(repmat(yp0(:,ii),N,1),D,N,mM);

    Kn = exp( -sum(bsxfun(@times,(xp - yp).^2,1./(len.^2))/2,1));
    K(:,ii)  = squeeze(Kn);

end

author: cortex-lab


text:
%JTRAJ Compute a joint space trajectory
%
% [Q,QD,QDD] = JTRAJ(Q0, QF, M) is a joint space trajectory Q (MxN) where the joint
% coordinates vary from Q0 (1xN) to QF (1xN).  A quintic (5th order) polynomial is used
% with default zero boundary conditions for velocity and acceleration.
% Time is assumed to vary from 0 to 1 in M steps.  Joint velocity and
% acceleration can be optionally returned as QD (MxN) and QDD (MxN) respectively.
% The trajectory Q, QD and QDD are MxN matrices, with one row per time step,
% and one column per joint.
%
% [Q,QD,QDD] = JTRAJ(Q0, QF, M, QD0, QDF) as above but also specifies
% initial QD0 (1xN) and final QDF (1xN) joint velocity for the trajectory.
%
% [Q,QD,QDD] = JTRAJ(Q0, QF, T) as above but the number of steps in the
% trajectory is defined by the length of the time vector T (Mx1).
%
% [Q,QD,QDD] = JTRAJ(Q0, QF, T, QD0, QDF) as above but specifies initial and
% final joint velocity for the trajectory and a time vector.
%
% Notes::
% - When a time vector is provided the velocity and acceleration outputs
%   are scaled assumign that the time vector starts at zero and increases
%   linearly.
%
% See also QPLOT, CTRAJ, SerialLink.jtraj.

% Copyright (C) 1993-2017, 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 [qt,qdt,qddt] = jtraj(q0, q1, tv, qd0, qd1)
    if length(tv) > 1
        tscal = max(tv);
        t = tv(:)/tscal;
    else
        tscal = 1;
        t = (0:(tv-1))'/(tv-1); % normalized time from 0 -> 1
    end

    q0 = q0(:);
    q1 = q1(:);

    if nargin == 3
        qd0 = zeros(size(q0));
        qd1 = qd0;
    elseif nargin == 5
        qd0 = qd0(:);
        qd1 = qd1(:);
    else
        error('incorrect number of arguments')
    end

    % compute the polynomial coefficients
    A = 6*(q1 - q0) - 3*(qd1+qd0)*tscal;
    B = -15*(q1 - q0) + (8*qd0 + 7*qd1)*tscal;
    C = 10*(q1 - q0) - (6*qd0 + 4*qd1)*tscal;
    E = qd0*tscal; % as the t vector has been normalized
    F = q0;

    tt = [t.^5 t.^4 t.^3 t.^2 t ones(size(t))];
    c = [A B C zeros(size(A)) E F]';

    qt = tt*c;

    % compute optional velocity
    if nargout >= 2
        c = [ zeros(size(A)) 5*A 4*B 3*C  zeros(size(A)) E ]';
        qdt = tt*c/tscal;
    end

    % compute optional acceleration
    if nargout == 3
        c = [ zeros(size(A))  zeros(size(A)) 20*A 12*B 6*C  zeros(size(A))]';
        qddt = tt*c/tscal^2;
    end

author: petercorke

repo: robotics-toolbox-matlab


text:
%FIGURE_Critical_Angle

MoI = logspace(-3,0,100);
Mu = [0,0.05,0.1,0.5,1,5,10,inf];
P.m = 1;
P.g = 1;
P.L = 1;

Data = cell(length(Mu),0);
for i=1:length(Mu)
    Data(i).idxF = false(size(MoI));
    Data(i).idxB = false(size(MoI));
    Data(i).fail = false(size(MoI));
    Data(i).th = zeros(size(MoI));
    P.u = Mu(i);
    disp(['Running u = ' num2str(P.u)]);
    for j=1:length(MoI)
        P.I = MoI(j);
        C = ToppleFromRest(P);
        if ~isempty(C.th)
            if strcmp(C.exit,'SlipForwards')
               Data(i).idxF(j) = true;
            elseif strcmp(C.exit,'SlipBackwards')
                Data(i).idxB(j) = true;
            end
            Data(i).th(j) = C.th;
        else
            Data(i).fail(j) = true;
        end
    end
end

style = {'k--','b--','r--','m--','g--',...
    'k-','b-','r-','m-','g-'...
    'k:','b:','r:','m:','g:'};
LINEWIDTH = 3;
FontSize.Title = 16;
FontSize.label = 12;
names = cell(1,length(Mu));
for i=1:length(Mu)
    names{i} = num2str(['u = ' num2str(Mu(i))]);
end

figH = figure(400); clf;
set(figH,'Name','CriticalAngleLim','NumberTitle','off')
IDX = false(2,length(Mu));
for i=1:length(Mu)
    sty = style{mod(i-1,length(style))+1};
    idx = Data(i).idxB;  IDX(1,i) = sum(idx)~=0;
    subplot(2,1,1);  hold on;
    semilogx(MoI(idx),Data(i).th(idx)*180/pi,sty,'LineWidth',LINEWIDTH);
    idx = Data(i).idxF;  IDX(2,i) = sum(idx)~=0;
    subplot(2,1,2);  hold on;
    semilogx(MoI(idx),Data(i).th(idx)*180/pi,sty,'LineWidth',LINEWIDTH);
end

subplot(2,1,1);
 title('Backwards Slip (falling forward)','FontSize',FontSize.Title)
    xlabel('Moment of Inertia','FontSize',FontSize.label)
    ylabel('Critical Angle (deg)','FontSize',FontSize.label)
    legend(names(IDX(1,:)),'Location','NorthEast');
    set(gca,'Xscale','log')

    subplot(2,1,2);
 title('Forwards Slip (falling forward)','FontSize',FontSize.Title)
    xlabel('Moment of Inertia','FontSize',FontSize.label)
    ylabel('Critical Angle (deg)','FontSize',FontSize.label)
    legend(names(IDX(2,:)),'Location','NorthWest');
    set(gca,'Xscale','log')

save('DATA_Critical_Fall.mat','Data');
save2pdf('../WriteUp/Figures/Critical_Fall.pdf',figH,600);

%Later on, compute the acceleration of the tip at the onset of sliding.

author: MatthewPeterKelly

repo: dscTutorials


text:
function [varargout] = isrow(varargin)

% ISROW is a drop-in replacement for the same function that was
% introduced in MATLAB R2010b.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% try to automatically remove incorrect compat folders from the path, see https://github.com/fieldtrip/fieldtrip/issues/899

if isempty(strfind(mfilename('fullpath'), matlabroot))
  % automatic cleanup does not work when FieldTrip is installed inside the MATLAB path, see https://github.com/fieldtrip/fieldtrip/issues/1527

  alternatives = which(mfilename, '-all');
  if ~iscell(alternatives)
    % this is needed for octave, see https://github.com/fieldtrip/fieldtrip/pull/1171
    alternatives = {alternatives};
  end

  keep = true(size(alternatives));
  for i=1:numel(alternatives)
    keep(i) = keep(i) && ~any(alternatives{i}=='@');  % exclude methods from classes
    keep(i) = keep(i) && alternatives{i}(end)~='p';   % exclude precompiled files
  end
  alternatives = alternatives(keep);

  if exist(mfilename, 'builtin') || any(strncmp(alternatives, matlabroot, length(matlabroot)) & cellfun(@isempty, strfind(alternatives, fullfile('private', mfilename))))
    % remove this directory from the path
    p = fileparts(mfilename('fullpath'));
    warning('removing "%s" from your path, see http://bit.ly/2SPPjUS', p);
    rmpath(p);
    % call the original MATLAB function
    if exist(mfilename, 'builtin')
      [varargout{1:nargout}] = builtin(mfilename, varargin{:});
    else
      [varargout{1:nargout}] = feval(mfilename, varargin{:});
    end
    return
  end

end % automatic cleanup of compat directories

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is where the actual replacement code starts
% function tf = isrow(x)

% deal with the input arguments
if nargin==1
  [x] = deal(varargin{1:1});
else
  error('incorrect number of input arguments')
end

tf = length(size(x))==2 && size(x,1)==1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% deal with the output arguments

varargout = {tf};


text:
function [result coeff] = RS_PF_OLS_Wald(y,x1,x2)
% Calculates OLS Wald-Test
% H0: beta=0 in the model
% y = x1*beta+x2*gamma+eps with HAC consistent var-covar
% If want to test all coefficients, simply use (y,x1) and do not include x2

q = 0;
n = size(y,1);
p = size(x1,2);

if nargin>2;
    q = size(x2,2);
    x = [x1,x2];
    R = [eye(p),zeros(p,q)];
else
    x = x1;
    R = eye(p);
end;

coeff = ((inv(x'*x))*(x'*y));
nlag  = round(n^(1/4));

% Compute Newey-West adjusted heteroscedastic-serial consistent
% least-squares regression
nwresult   = bear.RS_PF_nwest(y,x,nlag);
varbetahat = nwresult.vcv;

result = (R*coeff)'/(R*varbetahat*R')*R*coeff;

author: european-central-bank

repo: BEAR-toolbox


text:
% pop_select() - given an input EEG dataset structure, output a new EEG data structure
%                retaining and/or excluding specified time/latency, data point, channel,
%                and/or epoch range(s).
% Usage:
%   >> OUTEEG = pop_select(INEEG, 'key1', value1, 'key2', value2 ...);
%
% Graphic interface:
%   "Time range" - [edit box] RETAIN only the indicated epoch latency or continuous data
%                  time range: [low high] in ms, inclusive. For continuous data, several
%                  time ranges may be specified, separated by semicolons.
%                  Example: "5 10; 12 EEG.xmax" will retain the indicated
%                  stretches of continuous data, and remove data portions outside
%                  the indicated ranges, e.g. from 0 s to 5 s and from 10 s to 12 s.
%                  Command line equivalent: 'time' (or 'notime' - see below)
%   "Time range" - [checkbox] EXCLUDE the indicated latency range(s) from the data.
%                  For epoched data, it is not possible to remove a range of latencies
%                  from the middle of the epoch, so either the low and/or the high values
%                  in the specified latency range (see above) must be at an epoch boundary
%                  (EEG.xmin, EEGxmax).  Command line equivalent: [if checked] 'notime'
%   "Point range" - [edit box] RETAIN the indicated data point range(s).
%                  Same options as for the "Time range" features (above).
%                  Command line equivalent: 'point' (or 'nopoint' - see below).
%   "Point range" - [checkbox] EXCLUDE the indicated point range(s).
%                  Command line equivalent: [if checked] 'nopoint'
%   "Epoch range" - [edit box] RETAIN the indicated data epoch indices in the dataset.
%                  This checkbox is only visible for epoched datasets.
%                  Command line equivalent: 'trial' (or 'notrial' - see below)
%   "Epoch range" - [checkbox] EXCLUDE the specified data epochs.
%                   Command line equivalent: [if checked] 'notrial'
%   "Channel range" - [edit box] RETAIN the indicated vector of data channels
%                  Command line equivalent: 'channel' (or 'nochannel' - see below)
%   "Channel range" - [checkbox] EXCLUDE the indicated channels.
%                  Command line equivalent: [if checked] 'nochannel'
%   "..." - [button] select channels by name.
%   "Scroll dataset" - [button] call the eegplot() function to scroll the
%                  channel activities in a new window for visual inspection.
%                  Commandline equivalent: eegplot() - see its help for details.
% Inputs:
%   INEEG         - input EEG dataset structure
%
% Optional inputs
%   'time'        - [min max] in seconds. Epoch latency or continuous data time range
%                   to retain in the new dataset, (Note: not ms, as in the GUI text entry
%                   above). For continuous data (only), several time ranges can be specified,
%                   separated by semicolons. Example: "5 10; 12 EEG.xmax" will retain
%                   the indicated times ranges, removing data  outside the indicated ranges
%                   e.g. here from 0 to 5 s and from 10 s to 12 s. (See also, 'notime')
%   'notime'      - [min max] in seconds. Epoch latency or continuous dataset time range
%                   to exclude from the new dataset. For continuous data, may be
%                   [min1 max1; min2 max2; ...] to exclude several time ranges. For epoched
%                   data, the latency range must include an epoch boundary, as latency
%                   ranges in the middle of epochs cannot be removed from epoched data.
%   'point'       - [min max] epoch or continuous data point range to retain in the new
%                   dataset. For continuous datasets, this may be [min1 max1; min2 max2; ...]
%                   to retain several point ranges. (Notes: If both 'point'/'nopoint' and
%                   'time' | 'notime' are specified, the 'point' limit values take precedence.
%                   The 'point' argument was originally a point vector, now deprecated).
%   'nopoint'     - [min max] epoch or continuous data point range to exclude in the new dataset.
%                   For epoched data, the point range must include either the first (0)
%                   or the last point (EEG.pnts), as a central point range cannot be removed.
%   'trial'       - array of trial indices to retain in the new dataset
%   'notrial'     - array of trial indices to exclude from the new dataset
%   'sorttrial'   - ['on'|'off'] sort trial indices before extracting them (default: 'on').
%   'channel'     - vector of channel indices to retain in the new
%                   dataset. Can also be a cell array of channel names.
%   'nochannel'   - vector of channel indices to exclude from the new
%                   dataset. Can also be a cell array of channel names.
%   'newname'     - name for the new dataset (OUTEEG)
%
% Outputs:
%   OUTEEG        - new EEG dataset structure
%
% Note: This function performs a conjunction (AND) of all its optional inputs.
%       Using negative counterparts of all options, any logical combination is
%       possible.
%
% Author: Arnaud Delorme, CNL/Salk Institute, 2001; SCCN/INC/UCSD, 2002-
%
% see also: eeglab()

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

% 01-25-02 reformated help & license -ad
% 01-26-02 changed the format for events and trial conditions -ad
% 02-04-02 changed display format and allow for negation of inputs -ad
% 02-17-02 removed the event removal -ad
% 03-17-02 added channel info subsets selection -ad
% 03-21-02 added event latency recalculation -ad

function [EEG, com] = pop_select( EEG, varargin);

com = '';
if nargin < 1
    help pop_select;
    return;
end;

if nargin < 2
   geometry = { [1 1 1] [1 1 0.25 0.23 0.51] [1 1 0.25 0.23 0.51] [1 1 0.25 0.23 0.51] ...
           [1 1 0.25 0.23 0.51] [1] [1 1 1]};
   uilist = { ...
         { 'Style', 'text', 'string', 'Select data in:', 'fontweight', 'bold'  }, ...
         { 'Style', 'text', 'string', 'Input desired range', 'fontweight', 'bold'  }, ...
         { 'Style', 'text', 'string', 'on->remove these', 'fontweight', 'bold'  }, ...
         { 'Style', 'text', 'string', 'Time range [min max] (s)', 'fontangle', fastif(length(EEG)>1, 'italic', 'normal') }, ...
         { 'Style', 'edit', 'string', '', 'enable', fastif(length(EEG)>1, 'off', 'on') }, ...
         { }, { 'Style', 'checkbox', 'string', '    ', 'enable', fastif(length(EEG)>1, 'off', 'on') },{ }, ...
         ...
         { 'Style', 'text', 'string', 'Point range (ex: [1 10])', 'fontangle', fastif(length(EEG)>1, 'italic', 'normal') }, ...
         { 'Style', 'edit', 'string', '', 'enable', fastif(length(EEG)>1, 'off', 'on') }, ...
         { }, { 'Style', 'checkbox', 'string', '    ', 'enable', fastif(length(EEG)>1, 'off', 'on') },{ }, ...
         ...
         { 'Style', 'text', 'string', 'Epoch range (ex: 3:2:10)', 'fontangle', fastif(length(EEG)>1, 'italic', 'normal') }, ...
         { 'Style', 'edit', 'string', '', 'enable', fastif(length(EEG)>1, 'off', 'on') }, ...
         { }, { 'Style', 'checkbox', 'string', '    ', 'enable', fastif(length(EEG)>1, 'off', 'on') },{ }, ...
         ...
         { 'Style', 'text', 'string', 'Channel range' }, ...
         { 'Style', 'edit', 'string', '', 'tag', 'chans' }, ...
         { }, { 'Style', 'checkbox', 'string', '    ' }, ...
         { 'style' 'pushbutton' 'string'  '...', 'enable' fastif(isempty(EEG.chanlocs), 'off', 'on') ...
           'callback' 'tmpchanlocs = EEG(1).chanlocs; [tmp tmpval] = pop_chansel({tmpchanlocs.labels}, ''withindex'', ''on''); set(findobj(gcbf, ''tag'', ''chans''), ''string'',tmpval); clear tmp tmpchanlocs tmpval' }, ...
           { }, { }, { 'Style', 'pushbutton', 'string', 'Scroll dataset', 'enable', fastif(length(EEG)>1, 'off', 'on'), 'callback', ...
                          'eegplot(EEG.data, ''srate'', EEG.srate, ''winlength'', 5, ''limits'', [EEG.xmin EEG.xmax]*1000, ''position'', [100 300 800 500], ''xgrid'', ''off'', ''eloc_file'', EEG.chanlocs);' } {}};
   results = inputgui( geometry, uilist, 'pophelp(''pop_select'');', 'Select data -- pop_select()' );
   if length(results) == 0, return; end;

   % decode inputs
   % -------------
   args = {};
   if ~isempty( results{1} )
       if ~results{2}, args = { args{:}, 'time', eval( [ '[' results{1} ']' ] ) };
       else            args = { args{:}, 'notime', eval( [ '[' results{1} ']' ] ) }; end;
   end;

   if ~isempty( results{3} )
       if ~results{4}, args = { args{:}, 'point', eval( [ '[' results{3} ']' ] ) };
       else            args = { args{:}, 'nopoint', eval( [ '[' results{3} ']' ] ) }; end;
   end;

   if ~isempty( results{5} )
       if ~results{6}, args = { args{:}, 'trial', eval( [ '[' results{5} ']' ] ) };
       else            args = { args{:}, 'notrial', eval( [ '[' results{5} ']' ] ) }; end;
   end;

   if ~isempty( results{7} )
       [ chaninds chanlist ] = eeg_decodechan(EEG.chanlocs, results{7});
       if isempty(chanlist), chanlist = chaninds; end;
       if ~results{8}, args = { args{:}, 'channel'  , chanlist };
       else            args = { args{:}, 'nochannel', chanlist }; end;
   end;

else
    args = varargin;
end;

%----------------------------AMICA---------------------------------
if isfield(EEG.etc,'amica') && isfield(EEG.etc.amica,'prob_added')
    for index = 1:2:length(args)
       if strcmpi(args{index}, 'channel')
           args{index+1} = [ args{index+1} EEG.nbchan-(0:2*EEG.etc.amica.num_models-1)];

       end;

    end;
end;
%--------------------------------------------------------------------

% process multiple datasets
% -------------------------
if length(EEG) > 1
    [ EEG com ] = eeg_eval( 'pop_select', EEG, 'warning', 'on', 'params', args);
    return;
end;

if isempty(EEG.chanlocs), chanlist = [1:EEG.nbchan];
else                      chanlocs = EEG.chanlocs; chanlist = { chanlocs.labels };
end;
g = hlp_microcache('popselect',@finputcheck,args, { 'time'    'real'      []         []; ...
                        'notime'  'real'      []         []; ...
                        'trial'   'integer'   []         [1:EEG.trials]; ...
                        'notrial' 'integer'   []         []; ...
                        'point'   'integer'   []         []; ...
                        'nopoint' 'integer'   []         []; ...
                        'channel'   { 'integer' 'cell' }  []  chanlist;
                        'nochannel' { 'integer' 'cell' }   []  [];
                        'trialcond'   'integer'   []         []; ...
                        'notrialcond' 'integer'   []         []; ...
                        'sorttrial'   'string'    { 'on' 'off' } 'on' }, 'pop_select');
if isstr(g), error(g); end;

if strcmpi(g.sorttrial, 'on')
    g.trial = sort(setdiff( g.trial, g.notrial ));
else
    g.trial(ismember(g.trial,g.notrial)) = [];
    % still warn about & remove duplicate trials (may be removed in the future)
    [p,q] = unique(g.trial);
    if length(p) ~= length(g.trial)
        disp('Warning: trial selection contained duplicated elements, which were removed.');
    end
    g.trial = g.trial(sort(q));
end

if isempty(g.channel) && ~iscell(g.nochannel) && ~iscell(chanlist)
    g.channel = [1:EEG.nbchan];
end;

if iscell(g.channel) && ~iscell(g.nochannel) && ~isempty(EEG.chanlocs)
     noChannelAsCell = {};
     for nochanId = 1:length(g.nochannel)
         noChannelAsCell{nochanId} = EEG.chanlocs(g.nochannel(nochanId)).labels;
     end;
     g.nochannel =   noChannelAsCell;
end;

if strcmpi(g.sorttrial, 'on')
    g.channel = sort(setdiff( lower(g.channel), lower(g.nochannel) ));
else
    g.channel(ismember(lower(g.channel),lower(g.nochannel))) = [];
    % still warn about & remove duplicate channels (may be removed in the future)
    [p,q] = unique(g.channel);
    if length(p) ~= length(g.channel)
        disp('Warning: channel selection contained duplicated elements, which were removed.');
    end
    g.channel = g.channel(sort(q));
end

if ~isempty(EEG.chanlocs)
    if strcmpi(g.sorttrial, 'on')
        g.channel = eeg_decodechan(EEG.chanlocs, g.channel);
    else
        % we have to protect the channel order against changes by eeg_decodechan
        if iscell(g.channel)
            % translate channel names into indices
            [inds,names] = eeg_decodechan(EEG.chanlocs, g.channel);
            % and sort the indices back into the original order of channel names
            [tmp,I] = ismember(lower(g.channel),lower(names));
            g.channel = inds(I);
        end
    end
end;

if ~isempty(g.time) & (g.time(1) < EEG.xmin*1000) & (g.time(2) > EEG.xmax*1000)
   error('Wrong time range');
end;
if min(g.trial) < 1 | max( g.trial ) > EEG.trials
   error('Wrong trial range');
end;
if min(g.channel) < 1 | max( g.channel ) > EEG.nbchan
   error('Wrong channel range');
end;

if size(g.point,2) > 2,
    g.point = [g.point(1) g.point(end)];
    disp('Warning: vector format for point range is deprecated');
end;
if size(g.nopoint,2) > 2,
    g.nopoint = [g.nopoint(1) g.nopoint(end)];
    disp('Warning: vector format for point range is deprecated');
end;
if ~isempty( g.point )
    g.time = zeros(size(g.point));
    for index = 1:length(g.point(:))
        g.time(index) = eeg_point2lat(g.point(index), 1, EEG.srate, [EEG.xmin EEG.xmax]);
    end;
    g.notime = [];
end;
if ~isempty( g.nopoint )
    g.notime = zeros(size(g.nopoint));
    for index = 1:length(g.nopoint(:))
        g.notime(index) = eeg_point2lat(g.nopoint(index), 1, EEG.srate, [EEG.xmin EEG.xmax]);
    end;
    g.time = [];
end;
if ~isempty( g.notime )
    if size(g.notime,2) ~= 2
        error('Time/point range must contain 2 columns exactly');
    end;
    if g.notime(2) == EEG.xmax
        g.time = [EEG.xmin g.notime(1)];
    else
        if g.notime(1) == EEG.xmin
            g.time = [g.notime(2) EEG.xmax];
        elseif EEG.trials > 1
            error('Wrong notime range. Remember that it is not possible to remove a slice of time for data epochs.');
        end;
    end;
    if floor(max(g.notime(:))) > EEG.xmax || min(g.notime(:)) < EEG.xmin
        error('Time/point range out of data limits');
    end;
end;
if ~isempty(g.time)
    if size(g.time,2) ~= 2
        error('Time/point range must contain 2 columns exactly');
    end;
    for index = 1:length(g.time)
        if g.time(index) > EEG.xmax
            g.time(index) = EEG.xmax;
            disp('Upper time limits exceed data, corrected');
        elseif g.time(index) < EEG.xmin
            g.time(index) = EEG.xmin;
            disp('Lower time limits exceed data, corrected');
        end;
    end;
end;

% select trial values
%--------------------
if ~isempty(g.trialcond)
   try, tt = struct( g.trialcond{:} ); catch
      error('Trial conditions format error');
   end;
   ttfields = fieldnames (tt);
   for index = 1:length(ttfields)
        if ~isfield( EEG.epoch, ttfields{index} )
            error([ ttfields{index} 'is not a field of EEG.epoch' ]);
        end;
        tmpepoch = EEG.epoch;
	    eval( [ 'Itriallow  = find( [ tmpepoch(:).' ttfields{index} ' ] >= tt.' ttfields{index} '(1) );' ] );
	    eval( [ 'Itrialhigh = find( [ tmpepoch(:).' ttfields{index} ' ] <= tt.' ttfields{index} '(end) );' ] );
	    Itrialtmp = intersect(Itriallow, Itrialhigh);
	    g.trial = intersect( g.trial(:)', Itrialtmp(:)');
   end;
end;

if isempty(g.trial)
   error('Empty dataset, no trial');
end;
if length(g.trial) ~= EEG.trials
	fprintf('Removing %d trial(s)...\n', EEG.trials - length(g.trial));
end;
if length(g.channel) ~= EEG.nbchan
	fprintf('Removing %d channel(s)...\n', EEG.nbchan - length(g.channel));
end;

% For AMICA probabilities...
%-----------------------------------------------------
if isfield(EEG.etc, 'amica') && ~isempty(EEG.etc.amica) && isfield(EEG.etc.amica, 'v_smooth') && ~isempty(EEG.etc.amica.v_smooth) && ~isfield(EEG.etc.amica,'prob_added')
    if isfield(EEG.etc.amica, 'num_models') && ~isempty(EEG.etc.amica.num_models)
        if size(EEG.data,2) == size(EEG.etc.amica.v_smooth,2) && size(EEG.data,3) == size(EEG.etc.amica.v_smooth,3) && size(EEG.etc.amica.v_smooth,1) == EEG.etc.amica.num_models

            EEG = eeg_formatamica(EEG);

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

            [EEG com] = pop_select(EEG,args{:});

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

            EEG = eeg_reformatamica(EEG);
            EEG = eeg_checkamica(EEG);
            return;
        else
            disp('AMICA probabilities not compatible with size of data, probabilities cannot be rejected')

            disp('Resuming rejection...')
        end
    end

end
% ------------------------------------------------------

% recompute latency and epoch number for events
% ---------------------------------------------
if length(g.trial) ~= EEG.trials & ~isempty(EEG.event)
    if ~isfield(EEG.event, 'epoch')
        disp('Pop_epoch warning: bad event format with epoch dataset, removing events');
        EEG.event = [];
    else
		if isfield(EEG.event, 'epoch')
			keepevent = [];
			for indexevent = 1:length(EEG.event)
				newindex = find( EEG.event(indexevent).epoch == g.trial );
				if ~isempty(newindex)
					keepevent = [keepevent indexevent];
					if isfield(EEG.event, 'latency')
						EEG.event(indexevent).latency = EEG.event(indexevent).latency - (EEG.event(indexevent).epoch-1)*EEG.pnts + (newindex-1)*EEG.pnts;
					end;
					EEG.event(indexevent).epoch = newindex;
				end;
			end;
            diffevent = setdiff([1:length(EEG.event)], keepevent);
			if ~isempty(diffevent)
				disp(['Pop_select: removing ' int2str(length(diffevent)) ' unreferenced events']);
				EEG.event(diffevent) = [];
			end;
		end;
    end;
end;

% performing removal
% ------------------
if ~isempty(g.time) | ~isempty(g.notime)
    if EEG.trials > 1
        % select new time window
        % ----------------------
        try,   tmpevent = EEG.event;
               tmpeventlatency = [ tmpevent.latency ];
        catch, tmpeventlatency = [];
        end;
        alllatencies = 1-(EEG.xmin*EEG.srate); % time 0 point
        alllatencies = linspace( alllatencies, EEG.pnts*(EEG.trials-1)+alllatencies, EEG.trials);
        [EEG.data tmptime indices epochevent]= epoch(EEG.data, alllatencies, ...
                                                     [g.time(1) g.time(2)]*EEG.srate, 'allevents', tmpeventlatency);
        tmptime = tmptime/EEG.srate;
        if g.time(1) ~= tmptime(1) & g.time(2)-1/EEG.srate ~= tmptime(2)
            fprintf('pop_select(): time limits have been adjusted to [%3.3f %3.3f] to fit data points limits\n', tmptime(1), tmptime(2)+1/EEG.srate);
        end;
        EEG.xmin = tmptime(1);
        EEG.xmax = tmptime(2);
        EEG.pnts = size(EEG.data,2);
        alllatencies = alllatencies(indices);

        % modify the event structure accordingly (latencies and add epoch field)
        % ----------------------------------------------------------------------
        allevents = [];
        newevent = [];
        count = 1;
        if ~isempty(epochevent)
            newevent = EEG.event(1);
            for index=1:EEG.trials
                for indexevent = epochevent{index}
                    newevent(count)         = EEG.event(indexevent);
                    newevent(count).epoch   = index;
                    newevent(count).latency = newevent(count).latency - alllatencies(index) - tmptime(1)*EEG.srate + 1 + EEG.pnts*(index-1);
                    count = count + 1;
                end;
            end;
        end;
        EEG.event = newevent;

        % erase event-related fields from the epochs
        % ------------------------------------------
        if ~isempty(EEG.epoch)
            fn = fieldnames(EEG.epoch);
            EEG.epoch = rmfield(EEG.epoch,{fn{strmatch('event',fn)}});
        end;
    else
        if isempty(g.notime)
            if length(g.time) == 2 && EEG.xmin < 0
                disp('Warning: negative minimum time; unchanged to ensure correct latency of initial boundary event');
            end;
            g.notime = g.time';
            g.notime = g.notime(:);
            if g.notime(1) ~= 0, g.notime = [EEG.xmin g.notime(:)'];
            else                 g.notime = [g.notime(2:end)'];
            end;
            if g.time(end) == EEG.xmax, g.notime(end) = [];
            else                        g.notime(end+1) = EEG.xmax;
            end;

            for index = 1:length(g.notime)
                if g.notime(index) ~= 0  & g.notime(index) ~= EEG.xmax
                    if mod(index,2), g.notime(index) = g.notime(index) + 1/EEG.srate;
                    else             g.notime(index) = g.notime(index) - 1/EEG.srate;
                    end;
                end;
            end;
            g.notime = reshape(g.notime, 2, length(g.notime)/2)';
        end;

        nbtimes = length(g.notime(:));
        points = eeg_lat2point(g.notime(:)', ones(1,nbtimes), EEG.srate, [EEG.xmin EEG.xmax]);
        points = reshape(points, size(g.notime));
        EEG = eeg_eegrej(EEG, points);
    end
end;

% performing removal
% ------------------
if ~isequal(g.channel,1:size(EEG.data,1)) || ~isequal(g.trial,1:size(EEG.data,3))
    EEG.data  = EEG.data(g.channel, :, g.trial);
    if isfield(EEG,'srcpot') && ~isempty(EEG.srcpot)
        EEG.srcpot = EEG.srcpot(:, :, g.trial);  % TODO: this is an ugly hack
    end
end
if ~isempty(EEG.icaact), EEG.icaact = EEG.icaact(:,:,g.trial); end;
EEG.trials    = length(g.trial);
EEG.pnts      = size(EEG.data,2);
EEG.nbchan    = length(g.channel);
if ~isempty(EEG.chanlocs)
    EEG.chanlocs = EEG.chanlocs(g.channel);
end;
if ~isempty(EEG.epoch)
   EEG.epoch = EEG.epoch( g.trial );
end;
if ~isempty(EEG.specdata)
	if length(g.point) == EEG.pnts
   		EEG.specdata = EEG.specdata(g.channel, :, g.trial);
   	else
   		EEG.specdata = [];
   		fprintf('Warning: spectral data were removed because of the change in the numner of points\n');
   	end;
end;

% ica specific
% ------------
if ~isempty(EEG.icachansind)

    rmchans = setdiff( EEG.icachansind, g.channel ); % channels to remove

    % channel sub-indices
    % -------------------
    icachans = 1:length(EEG.icachansind);
    for index = length(rmchans):-1:1
        chanind           = find(EEG.icachansind == rmchans(index));
        icachans(chanind) = [];
    end;

    % new channels indices
    % --------------------
    count   = 1;
    newinds = [];
    for index = 1:length(g.channel)
        if any(EEG.icachansind == g.channel(index))
            newinds(count) = index;
            count          = count+1;
        end;
    end;
    EEG.icachansind = newinds;

else
    icachans = 1:size(EEG.icasphere,2);
end;

if ~isempty(EEG.icasphere)
    EEG.icasphere = EEG.icasphere(:,icachans);
end;
if ~isempty(EEG.icawinv)
    EEG.icawinv = EEG.icawinv(icachans,:);
end;
if ~isempty(EEG.specicaact)
    if length(g.point) == EEG.pnts
        EEG.specicaact = EEG.specicaact(icachans, :, g.trial);
    else
        EEG.specicaact = [];
        fprintf('Warning: spectral ICA data were removed because of the change in the numner of points\n');
    end;
end;

% for stats, can adapt remove the selected trials and electrodes
% in the future to gain time -----------------------------------
EEG.reject = struct('rejmanual',{[]});
EEG.stats = struct('jp',{[]});
EEG = eeg_checkset(EEG, 'eventconsistency');

% generate command
% ----------------
if nargout > 1
com = sprintf('EEG = pop_select( %s,%s);', inputname(1), vararg2str(args));
end

return;

% ********* OLD, do not remove any event any more
% ********* in the future maybe do a pack event to remove events not in the time range of any epoch

if ~isempty(EEG.event)
    % go to array format if necessary
    if isstruct(EEG.event), format = 'struct';
    else                     format = 'array';
    end;
    switch format, case 'struct', EEG = eventsformat(EEG, 'array'); end;

    % keep only events related to the selected trials
    Indexes = [];
    Ievent  = [];
    for index = 1:length( g.trial )
        currentevents = find( EEG.event(:,2) == g.trial(index));
        Indexes = [ Indexes ones(1, length(currentevents))*index ];
        Ievent  = union( Ievent, currentevents );
    end;
    EEG.event = EEG.event( Ievent,: );
    EEG.event(:,2) = Indexes(:);

    switch format, case 'struct', EEG = eventsformat(EEG, 'struct'); end;
end;

author: PatternRecognition

repo: OpenBMI


text:
function build_mmx(verbose)
% BUILD_MMX - compiles mmx() for different platforms and provides help
%            regarding compilation.
%
%  BUILD_MMX will try to compile, in this order, 3 different builds of mmx:
%  mmx_mkl_single    - linked to Intel's single-threaded MKL library (usually fastest)
%  mmx_mkl_multi     - linked to the multithreaded BLAS/LAPACK libraries that come
%                      with Matlab.
%  mmx_naive         - does not link to anything, uses simple C-loops.
%
%  The first time BUILD_MMX succeeds, it will compile again to 'mmx', so
%  that the mex-file mmx should be the fastest possible build on your
%  system.
%
%  BUILD_MMX has been tested on Win32, Win64, OSX, Linux 64
%

% %% FOR LINUX OR MAC SYSTEMS:
%
% To properly link to Intel's MKL, user needs to repackage their libraries
% into one single statically linked library. The instructions are as
% follows:
%
%
% Download Intel MKL for Linux here:
% http://software.intel.com/en-us/articles/non-commercial-software-download/
%
% Donwload Intel MKL for Mac here:
% https://registrationcenter.intel.com/RegCenter/AutoGen.aspx?ProductID=1518&AccountID=&EmailID=&ProgramID=&RequestDt=&rm=EVAL&lang=
%
% The Default installation directory for both Linux and Mac will be
% /opt/intel/
% with the MKL libraries in /opt/intel/mkl
%
% %% To build needed static Library
%    assuming default installation directory
%
% Run the following commands in Linux/Mac terminal:
%
% sudo -s
% cd /opt/intel/mkl/tools/builder
% cat blas_example_list > blas_lapack_list
% cat lapack_example_list >> blas_lapack_list
%
% For Linux 64 bit:
% make libintel64 interface=ilp64 export=blas_lapack_list name=libsingle_mkl_ilp64 threading=sequential
% For Linux 32 bit:
% make libia32 interface=lp64 export=blas_lapack_list name=libsingle_mkl_32 threading=sequential
%
% For Mac:
% make libuni interface=ilp64 export=blas_lapack_list name=libsingle_mkl_ilp64 threading=sequential
%
% A new libsingle_mkl_ilp64.so, libsingle_mkl_32.so, or
% libsingle_mkl_ilp64.dylib will appear.
% This needs to be copied to Matlab's external libraries directory.
%
% For Mac:
% cp libsingle_mkl_ilp64* MATLAB_ROOT/extern/lib/maci64
%
% For Linux 64 bit:
% cp libsingle_mkl_ilp64* MATLAB_ROOT/extern/lib/glnxa64
% For Linux 32 bit:
% cp libsingle_mkl_32* MATLAB_ROOT/extern/lib/glnx86
%
% Where MATLAB_ROOT is the installation directory of your Matlab.

if nargin == 0
   verbose = false;
end

clc

build_names  = {'mmx_mkl_single', 'mmx_mkl_multi','mmx_naive'};

built_mmx   = false;

arch        = computer('arch');

for b = 1:3
   name = build_names{b};

   [link, define]  = deal({});
   [inc_dir, link_dir, Cflags, Lflags]  = deal('');

   switch arch
      case {'win64','win32'}
         switch name
            case 'mmx_naive'
               define   = {'WIN_SYSTEM'};

            case 'mmx_mkl_multi'
               root     = matlabroot;
               if strcmp(arch,'win32')
                  inc_dir  = [root '\extern\lib\win32\microsoft'];
               else
                  inc_dir  = [root '\extern\lib\win64\microsoft'];
               end
               link     = {'libmwblas','libmwlapack'};
               define   = {'WIN_SYSTEM','USE_BLAS'};

            case 'mmx_mkl_single'
               root     = 'C:\Program Files (x86)\Intel\Composer XE 2011 SP1\mkl';
               inc_dir  = [root '\include'];
               if strcmp(arch,'win32')
                  link_dir  = [root '\lib\ia32'];
                  link     = {'mkl_intel_c','mkl_sequential','mkl_core'};
                  define   = {'WIN_SYSTEM','USE_BLAS','MKL_32'};
               else
                  link_dir  = [root '\lib\intel64'];
                  link     = {'mkl_intel_ilp64','mkl_sequential','mkl_core'};
                  define   = {'WIN_SYSTEM','USE_BLAS','MKL_ILP64'};
               end
         end
      case {'glnxa64','glnx86'}
         switch name
            case 'mmx_naive'
               link     = {'pthread'};
               define   = {'UNIX_SYSTEM'};
            case 'mmx_mkl_multi'
               if strcmp(arch,'glnx86')
               inc_dir  = [matlabroot '/extern/lib/glnx86'];
               else
               inc_dir  = [matlabroot '/extern/lib/glnxa64'];
               end
               link     = {'mwblas','mwlapack','pthread'};
               define   = {'UNIX_SYSTEM','USE_BLAS'};
            case 'mmx_mkl_single'
               root = '/opt/intel/mkl';
               inc_dir  = [ root '/include'];
               if strcmp(arch,'glnx86')
                link_dir  = [matlabroot '/extern/lib/glnx86'];
                link     = {'single_mkl_32','pthread'};
                define   = {'UNIX_SYSTEM', 'USE_BLAS', 'MKL_32'};
               else
                link_dir  = [matlabroot '/extern/lib/glnxa64'];
                link     = {'small_mkl_ilp64','pthread'};
                define   = {'UNIX_SYSTEM', 'USE_BLAS', 'MKL_ILP64'};
               end
         end
      case {'maci64'}
         switch name
            case 'mmx_naive'
               link     = {'pthread'};
               define   = {'UNIX_SYSTEM'};

            case 'mmx_mkl_multi'
               root     = matlabroot;
               inc_dir  = [root '/extern/lib/maci64'];
               link     = {'mwblas','mwlapack','pthread'};
               define   = {'UNIX_SYSTEM','USE_BLAS'};

            case 'mmx_mkl_single'
               root     = '/opt/intel/mkl';
               inc_dir  = [ root '/include'];
               link_dir = [matlabroot '/extern/lib/maci64'];
               link     = {'single_mkl_ilp64','pthread'};
               %link     = {'small_mkl_ilp64','pthread'};
               define   = {'UNIX_SYSTEM', 'USE_BLAS', 'MKL_ILP64'};
         end

      otherwise
         error unsupported_architecture
   end

   if ~isempty(link_dir)
      if strcmp(arch,'glnxa64') || strcmp(arch,'maci64')
         L_dir  = {['LDFLAGS="\$LDFLAGS -L' link_dir  ' ' Lflags '"']};
      else
         L_dir  = {['-L' link_dir]};
      end
   else
      L_dir  = {};
   end

   if ~isempty(inc_dir)
      if strcmp(arch,'glnxa64') || strcmp(arch,'maci64')
         I_dir  = {['CXXFLAGS="\$CXXFLAGS -I' inc_dir ' ' Cflags '"']};
      else
         I_dir  = {['-I' inc_dir]};
      end
   else
      I_dir  = {};
   end

   prefix   = @(pref,str_array) cellfun(@(x)[pref x],str_array,'UniformOutput',0);
   l_link   = prefix('-l',link);
   D_define = prefix('-D',define);

   if verbose
      verb  = {'-v'};
   else
      verb  = {};
   end

   try
      check_dir(link_dir, link)
      check_dir(inc_dir)
      clear(name)
      command = {verb{:}, I_dir{:}, L_dir{:}, l_link{:}, D_define{:}}; %#ok<*CCAT>
      fprintf('==========\nTrying to compile ''%s'', using \n',name);
      fprintf('%s, ',command{:})
      fprintf('\n')
      mex(command{:}, '-output', name, 'mmx.cpp');
      fprintf('Compilation of ''%s'' succeeded.\n',name);
      if ~built_mmx
         fprintf('Compiling again to ''mmx'' target using ''%s'' build.\n',name);
         mex(command{:}, '-output','mmx','mmx.cpp');
         built_mmx = true;
      end
   catch err
      fprintf('Compilation of ''%s'' failed with error:\n%s\n',name,err.message);
   end
end

function check_dir(dir,files)
if ~isempty(dir)
   here = cd(dir);
   if nargin == 2
      for i = 1:size(files)
         if isempty(ls(['*' files{i} '.*']))
            cd(here);
            error('could not find file %s', files{i});
         end
      end
   end
   cd(here);
end

author: NicolasBoumal

repo: manopt


text:
function [varargout] = plgndr(varargin)

% PLGNDR associated Legendre function
%
% y = plgndr(n,k,x) computes the values of the associated Legendre functions
% of degree N and order K
%
% implemented as MEX file

% the original implementation was based on "Numerical Recipes in C", version 2.0
% but has been replaced with an equvalent function from GNU Scientific Library

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

% compile the missing mex file on the fly
% remember the original working directory
pwdir = pwd;
pwdir_ressetter=onCleanup(@()cd(pwdir));

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc  = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

mexfullpath=fullfile(mexdir,sprintf('%s.%s',mexname,mexext()));
disp(mexfullpath);
has_mex_func=@()exist(mexfullpath,'file');

if ~has_mex_func()
  % try to compile the mex file on the fly
  ft_warning('trying to compile MEX file from %s', mexsrc);
  cd(mexdir);
  mex(mexsrc);

  if ~has_mex_func()
    ft_error('could not locate / compile MEX file for %s in %s', ...
                mexname, mexfullpath);
  end
end

% execute the mex file that was just created
funname   = mfilename;
funhandle = str2func(funname);
[varargout{1:nargout}] = funhandle(varargin{:});


text:
function cvt_movie ( )

%*****************************************************************************80
%
%% CVT_MOVIE generates and animates a CVT dataset.
%
%  Discussion:
%
%    This program is meant to be used interactively.  It's also
%    possible to prepare a simple input file beforehand and use it
%    in batch mode.
%
%    The program requests input values from the user:
%
%    * NDIM, the spatial dimension,
%    * N, the number of points to generate,
%    * SEED, a seed to use for random number generation;
%    * INIT, initialize the points:
%      ** file, by reading data from file;
%      ** 'GRID', picking points from a grid;
%      ** 'HALTON', from a Halton sequence;
%      ** 'RAND', using MATLAB's RAND function;
%      ** 'UNIFORM', using a simple uniform RNG;
%    * IT_MAX, the maximum number of iterations;
%    * IT_FIXED, the number of iterative steps to take
%      using a fixed set of sampling points.
%    * SAMPLE, how to conduct the sampling:
%      ** 'GRID', picking points from a grid;
%      ** 'HALTON', from a Halton sequence;
%      ** 'RAND', using MATLAB's RAND function;
%      ** 'UNIFORM', using a simple uniform RNG;
%    * SAMPLE_NUM, the number of sampling points;
%    * BATCH, the number of sampling points to create at one time;
%    * MOVIE_NAME is the name of the file in which the movie is stored;\n' );
%    * OUTPUT, a file into which the data can be stored.
%
%    To indicate that no further computations are desired, it is
%    enough to input a nonsensical value, such as -1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Qiang Du, Vance Faber, and Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review, Volume 41, 1999, pages 637-676.
%
  DEBUG = 1;

  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CVT_MOVIE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Generate and animate a CVT dataset.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program is meant to be used interactively.\n' );
  fprintf ( 1, '  It is also possible to prepare a simple input \n' );
  fprintf ( 1, '  file beforehand and use it in batch mode.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program requests input values from the user:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * NDIM, the spatial dimension,\n' );
  fprintf ( 1, '  * N, the number of points to generate,\n' );
  fprintf ( 1, '  * SEED, a seed to use for random number generation,\n' );
  fprintf ( 1, '  * INIT, initialize the points:\n' );
  fprintf ( 1, '    ** file, read data from a file;\n' );
  fprintf ( 1, '    ** ''GRID'', by picking points from a grid;\n' );
  fprintf ( 1, '    ** ''HALTON'', from a Halton sequence;\n' );
  fprintf ( 1, '    ** ''RAND'', using MATLAB''s RAND function;\n' );
  fprintf ( 1, '    ** ''UNIFORM'', using a simple uniform RNG;\n' );
  fprintf ( 1, '  * IT_MAX, the maximum number of iterations.\n' );
  fprintf ( 1, '  * IT_FIXED, the number of iterative steps to take\n' );
  fprintf ( 1, '    using a fixed set of sampling points.\n' );
  fprintf ( 1, '  * SAMPLE, how to conduct the sampling.\n' );
  fprintf ( 1, '    ** ''GRID'', by picking points from a grid;\n' );
  fprintf ( 1, '    ** ''HALTON'', from a Halton sequence;\n' );
  fprintf ( 1, '    ** ''RAND'', using MATLAB''s RAND function;\n' );
  fprintf ( 1, '    ** ''UNIFORM'', using a simple uniform RNG;\n' );
  fprintf ( 1, '  * SAMPLE_NUM, the number of sample points;\n' );
  fprintf ( 1, '  * BATCH, the number of sampling points to create at one time;\n' );
  fprintf ( 1, '  * OUTPUT, a file into which the data is stored.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  To indicate that no further computations are \n' );
  fprintf ( 1, '  desired, it is enough to input a nonsensical value, \n' );
  fprintf ( 1, '  such as -1.\n' );

  while ( 1 )

    fprintf ( 1, '  *\n' );
    fprintf ( 1, ' *\n' );
    fprintf ( 1, '*  Ready to generate a new dataset:\n' );
    fprintf ( 1, ' *\n' );
    fprintf ( 1, '  *\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  NDIM is the spatial dimension.\n' );
    fprintf ( 1, '  (Try ''2'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    ndim = [];
    ndim = input ( '  Enter NDIM:  ' );

    fprintf ( 1, '  Default NDIM = %12d\n', ndim );

    if ( ndim < 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of NDIM = %12d\n', ndim );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  N is the number of points to generate.\n' );
    fprintf ( 1, '  (Try ''25'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    n = [];
    n = input ( '  Enter N:  ' );

    fprintf ( 1, '  User input N = %12d\n', n );

    if ( n < 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of N = %12d\n', n );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  SEED is a seed for the random number generation.\n' );
%   fprintf ( 1, '  (Try ''123456789'' if you have no preference.)\n' );
%   fprintf ( 1, '  (Any value less than 0 terminates execution.)\n' );
%   seed = [];
%   seed = input ( '  Enter SEED:  ' );
    seed = 123456789;
    fprintf ( 1, '  Default SEED = %d\n', seed );

    if ( seed < 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of SEED = %12d\n', seed );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  INIT is the method of initializing the data:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  file       read data from a file;\n' );
    fprintf ( 1, '  ''GRID''     by picking points from a grid;\n' );
    fprintf ( 1, '  ''HALTON''   from a Halton sequence;\n' );
    fprintf ( 1, '  ''RAND''     using MATLAB''s RAND function;\n' );
    fprintf ( 1, '  ''UNIFORM''  using a simple uniform RNG;\n' );
    fprintf ( 1, ' \n' );
    fprintf ( 1, '  (Try ''RAND'' if you have no preference.)\n' );
    fprintf ( 1, '  (A blank value terminates execution).\n' );
    fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
    fprintf ( 1, ' \n' );

    init_string = [];
    init_string = input ( '  Enter INIT:  ' );

    fprintf ( 1, '  User input INIT = "%s".\n', init_string );

    if ( s_len_trim ( init_string ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of INIT \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    if ( s_eqi ( init_string, 'RAND'  ) )
      init = -1;
    elseif ( s_eqi ( init_string, 'RANDOM' ) )
      init_string = 'RAND';
      init = -1;
    elseif ( s_eqi ( init_string, 'UNIFORM' ) )
      init = 0;
    elseif ( s_eqi ( init_string, 'HALTON'  ) )
      init = 1;
    elseif ( s_eqi ( init_string, 'GRID'    ) )
      init = 2;
    elseif ( 0 < s_len_trim ( init_string ) )
      init = 3;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of INIT \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  IT_MAX is the maximum number of iterations.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  An iteration carries out the following steps:\n' );
    fprintf ( 1, '  * the Voronoi region associated with each\n' );
    fprintf ( 1, '    generator is estimated by sampling;\n' );
    fprintf ( 1, '  * the centroid of each Voronoi region is estimated.\n' );
    fprintf ( 1, '  * the generator is replaced by the centroid.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  If "enough" sampling points are used,\n' );
    fprintf ( 1, '  and "enough" iterations are taken, this process\n' );
    fprintf ( 1, '  will converge.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''50'' if you have no preference.)\n' );
    fprintf ( 1, '  (A negative value terminates execution).\n' );
    fprintf ( 1, '\n' );
    it_max = [];
    it_max = input ( '  Enter IT_MAX:  ' );

    fprintf ( 1, '  User input IT_MAX = %12d\n', it_max );

    if ( it_max < 0 )
      fprintf ( 1, ' \n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of IT_MAX = %12d\n', it_max );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  IT_FIXED is the number of consecutive iterations\n' );
    fprintf ( 1, '  to take with a fixed set of sample points.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Setting IT_FIXED to 1 means a new set of sample\n' );
    fprintf ( 1, '  points is generated on every iterative step;\n' );
    fprintf ( 1, '  Setting IT_FIXED equal to IT_MAX means a single set\n' );
    fprintf ( 1, '  of sample points is used for the entire iteration.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Any value between 1 and IT_MAX is reasonable.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try "%d" if you do not have a preference).\n', it_max );
    fprintf ( 1, '  (A 0 or negative value terminates execution).\n' );
    fprintf ( 1, '\n' );
    it_fixed = [];
    it_fixed = input ( '  Enter IT_FIXED:  ' );

    fprintf ( 1, '  User input IT_FIXED = %12d\n', it_fixed );

    if ( it_max < 0 )
      fprintf ( 1, ' \n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of IT_FIXED = %12d\n', it_fixed );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SAMPLE is the method of sampling the region:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  ''GRID''     by picking points from a grid;\n' );
    fprintf ( 1, '  ''HALTON''   from a Halton sequence;\n' );
    fprintf ( 1, '  ''RAND''     using MATLAB''s RAND function;\n' );
    fprintf ( 1, '  ''UNIFORM''  using a simple uniform RNG;\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''RAND'' if you have no preference.)\n' );
    fprintf ( 1, '  (A blank value terminates execution).\n' );
    fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
    fprintf ( 1, '\n' );

    sample_string = [];
    sample_string = input ( '  Enter SAMPLE:  ' );

    fprintf ( 1, '  User input SAMPLE = "%s".\n', sample_string );

    if ( s_len_trim ( sample_string ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of SAMPLE \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    if ( s_eqi ( sample_string, 'RAND'  ) )
      sample = -1;
    elseif ( s_eqi ( sample_string, 'RANDOM' ) )
      sample = -1;
      sample_string = 'RAND';
    elseif ( s_eqi ( sample_string, 'UNIFORM' ) )
      sample = 0;
    elseif ( s_eqi ( sample_string, 'HALTON'  ) )
      sample = 1;
    elseif ( s_eqi ( sample_string, 'GRID'    ) )
      sample = 2;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of SAMPLE \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SAMPLE_NUM is the number of sample points.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The Voronoi regions will be explored by generating\n' );
    fprintf ( 1, '  SAMPLE_NUM points.  For each sample point, the\n' );
    fprintf ( 1, '  nearest generator is found.  Using more points\n' );
    fprintf ( 1, '  gives a better estimate of these regions.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  SAMPLE_NUM should be much larger than N, the\n' );
    fprintf ( 1, '  number of generators.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''10000'' if you have no preference.)\n' );
    fprintf ( 1, '  (A zero or negative value terminates execution.)\n' );
    fprintf ( 1, '\n' );

    sample_num = [];
    sample_num = input ( '  Enter SAMPLE_NUM:  ' );

    fprintf ( 1, '  User input SAMPLE_NUM = %12d\n', sample_num );

    if ( sample_num <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of SAMPLE_NUM = %12d\n', sample_num );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  BATCH is the number of sample points to create\n' );
%   fprintf ( 1, '  at one time\n' );
%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  BATCH should be between 1 and SAMPLE_NUM.\n' );
%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  It is FASTER to set BATCH to SAMPLE_NUM;\n' );
%   fprintf ( 1, '  setting BATCH to 1 requires the least memory.\n' );
%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  (Try ''%d'' if you have no preference.)\n', ...
%     min ( 1000, sample_num ) );
%   fprintf ( 1, '  (A zero or negative value terminates execution.)\n' );
%   fprintf ( 1, '\n' );

%   batch = [];
%   batch = input ( '  Enter BATCH:  ' );
    batch = 1000;

    fprintf ( 1, '  Default BATCH = %12d\n', batch );

    if ( batch <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of BATCH = %12d\n', batch );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  MOVIE_NAME is the name of the file in which the movie is stored;\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''movie.avi'' if you have no preference.)\n' );
    fprintf ( 1, '  (A blank value terminates execution).\n' );
    fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
    fprintf ( 1, ' \n' );

    movie_name = [];
    movie_name = input ( '  Enter MOVIE_NAME:  ' );

    fprintf ( 1, '  User input MOVIE_NAME = %s\n', movie_name );

%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  OUTPUT is a file into which the data is stored;\n' );
%   fprintf ( 1, '\n' );
%   fprintf ( 1, '  (Try ''cvt.txt'' if you have no preference.)\n' );
%   fprintf ( 1, '  (A blank value terminates execution).\n' );
%   fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
%   fprintf ( 1, ' \n' );

%   file_out_name = [];
%   file_out_name = input ( '  Enter OUTPUT:  ' );

    file_out_name = 'cvt.txt';

    fprintf ( 1, '  Default OUTPUT = "%s".\n', file_out_name );

    if ( s_len_trim ( file_out_name ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of OUTPUT \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    if ( s_len_trim ( file_out_name ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CVT_MOVIE\n' );
      fprintf ( 1, '  The input value of OUTPUT \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end
%
%  Initialize the data.
%
    if ( init == 3 )
      r = data_read ( sample_string, ndim, n );
    else
      r = [];
    end

    seed_init = seed;
%
%  Initialize the data unless the user has already done that.
%
    if ( init ~= 3 )

      initialize = 1;

      [ r, seed ] = cvt_sample ( ndim, n, n, init, initialize, seed );

    end
%
%  If the initialization and sampling steps use the same random number
%  scheme, then the sampling scheme does not have to be initialized.
%
    if ( init == sample )
      initialize = 0;
    else
      initialize = 1;
    end

    if ( file_exist ( movie_name ) )
      file_delete ( movie_name )
    end

    if ( DEBUG )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Step    SEED          L2-Change     Energy\n' );
      fprintf ( 1, '\n' );
    end

    num_frames_per_second = 10;
    aviobj = avifile ( movie_name, 'fps', num_frames_per_second );
    box = [ 0.0, 0.0; 1.0, 0.0; 1.0, 1.0; 0.0, 1.0 ]';

    it_num = 0;
    it_diff = 0.0;
    seed_base = seed_init;
    seed = seed_init;

    while ( it_num < it_max )
%
%  Once every IT_FIXED steps, update the base value of SEED,
%  either to SEED_INIT before the first call to CVT_ITERATE, or to the
%  output value of SEED from the previous call to CVT_ITERATE.
%
%  Otherwise, reset the value of SEED to SEED_BASE.
%
      if ( mod ( it_num, it_fixed ) == 0 )
        seed_base = seed;
      else
        seed = seed_base;
      end

      it_num = it_num + 1;

      seed_init2 = seed;

      [ r, seed, it_diff, energy ] = cvt_iterate ( ndim, n, batch, sample, ...
        initialize, sample_num, seed, r );

      scatter ( r(1,1:n), r(2,1:n), [], 'r', 'filled' );
      line ( [ 0.0, 1.0, 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0, 1.0, 0.0 ] );
      axis ( [ -0.05, 1.05, -0.05, 1.05 ] );
      axis equal
%
%  Label the axes and the plot.
%
    xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
      'FontSize', 16 );

    ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
      'FontSize', 16, 'Rotation', 0 );

    it_string = sprintf ( 'Step %4d', it_num );

    title ( it_string, 'FontName', 'Helvetica', 'FontWeight', ...
      'bold', 'FontSize', 16 );

      frame = getframe ( gca );
      aviobj = addframe ( aviobj, frame );

      initialize = 0;

      if ( DEBUG )
        fprintf ( 1, '  %4d  %12d  %14e  %14e\n', ...
          it_num, seed_init2, it_diff, energy );
      end

    end

    aviobj = close ( aviobj );
%
%  Write the final data to a file.
%
    cvt_write ( ndim, n, batch, seed_init, seed, init_string, it_max, ...
      it_fixed, it_num, it_diff, energy, sample_string, sample_num, ...
      r, file_out_name );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  The data was written to the file "%s".\n', ...
      file_out_name );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Press RETURN to proceed.\n' );

    pause

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Final value of SEED = %d\n', seed );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = Spread(S1,S2)
    nsim = size(S1,1);
    nt = size(S1,2);
    A = S1(:,nt)./S1(:,1) - S2(:,nt)./S2(:,1);
    val = [max(A,[],2) zeros(nsim,1)];
    y = mean(max(val,[],2));
end

author: Sable

repo: mcbench-benchmarks


text:
function [vMc, mStd_McBoot, fStd_Mc, blo, bhi] = calc_BstMcPercInt(mCatalog,fBinning, alpha)
% [vProbMin, vMcBest, mMag_bstsamp, fStd_Mc, fConfLow, fConfUp] = calc_BstMcPerInt(mCatalog, fBinning)
%---------------------------------------------------------------------------
% Bootstrap EQ catalog and determine Mc confidence due to modelling with
% MLS fitting entire FMD distribution
%
% Incoming variables:
% mCatalog   : EQ catalog
% fBinning   : Magnitude binning interval
%
% Outgoing variables:
% vProbMin     : Vector of maximum likelihood scores
% vMcBest      : Best Mc estimate according to MLS
% mMag_bstsamp : Matrix of bootstrap samples of magnitudes
% fStd_Mc      : Standard deviation (assuming normal distribution
%
% J. Woessner: woessner@seismo.ifg.ethz.ch
% last update: 31.01.03

% Initialize
vMls = [];
vMc = [];
vStd_Mcboot = [];
vMls_boot = [];
vMc_boot = [];
mStd_McBoot = [];

nB = 100; % Bootstrap replicates to calculate confidence intervals

% Get magnitudes
vMags = mCatalog(:,6);

% Create bootstrap samples using bootstrap matlab toolbox
nSample = 50;
mMag_bstsamp = bootrsp(vMags,nSample);
!date
% First step: Estimate standard error of Mc original distribution
% Determine Mc uncertainty
for nSamp=1:nSample
    mCatalog(:,6) = mMag_bstsamp(:,nSamp);
    [mResult, fMls, fMc, fMu, fSigma, mDatPredBest, vPredBest] = calc_McCdfnormal(mCatalog, fBinning);
    vMls = [vMls; fMls];
    vMc =  [vMc; fMc];
end
fStd_Mc = std(vMc);
!date
% Second step: Estimate standard error for each Mc of a bootstrap sample
for i = 1:nB
    i
    mCatalog(:,6) = mMag_bstsamp(:,i);
    mMag_bst = bootrsp(mCatalog(:,6),nSample);
    % Bootstrap using the sample for estimating Mc
    for n= 1:nSample
        mCatalog(:,6) = mMag_bstsamp(:,nSample);
        [mResult, fMls, fMc, fMu, fSigma, mDatPredBest, vPredBest] = calc_McCdfnormal(mCatalog, fBinning);
        vMls_boot = [vMls_boot; fMls];
        vMc_boot =  [vMc_boot; fMc];
    end
    fStd_Mcboot = std(vMc_boot);
    vStd_Mcboot = [vStd_Mcboot; fStd_Mcboot];
    mStd_McBoot = [mStd_McBoot; vStd_Mcboot];
end

% Percentile interval
fExp = nB*alpha/2;
sbval = sort(mStd_Mcboot);
blo = sbval(k);
bhi = sbval(nB-k);

author: CelsoReyes


text:
function h = m_scatter(long,lat,varargin)
% M_SCATTER Scatter/bubble plot
%    M_SCATTER(LON,LAT,Y,S,C) displays colored circles at the locations specified
%    by the vectors LON and LAT (which must be the same size).
%
%    S determines the area of each marker (in points^2). S can be a
%    vector the same length a X and Y or a scalar. If S is a scalar,
%    MATLAB draws all the markers the same size. If S is empty, the
%    default size is used.
%
%    C determines the colors of the markers. When C is a vector the
%    same length as X and Y, the values in C are linearly mapped
%    to the colors in the current colormap. When C is a
%    length(X)-by-3 matrix, it directly specifies the colors of the
%    markers as RGB values. C can also be a color string. See ColorSpec.
%
%    M_SCATTER(LON,LAT) draws the markers in the default size and color.
%    M_SCATTER(LON,LAT,S) draws the markers at the specified sizes (S)
%    with a single color. This type of graph is also known as
%    a bubble plot.
%    M_SCATTER(...,M) uses the marker M instead of 'o'.
%    M_SCATTER(...,'filled') fills the markers.
%
%    H = M_SCATTER(...) returns handles to the scatter objects created.
%
%    Use M_PLOT for single color, single marker size scatter plots.
%

% RP 2016

global MAP_PROJECTION MAP_VAR_LIST

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

if nargin < 2
  help m_scatter
  return
end

[x,y] = m_ll2xy(long,lat);

h=scatter(x,y,varargin{:});

set(h,'tag','m_scatter');

if nargout == 0
  clear h
end

return

author: goGPS-Project


text:
function g = ipopt_callback_g(x,model)

global latest_x_g
global latest_G
global latest_g
x = x(:);

% Compute the nonlinear terms in the constraints and Jacobians for later
[g,geq,dg,dgeq] = fmincon_con_liftlayer(x,model);

% Append with linear constraints
g = [g;geq];
if ~isempty(model.A)
    g = [g;model.A*x - model.b];
end
if ~isempty(model.Aeq)
    g = [g;model.Aeq*x - model.beq];
end

% Append with Jacobians with linear terms
G = [dg';dgeq'];
if ~isempty(model.A)
    G = [G;model.A];
end
if ~isempty(model.Aeq)
    G = [G;model.Aeq];
end

% Save the Jacobian, and information about for which x it was computed
latest_G = sparse(G);
latest_x_g = x;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ handle,handle_test,handle_train] = ml_plot_cv_grid_states_regression(stats,parameterse,options)
%ML_PLOT_CV_GRID_STATES_REGRESSION Plots the results of grid search
% K-fold Cross Validation on regression functions
%
%   input -----------------------------------------------------------------
%
%       o stats     : struct,   multi-layer structure, see output of
%                               ml_get_cv_grid_states_regression.m
%
%       o options   : struct,   plot options.
%
%               options.title = 'title_figure';
%               options.metrics = cell array of metrics to plot.
%
%                   - options.metrics = {'mse','nmse'};
%
%   output ----------------------------------------------------------------
%
%       o handle : figure handle
%
%   Documentation --------------------------------------------------------.
%
%   These are the list of metrics provided by stats and can be added to the options.metrics list:
%
%   '1'  - mean squared error                               (mse)
%   '2'  - normalised mean squared error                    (nmse)
%   '3'  - root mean squared error                          (rmse)
%   '4'  - normalised root mean squared error               (nrmse)
%   '5'  - mean absolute error                              (mae)
%   '6'  - mean  absolute relative error                    (mare)
%   '7'  - coefficient of correlation                       (r)
%   '8'  - coefficient of determination                     (d)
%   '9'  - coefficient of efficiency                        (e)
%  '10'  - maximum absolute error                           (me)
%  '11'  - maximum absolute relative error                  (mre)
%
%
%% Input processing

title_name  = 'add a title in options.title';
metrics     = {};
num_metrics = 0;
para_name   = 'parameters';

if isfield(options,'title'),        title_name      = options.title;        end
if isfield(options,'metrics'),      metrics         = options.metrics;      end
if isfield(options,'para_name'),    para_name       = options.para_name;    end

%% Check input

if isempty(metrics)
    error('options.metrics = {} is empty, you should provide at least one metric to plot; ex: mse');
end

num_metrics = length(metrics);

%% Plot

handle = figure('Color', [1 1 1]);
hold on; box on; grid on;

handle_test  = zeros(1,num_metrics);
handle_train = zeros(1,num_metrics);

% 1 Parameter to do grid search on
if size(parameterse,1) == 1
    x_index      = 1:1:length(parameterse);

    for i=1:num_metrics

        metric_name      = metrics{i};

        [h_train,h_test] = plot_metric(x_index,stats,metric_name,num_metrics);

        handle_train(i)  = h_train;
        handle_test(i)   = h_test;

    end

    if num_metrics > 1
        legend(handle_test,metrics);
    else
        legend([h_train,h_test],'Train','Test');
    end

    title(title_name);

    ax            = gca;
    ax.XTick      = x_index;
    ax.XTickLabel = parameterse;
    xlabel(para_name);

    ylabel(metric_name);

% 2 Parameters to do grid search on
elseif size(parameterse,1) == 2
    for i=1:num_metrics
        colormap hot;
        x = parameterse(1,:);
        y = parameterse(2,:);

        h1 = subplot(1, 2, 1);
        z = stats.train.(metrics{i}).mean;
        contourf(x,y,z)
        title(h1,'Train Mean Accuracy')
        set(h1,'XTick',1:length(parameterse(1,:)));
        set(h1,'XTickLabel', parameterse(1,:));
        xlabel(h1,options.para_names{1});
        set(h1,'YTick',1:length(parameterse(2,:)));
        set(h1,'YTickLabel',parameterse(2,:));
        ylabel(h1,options.para_names{2});
        colorbar
        grid off
        axis square

        h2 = subplot(1, 2, 2);
        z = stats.test.(metrics{i}).mean;
        contourf(x,y,z)
        title(h2,'Test Mean Accuracy')
        set(h2,'XTick',1:length(parameterse(1,:)));
        set(h2,'XTickLabel', parameterse(1,:));
        xlabel(h2,options.para_names{1});
        set(h2,'YTick',1:length(parameterse(2,:)));
        set(h2,'YTickLabel',parameterse(2,:));
        ylabel(options.para_names{2});
        ylabel(h2,options.para_names{2});
        colorbar
        grid off
        axis square
        suptitle([title_name ' : ' metrics{i}])
    end

% 3 Parameters to do grid search on
elseif size(parameterse,1) == 3
    for i=1:num_metrics
        colormap hot;
        x = parameterse(3,:);
        y = parameterse(2,:);
        z = parameterse(1,:);
        [X, Y, Z] = meshgrid(x, y, z);
        V = stats.train.(metrics{i}).mean;

        h1 = subplot(1, 2, 1);
        slice(X, Y, Z, V, [], y, []);hold on;
        contourslice(X, Y, Z, V, [], y, []);
        shading interp
        set(gcf,'renderer', 'zbuffer');
        xlabel(options.para_names{3})
        ylabel(options.para_names{2})
        zlabel(options.para_names{1})
        title(strcat(title_name,'Train NMSE'),'Interpreter','Latex')
        cb1 = colorbar('southoutside');
        alpha(0.95)
        view(-56, 17)
        axis square

        h2 = subplot(1, 2, 2);
        V = stats.test.(metrics{i}).mean;

        slice(X, Y, Z, V, [], y, []);hold on;
        contourslice(X, Y, Z, V, [], y, []);
        shading interp
        set(gcf,'renderer', 'zbuffer');
        xlabel(options.para_names{3})
        ylabel(options.para_names{2})
        zlabel(options.para_names{1})
        title(strcat(title_name,'Test NMSE'),'Interpreter','Latex')
        cb2 = colorbar('southoutside');
        alpha(0.95)
        view(-56, 17)
        axis square

    end

end

end

function [h_train,h_test] = plot_metric(x_index,stats,metric_name,num_metrics)

if strcmp(metric_name,'mse')
    train_means = stats.train.mse.mean;
    train_stds  = stats.train.mse.std;
    test_means  = stats.test.mse.mean;
    test_stds   = stats.test.mse.std;

elseif strcmp(metric_name,'nmse')

    train_means = stats.train.nmse.mean;
    train_stds  = stats.train.nmse.std;
    test_means  = stats.test.nmse.mean;
    test_stds   = stats.test.nmse.std;

elseif strcmp(metric_name,'rmse')

    train_means = stats.train.rmse.mean;
    train_stds  = stats.train.rmse.std;
    test_means  = stats.test.rmse.mean;
    test_stds   = stats.test.rmse.std;

elseif strcmp(metric_name,'nrmse')

    train_means = stats.train.nrmse.mean;
    train_stds  = stats.train.nrmse.std;
    test_means  = stats.test.nrmse.mean;
    test_stds   = stats.test.nrmse.std;

elseif strcmp(metric_name,'mae')

    train_means = stats.train.mae.mean;
    train_stds  = stats.train.mae.std;
    test_means  = stats.test.mae.mean;
    test_stds   = stats.test.mae.std;

elseif strcmp(metric_name,'mare')

    train_means = stats.train.mare.mean;
    train_stds  = stats.train.mare.std;
    test_means  = stats.test.mare.mean;
    test_stds   = stats.test.mare.std;

elseif strcmp(metric_name,'r')

    train_means = stats.train.r.mean;
    train_stds  = stats.train.r.std;
    test_means  = stats.test.r.mean;
    test_stds   = stats.test.r.std;

elseif strcmp(metric_name,'d')

    train_means = stats.train.d.mean;
    train_stds  = stats.train.d.std;
    test_means  = stats.test.d.mean;
    test_stds   = stats.test.d.std;

elseif strcmp(metric_name,'e')

    train_means = stats.train.e.mean;
    train_stds  = stats.train.e.std;
    test_means  = stats.test.e.mean;
    test_stds   = stats.test.e.std;

elseif strcmp(metric_name,'me')

    train_means = stats.train.me.mean;
    train_stds  = stats.train.me.std;
    test_means  = stats.test.me.mean;
    test_stds   = stats.test.me.std;

elseif strcmp(metric_name,'mre')

    train_means = stats.train.mre.mean;
    train_stds  = stats.train.mre.std;
    test_means  = stats.test.mre.mean;
    test_stds   = stats.test.mre.std;

else
    error(['no such metric: ' metric_name]);
end

    [h_train,h_test] = plot_statisics(x_index,train_means,train_stds,test_means,test_stds,num_metrics);

end

function [h_train,h_test] = plot_statisics(x_index,train_means,train_stds,test_means,test_stds,num_metrics)
%
%   input -----------------------------------------------------------------
%
%       o x_index       : x values for x-y plot
%
%       o train_means   : (1 x N), set of means of k-CV
%       o train_stds    : (1 x N), set of standard deviations associated with
%                                  means.
%
%       o test_means    : same as train
%       o test_stds     : same as train
%
%       o color_train   : (1 x 3), RGB color of train plot
%       o color_test    : (1 x 3), RGB color of test plot
%
%       o num_metrics   : (1 x 1), number of metrics being plotted
%
%       o metric_name   : string, metrics name.
%

if num_metrics > 1

    h_train           = errorbar(x_index,train_means,train_stds,'--s');
    h_test            = errorbar(x_index,test_means,test_stds,'-s','Color',h_train.Color);
else
    h_train           = errorbar(x_index,train_means,train_stds,'-gs','Color', [0 1 0]);
    h_test            = errorbar(x_index,test_means,test_stds,'-rs','Color', [1 0 0]);
end

end

author: epfl-lasa


text:
function p_ = EntropyProg(p,A,b,Aeq,beq)
% This function computes the entropy-pooling change of measure, see
% "A. Meucci - Fully Flexible Views: Theory and Practice -
% The Risk Magazine, October 2008, p 97-102"
% available at www.symmys.com > Research > Working Papers

% Code by A. Meucci, November 2008
% Last version available at www.symmys.com > Teaching > MATLAB

K_=size(A,1);
K=size(Aeq,1);
A_=A';
b_=b';
Aeq_=Aeq';
beq_=beq';
x0=zeros(K_+K,1);
InqMat=-eye(K_+K); InqMat(K_+1:end,:)=[];
InqVec=zeros(K_,1);

options = optimset('GradObj','on','Hessian','on');
if ~K_
    v=fminunc(@nestedfunU,x0,options);
    p_=exp(log(p)-1-Aeq_*v);
else
    lv=fmincon(@nestedfunC,x0,InqMat,InqVec,[],[],[],[],[],options);
    l=lv(1:K_);
    v=lv(K_+1:end);
    p_=exp(log(p)-1-A_*l-Aeq_*v);
end

    function [mL g H] = nestedfunU(v)

        x=exp( log(p)-1-Aeq_*v );
        x=max(x,10^(-32));
        L=x'*(log(x)-log(p)+Aeq_*v)-beq_*v;
        mL=-L;

        g = [beq-Aeq*x];
        H = [Aeq*((x*ones(1,K)).*Aeq_)];  % Hessian computed by Chen Qing, Lin Daimin, Meng Yanyan, Wang Weijun
    end

    function [mL g H] = nestedfunC(lv)

        l=lv(1:K_);
        v=lv(K_+1:end);
        x=exp( log(p)-1-A_*l-Aeq_*v );
        x=max(x,10^(-32));
        L=x'*(log(x)-log(p))+l'*(A*x-b)+v'*(Aeq*x-beq);
        mL=-L;

        g = [b-A*x; beq-Aeq*x];
        H = [A*((x*ones(1,K_)).*A_)  A*((x*ones(1,K)).*Aeq_) % Hessian computed by Chen Qing, Lin Daimin, Meng Yanyan, Wang Weijun
            Aeq*((x*ones(1,K_)).*A_)   Aeq*((x*ones(1,K)).*Aeq_)];
    end
end


text:
%%*******************************************************
%% maxcut: MAXCUT problem.
%%
%% (P)    min  Tr C*X
%%        s.t.  diag(X) = e;
%%
%% C = -(diag(B*e)-B)/4.
%%
%% (dual problem)   max  e'*y
%%                  s.t. diag(y) + S = C.
%%-------------------------------------------------------
%%
%% [X,y,objval] = maxcut(B);
%%
%% B: weighted adjacency matrix of a graph.
%%
%% DSDP5.0
%% Copyright (c) 2004 by
%% S. Benson, Y. Ye
%% Last modified: Jan 2004
%%*******************************************************

 function [X,y,objval] = maxcut(B);

   if ~isreal(B); error('only real B allowed'); end;

   n = length(B); e = ones(n,1);
   CC = -(spdiags(B*e,0,n,n)-B)/4.0;
   b = e;

   AC = cell(1,3);
   AC{1,1}='SDP';
   AC{1,2}=n;
   nn=n*(n+1)/2;

   AAC=sparse(nn,n);
   for k = 1:n; AAC(:,k) = dsparse(k,k,1,n,n); end;
   AAC=[AAC sparse(dvec(CC))];
   AC{1,3}=AAC;
   y0 = -1.1*abs(CC)*e;
   OPTIONS=doptions;
   OPTIONS.gaptol=0.0001;
   OPTIONS.r0=0;
   OPTIONS.print=1;
   OPTIONS.rho=5;
   OPTIONS.zbar=norm(B,1);
   [STAT,y,Xv] = dsdp(b,AC,OPTIONS,y0);

   [objval,dobj,err] = derror(STAT,y,Xv,b,AC);
   X=dmat(Xv{1});
return;
%%=======================================================

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function r = issubfield(s, f)

% ISSUBFIELD tests for the presence of a field in a structure just like the standard
% Matlab ISFIELD function, except that you can also specify nested fields
% using a '.' in the fieldname. The nesting can be arbitrary deep.
%
% Use as
%   f = issubfield(s, 'fieldname')
% or as
%   f = issubfield(s, 'fieldname.subfieldname')
%
% This function returns true if the field is present and false if the field
% is not present.
%
% See also ISFIELD, GETSUBFIELD, SETSUBFIELD

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

%try
%  getsubfield(s, f);    % if this works, then the subfield must be present
%  r = true;
%catch
%  r = false;                % apparently the subfield is not present
%end

if isempty(f) || isempty(s)
  r = false;
else
  t = textscan(f,'%s','delimiter','.');
  t = t{1};
  r = true;
  for k = 1:numel(t)
    if isfield(s, t{k})
      s = s.(t{k});
    else
      r = false;
      return;
    end
  end
end


text:
function b = ContrastChange(img,level)

img = im2double(img)*255;
img= img*level;
b = uint8(img);
end

author: xialeiliu

repo: RankIQA

save_path: github-repos/MATLAB/xialeiliu-RankIQA


text:
function esvm_show_exemplar_frames(allmodels, N_PER_PAGE, dataset_params)
% Draw the initialized exemplar frames as a 1x3 row of 3 images
% Shows these 3 fields: input+gtbb+template, template mask+gtbb, HOG descriptor
% The visualization really shows what the template region is, and
% its relation to the ground-truth selection region
%
% Copyright (C) 2011-12 by Tomasz Malisiewicz
% All rights reserved.
%
% This file is part of the Exemplar-SVM library and is made
% available under the terms of the MIT license (see COPYING file).
% Project homepage: https://github.com/quantombone/exemplarsvm

if ~exist('dataset_params','var')
  dataset_params = [];
end

pinds = do_partition(1:length(allmodels),N_PER_PAGE);
for i = 1:length(pinds)
  models = allmodels(pinds{i});

  figure(i)
  clf
  N = length(models);
  for i = 1:N
    m = models{i};
    o = (i-1)*3;
    subplot(N,3,o+1)
    I = convert_to_I(m.I);
    imagesc(I)
    plot_bbox(m.model.bb(1,:),'',...
              [1 0 0], [0 1 0], 0 ,[1 3],...
              m.model.hg_size)
    plot_bbox(m.gt_box,'',[0 0 1])
    axis image
    axis off
    title(sprintf('Ex %s.%d %s',m.curid,m.objectid,m.cls))

    subplot(N,3,o+2)
    onimage = m.model.mask*1;
    onimage(onimage==0) = 2;
    colors = [1 0 0; 0 0 1];
    cim = colors(onimage(:),:);
    cim = reshape(cim,[size(m.model.mask,1) size(m.model.mask,2) 3]);
    imagesc(cim);
    fullimbox = [0 0 size(cim,2) size(cim,1)]+.5;
    xform = find_xform(m.model.bb(1,1:4), fullimbox);
    gtprime = apply_xform(m.gt_box,xform);
    plot_bbox(fullimbox,'',...
              [1 0 0], [0 1 0], 0 ,[1 3],...
              m.model.hg_size)
    plot_bbox(gtprime,'',[0 0 1])
    axis image
    axis off
    grid on
    [u,v] = find(m.model.mask);

    curselection = [min(v) min(u) max(v) max(u)];
    curos = getosmatrix_bb(curselection, gtprime);

    title(sprintf('%s: Template:[%d x %d] \ncuros=%.2f Mask: [%d x %d]',...
                  m.model.init_params.init_type,...
                  m.model.hg_size(1),m.model.hg_size(2),...
                  curos,range(u)+1,range(v)+1));

    subplot(N,3,o+3)
    hogim = HOGpicture(repmat(m.model.mask,[1 1 size(m.model.w,3)]).* ...
                       m.model.w);
    imagesc(hogim)
    axis image
    axis off
    grid on
    title('HOG features')
    drawnow
  end
end

author: quantombone

repo: exemplarsvm