text:
function h = conv(f, g)
%CONV   Convolution of BNDFUN objects.
%   H = CONV(F, G) produces the convolution of BNDFUN objects F and G:
%                     -
%                    /
%           H(x) =   |    F(t) G(x-t) dt,  x in [a + c, b + d]
%                    /
%                   -
%   where domain(F) is [a, b] and domain(G) is [c, d]. The integral is taken
%   over all t for which the integrand is defined: max(a, x - d) <= t <= min(b,
%   x - c).  The breakpoints of H are all pairwise sums of the breakpoints of F
%   and G.
%
%   Note that CONV only supports piecewise-smooth functions on bounded domains.
%
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.
%
% Nick Hale and Alex Townsend, 2014

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEVELOPER NOTE:
%
% For further details, see Hale and Townsend, "An algorithm for the convolution
% of Legendre series", SIAM Journal on Scientific Computing, Vol. 36, No. 3,
% pages A1207-A1220, 2014.
%
% In the following, it is assumed that the length of the domain of g is greater
% than the length of the domain of f. If this is not the case, then simply
% compute h = conv(g, f), which is equivalent. We assume f and g are polynomials
% of degree M and N, resepectively. If f and g are piecewise-defined, one can
% use the bilinearity of convolution and convolve each each the FUNs
% individually.
%
% The general convolution domain (for smooth functions) is as follows:
%           .___________________________
%          /|                  |      /
%        /  |                  |    /
%      /    |                  |  /
%    /______|__________________|/
%  a+c     b+c                a+d     b+d
%
% The triangular pieces at end are dealt with using a convolution theorem for
% Legendre polynomials, which leads to a convenient recurrence relation. The
% cost of this is O(m*n) and results in a polynomial of degree m+n. See
% EASYCONV() for details.
%
% Similarly one can show the interior rectangle results in a polynomial of
% degree n. This can be computed by patching with R = floor[(d-c) / (b-a)]
% parallelograms (although triangle Z is not used). Each parallelogram requires
% restricting g to a suitable subdomain, which costs O(n^2) operations. Total
% complexity ratio*(m*n + n*n)
%            ___________________________
%          /       /       /:     /   /
%        /       /       /  : Z /   /          <-- R patches
%      /       /       /    : /   /
%    /_______/_______/______/__ /
%  a+c     b+c             fl  a+d     b+d
%
% The final piece is computed via further parallelogram subdivision starting
% from the right (B), and a smaller subdivision in both f and g (C).
% Contributions from D and E are discarded, as they have already been counted
% above. Complexity O(m^2 + n^2)
%   ________________
%   : /E|C/:      /
%   :   /  :    /
%   : /D|B :  /
%   /___|__:/
%  a+c fl a+d     b+d
%
% Rather than make a BNDFUN corresponding to the each of the patches, we
% instead evaluate directly on a corresonding Chebyshev grid of appropriate
% size, which turns out to be far more efficient.
%
% Total complexity: O( R*(m*n + n*n) + m*m )
%
% However, in the case when R is too big we revert to a Clenshaw-Curtis
% quadrature-based approach to compute the convolution in the inner rectangle.
% This approach has a complexity O( (m + n)^3 ). We naively compare this with
% the big-O term above to determine which approach to use.

% [TODO]: It's possible this should be pushed further to the chebtech level.
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Return empty for an empty input:
if ( isempty(f) || isempty(g) )
    h = bndfun();
    return
end

% Extract the domain:
domF = f.domain;
a = domF(1);
b = domF(2);

domG = g.domain;
c = domG(1);
d = domG(2);

% Ensure that g is the signal (i.e., on the larger domain) and f is the filter:
if ( (b - a) > (d - c) )
    h = conv(g, f);
    return
end

% Useful things:
M = length(f);                               % Length of F
N = length(g);                               % Length of g
numPatches = floor((d - c) / (b - a));       % Number of patches required
x = chebpts(N, [b+c, a+d], 1);               % Chebyshev grid for interior piece
y = 0*x;                                     % Initialise values in interior
map = @(x, a, b) (x-a)/(b-a) - (b-x)/(b-a);  % Map from [a, b] --> [-1, 1]

% If there are too many patches then the HT approach is too slow. In such a case
% we resort to the standard quadrature-based approach (but still use the HT
% approach for the two triangular domains at the ends).
coeffsConvCost = numPatches*((M+N)*N);
quadConvCost = (M+N)^3;
if ( numPatches > 1 && coeffsConvCost > quadConvCost )
    h_left = conv(f, restrict(g, c+[0, b-a]));     % Left triangle
    h_right = conv(f, restrict(g, d-[(b-a) 0]));   % Right tirangle
    % Middle:
    [t, w] = legpts((M+N+5)/2, [a,b]);             % Legendre grid
    [tt, xx] = meshgrid(t,x);                      % Cheb/Leg grid to evaluate g
    ft = feval(f,t);                               % Evaluate f
    gxmt = feval(g, xx-tt);                        %    and g (expensive)
    y = gxmt*(w'.*ft);                             % Compute integral
    y = chebtech1.vals2coeffs(y);                  % Convert values to coeffs
    % Trim small coefficients:
    ay = abs(y); my = max(ay); loc = max(find(ay > 10*eps*my, 1, 'last'),1);
    if ( isempty(loc) ), loc = 1; end              % Deal with case when y = 0
    y = y(1:loc);
    data.domain = [b+c, a+d];
    h_mid = bndfun({[],y}, data);                  % Make h_mid from coeffs
    h = {h_left{1}, h_mid, h_right{end}};          % Combine three pieces
    h = fixMaps(h);                                % Ensure domain ends match
    return
end

% Trim small coefficients in f:
f_cheb = get(f, 'coeffs'); af = abs(f_cheb); mf = max(af);
loc = find(af > 10*eps*mf, 1, 'last'); if ( isempty(loc) ), loc = 1; end
f_cheb = f_cheb(1:loc);
f_leg = cheb2leg(f_cheb);                   % Legendre coefficients of f

% Restrict g:
doms = c + (b - a)*(0:numPatches);
g_restricted = restrict(g, doms);
doms = c + (b - a)*(0:numPatches);           % Subdomains
if ( ~iscell(g_restricted) )
    % If doms happened to be domain(g), restrict would return a cell.
    g_restricted = {g_restricted};
end

% Loop over the patches:
for k = 1:numPatches
                          %         _____
    dk = doms([k, k+1]);  %       /|    /
    dk_left  = a + dk(1); %     /  |  /
    dk_mid   = a + dk(2); %   /____|/
    dk_right = b + dk(2); %  dkl  dkm   dkr
    gk = g_restricted{k};                          % g on this subdomain
    gk = simplify(gk);                             % Simplify for efficiency
    gk_leg = cheb2leg(get(gk, 'coeffs'));          % Legendre coefficients
    [hLegL, hLegR] = easyConv(f_leg, gk_leg);      % Convolution on this domain

    % The left triangle for the kth patch:
    ind = (dk_left <= x) & (x < dk_mid); % Locate the grid values in [dkl, dkr]:
    if ( k == 1 ) % First piece:
        hLegL = leg2cheb(hLegL);                   % Cheb. coeffs of left tri.
        data.domain = [dk_left, dk_mid];
        h_left = bndfun({[], hLegL}, data);        % Make BNDFUN from coeffs
    else          % Subsequent left pieces
        z = map(x(ind), dk_left, dk_mid);          % Map grid points to [-1, 1]
        tmp = clenshawLegendre(z, hLegL);          % Evaluate via recurrence
        y(ind) = y(ind) + tmp;                     % Append
    end

    % The right triangle for the kth patch:
    if ( k < numPatches )                          % Not needed for final patch!
        % Locate the grid values in [dkl, dkr]:
        ind = (dk_mid <= x) & (x < dk_right);
        z = map(x(ind), dk_mid, dk_right);
        tmp = clenshawLegendre(z, hLegR);
        y(ind) = y(ind) + tmp;
    end

end

if ( abs((b-a)-(d-c)) < 10*eps(norm([a b c d], inf)) )
    % If there's only one patch, then we already have all the information reqd.
    hLegR = leg2cheb(hLegR);                        % Cheb coeffs of right tri.
    data.domain = d + [a b];
    h_right = bndfun({[], hLegR}, data);            % Make BNDFUN from coeffs
    h_mid = bndfun();

else
    % Final right right triangle:
    %  ________________
    %  : /E|C/:      /
    %  :   /  : A  /
    %  : /D|B :  /
    %  /___|__:/
    %     fl a+d     b+d

    finishLocation = a + c + numPatches*(b - a);    % Where patches got to. (fl)
    gk = restrict(g, d-[(b-a) 0]);                  % g on appropriate domain
    gk = simplify(gk);                              % Simplify for efficiency
    gk_leg = cheb2leg(get(gk, 'coeffs'));           % Legendre coeffs
    [hLegL, hLegR] = easyConv(f_leg, gk_leg);       % Conv on A and B
    hLegR = leg2cheb(hLegR);                        % Cheb coeffs on A
    data.domain = [d+a, d+b];
    h_right = bndfun({[], hLegR}, data);            % Make BNDFUN from coeffs

    % Remainder piece: (between fl and a+d)
    remainderWidth = d + a - finishLocation; % b+d-fl-(b-a)
    if ( remainderWidth > 0 )
        ind = finishLocation <= x;           % Discard D and E

        % B: (Coeffs were computed above)
        z = map(x(ind), d - b + 2*a, d + a); % Map grid to [-1, 1]
        tmp = clenshawLegendre(z, hLegL);    % Evaluate via recurrence
        y(ind) = tmp;                        % Store

        % C:
        domfk = b + [-remainderWidth, 0];               % Domain of fk
        domfk(1) = max(domfk(1), f.domain(1));          % Ensure domfk is a
        domfk(end) = min(domfk(end), f.domain(end));    %  valid subdomain
        fk = restrict(f, domfk);                        % Restrict f
        fk = simplify(fk);                              % Simplify f
        fk_leg = cheb2leg(get(fk, 'coeffs'));           % Legendre coeffs
        domgk = [finishLocation, d + a] - b;            % Domain of gk
        domgk(1) = max(domgk(1), g.domain(1));          % Ensure domgk is a
        domgk(end) = min(domgk(end), g.domain(end));    %  valid subdomain
        gk = restrict(g, domgk);                        % Restrict g
        gk = simplify(gk);                              % Simplify g
        gk_leg = cheb2leg(get(gk, 'coeffs'));           % Legendre coeffs
        [~, hLegR] = easyConv(fk_leg, gk_leg);          % Conv
        z = map(x(ind), finishLocation, d + a);         % Map to [-1, 1]
        tmp = clenshawLegendre(z, hLegR);               % Eval via recurrence
        y(ind) = y(ind) + tmp*remainderWidth/(b - a);   % Scale and append
    end
    % Convert values to coeffs (we don't want to construct a chebtech1)
    y = chebtech1.vals2coeffs(y);

    % Construct BNDFUN of the interior (rectangle) using coefficients:
    data.domain = [b+c, a+d];
    h_mid = bndfun({[], y}, data);

end

% h_mid can be empty so return the three or two pieces as a cell array:
if ( isempty(h_mid) )
    h = {h_left*(b-a)/2, h_right*(b-a)/2};
else
    h = {h_left*(b-a)/2, h_mid*(b-a)/2, h_right*(b-a)/2};
end
h = fixMaps(h);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function h = fixMaps(h)
% Ensure endpoints of BNDFUNs for adjacent subintervals match exactly. NB: This
% is a bit inefficient if consecutive endpoints get adjusted, as some BNDFUNs
% may have their maps changed twice, but the map-change operation is fast enough
% that this shouldn't matter.
for n = 2:1:numel(h)
    end_left = h{n-1}.domain(end);
    end_right = h{n}.domain(1);
    if ( end_left ~= end_right )
        hs_left = norm(h{n-1}.domain, Inf);  % hscale for left piece.
        hs_right = norm(h{n}.domain, Inf);   % hscale for right piece.

        % If there's a mismatch, it should be small because the BNDFUNs should
        % come out in order corresponding to consecutive adjacent subintervals.
        if ( abs(end_left - end_right) < 2*eps*max(hs_left, hs_right) )
            new_end = (end_left + end_right)/2;
            h{n-1} = changeMap(h{n-1}, [h{n-1}.domain(1) new_end]);
            h{n} = changeMap(h{n}, [new_end h{n}.domain(end)]);
        else
            % We should only get here if a programmer error made elsewhere
            % causes the BNDFUNs to get out of order.
            error('CHEBFUN:BNDFUN:conv:nonConsecutive', ...
                'Pieces do not belong to consecutive subintervals.');
        end
    end
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [gammaL, gammaR] = easyConv(alpha, beta)
% Convolution using Legendre expansions and the analoguous convolution theorem.
% See Hale and Townsend, "An algorithm for the convolution of Legendre series",
% SIAM Journal on Scientific Computing, Vol. 36, No. 3, pages A1207-A1220, 2014.

% Better computational efficiency is achieved when g has the lower degree:
if ( length(beta) > length(alpha) )
    tmp = alpha;
    alpha = beta;
    beta = tmp;
end

% Maximum degree of result:
MN = length(alpha) + length(beta);

% Pad to make length n + 1.
alpha = [ alpha ; zeros(MN - length(alpha), 1) ];

% S represents multiplication by 1/z in spherical Bessel space:
e = [[1 ; 1./(2*(1:(MN-1)).'+1)], [1 ; zeros(MN-1, 1)], -1./(2*(0:MN-1).'+1)];
S = spdiags(e, -1:1, MN, MN);

gammaL = rec(S, alpha, beta, -1); % Chebyshev coeffs for the left piece
S(1,1) = -1;                      % Update S
gammaR = rec(S, -alpha, beta, 1); % Chebyshev coeffs for the right piece

%%%%%%%%%%%%%%%%%%%%%%%%%%% MATRIX FREE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function gamma = rec(S, alpha, beta, sgn)
        % Compute the Legendre coefficients of the convolution on L/R piece.
        % See Theorem 4.1 of paper.

        % Initialise scl:
        N = length(beta);
        scl = 1./(2*(1:N).'-1);
        scl(2:2:end) = -scl(2:2:end);

        % First column of B:
        vNew = S*alpha;
        v = vNew;
        gamma = beta(1)*vNew;
        beta_scl = scl.*beta;
        beta_scl(1) = 0;
        gamma(1) = gamma(1) + vNew(1:N).'*beta_scl;

        % The scalar case is trivial:
        if ( length(beta) == 1 )
            return
        end

        % Second column of B:
        vNew = S*v + sgn*v;
        vOld = v;
        v = vNew;
        vNew(1) = 0;

        gamma = gamma + beta(2)*vNew;
        beta_scl = -beta_scl*((2 - 0.5)/(2 - 1.5));
        beta_scl(2) = 0;
        gamma(2) = gamma(2) + vNew(1:N).'*beta_scl;

        % Loop over remaining columns using recurrence:
        for n = 3:N
            vNew = (2*n-3) * (S * v) + vOld; % Recurrence
            vNew(1:n-1) = 0;                 % Zero terms
            gamma = gamma + vNew*beta(n);    % Append to g

            % Recurrence is unstable for j < k. Correct for upper-tri part:
            beta_scl = -beta_scl*((n-.5)/(n-1.5));
            beta_scl(n) = 0;
            gamma(n) = gamma(n) + vNew(1:N).'*beta_scl;

            vOld = v;
            v = vNew;

        end

        ag = abs(gamma);
        mg = max(ag);
        loc = find(ag > eps*mg, 1, 'last');
        gamma = gamma(1:loc);

    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function val = clenshawLegendre(x, alpha)
% Evaluate a Legendre expansion with coefficient alpha at x.

n = length(alpha);
b_old = 0;
b_cur = 0;
for k = (n-1):-1:1
  b_new = alpha(k+1) + (2*k + 1)/(k + 1)*x.*b_cur - (k + 1)/(k + 2)*b_old;
  b_old = b_cur;
  b_cur = b_new;
end
val = alpha(1) + x.*b_cur - .5*b_old;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function g = mlpOutputGrad(model, X)

% MLPOUTPUTGRAD Evaluate derivatives of mlp model outputs with respect to parameters.
% FORMAT
% DESC evaluates the derivates of a multi-layer perceptron's
% outputs with respect to the parameters of the multi-layer
% perceptron. Currently it simply wraps the NETLAB mlpderiv
% function.
% ARG model : the model for which the derivatives are to be
% computed.
% ARG X : the input data locations where the gradients are to be
% computed.
% RETURN g : the gradient of the outputs of the multi-layer
% perceptron with respect to each of the parameters. The size of
% the matrix is number of data x number of parameters x number of
% outputs of the model.
%
% SEEALSO : mlpCreate, mlpderiv
%
% COPYRIGHT : Neil D. Lawrence, 2006

% MLTOOLS

if length(model.hiddenDim) == 1
  g = mlpderiv(model, X);
else

  [Y, G, A] = mlpOut(model, X);
  gw = cell(1, length(model.w));
  for i = 1:length(gw)
    gw{i} = zeros(size(model.w{i}));
  end
  for i = 1:length(G);
    WdG{i} = (1-G{i}.*G{i})*w{i};
  end
  for k = 1:model.outputDim
    gw{end}(:, k) = Z{end}(:, k);
    gb{end} = 1;
    for i = length(model.w)-1:-1:1
      %gw{i} =
      error('Not yet implemented');
    end
  end

% Evaluate second-layer gradients.
gw2 = z'*deltas;
gb2 = sum(deltas, 1);

% Now do the backpropagation.
delhid = deltas*net.w2';
delhid = delhid.*(1.0 - z.*z);

% Finally, evaluate the first-layer gradients.
gw1 = x'*delhid;
gb1 = sum(delhid, 1);

g = [gw1(:)', gb1, gw2(:)', gb2];
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [afixed,sqnorm,ierr] = lsearch (afloat,L,D,Chi2,ncands)
%LSEARCH: Integer ambiguity resolution, search
%
% This routine finds the integer vector which is closest to a given
% float vector, in a least squares sence. This is the search-step in
% integer ambiguity resolution. It is best to perform this search only
% on ambiguities which have been decorrelated using LAMBDA.
%
% Input arguments:
%    afloat : Float ambiguities (\hat{a})
%    L      : LtDL-decomposition of the decorrelated
%    D        variance-covariance matrix of the ambiguities
%    Chi2   : Size of the search ellipsoid
%    ncands : Number of requested candidates
%
% Output arguments:
%    afixed : Estimated integers (matrix)
%    sqnorm : Corresponding squared norms (vector, sorted)
%    ierr   : Error code: 0: No errors found
%                         1: Not enough candidates found

% ----------------------------------------------------------------------
% File.....: lsearch.m
% Date.....: 19-MAY-1999
% Author...: Peter Joosten
%            Mathematical Geodesy and Positioning
%            Delft University of Technology
% ----------------------------------------------------------------------

% -------------------------------
% --- Initializing statements ---
% -------------------------------

Linv      = (L)^(-1);
Dinv      = 1./D;

True      = 1;
False     = 0;

n         = max(size(afloat));

right     = [zeros(n,1) ; Chi2];
left      = [zeros(n+1,1)];
dq        = [Dinv(2:n)./Dinv(1:n-1) 1/Dinv(n)];

cand_n    = False;
c_stop    = False;
endsearch = False;

ncan      = 0;

i         = n + 1;
iold      = i;
ierr      = 0;

afixed = zeros(n,ncands);
sqnorm = zeros(1,ncands);

% ----------------------------------
% --- Start the main search-loop ---
% ----------------------------------

while ~ (endsearch);

   i = i - 1;

   if iold <= i
      lef(i) = lef(i) + Linv(i+1,i);
   else
      lef(i) = 0;
      for j = i+1:n;
         lef(i) = lef(i) + Linv(j,i)*distl(j,1);
      end;
   end;
   iold = i;

   right(i)   = (right(i+1) - left(i+1)) * dq(i);
   reach      = sqrt(right(i));
   delta      = afloat(i) - reach - lef(i);
   distl(i,1) = ceil(delta) - afloat(i);

   if distl(i,1) > reach - lef(i)

%     ----------------------------------------------------
%     --- There is nothing at this level, so backtrack ---
%     ----------------------------------------------------

      cand_n = False;
      c_stop = False;

      while (~ c_stop) & (i < n);

         i = i + 1;
         if distl(i) < endd(i);
            distl(i) = distl(i) + 1;
            left(i)  = (distl(i) + lef(i)) ^ 2;
            c_stop   = True;
            if i == n; cand_n = True; end;
         end;

      end;

      if (i == n) & (~ cand_n); endsearch = True; end;

   else

%     ----------------------------
%     --- Set the right border ---
%     ----------------------------

      endd(i) = reach - lef(i) - 1;
      left(i) = (distl(i,1) + lef(i)) ^ 2;

   end

   if i == 1;

%     -------------------------------------------------------------------
%     --- Collect the integer vectors and corresponding               ---
%     --- squared distances, add to vectors "afixed" and "sqnorm" if: ---                             ---
%     --- * Less then "ncands" candidates found so far                ---
%     --- * The squared norm is smaller than one of the previous ones ---
%     -------------------------------------------------------------------

      t       = Chi2 - (right(1)-left(1)) * Dinv(1);
      endd(1) = endd(1) + 1;

      while distl(1) <= endd(1);

         if ncan < ncands;

            ncan             = ncan + 1;
            afixed(1:n,ncan) = distl + afloat;
            sqnorm(ncan)     = t;

         else

            [maxnorm,ipos] = max(sqnorm);
            if t < maxnorm;
               afixed(1:n,ipos) = distl + afloat;
               sqnorm(ipos)     = t;
            end;

         end;

         t       = t + (2 * (distl(1) + lef(1)) + 1) * Dinv(1);
         distl(1) = distl(1) + 1;

      end;

%     -------------------------
%     --- And backtrack ... ---
%     -------------------------

      cand_n = False;
      c_stop = False;

      while (~ c_stop) & (i < n);

         i = i + 1;

         if distl(i) < endd(i);
            distl(i) = distl(i) + 1;
            left(i)  = (distl(i) + lef(i)) ^ 2;
            c_stop   = True;
            if i == n; cand_n = True; end;
         end;

      end;

      if (i == n) & (~ cand_n); endsearch = True; end;

   end;

end;

% ----------------------------------------------------------------------
% --- Sort the resulting candidates, according to the norm
% ----------------------------------------------------------------------

tmp    = sortrows ([sqnorm' afixed']);
sqnorm = tmp(:,1)';
afixed = tmp(:,2:n+1)';

% ------------------------
% --- Check for errors ---
% ------------------------

if ncan < ncands; ierr = 1; end;

% ----------------------------------------------------------------------
% End of routine: lsearch
% ----------------------------------------------------------------------

author: goGPS-Project


text:
function x=v_usasi(n,fs)
%V_USASI generates N samples of USASI noise at sample frequency FS X=(N,FS)

% This routine is based on the USASI noise defined in [1] which was later
% reissued as [2]. USASI noise is intended to simulate the long-term average
% of typical audio program material. The routine does not currently implement
% the pulsation at 2.5Hz 12.5% duty cycle that is recommended by the standard.
% Also it should probably be scaled to a well-defined power.
%
%  [1] NRSC AM Reemphasis, Deemphasize, and Broadcast Audio Transmission Bandwidth Specifications,
%      EIA-549 Standard, Electronics Industries Association , July 1988.
%  [2] NRSC AM Reemphasis, Deemphasize, and Broadcast Audio Transmission Bandwidth Specifications,
%      NRSC-1-A Standard, Sept 2007, Online: http://www.nrscstandards.org/SG/NRSC-1-A.pdf

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: v_usasi.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<2 fs=8000; end
b=[1 0 -1];
a=poly(exp(-[100 320]*2*pi/fs));

x=v_randfilt(b,a,n);

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [varargout]=cmaperise(varargin)

% function [Cmapped,indMap]=cmaperise(C,cmap,clim)
%-------------------------------------------------------------------------
% This function creates RGB colors for the data C using the colormap cmap
% and the limits clim.
%
%
% Change log:
% 2019/06/25 Added variable input handling
%-------------------------------------------------------------------------

%% Parse input

switch nargin
    case 1
        C=varargin{1};
        cmap=[];
        clim=[];
    case 2
        C=varargin{1};
        cmap=varargin{2};
        clim=[];
    case 3
        C=varargin{1};
        cmap=varargin{2};
        clim=varargin{3};
end

if isempty(cmap)
    cmap=viridis(250);
end

if isempty(clim)
   clim=[min(C(:)) max(C(:))];
   if diff(clim)<eps(0)
       clim=[min(C(:))-1e-3 max(C(:))+1e-3];
   end
end

%% Use colormap to define RGB values

C=C(:); %Color data as column

p=(C-clim(1))./(clim(2)-clim(1));
p(p<0)=0;
p(p>1)=1;

p(isnan(p))=0; %Force nan data to 0

IND_cmap=round((p*(size(cmap,1)-1))+1);
Cmapped=cmap(IND_cmap,:);

%%
varargout{1}=Cmapped;
if nargout==2
    varargout{2}=IND_cmap;
end

end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function doublePendulumAnimate(sol,P)

%This function is used to run an animation of the double pendulum:

duration = sol.x(end) - sol.x(1);

tic;  %Start a timer
timeNow = 0;

h = gcf;

len = P.l1+P.l2;
axisVec = [-len,len,-len,len];

s1 = sqrt(P.m1);
s2 = sqrt(P.m2);
sMin = min(s1,s2);
s1 = s1/sMin;
s2 = s2/sMin;

while timeNow < duration

    zNow = deval(sol,timeNow);
    [p1,p2,g1,g2] = doublePendulumPosition(zNow,P);

    figure(h);
    plotFrame(timeNow,p1,p2,g1,g2,s1,s2);
    axis(axisVec); axis square; axis off;
    drawnow;

    pause(0.001);
    timeNow = toc;
end

end

function plotFrame(time,p1,p2,g1,g2,s1,s2)

linkOneX = [0;p1(1)];
linkOneY = [0;p1(2)];

linkTwoX = [p1(1);p2(1)];
linkTwoY = [p1(2);p2(2)];

clf;
hold on;

plot(0,0,'ks','MarkerSize',20,'LineWidth',3);

plot(linkOneX,linkOneY,'r-','LineWidth',6);
plot(linkTwoX,linkTwoY,'b-','LineWidth',6);

plot(0,0,'k.','MarkerSize',25);
plot(p1(1),p1(2),'k.','MarkerSize',25);
plot(p2(1),p2(2),'k.','MarkerSize',25);

plot(g1(1),g1(2),'ko','MarkerSize',round(15*s1),'LineWidth', 3);
plot(g2(1),g2(2),'ko','MarkerSize',round(15*s2),'LineWidth', 3);
plot(g1(1),g1(2),'kx','MarkerSize',round(15*s1),'LineWidth', 3);
plot(g2(1),g2(2),'kx','MarkerSize',round(15*s2),'LineWidth', 3);

title(sprintf('Time: %4.3f',time));

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function DEM_demo_OU
% DEM demo for linear deconvolution:  This demo considers the deconvolution
% of one of the simplest dynamical process; a random walk or Ornstein-
% Uhlenbeck process.  It shows how DEM can infer on the causes as stochastic
% innovations (c.f., Bayesian filtering) by exploiting temporal
% correlations.  Strictly speaking this is not a Ornstein-Uhlenbeck process
% because the innovations are themselves correlated and would normally be a
% Wiener process

% get a simple convolution model
%==========================================================================
M       = spm_DEM_M('OU');

% and generate data
%==========================================================================
N       = 64;                                 % length of data sequence
DEM     = spm_DEM_generate(M,N,{},{[] 8});

% invert model
%==========================================================================
DEM     = spm_DEM(DEM);

% overlay true values
%--------------------------------------------------------------------------
spm_DEM_qU(DEM.qU,DEM.pU)

author: spm


text:
% StackExchange Mathematics Q1909139
% https://math.stackexchange.com/questions/1909139
% Projection of a Symmetric Matrix onto the Probability Simplex
% References:
%   1.
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     17/06/2021
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

numRows = 4;

%% Generate Data

mY = randn(numRows, numRows);
% mY = mY + mY.';

vI = zeros(numRows * numRows, 1);
vI(1:(numRows + 1):(numRows ^ 2)) = 1;

hObjFun = @(mX) 0.5 * sum((mX(:) - mY(:)) .^ 2);

%% Solution by CVX

solverString = 'CVX';

% cvx_solver('SDPT3'); %<! Default, Keep numRows low
% cvx_solver('SeDuMi');
% cvx_solver('Mosek'); %<! Can handle numRows > 500, Very Good!
% cvx_solver('Gurobi');

hRunTime = tic();

cvx_begin('quiet')
    % cvx_precision('best');
    variable mX(numRows, numRows) symmetric;
    minimize( (0.5 * sum_square(mX(:) - mY(:))) );
    subject to
        % mX <In> symmetric(numRows); %<! Doesn't work for some reason
        trace(mX) == 1;
        mX >= 0;
cvx_end

runTime = toc(hRunTime);

DisplayRunSummary(solverString, hObjFun, mX, runTime, cvx_status);

sCvxSol.vXCvx     = mX(:);
sCvxSol.cvxOptVal = hObjFun(mX);

%% Solution by Orthogonal Projection onto the Intersection of Convex Sets

hP1 = @(vX) reshape((reshape(vX, numRows, numRows) + reshape(vX, numRows, numRows).') / 2, 1, numRows * numRows);
hP2 = @(vX) max(vX, 0); %<! Projection mX >= 0
hP3 = @(vX) (vX .* ~vI) + ((vX - ((sum(vX .* vI) - 1) / numRows)) .* vI); %<! Projection trace(mX) == 1

hRunTime = tic();
vX = OrthogonalProjectionOntoConvexSets({hP1; hP2; hP3}, mY(:), 100, 1e-6);
runTime = toc(hRunTime);

mX1 = reshape(vX, numRows, numRows);

DisplayRunSummary(solverString, hObjFun, reshape(vX, numRows, numRows), runTime);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
% Generate samples from the HHMM with the true params.

seed = 1;
rand('state', seed);
randn('state', seed);

discrete_obs = 0;

bnet = mk_square_hhmm(discrete_obs, 1);
Q1 = 1; Q2 = 2; Q3 = 3; F3 = 4; F2 = 5; Onode = 6;
Qnodes = [Q1 Q2 Q3]; Fnodes = [F2 F3];

for seqi=1:1
  evidence = sample_dbn(bnet, 'stop_test', 'is_F2_true_D3');
  clf
  plot_square_hhmm(evidence);
  %pretty_print_hhmm_parse(evidence, Qnodes, Fnodes, Onode, []);
  fprintf('sequence %d has length %d; press key to continue\n', seqi, size(evidence,2))
  pause
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [Psi, logQ] = RestrictedGibbsScan_Merge( Psi, data_struct, hyperparams, model, TargetPsi )

featIDs = Psi.activeFeatIDs;
objIDs = find( sum( Psi.F(:, featIDs ) , 2 ) > 0 )';

logQ = struct( 'z', 0, 'TS', 0,'theta', 0, 'F', 0, 'all', 0 );

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

[Psi.stateSeq, logQ.z] = sampleStateSeq_RGS( Psi.F, Psi.stateSeq, Psi.TS, Psi.theta, data_struct, model, objIDs, TargetPsi );
[Psi.TS, logQ.TS] = sampleTransParams_RGS( Psi.F, Psi.stateSeq, Psi.TS, hyperparams, model, objIDs, TargetPsi);
[Psi.theta, logQ.theta] = sampleTheta_RGS( Psi.F, Psi.stateSeq, Psi.theta, data_struct, model, featIDs , TargetPsi);

logQ.all = logQ.F + logQ.theta + logQ.TS + logQ.z;

author: michaelchughes

repo: NPBayesHMM


text:
function [ loadFiles, loadDirs ] = findFilesDirs( imagePath)
% recursive search
%
% (c) Marc Fischer, Thomas Kuestner
% ---------------------------------------------------------------------

[filesInPath, dirsInPath] = rdir(imagePath);
% filter files (just use *.mat and *.dat)
allFiles = cell(1,length(filesInPath));
allDirs = cell(1,length(dirsInPath));
foundFile = false(1,length(filesInPath));
if ~isempty(dirsInPath)
    nDirsInPath = length(dirsInPath);
    for k = 1:nDirsInPath
        allDirs{k} = dirsInPath.fullpath;
    end;
else
    nDirsInPath = 1;
end;
foundDir = false(1,length(nDirsInPath));

for i=1:length(filesInPath)
    allFiles{i} = filesInPath(i).fullpath;
        tmp = regexp(allFiles{i},'\w+x\d*\w+\.[md]at', 'once');
        tmp2 = regexp(allFiles{i},'\w+_rec\w+', 'once');
        if (~isempty(tmp) && isempty(tmp2))
            foundFile(i) = true;
        end;
end;
for k=1:nDirsInPath; % for DICOM only check for dirs, since fReadDICOM searches a whole folder -> if folder contains *.IMA-file, set found(k) true
    for i=1:length(filesInPath)
        if ~isempty(allDirs)
            tmp = regexp(allFiles{i},'\w+.IMA', 'once');
            if(~isempty(tmp))
                dirString = regexptranslate('escape', allDirs{k});
                tmp_dir = regexp(allFiles{i},[dirString '\w*'],'once');
                if(~isempty(tmp_dir))
                    foundDir(k) = true;
                end;
            end;
        else
            allDirs{k} = imagePath;
            allFiles{i} = filesInPath(i).fullpath;
            tmp = regexp(allFiles{i},'\w+.IMA', 'once');
            if(~isempty(tmp))
                dirString = regexptranslate('escape', allDirs{k});
                regexp(allFiles{i},[dirString '\w*']);
                foundDir(k) = true;
            end;
        end;
    end;
end

loadFiles = allFiles(foundFile);
loadDirs = allDirs(foundDir);

% remove same file names (only needed for mat/dat)
helper = cellfun(@(x) x(1:end-4), loadFiles, 'UniformOutput', false); % removes file extension in string
[~, idxloadFiles, ~] = unique(helper); % sorts out double entries
loadFiles = loadFiles(idxloadFiles);

end

author: thomaskuestner


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Setting up advection-diffusion solver
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [C,laplacian_C] = please_Update_Adv_Diff_Concentration_Flux_Limiter_FV(C,dt,dx,dy,uX,uY,k)

% C:     concentration
% dt:    time-step
% dx,dy: spatial steps in x and y, respectively
% uX:    x-Component of Velocity
% uY:    y-Component of Velocity
% k:     diffusion coefficient

% Compute Fluxes (Note: these calculations could be parallalized)
selection = 'superbee';
Fx = give_Necessary_Fluxes(C,dx,uX,'x',selection,dt); % Fluxes in x
Fy = give_Necessary_Fluxes(C,dy,uY,'y',selection,dt); % Fluxes in y

% "forward difference of fluxes in x"
F2x = [Fx(:,2:end) Fx(:,1)];
F1x = [Fx(:,end) Fx(:,1:end-1)];
diffX = 0.5/dx*( F2x - F1x );

% "forward differences of fluxes in y"
F2y = [Fy(2:end,:); Fy(1,:)];
F1y = [Fy(end,:); Fy(1:end-1,:)];
diffY = 0.5/dy*( F2y - F1y );

diffY=0;

% Compute 2nd Derivative Terms
Cxx = DD(C,dx,'x');
Cyy = DD(C,dy,'y');

% Forms Laplacian
laplacian_C = Cxx+Cyy;

% UPWIND
C = C - dt * ( diffX + diffY ) + dt*( k*laplacian_C );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Computes derivative based on sign of Velocity, u, using UPWIND
% approach
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C_z = give_Necessary_Fluxes(C,dz,uZ,string,selection,dt)

C_z = zeros(size(C));
len = length(uZ(:,1));
signs = sign(uZ);

if strcmp(string,'x')

    %For periodicity on ends w/ UPWIND
    for i=1:len

        %left side of grid
        if signs(i,1) <= 0
            r = ( C(i,2) - C(i,1) ) / ( C(i,1) - C(i,len) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(i,1) =  uZ(i,1)*C(i,2) ;%  + 0.5*abs( uZ(i,1) )*( 1 - abs( 0.5*dt*uZ(i,1)/dz ) )*phi*( C(i,2) - C(i,len) );
        else
            r = ( C(i,1) - C(i,len) ) / ( C(i,2) - C(i,1) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(i,1) =  uZ(i,1)*C(i,len) ;%+ 0.5*abs( uZ(i,1) )*( 1 - abs( 0.5*dt*uZ(i,1)/dz ) )*phi*( C(i,2) - C(i,len) );
        end

        %right side of grid
        if signs(len,1) <= 0
            r = ( C(i,1) - C(i,len) ) / ( C(i,len) - C(i,len-1) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(i,len) =  uZ(i,len)*C(i,1)  ;%   + 0.5*abs( uZ(i,len) )*( 1 - 0.5*abs( dt*uZ(i,len)/dz ) )*phi*( C(i,1) - C(i,len-1) );
        else
            r = ( C(i,len) - C(i,len-1) ) / ( C(i,1) - C(i,len) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(i,len) =  uZ(i,len)*C(i,len-1);% + 0.5*abs( uZ(i,len) )*( 1 - 0.5*abs( dt*uZ(i,len)/dz ) )*phi*( C(i,1) - C(i,len-1) );
        end

    end
    %Standard Upwind
    for i=1:len
        for j=2:len-1
            if signs(i,j) <= 0
                r = ( C(i,j+1) - C(i,j) ) / ( C(i,j) - C(i,j-1) );
                phi = please_Give_Flux_Limiter(r,selection);
                C_z(i,j) = uZ(i,j)*C(i,j+1);% + 0.5*abs( uZ(i,j) )*( 1 - 0.5*abs( dt*uZ(i,j)/dz ) )*phi*( C(i,j+1) - C(i,j-1) );
            else
                r = ( C(i,j) - C(i,j-1) ) / ( C(i,j+1) - C(i,j) );
                phi = please_Give_Flux_Limiter(r,selection);
                C_z(i,j) = uZ(i,j)*C(i,j-1);% + 0.5*abs( uZ(i,j) )*( 1 - 0.5*abs( dt*uZ(i,j)/dz ) )*phi*( C(i,j+1) - C(i,j-1) );
            end
        end
    end

    % Ends x-Direction calculation %

elseif strcmp(string,'y')

    %For periodicity on ends w/ UPWIND
    for i=1:len

        %bottom of grid
        if signs(1,i) <= 0
            r = ( C(2,i) - C(1,i) ) / ( C(1,i) - C(len,i) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(1,i) =  uZ(1,i)*C(2,i)   + 0.5*abs( uZ(1,i) )*( 1 - abs( 0.5*dt*uZ(1,i)/dz ) )*phi*( C(2,i) - C(len,i) );
        else
            r = ( C(1,i) - C(len,i) ) / ( C(2,i) - C(1,i) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(1,i) =  uZ(1,i)*C(len,i) + 0.5*abs( uZ(1,i) )*( 1 - abs( 0.5*dt*uZ(1,i)/dz ) )*phi*( C(2,i) - C(len,i) );
        end

        %top of grid
        if signs(len,1) <= 0
            r = ( C(1,i) - C(len,i) ) / ( C(len,i) - C(len-1,i) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(len,i) =  uZ(len,i)*C(1,i)       + 0.5*abs( uZ(len,i) )*( 1 - abs( 0.5*dt*uZ(len,i)/dz ) )*phi*( C(1,i) - C(len-1,i) );
        else
            r = ( C(len,i) - C(len-1,i) ) / ( C(1,i) - C(len,i) );
            phi = please_Give_Flux_Limiter(r,selection);
            C_z(len,i) =  uZ(len,i)*C(len-1,i)   + 0.5*abs( uZ(len,i) )*( 1 - abs( 0.5*dt*uZ(len,i)/dz ) )*phi*( C(1,i) - C(len-1,i) );
        end

    end

    %Standard Upwind
    for i=1:len
        for j=2:len-1
            if signs(j,i) <= 0
                r = ( C(j+1,i) - C(j,i) ) / ( C(j,i) - C(j-1,i) );
                phi = please_Give_Flux_Limiter(r,selection);
                C_z(j,i) = uZ(j,i)*C(j+1,i) + 0.5*abs( uZ(j,i) )*( 1 - abs( 0.5*dt*uZ(j,i)/dz ) )*phi*( C(j+1,i) - C(j-1,i) );
            else
                r = ( C(j,i) - C(j-1,i) ) / ( C(j+1,i) - C(j,i) );
                phi = please_Give_Flux_Limiter(r,selection);
                C_z(j,i) = uZ(j,i)*C(j-1,i) + 0.5*abs( uZ(j,i) )*( 1 - abs( 0.5*dt*uZ(j,i)/dz ) )*phi*( C(j+1,i) - C(j-1,i) );
            end
        end
    end

    % Ends y-Direction calculation %

else

    fprintf('\n\n\n ERROR IN FUNCTION FOR COMPUTING FLUX LIMITERS\n');

end

clear signs; clear len;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Finds CENTERED finite difference approximation to 2ND
% DERIVATIVE in z direction, specified by input and 'string'
% Note: It automatically accounts for periodicity of the domain.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u_zz = DD(u,dz,string)

% u:      velocity
% dz:     spatial step in "z"-direction
% string: specifies which 2ND derivative to take (to enforce periodicity)

len = length(u(:,1));

if strcmp(string,'x')

    %For periodicity on ends
    u_zz(:,1) =  ( u(:,2) - 2*u(:,1)   + u(:,len) )   / (dz^2);
    u_zz(:,len)= ( u(:,1) - 2*u(:,len) + u(:,len-1) ) / (dz^2);

    %Standard Upwind Scheme (Centered Difference)
    for j=2:len-1
        u_zz(:,j) = ( u(:,j+1) - 2*u(:,j) + u(:,j-1) ) / (dz^2);
    end

elseif strcmp(string,'y')

    %For periodicity on ends
    u_zz(1,:) =  ( u(2,:) - 2*u(1,:)   + u(len,:) )   / (dz^2);
    u_zz(len,:)= ( u(1,:) - 2*u(len,:) + u(len-1,:) ) / (dz^2);

    %Standard Upwind Scheme (Centered Difference)
    for j=2:len-1
        u_zz(j,:) = ( u(j+1,:) - 2*u(j,:) + u(j-1,:) ) / (dz^2);
    end

else

    fprintf('\n\n\n ERROR IN FUNCTION DD FOR COMPUTING 2ND DERIVATIVE\n');
    fprintf('Need to specify which desired derivative, x or y.\n\n\n');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes flux limiter with choice of which one
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function phi = please_Give_Flux_Limiter(r,selection)

if strcmp(selection,'superbee')
    max1 = max( min(1,2*r),min(2,r) );
    phi = max(0,max1);
elseif strcmp(selection,'vanLeer')
    phi = ( r + abs(r) ) / ( 1 + abs(r) );
else
    fprintf('\n\n');
    error('NEED TO CHOOSE AN APPROPRIATE FLUX LIMITER');
end

author: nickabattista


text:
function [imgCrop, scaling, seg, warn, valid] = LMobjectnormalizedcrop(img, annotation, j, b, height, width)
%
% Crop object from image using a normalized frame
% [imgCrop, scaling] = LMobjectnormalizedcrop(img, annotation, j, [bh bw], height, width)
% [imgCrop, scaling] = LMobjectnormalizedcrop(img, [xmin ymin xmax ymax], j, [bh bw], height, width)
%  extract object index j from the annotation.
%
%   ------------
%   |    bh    |   b = [bh bw] = number of boundary pixels
%   |   ----   |   h = height inner bounding box
%   |bw |  |   |   w = width inner bounding box
%   |   |h |   |
%   |   ----   |
%   |    w     |
%   ------------
%

if length(b) == 1
    b = b *[1 1];
end
bh = b(1);
bw = b(2);

if isstruct(annotation)
    [X,Y,t] = getLMpolygon(annotation.object(j).polygon);
else
    xmin = annotation(1); % [xmin ymin xmax ymax]
    ymin = annotation(2); % [xmin ymin xmax ymax]
    xmax = annotation(3); % [xmin ymin xmax ymax]
    ymax = annotation(4); % [xmin ymin xmax ymax]
    X = [xmin xmin xmax xmax];
    Y = [ymin ymax ymax ymin];
    t = 1;
end
numFrames = length(t);

if nargin < 6
    % decide which dimension is larger and use
    Dx = max(X)-min(X);
    Dy = max(Y)-min(Y);
    if Dy>Dx
        width = round(height*Dx/Dy);
    else
        % if Dx>Dy,
        h = round(height*Dy/Dx);
        width = height;
        height = h;
    end
end

imgCrop = zeros([height+2*bh width+2*bw size(img,3) numFrames], 'uint8');
seg = zeros([height+2*bh width+2*bw numFrames], 'uint8');
valid = zeros([height+2*bh width+2*bw numFrames], 'uint8');
warn = zeros(1, numFrames);
for f = t
    if numFrames>1
        i = find(t==f); i = i(1);
        x = X(:,i);
        y = Y(:,i);
    else
        i = 1;
        x = X;
        y = Y;
    end

    %  bb = boundingbox = [xmin ymin xmax ymax]
    bb = [min(x) min(y) max(x) max(y)];

    % 1) Resize image so that object is normalized in size
    scaling = min(height/(bb(4)-bb(2)), width/(bb(3)-bb(1)));
    if numFrames>1
        [foo, I] = LMimscale([], img(:,:,:,f+1), scaling, 'bilinear');
    else
        [foo, I] = LMimscale([], img, scaling, 'bilinear');
    end
    bb = bb*scaling; x = x*scaling; y = y*scaling;

    % 2) pad image (just to make sure)
    margin = 5+max(width,height)+ceil(max(bh,bw)+1 + max([0 -bb(1) -bb(2) bb(4)-size(I,1) bb(3)-size(I,2)]));
    [foo, I] = LMimpad(foo, single(I), [size(I,1)+2*margin size(I,2)+2*margin], NaN);
    bb = bb+margin; x = x+margin; y = y+margin;

    % 2) Crop result
    cx = fix((bb(3)+bb(1))/2); cy = fix((bb(4)+bb(2))/2);

    bb = round([cx-width/2-bw+1 cy-height/2-bh+1 cx+width/2+bw cy+height/2+bh]);

    x = x - bb(1)+2; y = y - bb(2)+2;

    % Image crop:
    crop = I(bb(2):bb(4), bb(1):bb(3), :);
    imgCrop(:,:,:,i) = uint8(crop);

    [xx,yy] = meshgrid(1:size(imgCrop,2),1:size(imgCrop,1));
    seg(:,:,i) = uint8(255*double(inpolygon(xx, yy, x, y)));

    valid(:,:,i) = 1-isnan(sum(crop,3));

    % warnings by crops of scales
    warn(i) = (scaling>1) || (sum(isnan(crop(:)))>500);
end

author: CSAILVision

repo: LabelMeToolbox


text:
% demos for ch06

%% Kernel regression with gaussian kernel
clear; close all;
n = 100;
x = linspace(0,2*pi,n);   % test data
t = sin(x)+rand(1,n)/2;
model = knReg(x,t,1e-4,@knGauss);
[y,s] = knRegPred(model,x);
plotCurveBar(x,y,s);
hold on;
plot(x,t,'o');
hold off;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function [F,V,faceBoundaryMarker]=quadBox(boxDim,boxEl)

%%
F=[]; V=[]; faceBoundaryMarker=[];

[Xtb,Ytb]=meshgrid(0:boxEl(1),0:boxEl(2));

[Ylr,Zlr]=meshgrid(0:boxEl(2),0:boxEl(3));
[f,v]=surf2patch(0*ones(size(Ylr)),Ylr,Zlr); %Left
f=fliplr(f);
F=[F;f+size(V,1)]; V=[V;v;]; faceBoundaryMarker=[faceBoundaryMarker; 1*ones(size(f,1),1)];

[f,v]=surf2patch(boxEl(1)*ones(size(Ylr)),Ylr,Zlr); %Right
F=[F;f+size(V,1)]; V=[V;v;]; faceBoundaryMarker=[faceBoundaryMarker; 2*ones(size(f,1),1)];

[Yfb,Zfb]=meshgrid(0:boxEl(1),0:boxEl(3));
[f,v]=surf2patch(Yfb,0*ones(size(Yfb)),Zfb); %Front
F=[F;f+size(V,1)]; V=[V;v;]; faceBoundaryMarker=[faceBoundaryMarker; 3*ones(size(f,1),1)];

[f,v]=surf2patch(Yfb,boxEl(2)*ones(size(Yfb)),Zfb); %Back
f=fliplr(f);
F=[F;f+size(V,1)]; V=[V;v;]; faceBoundaryMarker=[faceBoundaryMarker; 4*ones(size(f,1),1)];

[f,v]=surf2patch(Xtb,Ytb,0*ones(size(Xtb))); %Bottom
f=fliplr(f);
F=[F;f+size(V,1)]; V=[V;v;]; faceBoundaryMarker=[faceBoundaryMarker; 5*ones(size(f,1),1)];

[f,v]=surf2patch(Xtb,Ytb,boxEl(3)*ones(size(Xtb))); %Top
F=[F;f+size(V,1)]; V=[V;v;]; faceBoundaryMarker=[faceBoundaryMarker; 6*ones(size(f,1),1)];

%%
% Remove double nodes by merging
[F,V]=mergeVertices(F,V);

%%
% Scale coordinates
maxV=max(V,[],1);
V=V./maxV(ones(size(V,1),1),:);
V=V.*boxDim(ones(size(V,1),1),:);
V=V-boxDim(ones(size(V,1),1),:)/2;

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


text:
function s = trunc_singular(s, tol, relative, maxrank)
% REL_TRUNC_SINGULAR Helper routine to truncate singular values

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

    if ~exist('relative','var'),    relative = true;     end
    if ~exist('maxrank','var'), maxrank = length(s); end

    summ = cumsum(s.^2,'reverse');

    if relative
        s = find(summ > tol^2, 1, 'last');
        if isempty(s), s = 1; end
    else
        s = find(summ > tol^2*summ(1), 1, 'last');
        if isempty(s), s = 1; end
    end

    s = min([s, maxrank, length(s)]);

end

author: NicolasBoumal

repo: manopt


text:
function source = nmt_polaritytweak(cfg,source)

inside_idx = find(source.inside);

s = cell2mat(source.avg.mom(inside_idx));
ori = cell2mat(source.avg.ori(inside_idx));

if(isfield(cfg,'toilim'))
    tsel = dsearchn(source.time',cfg.toilim');
    tsel = tsel(1):tsel(2);
else
    [dum,tsel]=max(max(abs(s)));
end

% flip based on polarity of voxel with maximum power in desired time window
[dum,b]=max(nmt_rownorm(s(:,tsel)));
flipper=sign(s(b,tsel)*s(:,tsel)').*speye(size(s,1));
s = flipper*s;
ori = ori*flipper;

for ii=1:size(s,1)
    source.avg.mom{inside_idx(ii)} = s(ii,:);
    source.avg.ori{inside_idx(ii)} = ori(:,ii);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [data, warns] = psse_parse(records, sections, verbose, rev)
%PSSE_PARSE  Parses the data from a PSS/E RAW data file.
%   DATA = PSSE_PARSE(RECORDS, SECTIONS)
%   DATA = PSSE_PARSE(RECORDS, SECTIONS, VERBOSE)
%   DATA = PSSE_PARSE(RECORDS, SECTIONS, VERBOSE, REV)
%   [DATA, WARNINGS] = PSSE_PARSE(RECORDS, SECTIONS, ...)
%
%   Parses the data from a PSS/E RAW data file (as read by PSSE_READ)
%   into a struct.
%
%   Inputs:
%       RECORDS : cell array of strings, corresponding to the lines
%                 in the RAW file
%       SECTIONS : struct array with indices marking the beginning
%                  and end of each section, and the name of the
%                  section, fields are:
%           first   : index into RECORDS of first line of section
%           last    : index into RECORDS of last line of section
%           name    : name of the section, as extracted from the
%                     END OF ... DATA comments
%       VERBOSE      :  1 (default) to display progress info, 0 otherwise
%       REV          :  (optional) assume the input file is of this
%                       PSS/E revision number, attempts to determine
%                       REV from the file by default
%
%   Output(s):
%       DATA :  a struct with the following fields, each with two
%               sub-fields, 'num' and 'txt' containing the numeric and
%               text data read from the file for the corresponding section
%           id
%           bus
%           load
%           gen
%           shunt
%           branch
%           trans2
%           trans3
%           area
%           twodc
%           swshunt
%       WARNINGS :  cell array of strings containing accumulated
%                   warning messages
%
%   See also PSSE2MPC, PSSE_READ, PSSE_PARSE_SECTION, PSSE_PARSE_LINE

%   MATPOWER
%   Copyright (c) 2014-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   Based on mpreadraw.m, written by: Yujia Zhu, Jan 2014, yzhu54@asu.edu.
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% default args
if nargin < 4
    rev = 0;
    if nargin < 3
        verbose = 0;
    end
end
defaultrev = 23;

%% inititialize section counter
s = 1;
warns = {};

%%-----  case identification data  -----
%% In some files, SBASE is followed by a comment with the REV number such as:
%%  / PSS/E-29.0
%%  / PSS(tm)E-30 RAW
if verbose
    if rev
        fprintf('Forcing interpretation as PSS/E revision %d\n', rev);
    else
        fprintf('Attempting to determine PSS/E revision from content.\n');
    end
    fprintf('Parsing case identification data ...');
end
if rev
    warns{end+1} = sprintf('Conversion explicitly using PSS/E revision %d', rev);
end
[d, c] = psse_parse_line(records{1}, 'dfdfff');
nn = length(d);
data.id.IC = d{1};
if isempty(d{2}) || d{2} <= 0
    error('ERROR: Probable corrupt file, unable to read a valid SBASE value from the first line.');
else
    data.id.SBASE = d{2};
end
if ~isempty(d{3})
    data.id.REV = d{3};
else    %% attempt to extract revision from comment
    tmp = regexp(c, 'PSS(/|\(tm\))E-(?<rev>\d+)', 'tokens');
    if ~isempty(tmp) && size(tmp{1}, 2) == 2
        data.id.REV = str2num(tmp{1}{2});
    else
        data.id.REV = 0;
    end
end
if ~isempty(d{4})
    data.id.XFRRAT = d{4};
else
    data.id.XFRRAT = 0;
end
if ~isempty(d{5})
    data.id.NXFRAT = d{5};
else
    data.id.NXFRAT = 0;
end
if ~isempty(d{6})
    data.id.BASFRQ = d{6};
else
    data.id.BASFRQ = 0;
end
data.id.comment0 = c;
data.id.comment1 = records{2};
data.id.comment2 = records{3};
if verbose
    if rev
        if data.id.REV
            fprintf('.. override detected rev %2d w/%2d ... done.\n', data.id.REV, rev);
        else
            fprintf('...... unknown rev, using rev %2d ... done.\n', rev);
        end
    else
        if data.id.REV
            fprintf('......... rev %2d format detected ... done.\n', data.id.REV);
        else
            fprintf('...... unknown rev, using rev %2d ... done.\n', defaultrev);
        end
    end
end
if ~rev
    if data.id.REV
        rev = data.id.REV;      %% use detected value
    else
        rev = defaultrev;       %% none detected, use default value
        data.id.REV = defaultrev;
        warns{end+1} = sprintf('Unknown REV, using REV %2d format.', defaultrev);
    end
else
    data.id.REV = rev;          %% use override value
end
if isempty(data.id.IC) || data.id.IC ~= 0
    warns{end+1} = sprintf('IC = %d indicates that this may be a change case, rather than base case\n         PSSE2MPC is NOT designed to handle change cases.', data.id.IC);
    if verbose
        fprintf('WARNING : %s\n', warns{end});
    end
end
s = s + 1;

%%-----  bus data  -----
if rev < 24         %% includes load data
    [data.bus, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
        'bus', 'ddffffdffsfd');
elseif rev < 31     %% includes fixed shunt data, load separate
    [data.bus, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
        'bus', 'dsfdffddffd');
else                %% fixed shunt and load data separate
    [data.bus, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
        'bus', 'dsfddddffff..');
%       'bus', 'dsfddddffffff');
end
s = s + 1;

%%-----  load data  -----
if rev >= 24
    [data.load, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
        'load', 'd.d..ffffff...');
%       'load', 'dsdddffffffddd');
    s = s + 1;
end

%%-----  fixed shunt data  -----
if rev > 30     %% fixed shunt data is included in bus data for rev <= 30
    [data.shunt, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
        'fixed shunt', 'd.dff');
%       'fixed shunt', 'dsdff');
    s = s + 1;
end

%%-----  generator data  -----
[data.gen, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
    'generator', 'd.fffff.f.....d.ff...........');
%   'generator', 'dsfffffdffffffdfffdfdfdfdfsdf');
s = s + 1;

%%-----  branch data  -----
if rev <= 27   %% includes transformer ratio, angle
    [data.branch, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
        'branch', 'dd.ffffffffffffd');
%       'branch', 'dddffffffffffffd');
elseif rev < 31
    [data.branch, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
        'branch', 'dd.ffffffffffd');
%       'branch', 'ddsffffffffffdfdfdfdfdf');
else
    [data.branch, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
        'branch', 'dd.ffffffffffd');
%       'branch', 'ddsffffffffffddfdfdfdfdf');
end
s = s + 1;

%%-----  skip transformer adjustment data  -----
if rev <= 27
    [s, warns] = psse_skip_section(warns, sections, s, verbose, 'transformer adjustment');
end

%%-----  transformer data  -----
if rev > 27
    %% PSS/E stores two winding and three winding transformer data in the same
    %% section in RAW file. We read in 2 passes, first pass determines the type of
    %% each, second pass reads the data.
    label = 'transformer';
    if ~isempty(sections(s).name) && ~strcmpi(label, sections(s).name)
        if verbose > 1
            fprintf('-----  WARNING:  Expected section labeled: ''%s''\n', upper(label));
            fprintf('-----            Found section labeled:    ''%s''\n', sections(s).name);
        end
    end

    %% Step 1 : Count and collect transformer types
    if verbose
        fprintf('Analyzing transformer types ...');
    end

    %% estimate max number of transformers by number of lines in section
    nt2 = round((sections(s).last - sections(s).first + 1) / 4);
    nt3 = round((sections(s).last - sections(s).first + 1) / 5);

    %% initialize indices for transformer types
    idx2 = zeros(nt2, 1);
    idx3 = zeros(nt3, 1);

    %% set up counters
    i = sections(s).first;      %% initialize record index
    i2 = 0;
    i3 = 0;

    while i <= sections(s).last
        %% determine transformer type
        pat = '[^''",\s/]+\s*(,|\s)\s*[^''",\s/]+\s*(,|\s)\s*([^''",\s/]+)';
        m = regexp(records{i}, pat, 'tokens', 'once');
        if length(m) ~= 3
            disp(m);
            error('m should be length 3');
        end
        if length(m{3}) == 1 && m{3}(1) == '0'  %% two-winding
            i2 = i2 + 1;
            idx2(i2) = i;
            i = i + 4;
        else                                    %% three-winding
            i3 = i3 + 1;
            idx3(i3) = i;
            i = i + 5;
        end
    end
    nt2 = i2;
    nt3 = i3;

    if verbose
        str = sprintf(' %d(%d) two(three)-winding.', nt2, nt3);
        spacers = repmat('.', 1, 36-length(str));
        fprintf('%s %s ... done.\n', spacers, str);
    end

    %% trim index vectors down to size
    idx2 = idx2(1:nt2);
    idx3 = idx3(1:nt3);

    %% parse record 1 (cols 1-20)
    [t2_1, warns] = psse_parse_section(warns, records(idx2), verbose, ...
        '2-winding transformers (1)', 'dd..ddd....d........');
%       '2-winding transformers (1)', 'dddsdddffdsddfdfdfdf');
    [t3_1, warns] = psse_parse_section(warns, records(idx3), verbose, ...
        '3-winding transformers (1)', 'ddd.ddd....d........');
%       '3-winding transformers (1)', 'dddsdddffdsddfdfdfdf');

    %% two-winding
    %% parse record 2 (cols 21-23)
    [t2_2, warns] = psse_parse_section(warns, records(idx2+1), verbose, ...
        '2-winding transformers (2)', 'fff');

    %% parse record 3 (cols 24-39)
    %% parse up to CX1, should warn if CNXA1 is present and non-zero
    [t2_3, warns] = psse_parse_section(warns, records(idx2+2), verbose, ...
        '2-winding transformers (3)', 'ffffff..........');
%       '2-winding transformers (3)', 'ffffffddffffddff');

    %% parse record 4 (cols 40-41)
    [t2_4, warns] = psse_parse_section(warns, records(idx2+3), verbose, ...
        '2-winding transformers (4)', 'ff');

    %% three-winding
    %% parse record 2 (cols 21-31)
    [t3_2, warns] = psse_parse_section(warns, records(idx3+1), verbose, ...
        '3-winding transformers (2)', 'fffffffffff');
%       '3-winding transformers (2)', 'fffffffffff');

    %% parse record 3 (cols 32-47)
    %% parse up to CX1, should warn if CNXA1 is present and non-zero
    [t3_3, warns] = psse_parse_section(warns, records(idx3+2), verbose, ...
        '3-winding transformers (3)', 'ffffff..........');
%       '3-winding transformers (3)', 'ffffffddffffddff');

    %% parse record 4 (cols 48-63)
    %% parse up to CX2
    [t3_4, warns] = psse_parse_section(warns, records(idx3+3), verbose, ...
        '3-winding transformers (4)', 'ffffff..........');
%       '3-winding transformers (4)', 'ffffffddffffddff');

    %% parse record 5 (cols 64-79)
    %% parse up to CX3
    [t3_5, warns] = psse_parse_section(warns, records(idx3+4), verbose, ...
        '3-winding transformers (5)', 'ffffff..........');
%       '3-winding transformers (5)', 'ffffffddffffddff');

    %% assemble two-winding transformer records
    data.trans2.num = [t2_1.num(:, 1:20) t2_2.num(:, 1:3) t2_3.num(:, 1:16) t2_4.num(:, 1:2)];
    data.trans2.txt = [t2_1.txt(:, 1:20) t2_2.txt(:, 1:3) t2_3.txt(:, 1:16) t2_4.txt(:, 1:2)];

    %% assemble three-winding transformer records
    data.trans3.num = [t3_1.num(:, 1:20) t3_2.num(:, 1:11) t3_3.num(:, 1:16) t3_4.num(:, 1:16) t3_5.num(:, 1:16)];
    data.trans3.txt = [t3_1.txt(:, 1:20) t3_2.txt(:, 1:11) t3_3.txt(:, 1:16) t3_4.txt(:, 1:16) t3_5.txt(:, 1:16)];

    % if verbose
    %     fprintf('%s\n', upper(label));
    %     fprintf('%s\n', sections(s).name);
    % end
    s = s + 1;
end

%%-----  area interchange data  -----
[data.area, warns] = psse_parse_section(warns, records, sections, s, verbose, 'area', 'ddffs');
s = s + 1;

%%-----  two-terminal DC transmission line data  -----
label = 'two-terminal DC';
if ~isempty(sections(s).name) && ~strcmpi(label, sections(s).name)
    if verbose > 1
        fprintf('-----  WARNING:  Expected section labeled: ''%s''\n', upper(label));
        fprintf('-----            Found section labeled:    ''%s''\n', sections(s).name);
    end
end
idx = sections(s).first:3:sections(s).last;
if rev < 31
    [dc1, warns] = psse_parse_section(warns, records(idx), verbose, ...
        'two-terminal DC (1)', '.d.ff.......');
%       'two-terminal DC (1)', 'ddffffffsfdf');
    [dc2, warns] = psse_parse_section(warns, records(idx+1), verbose, ...
        'two-terminal DC (2)', 'd.ff.............');
%       'two-terminal DC (2)', 'ddffffffffffdddsf');
    [dc3, warns] = psse_parse_section(warns, records(idx+2), verbose, ...
        'two-terminal DC (3)', 'd.ff.............');
%       'two-terminal DC (3)', 'ddffffffffffdddsf');
else
    [dc1, warns] = psse_parse_section(warns, records(idx), verbose, ...
        'two-terminal DC (1)', '.d.ff.......');
%       'two-terminal DC (1)', 'sdffffffsfdf');
    [dc2, warns] = psse_parse_section(warns, records(idx+1), verbose, ...
        'two-terminal DC (2)', 'd.ff.............');
%       'two-terminal DC (2)', 'ddffffffffffdddDf');
    [dc3, warns] = psse_parse_section(warns, records(idx+2), verbose, ...
        'two-terminal DC (3)', 'd.ff.............');
%       'two-terminal DC (3)', 'ddffffffffffdddDf');
end

%% assemble two-terminal DC transmission line
data.twodc.num = [dc1.num dc2.num dc3.num];
data.twodc.txt = [dc1.txt dc2.txt dc3.txt];
% if verbose
%     fprintf('%s\n', upper(label));
%     fprintf('%s\n', sections(s).name);
% end
s = s + 1;

%%-----  skip voltage source converter data  -----
if rev > 28
    [s, warns] = psse_skip_section(warns, sections, s, verbose, 'voltage source converter');
end

%%-----  switched shunt data  -----
if rev < 31
    %% parse up to B1
    if rev <= 27
        [data.swshunt, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
            'switched shunt', 'd....f');
%           'switched shunt', 'ddffdfdfdfdfdfdfdfdfdf');
    elseif rev <= 29
        [data.swshunt, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
            'switched shunt', 'd.....f');
%           'switched shunt', 'ddffdsfdfdfdfdfdfdfdfdf');
    else    %%  rev == 30
        [data.swshunt, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
            'switched shunt', 'd......f');
%           'switched shunt', 'ddffdfsfdfdfdfdfdfdfdfdf');
    end
    s = s + 1;
end

%%-----  skip impedance correction data  -----
[s, warns] = psse_skip_section(warns, sections, s, verbose, 'impedance correction');

%%-----  skip multi-terminal DC data  -----
[s, warns] = psse_skip_section(warns, sections, s, verbose, 'multi-terminal DC');

%%-----  skip multi-section line data  -----
[s, warns] = psse_skip_section(warns, sections, s, verbose, 'multi-section line');

%%-----  skip zone data  -----
[s, warns] = psse_skip_section(warns, sections, s, verbose, 'zone');

%%-----  skip inter-area transfer data  -----
[s, warns] = psse_skip_section(warns, sections, s, verbose, 'inter-area transfer');

%%-----  skip owner data  -----
if rev > 24
    [s, warns] = psse_skip_section(warns, sections, s, verbose, 'owner');
end

%%-----  skip FACTS control device data  -----
if rev > 25
    [s, warns] = psse_skip_section(warns, sections, s, verbose, 'FACTS control device');
end

%%-----  switched shunt data  -----
if rev > 30
    %% parse up to B1
    if rev < 32
        [data.swshunt, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
            'switched shunt', 'd......f');
%           'switched shunt', 'ddffdfsfdfdfdfdfdfdfdfdf');
    else
        [data.swshunt, warns] = psse_parse_section(warns, records, sections, s, verbose, ...
            'switched shunt', 'd........f');
%           'switched shunt', 'ddddffdfsfdfdfdfdfdfdfdfdf');
    end
    s = s + 1;
end

%%-----  skip GNE device data  -----
if rev > 31
    [s, warns] = psse_skip_section(warns, sections, s, verbose, 'GNE device');
end

%%-----  skip induction machine data  -----
if rev > 32
    [s, warns] = psse_skip_section(warns, sections, s, verbose, 'induction machine');
end

%%-----  check for extra sections  -----
if s <= length(sections)
    warns{end+1} = sprintf('Found %d additional section(s)', length(sections)-s+1);
    if verbose > 1
        fprintf('-----  WARNING:   Found %d additional section(s):\n', length(sections)-s+1);
    end
end
while s <= length(sections)
    n = sections(s).last - sections(s).first + 1;
    if n
        str = sprintf('with %d line(s)', n);
    else
        str = sprintf('(empty)');
    end
    if isempty(sections(s).name)
        warns{end+1} = sprintf('  unlabeled section %s', str);
        if verbose > 1
            fprintf('-----            unlabeled section %s\n', str);
        end
    else
        warns{end+1} = sprintf('  ''%s DATA'' %s', sections(s).name, str);
        if verbose > 1
            fprintf('-----            ''%s DATA'' %s\n', sections(s).name, str);
        end
    end
    s = s + 1;
end

%%---------------------------------------------------------------------
function [s, warns] = psse_skip_section(warns, sections, s, verbose, label)
%PSSE_SKIP_SECTION  Skips over a section without extracting any data.
%   [SIDX, WARNINGS] = PSSE_SKIP_SECTION(WARNINGS, SECTIONS, SIDX, VERBOSE, LABEL)

if s > length(sections)
    if verbose
        spacers = repmat('.', 1, 58-length(label));
        fprintf('No %s data read %s done.\n', label, spacers);
    end
else
    nr = sections(s).last - sections(s).first + 1;
    if nr > 1
        ss = 'lines';
    else
        ss = 'line';
    end
    if nr
        warns{end+1} = sprintf('Skipped %d %s of %s data.', nr, ss, label);
    end
    if ~isempty(sections(s).name) && ~strcmp(upper(label), sections(s).name)
        warns{end+1} = sprintf('Section label mismatch, found ''%s'', expected ''%s''', ...
            sections(s).name, upper(label));
        if verbose
            fprintf('-----  WARNING:  Found section labeled:    ''%s''\n', sections(s).name);
            fprintf('-----            Expected section labeled: ''%s''\n', upper(label));
        end
    end
    if verbose && nr
        spacers = repmat('.', 1, 47-length(ss)-length(label));
        fprintf('Skipping%6d %s of %s data %s done.\n', nr, ss, label, spacers);
    end
    s = s + 1;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function[h]=regionplot(varargin)
%REGIONPLOT  Plots a simple box indicating a latitude / longitude region.
%
%   REGIONPLOT(REGION) draws a latitude / longitue box indicated by REGION
%   on the current axes.
%
%   REGION is an array with the format [WEST EAST SOUTH NORTH]. Longitudes
%   may either be specified on the interval [-180, 180] or on [0, 360].
%
%   By default, longitudes will be converted to [-180,180] unless the
%   region overlaps the dateline, in which case [0,360] will be used.
%
%   The current axes will be held as the box is plotted, then returned to
%   its previous HOLD state.
%
%   REGION can also be a cell array of N arrays, each with the format
%   REGION{N}=[WEST EAST SOUTH NORTH].  All N regions will then be plotted.
%
%   REGIONPLOT is compatible with the M_MAP toolbox, as described below.
%
%   H=REGIONPLOT outputs the line handle H to the plotted box.
%   __________________________________________________________________
%
%   Additional options
%
%   H=REGIONPLOT returns the graphics handle to the plotted box.
%
%   REGIONPLOT(REGION,STY) uses style string STY, in the LINESTYLE format,
%   with a default value of STY='k'.
%
%   REGIONPLOT(REGION,'M_MAP') will work with Rich Pawlowicz's M_MAP
%   package by calling M_PLOT.
%
%   The additional input parameters can be input together in either order,
%   i.e. REGIONPLOT(REGION,STY,'M_MAP') or REGIONPLOT(REGION,'M_MAP',STY).
%   __________________________________________________________________
%
%   See also INREGION, TRAJEXTRACT, TOPOPLOT.
%
%   'regionplot --f' generates some sample figures.
%
%   Usage: regionplot(region);
%          h=regionplot(region,'m_map');
%          h=regionplot(region,'m_map','2k');
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2013--2016 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--f')
    type makefigs_regionplot
    makefigs_regionplot;
    return
end
region=varargin{1};
varargin=varargin(2:end);
sty='k';
str='matlab';
for i=1:2
    if length(varargin)>=1
        if ischar(varargin{1})
            if strfind(lower(varargin{1}),'m_m')
                str=lower(varargin{1});
            else
                sty=varargin{1};
            end
            varargin=varargin(2:end);
        end
    end
end

if ~iscell(region)
    h=regionplot_one(region,str,sty);
    set(h,'visible','on')
    %set(h(~isnan(h)),'visible','on')
else
    bhold=ishold;
    set(gca,'visible','off')
    hold on
    storestate=get(gcf,'BackingStore');
    set(gcf,'BackingStore','off')
    for i=1:length(region)
         h(i)=regionplot_one(region{i},str,sty);
    end
    linestyle(h,'default');
    set(h,'visible','on')
    %set(h(~isnan(h)),'visible','on')

    set(gcf,'BackingStore',storestate)
    if ~bhold
        hold off
    end
     set(gca,'visible','on')
end

if nargout==0
    clear h
end

function[h]=regionplot_one(region,str,sty)

west=region(1);
east=region(2);
south=region(3);
north=region(4);

ax=axis;
if ax(1)>180 || deg180(west)>deg180(east)
    west=deg360(west);
    east=deg360(east);
else
    %west=deg180(west);
    %east=deg180(east);
end

bhold=ishold;
hold on
if strcmpi(str(1:3),'m_m')
     h=m_plot([west east east west west],[south south north north south],'visible','off');
else
     h=plot([west east east west west],[south south north north south],'visible','off');
end
linestyle(h,sty);
if ~bhold,
    hold off
end

%use ebasnfloats
%%Region shifted to be greater than 180
%region=[-30+250 -21+250 24 35];
%%Region shifted to overlap dateline
%region=[-30+205 -21+205 24 35];
%figure,cellplot(celladd(205,lon),lat),latratio,regionplot(region),axis tight
%title('Floats from Eastern Basin experiment')
%xlabel('Longitude'),ylabel('Latitude'),boxon
%ax=axis;

%reporttest('REGIONPLOT',aresame())

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% function P = nodal2center(y,m,mex)
%
% transfers a nodal grid to a cell-centered grid
% if nargin==1, builds P explicitely, else return results of P(y), matrix free; endif
% depending on numel(Y), the matrix free version also handles P'(Y)
%
% Input:
%   y        input points,  nodal or cell-centered
%   m        number of discretization points
%   mex      bool (use mex code for matrix-free implementation)
%
% Output:
%   P        the projection matrix P, if nargin == 2
%            P*y,                     if y is nodal
%            P'*y,                    if y is cell-centered
%==============================================================================

function P = nodal2center(y,m,mex)

if nargin == 0 % help and minimal example
  help(mfilename);
  runMinimalExample;
  P = 'endMinimalExample';
  return
end

Js = {[1 2 3], [2 1 3], [3 1 2]}; % different permutations of dimensions

if nargin==1, m = y; end
dim = length(m);

% Here starts the matrix-based code
if nargin==1
  % ----------------------------------
  % build the matrix P and return it
  % ----------------------------------
  av = @(i) spdiags(ones(m(i),1)*[ 1,1],[0,1],m(i),m(i)+1)/2;
  zero = sparse(prod(m),prod(m+1));
  switch dim
    case 2
      A = kron(av(2), av(1));
      P = sparse([A zero; zero A]);
    case 3
      A = kron(av(3), kron(av(2), av(1)));
      P = sparse([A zero zero; zero A zero; zero zero A]);
    otherwise
      error('Dimension must be either 2 or 3.')
  end
  return; % end of story
end

% Here starts the MEX based code
if (dim > 1) && exist(fullfile(FAIRpath,'kernel',[mfilename,'C.',mexext])) == 3,
  if numel(y) == length(m)*prod(m)
    % cell-centered ->  nodal
    status = true;
  else
    % nodal -> cell-centered
    status = false;
  end;
  P = nodal2centerC(y, m, dim, status);
  return;
end;

% Here starts the matrix-free code

if numel(y) == length(m)*prod(m)
  % -----------------------
  % cell-centered ->  nodal
  % -----------------------
  y = reshape(y,prod(m),[]);
  Z = zeros(prod(m+1),size(y,2));         % allocate memory
  for i=1:size(y,2)                       % run over all components of Y
    yi = reshape(y(:,i),m);
    for j=1:dim                           % run over all dimensions in yi
      % J = [j,setdiff(1:dim,j)];         % make j-th dimension first
      J = Js{j};                          % make j-th dimension first
      yi = permute(yi,J);

      zi = yi([1 1:end],:,:);
      zi(2:end-1,:,:) = zi(2:end-1,:,:) + yi(2:end,:,:);

      yi = .5* ipermute(zi,J);            % undo permutation
    end
    Z(:,i) = yi(:);                       % store i-th component
  end
else
  % -----------------------
  % nodal -> cell-centered
  % -----------------------
  % MATLAB implementation
  y = reshape(y,prod(m+1),[]);
  Z = zeros(prod(m),size(y,2));           % allocate memory
  for i=1:size(y,2)                       % run over all components of Y
    if dim==1                             % 1D case is simplr and does not admit oermutations
      P = 0.5*(y(1:end-1)+y(2:end));
      return
    else
      yi = reshape(y(:,i),m+1);           % reorganize Y
    end
    for j=1:dim                           % run over all dimensions in yi
      % J = [j,setdiff(1:dim,j)];         % make j-th dimension first
      J = Js{j};                          % make j-th dimension first
      yi = permute(yi,J);
      % average nodal to center
      yi = 0.5*(yi(1:end-1,:,:)+yi(2:end,:,:));
      yi = ipermute(yi,J);                % undo permutation
    end
    Z(:,i) = reshape(yi,[],1);            % store i-th component
  end

end
P = reshape(Z,[],1);                      % reshape
%------------------------------------------------------------------------------

function runMinimalExample

omega = [0 2 0 1];
m     = [8,7];
yn    = getNodalGrid(omega,m);
xc    = nodal2center(yn,m);
xc    = reshape(xc,[],2);

FAIRfigure(1); clf;
subplot(2,1,1); spy(nodal2center([4,5,6])); title('spy(nodal2gcenter operator)');
subplot(2,1,2); plotGrid(yn,omega,m,'color','b'); hold on; plot(xc(:,1),xc(:,2),'rx');
title(mfilename);

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


text:
function [g,tfr]=pgauss(L,varargin)
%PGAUSS  Sampled, periodized Gaussian
%   Usage: g=pgauss(L);
%          g=pgauss(L,tfr);
%          g=pgauss(L,...);
%          [g,tfr]=pgauss( ... );
%
%   Input parameters:
%         L    : Length of vector.
%         tfr  : ratio between time and frequency support.
%
%   Output parameters:
%         g    : The periodized Gaussian.
%
%   `pgauss(L,tfr)` computes samples of a periodized Gaussian. The function
%   returns a regular sampling of the periodization of the function
%   $exp(-pi*(x.^2/tfr))$.
%
%   The $l^2$ norm of the returned Gaussian is equal to 1.
%
%   The parameter *tfr* determines the ratio between the effective support
%   of *g* and the effective support of the DFT of *g*. If $tfr>1$ then *g*
%   has a wider support than the DFT of *g*.
%
%   `pgauss(L)` does the same setting *tfr=1*.
%
%   `[g,tfr] = pgauss( ... )` will additionally return the time-to-frequency
%   support ratio. This is useful if you did not specify it (i.e. used the
%   `'width'` or `'bw'` flag).
%
%   The function is whole-point even. This implies that `fft(pgauss(L,tfr))`
%   is real for any *L* and *tfr*. The DFT of *g* is equal to
%   `pgauss(L,1/tfr)`.
%
%   In addition to the `'width'` flag, `pgauss` understands the following
%   flags at the end of the list of input parameters:
%
%     'fs',fs     Use a sampling rate of *fs* Hz as unit for specifying the
%                 width, bandwidth, centre frequency and delay of the
%                 Gaussian. Default is *fs=[]* which indicates to measure
%                 everything in samples.
%
%     'width',s   Set the width of the Gaussian such that it has an
%                 effective support of *s* samples. This means that
%                 approx. 96% of the energy or 79% of the area
%                 under the graph is contained within *s* samples.
%                 This corresponds to -6dB or to width at the
%                 half of the height.
%                 This is equivalent to calling `pgauss(L,pi*s^2/4L*log(2))`.
%
%     'atheight',ah  Used only in conjuction with 'width'. Forces the
%                    Gaussian to width *s* at the *ah* fraction of the
%                    height.
%
%     'bw',bw     As for the `'width'` argument, but specifies the width
%                 in the frequency domain. The bandwidth is measured in
%                 normalized frequencies, unless the `'fs'` value is given.
%
%     'cf',cf     Set the centre frequency of the Gaussian to *fc*.
%
%     'wp'        Output is whole point even. This is the default.
%
%     'hp'        Output is half point even, as most Matlab filter
%                 routines.
%
%     'delay',d   Delay the output by *d*. Default is zero delay.
%
%   In addition to these parameteres, `pgauss` accepts any of the flags
%   from |setnorm|. The output will be normalized as specified.
%
%   If this function is used to generate a window for a Gabor frame, then
%   the window giving the smallest frame bound ratio is generated by
%   `pgauss(L,a*M/L)`.
%
%   Examples:
%   ---------
%
%   This example creates a Gaussian function, and demonstrates that it is
%   its own Discrete Fourier Transform:::
%
%     g=pgauss(128);
%
%     % Test of DFT invariance: Should be close to zero.
%     norm(g-dft(g))
%
%   The next plot shows the Gaussian in the time domain:::
%
%     plot(fftshift(pgauss(128)));
%
%   The next plot shows the Gaussian in the frequency domain on a log scale:::
%
%     magresp(pgauss(128),'dynrange',100);
%
%   The next plot shows the Gaussian in the time-frequency plane:::
%
%     sgram(pgauss(128),'tc','nf','lin');
%
%   See also:  dgtlength, psech, firwin, pbspline, setnorm
%
%   Demos:  demo_pgauss
%
%   References: mazh93

% AUTHOR : Peter L. Søndergaard.

%   First reference on this found in mazh93 eq. 63

if nargin<1
  error('Too few input parameters.');
end;

if (prod(size(L,1))~=1 || ~isnumeric(L))
  error('L must be a scalar');
end;

if rem(L,1)~=0
  error('L must be an integer.')
end;

% Define initial value for flags and key/value pairs.
definput.import={'setnorm'};

definput.flags.centering={'wp','hp'};
definput.flags.delay={'nodelay','delay'};
definput.flags.width={'tfr','width','bw'};

definput.keyvals.tfr=1;
definput.keyvals.delay=0;
definput.keyvals.width=0;
definput.keyvals.fs=[];
definput.keyvals.cf=0;
definput.keyvals.bw=0;
definput.keyvals.atheight=[];

[flags,keyvals,tfr]=ltfatarghelper({'tfr'},definput,varargin);

if (prod(size(tfr,1))~=1 || ~isnumeric(tfr))
  error('tfr must be a scalar.');
end;

if ~isempty(keyvals.atheight) && ~flags.do_width
    error(['%s: Param. ''atheight'' must be used together with param.',...
           ' ''width''. '],upper(mfilename));
end

if isempty(keyvals.atheight)
    keyvals.atheight = 0.5;
end

if keyvals.atheight >= 1 || keyvals.atheight <=0
    error('%s: Param. ''atheight'' must be in the range ]0,1[.',...
           upper(mfilename));
end

fs=keyvals.fs;

if flags.do_wp
  cent=0;
else
  cent=0.5;
end;

if isempty(fs)

  if flags.do_width
     tfr=pi/(4*log(1/keyvals.atheight))*keyvals.width^2/L;
  end;

  if flags.do_bw
     tfr=L/(keyvals.bw*L/2)^2;
  end;

  delay_s=keyvals.delay;
  cf_s   =keyvals.cf;
else

  if flags.do_width
    tfr=(keyvals.width*fs)^2/L;
  end;

  if flags.do_bw
    tfr=L/(keyvals.bw/fs*L)^2;
  end;

  delay_s=keyvals.delay*fs;
  cf_s   =keyvals.cf/fs*L;
end;

g=comp_pgauss(L,tfr,cent-delay_s,cf_s);

g=setnorm(g,flags.norm);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [] = writefort11( fort11dat, finame )
%
%

if ( nargin == 1 )
  finputname = 'fort.11_1' ;
else
  finputname = finame ;
end

fid = fopen(finputname,'w') ;

fprintf( fid, '%s\n', fort11dat.DataTitle ) ;
%fprintf( fid, '%s\n', fort11dat.DataSubTitle ) ;
fprintf( fid, '%d\n', fort11dat.DTIMINC ) ;

fprintf( fid, '%d\n', fort11dat.NumOfNodes ) ;

if iscell(fort11dat.Val)
    ValL = length(fort11dat.Val);
    valpernode = size(fort11dat.Val{1},1);
else
    % Write user-defined values
    [valpernode,~,ValL] = size( fort11dat.Val ) ;
end

% Looping over the number of times
for t = 1:ValL

    if iscell(fort11dat.Val)
        val = fort11dat.Val{t};
    else
        val = fort11dat.Val(:,:,t) ;
    end
    % Get a proper format
    str = '%d' ;
    for ll = 1: valpernode - 1
        str = [str ' %15.9e'] ;
    end
    str = [str '\n' ] ;

    fprintf( fid, str, val ) ;
end

fclose(fid) ;
%EOF
end

author: CHLNDDEV


text:
function tree_output = eval_cartree(Data,RETree)

%Evaluates the tree output

tree_output = mx_eval_cartree(Data,RETree.nodeCutVar,RETree.nodeCutValue,RETree.childnode,RETree.nodelabel);

author: Sable

repo: mcbench-benchmarks


text:
function [label, model, llh] = mixGaussEm(X, init)
% Perform EM algorithm for fitting the Gaussian mixture model.
% Input:
%   X: d x n data matrix
%   init: k (1 x 1) number of components or label (1 x n, 1<=label(i)<=k) or model structure
% Output:
%   label: 1 x n cluster label
%   model: trained model structure
%   llh: loglikelihood
% Written by Mo Chen (sth4nth@gmail.com).
%% init
fprintf('EM for Gaussian mixture: running ... \n');
tol = 1e-6;
maxiter = 500;
llh = -inf(1,maxiter);
R = initialization(X,init);
for iter = 2:maxiter
    [~,label(1,:)] = max(R,[],2);
    R = R(:,unique(label));   % remove empty clusters
    model = maximization(X,R);
    [R, llh(iter)] = expectation(X,model);
    if abs(llh(iter)-llh(iter-1)) < tol*abs(llh(iter)); break; end;
end
llh = llh(2:iter);

function R = initialization(X, init)
n = size(X,2);
if isstruct(init)  % init with a model
    R  = expectation(X,init);
elseif numel(init) == 1  % random init k
    k = init;
    label = ceil(k*rand(1,n));
    R = full(sparse(1:n,label,1,n,k,n));
elseif all(size(init)==[1,n])  % init with labels
    label = init;
    k = max(label);
    R = full(sparse(1:n,label,1,n,k,n));
else
    error('ERROR: init is not valid.');
end

function [R, llh] = expectation(X, model)
mu = model.mu;
Sigma = model.Sigma;
w = model.w;

n = size(X,2);
k = size(mu,2);
R = zeros(n,k);
for i = 1:k
    R(:,i) = loggausspdf(X,mu(:,i),Sigma(:,:,i));
end
R = bsxfun(@plus,R,log(w));
T = logsumexp(R,2);
llh = sum(T)/n; % loglikelihood
R = exp(bsxfun(@minus,R,T));

function model = maximization(X, R)
[d,n] = size(X);
k = size(R,2);
nk = sum(R,1);
w = nk/n;
mu = bsxfun(@times, X*R, 1./nk);

Sigma = zeros(d,d,k);
r = sqrt(R);
for i = 1:k
    Xo = bsxfun(@minus,X,mu(:,i));
    Xo = bsxfun(@times,Xo,r(:,i)');
    Sigma(:,:,i) = Xo*Xo'/nk(i)+eye(d)*(1e-6);
end

model.mu = mu;
model.Sigma = Sigma;
model.w = w;

function y = loggausspdf(X, mu, Sigma)
d = size(X,1);
X = bsxfun(@minus,X,mu);
[U,p]= chol(Sigma);
if p ~= 0
    error('ERROR: Sigma is not PD.');
end
Q = U'\X;
q = dot(Q,Q,1);  % quadratic term (M distance)
c = d*log(2*pi)+2*sum(log(diag(U)));   % normalization constant
y = -(c+q)/2;


text:
function tc_saveParfiles(tc);
% tc_saveParfiles(tc);
%
% Time Course UI: save parfiles specified in tc struct.
%
% If an onset shift is specified in the 'onsetDelta' param, this is
% applied to the parfiles before saving.
%
% 07/04 ras.
global HOMEDIR;

dlgTitle = 'Save parfiles...';
for i = 1:length(tc.params.parfiles)
    prompt{i} = sprintf('Scan %i parfile name:',tc.params.scans(i));
    defaults{i} = tc.params.parfiles{i};
end

answer = inputdlg(prompt,dlgTitle,1,defaults);

scans = tc.params.scans;

for i = 1:length(answer)
    ind = find(tc.trials.run==i);
    par.cond = tc.trials.cond(ind);
    D = tc.params.onsetDelta;
    if i==1
        par.onset = tc.trials.onsetSecs(ind) + D;
    else
        par.onset = tc.trials.onsetSecs(ind) + D - tc.trials.onsetSecs(ind(1));
    end
    par.label = tc.trials.label(ind);
    par.color = cell(size(par.label));

    % use the currently-assigned cond names for the first label
    % for each condition
    for c = unique(par.cond)
        subind = find(par.cond==c);
        par.label{subind(1)} = tc.trials.condNames{find(tc.trials.condNums==c)};
        par.color{subind(1)} = tc.trials.condColors{find(tc.trials.condNums==c)};
    end

    parPath = fullfile(parfilesDir(initHiddenInplane), answer{i});
    writeParfile(par, parPath);
end

% because we applied the onset shift in the parfiles, we can set the
% onsetDelta param back to 0, and save it, so we don't shift twice:
if tc.params.onsetDelta ~= 0
    tc.params.onsetDelta = 0;
    hI = initHiddenInplane(tc.params.dataType, tc.params.scans(1));
    er_setParams(hI, tc.params);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [x, y, z, it]=icosahedron_nodes(f, r)
% % Syntax;
% %
% % [x, y, z, it]=icosahedron_nodes(f, r);
% %
% % ***********************************************************
% %
% % Description
% %
% % Program makes the master triangle for an icosahedron of
% % f frequency of type 2.
% %
% % ***********************************************************
% %
% % Input Variables
% %
% % f is the frquency of the icosahedron.
% % f=2 is similar to a soccer ball.
% %
% % r is the radius of the icosahedron.
% %
% % ***********************************************************
% %
% % Output Variables
% %
% % x, y, and z are vectors with the rectangular coordinates for the nodes
% %             of a hemispherical icosahedron.
% %
% % it is the matrix of coordinate indices for the rectangular coordinate
% %    vectors.
% %
% % ***********************************************************
% %
%
% Example
%
% f=2;  % 2 frequency icosahedron.
%
% r=4;  % 4 meter radius.
%
% [x2, y2, z2, it]=icosahedron_nodes(f, r);
%
% %
% % ***********************************************************
% %
% % This program was written by Edward L. Zechmann
% %
% %     date     January 2007
% %
% % modified 11   March   2008  added examples
% %
% % ***********************************************************
% %
% % Feel free to modify this code.
% %

%calculate the coordinate indices
% i1 -> x
% i2 -> y
% i3 -> z

num_pts=0.5*((f+1)^2+(f+1));

it=zeros(3,num_pts);

for e1=1:(f+1);
    i1=f-e1+1;
    for e2=1:(e1);
        i2=(e2-1);
        i3=(f-i1-i2);

        a=e1+1;
        b=f+1;
        if a <= f+1;
            pt=0.5*(b^2+b-a^2+a)+i2+1;
        else
            pt=i2+1;
        end

        it(1, pt)=i1;
        it(2, pt)=i2;
        it(3, pt)=i3;
    end
end

x1=1:num_pts;
y1=1:num_pts;
z1=1:num_pts;

% Coordinate formulas for base pentagon
% Icosahedron class 1, type 2, frequency f
tau=0.5*(1+sqrt(5));

for e1=1:num_pts;

    i1=it(1, e1);
    i2=it(2, e1);
    i3=it(3, e1);

    x1(e1)=i1*sin(72/180*pi);
    y1(e1)=i2+i1*cos(72/180*pi);
    z1(e1)=f/2+i3/tau;

end

[rho, theta, phi]=spherical_angle_ed(x1, y1, z1);
[x, y, z]=spherical_to_rectangular(r, theta, phi);