text:
function terms = strsplit_re(s, delimiter)
%STRSPLIT_RE Splits a string into multiple terms with regex delimiter
%
%   terms = strsplit_re(s, delimiter);
%       splits the string s into multiple terms with delimiter specified
%       by a regular expression.
%
%       The extracted terms are returned in form of a cell array of
%       strings.
%
%   Remarks
%   -------
%       - If there are two consecutive non-whitespace delimiters, it is
%         regarded that there is an empty-string term between them.
%
%   Examples
%   --------
%       % split a string by whitespaces
%       ts = strsplit_re('I am using MATLAB', '\s+');
%       ts <- {'I', 'am', 'using', 'MATLAB'}
%
%       % extract terms delimited by commas
%       ts1 = strsplit_re('Apple, Orange , Banana', ',');
%       ts1 <- {'Apple', ' Orange ', ' Banana'}
%
%       ts2 = strsplit_re('Apple, Orange , Banana', '\s*,\s*');
%       ts2 <- {'Apple', 'Orange', 'Banana'}
%
%       % consecutive delimiters results in empty terms
%       ts = strsplit_re(',Apple,,Orange,Banana', ',');
%       ts <- {'', 'Apple', '', 'Orange', 'Banana'}
%

%   History
%   -------
%       - Created by Dahua Lin, on Oct 9, 2008
%

%% parse and verify input arguments

assert(ischar(s) && ndims(s) == 2 && size(s,1) <= 1, ...
    'strsplit_re:invalidarg', ...
    'The first input argument should be a char string.');

d = delimiter;
assert(ischar(d) && ndims(d) == 2 && size(d,1) == 1 && ~isempty(d), ...
    'strsplit_re:invalidarg', ...
    'The delimiter should be a non-empty char string.');

%% main

[sp, ep] = regexp(s, d, 'start', 'end');

if ~isempty(sp)
    nt = numel(sp) + 1;
    terms = cell(1, nt);
    p = 1;
    for i = 1 : nt-1
        if p < sp(i)
            terms{i} = s(p:sp(i)-1);
        else
            terms{i} = '';
        end
        p = ep(i) + 1;
    end
    if p <= length(s)
        terms{nt} = s(p:end);
    else
        terms{nt} = '';
    end
else
    terms = {s};
end

author: Sable

repo: mcbench-benchmarks


text:
function value = p23_exact ( dim_num )

%*****************************************************************************80
%
%% P23_EXACT returns the exact integral for problem 23.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Output, real VALUE, the exact value of the integral.
%
  c = 0.0;
  c = p23_r8 ( 'G', 'C', c );

  e = [];
  e = p23_i4vec ( 'G', 'E', dim_num, e );

  value = c;
  for i = 1: dim_num
    value = value * gamma ( e(i) + 1 );
  end

  value = value / gamma ( sum ( e(1:dim_num) ) + 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef OptimalExponentForMachingVigdergauz < handle

    properties (Access = private)
        sample
        rho
        txi
        rhoV
        txiV
        errorV
        qV
        mxV
        myV
    end

    methods (Access = public)

        function obj = OptimalExponentForMachingVigdergauz()
            obj.init()
            obj.compute()
            obj.save();
        end
    end

    methods (Access = private)

        function init(obj)
            obj.createSample();
        end

        function createSample(obj)
            s.type = 'FromMxMy';
            obj.sample = SamplePointsCreatorForOptimalExponentComputer.create(s);
            obj.sample.compute();
        end

        function compute(obj)
            for ipoint = 1: length(obj.sample.rhoV)
                obj.rho = obj.sample.rhoV(ipoint);
                obj.txi = obj.sample.txiV(ipoint);
                [qOpt,error] = obj.computeOptimalExponent();
                obj.qV(ipoint) = qOpt;
                obj.errorV(ipoint) = error;
                obj.mxV(ipoint) = obj.computeMx(obj.txi,obj.rho,qOpt);
                obj.myV(ipoint) = obj.computeMy(obj.txi,obj.rho,qOpt);
            end
        end

        function [qOpt,error] = computeOptimalExponent(obj)
            s.rho = obj.rho;
            s.txi = obj.txi;
            s.savingFrames = true;
            optimalExponent = OptimalExponentClosestToVigergauz(s);
            optimalExponent.compute();
            qOpt = optimalExponent.qOpt;
            error = optimalExponent.error;
        end

        function save(obj)
           d.mx = obj.mxV;
           d.my = obj.myV;
           d.rho = obj.rhoV;
           d.txi = obj.txiV;
           d.q  = obj.qV;
           d.error = obj.errorV;
           fN = 'OptimalSuperEllipseExponentMatchingVigdergauz';
           pD = 'Topology Optimization/Vademecums';
           file2SaveName = [pD,'/',fN,'.mat'];
           save(file2SaveName,'d');
        end

    end

    methods (Access = private, Static)

        function mx = computeMx(txi,rho,q)
            mx = SuperEllipseParamsRelator.mx(txi,rho,q);
        end

        function my = computeMy(txi,rho,q)
            my = SuperEllipseParamsRelator.my(txi,rho,q);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function img = discreteCurve(varargin)
%DISCRETECURVE Discretize a planar curve
%
%   IMG = discreteCurve(DIM, CURVE, WIDTH)
%   DIM is the size of image, with the format [x0 dx x1;y0 dy y1]
%   CURVE is a series of points describing the curve
%   WIDTH is the max distance between pixel centers and points of the
%   curve.
%
%   IMG = discreteCurve(LX, LY, ...);
%   Specifes the pixels coordinates with the two row vectors LX and LY.
%
%   Example
%   % creates a ring
%   circle = circleAsPolygon([25 25 15], 120);
%   img = discreteCurve([1 1 50;1 1 50], circle, 3);
%   imshow(img);
%
%   See Also
%   imShapes, discretePolyline, discretePolygon

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2007-03-19
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

%   HISTORY
%   19/06/2007: update doc
%   04/03/2009: use meshgrid
%   29/05/2009: use more possibilities for specifying grid

% compute coordinate of image voxels
[lx, ly, varargin] = parseGridArgs(varargin{:});
[x, y]   = meshgrid(lx, ly);

% get polyline vertex coordinates
curve = varargin{1};
varargin(1) = [];

% determines width of the polyline
width = 2;
if ~isempty(varargin)
    width = varargin{1};
end

try
    % try with vectorized version (greedy !!!)
    dist = reshape(minDistancePoints([x(:) y(:)], curve), size(x));
catch
    % if not enough memory, use loop instead
    dist = zeros(size(x));
    for i = 1:length(lx)
        dist(:, i) = minDistancePoints([x(:,i) y(:,i)], curve);
    end
end

% create image : simple threshold over 2 dimensions
img = abs(dist) < width;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [ n_data, n, x, fx ] = he_polynomial_values ( n_data )

%*****************************************************************************80
%
%% HE_POLYNOMIAL_VALUES: tabulated values of He(i,x).
%
%  Discussion:
%
%    He(i,x) represents the probabilist's Hermite polynomial.
%
%    In Mathematica, the function can be evaluated by:
%
%      He(n,x) = HermiteH[n,x/Sqrt[2]] / Sqrt [ 2^n ]
%
%  First terms:
%
%   1
%   X
%   X^2  -  1
%   X^3  -  3 X
%   X^4  -  6 X^2 +   3
%   X^5  - 10 X^3 +  15 X
%   X^6  - 15 X^4 +  45 X^2 -   15
%   X^7  - 21 X^5 + 105 X^3 -  105 X
%   X^8  - 28 X^6 + 210 X^4 -  420 X^2 +  105
%   X^9  - 36 X^7 + 378 X^5 - 1260 X^3 +  945 X
%   X^10 - 45 X^8 + 630 X^6 - 3150 X^4 + 4725 X^2 - 945
%
%  Recursion:
%
%    He(0,X) = 1,
%    He(1,X) = X,
%    He(N,X) = X * He(N-1,X) - (N-1) * He(N-2,X)
%
%  Norm:
%
%    Integral ( -oo < X < +oo ) exp ( - 0.5 * X^2 ) * He(M,X) He(N,X) dX
%    = sqrt ( 2 * pi ) * N! * delta ( N, M )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0
%    before the first call.  On each call, the routine increments N_DATA by 1,
%    and returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the order of the polynomial.
%
%    Output, real X, the point where the polynomial is evaluated.
%
%    Output, real FX, the value of the function.
%
  n_max = 18;

  fx_vec = [ ...
    1.000000000000000E+00, ...
    5.000000000000000E+00, ...
    24.00000000000000E+00, ...
    110.0000000000000E+00, ...
    478.0000000000000E+00, ...
    1950.000000000000E+00, ...
    7360.000000000000E+00, ...
    25100.00000000000E+00, ...
    73980.00000000000E+00, ...
    169100.0000000000E+00, ...
    179680.0000000000E+00, ...
   -792600.0000000000E+00, ...
   -5939480.000000000E+00, ...
    0.000000000000000E+00, ...
    6.281250000000000E+00, ...
    6.000000000000000E+00, ...
    18.00000000000000E+00, ...
    90150.00000000000E+00 ];

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

  x_vec = [ ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    5.0E+00, ...
    0.0E+00, ...
    0.5E+00, ...
    1.0E+00, ...
    3.0E+00, ...
    1.0E+01 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    x = 0.0;
    fx = 0.0;
  else
    n = n_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function y = base64decode(x)
%BASE64DECODE Perform base64 decoding on a string.
%
%   BASE64DECODE(STR) decodes the given base64 string STR.
%
%   Any character not part of the 65-character base64 subset set is silently
%   ignored.  Characters occuring after a '=' padding character are never
%   decoded.
%
%   STR doesn't have to be a string.  The only requirement is that it is a
%   vector containing values in the range 0-255.
%
%   If the length of the string to decode (after ignoring non-base64 chars) is
%   not a multiple of 4, then a warning is generated.
%
%   This function is used to decode strings from the Base64 encoding specified
%   in RFC 2045 - MIME (Multipurpose Internet Mail Extensions).  The Base64
%   encoding is designed to represent arbitrary sequences of octets in a form
%   that need not be humanly readable.  A 65-character subset ([A-Za-z0-9+/=])
%   of US-ASCII is used, enabling 6 bits to be represented per printable
%   character.
%
%   See also BASE64ENCODE.

%   Author:      Peter John Acklam
%   Time-stamp:  2004-09-20 08:20:50 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % check number of input arguments
   error(nargchk(1, 1, nargin));

   % remove non-base64 chars
   x = x (   ( 'A' <= x & x <= 'Z' ) ...
           | ( 'a' <= x & x <= 'z' ) ...
           | ( '0' <= x & x <= '9' ) ...
           | ( x == '+' ) | ( x == '=' ) | ( x == '/' ) );

   if rem(length(x), 4)
      warning('Length of base64 data not a multiple of 4; padding input.');
   end

   k = find(x == '=');
   if ~isempty(k)
      x = x(1:k(1)-1);
   end

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Now perform the following mapping
   %
   %   A-Z  ->  0  - 25
   %   a-z  ->  26 - 51
   %   0-9  ->  52 - 61
   %   +    ->  62
   %   /    ->  63

   y = repmat(uint8(0), size(x));
   i = 'A' <= x & x <= 'Z'; y(i) =    - 'A' + x(i);
   i = 'a' <= x & x <= 'z'; y(i) = 26 - 'a' + x(i);
   i = '0' <= x & x <= '9'; y(i) = 52 - '0' + x(i);
   i =            x == '+'; y(i) = 62 - '+' + x(i);
   i =            x == '/'; y(i) = 63 - '/' + x(i);
   x = y;

   nebytes = length(x);         % number of encoded bytes
   nchunks = ceil(nebytes/4);   % number of chunks/groups

   % add padding if necessary
   if rem(nebytes, 4)
      x(end+1 : 4*nchunks) = 0;
   end

   x = reshape(uint8(x), 4, nchunks);
   y = repmat(uint8(0), 3, nchunks);            % for the decoded data

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Rearrange every 4 bytes into 3 bytes
   %
   %    00aaaaaa 00bbbbbb 00cccccc 00dddddd
   %
   % to form
   %
   %    aaaaaabb bbbbcccc ccdddddd

   y(1,:) = bitshift(x(1,:), 2);                    % 6 highest bits of y(1,:)
   y(1,:) = bitor(y(1,:), bitshift(x(2,:), -4));    % 2 lowest bits of y(1,:)

   y(2,:) = bitshift(x(2,:), 4);                    % 4 highest bits of y(2,:)
   y(2,:) = bitor(y(2,:), bitshift(x(3,:), -2));    % 4 lowest bits of y(2,:)

   y(3,:) = bitshift(x(3,:), 6);                    % 2 highest bits of y(3,:)
   y(3,:) = bitor(y(3,:), x(4,:));                  % 6 lowest bits of y(3,:)

   % remove padding
   switch rem(nebytes, 4)
      case 2
         y = y(1:end-2);
      case 3
         y = y(1:end-1);
   end

   % reshape to a row vector and make it a character array
   y = char(reshape(y, 1, numel(y)));

author: jianxiongxiao

repo: ProfXkit


text:
% like mgram1, except we use a durational HMM instead of an HHMM2

past = 0;

words = {'the', 't', 'h', 'e'};
data = 'the';
nwords = length(words);
word_len = zeros(1, nwords);
word_prob = normalise(ones(1,nwords));
word_logprob = log(word_prob);
for wi=1:nwords
  word_len(wi)=length(words{wi});
end
D = max(word_len);

alphasize = 26*2;
data = letter2num(data);
T = length(data);

% node numbers
W = 1; % top level state = word id
L = 2; % bottom level state = letter position within word
F = 3;
O = 4;

ss = 4;
intra = zeros(ss,ss);
intra(W,[F L O])=1;
intra(L,[O F])=1;

inter = zeros(ss,ss);
inter(W,W)=1;
inter(L,L)=1;
inter(F,[W L O])=1;

% node sizes
ns = zeros(1,ss);
ns(W) = nwords;
ns(L) = D;
ns(F) = 2;
ns(O) = alphasize;
ns2 = [ns ns];

% Make the DBN
bnet = mk_dbn(intra, inter, ns, 'observed', O);
eclass = bnet.equiv_class;

% uniform start distrib over words, uniform trans mat
Wstart = normalise(ones(1,nwords));
Wtrans = mk_stochastic(ones(nwords,nwords));

% always start in state d = length(word) for each bottom level HMM
Lstart = zeros(nwords, D);
for i=1:nwords
  l = length(words{i});
  Lstart(i,l)=1;
end

% make downcounters
RLtrans = mk_rightleft_transmat(D, 0); % 0 self loop prob
Ltrans = repmat(RLtrans, [1 1 nwords]);

% Finish when downcoutner = 1
Fprob = zeros(nwords, D, 2);
Fprob(:,1,2)=1;
Fprob(:,2:end,1)=1;

% Define CPDs for slice
bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', Wstart);
bnet.CPD{eclass(L,1)} = tabular_CPD(bnet, L, 'CPT', Lstart);
bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', Fprob);

% Define CPDs for slice 2
bnet.CPD{eclass(W,2)} = hhmmQ_CPD(bnet, W+ss, 'Fbelow', F, 'startprob', Wstart,  'transprob', Wtrans);
bnet.CPD{eclass(L,2)} = hhmmQ_CPD(bnet, L+ss, 'Fself', F, 'Qps', W+ss, 'startprob', Lstart, 'transprob', Ltrans);

if 0
% To test it is generating correctly, we create an artificial
% observation process that capitalizes at the start of a new segment
% Oprob(Ft-1,Qt,Dt,Yt)
Oprob = zeros(2,nwords,D,alphasize);
Oprob(1,1,3,letter2num('t'),1)=1;
Oprob(1,1,2,letter2num('h'),1)=1;
Oprob(1,1,1,letter2num('e'),1)=1;
Oprob(2,1,3,letter2num('T'),1)=1;
Oprob(2,1,2,letter2num('H'),1)=1;
Oprob(2,1,1,letter2num('E'),1)=1;
Oprob(1,2,1,letter2num('a'),1)=1;
Oprob(2,2,1,letter2num('A'),1)=1;
Oprob(1,3,1,letter2num('b'),1)=1;
Oprob(2,3,1,letter2num('B'),1)=1;
Oprob(1,4,1,letter2num('c'),1)=1;
Oprob(2,4,1,letter2num('C'),1)=1;

% Oprob1(Qt,Dt,Yt)
Oprob1 = zeros(nwords,D,alphasize);
Oprob1(1,3,letter2num('t'),1)=1;
Oprob1(1,2,letter2num('h'),1)=1;
Oprob1(1,1,letter2num('e'),1)=1;
Oprob1(2,1,letter2num('a'),1)=1;
Oprob1(3,1,letter2num('b'),1)=1;
Oprob1(4,1,letter2num('c'),1)=1;

bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', Oprob);
bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', Oprob1);

evidence = cell(ss,T);
%evidence{W,1}=1;
sample = cell2num(sample_dbn(bnet, 'length', T, 'evidence', evidence));
str = num2letter(sample(4,:))
end

[log_obslik, obslik, match] = mk_mgram_obslik(lower(data), words, word_len, word_prob);
% obslik(j,t,d)
softCPDpot = cell(ss,T);
ens = ns;
ens(O)=1;
ens2 = [ens ens];
for t=2:T
  dom = [F W+ss L+ss O+ss];
  % tab(Ft-1, Q2, Dt)
  tab = ones(2, nwords, D);
  if past
    tab(1,:,:)=1; % if haven't finished previous word, likelihood is 1
    tab(2,:,:) = squeeze(obslik(:,t,:)); % otherwise likelihood of this segment
  else
    for d=1:max(1,min(D,T+1-t))
      tab(2,:,d) = squeeze(obslik(:,t+d-1,d));
    end
  end
  softCPDpot{O,t} = dpot(dom, ens2(dom), tab);
end
t = 1;
dom = [W L O];
% tab(Q2, Dt)
tab = ones(nwords, D);
if past
  tab = squeeze(obslik(:,t,:));
else
  for d=1:min(D,T-t)
    tab(:,d) = squeeze(obslik(:,t+d-1,d));
  end
end
softCPDpot{O,t} = dpot(dom, ens(dom), tab);

%bnet.observed = [];
% uniformative observations
%bnet.CPD{eclass(O,2)} = tabular_CPD(bnet, O+ss, 'CPT', mk_stochastic(ones(2,nwords,D,alphasize)));
%bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', mk_stochastic(ones(nwords,D,alphasize)));

engine = jtree_dbn_inf_engine(bnet);
evidence = cell(ss,T);
% we add dummy data to O to force its effective size to be 1.
% The actual values have already been incorporated into softCPDpot
evidence(O,:) = num2cell(ones(1,T));
[engine, ll_dbn] = enter_evidence(engine, evidence, 'softCPDpot', softCPDpot);

%evidence(F,:) = num2cell(2*ones(1,T));
%[engine, ll_dbn] = enter_evidence(engine, evidence);

gamma = zeros(nwords, T);
for t=1:T
  m = marginal_nodes(engine, [W F], t);
  gamma(:,t) = m.T(:,2);
end

gamma

xidbn = zeros(nwords, nwords);
for t=1:T-1
  m = marginal_nodes(engine, [W F W+ss], t);
  xidbn = xidbn + squeeze(m.T(:,2,:));
end

% thee
% xidbn(1,4)  = 0.9412  the->e
% (2,3)=0.0588 t->h
% (3,4)=0.0588 h-e
% (4,4)=0.0588 e-e

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [irf_record]=irf(beta_gibbs,It,Bu,IRFperiods,n,m,p,k)

% function [irf_record]=irf(beta_gibbs,It,Bu,IRFperiods,n,m,p,k)
% runs the gibbs sampler to obtain draws from the posterior distribution of IRFs
% inputs:  - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - integer 'It': total number of iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'Bu': number of burn-in iterations of the Gibbs sampler (defined p 28 of technical guide)
%          - integer 'IRFperiods': number of periods for IRFs
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
% outputs: -  cell 'irf_record': record of the gibbs sampler draws for the IRFs

% this function implements algorithm 2.2.1

% create the cell aray that will store the values from the simulations
irf_record=cell(n,n);

Bgibbs=reshape(beta_gibbs,k,n,It-Bu);

% deal with shocks in turn
for ii=1:n

   % step 1: repeat the simulation process a number of times equal to the number of Gibbs iterations
   for kk=1:It-Bu

   % step 3: draw beta from its posterior distribution
   B=squeeze(Bgibbs(:,:,kk));
   % create a matrix of zeros of dimension p*n
   Y=zeros(p,n);
   %  step 2: set the value of the last row, column i, equal to 1
   Y(p,ii)=1;

% if prior~=61
      % step 4: for each iteration kk, repeat the algorithm for periods T+1 to T+h
      for jj=1:IRFperiods-1
      % use the function lagx to obtain a matrix temp, containing the endogenous regressors
      temp=bear.lagx(Y,p-1);
      % define the vector X
      X=[temp(end,:) zeros(1,m)];
      % obtain the predicted value for T+jj
      yp=X*B;
      % concatenate yp at the top of Y
      Y=[Y;yp];
      % repeat until values are obtained for T+h
      end
% elseif prior==61 %for prior=61 mean adjusted model
% 	  % for each iteration kk, repeat the algorithm for periods T+1 to T+h
%       for jj=1:IRFperiods
%       % use the function lagx to obtain the matrix X
%       X=bear.lagx(Y,p-1);
%       X=X(end,:);
%       % obtain predicted value for T+jj
%       yp=X*B;
%       % concatenate yp at the top of Y
%       Y=[Y;yp];
%       % repeat until values are obtained for T+h
%       end
% end

   % step 5: record the results from current iteration in cell irf_record
      % loop over variables
      for jj=1:n
      % consider column jj of matrix 'Y' and trim the (p-1) initial periods: what remains is the series of IRFs for period T to period T+h-1, for variable jj
      temp=Y(p:end,jj);
      % record these values in the corresponding matrix of irf_record
      irf_record{jj,ii}(kk,:)=temp';
      end

   % then go for next iteration
   end

% conduct the same process with shocks in other variables
end

author: european-central-bank

repo: BEAR-toolbox


text:
%% STEP 0: Calculate Distortion Parameters
% This is a main script to perform the following steps:
% 1) Estimate camera parameters from checkerboard images
% 2) Use these parameters to correct the distortion
% 3) Plot the camera parameters before and after the correction
% 4) Save the parameters to correct the distortion from image points in Step 1p and Step 3.

%%
clear all; close all

fs=get(0, 'DefaultUIControlFontSize');
set(0, 'DefaultUIControlFontSize', 10);

%% NEW or REPEAT?
% select repeat if you already ran analysis with these same images, but want to repeat with a different distortion model
% repeat is faster because the point detection is skipped (which is the slowest part in this code)

repeatButton = questdlg('New analysis of repeat calibration with a different model?', 'New analysis of repeat calibration with a different model?', 'New', 'Repeat', 'New');
switch repeatButton
    case 'Repeat'
        repeatLogic=true(1);
    case 'New'
        repeatLogic=false(1);
end

%% CHOOSE PATHS OPTIONS

if repeatLogic % if a repeated analysis, request to point to the existing parameters
    initialPath=pwd;
    folder_paths = uipickfiles('FilterSpec',initialPath,'Prompt','Select one or multiple cameraCBparameters files');
else
    % initial image path. Change this path if you want the UI to start with a specific path. Otherwise leave [] or pwd.
    imagePath=pwd;
    % select the folder containing the checkerboard images (if imagePath=[] then the initial path is the current path)
    folder_paths = uipickfiles('FilterSpec',imagePath,'Prompt','Select one or multiple folders containing checkerboard images for analysis');
end

% save camera parameters? y/n. if y, choose save path and overwrite options
[saveCameraParametersLogic,savePath]=QsaveCameraParameters(folder_paths);

% save undistorted images? y/n. if yes, warn for overwriting. undistorted images will be saved in savePath
[saveUndistortedImagesLogic,overWriteUDimagesLogic]=QsaveUndistortedImages;

% figures path to save the plotted results
figuresPath=fullfile(savePath, 'figures');
warning('off','MATLAB:MKDIR:DirectoryExists');
mkdir(figuresPath);

%% SELECT CHECKERBOARD PARAMETERS

if repeatLogic % if repeated analysis, extract the checkerboard parameters from the files
    cameraCBparameters1=load(folder_paths{1});
    Nrows=cameraCBparameters1.cameraCBparameters.boardSize(1);
    Ncols=cameraCBparameters1.cameraCBparameters.boardSize(2);
    squareSize=cameraCBparameters1.cameraCBparameters.squareSize;
else % new analysis. ask for checkerboard parameters
    %initial parameters: these are the default parameters which appear in the dialog box. Change them if you want other numbers to appear instead
    Nrows=15; % Number of black rows (should be uneven)
    Ncols=20; % Number of black columns (should be even)
    squareSize=10; %[mm]
    % dialog box
    answer = inputdlg({'Enter number of rows (odd):','Enter number of columns (even):','Enter square size [mm]:'},'Input',[1,50],{num2str(Nrows),num2str(Ncols),num2str(squareSize)});
    % extract answers
    Nrows=str2double(answer{1});
    Ncols=str2double(answer{2});
    squareSize=str2double(answer{3});
end

%% SELECT DISTORTION MODEL

% dialog box for selecting the distortion model.
% The default is the full model (3 radial parameters, 2 tangential (1 stands for true), and a skew parameter (1 = true).
answer = inputdlg({'Enter number of radial distortion coefficients (2 or 3):',...
    'Estimate tangential distortion? (1 or 0 for yes/no):',...
    'Estimate skew? (1 or 0 for yes/no)'},...
    'Input',[1,70],{'3','1','1'});
optStruct=struct;
optStruct.NumRadialDistortionCoefficients=str2num(answer{1});
optStruct.EstimateTangentialDistortion=logical(str2num(answer{2}));
optStruct.EstimateSkew=logical(str2num(answer{3}));

%% compute and plot camera parameters for each camera

Ncam=numel(folder_paths); % number of cameras in this analysis
cameraCBparametersAllCams=cell(Ncam,1); % assign cell array for all camera parmaters

for ic=1:Ncam % loop over all cameras

    % Extract images and info
    if repeatLogic % if repeated analysis, extract cameraCBparameters and image info from files
        cameraCBparameters=load(folder_paths{ic});
        cameraCBparameters=cameraCBparameters.cameraCBparameters;
        CBimagesInfo=cameraCBparameters.imagesInfo;
        I1=imread(CBimagesInfo.imageFileNames{1}); % read first image from path
        I=zeros([size(I1) CBimagesInfo.Nimages],'uint8'); % reallocate image var
        for ip=1:CBimagesInfo.Nimages
            I(:,:,:,ip)=imread(CBimagesInfo.imageFileNames{ip}); % load all images into I
        end
        CBimagesInfo.I=I;
        % plot all images in one figure
        plotAllCameraImages(CBimagesInfo);
        % re-calculate the distortion parameters with the selected model
        hm=msgbox(['Please wait while computing distortion parameters for camera ' num2str(CBimagesInfo.icam)]);
        cameraCBparameters=RecalculateCBcalibrationParameters(cameraCBparameters,optStruct);
        delete(hm);
    else
        % if New, extract only image info
        CBimagesInfo=extractImagesInfo(folder_paths{ic});
        % plot all images in one figure
        plotAllCameraImages(CBimagesInfo);
        % Extract images, Detect the checkerboard points, calculate camera paramters, and save a structure containing all necessary parameters
        set(0, 'DefaultUIControlFontSize', 11);
        hm=msgbox(['Please wait while computing distortion parameters for camera ' num2str(CBimagesInfo.icam)]);
        cameraCBparameters=calculateCBcalibrationParameters(CBimagesInfo,squareSize,optStruct);
        delete(hm);
            % check if detected boardsize matches entered values
            if (cameraCBparameters.boardSize(1)~=Nrows) || (cameraCBparameters.boardSize(2)~=Ncols)
                error('Detected number of columns or rows does not match entered values');
            end
    end

    % plot camera parameters and reorojection errors before and after distortion correction
    plot_camera_parameters_2tabs(cameraCBparameters);
        if saveCameraParametersLogic
            savefig(fullfile(figuresPath,[ 'params_cam' num2str(CBimagesInfo.icam)]));
        end

    % plot reprojected points vs. true points and straight lines on each image
    plot_reprojectVSreal_points(CBimagesInfo,cameraCBparameters);

    % undistort images and save if required
    hm=msgbox(['Please wait while correcting distortion from images of ' num2str(CBimagesInfo.icam)]);
    undistortImagesSavePlot(CBimagesInfo,cameraCBparameters,saveUndistortedImagesLogic,overWriteUDimagesLogic,savePath)
    delete(hm);

    % remove images with high errors
    answer=questdlg('Do you want to remove some of the images with higher errors?','Remove images?','Yes','No','No');
    switch answer
        case 'Yes'
            % reprojection errors per image
            Error=cameraCBparameters.cameraParameters.ReprojectionErrors;
            ErrorMgn=squeeze(sqrt(Error(:,1,:).^2+Error(:,2,:).^2));
            meanError=squeeze(mean(ErrorMgn,1));
            prctile80Error=prctile(meanError,80);
            % select threshold
            answer=inputdlg('Select threshold (mean error in pixels) for removing images','select threshold',1,{num2str(prctile80Error)});
            threshold=str2num(answer{1});
            indImgHighError=cameraCBparameters.imagesUsed(find(meanError>threshold));
            if length(indImgHighError)>0
                msgbox([ num2str(length(indImgHighError)) ' images will be removed from the calibration parameters calculation']);
                fileNamesImgHighError=cameraCBparameters.imagesInfo.imageFileNames(indImgHighError)';
                % delete high error images
                CBimagesInfo.I(:,:,:,indImgHighError)=[];
                CBimagesInfo.imageFileNames(indImgHighError)=[];
                CBimagesInfo.Nimages=numel(CBimagesInfo.imageFileNames);
                % plot all images in one figure
                plotAllCameraImages(CBimagesInfo);
                % Extract images, Detect the checkerboard points, calculate camera paramters, and save a structure containing all necessary parameters
                hm=msgbox(['Please wait while re-computing distortion parameters for camera ' num2str(CBimagesInfo.icam)]);
                cameraCBparameters=calculateCBcalibrationParameters(CBimagesInfo,squareSize,optStruct);
                delete(hm);
                % plot camera parameters and reorojection errors before and after distortion correction
                plot_camera_parameters_2tabs(cameraCBparameters);
                if saveCameraParametersLogic
                    savefig(fullfile(figuresPath,[ 'params_cam', num2str(CBimagesInfo.icam)]));
                end
            end
    end

    % save camera parameters into the cell array of all cameras
    cameraCBparametersAllCams{ic}=cameraCBparameters;
    % save parameters into savePath
    if saveCameraParametersLogic
        save(fullfile(savePath, ['cameraCBparameters_cam_', num2str(cameraCBparameters.icam)]),'cameraCBparameters');
    end

end

% save cell array containing the camera parameters for all cameras in this analysis
if saveCameraParametersLogic
    save(fullfile(savePath, 'cameraCBparametersAllCams'),'cameraCBparametersAllCams');
end

%% plot camera instrinsic statistics for all cameras (if more than 1 camera)

Ncam=numel(cameraCBparametersAllCams);
if Ncam>1
    plotButton = questdlg('Plot intrinsic parameters statistics for all cameras?', 'Plot?', 'Yes', 'No', 'Yes');
    if strcmp(plotButton,'Yes')
        plotIntrinsicStatsAll(cameraCBparametersAllCams);
                if saveCameraParametersLogic
                    savefig(fullfile(figuresPath, 'IntrinsicStats'));
                end
    end
end

%% finish

h=msgbox('STEP0 is completed');
h.CurrentAxes.Children.FontSize=11;

set(0, 'DefaultUIControlFontSize', fs);

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function res = chantype(this, varargin)
% Method for setting/getting channel types
% FORMAT chantype(this, ind, type)
%   ind - channel index
%   type - type (string: 'EEG', 'MEG', 'LFP' etc.)
%
% FORMAT chantype(this, ind), chantype(this)
% Sets channel types to default using Fieldtrip channelselection
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: chantype.m 5933 2014-03-28 13:22:28Z vladimir $

if this.montage.Mind==0
    res = getset(this, 'channels', 'type', varargin{:});
else
    if nargin == 3
        this.montage.M(this.montage.Mind) = getset(this.montage.M(this.montage.Mind), 'channels', 'type', varargin{:});
        res = this;
    else
        res = getset(this.montage.M(this.montage.Mind), 'channels', 'type', varargin{:});
    end
end

author: spm


text:
function m = cell2mat(c)
%CELL2MAT Convert the contents of a cell array into a single matrix.
%   M = CELL2MAT(C) converts a multidimensional cell array with contents of
%   the same data type into a single matrix. The contents of C must be able
%   to concatenate into a hyperrectangle. Moreover, for each pair of
%   neighboring cells, the dimensions of the cell's contents must match,
%   excluding the dimension in which the cells are neighbors. This constraint
%   must hold true for neighboring cells along all of the cell array's
%   dimensions.
%
%   The dimensionality of M, i.e. the number of dimensions of M, will match
%   the highest dimensionality contained in the cell array.
%
%   CELL2MAT is not supported for cell arrays containing cell arrays or
%   objects.
%
%	Example:
%	   C = {[1] [2 3 4]; [5; 9] [6 7 8; 10 11 12]};
%	   M = cell2mat(C)
%
%	See also MAT2CELL, NUM2CELL

% Copyright 1984-2006 The MathWorks, Inc.
% $Revision: 1.10.4.7 $  $Date: 2006/11/11 22:43:58 $

% Error out if there is no input argument
if nargin==0
    error('MATLAB:cell2mat:NoInputs',['No input argument specified. ' ...
        'There should be exactly one input argument.'])
end
% short circuit for simplest case
elements = numel(c);
if elements == 0
    m = [];
    return
end
if elements == 1
    if isnumeric(c{1}) || ischar(c{1}) || islogical(c{1}) || isstruct(c{1})
        m = c{1};
        return
    end
end
% Error out if cell array contains mixed data types
cellclass = class(c{1});
ciscellclass = cellfun('isclass',c,cellclass);
if ~all(ciscellclass(:))
    error('MATLAB:cell2mat:MixedDataTypes', ...
        'All contents of the input cell array must be of the same data type.');
end

% Error out if cell array contains any cell arrays or objects
ciscell = iscell(c{1});
cisobj = isobject(c{1});
if cisobj || ciscell
    error('MATLAB:cell2mat:UnsupportedCellContent',['Cannot support cell ' ...
        'arrays containing cell arrays or objects.']);
end

% If cell array of structures, make sure the field names are all the same
if isstruct(c{1})
    cfields = cell(elements,1);
    for n=1:elements
        cfields{n} = fieldnames(c{n});
    end
    % Perform the actual field name equality test
    if ~isequal(cfields{:})
        error('MATLAB:cell2mat:InconsistentFieldNames', ...
            ['The field names of each cell array element must be consistent ' ...
            'and in consistent order.'])
    end
end

% If cell array is 2-D, execute 2-D code for speed efficiency
if ndims(c) == 2
    rows = size(c,1);
    cols = size(c,2);
    if (rows < cols)
        m = cell(rows,1);
        % Concatenate one dim first
        for n=1:rows
            m{n} = cat(2,c{n,:});
        end
        % Now concatenate the single column of cells into a matrix
        m = cat(1,m{:});
    else
        m = cell(1, cols);
        % Concatenate one dim first
        for n=1:cols
            m{n} = cat(1,c{:,n});
        end
        % Now concatenate the single column of cells into a matrix
        m = cat(2,m{:});
    end
    return
end

csize = size(c);
% Treat 3+ dimension arrays

% Construct the matrix by concatenating each dimension of the cell array into
%   a temporary cell array, CT
% The exterior loop iterates one time less than the number of dimensions,
%   and the final dimension (dimension 1) concatenation occurs after the loops

% Loop through the cell array dimensions in reverse order to perform the
%   sequential concatenations
for cdim=(length(csize)-1):-1:1
    % Pre-calculated outside the next loop for efficiency
    ct = cell([csize(1:cdim) 1]);
    cts = size(ct);
    ctsl = length(cts);
    mref = {};

    % Concatenate the dimension, (CDIM+1), at each element in the temporary cell
    %   array, CT
    for mind=1:prod(cts)
        [mref{1:ctsl}] = ind2sub(cts,mind);
        % Treat a size [N 1] array as size [N], since this is how the indices
        %   are found to calculate CT
        if ctsl==2 && cts(2)==1
            mref = {mref{1}};
        end
        % Perform the concatenation along the (CDIM+1) dimension
        ct{mref{:}} = cat(cdim+1,c{mref{:},:});
    end
    % Replace M with the new temporarily concatenated cell array, CT
    c = ct;
end

% Finally, concatenate the final rows of cells into a matrix
m = cat(1,c{:});

repo: REST


text:
 function val = byte_to_uint32(buff)
%function val = byte_to_uint32(buff)
% convert 4-byte buffer into uint32

try
	val = typecast(buff, 'uint32')
catch
	u = @(x) uint32(x);
	buff = u(buff);
	s = u(256);
	val = buff(4) + s * (buff(3) + s * (buff(2) + s * buff(1)));
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [Pfa_asv, Pmiss_asv, Pmiss_spoof_asv] = obtain_asv_error_rates(tar_asv, non_asv, spoof_asv, thresh_asv)

Ntar_asv    = length(tar_asv);
Nnon_asv    = length(non_asv);
Nspoof_asv  = length(spoof_asv);

% Obtain ASV false alarm and miss rates
Pfa_asv     = sum(non_asv >= thresh_asv)./Nnon_asv;
Pmiss_asv   = sum(tar_asv <  thresh_asv)./Ntar_asv;
if isempty(spoof_asv)
    Pmiss_spoof_asv = [];
else
    Pmiss_spoof_asv = sum(spoof_asv <  thresh_asv)./Nspoof_asv;
end

author: nesl


text:
function thresholdOpt = learnThetaMS(params,scale)

fprintf('Learning theta_MS for scale = %d \n',scale);
bestScoreThreshold = -inf;
level = find(params.MS.scale == scale);

for idxThr = 1:length(params.MS.domain(level,:)) %loop over the possible threshold values

    threshold = params.MS.domain(level,idxThr);
    struct = load([params.trainingImages 'structGT.mat']);
    structGT= struct.structGT;%training GT

    scoreThreshold = 0;

    for idxImgGT = 1:length(structGT) %for every img in GT

        img = imread([params.trainingImages '/' structGT(idxImgGT).img]);

        saliencyMAP = saliencyMap(img,params.MS.filtersize,scale);%compute the saliency map - for the current scale
        thrmap = saliencyMAP >= threshold;
        salmap = saliencyMAP .* thrmap;
        thrmapIntegralImage = computeIntegralImage(thrmap);
        salmapIntegralImage =  computeIntegralImage(salmap);
        scoreScale = slidingWindowComputeScore(double(saliencyMAP), scale, 1, 1, threshold, salmapIntegralImage, thrmapIntegralImage);

        [xmin ymin xmax ymax score] = nms4d(double(scoreScale),scale,scale,params.MS.sizeNeighborhood);%non maximum supression

        indexPositive = find(score > 0);
        xmin = xmin(indexPositive);
        ymin = ymin(indexPositive);
        xmax = xmax(indexPositive);
        ymax = ymax(indexPositive);
        score = score(indexPositive);

        for idxObject = 1:size(structGT(idxImgGT).boxes,1) %for every annotated object

            gtBox = structGT(idxImgGT).boxes(idxObject,:);
            [height width ~] = size(img);
            gtBox(1) = gtBox(1)*scale/width;
            gtBox(3) = gtBox(3)*scale/width;
            gtBox(2) = gtBox(2)*scale/height;
            gtBox(4) = gtBox(4)*scale/height;

            maxPascalScore = 0;
            for w = 1:length(score)
                nmsBox = [xmin(w) ymin(w) xmax(w) ymax(w)];
                pascalScore = computePascalScore(gtBox,nmsBox);
                if maxPascalScore < pascalScore
                    maxPascalScore = pascalScore;
                end
            end
            scoreThreshold = scoreThreshold + maxPascalScore;
        end

    end

    if bestScoreThreshold < scoreThreshold
        bestScoreThreshold = scoreThreshold;
        thresholdOpt = threshold;
    end
    fprintf('Best current theta_MS for scale = %d is %f \n',scale,thresholdOpt)
end

end

function saliencyMAP = saliencyMap(inImg,filtersize,scale)

inImg = im2double(rgb2gray(inImg));
inImg = imresize(inImg,[scale,scale],'bilinear');

%Spectral Residual
myFFT = fft2(inImg);
myLogAmplitude = log(abs(myFFT));
myPhase = angle(myFFT);
mySmooth = imfilter(myLogAmplitude,fspecial('average',filtersize),'replicate');
mySpectralResidual = myLogAmplitude-mySmooth;
saliencyMAP = abs(ifft2(exp(mySpectralResidual+1i*myPhase))).^2;

%After Effect
saliencyMAP = imfilter(saliencyMAP,fspecial('disk',filtersize));
saliencyMAP = mat2gray(saliencyMAP);
end

author: Cloud-CV

repo: object-proposals


text:
function jac = p36_jac ( neqn, t, y )

%*****************************************************************************80
%
%% P36_JAC evaluates the jacobian for problem p36.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the jacobian.
%
%    Output, real JAC(NEQN,NEQN), the jacobian matrix.
%
  jac = zeros ( neqn, neqn );

  jac(1,1) =  0.0;
  jac(1,2) =  1.0;

  jac(2,1) = 1.0 - 3.0 * y(1).^2;
  jac(2,2) = 0.0;

  return
end


text:
function [fiberspool, fiberDiameter] =dtiCullFibers(fg, dt6filename, Tt, distanceCrit, fiberLenCrit, averLACrit)
% ER's idea of how to 'cull' fibers - Use with caution ...
%
% [fiberspool, fiberDiameter] = dtiCullFibers(fg, dt6filename, Tt, distanceCrit, fiberLenCrit, averLACrit)
%
% fg: Fiber group
% dt6filename:
% Tt:
% distanceCrit:  Separation between fibers
% fiberLenCrit: (mm)
% averLACrit:  (mm)
%
% Fibers are not removed (culled) from a fiber group if
%
% (1) Length of a trajectory > 10.0mm
% (2) Average linear anisotropy along fiber trajectory > 0.1
% (3) Fibers meeting these criteria have a between fiber distance >1 mm (on average).
%     Default threshold value is 1 mm;
%
% Although arguable, it is natural to set Tt to resampled data voxel size
% (as Zhang 2003 did), and . Zhang et al. (2003) IEEE Transactions on
% Vizualization and Computer Graphics 9(4)
%
%11/2008: ER wrote it
%
% (c) Stanford VISTA Team

totalNfibers=size(fg.fibers, 1);

if (~exist('fiberLenCrit','var') || isempty(fiberLenCrit))
    fiberLenCrit=10; %in mm
end

if(~exist('averLACrit','var') || isempty(averLACrit))
    averLACrit=0.1;%arbitrary
end

if(~exist('distanceCrit','var') || isempty(distanceCrit))
    distanceCrit = 1;% %distanceCrit=1.7*2;%for voxels of 2; 1.7*voxel Or do I mean voxels as in opts.stepSizeMm?
end

% What is Tt?
if(~exist('Tt','var') || isempty(Tt))
    Tt=1; %Tt=2;
end

%Deal with files with no subgroup field
if  ~isfield(fg, 'subgroup') || isempty(fg.subgroup)
    fg.subgroup=repmat(1, [1 size(fg.fibers, 1)]);
end

% (1) Fiber length: infer from number of points; length(fg.fibers{ii}). On
%average the distance between the neighbouring points is....opts.stepSizeMm= 1;

%Create a FG which consists of only fibers that are longer than fiberLenCrit
fg1=fg;
fiberLen = cellfun('length', fg.fibers);
fg1.fibers=fg.fibers(fiberLen>fiberLenCrit);
if ~isempty(fg.seeds)
    fg1.seeds=fg.seeds(fiberLen>fiberLenCrit, :);
end
fg1.subgroup=fg.subgroup(fiberLen>fiberLenCrit);

fg1.name=[fg.name '_>10mm'];
shortNfibers=size(fg.fibers, 1)-size(fg1.fibers, 1);
clear fg;
display(['Eliminated ' num2str(shortNfibers) ' fibers that were shorter than ' num2str(fiberLenCrit) ' mm']);

%(2) Linear anisotropy is defined in Westin, Peled, Gubjartsson, Kikinis,
%Jolesc (ISMRM, 1997): cl=(l1-l2)/(l1+l2+l3)
valName ='shapes'; interpMethod='trilin'; %(triplet of values indicating linearity, planarity and spherisity)
dt      = dtiLoadDt6(dt6filename);

%Create a FG which consists only of long fibers that have average LA of at
%least averLACrit
fg2 = fg1; fg2.seeds=[]; fg2.fibers=[]; fg2.subgroup=[];
fg2.name = [fg1.name '_la>' num2str(averLACrit)];

%Need to break the FG into pieces of 5000 so that dtiGetValFromFibers does not
%choke on it.
stepsize=5000;

fgtemp.seeds=[];

fprintf(1, ['Computing linear anisotropy for fibers (out of ' num2str(size(fg1.fibers, 1)) '): ']);
for group_of_fibers=1:stepsize:size(fg1.fibers, 1)

    fprintf(1, [num2str(group_of_fibers) ' to ' num2str(min(group_of_fibers+stepsize-1, size(fg1.fibers, 1))) ' ... ']);
    if ~isempty(fg1.seeds)
        fgtemp.seeds=fg1.seeds(group_of_fibers:min(group_of_fibers+stepsize-1, size(fg1.seeds, 1)), :);
    end
    fgtemp.fibers=fg1.fibers(group_of_fibers:min(group_of_fibers+stepsize-1, size(fg1.fibers, 1)));
    fgtemp.subgroup=fg1.subgroup(group_of_fibers:min(group_of_fibers+stepsize-1, size(fg1.fibers, 1)));

    val = dtiGetValFromFibers(dt.dt6, fgtemp, inv(dt.xformToAcpc), valName, interpMethod);
    shapes=cellfun(@mean, val, 'UniformOutput', 0);
    shapesC=vertcat(shapes{:});
    %shapesC(:, 1) %gives you linear anisotropy

    fg2.fibers=vertcat(fg2.fibers, fgtemp.fibers(shapesC(:, 1)>averLACrit));
    if ~isempty(fgtemp.seeds)
        fg2.seeds=vertcat(fg2.seeds, fgtemp.seeds(shapesC(:, 1)>averLACrit, :));
    end
    fg2.subgroup=horzcat(fg2.subgroup, fgtemp.subgroup(shapesC(:, 1)>averLACrit));

end
fprintf(1, '\n');
lowanisotropyNfibers=size(fg1.fibers, 1)-size(fg2.fibers, 1);
fprintf('Eliminated %d fibers that had linear anisotropy smaller than %.3f\n',...
    lowanisotropyNfibers, averLACrit);
clear fg1 fgtemp;

%Now, this approach has to be tested for robustness towards the order in
%which the fibers are tested. If a fiber B is considered too close to fiber
%A, fiber B is dropped (there is no creating an average representation of
%the two, as the latter might lead to weird drifft effects swipping
%streamsheets

% Shuffle the fibers to ensure the lack of effects that comes from
% systematic seed placement.
fg3 = dtiShuffleFibers(fg2);
clear fg2;

% This will be the selected fibers, I guess. (BW).  This method of creating
% the pool is odd and special cased.  It should be a general function.
fiberspool = fg3;
fiberspool.name = sprintf('%s-culled',fg3.name);
fiberspool.fibers = [];
fiberspool.seeds  = [];
fiberspool.subgroup = [];

% Always put one fiber (the first one) into the pool
fiberspool.fibers{1} = fg3.fibers{1};
if ~isempty(fg3.seeds)
    fiberspool.seeds = fg3.seeds(1, :);
end
% I guess for subgropu, too
fiberspool.subgroup = fg3.subgroup(1);

% Should be a mrvWaitbar?  It steps in 5000 fiber chunks.  Not sure why.
fprintf('Distance testing: %d fibers\n', size(fg3.fibers, 1));
stepSize = 5000;
for i=1:stepSize:size(fg3.fibers, 1) % Step over 5000 fibers at a time.

    % There is an index 'i', an index 'ii', and an index 'jj'. We can be
    % clearer. We compare the fibers in fg3 with each fibers in the pool We
    % add fibers from fg3 to the pool if they are not close to anything
    % already in the pool.
    for ii=i:min(i+(stepSize-1), size(fg3.fibers, 1))
        inpool = 1;   %Assume ii fiber in fg3 is in the keep pool
        for jj = 1:size(fiberspool.fibers, 1)

            % Compare the iith fiber in fg3 to each of the fibers in the
            % pool
            %
            % InterfiberZhangDistance should be renamed fgDistance.  That
            % function should have a switch to define which algorithm is
            % used, e.g. Zhang. The Dt level should be explained better.
            % Also, I think this algorithm is supposed to be using Dt + Tt,
            % not Dt.
            [Dt avgDt] = InterfiberZhangDistance(fg3.fibers{ii}, fiberspool.fibers{jj}, Tt);

            if Dt <= distanceCrit
                inpool=0; % Reject this fiber
                break
            end

        end

        % Add the new fiber to the pool.  It will be compared to the
        % remaining fibers in fg3.
        if inpool==1
            fiberspool.fibers=vertcat(fiberspool.fibers, fg3.fibers{ii});
            if ~isempty(fg3.seeds)
                fiberspool.seeds=vertcat(fiberspool.seeds, fg3.seeds(ii, :));
            end
            fiberspool.subgroup=horzcat(fiberspool.subgroup, fg3.subgroup(ii));
        end

    end
    fprintf('Distance testing: %d fibers of %d\n', ii, size(fg3.fibers, 1))
end

fprintf('Out of total %d fibers %d fibers remain\n',...
    totalNfibers,size(fiberspool.fibers, 1));
fprintf('Eliminated %d fibers because they were shorter than %d mm\n',...
    shortNfibers,fiberLenCrit);
fprintf('Eliminated %d fibers because linear anisotropy less than %d\n',...
    lowanisotropyNfibers, averLACrit);
fprintf('Eliminated %d fibers because they were adjacent to others\n',...
    size(fg3.fibers,1) - size(fiberspool.fibers, 1)) ;

% Add a comment
fg.cullingParams={'Tt', Tt, 'aboveThresDistanceCrit', distanceCrit};
fiberDiameter = Tt + distanceCrit;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function res = norm( x, safe )
    %NORM Norm of a TT/MPS tensor.
    %   norm(X) computes the Frobenius norm of the TT/MPS tensor X.
    %
    %   norm(X, SAFE) with SAFE=true computes the Frobenius norm of the TT/MPS tensor X
    %       with reorthogonalization to increase the accuracy
    %
    %   See also INNERPROD

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    if ~exist('safe','var')
        safe = true;
    end

    if safe
        x = orthogonalize(x, x.order );
        res = norm( x.U{end}(:) );
    else
        res = sqrt(innerprod( x, x ));

        if res < 1e-7
            x = orthogonalize(x, x.order );
            res = norm( x.U{end}(:) );
        end
    end

end

author: NicolasBoumal

repo: manopt


text:
function p4 = angle_half_2d ( p1, p2, p3 )

%*****************************************************************************80
%
%% ANGLE_HALF_2D finds half an angle in 2D.
%
%  Discussion:
%
%    The original angle is defined by the sequence of points P1, P2 and P3.
%
%    The point P4 is calculated so that:
%
%      (P1,P2,P4) = (P1,P2,P3) / 2
%
%        P1
%        /
%       /   P4
%      /  .
%     / .
%    P2--------->P3
%
%    Thanks to Cesar Fraga Bobis for pointing out a typographical error in
%    a previous version of this routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), P3(2,1), points defining the angle.
%
%    Input, real P4(2,1), a point defining the half angle.
%    The vector P4 - P2 will have unit norm.
%
  p4(1:2,1) = 0.5 * ( ...
      ( p1(1:2,1) - p2(1:2,1) ) / sqrt ( sum ( ( p1(1:2,1) - p2(1:2,1) ).^2 ) ) ...
    + ( p3(1:2,1) - p2(1:2,1) ) / sqrt ( sum ( ( p3(1:2,1) - p2(1:2,1) ).^2 ) ) );

  p4(1:2,1) = p2(1:2,1) + p4(1:2,1) / sqrt ( sum ( ( p4(1:2,1) ).^2 ) );

  return
end


text:
function [ebsd,interface,options] = loadEBSD(fname,varargin)

warning('loadEBSD is depreciated. Please use instead EBSD.load');
[ebsd,interface,options] = EBSD.load(fname,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function [stat, cfg] = ft_statistics_mvpa(cfg, dat, design)

% FT_STATISTICS_MVPA performs multivariate pattern classification
% on the data. If the data has not been averaged over time, classification
% is performed separately for every time point. Additionally, searchlight
% analysis (classification for each channel/voxel
% separately) or time x time generalisation can be performed.
%
% This function should not be called directly, instead
% you should call the function that is associated with the type of
% data on which you want to perform the test.
%
% Use as
%   stat = ft_timelockstatistics(cfg, data1, data2, data3, ...)
%   stat = ft_freqstatistics    (cfg, data1, data2, data3, ...)
%   stat = ft_sourcestatistics  (cfg, data1, data2, data3, ...)
% where the data is obtained from FT_TIMELOCKANALYSIS, FT_FREQANALYSIS
% or FT_SOURCEANALYSIS respectively, or from FT_TIMELOCKGRANDAVERAGE,
% FT_FREQGRANDAVERAGE or FT_SOURCEGRANDAVERAGE respectively.
%
% The configuration can contain
%   cfg.searchlight     = 'yes' or 'no', performs searchlight analysis
%                         (default 'no'). More information see below
%   cfg.timextime       = 'yes' or 'no', performs time x time
%                         generalisation. In other words, the classifier
%                         is trained at each time point and tested at
%                         every time point. The result is a time x time
%                         matrix of classification performance.
%                         (default 'no')
%                         Note that searchlight and timextime cannot be
%                         run simultaneously (at least one option needs
%                         to be set to 'no').
%   cfg.mvpa            = structure that contains detailed options for the
%                         MVPA procedure. See
%                         https://github.com/treder/MVPA-Light for more
%                         details.
%
%
%   cfg.mvpa.classifier      = 'lda'          Regularised linear discriminant analysis
%                                        (LDA) (for two classes)
%                         'multiclass_lda' LDA for more than two classes
%                         'logreg'       Logistic regression
%                         'svm'          Support Vector Machine (SVM)
%                         'ensemble'     Ensemble of classifiers. Any of the other
%                                        classifiers can be used as a learner.
%                         'kernel_fda'   Kernel Fisher Discriminant Analysis
%   cfg.mvpa.metric          = string, performance metric. Possible metrics:
%                         accuracy auc tval dval confusion precision
%                         recall f1
%   See https://github.com/treder/MVPA-Light for an overview of all
%   classifiers and metrics.
%
%   cfg.mvpa.param           = struct, structure with hyperparameters for the
%                         classifier (see HYPERPARAMETERS below)
%
%   cfg.mvpa.balance         = string, for imbalanced data that does not have
%                         the same number of instances in each class
%                         'oversample'     oversamples the minority classes
%                         'undersample'    undersamples the minority classes
%                         such that all classes have the same number of
%                         instances. Note that undersampling is at the
%                         level of the repeats, whereas we oversampling occurs within each
%                         training set (for an explanation see mv_balance_classes).
%                         You can also give an integer number for undersampling.
%                         The samples will be reduced to this number. Note that
%                         concurrent over/undersampling (oversampling of the
%                         smaller class, undersampling of the larger class) is not
%                         supported at the moment
%  cfg.mvpa.replace          = bool, if balance is set to 'oversample' or 'undersample',
%                         replace determines whether data is drawn with
%                         replacement (default 1)
%  cfg.mvpa.normalise        = string, normalises the data across samples, for each time point
%                         and each feature separately, using 'zscore' or 'demean'
%                         (default 'zscore'). Set to 'none' or [] to avoid normalisation.
%  cfg.mvpa.feedback         = 'yes' or 'no', whether or not to print feedback on the console (default 'yes')
%
% To obtain a realistic estimate of classification performance,
% cross-validation is used. It is controlled by the following parameters:
%   cfg.mvpa.cv              = string, cross-validation type, either 'kfold', 'leaveout'
%                         or 'holdout'. If 'none', no cross-validation is
%                         used and the classifier is tested on the training
%                         set. (default 'kfold')
%   cfg.mvpa.k               = number of folds in k-fold cross-validation (default 5)
%   cfg.mvpa.repeat          = number of times the cross-validation is repeated
%                         with new randomly assigned folds (default 5)
%   cfg.mvpa.p               = if cfg.cv is 'holdout', p is the fraction of test
%                         samples (default 0.1)
%   cfg.mvpa.stratify        = if 1, the class proportions are approximately
%                         preserved in each test fold (default 1)
%
%
% More information about each classifier is found in the documentation of
% MVPA-Light (github.com/treder/MVPA-Light/).
%
% HYPERPARAMETERS:
% Each classifier comes with its own set of hyperparameters, such as
% regularisation parameters and the kernel. Hyperparameters can be set
% using the cfg.param substruct. For instance, in LDA, cfg.param.lambda =
% 'auto' sets the lambda regularisation parameter.
%
% The specification of the hyperparameters is found in the training function
% for each classifier at github.com/treder/MVPA-Light/tree/master/classifier
% If a hyperparameter is not specified, default values are used.
%
% SEARCHLIGHT ANALYSIS:
% Classification using a feature searchlight approach can highlight which
% feature(s) are informative. To this end, classification is performed on
% each feature separately. However, neighbouring features can enter the
% classification together when specified. Optional parameters:
%
% cfg.mvpa.neighbours   = neighbourhood structure, see FT_PREPARE_NEIGHBOURS
%                         Alternatively, a [features x features] matrix specifying
%                         which features are neighbours of each other. This
%                         matrix can be either
%                         a GRAPH consisting of 0's and 1's. A 1 in the
%                         (i,j)-th element signifies that feature i and feature j
%                         are neighbours, and a 0 means they are not neighbours
%                         or a DISTANCE MATRIX, where larger values mean larger distance.
%                         If no matrix is provided, every feature is only neighbour
%                         to itself and classification is performed for each feature
%                         separately.
% cfg.mvpa.size         = size of the 'neighbourhood' of a feature.
%                         number of steps taken through the neighbourhood
%                         to include neighbours
%                         0: only the feature itself is considered (no neighbours)
%                         1: the feature and its immediate neighbours
%                         2: the feature, its neighbours, and its neighbours'
%                         neighbours
%                         3+: neighbours of neighbours of neighbours etc
%                         (default 1)
%                         if cfg.neighbours is a distance matrix, size defines the number of
%                         neighbouring features that enter the classification
%                         0: only the feature itself is considered (no neighbours)
%                         1: the feature and its first closest neighbour
%                            according to the distance matrix
%                         2+: the 2 closest neighbours etc.
%
% -- TODO: for time x time generalisation, in MVPA light we can use two
% different datasets (one for training the classifier, the other one for
% testing). This could be realised eg using extra fields in cfg such as
% cfg.X2 for dataset 2 and cfg.design2 for the second design matrix
%
% Returns:
%   stat        = struct with results. the .metric field contains the
%                   requested metrics
%

% Copyright (C) 2019, Matthias Treder
%
% 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$

ft_hastoolbox('mvpa-light', 1);

% do a sanity check on the input data
assert(isnumeric(dat),    'this function requires numeric data as input, you probably want to use FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS instead');
assert(isnumeric(design), 'this function requires numeric data as input, you probably want to use FT_TIMELOCKSTATISTICS, FT_FREQSTATISTICS or FT_SOURCESTATISTICS instead');

% cfg: set defaults
cfg.searchlight     = ft_getopt(cfg, 'searchlight', 'no');
cfg.timextime       = ft_getopt(cfg, 'timextime',   'no');
cfg.mvpa            = ft_getopt(cfg, 'mvpa',        []);
cfg.mvpa.metric     = ft_getopt(cfg.mvpa, 'metric',   'accuracy');
cfg.mvpa.feedback   = ft_getopt(cfg.mvpa, 'feedback',   'yes');
cfg.mvpa.neighbours = ft_getopt(cfg.mvpa, 'neighbours', []);

% flip dimensions such that the number of trials comes first
dat = dat';

% if cfg.dim has two entries which are non-singleton then the data has been
% 3D before. We must reshape it from 2D to 3D to run it in MVPA-Light
data_is_3D = (numel(cfg.dim) > 1 && all(cfg.dim > 1));  % checks whether data was 3D before being reshaped
if data_is_3D
  dat = reshape(dat, size(dat,1), cfg.dim(1), cfg.dim(2));
end

y = design;

%% perform sanity checks on parameters
if istrue(cfg.timextime) && istrue(cfg.searchlight)
  ft_error('you should not set timextime = ''yes'' and searchlight = ''yes'' simultaneously')
end

% timextime = 'yes' but data is not 3D we should change timextime to 'no'
if istrue(cfg.timextime) && ~data_is_3D
  ft_warning('timextime = ''yes'' but data has no time dimension, setting timextime = ''no''');
  cfg.timextime = 'no';
end

label = [];
dim   = [];
dimord = [];

%% Call MVPA-Light

if strcmp(cfg.searchlight, 'yes')
  % --- searchlight analysis ---

  if isstruct(cfg.mvpa.neighbours)
    cfg.mvpa.neighbours = channelconnectivity(struct('neighbours',cfg.mvpa.neighbours, 'channel', {cfg.channel}));
  end
  [perf, result] = mv_searchlight(cfg.mvpa, dat, y);

  % this preserves any spatial dimension, so no adjustment is done to a
  % channel list, if present
  if isfield(cfg, 'channel')
    label = cfg.channel;
  end

  if isfield(cfg, 'dim')
    dim = cfg.dim;
  end

elseif strcmp(cfg.timextime, 'yes')
  % --- time x time generalisation ---
  [perf, result] = mv_classify_timextime(cfg.mvpa, dat, y);

  % this does note preserve any spatial dimension, so label should be
  % adjusted
  label = squeezelabel(label, cfg);
  dim   = squeezedim(dim, cfg);
  dimord = 'time_time';

elseif data_is_3D
  % --- classification across time ---
  [perf, result] = mv_classify_across_time(cfg.mvpa, dat, y);

  % this does note preserve any spatial dimension, so label should be
  % adjusted
  label = squeezelabel(label, cfg);
  dim   = squeezedim(dim, cfg);

else
  % --- data has no time dimension, perform only cross-validation ---
  [perf, result] = mv_crossvalidate(cfg.mvpa, dat, y);

  % this does note preserve any spatial dimension, so label should be
  % adjusted
  label = squeezelabel(label, cfg);
  dim   = squeezedim(dim, cfg);

end

%% setup stat struct
stat = [];
if ~iscell(cfg.mvpa.metric), cfg.mvpa.metric = {cfg.mvpa.metric}; end
if ~iscell(perf),            perf            = {perf};            end
for mm=1:numel(perf)

  % Performance metric
  stat.(cfg.mvpa.metric{mm}) = perf{mm};

  % Std of performance
  if iscell(result.perf_std)
    stat.([cfg.mvpa.metric{mm} '_std']) = result.perf_std{mm};
  else
    stat.([cfg.mvpa.metric{mm} '_std']) = result.perf_std;
  end
end

% return the MVPA-Light result struct as well
stat.mvpa = result;

if ~isempty(label)
  stat.label = label;
end

if ~isempty(dim)
  stat.dim = dim;
end

if ~isempty(dimord)
  stat.dimord = dimord;
end

function label = squeezelabel(label, cfg)

if isfield(cfg, 'channel')
  label = sprintf('combined(%s)', sprintf('%s',cfg.channel{:}));
end

function dim = squeezedim(dim, cfg)

if isfield(cfg, 'dim')
  dim = cfg.dim;
  dim(1) = 1;
end


text:
function [] = aps_modis(start_step,end_step)
% [] = aps_modis(start_step,end_step)
% Computes the MODIS delays based on the OSCAR service provied by JPL
%
%     Copyright (C) 2015  Bekaert David - University of Leeds
%     Email: eedpsb@leeds.ac.uk or davidbekaert.com
%
%     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.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
% By David Bekaert - University of Leeds
% May 2014
%
% Modifications:
% DB	05/2014 	Start from the MERIS scrip the program the function
% DB    07/2014     Have re-visualization of sounding data when already
%                   processed
% DB    08/2014     Include SAR varying conversion factors
% DB    08/2014     Include option to use MODIS recalibrated data
% DB    09/2014     Check if conversion factors are known before
% DB    04/2015     Fix in case filename grd is not starting with OSCAR_Modis*grd
% DB    02/2016     Include aps_systemcall and suppress filename output.

%% loading the StaMPS variables
curdir = pwd;

if start_step==0
   % get the data from the OSCAR webite (JPL)
   aps_modis_Python
end

if start_step==1
     % Compute the MERIS conversion factors from the sounding data
    fprintf('\n\nStep 1: Compute the spectrometer conversion factors from the sounding data \n')

    % getting parameters from parms_aps file
    sounding_dir = getparm_aps('sounding_dir',1);
    start_date = getparm_aps('sounding_start_date',1);
    end_date = getparm_aps('sounding_end_date',1);
    time_stamp = getparm_aps('sounding_time_stamp',1);
    sounding_ifg_dates = getparm_aps('sounding_ifg_dates',1);
    start_year_str = start_date(1:4);
    end_year_str = end_date(1:4);
    start_str = start_date(5:6);
    end_str = end_date(5:6);
    time_stamp_str = [];
    for k=1:size(time_stamp,1)
        if k>1
            time_stamp_str = [time_stamp_str '_' time_stamp(k,:)];
        else
            time_stamp_str = [time_stamp(k,:)];
        end
    end

     % see if the filename already exist
    if strcmp(sounding_ifg_dates,'y')
        savename = [sounding_dir filesep 'Spectrometer' filesep 'Spectrometer_sensitivity_SAR_dates_1month_' time_stamp_str 'Hr_' num2str(start_year_str) start_str '_' num2str(end_year_str) end_str '.mat'];
    else
        savename = [sounding_dir filesep 'Spectrometer' filesep 'Spectrometer_sensitivity_' time_stamp_str 'Hr_' num2str(start_year_str) start_str '_' num2str(end_year_str) end_str '.mat'];
    end

    if exist(savename,'file')==2
        % check if user wants to visualize or re-run the data
        str = '';
        while strcmpi(str,'y')~=1 && strcmpi(str,'n')~=1
            str = input(['This file already exist, do you want to re-process the data? [y/n] \n'],'s');
        end

        if strcmpi(str,'y')
            rerun_flag =1;
        else
            sounding_spectrometer_sens_display
            rerun_flag =0;
        end
    else
       rerun_flag =1;
    end

    if rerun_flag ==1;
        sounding_spectrometer_sens;
    end
end

if start_step<=2 && end_step >=2
    % SAR delays using MODIS data
    fprintf('\n\nStep 2: Compute MODIS tropospheric delay for individual dates\n')

    % generating a file list of the files to be processed:
    modis_datapath = getparm_aps('modis_datapath',1);
    command_str = 'echo files > modis_batch_file.txt';
    command_str2 = ['ls -d ' modis_datapath filesep '2*' filesep '*.grd >> modis_batch_file.txt'];

    aps_systemcall(command_str);
    aps_systemcall(command_str2);

    % running the MODIS computation on the file list
    aps_modis_SAR('modis_batch_file.txt')
end
if start_step<=3 && end_step >=3
    fprintf('\n\nStep 3: Computes MODIS tropospheric delay for interferograms\n')
    aps_modis_InSAR
end

cd(curdir)

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function H = hopkins (dataIn,xlim1, xlim2, ylim1, ylim2, m)
% HOPKINS calculates H for Hopkins statistics estimnation
%H = hopkins (dataIn,xlim1, xlim2, ylim1, ylim2, m)

sizeData = size(dataIn);
nPoints = sizeData(1); %number of datapoints
sizeVec = [xlim2-xlim1, ylim2 - ylim1];

%generate m random sampling points (not data points)
sp(:,1) = xlim1 + sizeVec(1)*rand(m,1);
sp(:,2) = ylim1 + sizeVec(2)*rand(m,1);

%select m random data points:
pointIndex = [];
sdp = [];
ii = 1;
while ii <= m
    rn = round(nPoints * rand);
    if ~any(pointIndex == rn)
        if ~rn==0
            pointIndex(ii) = rn;
            sdp(ii,:) = dataIn (rn,:);
            ii = ii + 1;
        end
    end
end

% minimal distance between random sampling points and data points:
U = sum(min(dist2 (dataIn, sp), [], 1)); %minimum distances

leaveIndex = removerows((1:nPoints)',pointIndex); %excludes data points in sdp
dataR = dataIn(leaveIndex,:);
W = sum(min(dist2 (dataR, sdp), [], 1));

H = U/(U+W);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
clear
clc
im = testmat(109,105);
a = 0:17:179;
[W5,p5,t5] = build5(im,a);
[W7,p7,t7] = build7(im,a);
[W8,p8,t8] = build8(im,a);
fprintf(['\nElapsed time:\nbuild5 %f seconds'...
    '\nbuild7 %f seconds\nbuild8 %f seconds\r'],...
    t5,t7,t8);
%Variance
% var_78 = var(W7(:)-W8(:));
% fprintf('\nVariance8: %f\r',var_78);

author: Sable


text:
function [ fHdl ] = plot( obj, subjects )
% Plot cluster and subject-level estimation result from HUGE model.
%
% INPUTS:
%   obj - A tapas_Huge object containing estimation results.
%
% OPTIONAL INPUTS:
%   subjects - A vector containing indices of subjects for which to plot
%              detailed results.
%
% OUTPUTS:
%   fHdl - Handle of first figure.
%
% See also tapas_Huge.ESTIMATE
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch)
% Copyright (C) 2019 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This software is provided "as is", without warranty of any kind, express
% or implied, including, but not limited to the warranties of
% merchantability, fitness for a particular purpose and non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

if nargout > 0
    fHdl = figure( );
else
    figure( );
end
if nargin < 2
    subjects = [];
end

tickLabels = obj.parse_labels( obj.dcm, obj.labels, obj.idx );

%% assignments/boxplot
% figure;
if obj.K > 1
    % subject assignment
    bar(obj.posterior.q_nk,'stacked');
    axis([0 obj.N+1 0 1])
    title('assignments')
    ylabel('q_{nk}')
    xlabel('subject index')
    title('Assignment');
else
    % boxplot of MAP estimates of DCM parameters
    hold on
    line([0 obj.idx.P_c + obj.idx.P_h], [0 0], 'color', 'k')
    boxplot(obj.posterior.mu_n)
    ylabel('\mu_n')
    title('Empirical Bayes');
    set(gca,'XTick',1:obj.idx.P_c + obj.idx.P_h, 'XTickLabelRotation', 60, ...
        'XTickLabel', tickLabels, 'TickLabelInterpreter', 'none');
end

%% cluster
figure
hold on
% plot posterior cluster mean and 95% marginal credible intervals
legends = cell(obj.K,1);
for k = 1:obj.K
    xOffset = ((k-1)/max(1,(obj.K-1)) - .5)/4;
    switch obj.posterior.method
        case 'VB'
            clMean = obj.posterior.m(k,:);
            clStd = sqrt(diag(obj.posterior.S(:,:,k))'/...
                (obj.posterior.tau(k)*...
                (obj.posterior.nu(k) - obj.idx.P_c + 1)));
            s = tinv(1-0.025, obj.posterior.nu(k));
            errorbar((1:obj.idx.P_c) + xOffset, clMean, s*clStd, 'd');
        case 'MH'
            clMean = obj.posterior.mean.mu(k,:);
            [~,i1] = min(abs(obj.posterior.quantile.levels - .025));
            [~,i2] = min(abs(obj.posterior.quantile.levels - .975));
            neg = clMean - obj.posterior.quantile.mu(k,:,i1);
            pos = obj.posterior.quantile.mu(k,:,i2) - clMean;
            errorbar((1:obj.idx.P_c) + xOffset, clMean, neg, pos, 'd');
    end
    legends{k} = ['cluster ' num2str(k)];
end
line([0 obj.idx.P_c+1], [0 0], 'color', 'k')
xlim([0 obj.idx.P_c+1])
ylabel('\mu_k')
set(gca,'XTick',1:obj.idx.P_c, 'XTickLabelRotation', 60, ...
    'XTickLabel', tickLabels(1:obj.idx.P_c), 'TickLabelInterpreter', 'none');
legend(legends)
title('Clusters')

%% DCM
% plot 25 samples from posterior over (noise free) BOLD response
nSmp = 25;
for n = subjects(:)'
    figure
    hold on
    % draw samples from posterior over DCM parameters
    switch obj.posterior.method
        case 'VB'
            postMean = obj.posterior.mu_n(n,:);
            postStd = chol(obj.posterior.Sigma_n(:,:,n));
            postSmp = randn(nSmp,obj.idx.P_c + obj.idx.P_h);
            postSmp = bsxfun(@plus, postSmp*postStd, postMean);
        case 'MH'
            nTrace = length(obj.trace.smp);
            nSmp = min(nSmp, nTrace);
            idx = randsample(nTrace, nSmp);
            tmp = [reshape([obj.trace.smp(idx).theta_c], obj.N, obj.idx.P_c, []), ...
                reshape([obj.trace.smp(idx).theta_h], obj.N, obj.idx.P_h, [])];
            postSmp = permute(tmp(n,:,:), [3 2 1]);
    end
    legends = {'measured'};
    plot(obj.data(n).bold(:),'k')
    % plot ground truth if available
    if ~isempty(obj.model)
        [ ~, epsilon ] = obj.gen_bold( n, obj.model.theta(n,:) );
        plot(obj.data(n).bold(:) - epsilon(:),'r')
        legends = [legends, {'ground truth'}]; %#ok<AGROW>
    end
    legends = [legends ,{'posterior samples'}]; %#ok<AGROW>
    for iSmp = 1:nSmp
        [ ~, epsilon ] = obj.gen_bold( n, postSmp(iSmp,:) );
        plot(obj.data(n).bold(:) - epsilon(:),'b')
    end
    legend(legends);
    title(['Subject ' num2str(n)]);
    ylabel('BOLD')
    xlabel('sample index')
end

end

author: translationalneuromodeling

repo: tapas


text:
function b = isint(x)
%ISINT  True for integers.
%
%   ISINT(X) returns 1's where the elements of X are integers and 0's where
%   they are not.  For example, ISINT([ 2e9 pi 3+5i NaN ]) is [ 1 0 1 0 ].
%
%   See also ISEVEN, ISODD.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-04-12 14:30:59 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   error(nargchk(1, 1, nargin));
   if ~isnumeric(x)
      error('Argument must be a numeric array.');
   end

   cls = class(x);                      % class of input argument
   if isempty(x)
      b = feval(cls, x);                % return empty array of same class
   else
      switch cls
         case 'double'
            b = x == round(x);
         case 'single'
            % "mod" is not defined for class "single"; so convert input to
            % double, compare, and convert back
            d = double(x);
            b = single(d == round(d));
         case {'int8', 'int16', 'int32', 'int64', ...
               'uint8', 'uint16', 'uint32', 'uint64'}
            % return an array of ones of the same class as input argument
            b = logical(repmat(feval(cls, 1), size(x)));
         otherwise
            error('Argument is of unrecognized class.');
      end
   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function BandTrader(sys)
%BANDTRADER - Band Trading Strategy
%
%  This MATLAB function implements a simple band trading strategy. A band
%  consists of two lines that form the upper and lower boundaries of the band.
%  The upper and lower boundaries are used to to enter and exit trades.
%  For example, if prices fall below the lower boundary a buy signal is
%  generated.
%
%  Requirements:
%  - MATLAB 2012b
%  - TA Developer Toolbox (http://www.tadeveloper.com)
%
%  Please refer to http://www.tadeveloper.com/docs/80-matlab-algo-trading
%  for a full description of the strategy
%

% custom trading parameters
if sys.TradingParameters.isKey('movingAverage'),
    movingAverage = sys.TradingParameters('movingAverage');
else
    movingAverage = 10; % default value if unset
end;
if sys.TradingParameters.isKey('bandWidth'),
    bandWidth = sys.TradingParameters('bandWidth');
else
    bandWidth = 2; % default value if unset
end;

% Construct the center line
Average = (sys.High + sys.Low)/2;
CenterLine = Ema(Average, movingAverage);
PlotDataSeries(sys, CenterLine, 'blue', 'solid', 2, '');

% Construct the bands
AverageTrueRange = Atr(sys.BarData, movingAverage);
UpperBand = CenterLine + bandWidth * AverageTrueRange;
LowerBand = CenterLine - bandWidth * AverageTrueRange;
PlotDataSeries(sys, UpperBand, 'magenta', 'solid', 2, '');
PlotDataSeries(sys, LowerBand, 'magenta', 'solid', 2, '');

% Entry and Exit Signals
LongSignal = sys.Close < LowerBand;
ShortSignal = sys.Close > UpperBand;
AddLongSignal(sys, LongSignal);
AddShortSignal(sys, ShortSignal);

end

repo: mcbench-benchmarks


function [acc, pred] = CrossValidateKNN(y, X, tCL, k, knn_size);
% [acc, pred] = CrossValidateKNN(y, X, tCL, k, knn_size);
%
% Cross-validation for evaluating the k-nearest neighbor classifier with
% a learned metric.  Performs k-fold cross validation, training on the
% training fold and evaluating on the test fold
%
% y: (n x 1) true labels
%
% X: (n x m) data matrix
%
% tCL: Metric learning algorithm that takes in true labels as first
% argument, and data as a second
%
% k: Number of cross-validated folds
%
% knn_size: size of nearest neighbor window
%
% Returns
% acc: cross-validated accuracy
% pred: predictions on test set for each row in X

[n,m] = size(X);
if (n ~= length(y)),
   disp('ERROR: num rows of X must equal length of y');
   return;
end

%permute the rows of X and y
rp = randperm(n);
y = y(rp);
X = X(rp, :);

pred = zeros(n,1);
for (i=1:k),
   test_start = ceil(n/k * (i-1)) + 1;
   test_end = ceil(n/k * i);

   yt = [];
   Xt = zeros(0, m);
   if (i > 1);
       yt = y(1:test_start-1);
       Xt = X(1:test_start-1,:);
   end
   if (i < k),
       yt = [yt; y(test_end+1:length(y))];
       Xt = [Xt; X(test_end+1:length(y), :)];
   end

   nt = length(yt);
   yt = yt(1:nt);
   Xt = Xt(1:nt, :);

   %train model
   M = feval(tCL, yt, Xt);

   %evaluate model
   XT = X(test_start:test_end, :);
   yT =  y(test_start:test_end);
   pred(test_start:test_end) = KNN(yt, Xt, sqrtm(M), knn_size, XT);
end
acc = sum(pred==y)/n;


text:
% This file is designed for customer testing. Users can put their target
% images into the folder "dataset/test/images". This file supports testing
% images with arbitrary resolution

function testGAIC_qualitative_customer()

imDir = dir(fullfile('dataset','test','images','*.*'));
imDir = imDir(3:end);
imList = fullfile('dataset','test','images',{imDir.name});

load(['imdb_GAIC1000.mat'],'imdb');
gt_scores = cat(2,imdb.bbox.gt_scores{imdb.images.set==1});
gt_scores_means = mean(gt_scores);
gt_scores_stds = std(gt_scores);

netStruct = load(fullfile('data','pretrained_models','net-epoch-37.mat'),'net');
net = dagnn.DagNN.loadobj(netStruct.net) ;
net.mode = 'test' ;
net.move('gpu') ;
probVarI = net.getVarIndex('predcls');
net.vars(probVarI).precious = 1;
minScale = net.meta.minScale;

for aspect_ratios = {'1_1','4_3','16_9','2_3','4_5','5_7'}

    for i = 1:numel(imList)
        try
            img = imread(imList{i});
        catch
            continue
        end

        [x1,x2,x3] = size(img);

        imre = single(imresize(img,minScale/min([x1,x2]),'bilinear'));
        imre = bsxfun(@minus,imre,net.meta.normalization.averageImage);

        [r1,r2,r3] = size(imre);
        r1 = 32*round(r1/32);
        r2 = 32*round(r2/32);
        imre = imresize(imre,[r1,r2],'bilinear');

        boxes = generateBoxes(imre,aspect_ratios{1});
        if isempty(boxes)
            continue;
        end

        inputs = {'input', gpuArray(imre), 'rois', gpuArray(single([ones(1,size(boxes,2));boxes]))} ;

        net.eval(inputs) ;
        preds = squeeze(gather(net.vars(probVarI).value)) ;
        preds = preds * gt_scores_stds + gt_scores_means;
        [~,id_preds] = sort(preds,'descend');

        predBox = boxes(:,id_preds(1));

        scale1 = x1/r1;
        scale2 = x2/r2;
        boxes_s(1,:) = max(floor(predBox(1,:) * scale1),1);
        boxes_s(2,:) = max(floor(predBox(2,:) * scale2),1);
        boxes_s(3,:) = min(ceil(predBox(3,:) * scale1),x1);
        boxes_s(4,:) = min(ceil(predBox(4,:) * scale2),x2);

        outDir = fullfile('dataset','test','result');
        mkdir(outDir);

        imwrite(img(boxes_s(1):boxes_s(3),boxes_s(2):boxes_s(4),:),...
            fullfile(outDir,[imDir(i).name(1:end-4) '_preds_' aspect_ratios{1} '.jpg']));

    end

end

function boxes = generateBoxes(im,aspect_ratios)

switch aspect_ratios
    case '1_1'
        boxes = generateBoxes_1_1(im);
    case '4_3'
        boxes = generateBoxes_4_3(im);
    case '16_9'
        boxes = generateBoxes_16_9(im);
    case '2_3'
        boxes = generateBoxes_2_3(im);
    case '4_5'
        boxes = generateBoxes_4_5(im);
    case '5_7'
        boxes = generateBoxes_5_7(im);
    case '1_2'
        boxes = generateBoxes_1_2(im);
    case '1_3'
        boxes = generateBoxes_1_3(im);
    case '1_4'
        boxes = generateBoxes_1_4(im);
    case '2_1'
        boxes = generateBoxes_2_1(im);
    case '3_1'
        boxes = generateBoxes_3_1(im);
    case '4_1'
        boxes = generateBoxes_4_1(im);
end
boxes = round(boxes);

function boxes = generateBoxes_16_9(im)

boxes = [];
[s1,s2,s3] = size(im);
id = floor(max(s1,s2)/8);
out_Dim = [16*(14:id)',9*(14:id)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);

cnt = 0;
if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:16:(s1-out_Dim(scale,2)+1);
    step_y = 1:9:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_4_3(im)

boxes = [];
[s1,s2,s3] = size(im);
id = floor(max(s1,s2)/8);
out_Dim = [8*(25:id)',6*(25:id)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);

cnt = 0;
if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:16:(s1-out_Dim(scale,2)+1);
    step_y = 1:12:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_1_1(im)

boxes = [];
[s1,s2,s3] = size(im);

out_Dim = [4*(40:2:64)',4*(40:2:64)'];

cnt = 0;
if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:12:(s1-out_Dim(scale,2)+1);
    step_y = 1:12:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_2_3(im)

boxes = [];

[s1,s2,s3] = size(im);
id = floor(max(s1,s2)/8);
out_Dim = [8*(18:id)',12*(18:id)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);

cnt = 0;
if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:8:(s1-out_Dim(scale,2)+1);
    step_y = 1:12:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_4_5(im)

boxes = [];

[s1,s2,s3] = size(im);
id = floor(max(s1,s2)/8);
out_Dim = [8*(20:id)',10*(20:id)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);

if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2

for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end

end

cnt = 0;
for scale = 1:size(out_Dim,1)
    step_x = 1:8:(s1-out_Dim(scale,2)+1);
    step_y = 1:10:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_5_7(im)

boxes = [];

[s1,s2,s3] = size(im);
id = floor(max(s1,s2)/8);
out_Dim = [10*(15:id)',14*(15:id)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);
cnt = 0;

if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:8:(s1-out_Dim(scale,2)+1);
    step_y = 1:14:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_1_2(im)

boxes = [];

[s1,s2,s3] = size(im);
% id = floor(max(s1,s2)/8);
out_Dim = [6*(18:100)',12*(18:100)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);

cnt = 0;
if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:8:(s1-out_Dim(scale,2)+1);
    step_y = 1:16:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_1_3(im)

boxes = [];

[s1,s2,s3] = size(im);
out_Dim = [4*(18:100)',12*(18:100)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);

cnt = 0;
if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:8:(s1-out_Dim(scale,2)+1);
    step_y = 1:16:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_1_4(im)

boxes = [];

[s1,s2,s3] = size(im);
out_Dim = [3*(18:100)',12*(18:100)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);

cnt = 0;
if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:8:(s1-out_Dim(scale,2)+1);
    step_y = 1:16:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_2_1(im)

boxes = [];

[s1,s2,s3] = size(im);
out_Dim = [12*(18:100)',6*(18:100)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);

cnt = 0;
if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:16:(s1-out_Dim(scale,2)+1);
    step_y = 1:8:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_3_1(im)

boxes = [];

[s1,s2,s3] = size(im);
out_Dim = [12*(18:100)',4*(18:100)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);

cnt = 0;
if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:16:(s1-out_Dim(scale,2)+1);
    step_y = 1:8:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

function boxes = generateBoxes_4_1(im)

boxes = [];

[s1,s2,s3] = size(im);
out_Dim = [12*(18:100)',3*(18:100)'];
out_Dim = out_Dim(out_Dim(:,1)<=s2,:);
out_Dim = out_Dim(out_Dim(:,2)<=s1,:);

cnt = 0;
if out_Dim(end,1)*out_Dim(end,2)>0.4*s1*s2
for scale = 1:size(out_Dim,1)
    if out_Dim(1,1)*out_Dim(1,2)<0.4*s1*s2
        out_Dim = out_Dim(2:end,:);
    end
end
end

for scale = 1:size(out_Dim,1)
    step_x = 1:16:(s1-out_Dim(scale,2)+1);
    step_y = 1:8:(s2-out_Dim(scale,1)+1);
    for x1 = step_x
        for y1 = step_y
            cnt = cnt + 1;
            boxes(1,cnt) = x1;
            boxes(2,cnt) = y1;
            boxes(3,cnt) = x1+out_Dim(scale,2)-1;
            boxes(4,cnt) = y1+out_Dim(scale,1)-1;
        end
    end
end
boxes = round(boxes);

author: HuiZeng

repo: Grid-Anchor-based-Image-Cropping


text:
function colormap_centered(hax)

if nargin < 1
  hax = gca();
end

for n=1:length(hax)
  cl = get(hax(n), 'clim');
  l = max(abs(cl));
  cl = [-l l];
  set(hax(n), 'clim', cl);
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [us] = ns2us(ns)
% Convert time from nanoseconds to microseconds.
% Chad Greene 2012
us = ns*0.001 ;


text:
function scale = p23_scale ( neqn )

%*****************************************************************************80
%
%% P23_SCALE returns scale factors for problem p23.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = [ 1.0; 1.0 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% generate new location for the K users
close all
clear all

K=4;
%%
Pt=zeros(K,2);
%%
Lroom=200;
Wroom=30;
k1=[1,0];
k2=[0,1];
R=10;
%%
for k0=1:K
    r=rand(1,1)*R;
    theta=rand(1,1)*2*pi;
    px=r*cos(theta);
    py=r*sin(theta);
    pt=[Lroom,Wroom]+px*k1+py*k2;
    Pt(k0,:)=pt;
end
%%
save('user_location.mat','K','Pt','Lroom','Wroom','R');
figure
plot(Pt(:,1),Pt(:,2),'ro');
xlim([Lroom-R,Lroom+R]);ylim([Wroom-R,Wroom+R]);
hold on
theta=linspace(0,1,100).*2.*pi;
hold on
plot(Lroom+R*cos(theta),Wroom+R*sin(theta),'r.')

author: guohuayan