text:
function [xo,Ot,nS]=buscarnd(S,x0,ip,nOt,samples,Lb,Ub,problem,tol,mxit,R,red,mem)
%   Unconstrained global optimization using adaptive random search.
%
%   [xo,Ot,nS]=buscarnd(S,x0,ip,nOt,samples,Lb,Ub,problem,tol,mxit,R,red,mem)
%
%   S: objective function
%   x0: initial point
%   ip: (0): no plot (default), (>0) plot figure ip with pause, (<0) plot figure ip
%   nOt: maximum number of optimal points (default = 1)
%   samples: number of samples per stage (default = 3*size(x0(:),1))
%   Lb, Ub: lower and upper bound vectors to plot (default = x0*(1+/-2))
%   problem: (-1): minimum (default), (1): maximum
%   tol: tolerance (default = 1e-4)
%   mxit: maximum number of stages (default = 50*(1+4*~(ip>0)))
%   R: axis vector of the hyperellipse centered in x0 (default = max(0.1*abs(x0+~x0),1))
%   red: factor to reduce the size of the axis vector (0,1) (default = 0.2)
%   mem: number of stored points during the exploration phase [0, min(18,samples-1)]
%   xo: matrix of optimal points (column vectors)
%   Ot: vector of optimal values of S
%   nS: number of objective function evaluations

%   Copyright (c) 2001 by LASIM-DEQUI-UFRGS
%   $Revision: 1.0 $  $Date: 2001/07/10 19:40:05 $
%   Argimiro R. Secchi (arge@enq.ufrgs.br)
%
%   Based on the algorithm of the same author written in C
%   for constrained global optimization in 1989/09/29 and published as:
%   A.R. Secchi and C.A. Perlingeiro, "Busca Aleatoria Adaptativa",
%   in Proc. of XII Congresso Nacional de Matematica Aplicada e
%   Computacional, Sao Jose do Rio Preto, SP, pp. 49-52 (1989).
%
% Modified by Giovani Tonel(giovani.tonel@ufrgs.br) on September 2006

rand('state',sum(100*clock)); 	%rand('state',sum(100*clock)) resets it to
													%		         a different state each time.

if nargin < 2,
   error('buscarnd requires two input arguments');
 end
 if nargin < 3 | isempty(ip),
   ip=0;
 end
 if nargin < 4 | isempty(nOt),
   nOt=1;
 end
 if nargin < 5 | isempty(samples),
   samples=3*size(x0(:),1);
 end
 if nargin < 6 | isempty(Lb),
   Lb=-x0-~x0;
 end
 if nargin < 7 | isempty(Ub),
   Ub=2*x0+~x0;
 end
 if nargin < 8 | isempty(problem),
   problem=-1;
 end
 if nargin < 9 | isempty(tol),
   tol=1e-4;
 end
 if nargin < 10 | isempty(mxit),
   mxit=50*(1+4*~(ip>0));
 end
 if nargin < 11 | isempty(R),
   R=max(0.1*abs(x0+~x0),1);
 end
 if nargin < 12 | isempty(red) | red <= 0 | red >= 1,
   red=0.2;
 end
 if nargin < 13 | isempty(mem),
   mem=min(18,samples-1);
 end
 if mem >= samples,
   mem=samples-1;
 end
                    % initialization
 x0=x0(:);
 R=abs(R(:));
 n=size(x0,1);
 distrib=1;          % distr =  1  --> exploration phase
                     % distr >  1  --> progression phase
 asym1=-1*ones(n,1); % asym1 =  1  --> asymmetric distribution to reduce xo
                     % asym1 = -1  --> asymmetric distribution to increase xo
 asym2=2*ones(n,1);  % asym2 =  2  --> symmetric distribution
                     % asym2 = 1.5 --> asymmetric distribution
 metric=0.5*norm(R);
 a0=0.5*metric;
 a2=0.1*metric;
 a3=10*tol;
 n3=1;
 holdm=floor(1/red+0.5)*n;
 nhold=0;            % counter to keep distribution constant during holdm times

 top=0;
 if mem > 0,
   n4=10*mem;
   idx=[1:n4+1];
   last=n4;
   first=0;
   idx(last+1)=first;
   mem_x=zeros(n,n4);
   mem_S=-ones(1,n4)*inf;
 end

%     Criterion of axis vector reduction
%
%     R_red = R*0.5^(distrib-1)
%
%     or norm(R_red) = norm(R)*0.5^(distrib-1)
%
%     then if the criterion is norm(R_red) < tol = 10^(-asc)
%
%     the distrib is limited by
%
%     distrib = 1 + (asc + log10(norm(R))) / log10(2)

 lim_distrib=floor(0.5*(1+log10(metric/tol)/log10(2.0)));

 if ip & n == 2,
   figure(abs(ip));
   [X1,X2]=meshgrid(Lb(1):(Ub(1)-Lb(1))/20:Ub(1),Lb(2):(Ub(2)-Lb(2))/20:Ub(2));
   [n1,n2]=size(X1);
   f=zeros(n1,n2);
   for i=1:n1,
    for j=1:n2,
      f(i,j)=feval(S,[X1(i,j);X2(i,j)]);
    end
   end
   mxf=max(max(f));
   mnf=min(min(f));
   df=mnf+(mxf-mnf)*(2.^(([0:10]/10).^2)-1);
   [v,h]=contour(X1,X2,f,df); hold on;
%   clabel(v,h);
   h1=plot(x0(1),x0(2),'ro');
   legend(h1,'start point');

   if ip > 0,
     disp('Pause: hit any key to continue...'); pause;
   end
 end

 xOt=[];
 Ot=[];
 xo=x0;
 yo=feval(S,x0)*problem;
 nS=1;
 it=0;
 opt=0;
 x=zeros(n,samples);
 y=zeros(1,samples);

 while it < mxit & opt < nOt,
  l3=0;
  it=it+1;
  for j=1:samples,   % sampling
    a=min(max(rand(n,1),0.05),0.95);
    x(:,j)=xo+(R.*asym1.*(asym2.*a-1).^distrib)/distrib;
    y(j)=feval(S,x(:,j))*problem;
    nS=nS+1;
  end

  l4=mem & distrib == 1;  % sorting the samples
  if samples > 1,
    if ~l4,
      [y(1),i]=max(y);
      x(:,1)=x(:,i);
    else
      [y,i]=sort(-y); y=-y;
      x=x(:,i);
    end
  end

  if l4,                  % memorize samples
    n2=top;
    l1=0;
    for j=1:samples,
      n1=0;
      if opt,
        for i=1:opt
          if norm(xOt(:,i)-x(:,j)) < a2,
            n1=1;
            if j == 1,
              l1=1;
            end
            break;
          end
        end
      end

      if ~n1,
	for i=1:j-1
	   if norm(x(:,j)-x(:,i)) < metric,
	     n1=1;
	     break;
	   end
	end

	if ~n1,
	  k=idx(1);
	  for i=1:top
             if norm(mem_x(:,k)-x(:,j)) < metric,
               n1=1;
               break;
             else
               k=idx(k+1);
             end
          end

	  if ~n1,
            first=idx(first+1);
            if ~first,              % expand memory
              mem_x=[mem_x,zeros(n,n4)];
              mem_S=[mem_S,-ones(1,n4)*inf];
              idx=[idx,[n2+1:n2+n4+1]];
              first=n2;
              idx(last+1)=first;
              last=n2+n4;
              idx(last+1)=0;
            end
            k=first;
            n2=n2+1;
          end

          if ~n1 | y(j) > mem_S(k),
            mem_x(:,k)=x(:,j);
            mem_S(:,k)=y(j);
          end

	  if l1,
            l1=0;
            x(:,1)=x(:,j);
            y(1)=y(j);
          end
        end
      end

      if n2-top == mem,
        break;
      end
    end
    top=n2;
  end
                  % analysis of best sampled point
  a1=norm(x(:,1)-xo)/(0.1+norm(xo))+abs(y(1)-yo)/(0.1+abs(yo));
  l1=y(1) > yo;
  if l1,
    if norm(x(:,1)-xo) > a0,
      nhold=0;
      n3=1;
    end

    z=xo; xo=x(:,1); x(:,1)=z;
    a4=yo; yo=y(1); y(1)=a4;

    if ip & n == 2,
      plot([x(1) xo(1)],[x(2) xo(2)],'r');
      if ip > 0,
        disp('Pause: hit any key to continue...'); pause;
      end
    end
  end

  if a1 < tol & distrib >= lim_distrib,
    l2=1;
    if opt,
      for i=1:opt,
        if norm(xOt(:,i)-xo) < a3,
          l2=0;
          break;
        end
      end
    end

    if l2,
      opt=opt+1;
      Ot=[Ot,yo];
      xOt=[xOt,xo];
      if ip & n == 2,
        h2=plot(xo(1,:),xo(2,:),'r*');
        if opt == 1,
          legend([h1,h2],'start point','optimum');
        end
      end
    end

        % rescue another point from memory
    if ~top,
      if opt,
        xo=xOt;
        yo=Ot;
      end
      break;
    end

    if l2 & ~l1,
      n2=0;
      for j=1:top,
        k=idx(n2+1);
        if norm(mem_x(:,k)-xo) < a2,
          top=top-1;
          idx(last+1)=k;
          last=k;
          idx(n2+1)=idx(k+1);
          idx(k+1)=0;
          if k == first,
            first=n2;
          end
        else
          n2=k;
        end
      end
    end

    if ~top,
      if opt,
        xo=xOt;
        yo=Ot;
      end
      break;
    end

    k=idx(1);
    xo=mem_x(:,k);
    yo=mem_S(k);
    idx(last+1)=k;
    last=k;
    idx(1)=idx(k+1);
    idx(k+1)=1;
    top=top-1;

    if k == first,
      first=1;
    end

    l3=1;

    if ip & n == 2 & opt < nOt,
      plot(xo(1),xo(2),'ro');
    end
  elseif opt & (~l1 | (l1 & ~l4))
        n1=0;
	for i=1:opt,
          if norm(xOt(:,i)-xo) < a2,
            n1=1;
            break;
          end
        end

	if n1,
          if ~top,
            xo=xOt;
            yo=Ot;
            break;
          end

          k=idx(1);
          xo=mem_x(:,k);
          yo=mem_S(k);
          idx(last+1)=k;
          last=k;
          idx(1)=idx(k+1);
          idx(k+1)=1;
          top=top-1;

          if k == first,
            first=1;
          end

	  l3=1;

          if ip & n == 2 & opt < nOt,
            plot(xo(1),xo(2),'ro');
          end
        end
  end
                 % adjust search direction and criterion
  if l3,
    n3=1;
    distrib=1;
    nhold=0;
  elseif (a1 < a2 | ~l1) & (nhold+n3 > holdm),
        n3 = n3+~mod(distrib,3);
        distrib=distrib+2;
	nhold=0;
  else
    nhold=nhold+n3;
  end

  for i=1:n,
    if l3 | (~l3 & x(i,1) == xo(i)),
      asym2(i)=2;
    else
      asym2(i)=1.5;
      if x(i,1) < xo(i),
        asym1(i)=1;
      else
        asym1(i)=-1;
      end
    end
  end
 end

 if it == mxit,
   %disp('Warning Buscarnd: reached maximum number of stages!');
   if opt,
     yo=Ot;
     xo=xOt;
   end
 end

 if opt == nOt,
   yo=Ot;
   xo=xOt;
 end

 if opt > 1,
   [yo,i]=sort(-yo); yo=-yo;
   xo=xo(:,i);
 end

 Ot=yo*problem;


text:
function Show_EWT(ewt,f,rec)

% ====================================================================
% function Show_EWT(f,ewt)
%
% This function plots the successive filtered components (low scale
% first and then wavelets scales). The original and
% reconstructed signals are plotted on a different graph.
% If f and rec are provided, it also plot the original and reconstructed
% signals on a separate figure
%
% Inputs:
%   -ewt: EWT components
%   -f: input signal  (OPTIONNAL)
%   -rec: reconstructed signal (OPTIONNAL)
%
% Author: Jerome Gilles
% Institution: UCLA - Department of Mathematics
% Year: 2012
% Version: 1.0
% =====================================================================

%These lines plot the EWT components
figure;
x=0:1/length(ewt{1}):(length(ewt{1})-1)/length(ewt{1});
l=1;
if length(ewt)>6
    lm=6;
else
    lm=length(ewt);
end

for k=1:length(ewt)
   hold on; subplot(lm,1,l); plot(x,ewt{k});
   if mod(k,6) == 0
        figure;
        l=1;
   else
    l=l+1;
   end
end

%These lines plot f and its reconstruction
if nargin>1
    figure;
    subplot(2,1,1);plot(x,f);
    subplot(2,1,2);plot(x,rec);
end

author: Sable

repo: mcbench-benchmarks


text:
function RMC = fmri_mcorrestriction(RM,TR,nHEst,Delta,Tau)
%
% RMC = fmri_mcorrestriction(RM,TR,nHEst,Delta,Tau)
%
% Creates a restriction matrix for correlating with a
% hemodynmic response function (HRF).  The HRF is
% computed using a Gamma function with parameters
% Delta and Tau over the estimation time window.
%
% RM is a restriction matrix created by CreateRM.
%
% See also: CreateRM(), HemoDyn()
%
% $Id: fmri_mcorrrestriction.m,v 1.1 2004/03/11 01:30:16 sayres Exp $

t = TR*[0:nHEst-1]';
hHDIR = fmri_hemodyn(t,Delta,Tau);
nCond = size(RM,2)/nHEst; % excluding fix %
q = repmat(hHDIR',size(RM,1),nCond);
RMC = RM .* q;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function uninstall
% AnalogDevices.uninstall removes AnalogDevices HDL BSPs

% Copyright 2015 MathWorks, Inc. All Rights Reserved.

    AnalogDevices.install(1);
end

author: analogdevicesinc


text:
function drawFace(p)
%
% head
% eyes
% nose
% mouth
% eyebrows
% hair
%
% 1 - head width    0 1
% 2 - head height   0 1
% 3-5 head color    0 1
%
% 6     y eye level
% 7     eye spacing
% 8     eye line length
% 9     left eye line width
% 10    right eye line width
%
% 11    nose to eye dist
% 12    nose length
% 13    nose line width
%
% 14    mouth width
% 15    smile scale factor
% 16    smile displacement
% 17    smile line width
%
% 18 x
% 19 y

cla

% ------------------------------------------
% * draw head
% ------------------------------------------
[X,Y,Z]=ellipsoid(p(18),p(19),0,p(1),p(2),1);
surf(X,Y,Z,'EdgeColor','None','FaceColor',[p(3) p(4) p(5)])
%axis off; set(gcf,'Color','w');axis([-1 1 -1 1])
hold on

% ------------------------------------------
% * draw eyes
% ------------------------------------------
y = p(6) * p(2) + p(19);
x = p(18);
s = (p(7) * p(1)) ./ 2;
plot3(x+[-p(8) * p(1) - s -s] ,[y y],[2 2],'k','LineWidth',p(9)*8)
plot3(x+[s (p(8)*p(1)+s)],[y y],[2 2],'k','LineWidth',p(10)*8)

% ------------------------------------------
% * draw nose
% ------------------------------------------
nosestart = y - .5*(p(11) * p(2));
noseend = max(y,p(12) * (p(2) - nosestart));
plot3(x+[0 0],[nosestart noseend],[2 2],'k','LineWidth',p(13)*8)

% ------------------------------------------
% * draw mouth
% ------------------------------------------
x2 = x; % displacement
wid = (p(14) * p(1)) ./ 2;
x = -wid:.05.*p(1):wid;
ys = p(15) * 2 * (x) .^ 2;
ys = y+ ys - (p(16)*p(2));
plot3(x+x2,ys,2*ones(size(x)),'k','LineWidth',p(17)*5)

drawnow

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function y = interpftn(x, sz, win)
%INTERPFTN Resample data using Fourier interpolation.
%
% DESCRIPTION:
%       interpftn resamples an N-D matrix to the size given in sz using
%       Fourier interpolation.
%
% USAGE:
%       y = interpftn(x, sz)
%       y = interpftn(x, sz, win)
%
% INPUTS:
%       x       - matrix to interpolate
%       sz      - new size
%
% OPTIONAL INPUTS:
%       win     - name of windowing function to use
%
% OUTPUTS:
%       y       - resampled matrix
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 20th October 2010
%       last update - 20th December 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox
%
% See also interpft

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% check inputs for window
if nargin == 2
    win = 0;
end

% extract the size of the input matrix
x_sz = size(x);

% check enough coefficients have been given
if sum(x_sz ~= 1) ~= numel(sz)
    error('the number of scaling coefficients must equal the number of dimensions in x');
end

% compute the interpolation over each matrix dimension using the inbuilt
% interpft (dimensions with no interpolation required are skipped)
y = x;
for p = 1:numel(sz)
    if sz(p) ~= x_sz(p)
        y = interpft(y, sz(p), p, win);
    end
end

function y = interpft(x, ny, dim, win)
%INTERPFT 1-D interpolation using FFT method.
%   Y = INTERPFT(X,N) returns a vector Y of length N obtained
%   by interpolation in the Fourier transform of X.
%
%   If X is a matrix, interpolation is done on each column.
%   If X is an array, interpolation is performed along the first
%   non-singleton dimension.
%
%   INTERPFT(X,N,DIM) performs the interpolation along the
%   dimension DIM.
%
%   Assume x(t) is a periodic function of t with period p, sampled
%   at equally spaced points, X(i) = x(T(i)) where T(i) = (i-1)*p/M,
%   i = 1:M, M = length(X).  Then y(t) is another periodic function
%   with the same period and Y(j) = y(T(j)) where T(j) = (j-1)*p/N,
%   j = 1:N, N = length(Y).  If N is an integer multiple of M,
%   then Y(1:N/M:N) = X.
%
%   Example:
%      % Set up a triangle-like signal signal to be interpolated
%      y  = [0:.5:2 1.5:-.5:-2 -1.5:.5:0]; % equally spaced
%      factor = 5; % Interpolate by a factor of 5
%      m  = length(y)*factor;
%      x  = 1:factor:m;
%      xi = 1:m;
%      yi = interpft(y,m);
%      plot(x,y,'o',xi,yi,'*')
%      legend('Original data','Interpolated data')
%
%   Class support for data input x:
%      float: double, single
%
%   See also INTERP1.

%   Robert Piche, Tampere University of Technology, 10/93.
%   Copyright 1984-2006 The MathWorks, Inc.
%   $Revision: 5.15.4.4 $  $Date: 2006/12/15 19:27:56 $

% y = zeros(ny, 1);
% y(1:length(x)) = fftshift(fft(x));
% y = real(ifft(ifftshift(y)));

% Bug fix, additional documentation, and win option
% author        - Bradley Treeby
% date          - 20th October 2010
% last update   - 1st April 2011

% if nargin==2,
%
%     % this forces x to be a column vector, nshifts is so it can be changed
%     % back if necessary
%     [x, nshifts] = shiftdim(x);
%
%     % check if x is a scalar, and return a row vector if so
%     if isscalar(x)
%       nshifts = 1;
%     end
%
% elseif nargin==3,
%
%     % if the dimension to take the fft over is given, rearrange the matrix
%     % so this dimension is the first dimension
%     perm = [dim:max(length(size(x)),dim) 1:dim-1];
%     x = permute(x,perm);
% end

% rearrange the matrix so the fft is taken over the first dimension
perm = [dim:max(length(size(x)),dim) 1:dim-1];
x = permute(x,perm);

siz = size(x);
[m, n] = size(x);
if ~isscalar(ny)
  error('MATLAB:interpft:NonScalarN', 'N must be a scalar.');
end

% if necessary, increase ny by an integer multiple to make ny > m
% here m is the length of the dimension over which the fft is being taken,
% and ny is the number of points to use including zero padding
if ny > m
   incr = 1;
else
    % if ny is zero, return nothing
    if ny == 0
       y = [];
       return
    end
    % otherwise make ny bigger to allow downsampling from integer points
    incr = floor(m/ny) + 1;
    ny = incr*ny;
end

% take the fft of the input function
a = fft(x,[],1);

% calculate the nyquist frequency
nyqst = ceil((m+1)/2);

% window the Fourier coefficient if win is given
if win
    % create window using getWin
    wind = getWin(siz(1), win);

    % shift the window
    wind = circshift(wind, nyqst);

    % create repmat variable
    siz_rep = siz;
    siz_rep(1) = 1;

    % repeat
    wind = repmat(wind, siz_rep);

    % apply the window
    a = a.*wind;

end

% zero pad with the zeros in the middle
b = [a(1:nyqst,:) ; zeros(ny-m, prod(siz(2:end))) ; a(nyqst+1:m,:)];

if rem(m, 2) == 0
    % if the sequence has an even number of points, make the sequence
    % symmetric, e.g., turn
    % EP P P N   0 0 0 0 0 0   P P
    % into
    % EP P P N/2 0 0 0 0 0 N/2 P P
    % where EP = end point P = points, N = nyquist point
    b(nyqst,:) = b(nyqst,:)/2;
    b(nyqst+ny-m,:) = b(nyqst,:);
end

% take the inverse FFT
y = ifft(b,[],1);

% if the input was real, throw away any residual complex bits
if isreal(x)
    y = real(y);
end

% make sure the amplitudes are correct by scaling by the new length
y = y * ny / m;

% this gets a downsampled version from an upsampled version
y = y(1:incr:ny,:);  % Skip over extra points when oldny <= m.

% reshape
[y_length, num_signals] = size(y);
y = reshape(y, [y_length, siz(2:end)]);

% sort out the resizing
y = ipermute(y,perm);

% if nargin==2,
%   y = reshape(y,[ones(1,nshifts) size(y,1) siz(2:end)]);
% elseif nargin==3,
%   y = ipermute(y,perm);
% end

repo: Aero-matlab


text:
%% Copyright (C) 2016 Lagu
%% Copyright (C) 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym frac (@var{x})
%% Return the fractional part of a symbolic expression.
%%
%% Examples:
%% @example
%% @group
%% y = frac(sym(3)/2)
%%   @result{} y = (sym) 1/2
%%
%% syms x
%% rewrite(frac(x), 'floor')
%%   @result{} ans = (sym) x - ⌊x⌋
%% @end group
%% @end example
%%
%% @seealso{@@sym/ceil, @@sym/floor, @@sym/fix, @@sym/round}
%% @end defmethod

function y = frac(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('frac', x);
end

%!test
%! f1 = frac(sym(11)/10);
%! f2 = sym(1)/10;
%! assert (isequal (f1, f2))

%!test
%! d = sym(-11)/10;
%! c = sym(9)/10;
%! assert (isequal (frac (d), c))

%!test
%! d = sym(-19)/10;
%! c = sym(1)/10;
%! assert (isequal (frac (d), c))

repo: octsympy


text:
function [shape] = ft_transform_headshape(transform, shape)

% FT_TRANSFORM_HEADSHAPE applies a homogenous coordinate transformation to a
% structure with headshape and fiducial information.
%
% Use as
%   shape = ft_transform_headshape(transform, shape)
%
% See also FT_READ_HEADSHAPE

% Copyright (C) 2008, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: ft_transform_headshape.m 946 2010-04-21 17:51:16Z roboos $

if any(transform(4,:) ~= [0 0 0 1])
  error('invalid transformation matrix');
end

if isfield(shape, 'pnt') && ~isempty(shape.pnt)
  % this also works if the structure describes electrode or gradiometer positions instead of a headshape
  shape.pnt = apply(transform, shape.pnt);
end

if isfield(shape, 'ori')
  % gradiometer coil orientations should only be rotated and not translated
  rotation = eye(4);
  rotation(1:3,1:3) = transform(1:3,1:3);
  if abs(det(rotation)-1)>10*eps
    error('only a rigid body transformation without rescaling is allowed for MEG sensors');
  end
  % apply the rotation to the coil orientations
  shape.ori = apply(rotation, sens.ori);
end

if isfield(shape, 'fid') && isfield(shape.fid, 'pnt')
  % apply the same transformation on the fiducials
  shape.fid.pnt = apply(transform, shape.fid.pnt);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that applies the homogenous transformation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [new] = apply(transform, old)
old(:,4) = 1;
new = old * transform';
new = new(:,1:3);

author: PatternRecognition

repo: OpenBMI


text:
% Calibrate camera counts for the number of detected photons.
pathfile = './';
dark = readtimeseries([pathfile 'R_20110413/Dark_2/img_000000*__000.tif']);
bright = readtimeseries([pathfile 'R_20110413/Bright_2/img_000000*__000.tif']);
sv= size(bright);
out_cal = cal_readnoise(bright, dark(0:sv(1)-1,0:sv(2)-1,:));

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function x = beta_inv(p, a, b)
% PURPOSE: inverse of the cdf (quantile) of the beta(a,b) distribution
%--------------------------------------------------------------
% USAGE: x = beta_inv(p,a,b)
% where:   p = vector of probabilities
%          a = beta distribution parameter, a = scalar
%          b = beta distribution parameter  b = scalar
% NOTE: mean [beta(a,b)] = a/(a+b), variance = ab/((a+b)*(a+b)*(a+b+1))
%--------------------------------------------------------------
% RETURNS: x at each element of p for the beta(a,b) distribution
%--------------------------------------------------------------
% SEE ALSO: beta_d, beta_pdf, beta_inv, beta_rnd
%--------------------------------------------------------------

% Anders Holtsberg, 18-11-93
% Copyright (c) Anders Holtsberg
% documentation modified by LeSage to
% match the format of the econometrics toolbox

if (nargin ~= 3)
    error('Wrong # of arguments to beta_inv');
end

if any(any((a<=0)|(b<=0)))
   error('beta_inv parameter a or b is nonpositive');
end
if any(any(abs(2*p-1)>1))
   error('beta_inv: A probability should be 0<=p<=1');
end

x = a ./ (a+b);
dx = 1;
while any(any(abs(dx)>256*eps*max(x,1)))
   dx = (betainc(x,a,b) - p) ./ beta_pdf(x,a,b);
   x = x - dx;
   x = x + (dx - x) / 2 .* (x<0);
end

author: ambropo

repo: VAR-Toolbox

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


text:
function experiments = stack_vottir2015()

set_global_variable('bundle', 'http://www.cvl.isy.liu.se/research/datasets/ltir/version1.0/ltir_v1_0_8bit.zip');
set_global_variable('legacy_rasterization', true);
set_global_variable('bounded_overlap', false);

baseline.name = 'baseline';
baseline.converter = [];
baseline.type = 'supervised';
baseline.tags = {'camera_motion', 'illum_change', 'occlusion', 'size','motion', 'empty'};

baseline.parameters.repetitions = 15;
baseline.parameters.burnin = 10;
baseline.parameters.skip_initialize = 5;
baseline.parameters.failure_overlap = 0;

experiments = {baseline};

end

author: votchallenge

repo: toolkit-legacy


text:
%MAIN_singlePendulum.m
%
% This script runs a simple simulation of a single pendulum

% Use: EoM_Single_Pendulum to write the equations of motion

m = 1.0;  % (kg) pendulum mass
g = 9.81; % (m/s^2) gravity
l = 1.0; % (m) pendulum length

tSpan = [0,10]; %Simulation time interval

th0 = (pi/180)*(-178);  %Initial angle
w0 = 0;  %Initial angular rate
z0 = [th0;w0];

userFunc = @(t,z)singlePendulumRhs(t,z,g,l);

options = odeset(...
    'AbsTol',1e-8,...
    'RelTol',1e-8,...
    'Vectorized','on');

% Run the simulation!
sol = ode45(userFunc,tSpan,z0,options);

% Break apart solution for plotting
nPlot = 1000;
time = linspace(tSpan(1),tSpan(2),nPlot);
z = deval(sol,time); %Evaluate solution from ode45 at points in time
th = z(1,:);
w = z(2,:);
[energy, kinetic, potential] = singlePendulumEnergy(th,w,m,g,l);

% Plotting!

figure(1111);clf;

subplot(3,1,1);
plot(time,th,'k-','LineWidth',2);
xlabel('time (s)')
ylabel('angle (rad)');

subplot(3,1,2);
plot(time,w,'k-','LineWidth',2);
xlabel('time (s)')
ylabel('angle rate (rad/s)');

subplot(3,1,3);  hold on;
plot(time,energy,'k-','LineWidth',3);
plot(time,kinetic,'r-','LineWidth',2);
plot(time,potential,'b-','LineWidth',2);
xlabel('time (s)')
ylabel('energy (J)');
legend('total','kinetic','potential');

author: MatthewPeterKelly

repo: dscTutorials


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

%%
clear; close all; clc;

%% Syntax
% |[Vr]=rotate_vertices(V,R_mat);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [centers, mincenter, mindist, lower, computed] = anchors(firstcenter,k,data)
% choose k centers by the furthest-first method
% URL: http://cseweb.ucsd.edu/~elkan/fastkmeans.html

[n,dim] = size(data);
centers = zeros(k,dim);
lower = zeros(n,k);
mindist = Inf*ones(n,1);
mincenter = ones(n,1);
computed = 0;
centdist = zeros(k,k);

for j = 1:k
    if j == 1
        newcenter = firstcenter;
    else
        [maxradius,i] = max(mindist);
        newcenter = data(i,:);
    end

    centers(j,:) = newcenter;
    centdist(1:j-1,j) = calcdist(centers(1:j-1,:),newcenter);
    centdist(j,1:j-1) = centdist(1:j-1,j)';
    computed = computed + j-1;

    inplay = find(mindist > centdist(mincenter,j)/2);
    newdist = calcdist(data(inplay,:),newcenter);
    computed = computed + size(inplay,1);
    lower(inplay,j) = newdist;

    move = find(newdist < mindist(inplay));
    shift = inplay(move);
    mincenter(shift) = j;
    mindist(shift) = newdist(move);
end

author: adikhosla

repo: feature-extraction


text:
function area = polygonArea(poly, varargin)
%POLYGONAREA Compute the signed area of a polygon
%
%   A = polygonArea(POINTS);
%   Compute area of a polygon defined by POINTS. POINTS is a N-by-2 array
%   of double containing coordinates of vertices.
%
%   Vertices of the polygon are supposed to be oriented Counter-Clockwise
%   (CCW). In this case, the signed area is positive.
%   If vertices are oriented Clockwise (CW), the signed area is negative.
%
%   If polygon is self-crossing, the result is undefined.
%
%   Examples
%     % compute area of a simple shape
%     poly = [10 10;30 10;30 20;10 20];
%     area = polygonArea(poly)
%     area =
%         200
%
%     % compute area of CW polygon
%     area2 = polygonArea(poly(end:-1:1, :))
%     area2 =
%         -200
%
%     % Computes area of a paper hen
%     x = [0 10 20  0 -10 -20 -10 -10  0];
%     y = [0  0 10 10  20  10  10  0 -10];
%     poly = [x' y'];
%     area = polygonArea(poly)
%     area =
%        400
%
%     % Area of unit square with 25% hole
%     pccw = [0 0; 1 0; 1 1; 0 1];
%     pcw = pccw([1 4 3 2], :) * .5 + .25;
%     polygonArea ([pccw; nan(1,2); pcw])
%     ans =
%        0.75
%
%   References
%   algo adapted from P. Bourke web page
%   http://paulbourke.net/geometry/polygonmesh/
%
%   See also:
%   polygons2d, polygonCentroid, polygonSecondAreaMoments, triangleArea
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 05/05/2004.
%

%   HISTORY
%   25/04/2005: add support for multiple polygons
%   12/10/2007: update doc

%% Process special cases

% in case of polygon sets, computes the sum of polygon areas
if iscell(poly)
    area = 0;
    for i = 1:length(poly)
        area = area + polygonArea(poly{i});
    end
    return;
end

% check there are enough points
if size(poly, 1) < 2
    area = 0;
    return;
end

% case of polygons with holes -> computes the sum of areas
if any(isnan(poly))
    area = sum(polygonArea(splitPolygons(poly)));
    return;
end

%% Process single polygons or single rings

% extract coordinates
if nargin == 1
    % polygon given as N-by-2 array
    px = poly(:, 1);
    py = poly(:, 2);

elseif nargin == 2
    % poylgon given as two N-by-1 arrays
    px = poly;
    py = varargin{1};
end

% indices of next vertices
N = length(px);
iNext = [2:N 1];

% compute area (vectorized version)
area = sum(px .* py(iNext) - px(iNext) .* py) / 2;

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [p] = eeg_linfit(p)

% eeg_linfit - returns slope/intercept of linear fit
%
% USEAGE: [p] = eeg_linfit(p)
%
% p is the eeg_toolbox struct.  For this function, the fields
% required are:
%
%   p.volt.data - voltage data matrix (Npoints,Nelec)
%   p.volt.timeArray - voltage sample points (msec)
%
% Returns slope and intercept matrices (Npoints,Nelec)
% for linear fit of each electrode into fields of p:
%
%   p.volt.fitslope
%   p.volt.fitintercept
%
% The column vectors of slope/intercept hold constant
% values.  So, the linear fit data can be generated by:
%
% y = p.volt.fitslope .* p.volt.timeArray + p.volt.fitintercept;
%
% Where necessary, the p.volt.timeArray is replicated
% across columns for each electrode to enable this calculation
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:52 $

% Licence:  GNU GPL, no implied or express warranties
% History:  07/00, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% volt.timeArray is essentially a column vector, but
% operations in the eeg_toolbox may replicate across N columns
time = p.volt.timeArray(:,1);

slope = zeros(1,size(p.volt.data,2));
intercept = slope;

for elec = 1:size(p.volt.data,2),

    y = polyfit(time,p.volt.data(:,elec),1);
    slope(1,elec) = y(1);
    intercept(1,elec) = y(2);
end

% make sure that volt.timeArray is same size as volt.data
if ~isequal(size(p.volt.timeArray),size(p.volt.data)),
    p.volt.timeArray = repmat(p.volt.timeArray(:,1),1,size(p.volt.data,2));
end

p.volt.fitslope = repmat(slope,size(p.volt.data,1),1);
p.volt.fitintercept = repmat(intercept,size(p.volt.data,1),1);

return


text:
function ori = project2FundamentalRegion(ori,varargin)
% projects orientation to a fundamental region
%
% Syntax
%   ori = project2FundamentalRegion(ori,rot_ref)
%
% Input
%  ori     - @orientation
%  ori_ref - reference @rotation
%
% Output
%  ori     - @orientation
%  omega   - rotational angle to reference rotation
%

if isempty(ori.SS) || ismember(ori.SS.id, [1,2])
  ori = project2FundamentalRegion@quaternion(ori,ori.CS,varargin{:});
else
  if ori.antipodal, ap = {'antipodal'}; else, ap = {}; end
  CS = ori.CS; SS = ori.SS;
  ori = project2FundamentalRegion@quaternion(ori,CS,SS,ap{:},varargin{:});

  % ensure result is of type orientation
  if ~isa(ori,'orientation')
    ori = orientation(ori,CS,SS,ap{:});
  end

end

author: mtex-toolbox

repo: mtex

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


text:
function truncated_normal_b_pdf_values_test ( )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_B_PDF_VALUES_TEST tests TRUNCATED_NORMAL_B_PDF_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRUNCATED_NORMAL_B_PDF_VALUES_TEST:\n' );
  fprintf ( 1, '  TRUNCATED_NORMAL_B_PDF_VALUES stores values of\n' );
  fprintf ( 1, '  the PDF of the Normal Distribution truncated to (-oo,B].\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '            MU         SIGMA             B             X               FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, mu, sigma, b, x, fx ] = ...
      truncated_normal_b_pdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %12f  %12f  %12f  %24.16f\n', ...
      mu, sigma, b, x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fx2 = p19_fx2 ( x )

%*****************************************************************************80
%
%% P19_FX2 evaluates the second derivative of the function for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the point at which F is to be evaluated.
%
%    Output, real FX2(*), the value of the second derivative at X.
%
  arg = - ( 10.0 - 30.0 * x ).^2;

  fx2 = - 10000.0 * cos ( 100.0 * x ) ...
    + 14400.0 * exp ( arg ) * ( 10.0 - 30.0 * x ) / sqrt ( pi );

  return
end


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

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

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

      inside = 1;

    end

  end

  return
end


text:
% SYNTAX:
% yAvgStd = hmrE_RunAvgStd2_oldversion_nirs(yAvgRuns, ySum2Runs, mlActRuns, nTrialsRuns)
%
% UI NAME:
% Run_Average_Standard_Deviation
%
% DESCRIPTION:
% Calculate avearge HRF standard deviation of all runs for one subject.
%
% INPUTS:
% yAvgRuns:
% ySum2Runs:
% mlActRuns:
% nTrialsRuns:
%
% OUTPUTS:
% yAvgStdOut: the standard deviation across runs
%
% USAGE OPTIONS:
% Run_Average_Standard_Deviation_on_Concentration_Data:  dcAvgStd  = hmrE_RunAvgStd2_oldversion_nirs(dcAvgRuns, dcSum2Runs, mlActRuns, nTrialsRuns)
% Run_Average_Standard_Deviation_on_Delta_OD_Data:       dodAvgStd = hmrE_RunAvgStd2_oldversion_nirs(dodAvgRuns, dodSum2Runs, mlActRuns, nTrialsRuns)
%

function yAvgStdOut = hmrE_RunAvgStd2_oldversion_nirs(yAvgRuns, ySum2Runs, mlActRuns, nTrialsRuns)

yAvgStdOut = DataClass().empty();

nDataBlks = length(yAvgRuns{1});
nTrials_tot = cell(nDataBlks,1);

for iBlk = 1:nDataBlks

    grp1 = [];

    for iRun = 1:length(yAvgRuns)

        yAvgStdOut(iBlk) = DataClass();

        yAvg    = yAvgRuns{iRun}(iBlk).GetDataTimeSeries('reshape');
        if isempty(yAvg)
            continue;
        end

        yAvgStd = zeros(size(yAvg,1), size(yAvg,2), size(yAvg,3), size(yAvg,4));
        ySum2   = ySum2Runs{iRun}(iBlk).GetDataTimeSeries('reshape');
        tHRF    = yAvgRuns{iRun}(iBlk).GetTime();
        nTrials = nTrialsRuns{iRun}{iBlk};
        if isempty(mlActRuns{iRun})
            mlActRuns{iRun} = cell(length(nDataBlks),1);
        end

        %
        datatype  = yAvgRuns{iRun}(iBlk).GetDataTypeLabel();
        if strncmp(datatype{1}, 'HRF Hb', length('HRF Hb'))
            ml    = yAvgRuns{iRun}(iBlk).GetMeasListSrcDetPairs('reshape');
        elseif strcmp(datatype{1}, 'HRF dOD')
            ml    = yAvgRuns{iRun}(iBlk).GetMeasList('reshape');
        end
        if isempty(mlActRuns{iRun}{iBlk})
            mlActRuns{iRun}{iBlk} = ones(size(ml,1),1);
        end
        mlAct = mlActRuns{iRun}{iBlk}(1:size(ml,1));

        nCond = size(nTrials,2);
        yAvgStdOut(iBlk).SetTime(tHRF);

        if strcmp(datatype{1}, 'HRF dOD')

            if isempty(grp1)
                grp1 = zeros(size(yAvg,1), size(yAvg,2), nCond);
                grp1Sum2 = zeros(size(yAvg,1), size(yAvg,2), nCond);
                nTrials_tot{iBlk} = zeros(size(yAvg,2), nCond);
            end

            lstChInc = find(mlAct==1);
            for iC = 1:nCond
                nT = nTrials(iC);
                if nT>0
                    if iRun==1
                        grp1(:,lstChInc,iC) = yAvg(:,lstChInc,iC) * nT;
                        grp1Sum2(:,lstChInc,iC) = ySum2(:,lstChInc,iC);
                        nTrials_tot{iBlk}(lstChInc,iC) = nT;
                    else
                        for iCh=1:length(lstChInc) %size(yAvg,2)
                            % Make sure 3rd arg to interp1 is column vector to guarauntee interp1 output is column vector
                            % which matches grp1 dimensions when adding the two.
                            grp1(:,lstChInc(iCh),iC) = grp1(:,lstChInc(iCh),iC) + interp1(tHRF,yAvg(:,lstChInc(iCh),iC),tHRF(:)) * nT;
                            grp1Sum2(:,lstChInc(iCh),iC) = grp1Sum2(:,lstChInc(iCh),iC) + interp1(tHRF,ySum2(:,lstChInc(iCh),iC),tHRF(:));
                            nTrials_tot{iBlk}(lstChInc(iCh),iC) = nTrials_tot{iBlk}(lstChInc(iCh),iC) + nT;
                        end
                    end
                end
            end

            if ~isempty(grp1)
                for iC = 1:size(grp1,3)
                    for iCh = 1:size(grp1,2)
                        yAvg(:,iCh,iC) = grp1(:,iCh,iC) / nTrials_tot{iBlk}(iCh,iC);

                        % We want to distinguish between no trials and 1 trial:
                        % If there are no trials, we have no HRF data and no std which
                        % the first case will calculate as opposed to one trial (2nd case)
                        % where we have all zeros.
                        if(nTrials_tot{iBlk}(iCh,iC)~=1)
                            yAvgStd(:,iCh,iC) = ( (1/(nTrials_tot{iBlk}(iCh,iC)-1))*grp1Sum2(:,iCh,iC) - (nTrials_tot{iBlk}(iCh,iC)/(nTrials_tot{iBlk}(iCh,iC)-1))*(grp1(:,iCh,iC) / nTrials_tot{iBlk}(iCh,iC)).^2).^0.5 ;
                        else
                            yAvgStd(:,iCh,iC) = zeros(size(grp1Sum2(:,iCh,iC)));
                        end

                        %%%% Snirf stuff: Once we get to the last run, we've accumulated our averages.
                        %%%% Now we can set channel descriptors for avg and standard deviation
                        if iRun == length(yAvgRuns)
                            yAvgStdOut(iBlk).AddChannelDod(ml(iCh,1), ml(iCh,2), ml(iCh,4), iC);
                        end
                    end

                    %%%% Snirf stuff: Once we get to the last run, we've accumulated our averages.
                    %%%% Now we can set channel descriptors for avg and standard deviation
                    if iRun == length(yAvgRuns)
                        yAvgStdOut(iBlk).AppendDataTimeSeries(yAvgStd(:,:,iC));
                    end
                end
            end

        elseif strncmp(datatype{1}, 'HRF Hb', length('HRF Hb'))

            if isempty(grp1)
                grp1 = zeros(size(yAvg,1), size(yAvg,2), size(yAvg,3), nCond);
                grp1Sum2 = zeros(size(yAvg,1), size(yAvg,2), size(yAvg,3), nCond);
                nTrials_tot{iBlk} = zeros(size(yAvg,3), nCond);
            end

            lstChInc = find(mlAct==1);
            for iC = 1:1:nCond
                nT = nTrials(iC);
                if nT>0
                    if iRun==1
                        grp1(:,:,lstChInc,iC) = yAvg(:,:,lstChInc,iC) * nT;
                        grp1Sum2(:,:,lstChInc,iC) = ySum2(:,:,lstChInc,iC);
                        nTrials_tot{iBlk}(lstChInc,iC) = nT;
                    else
                        for iCh=1:length(lstChInc) %size(yAvg,3)
                            for iHb=1:size(yAvg,2)
                                % Make sure 3rd arg to interp1 is column vector to guarauntee interp1 output is column vector
                                % which matches grp1 dimensions when adding the two.
                                grp1(:,iHb,lstChInc(iCh),iC) = grp1(:,iHb,lstChInc(iCh),iC) + interp1(tHRF,yAvg(:,iHb,lstChInc(iCh),iC),tHRF(:)) * nT;
                                grp1Sum2(:,iHb,lstChInc(iCh),iC) = grp1Sum2(:,iHb,lstChInc(iCh),iC) + interp1(tHRF,ySum2(:,iHb,lstChInc(iCh),iC),tHRF(:));
                            end
                            nTrials_tot{iBlk}(lstChInc(iCh),iC) = nTrials_tot{iBlk}(lstChInc(iCh),iC) + nT;
                        end
                    end
                end
            end

            if ~isempty(grp1)
                for iC = 1:size(grp1,4)
                    for iCh = 1:size(grp1,3)
                        yAvg(:,:,iCh,iC) = grp1(:,:,iCh,iC) / nTrials_tot{iBlk}(iCh,iC);

                        % We want to distinguish between no trials and 1 trial:
                        % If there are no trials, we have no HRF data and no std which
                        % the first case will calculate as opposed to one trial (2nd case)
                        % where we have all zeros.
                        if(nTrials_tot{iBlk}(iCh,iC)~=1)
                            yAvgStd(:,:,iCh,iC) = ( (1/(nTrials_tot{iBlk}(iCh,iC)-1))*grp1Sum2(:,:,iCh,iC) - (nTrials_tot{iBlk}(iCh,iC)/(nTrials_tot{iBlk}(iCh,iC)-1))*(grp1(:,:,iCh,iC) / nTrials_tot{iBlk}(iCh,iC)).^2).^0.5 ;
                        else
                            yAvgStd(:,:,iCh,iC) = zeros(size(grp1Sum2(:,:,iCh,iC)));
                        end

                        %%%% Snirf stuff: Once we get to the last run, we've accumulated our averages.
                        %%%% Now we can set channel descriptors for avg and standard deviation
                        if iRun == length(yAvgRuns)
                            yAvgStdOut(iBlk).AddChannelHbO(ml(iCh,1), ml(iCh,2), iC);
                            yAvgStdOut(iBlk).AddChannelHbR(ml(iCh,1), ml(iCh,2), iC);
                            yAvgStdOut(iBlk).AddChannelHbT(ml(iCh,1), ml(iCh,2), iC);
                        end
                    end

                    %%%% Snirf stuff: Once we get to the last run, we've accumulated our averages.
                    %%%% Now we can set channel descriptors for avg and standard deviation
                    if iRun == length(yAvgRuns)
                        yAvgStdOut(iBlk).AppendDataTimeSeries(yAvgStd(:,:,:,iC));
                    end
                end
            end
        end
    end
end

author: BUNPC


text:
function y = sum(varargin)
% Fixes special case sum([],1), sum([],1,...), sum('',1), sum('',1,...)
% Fixes special case sum(sym([])), sum(sym([]),...) (without specifying dimension)
% Adds 'omitnan'  functionality
if nargin>1 && isequal(size(varargin{1}), [0 0]) && isequal(varargin{2}, 1)
    y = zeros(1,0);
elseif ( nargin==1 && isa(varargin{1}, 'sym') && isequal(size(varargin{1}), [0 0]) ) || ...
        ( nargin>1  && isa(varargin{1}, 'sym') && isequal(size(varargin{1}), [0 0]) && ~isnumeric(varargin{2}) )
    y = 0;
elseif nargin>1 && strcmp(varargin{end},'omitnan')
    y = nansum(varargin{1:end-1});
else
    if ~isa(varargin{1},'sym')
        y = builtin('sum', varargin{:});
    else
        y = builtin('@sym/sum', varargin{:});
    end
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function v = gradient(f)
%GRADIENT Gradient of a BALLFUN in cartesian coordinates.
%   GRADIENT(F) is the gradient of the BALLFUN F expressed in
%   cartesian coordinates.
%
% See also DIV, CURL

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

Vx = diff(f, 1);
Vy = diff(f, 2);
Vz = diff(f, 3);
v = ballfunv(Vx, Vy, Vz);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function drawSecPot(dietPath, rxnmax, startMet, levels, cutoff)
% Draws part of the community structure when maximized with the given
% reaction
%
% INPUTS:
%	dietPath          char with path of the model
%   rxnmax            char with name of the rxn that gets maximized with fastFVA
%   startMet          char with starting metabolite for visualizing
%   levels            number with levels that are connected to the starting
%                     metabolite
%   cutoff            number of fluxes that are considered for drawing,
%                     only fluxes>abs(cutoff) are considered

global CBT_LP_SOLVER
if isempty(CBT_LP_SOLVER)
    initCobraToolbox
end

microbiota_model=load(dietPath);
modelF=fieldnames(microbiota_model);
model=microbiota_model.(modelF{1});

% parameters for fastFVA
cpxControl.PARALLELMODE = 1;
cpxControl.THREADS = 1;
cpxControl.AUXROOTTHREADS = 2;

% calculating max secretion potential
%(note: if diet constraint present max. secretion flux = max(EX_rxn[fe]) - min(EX_rxn[d]))
[minFlux, maxFlux, optsol, ret, fbasol, fvamin, fvamax]  = fastFVA(model,99.99,'max',{},rxnmax,'A',cpxControl);

% all the metabolites and reactions connected to the starting
startMetnr = find(string(model.mets) == char(startMet));

[a,b] = inorder(startMetnr, model, [], [], levels(1));

% names of the reactions that are of interest
rxoi = model.rxns(b);

%  filter reactions and fluxes that are greater than abs(cutoff)
rxngz = [];
flxgz = [];
n = length(rxoi)
for i=1:n
    progress = i/n;
    sprintf('Progress %f ', progress)
    if(abs(fvamax(find(string(model.rxns)== char(rxoi(i)), 1))) > cutoff)
        rxngz = [rxngz, rxoi(i)];
        %flxgz = [flxgz, fvamax(find(string(model.rxns)== char(rxoi(i)), 1))];
    end
end

%flxgz = flxgz';

diffmets = setdiff(model.mets, model.mets(a));
[Involved_mets, Dead_ends] = draw_by_rxn(model, rxngz, 'true', 'struc', {''}, diffmets, fvamax);

author: opencobra

repo: cobratoolbox


text:
function overlaps = bbox_overlap(boxes1,boxes2)
% Copyright (C) 2016 Hakan Bilen.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).
x11 = boxes1(:,1);
y11 = boxes1(:,2);
x12 = boxes1(:,3);
y12 = boxes1(:,4);

x21 = boxes2(:,1);
y21 = boxes2(:,2);
x22 = boxes2(:,3);
y22 = boxes2(:,4);

N1 = size(boxes1,1);
N2 = size(boxes2,1);

area1 = (x12-x11+1) .* (y12-y11+1);
area2 = (x22-x21+1) .* (y22-y21+1);

overlaps = zeros(N1,N2);

for i=1:N1

  xx1 = max(x11(i), x21);
  yy1 = max(y11(i), y21);
  xx2 = min(x12(i), x22);
  yy2 = min(y12(i), y22);

  w = max(0.0, xx2-xx1+1);
  h = max(0.0, yy2-yy1+1);

  inter = w.*h;
  overlaps(i,:) = inter ./ (area1(i) + area2 - inter);
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function hole_num = p06_hole_num ( )

%*****************************************************************************80
%
%% P06_HOLE_NUM counts the holes in problem 06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer HOLE_NUM, the number of holes.
%
  hole_num = 1;

  return
end


text:
function [hold_state, cax, next] = axes(majors,ticks)
% TERNARY.AXES create ternary axis
%   [hold_state, cax, next]  = TERNARY.AXES(MAJORS) creates a ternary axis system using the system
%   defaults and with MAJORS major tickmarks with percentage labels at tick
%   marks. Returns the hold state of the plot.
%
% [hold_state, cax, next]  = TERNARY.AXES(MAJORS,'fractions')  returns the
% axes with fraction of unit at tick marks.

% Author: Carl Sandrock 20050211
% Modifications: FVA:

% To Do

% Modifications

% Modifiers
% (CS) Carl Sandrock

percentage = 1;
if nargin > 1
    switch ticks
        case 'fraction'
            percentage = 0;
        otherwise
            percentage = 1;
    end
end

%TODO: Get a better way of offsetting the labels
xoffset = 0.01;
yoffset = 0.02;

% get hold state
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;

% get x-axis text color so grid is in same color
tc = get(cax,'xcolor');
ls = get(cax,'gridlinestyle');

% Hold on to current Text defaults, reset them to the
% Axes' font attributes so tick marks use them.
fAngle  = get(cax, 'DefaultTextFontAngle');
fName   = get(cax, 'DefaultTextFontName');
fSize   = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
fUnits  = get(cax, 'DefaultTextUnits');

set(cax, 'DefaultTextFontAngle',  get(cax, 'FontAngle'), ...
    'DefaultTextFontName',   get(cax, 'FontName'), ...
    'DefaultTextFontSize',   get(cax, 'FontSize'), ...
    'DefaultTextFontWeight', get(cax, 'FontWeight'), ...
    'DefaultTextUnits','data')

% only do grids if hold is off
if ~hold_state
	%plot axis lines
	hold on;
	plot ([0 1 0.5 0],[0 0 sin(1/3*pi) 0], 'color', tc, 'linewidth',1,...
                   'handlevisibility','off');
	set(gca, 'visible', 'off');

    % plot background if necessary
    if ~ischar(get(cax,'color')),
       patch('xdata', [0 1 0.5 0], 'ydata', [0 0 sin(1/3*pi) 0], ...
             'edgecolor',tc,'facecolor',get(gca,'color'),...
             'handlevisibility','off');
    end

	% Generate labels
	majorticks = linspace(0, 1, majors + 1);
	majorticks = majorticks(1:end-1);
    if percentage
        labels = num2str(majorticks'*100);
    else
        labels = num2str(majorticks','%1.1g');%Doesn't work! I don't want it to put a leading zero, that's all!
    end

    zerocomp = zeros(size(majorticks)); % represents zero composition

	% Plot right labels (no c - only b a)
	[lxc, lyc] = ternary.coords(1-majorticks, majorticks, zerocomp);
	text(lxc, lyc, [repmat('  ', length(labels), 1) labels]);

	% Plot bottom labels (no b - only a c)
	[lxb, lyb] = ternary.coords(majorticks, zerocomp, 1-majorticks); % fB = 1-fA
	text(lxb, lyb, labels, 'VerticalAlignment', 'Top');

	% Plot left labels (no a, only c b)
	[lxa, lya] = ternary.coords(zerocomp, 1-majorticks, majorticks);
	%text(lxa-xoffset, lya, labels, 'HorizontalAlignment','right');
	text(lxa-xoffset, lya, labels, 'HorizontalAlignment','right');

	nlabels = length(labels)-1;
	for i = 1:nlabels
        plot([lxa(i+1) lxb(nlabels - i + 2)], [lya(i+1) lyb(nlabels - i + 2)], ls, 'color', tc, 'linewidth',1,...
           'handlevisibility','off');
        plot([lxb(i+1) lxc(nlabels - i + 2)], [lyb(i+1) lyc(nlabels - i + 2)], ls, 'color', tc, 'linewidth',1,...
           'handlevisibility','off');
        plot([lxc(i+1) lxa(nlabels - i + 2)], [lyc(i+1) lya(nlabels - i + 2)], ls, 'color', tc, 'linewidth',1,...
           'handlevisibility','off');
    end
end%if ~hold_state

% Reset defaults
set(cax, 'DefaultTextFontAngle', fAngle , ...
    'DefaultTextFontName',   fName , ...
    'DefaultTextFontSize',   fSize, ...
    'DefaultTextFontWeight', fWeight, ...
    'DefaultTextUnits', fUnits );
return%[hold_state, cax, next]

author: Sable

repo: mcbench-benchmarks


text:
%klists ' '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros lists.pane file
%
% Parameters:
%
% Example: klists( {})
%
% Khoros helpfile follows below:

function varargout = klists(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = klists(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={};
maxval={};
minval={};
istoggle=[];
was_set=istoggle * 0;
paramtype={};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'lists" '],Inputs,narglist);


text:
function x = prox_cvx(v, lambda, f, l, u)
% PROX_CVX    The proximal operator of a generic function.
%
%   prox_cvx(v,lambda,f) is the proximal operator of f with parameter lambda
%   evaluated at v. Here, f is a closure, possibly described by an anonymous
%   function. It is also possible to provide lower and upper bounds, so
%   prox_cvx(v,f,l,u) is the prox operator of f + I_[l,u].
%
%   For example,
%
%     prox_cvx(v, lambda, (@(x) norm(x,1))
%
%   is equivalent to the soft thresholding operator.
%
%   WARNING: This is a *very* inefficient way of evaluating a proximal
%   operator. This function is mainly for rapid prototyping and
%   testing custom implementations of particular proximal operators.
%   It may also be useful when suffering from extreme laziness.

    if ~exist('l', 'var') || isnan(l)
        l = -Inf;
    end

    if ~exist('u', 'var') || isnan(u)
        u = Inf;
    end

    L = 1/(2*lambda);
    [m n] = size(v);

    if min(m,n) == 1
        cvx_begin quiet
            variable x(max(m,n))
            minimize(f(x) + L*sum_square(x - v))
            subject to
                v >= l;
                v <= u;
        cvx_end
    elseif m == n && m > 1
        cvx_begin quiet
            variable x(n,n) symmetric
            minimize(f(x) + L*pow_pos(norm(x - v, 'fro'), 2))
            subject to
                x == semidefinite(n)
        cvx_end
    else
        x = nan(size(v));
    end
end

author: cvxgrp

repo: proximal

save_path: github-repos/MATLAB/cvxgrp-proximal


text:
function [ x, seed ] = gompertz_sample ( a, b, seed )

%*****************************************************************************80
%
%% GOMPERTZ_SAMPLE samples the Gompertz PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    1 < A, 0 < B.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ cdf, seed ] = r8_uniform_01 ( seed );

  x = gompertz_cdf_inv ( cdf, a, b );

  return
end


text:
function [mh,mw]=tropmapf(time,pos,azel)

if pos(3)<-1000||pos(3)>20000
    mh=0; mw=0; return;
end

% hydro-ave-a,b,c, hydro-amp-a,b,c, wet-a,b,c at latitude 15,30,45,60,75
coef=[ 1.2769934E-3, 1.2683230E-3, 1.2465397E-3, 1.2196049E-3, 1.2045996E-3;...
       2.9153695E-3, 2.9152299E-3, 2.9288445E-3, 2.9022565E-3, 2.9024912E-3;...
       62.610505E-3, 62.837393E-3, 63.721774E-3, 63.824265E-3, 64.258455E-3;...

       0.0000000E-0, 1.2709626E-5, 2.6523662E-5, 3.4000452E-5, 4.1202191E-5;...
       0.0000000E-0, 2.1414979E-5, 3.0160779E-5, 7.2562722E-5, 11.723375E-5;...
       0.0000000E-0, 9.0128400E-5, 4.3497037E-5, 84.795348E-5, 170.37206E-5;...

       5.8021897E-4, 5.6794847E-4, 5.8118019E-4, 5.9727542E-4, 6.1641693E-4;...
       1.4275268E-3, 1.5138625E-3, 1.4572752E-3, 1.5007428E-3, 1.7599082E-3;...
       4.3472961E-2, 4.6729510E-2, 4.3908931E-2, 4.4626982E-2, 5.4736038E-2];

% height correction
aht=[2.53E-5, 5.49E-3, 1.14E-3];
ah=zeros(3,1); aw=zeros(3,1);

if azel(2)<=0,mh=0; mw=0; return;end

lat=pos(1)*180/pi;
if lat<0,yy=0.5;else,yy=0;end
y=(time2doy(time)-28.0)/365.25+yy;
cosy=cos(2.0*pi*y);
lat=abs(lat);

for i=1:3
    ah(i)=interpc(coef(i,:),lat)-interpc(coef(i+3,:),lat)*cosy;
    aw(i)=interpc(coef(i+6,:),lat);
%     j=fix(lat/15);
%     ah(i)=(coef(i,j)+(coef(i,j+1)-coef(i,j))*(lat/15-j))-...
%           (coef(i+3,j)+(coef(i+3,j+1)-coef(i+3,j))*(lat/15-j))*cosy;
%     aw(i)=(coef(i+6,j)+(coef(i+6,j+1)-coef(i+6,j))*(lat/15-j));
end

el=azel(2);
sinel=sin(el);

mh1=(1+ah(1)/(1+ah(2)/(1+ah(3))))/(sinel+(ah(1)/(sinel+ah(2)/(sinel+ah(3)))));
mh2=(1+aht(1)/(1+aht(2)/(1+aht(3))))/(sinel+(aht(1)/(sinel+aht(2)/(sinel+aht(3)))));
mh=mh1+(1/sinel-mh2)*pos(3)/1000;

mw=(1+aw(1)/(1+aw(2)/(1+aw(3))))/(sinel+(aw(1)/(sinel+ah(2)/(sinel+aw(3)))));

return

repo: GINav


text:
function[kappa,lambda,theta,phi,alpha,beta]=ellparams(varargin)
%ELLPARAMS  Ellipse parameters of a modulated bivariate or trivariate oscillation.
%
%   [KAPPA,LAMBDA,THETA,PHI]=ELLPARAMS(X,Y) where X and Y are analytic
%   signals, returns the parameters of the complex-valued signal
%   Z=REAL(X)+i REAL(Y), expressed as a modulated ellipse.
%
%   Here KAPPA is the RMS ellipse amplitude, LAMBDA is the linearity,
%   THETA is the orientation, and PHI is the instantaneous orbital phase.
%
%   ELLPARAMS(M), where M is matrix with two columns, also works.
%
%   See Lilly and Gascard (2006) and Lilly and Olhede (2010a) for details.
%
%   ELLPARAMS is inverted by ELLSIG, which returns the X and Y signals
%   given the ellipse parameters.
%
%   ELLPARAMS(...,DIM) performs the analysis with time running along
%   dimension DIM, as opposed to the default behavior of DIM=1.
%
%   ELLPARAMS also works if X and Y are cell arrays with each cell holding
%   a different analytic signal. The output will then also be cell arrays.
%   _______________________________________________________________________
%
%   Trivariate signals
%
%   ELLPARAMS also works for trivariate signals, which can be expressed as
%   a modulated ellipse in three dimensions.
%
%   [KAPPA,LAMBDA,THETA,PHI,ALPHA,BETA]=ELLPARAMS(X,Y,Z), where X, Y, and Z
%   are all analytic signals, also returns the zenith angle ALPHA and the
%   azimuth angle BETA in addition to the other ellipse parameters.
%
%   ELLPARAMS(M), where M is matrix with three columns, also works.
%
%   See Lilly (2010) for details on the trivariate case.
%   __________________________________________________________________
%
%   'ellparams --t' runs a test.
%
%   See also ELLSIG, ELLBAND, ELLDIFF, ELLVEL, ELLRAD, KL2AB, AB2KL.
%
%   Usage: [kappa,lambda,theta,phi]=ellparams(x,y);
%          [kappa,lambda,theta,phi]=ellparams(x,y,dim);
%          [kappa,lambda,theta,phi,alpha,beta]=ellparams(x,y,z);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2009--2018 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
    ellparams_test,normvect_test,return
end

[na,k,l,theta,phi,alpha,beta]=vempty;

if ~iscell(varargin{end})&&length(varargin{end})==1
    dim=varargin{end};
    varargin=varargin(1:end-1);
else
    dim=1;
end

[z,kappa,lambda,theta,phi,alpha,beta]=vempty;
if ~isempty(varargin{1})
    if ~iscell(varargin{1})
        x=varargin{1};
        y=varargin{2};
        if length(varargin)==3
            z=varargin{3};
        end
        [kappa,lambda,theta,phi,alpha,beta]=ellparams_one(x,y,z,dim);
    else
        x=varargin{1};
        y=varargin{2};
        for i=1:length(x)
            if ~isempty(x{i})
                if length(varargin)==2
                    [kappa{i,1},lambda{i,1},theta{i,1},phi{i,1}]=ellparams_one(x{i},y{i},[],dim);
                else
                    z=varargin{3};
                    [kappa{i,1},lambda{i,1},theta{i,1},phi{i,1},alpha{i,1},beta{i,1}]=...
                        ellparams_one(x{i},x{i},z{3}{i},dim);
                end
            else
                [kappa{i,1},lambda{i,1},theta{i,1},phi{i,1},alpha{i,1},beta{i,1}]=vempty;
            end
        end
    end
end

function[kappa,lambda,theta,phi,alpha,beta]=ellparams_one(x,y,z,dim)

[alpha,beta]=vempty;

if isempty(z)
    [kappa,lambda,theta,phi]=ellconv_xy2kl(abs(x),abs(y),angle(x),angle(y),dim);
else
    [nx,ny,nz]=normvect(x,y,z);
    warning('off','MATLAB:log:logOfZero')
    alpha=imag(log(sqrt(-1)*nx-ny));
    warning('on','MATLAB:log:logOfZero')
    beta=imag(log(nz+sqrt(-1)*sqrt(nx.^2+ny.^2)));
    [x,y,z]=vectmult(jmat3(-alpha,3),x,y,z);
    [x,y,z]=vectmult(jmat3(-beta,1),x,y,z);
    [kappa,lambda,theta,phi]=ellconv_xy2kl(abs(x),abs(y),angle(x),angle(y),dim);
end

function[kappa,lambda,theta,phi]=ellconv_xy2kl(X,Y,phix,phiy,dim)
%phia=double(phix+phiy+pi/2)/2;
%phid=double(phix-phiy-pi/2)/2;
%P=double(frac(1,2)*sqrt(squared(X)+squared(Y)+2.*X.*Y.*cos(2*phid)));
%N=double(frac(1,2)*sqrt(squared(X)+squared(Y)-2.*X.*Y.*cos(2*phid)));

phia=(phix+phiy+pi/2)/2;
phid=(phix-phiy-pi/2)/2;

P=frac(1,2)*sqrt(squared(X)+squared(Y)+2.*X.*Y.*cos(2*phid));
N=frac(1,2)*sqrt(squared(X)+squared(Y)-2.*X.*Y.*cos(2*phid));

phip=unwrap(phia+imlog(X.*rot(phid)+Y.*rot(-phid)),[],dim);
phin=unwrap(phia+imlog(X.*rot(phid)-Y.*rot(-phid)),[],dim);

kappa=sqrt(P.^2+N.^2);
lambda=frac(2*P.*N.*sign(P-N),P.^2+N.^2);

%For vanishing linearity, put in very small number to have sign information
lambda(lambda==0)=sign(P(lambda==0)-N(lambda==0))*(1e-10);

theta=phip/2-phin/2;
phi=  phip/2+phin/2;

theta=unwrap(theta,[],dim);
phi=unwrap(phi,[],dim);

lambda=real(lambda);

function[nx,ny,nz]=normvect(x,y,z)
%NORMVECT  Unit normal vector to the ellipse plane in three dimensions.
%
%   [NX,NY,NZ]=NORMVECT(X,Y,Z) returns the three components of the unit
%   normal vector to the plane containing the ellipse specified by the
%   three complex-valued arrays X, Y, and Z.
%
%   In vector notation the normal vector is defined as N=IMAG(X) x REAL(X),
%   where ``x'' is the vector cross product, and the unit normal is N/||N||.
%
%   The input arrays and output arrays are all the same size.
%
%   See Lilly (2010) for details.
%
%   'normvect --t' runs a test.
%
%   Usage: [nx,ny,nz]=normvect(x,y,z);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2010 J.M. Lilly --- type 'help jlab_license' for details

nx= (imag(y).*real(z)-imag(z).*real(y));
ny=-(imag(x).*real(z)-imag(z).*real(x));
nz= (imag(x).*real(y)-imag(y).*real(x));

denom=sqrt(nx.^2+ny.^2+nz.^2);
nx=frac(nx,denom);
ny=frac(ny,denom);
nz=frac(nz,denom);

function[]=normvect_test

load solomon
use solomon

[x,y,z]=anatrans(x./1e4,y./1e4,z./1e4);
[nx,ny,nz]=normvect(x,y,z);

dot=nx.*x+ny.*y+nz.*z;

reporttest('NORMVECT parallel part of X_+ vanishes, Solomon Islands',allall(abs(dot)<1e-10))

%Choose central part where signal is elliptical
vindex(x,y,z,7000:10000,1);

[ax,omx,upx]=instmom(x);
[ay,omy,upy]=instmom(y);
[az,omz,upz]=instmom(z);

dx=x.*(upx+sqrt(-1)*omx);
dy=y.*(upy+sqrt(-1)*omy);
dz=z.*(upz+sqrt(-1)*omz);

[nx,ny,nz]=normvect(x,y,z);

dot=nx.*dx+ny.*dy+nz.*dz; %Parallel part of derivative

[dnx,dny,dnz]=vdiff(nx,ny,nz,1);

dot2=-(dnx.*x+dny.*y+dnz.*z);

err=abs(dot-dot2).^2./abs(dot).^2;
err=flipud(sort(err));
err=err(60:end);

reporttest('NORMVECT derivative of parallel part matches, Solomon Islands (removing worst outliers)',allall(err<0.05))

function[]=ellparams_test

t=(0:1:925)';
kappa=3*exp(2*0.393*(t/1000-1));
lambda=0.5+0*kappa;
phi=(t/1000*5)*2*pi;
theta=pi/4+phi./14.45;

beta=pi/6+phi./14.45*lambda(1);
alpha=pi/6-phi./14.45*2*lambda(1)*sqrt(2);

[x,y,z]=ellsig(kappa,lambda,theta,phi,alpha,beta);
[kappa2,lambda2,theta2,phi2,alpha2,beta2]=ellparams(x,y,z);

x1=[kappa lambda theta phi alpha beta];
x2=[kappa2 lambda2 theta2 phi2 alpha2 beta2];
reporttest('ELLPARAMS rapidly changing trivariate ellipse',aresame(x1,x2,1e-8))

% %/*******************************************************************
% %Flip ellipse parameters if unit normal is pointing radially inwards
% %Components of unit normal vector to surface of earth
% [nx,ny,nz]=normvect(xr,yr,zr);
%
% %Replicate x, y, and z along columns
% xmat=vrep(x,size(nx,2),2)./radearth;
% ymat=vrep(y,size(ny,2),2)./radearth;
% zmat=vrep(z,size(nz,2),2)./radearth;
%
% %Projection of normal vector to plane onto normal to sphere
% proj=xmat.*nx+ymat.*ny+zmat.*nz;
%
% bool=(proj<0);
% theta(bool)=-theta(bool);
% lambda(bool)=-lambda(bool);
% beta(bool)=pi+beta(bool);
% nx(bool)=-nx(bool);
% ny(bool)=-ny(bool);
% nz(bool)=-nz(bool);
%
% if length(find(bool))>0
%     [xr2,yr2,zr2]=ellsig(kappa,lambda,theta,phi,alpha,beta);
%     tol=1e-6;
%     reporttest('ELLIPSEXTRACT adjustment for sign of normal vector',aresame(xr2,xr,tol)&&aresame(yr2,yr,tol)&&aresame(zr2,zr,tol))
% end
%
% %dev=1-abs(proj);
% %figure,plot(dev)
% [latn,lonn]=xyz2latlon(nx*radearth,ny*radearth,nz*radearth);
% %\*******************************************************************

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
clc;clear all;close all;
%***********************************************%
% This code runs on the VeRi-776 dataset.
% This code uses post fuse the results of CNN + CN + SIFT.
% We use the mAP and hit-1/5 rate as evaluation
%***********************************************%
%% load query hist and test hist
CNN_featureLen = 431;
CN_featureLen = 5600;
SIFT_featureLen = 10000;
nQuery = 1678;
nTest = 11579;

dist_SIFT = zeros(nTest, nQuery, 'double');
dist_CN = zeros(nTest, nQuery, 'double');
dist_CNN = zeros(nTest, nQuery, 'double');

fidin = fopen('dist_SIFT_776.txt');
for i = 1:nTest
    dist_line = fgetl(fidin);
    dist_line = str2num(dist_line);
    dist_SIFT(i, 1:length(dist_line)) = dist_line;
end
fclose(fidin);
fidin = fopen('dist_CN_776.txt');
for i = 1:nTest
    dist_line = fgetl(fidin);
    dist_line = str2num(dist_line);
    dist_CN(i, 1:length(dist_line)) = dist_line;
end
fclose(fidin);
fidin = fopen('dist_CNN_776.txt');
for i = 1:nTest
    dist_line = fgetl(fidin);
    dist_line = str2num(dist_line);
    dist_CNN(i, 1:length(dist_line)) = dist_line;
end
fclose(fidin);

dist = dist_CNN*0.7 + dist_CN*0.2 + dist_SIFT*0.1;

%% load ground truth index
maxgt = 256;
gt_index =  zeros(nQuery, maxgt);
fidin = fopen('gt_index_776.txt');

for i = 1:nQuery
    gt_index_line = fgetl(fidin);
    gt_line = str2num(gt_index_line);
    for j = 1:size(gt_line, 2)
       gt_index(i, j) = gt_line(j);
    end
end

maxjk = 256;
jk_index = zeros(nQuery, maxjk);
fidin = fopen('jk_index_776.txt');
for i = 1:nQuery
    jk_index_line = fgetl(fidin);
    jk_line = str2num(jk_index_line);
    for j = 1:size(jk_line, 2)
       jk_index(i, j) = jk_line(j);
    end
end
%% search the database and calcuate re-id accuracy
ap = zeros(nQuery, 1); % average precision
CMC = zeros(nQuery, nTest);
r1 = 0; % rank 1 precision with single query
for k = 1:nQuery
    k
    % load groud truth for each query (good and junk)
      good_index = reshape(gt_index(k,:), 1, []);
      good_index = good_index(good_index ~= 0);
      junk_index = reshape(jk_index(k,:), 1, []);
      junk_index = junk_index(junk_index ~= 0);
    tic
    %score = dist(:, k);
    %score = dist_null(:, k);
    score = dist_LOMO(:, k);
    [~, index] = sort(score, 'ascend');  % single query
    [ap(k), CMC(k, :)] = compute_AP(good_index, junk_index, index);% compute AP for single query

end
CMC = mean(CMC);
%% print result
fprintf('single query:  mAP = %f,	r1 precision = %f,	r5 precision = %f\r\n', mean(ap), CMC(1), CMC(5));
%% plot CMC curves
figure;
s = 50;
CMC_curve = CMC;
plot(1:s, CMC_curve(:, 1:s));

author: JDAI-CV

repo: VeRidataset

save_path: github-repos/MATLAB/JDAI-CV-VeRidataset


text:
function el = export_contour_image_sequence(args)
%Subfunction to handle contour_image sequences within the
%structure_set module.  Uses the same layout and principle as the parent
%function.
%
%   This function takes a scanInfo object.
%
%JRA 06/23/06
%NAV 07/19/16 updated to dcm4che3
%   replaced ml2dcm_Element to data2dcmElement
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Init output element to empty.
el = [];

%Unpack input data.
tag                 = args.tag;
%scanInfo            = args.data{1};
contourS            = args.data{1};
template            = args.template;

switch tag
    case 528720     %0008,1150  Referenced SOP Class UID
        data = contourS.SOP_Class_UID;
        el = data2dcmElement(data, tag);

    case 528725     %0008,1155  Referenced SOP Instance UID
        data = contourS.SOP_Instance_UID;
        el = data2dcmElement(data, tag);

    case 528736     %0008,1160  Referenced Frame Number
        %Currently not implemented.

    otherwise
        warning(['No methods exist to populate DICOM structure_set module''s contour_image_sequence field ' dec2hex(tag,8) '.']);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef Feature2 < handle

    properties
        feature
        descriptor
    end

    methods
        function f = Feature2(ff, dd)
            if nargin == 0
                f.feature = rand(100, 3);
                f.descriptor = rand(100, 20);
            elseif nargin == 1
                f.feature =ff;
                f.descriptor = [];
            elseif nargin == 2
                f.feature =ff;
                f.descriptor = dd;
            end
        end

        function s = char(f)
            s = '';
            for i=1:numrows(f.feature)
                s = char(s, sprintf('%f %f', f.feature(i,1:2)));
            end
        end

        function display(f)
            disp( f.char() );
        end

        function z = bob(f, a,b,c)
            fprintf('nargin=%d, a=%f\n', nargin, a);
            z = a;
        end

        function out = end(f,k,n)
            k,n
            %out = end(A.feature,k,n);
            if k == 1
                out = numrows(f.feature);
            else
                out = 1;
            end
        end

        function s = size(f)
            s = [numrows(f.feature) 1];
        end

        function n = length(f)
            n = numrows(f.feature);
        end
        %{
        function n = numel(f)
            n = numrows(f.feature);
        end
        %}
    end % methods
end % classdef

author: petercorke

repo: machinevision-toolbox-matlab


text:
function i4vec_sort_heap_d_test ( )

%*****************************************************************************80
%
%% I4VEC_SORT_HEAP_D_TEST tests I4VEC_SORT_HEAP_D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 20;
  b = 0;
  c = 3 * n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_SORT_HEAP_D_TEST\n' );
  fprintf ( 1, '  For a vector of integers,\n' );
  fprintf ( 1, '  I4VEC_SORT_HEAP_D descending sorts.\n' );

  seed = 123456789;

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  i4vec_print ( n, a, '  Unsorted:' );

  a = i4vec_sort_heap_d ( n, a );

  i4vec_print ( n, a, '  Descending sorted:' );

  return
end