text:
%REAL2RGB  Converts a real-valued matrix into a truecolor image
%
% Examples:
%   B = real2rgb(A, cmap);
%   B = real2rgb(A, cmap, lims);
%   [B lims map] = real2rgb(...);
%
% This function converts a real-valued matrix into a truecolor image (i.e.
% double array with values between 0 and 1) using the colormap specified
% (either user-defined or the name of a colormap function). The output
% image is suitable for display using IMAGE or IMSHOW, exporting using
% IMWRITE, texture mapping a surface etc.
%
% Colormaps specified by name, e.g. 'hot', can be reversed ('-hot'), made
% to convert linearly to grayscale when printed on a black & white printer
% ('hot*'), or both ('-hot*').
%
% Value limits and a colormap table can be output, for use generating the
% correct colorbar, e.g.:
%   [B lims map] = real2rgb(peaks(256), '-hot*');
%   hIm = imshow(B);
%   set(gcf, 'Colormap', map);
%   set(gca, 'CLim', lims);
%   set(hIm, 'CDataMapping', 'scaled');
%   colorbar;
%
% IN:
%   A - MxN real matrix.
%   cmap - JxK user-defined colormap, or a string indicating the name
%          of the colormap to be used. K = 3 or 4. If K == 4 then
%          cmap(1:end-1,4) contains the relative widths of the bins between
%          colors. If cmap is a colormap function name then the prefix '-'
%          indicates that the colormap is to be reversed, while the suffix
%          '*' indicates that the colormap bins are to be rescaled so that
%          each bin produces the same change in gray level, such that the
%          colormap converts linearly to grayscale when printed in black
%          and white.
%   lims - 1x2 array of saturation limits to be used on A. Default:
%          [min(A(:)) max(A(:))].
%
% OUT:
%   B - MxNx3 truecolor image.
%   lims - 1x2 array of saturation limits used on A. Same as input lims, if
%          given.
%   map - 256x3 colormap similar to that used to generate B.

% Copyright: Oliver Woodford, 2009-2011

% Thank you to Peter Nave for reporting a bug whereby colormaps larger than
% 256 entries long are returned.

function [B, lims, map] = real2rgb(A, cmap, lims)
% Don't do much if A is wrong size
[y, x, c] = size(A);
if c > 1
    error('A can only have 2 dimensions');
end
if y*x*c == 0
    % Create an empty array with the correct dimensions
    B = zeros(y, x, (c~=0)*3);
    return
end
if nargin < 2
    cmap = 'parula';
end

% Generate the colormap
if ischar(cmap)
    % If map starts with a '-' sign, invert the colormap
    reverseMap = cmap(1) == '-';
    % If the map ends with a '*', attempt to make map convert linearly to
    % grayscale
    grayMap = cmap(end) == '*';
    % Extract the map name
    cmap = lower(cmap(reverseMap+1:end-grayMap));
    % Load the map
    try
        % Check for a concise table first
        map = feval(cmap, Inf);
    catch
        map = [];
    end
    if invalid_map(map)
        try
            % Just load a large table
            map = feval(cmap, 256);
        catch
            error('Colormap ''%s'' not found', cmap);
        end
        if invalid_map(map)
            error('Invalid colormap');
        end
    end
    if reverseMap
        % Reverse the map
        map = map(end:-1:1,:);
        if size(map, 2) == 4
            % Shift up the bin lengths
            map(1:end-1,4) = map(2:end,4);
        end
    end
    if grayMap && size(map, 1) > 2
        % Ensure the map converts linearly to grayscale
        map(1:end-1,4) = abs(diff(map(:,1:3) * [0.299; 0.587; 0.114]));
    end
else
    % Table-based colormap given
    map = cmap;
end

% Only work with real doubles
B = reshape(double(real(A)), y*x, c);

% Compute limits and scaled values
maxInd = 1 + (size(map, 1) - 2) * (size(map, 2) ~= 4);
if nargin < 3
    lims = [];
end
[B, lims] = rescale(B, lims, [0 maxInd]);

% Compute indices and offsets
if size(map, 2) == 4
    % Non-linear colormap
    bins = map(1:end-1,4);
    cbins = cumsum(bins);
    bins(bins==0) = 1;
    bins = cbins(end) ./ bins;
    cbins = [0; cbins(1:end-1) ./ cbins(end); 1+eps];
    [ind, ind] = histc(B, cbins);
    B = (B - cbins(ind)) .* bins(ind);
    clear bins cbins
else
    % Linear colormap
    ind = min(floor(B), maxInd-1);
    B = B - ind;
    ind = ind + 1;
end

% Compute the output image
try
    B = bsxfun(@times, map(ind,1:3), 1 - B) + bsxfun(@times, map(ind+1,1:3), B);
catch
    % If no bsxfun
    B = B(:,[1 1 1]);
    B = map(ind,1:3) .* (1 - B) + map(ind+1,1:3) .* B;
end
B = min(max(B, 0), 1); % Rounding errors can make values slip outside bounds
B = reshape(B, y, x, 3);

if nargout > 2 && (size(map, 1) ~= 256 || size(map, 2) == 4)
    % Generate the colormap (for creating a colorbar with)
    map = reshape(real2rgb(0:255, map, [0 255]), 256, 3);
end
return

function notmap = invalid_map(map)
notmap = isempty(map) || ndims(map) ~= 2 || size(map, 1) < 1 || size(map, 2) < 3 || size(map, 2) > 4 || ~all(reshape(map(:,1:3) >= 0 & map(:,1:3) <= 1, [], 1));

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
% SOSDEMO9p --- Matrix SOS decomposition
% Section 3.9 of SOSTOOLS User's Manual
%
clear; echo on;
pvar x1 x2 x3;
vartable = [x1, x2, x3];
% =============================================
% First, initialize the sum of squares program
prog = sosprogram(vartable);   % No decision variables.

% =============================================
% Consider the following candidate sum of squares matrix P(x)

P = [x1^4+x1^2*x2^2+x1^2*x3^2 x1*x2*x3^2-x1^3*x2-x1*x2*(x2^2+2*x3^2);
     x1*x2*x3^2-x1^3*x2-x1*x2*(x2^2+2*x3^2) x1^2*x2^2+x2^2*x3^2+(x2^2+2*x3^2)^2];

% Test if P(x1,x2,x3) is an SOS matrix and return H so that P = H.'*H
[Q,Z,H] = findsos(P);

% Verify that P - H'*H = 0 and P- kron(I,Z)'*Q*kron(I,Z)= 0 to within
% numerical tolerance.
tol = 1e-8;
cleanpoly(P - H'*H ,tol)
cleanpoly(P- kron(eye(2),Z)'*Q*kron(eye(2),Z), tol)

% Verify that Q >=0
min(eig(Q))

% =============================================
% Program is feasible, thus P(x1,x2,x3) is an SOS matrix.

echo off;

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
%typ=0->ecef 2 nav || typ=1->nav 2 ecef (nav=Llh in rad)
% from my observation, ecef2geo is almost as accurate as Borkowski and
% Hirvonen iterative scheme, though no so accurate in height, 1e-4 m.
function [out]=ecef2geo(inp,typ)
%Earth Defs.
E_SQR=0.00669437999014;
SM_AXIS=6378137;
FLATTENING=0.00335281066475;

if (typ==0) %ecef2Llh
    X=inp(1);
    Y=inp(2);
    Z=inp(3);

    p=sqrt(X^2+Y^2);
    b=SM_AXIS*(1-FLATTENING);
    ep=((SM_AXIS^2)-(b^2))/b^2;
    theta=atan2(Z*SM_AXIS,p*b);

    l=atan2(Y,X);
    L=atan2((Z+ep*b*(sin(theta))^3),(p-E_SQR*SM_AXIS*(cos(theta))^3));
    Re=SM_AXIS/sqrt(1-E_SQR*sin(L)^2);
    h=p/cos(L)-Re;

    out=[L;l;h];
else %Llh2ecef
    sL=sin(inp(1));
    cL=cos(inp(1));
    sl=sin(inp(2));
    cl=cos(inp(2));
    h=inp(3);

    Re=SM_AXIS/(sqrt(1.0-E_SQR*sL*sL));
    X=(Re+h)*cL*cl;
    Y=(Re+h)*cL*sl;
    Z=(Re*(1-E_SQR)+h)*sL;
    out=[X;Y;Z];
end

repo: ekfmonoslam


text:
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
%   [all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
%   logisitc regression classifiers and returns each of these classifiers
%   in a matrix all_theta, where the i-th row of all_theta corresponds
%   to the classifier for label i

% Some useful variables
m = size(X, 1);
n = size(X, 2);

% You need to return the following variables correctly
all_theta = zeros(num_labels, n + 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
%               logistic regression classifiers with regularization
%               parameter lambda.
%
% Hint: theta(:) will return a column vector.
%
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell use
%       whether the ground truth is true/false for this class.
%
% Note: For this assignment, we recommend using fmincg to optimize the cost
%       function. It is okay to use a for-loop (for c = 1:num_labels) to
%       loop over the different classes.
%
%       fmincg works similarly to fminunc, but is more efficient when we
%       are dealing with large number of parameters.
%
% Example Code for fmincg:
%
%     % Set Initial theta
%     initial_theta = zeros(n + 1, 1);
%
%     % Set options for fminunc
%     options = optimset('GradObj', 'on', 'MaxIter', 50);
%
%     % Run fmincg to obtain the optimal theta
%     % This function will return theta and the cost
%     [theta] = ...
%         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
%                 initial_theta, options);
%

for c = 1:num_labels
     initial_theta = zeros(n + 1, 1);
     options = optimset('GradObj', 'on', 'MaxIter', 50);
     [theta] = fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), initial_theta, options);
     all_theta(c, :) = theta';
end;

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

end

author: zsiciarz

repo: ml-coursera


text:
%L0SMOOTH  Global image smoothing via L0 gradient minimization
%
%     dst = cv.l0Smooth(src)
%     dst = cv.l0Smooth(src, 'OptionName',optionValue, ...)
%
% ## Input
% * __src__ source image for filtering with unsigned 8-bit or signed 16-bit or
%   floating-point depth.
%
% ## Output
% * __dst__ destination image.
%
% ## Options
% * __Lambda__ parameter defining the smooth term weight. default 0.02
% * __Kappa__ parameter defining the increasing factor of the weight of the
%   gradient data term. default 2.0
%
% For more details about L0 Smoother, see the original paper [xu2011image].
%
% ## References
% [xu2011image]:
% > Li Xu, Cewu Lu, Yi Xu, and Jiaya Jia. "Image smoothing via L0 gradient
% > minimization". In ACM Transactions on Graphics (TOG), volume 30,
% > page 174. ACM, 2011.
%
% See also: cv.jointBilateralFilter
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function hydro = radiationIRFSS(hydro,Omax,R2t)
% Calculates the state space (SS) realization of the normalized radiation IRF.
% If this function is used, it must be implemented after the radiationIRF function.
%
% Default parameters can be used by inputting [].
% See ``WEC-Sim\examples\BEMIO`` for examples of usage.
%
% Parameters
% ----------
%     hydro : struct
%         Structure of hydro data
%
%     Omax : integer
%         Maximum order of the SS realization, the default is 10
%
%     R2t : float
%         :math:`R^2` threshold (coefficient of determination) for the SS
%         realization, where :math:`R^2` may range from 0 to 1, and the
%         default is 0.95
%
% Returns
% -------
%     hydro : struct
%         Structure of hydro data with radiation IRF state space
%         coefficients
%

p = waitbar(0,'Calculating state space radiation IRFs...');  % Progress bar

% Set defaults if empty
if isempty(Omax)==1;    Omax = 10;  end
if isempty(R2t)==1;     R2t = 0.95; end

t = hydro.ra_t;
dt = t(2)-t(1);
hydro.ss_A = zeros(sum(hydro.dof),sum(hydro.dof),Omax,Omax);
hydro.ss_B = zeros(sum(hydro.dof),sum(hydro.dof),Omax,1);
hydro.ss_C = zeros(sum(hydro.dof),sum(hydro.dof),1,Omax);
hydro.ss_D = zeros(sum(hydro.dof),sum(hydro.dof),1);
hydro.ss_K = zeros(sum(hydro.dof),sum(hydro.dof),length(t));
hydro.ss_conv = zeros(sum(hydro.dof),sum(hydro.dof));
hydro.ss_R2 = zeros(sum(hydro.dof),sum(hydro.dof));
hydro.ss_O = zeros(sum(hydro.dof),sum(hydro.dof));

for i=1:sum(hydro.dof)
    for j=1:sum(hydro.dof)

        K = squeeze(hydro.ra_K(i,j,:));
        R2i = norm(K-mean(K));  % Initial R2

        %Hankel Singular Value Decomposition
        O = 2;  % Initial state space order
        y = dt*K;
        h = hankel(y(2:end));
        [u,svh,v] = svd(h);
        svh=diag(svh);

        while R2i ~= 0.0
            u1 = u(1:length(K)-2,1:O);
            v1 = v(1:length(K)-2,1:O);
            u2 = u(2:length(K)-1,1:O);
            sqs = sqrt(svh(1:O));
            ubar = u1.'*u2;

            a = ubar.*((1./sqs)*sqs.');
            b = v1(1,:).'.*sqs;
            c = u1(1,:).*sqs.';
            d = y(1);

            iidd = inv(dt/2*(eye(O)+a));  % (T/2*I+T/2*A)^{-1} = 2/T(I+A)^{-1}
            ac = (a-eye(O))*iidd;         % (A-I)2/T(I+A)^{-1} = 2/T(A-I)(I+A)^{-1}
            bc = dt*(iidd*b);             % (T/2+T/2)*2/T(I+A)^{-1}B = 2(I+A)^{-1}B
            cc = c*iidd;                  % C*2/T(I+A)^{-1} = 2/T(I+A)^{-1}
            dc = d-dt/2*((c*iidd)*b);     % D-T/2C (2/T(I+A)^{-1})B = D-C(I+A)^{-1})B
            for k=1:length(t)
                ss_K(k) = ((cc*expm(ac*dt*(k-1)))*bc);  % Calc SS IRF approx
            end

            R2 = 1-(norm(K-ss_K.')/R2i)^2;  % Calc R2 for SS IRF approx
            if R2 >= R2t  status = 1;  break  % R2 threshold
            elseif O == Omax  status = 2;  break  % Max SS order threshold
            else  O = O+1;  end  % Increase state space order
        end
        if R2i ~= 0.0
            hydro.ss_A(i,j,1:O,1:O) = ac;
            hydro.ss_B(i,j,1:O,1) = bc;
            hydro.ss_C(i,j,1,1:O) = cc;
            hydro.ss_D(i,j) = dc;
            hydro.ss_K(i,j,:) = ss_K;
            hydro.ss_conv(i,j) = status;
            hydro.ss_R2(i,j) = R2;
            hydro.ss_O(i,j) = O;
        end
    end
    waitbar(i/(sum(hydro.dof)))
end
close(p)

end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [v, g] = fv_obj_func(w, num_threads)
% Get the object function value and gradient at w.
%   [v, g] = fv_obj_func(w, num_threads)
%
% Return values
%   v             Objective function value f(w)
%   g             Objective function gradient \nabla f(w)
%
% Arguments
%   w             Gradient and function evaluation point
%   num_threads   Number of worker threads to use for computing the gradient

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

[v, g] = fv_cache('gradient', w, num_threads);

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function [mtv,xx1,xx2] = vbmc_mtv(vp1,vp2,Ns)
%VBMC_MTV Marginal Total Variation distances between two variational posteriors.
%   MTV = VBMC_MTV(VP1,VP2) returns an estimate of the marginal total
%   variation distances between two variational posterior distributions VP1
%   and VP2. MTV is a D-element vector whose elements are the total variation
%   distance between the marginal distributions of VP1 and VP2, for each
%   coordinate dimension.
%
%   The total variation distance between two densities p1 and p2 is:
%       TV(p1, p2) = 1/2 \int | p1(x) - p2(x) | dx
%
%   MTV = VBMC_MTV(VP1,VP2,NS) uses NS random draws to estimate the MTV
%   (default NS=1e5).
%
%   [MTV,XX1,XX2] = VBMC_MTV(...) returns NS samples from the variational
%   posteriors VP1 and VP2 as, respectively, NS-by-D matrices XX1 and XX2,
%   where D is the dimensionality of the problem.
%
%   VP1 and/or VP2 can be N-by-D matrices of samples from variational
%   posteriors (they do not need have the same number of samples).
%
%   See also VBMC, VBMC_KLDIV, VBMC_PDF, VBMC_RND, VBMC_DIAGNOSTICS.

if nargin < 3 || isempty(Ns); Ns = 1e5; end

% This was removed because the comparison *has* to be in original space,
% given that the transform might change for distinct variational posteriors
% if nargin < 4 || isempty(origflag); origflag = true; end
origflag = true;

if vbmc_isavp(vp1)
    xx1 = vbmc_rnd(vp1,Ns,origflag,1);
    lb1 = vp1.trinfo.lb_orig;
    ub1 = vp1.trinfo.ub_orig;
else
    xx1 = vp1;
    lb1 = -Inf(1,size(vp1,2));
    ub1 = Inf(1,size(vp1,2));
end
if vbmc_isavp(vp2)
    xx2 = vbmc_rnd(vp2,Ns,origflag,1);
    lb2 = vp2.trinfo.lb_orig;
    ub2 = vp2.trinfo.ub_orig;
else
    xx2 = vp2;
    lb2 = -Inf(1,size(vp2,2));
    ub2 = Inf(1,size(vp2,2));
end

D = size(xx1,2);
nkde = 2^13;
mtv = zeros(1,D);

% Set bounds for kernel density estimate
lb1_xx = min(xx1); ub1_xx = max(xx1);
range1 = ub1_xx - lb1_xx;
lb1 = max(lb1_xx-range1/10,lb1);
ub1 = min(ub1_xx+range1/10,ub1);

lb2_xx = min(xx2); ub2_xx = max(xx2);
range2 = ub2_xx - lb2_xx;
lb2 = max(lb2_xx-range2/10,lb2);
ub2 = min(ub2_xx+range2/10,ub2);

% Compute marginal total variation
for i = 1:D
    [~,yy1,x1mesh] = kde1d(xx1(:,i),nkde,lb1(i),ub1(i));
    yy1 = yy1/(qtrapz(yy1)*(x1mesh(2)-x1mesh(1))); % Ensure normalization

    [~,yy2,x2mesh] = kde1d(xx2(:,i),nkde,lb2(i),ub2(i));
    yy2 = yy2/(qtrapz(yy2)*(x2mesh(2)-x2mesh(1))); % Ensure normalization

    f = @(x) abs(interp1(x1mesh,yy1,x,'spline',0) - interp1(x2mesh,yy2,x,'spline',0));
    bb = sort([x1mesh([1,end]),x2mesh([1,end])]);
    for j = 1:3
        xx_range = linspace(bb(j),bb(j+1),1e5);
        mtv(i) = mtv(i) + 0.5*qtrapz(f(xx_range))*(xx_range(2)-xx_range(1));
    end
end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function supported_labels = getSupportedLabels()
% auxilary function

supported_labels={'SLC','SEG','REP','AVG','SET','ECO','PHS','LIN','PAR','NAV','REV','SMS', ... # these label values are copied to the corresponding MDH field on Siemens
        'PMC', ... # for MoCo/PMC Pulseq version to recognize blocks that can be prospectively corrected for motion
        'NOROT','NOPOS','NOSCL', ... # instruct the interpreter to ignore the position, rotation or scaling of the FOV specified on the UI
        'ONCE' ... # a 3-state flag that instructs the interpreter to alter the sequence when executing multiple repeats as follows: blocks with ONCE==0 are executed on every repetition; ONCE==1: only the first repetition; ONCE==2: only the last repetition
    };

end

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function value = year_length_months_hebrew ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_MONTHS_HEBREW returns the number of months in a Hebrew year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of months in the year.
%
  if ( year_is_embolismic_hebrew ( y ) )
    value = 13;
  else
    value = 12;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
filename='Cantilever_triangle_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function pf = loadPoleFigure_labotex(fname,varargin)
% load labotex pole figure data
%
% Input
% fname - file name
%
% Output
% pf    - @PoleFigure
%
% See also
% ImportPoleFigureData loadPoleFigure

% open the file
fid = efopen(fname);

% read header

try

  % first line --> comments
  comment = fgetl(fid);
  % skip next lines
  fgetl(fid);
  s = fgetl(fid);
  assert(strncmpi('Structure Code',s,14));

  % crystal symmetry
  s = fgetl(fid);
  c = sscanf(s,'%d%f%f%f%f%f%f');
  spacegroup =  {'C1','C2','D2','C4','D4','T','O','C3','D3','C6','D6'};
  cs = crystalSymmetry(spacegroup{c(1)},c(2:4),c(5:7)*degree);

  % number of pole figures
  s = fgetl(fid);
  npf = sscanf(s,'%d');

  % skip one line
  s = fgetl(fid);

  % hr information
  for n = 1:npf
    s = fgetl(fid);
    hr = sscanf(s,'%f%f%f%f%f%f%f%f%d%d%d%d%d');
    theta = (hr(2):hr(4):hr(3))*degree;
    rho = (hr(5):hr(7):hr(6))*degree;
    r{n} = regularS2Grid('theta',theta,'rho',rho,'antipodal'); %#ok<AGROW>
    h{n} = Miller(hr(9),hr(10),hr(11),cs); %#ok<AGROW>
  end

  for n = 1:npf

    % read data
    tmp = textscan(fid,'%f',length(r{n}));
    d{n} = reshape(double(tmp{1}),size(r{n}));

  end

  % construct pole figure object
  pf = PoleFigure(h,r,d,cs,'comment',comment,varargin{:});

catch %#ok<CTCH>
  interfaceError(fname,fid);
end

fclose(fid);

author: mtex-toolbox

repo: mtex

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


text:
%SECS2HMS - converts a time in seconds to a string giving the time in hours, minutes and second
%Usage TIMESTRING = SECS2HMS(TIME)]);
%Example 1: >> secs2hms(7261)
%>> ans = 2 hours, 1 min, 1.0 sec
%Example 2: >> tic; pause(61); disp(['program took ' secs2hms(toc)]);
%>> program took 1 min, 1.0 secs

function time_string=secs2hms(time_in_secs)
    time_string='';
    nhours = 0;
    nmins = 0;
    if time_in_secs >= 3600
        nhours = floor(time_in_secs/3600);
        if nhours > 1
            hour_string = ' hours, ';
        else
            hour_string = ' hour, ';
        end
        time_string = [num2str(nhours) hour_string];
    end
    if time_in_secs >= 60
        nmins = floor((time_in_secs - 3600*nhours)/60);
        if nmins > 1
            minute_string = ' mins, ';
        else
            minute_string = ' min, ';
        end
        time_string = [time_string num2str(nmins) minute_string];
    end
    nsecs = time_in_secs - 3600*nhours - 60*nmins;
    time_string = [time_string sprintf('%2.1f', nsecs) ' secs'];
end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
%         I=imread('test1.bmp');
%         mapping=getmapping(8,'u2');
%         H1=lbp(I,1,8,mapping,'h'); %LBP histogram in (8,1) neighborhood %using uniform patterns
%         subplot(1,1,1),stem(H1);
      I=imread('test1.bmp');
      row=size(I,1);
      col=size(I,2);
      B=mat2cell(I,[row/2 row/2],[col/2 col/2]);
      mapping=getmapping(8,'u2');
      H.a=0;
      for i=1:4
      H1=lbp(B{i},1,8,mapping,'h'); %LBP histogram in (8,1) neighborhood %using uniform patterns
      H.hist{i}=H1;
      end
      subplot(2,2,1),stem(H.hist{1});
      subplot(2,2,2),stem(H.hist{2});
      subplot(2,2,3),stem(H.hist{3});
      subplot(2,2,4),stem(H.hist{4});

author: xiaomangui

repo: LBP-DBN-face-recognition


text:
function [entropy, samples, predicted_label]= get_recognition_entropy(model, tsdf_data, param)
% recognition entropy of the incomplete TSDF data.
% tsdf_data: contains surface voxels(+1), empty space(0), and unobserved
% space(-1).
% param: testing parameters. (sampling epochs, number of particles,
% earlyStop or not).
% predicted_label: p(y|tsdf).
% entropy: entropy of predicted_label.
% samples: completed samples of tsdf.

if ~exist('param','var')
    param.nparticles = 32;
    param.epochs = 20;
    param.earlyStop = true;
end

nparticles = param.nparticles;

sample_param = [];
sample_param.epochs = param.epochs;
sample_param.nparticles = nparticles;
sample_param.gibbs_iter = 1;
sample_param.earlyStop = param.earlyStop;

batch_data = repmat(permute(tsdf_data,[4,1,2,3]),nparticles,1); % running n chains altogether
mask = batch_data < 0;

[samples, predicted_label] = rec_completion_test(model, batch_data, mask, 0, sample_param);

% evaluate the entropy
predicted_label = mean(predicted_label,1);
temp = predicted_label(predicted_label>0);
entropy = -sum(temp .* log(temp));

author: zhirongw


text:
function x = i4vec_index_order ( n, x, indx )

%*****************************************************************************80
%
%% I4VEC_INDEX_ORDER sorts an I4VEC using an index vector.
%
%  Discussion:
%
%    The index vector itself is not modified.  Therefore, the pair
%    (X,INDX) no longer represents an index sorted vector.  If this
%    relationship is to be preserved, then simply set INDX(1:N)=(1:N).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the size of the current list.
%
%    Input, integer X(N), the list.
%
%    Input, integer INDX(N), the sort index of the list.
%
%    Output, integer X(N), the list has been sorted.
%
  y(1:n) = x(indx(1:n));
  x(1:n) = y(1:n);

  return
end


text:
function [CE,PC] = crust(P)
  % CRUST Given a set of points P determine the set of facets CE that form the
  % "crust" (assuming points are a sampling of some underlying curve/surface).
  % Implements "A New Voronoi-Based Surface Reconstruction Algorithm" [Amenta
  % et al. 1998].
  %
  % Inputs:
  %   P  #P by dim list of input points
  % Outputs:
  %   CE  #CE by dim list of facets indexing P
  %
  dim = size(P,2);
  switch dim
  case 3
    % As of MATLAB 2017a, much faster than calling tetgen
    Tri = delaunayTriangulation(P);
    [VV,VR] = Tri.voronoiDiagram;
    n = size(P,1);
    I = zeros(n,2);
    %% This loop is the bottleneck ~50% but actually scales reasonable well with
    %% input size (as expected since it is O(n) and the delaunay is O(n log n) )
    %for vi = 1:n
    %  VRi = VR{vi};
    %  VVi = VV(VRi,:);
    %  D = pdist2(VVi,P(vi,:));
    %  Ii = [0 0];
    %  [~,Ii(1)] = max(D);
    %  D = D.*sign(sum((VVi(Ii(1),:)-P(vi,:)).*(VVi-P(vi,:)),2));
    %  [~,Ii(2)] = min(D);
    %  I(vi,:) = VRi(Ii);
    %end

    % Vectorized version of for loop
    counts = cellfun(@(r) numel(r),VR);
    % Creepy way of getting indices repeated variable number of times
    J = cumsum([1;sparse(cumsum(counts),1,1)])';
    J = J(1:end-1);
    VRA = [VR{:}];
    U = VV(VRA,:) - P(J,:);
    D = normrow(U);
    S = sparse(VRA,J,D,size(VV,1),size(P,1));
    [~,I(:,1)] = max(S);
    Y = VV(I(:,1),:)-P;
    D = sign(sum(U.*Y(J,:),2)).*D;
    S = sparse(VRA,J,D,size(VV,1),size(P,1));
    [~,I(:,2)] = minnz(S);

    I = unique(I(:));
    C = VV(I,:);
    C = C(~any(isnan(C)|isinf(C),2),:);
    PC = [P;C];
    Tri = delaunayTriangulation(PC);
    FF = Tri.ConnectivityList;
    E = [ ...
      FF(:,[4 2 3]); ...
      FF(:,[3 1 4]); ...
      FF(:,[2 4 1]); ...
      FF(:,[1 3 2])];
    E = unique(sort(E,2),'rows');
    CE = E(all(E<=size(P,1),2),:);
  case 2
    % As of at least MATLAB 2017a, delaunayTriangulation is actually faster
    % than calling triangle
    %
    %[~,F] = triangle(P,[],[],'Flags','-c');
    %F = Tri.ConnectivityList;
    %[~,C] = circumradius(P,F);
    %C = C(~any(isnan(C),2),:);
    %C = unique(C,'rows');
    Tri = delaunayTriangulation(P);
    C = Tri.circumcenter;
    PC = [P;C];
    %[~,FF] = triangle(PC,[],[],'Flags','-c');
    Tri = delaunayTriangulation(PC);
    FF = Tri.ConnectivityList;
    E = edges(FF);
    CE = E(all(E<=size(P,1),2),:);
  end
end

author: alecjacobson

repo: gptoolbox


text:
function y = cvx_isconcave( x, full )
error( nargchk( 1, 2, nargin ) );
y = cvx_vexity( x );
if nargin < 2,
    y = nnz( y > 0 ) == 0;
else
    y = cvx_reshape( y <= 0, x.size_ );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function g = p03_gx ( x )

%*****************************************************************************80
%
%% P03_GX evaluates the underlying function for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(2), the point at which the function is to
%    be evaluated.
%
%    Output, real G(2), the value of the function at X.
%
  g(1) = 10000.0 * x(1) * x(2) - 1.0;
  g(2) = exp ( - x(1) ) + exp ( - x(2) ) - 1.0001;

  return
end


text:
%% Contours in OpenCV
%
% In the sample, we will:
%
% * Understand what contours are.
% * Learn about the hierarchy of contours, i.e. the parent-child relationship
%   in contours.
% * Learn to find and draw contours using the functions |cv.findContours| and
%   |cv.drawContours|.
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/d4/d73/tutorial_py_contours_begin.html>
% * <https://docs.opencv.org/3.2.0/d9/d8b/tutorial_py_contours_hierarchy.html>
%

%% Contours
%
% Contours can be explained simply as a curve joining all the continuous
% points (along the boundary), having same color or intensity. The contours
% are a useful tool for shape analysis and object detection and recognition.
%
% * For better accuracy, use binary images. So before finding contours, apply
%   threshold or canny edge detection.
% * In OpenCV, finding contours is like finding white object from black
%   background. So remember, object to be found should be white and background
%   should be black.
%
% There are three arguments in |cv.findContours| function, first one is source
% image, second is contour retrieval mode, third is contour approximation
% method. And it outputs the contours and hierarchy. The output contours are a
% cell array of all the contours in the image. Each individual contour is a
% cell array of (x,y) coordinates of boundary points of the object.
%

%% Drawing Contours
%
% To draw the contours, |cv.drawContours| function is used. It can also be
% used to draw any shape provided you have its boundary points. Its first
% argument is source image, second argument is the contours which should be
% passed as a cell array, third argument is index of contours (useful when
% drawing individual contour, or -1 to draw all contours). Remaining arguments
% are color, thickness, etc.
%
% Let's see how to find and draw contours of a binary image:
%
%   img = imread('image.jpg')
%   gray = cv.cvtColor(img, 'RGB2GRAY');
%   bw = cv.threshold(gray, 'Otsu');
%   [contours,hierarchy] = cv.findContours(bw, 'Mode','Tree', 'Method','Simple');
%   img = cv.drawContours(img, contours, 'Hierarchy',hierarchy, 'Color',[0 255 0]);
%

%% Contour Approximation Method
%
% This is the third argument in |cv.findContours| function. What does it
% denote actually?
%
% We previously explained that contours are the boundaries of a shape with the
% same intensity. It stores the (x,y) coordinates of the boundary of a shape.
% But does it store all the coordinates? That is specified by this contour
% approximation method.
%
% If you set |Method = None|, all the boundary points are stored. But actually
% do we need all the points? For example, you found the contour of a straight
% line. Do you need all the points on the line to represent that line? No, we
% need just two end points of that line. This is what |Method = Simple| does.
% It removes all redundant points and compresses the contour, thereby saving
% memory.
%
% Below image of a rectangle demonstrate this technique. Just draw a circle on
% all the coordinates in the contour array (drawn in blue color). First image
% shows points what you got with |None| (734 points) and second image shows
% the one with |Simple| (only 4 points). See how much memory it saves!
%
% <<https://docs.opencv.org/3.2.0/none.jpg>>
%

% image with a simple rectangular shape
img = zeros(200, 300, 'uint8');
img = cv.rectangle(img, [50 50 200 100], 'Color',255, 'Thickness','Filled');

% find contour w/o approximation
contours = cv.findContours(img, 'Method','None');
cntr1 = cat(1, contours{1}{:});

% find contour w/ approximation
contours = cv.findContours(img, 'Method','Simple');
cntr2 = cat(1, contours{1}{:});

% draw both contours points
out = repmat(img, 1, 1, 3);
out1 = cv.circle(out, cntr1, 3, 'Color',[255 0 0], 'Thickness','Filled');
out2 = cv.circle(out, cntr2, 3, 'Color',[255 0 0], 'Thickness','Filled');

% compare results
subplot(121), imshow(out1)
title({'Method = None', sprintf('%d points', size(cntr1,1))})
subplot(122), imshow(out2)
title({'Method = Simlpe', sprintf('%d points', size(cntr2,1))})

%% Contours Hierarchy
%
% Normally we use the |cv.findContours| function to detect objects in an
% image, right? Sometimes objects are in different locations. But in some
% cases, some shapes are inside other shapes. Just like nested figures. In
% this case, we call outer one as *parent* and inner one as *child*. This way,
% contours in an image has some relationship to each other. And we can specify
% how one contour is connected to each other, like, is it child of some other
% contour, or is it a parent etc. Representation of this relationship is
% called the *Hierarchy*.
%
% Consider an example image below :
%
% <<https://docs.opencv.org/3.2.0/hierarchy.png>>
%
% In this image, there are a few shapes numbered from *0-5*. *2* and *2a*
% denote the external and internal contours of the outermost box.
%
% Here, contours *0,1,2* are _external_ or _outermost_. We can say, they are
% in *hierarchy-0* or simply they are in _same hierarchy level_.
%
% Next comes *contour-2a*. It can be considered as a _child_ of *contour-2*
% (or in opposite way, *contour-2* is parent of *contour-2a*). So let it be in
% *hierarchy-1*. Similarly *contour-3* is child of *contour-2* and it comes in
% next hierarchy. Finally contours *4,5* are the children of *contour-3a*, and
% they come in the last hierarchy level. From the way the boxes are numbered,
% we would say *contour-4* is the first child of *contour-3a* (It can be
% *contour-5* also).
%
% We mentioned these things to understand terms like _same hierarchy level_,
% _external contour_, _child contour_, _parent contour_, _first child_, etc.
% Now let's get into OpenCV.
%

%%
% So each contour has its own information regarding what hierarchy it is, who
% is its child, who is its parent, etc. OpenCV represents it as an array of
% four values: |[Next, Previous, First_Child, Parent]|.
%
% * |Next| denotes next contour at the same hierarchical level
%
% For example, take *contour-0* in our picture. Who is next contour in its
% same level? It is *contour-1*. So simply put |Next = 1|. Similarly for
% *contour-1*, next is *contour-2*. So |Next = 2|.
%
% What about *contour-2*? There is no next contour in the same level. So
% simply, put |Next = -1|. What about *contour-4*? It is in same level with
% *contour-5*. So its next contour is *contour-5*, so |Next = 5|.
%
% * |Previous| denotes previous contour at the same hierarchical level
%
% It is same as above. Previous contour of *contour-1* is *contour-0* in the
% same level. Similarly for *contour-2*, it is *contour-1*. And for
% *contour-0*, there is no previous, so put it as |-1|.
%
% * |First_Child| denotes its first child contour
%
% There is no need of any explanation. For *contour-2*, child is *contour-2a*.
% So it gets the corresponding index value of *contour-2a*. What about
% *contour-3a*? It has two children. But we take only first child. And it is
% *contour-4*. So |First_Child = 4| for *contour-3a*.
%
% * |Parent| denotes index of its parent contour
%
% It is just opposite of |First_Child|. Both for *contour-4* and *contour-5*,
% parent contour is *contour-3a*. For *contour-3a*, it is *contour-3* and so
% on.
%
% Note: If there is no child or parent, that field is taken as |-1|.
%

%% Contour Retrieval Mode
%
% So now we know about the hierarchy style used in OpenCV, we can check into
% Contour Retrieval Modes in OpenCV with the help of same image given above.
%

%%
img1 = zeros(300, 450, 'uint8');
img1 = cv.line(img1, [30 10], [100 10], 'Color',255, 'Thickness',4);
img1 = cv.line(img1, [120 15], [430 15], 'Color',255, 'Thickness',4);
img1 = cv.rectangle(img1, [20 60 400 220], 'Color',255, 'Thickness',2);
img1 = cv.rectangle(img1, [50 100 340 150], 'Color',255, 'Thickness',2);
img1 = cv.rectangle(img1, [90 130 140 100], 'Color',255, 'Thickness','Filled');
img1 = cv.rectangle(img1, [260 155 100 50], 'Color',255, 'Thickness','Filled');
figure, imshow(img1), title('Image 1')

%%
img2 = zeros(300, 450, 'uint8');
img2 = cv.line(img2, [30 10], [100 10], 'Color',255, 'Thickness',4);
img2 = cv.line(img2, [120 15], [430 15], 'Color',255, 'Thickness',4);
img2 = cv.rectangle(img2, [20 60 400 220], 'Color',255, 'Thickness',2);
img2 = cv.rectangle(img2, [50 100 340 150], 'Color',255, 'Thickness',2);
img2 = cv.rectangle(img2, [90 130 140 100], 'Color',255, 'Thickness','Filled');
img2 = cv.rectangle(img2, [105 150 100 20], 'Color',0, 'Thickness','Filled');
img2 = cv.rectangle(img2, [105 190 100 30], 'Color',0, 'Thickness','Filled');
figure, imshow(img2), title('Image 2')

%% 1. |Mode = List|
%
% This is the simplest of the four flags (from explanation point of view). It
% simply retrieves all the contours, but doesn't create any parent-child
% relationship. Parents and kids are equal under this rule, and they are just
% contours, i.e they all belongs to same hierarchy level.
%
% So here, 3rd and 4th term in hierarchy array is always |-1|. But obviously,
% |Next| and |Previous| terms will have their corresponding values. Just check
% it yourself and verify it.
%
% Below is the result we get, and each row is hierarchy details of
% corresponding contour. For example, first row corresponds to contour 0.
% Next contour is contour 1. So |Next = 1|. There is no previous contour, so
% |Previous = -1|. And the remaining two, as told before, it is |-1|.

[~,hierarchy] = cv.findContours(img1, 'Mode','List');
hierarchy = cat(1, hierarchy{:});
display(hierarchy)

%%
% This is the good choice to use in your code, if you are not using any
% hierarchy features.
%

%% 2. |Mode = External|
%
% If you use this flag, it returns only extreme outer flags. All child
% contours are left behind. We can say, under this law, only the eldest in
% every family is taken care of. It doesn't care about other members of the
% family.
%
% So, in our image, how many extreme outer contours are there? i.e at
% *hierarchy-0* level?. Only 3, i.e contours *0,1,2*, right? Now try to find
% the contours using this flag. Here also, values given to each element is
% same as above. Compare it with above result. Below is what we get:

[~,hierarchy] = cv.findContours(img1, 'Mode','External');
hierarchy = cat(1, hierarchy{:});
display(hierarchy)

%%
% You can use this flag if you want to extract only the outer contours. It
% might be useful in some cases.
%

%% 3. |Mode = CComp|
%
% This flag retrieves all the contours and arranges them to a 2-level
% hierarchy, i.e external contours of the object (its boundary) are placed in
% *hierarchy-1*. And the contours of holes inside object (if any) is placed in
% *hierarchy-2*. If any object inside it, its contour is placed again in
% *hierarchy-1* only. And its hole in *hierarchy-2* and so on.
%
% Just consider the image of a "big white zero" on a black background. Outer
% circle of zero belongs to first hierarchy, and inner circle of zero belongs
% to second hierarchy.
%
% We can explain it with a simple image. Here the order of contours are
% labelled in red color and the hierarchy they belongs to, in green color
% (either 1 or 2). The order is same as the order OpenCV detects contours.
%
% <<https://docs.opencv.org/3.2.0/ccomp_hierarchy.png>>
%
% So consider first contour, i.e *contour-0*. It is *hierarchy-1*. It has two
% holes, contours *1,2*, and they belong to *hierarchy-2*. So for *contour-0*,
% Next contour in same hierarchy level is *contour-3*. And there is no
% previous one. And its first is child is *contour-1* in *hierarchy-2*. It has
% no parent, because it is in *hierarchy-1*. So its hierarchy array is
% |[3,-1,1,-1]|.
%
% Now take *contour-1*. It is in *hierarchy-2*. Next one in same hierarchy
% (under the parenthood of *contour-1*) is *contour-2*. No previous one. No
% child, but parent is *contour-0*. So array is |[2,-1,-1,0]|.
%
% Similarly *contour-2*. It is in *hierarchy-2*. There is not next contour in
% same hierarchy under *contour-0*. So no Next. Previous is *contour-1*. No
% child, parent is *contour-0*. So array is |[-1,1,-1,0]|.
%
% *contour-3*: Next in *hierarchy-1* is *contour-5*. Previous is *contour-0*.
% Child is *contour-4* and no parent. So array is |[5,0,4,-1]|.
%
% *contour-4*: It is in *hierarchy-2* under *contour-3* and it has no sibling.
% So no next, no previous, no child, parent is *contour-3*. So array is
% |[-1,-1,-1,3]|.
%
% Remaining you can fill up. This is the final answer we get:
%

[~,hierarchy] = cv.findContours(img2, 'Mode','CComp');
hierarchy = cat(1, hierarchy{:});
display(hierarchy)

%% 4. |Mode = Tree|
%
% It retrieves all the contours and creates a full family hierarchy list. It
% even tells, who is the grandpa, father, son, grandson and even beyond...
%
% For examle, take the above image, run it with |Mode = Tree|, reorder the
% contours as per the result given by OpenCV and analyze it. Again, red
% letters give the contour number and green letters give the hierarchy order.
%
% <<https://docs.opencv.org/3.2.0/tree_hierarchy.png>>
%
% Take *contour-0*: It is in *hierarchy-0*. Next contour in same hierarchy is
% *contour-7*. No previous contours. Child is *contour-1*. And no parent. So
% array is |[7,-1,1,-1]|.
%
% Take *contour-1*: It is in *hierarchy-1*. No contour in same level. No
% previous one. Child is *contour-2*. Parent is *contour-0*. So array is
% |[-1,-1,2,0]|.
%
% The remaining is similar. Below is the full answer:
%

[~,hierarchy] = cv.findContours(img2, 'Mode','Tree');
hierarchy = cat(1, hierarchy{:});
display(hierarchy)


text:
function [ y2, j2, f2 ] = ymdf_to_yjf_common ( y1, m1, d1, f1 )

%*****************************************************************************80
%
%% YMDF_TO_YJF_COMMON converts from YMDF to YJF form in the Common calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the YMDF date.
%
%    Output, integer Y2, J2, real F2, YJF date.
%

%
%  Check the date.
%
  [ y1, m2, d1, f1, ierror ] = ymdf_check_common ( y1, m1, d1, f1 );

  if ( ierror ~= 0 )
    y2 = 0;
    j2 = 0;
    f2 = 0.0;
    return
  end

  y2 = y1;
  j2 = d1;
  f2 = f1;
%
%  Add in the days of the elapsed months.
%
  for m = 1 : m1 - 1
    j2 = j2 + month_length_common ( y2, m );
  end

  return
end


text:
function poly_cof = r8poly_basis_1 ( ival, ntab, xtab )

%*****************************************************************************80
%
%% R8POLY_BASIS_1 computes the I-th Lagrange basis polynomial in standard form.
%
%  Discussion:
%
%    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
%    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
%    XTAB(J) for J not equal to I, and 1 when J is equal to I.
%
%    The Lagrange basis polynomials have the property that the interpolating
%    polynomial through a set of NTAB data points (XTAB,YTAB) may be
%    represented as
%
%      P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X)
%
%    Higher order interpolation at selected points may be accomplished
%    using repeated X values, and scaled derivative values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer IVAL, the index of the desired Lagrange basis polynomial.
%    IVAL should be between 1 and NTAB.
%
%    Input, integer NTAB, the number of data points XTAB.
%
%    Input, real XTAB(NTAB), the X values upon which the
%    Lagrange basis polynomial is to be based.
%
%    Output, real POLY_COF(1:NTAB), the polynomial coefficients for the IVAL-th
%    Lagrange basis polynomial.
%

%
%  Check IVAL.
%
  if ( ival < 1 | ntab < ival )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8POLY_BASE_1 - Fatal error!\n' );
    fprintf ( 1, '  IVAL must be between 1 and %d\n', ntab );
    fprintf ( 1, '  but your value is %d\n', ival );
    error ( 'R8POLY_BASE_1 - Fatal error!' );
  end
%
%  Initialize POLY_COF to the IVAL-th column of the identity matrix.
%
  poly_cof(1:ntab) = 0.0;
  poly_cof(ival) = 1.0;
%
%  Compute the divided difference table for the IVAL-th Lagrange basis
%  polynomial.
%
  poly_cof = data_to_dif ( ntab, xtab, poly_cof );
%
%  Convert the divided difference table coefficients to standard polynomial
%  coefficients.
%
  poly_cof = dif_to_r8poly ( ntab, xtab, poly_cof );

  return
end


text:
function g = tensor_creation(r, theta)

% tensor_creation - generate tensor list from polar coordinates
%
%   g = tensor_creation(r, theta)
%
%   r in [0,1] x [0,2*pi[
%
%   a = r(i)*cos(theta(i))+1)/2;
%   b = r(i)*sin(theta(i))/2;
%   g(i,:,:) = [a b; b 1-a];
%
%   Copyright (c) Gabriel Peyre

% r = r(:); theta = theta(:);
if min(r(:))<0 || max(r(:))>1 || min(theta(:))<0 || max(theta(:))>2*pi
%    warning('Tensor parameterization error');
end

g = zeros(2,2,size(theta,1),size(theta,2));

a = (r.*cos(theta)+1)/2;
b = r.*sin(theta)/2;
resh = @(x)reshape(x, [1 1 size(x,1) size(x,2)]);
a = resh(a); b = resh(b);

g(1,1,:,:) = a;
g(2,2,:,:) = 1-a;
g(1,2,:,:) = b;
g(2,1,:,:) = b;

end

author: optimaltransport


text:
function [spectrum, ntaper, freqoi] = ft_specest_irasa(dat, time, varargin)

% FT_SPECEST_IRASA separates the fractal components from the orginal power spectrum
% using Irregular-Resampling Auto-Spectral Analysis (IRASA)
%
% Use as
%   [spectrum, ntaper, freqoi] = ft_specest_irasa(dat, time, ...)
% where the input arguments are
%   dat        = matrix of chan*sample
%   time       = vector, containing time in seconds for each sample
% and the output arguments are
%   spectrum   = matrix of taper*chan*freqoi of fourier coefficients
%   ntaper     = vector containing number of tapers per element of freqoi
%   freqoi     = vector of frequencies in spectrum
%
% Optional arguments should be specified in key-value pairs and can include
%   freqoi     = vector, containing frequencies of interest
%   output     = string, indicating type of output ('fractal' or 'original', default 'fractal')
%   pad        = number, total length of data after zero padding (in seconds)
%   padtype    = string, indicating type of padding to be used, can be 'zero', 'mean', 'localmean', 'edge', or 'mirror' (default = 'zero')
%   polyorder  = number, the order of the polynomial to fitted to and removed from the data prior to the Fourier transform (default = 0, which removes the DC-component)
%   verbose    = boolean, output progress to console (0 or 1, default 1)
%
% This implements: Wen.H. & Liu.Z.(2016), Separating fractal and oscillatory components in the power
% spectrum of neurophysiological signal. Brain Topogr. 29(1):13-26. The source code accompanying the
% original paper is avaible from https://purr.purdue.edu/publications/1987/1
%
% For more information about the difference between the current and previous version and how to use this
% function, please see https://www.fieldtriptoolbox.org/example/irasa/
%
% See also FT_FREQANALYSIS, FT_SPECEST_MTMFFT, FT_SPECEST_MTMCONVOL, FT_SPECEST_TFR, FT_SPECEST_HILBERT, FT_SPECEST_WAVELET

% Copyright (C) 2019-2020, Rui Liu, Arjen Stolk
%
% 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$

% these are for speeding up computation of tapers on subsequent calls
persistent previous_argin previous_tap

% get the optional input arguments
freqoi    = ft_getopt(varargin, 'freqoi', 'all');
output    = ft_getopt(varargin, 'output', 'fractal');
pad       = ft_getopt(varargin, 'pad');
padtype   = ft_getopt(varargin, 'padtype', 'zero');
polyorder = ft_getopt(varargin, 'polyorder', 1);
verbose   = ft_getopt(varargin, 'verbose', true);
fbopt     = ft_getopt(varargin, 'feedback');

% check output option
if ~strcmp(output, {'fractal','original'})
  ft_error('The current version ft_specest_irasa outputs ''fractal'' or ''original'' power only. For more information about the update, see https://www.fieldtriptoolbox.org/example/irasa/');
end

% this does not work on integer data
dat = cast(dat, 'double');
if ~isa(dat, 'double') && ~isa(dat, 'single')
  dat = cast(dat, 'double');
end

% dat param
[nchan,ndatsample] = size(dat);
fsample = 1./mean(diff(time));

% subset param
subset_nsample = 2^floor(log2(ndatsample*0.9)); % the number of sub-segments in samples is the power of 2 that does not exceed 90% of input data
subset_num     = 10;                            % the number of sub-segments
subset_dist    = floor((ndatsample  - subset_nsample)/(subset_num - 1)); % distance between sub-segements in sample, to evenly distribute the sub-subsegments within the total length of input data

% resampling ratio
hset  = 1.1:0.05:1.9;
nhset = length(hset);

% remove polynomial fit from the data -> default is demeaning
if polyorder >= 0
  dat = ft_preproc_polyremoval(dat, polyorder, 1, ndatsample);
end

% zero padding
if isempty(pad)
  % if no padding is specified this is set equal to the current data length
  pad = subset_nsample/fsample;
end
if round(pad * fsample) < subset_nsample
  ft_error('the padding that you specified is shorter than the data');
end
endnsample = round(pad * fsample);  % total number of samples of padded data
endtime    = pad;                   % total time in seconds of padded data

% set freqboi and freqoi
freqoiinput = freqoi;
if isnumeric(freqoi) % if input is a vector
  freqboi   = round(freqoi ./ (fsample ./ endnsample)) + 1; % is equivalent to: round(freqoi .* endtime) + 1;
  freqboi   = unique(freqboi);
  freqoi    = (freqboi-1) ./ endtime; % boi - 1 because 0 Hz is included in fourier output
elseif strcmp(freqoi, 'all') % if input was 'all'
  freqboilim = round([0 fsample/2] ./ (fsample ./ endnsample)) + 1;
  freqboi    = freqboilim(1):1:freqboilim(2);
  freqoi     = (freqboi-1) ./ endtime;
end
nfreqboi = length(freqboi);
nfreqoi  = length(freqoi);

% throw a warning if input freqoi is different from output freqoi
if isnumeric(freqoiinput)
  % check whether padding is appropriate for the requested frequency resolution
  rayl = 1/endtime;
  if any(rem(freqoiinput,rayl)) % not always the case when they mismatch
    ft_warning('padding not sufficient for requested frequency resolution, for more information please see the FAQs on www.ru.nl/neuroimaging/fieldtrip');
  end
  if numel(freqoiinput) ~= numel(freqoi) % freqoi will not contain double frequency bins when requested
    ft_warning('output frequencies are different from input frequencies, multiples of the same bin were requested but not given');
  else
    if any(abs(freqoiinput-freqoi) >= eps*1e6)
      ft_warning('output frequencies are different from input frequencies');
    end
  end
end

% determine whether tapers need to be recomputed
current_argin = {output, time, endnsample, freqoi}; % reasoning: if time and endnsample are equal, it's the same length trial, if the rest is equal then the requested output is equal
if isequal(current_argin, previous_argin)
  % don't recompute tapers
  tap = previous_tap;
else
  if strcmp(output, 'fractal')
    % (re)compute tapers, 1:mid are upsample tapers, mid+1:end are downsample tapers
    tap = cell(nhset*2,1);
    for ih = 1:nhset
      [n, d] = rat(hset(ih)); % n > d
      tmp = hanning(size(resample(zeros(subset_nsample,nchan), n, d),1))';
      tap{ih,1} = tmp./norm(tmp, 'fro');% for upsampled subsets
      tmp = hanning(size(resample(zeros(subset_nsample,nchan), d, n),1))';
      tap{ih+nhset,1} = tmp./norm(tmp, 'fro');% for downsampled subsets
    end
  elseif strcmp(output, 'original')
    tap = hanning(subset_nsample)';
    tap = tap./norm(tap, 'fro');
  end
end

% set ntaper
if strcmp(output, 'fractal')
  ntaper = repmat(size(tap,2),nfreqoi,1); % pretend there's only one taper
elseif strcmp(output, 'original')
  ntaper = repmat(size(tap,1),nfreqoi,1);
end

% feedback of computing progress
if isempty(fbopt)
  fbopt.i = 1;
  fbopt.n = 1;
end
str = sprintf('nfft: %d samples, datalength: %d samples, %d tapers',endnsample,ndatsample,ntaper(1));
st  = dbstack;
if length(st)>1 && strcmp(st(2).name, 'ft_freqanalysis')
  % specest_mtmfft has been called by ft_freqanalysis, meaning that ft_progress has been initialised
  ft_progress(fbopt.i./fbopt.n, ['processing trial %d/%d ',str,'\n'], fbopt.i, fbopt.n);
elseif verbose
  fprintf([str, '\n']);
end

% compute irasa or fft
spectrum = cell(ntaper(1),1);
for itap = 1:ntaper(1)
  %%%%%%%% IRASA %%%%%%%%%%
  if strcmp(output,'fractal')
    pow  = zeros(nchan,nfreqboi,nhset);
    for ih = 1:nhset % loop resampling ratios hset
      upow = zeros(nchan,nfreqboi);
      dpow = zeros(nchan,nfreqboi);
      [n, d] = rat(hset(ih)); % n > d
      for k = 0 : subset_num-1 % loop #subset
        % pair-wised resampling
        subset_start = subset_dist*k + 1;
        subset_end = subset_start+subset_nsample-1;
        subdat = dat(:,subset_start:subset_end);
        udat = resample(subdat', n, d)'; % upsample
        ddat = resample(subdat', d, n)'; % downsample

        % compute auto-power of resampled data
        % upsampled
        postpad = ceil(round(pad * fsample) - size(udat,2));
        tmp = fft(ft_preproc_padding(bsxfun(@times,udat,tap{ih,1}), padtype, 0, postpad), endnsample, 2);
        tmp = tmp(:,freqboi);
        tmp = tmp .* sqrt(2 ./ endnsample);
        tmp = abs(tmp).^2;
        upow = upow + tmp;
        % downsampled
        postpad = ceil(round(pad * fsample) - size(ddat,2));
        tmp = fft(ft_preproc_padding(bsxfun(@times,ddat,tap{ih+nhset,1}), padtype, 0, postpad), endnsample, 2);
        tmp = tmp(:,freqboi);
        tmp = tmp .* sqrt(2 ./ endnsample);
        tmp = abs(tmp).^2;
        dpow = dpow + tmp;
      end % loop #subset

      % average across subsets for noise filtering
      upow = upow / subset_num;
      dpow = dpow / subset_num;

      % geometric mean for relocating oscillatory component
      pow(:,:,ih) = sqrt(upow.*dpow);
    end % loop resampling ratios hset

    % median across resampling factors
    spectrum{itap} = median(pow,3);

    %%%%%%%% FFT %%%%%%%%%%
  elseif strcmp(output,'original')
    pow  = zeros(nchan,nfreqboi);
    for k = 0 : subset_num-1 % loop #subset
      % pair-wised resampling
      subset_start = subset_dist*k + 1;
      subset_end = subset_start+subset_nsample-1;
      subdat = dat(:,subset_start:subset_end);
      % compute auto-power
      postpad = ceil(round(pad * fsample) - size(subdat,2));
      tmp = fft(ft_preproc_padding(bsxfun(@times,subdat,tap(itap,:)), padtype, 0, postpad), endnsample, 2);
      tmp = tmp(:,freqboi);
      tmp = tmp .* sqrt(2 ./ endnsample);
      tmp = abs(tmp).^2;
      pow = pow + tmp;
    end % loop #subset

    % average across subsets for noise filtering
    spectrum{itap} = pow / subset_num;

  end % if output
end % for taper

spectrum = reshape(vertcat(spectrum{:}),[nchan ntaper(1) nfreqboi]); % collecting in a cell-array and later reshaping provides significant speedups
spectrum = permute(spectrum, [2 1 3]);

% remember the current input arguments, so that they can be
% reused on a subsequent call in case the same input argument is given
previous_argin = current_argin;
previous_tap   = tap;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function TCnew = align(TC,varargin)

%ALIGN Align and trim traces in threecomp object.
% TC = ALIGN(TC) resample traces so that samples from each component align
% in time. As needed, this may include: trimming the start or end times;
% changing the sample rate; resampling traces to eliminate subsample
% offsets. If start and end times do not overlap for all traces, the
% waveforms returned will be of zero length.
%
% TC = ALIGN(TC,ALIGNFREQ) specifies an explicit frequency for the
% resampling. If no frequency is given, then traces are resampled to same
% frequency as the Z component in each threecomp object.
%
% If a trace contains fewer than 25 samples, that threecomp element will
% not be aligned and the resulting traces will be of zero length (empty).

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

GETFREQ = 1;
if ~isempty(varargin)
   alignFreq = varargin{1};
   if ~isreal(alignFreq)
      error('Threecomp:align:resampling rate must be a real number');
   end
   GETFREQ = 0;
end

TCnew = TC;
for n = 1:length(TC)
    if GETFREQ
       alignFreq = get(TC(n).traces(1),'FREQ');
    end
    TCnew(n) = do_one(TC(n),alignFreq);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Process a single threecomp trace tuple
function TC = do_one(TC,alignFreq)

% CHECK START/END TIMES AND FREQUENCIES
SAMETIMES = 1;
SAMEFREQUENCY = 1;
SUFFICIENTLENGTH = 1;
startTime = get(TC.traces,'START_MATLAB');
endTime = get(TC.traces,'END_MATLAB');
freq = get(TC.traces,'FREQ');
dataLength = get(TC.traces,'DATA_LENGTH');
if ( startTime(2)~=startTime(1) ) || ( startTime(3)~=startTime(1) )
    SAMETIMES = 0;
end
if ( endTime(2)~=endTime(1) ) || ( endTime(3)~=endTime(1) )
    SAMETIMES = 0;
end
if ( freq(2)~=freq(1) ) || ( freq(3)~=freq(1) )
    SAMEFREQUENCY = 0;
end
if min(dataLength)<25
   SUFFICIENTLENGTH = 0;
end

if ~SUFFICIENTLENGTH
    TC.traces = set(TC.traces,'DATA',[]);
else

    % ALIGN IF NECESSARY
    if ~SAMETIMES || ~SAMEFREQUENCY
        TC.traces = align(TC.traces , startTime , alignFreq );
    end

    % CROP TRACES IF NECESSARY
    if ~SAMETIMES || ~SAMEFREQUENCY
        samples = floor(get(TC.traces(1),'FREQ')*86400*(min(endTime) - max(startTime)));
        TC.traces = extract(TC.traces, 'TIME&SAMPLES', max(startTime), samples);
    end

end

author: geoscience-community-codes

repo: GISMO


text:
function basis_11_t4_test ( )

%*****************************************************************************80
%
%% BASIS_11_T4_TEST verifies BASIS_11_T4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  node_num = 4;

  t = [ ...
    2.0, 0.0; ...
    4.0, 3.0; ...
    0.0, 4.0; ...
    0.0, 0.0 ]';
%
%  The node associated with the fourth basis function is the centroid.
%
  t(1,4) = sum ( t(1,1:3) ) / 3.0;
  t(2,4) = sum ( t(2,1:3) ) / 3.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BASIS_11_T4_TEST:\n' );
  fprintf ( 1, '  Verify basis functions for element T4.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Physical Nodes:\n' );
  fprintf ( 1, '\n' );
  for j = 1 : node_num
    fprintf ( 1, '  %8d  %7f  %7f\n', j, t(1:2,j) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The basis function values at basis nodes\n' );
  fprintf ( 1, '  should form the identity matrix.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : node_num
    for j = 1 : node_num
      [ phi(i,j), dphidx(i,j), dphidy(i,j) ] = basis_11_t4 ( t, i, t(1:2,j) );
    end
  end

  for i = 1 : node_num
    fprintf ( 1, '  %7f  %7f  %7f  %7f\n', phi(i,1:node_num) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The X and Y derivatives should sum to 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      dPhidX sum    dPhidY sum\n' );
  fprintf ( 1, '\n' );

  for j = 1 : node_num

    sum_x = sum ( dphidx(1:node_num,j) );
    sum_y = sum ( dphidy(1:node_num,j) );
    fprintf ( 1, '  %14e  %14e\n', sum_x, sum_y );

  end

  return
end


text:
function [c] = tricircumcenter3d(node,elem)

% http://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html

v1 = node(elem(:,1),:);
v2 = node(elem(:,2),:);
v3 = node(elem(:,3),:);

ba = v2 - v1;
ca = v3 - v1;

baLength = sum(ba.^2,2);
caLength = sum(ca.^2,2);

crossbc = cross(ba,ca,2);

d = 0.5 ./sum(crossbc.^2,2);

xc = ((baLength.*ca(:,2) - caLength.*ba(:,2)).*crossbc(:,3) - (baLength.*ca(:,3) - caLength.*ba(:,3)).*crossbc(:,2)).*d;
yc = ((baLength.*ca(:,3) - caLength.*ba(:,3)).*crossbc(:,1) - (baLength.*ca(:,1) - caLength.*ba(:,1)).*crossbc(:,3)).*d;
zc = ((baLength.*ca(:,1) - caLength.*ba(:,1)).*crossbc(:,2) - (baLength.*ca(:,2) - caLength.*ba(:,2)).*crossbc(:,1)).*d;

c = [xc, yc, zc];

c = v1 + c;

repo: ifem


text:
%% <../index.html MATLAB xUnit Test Framework>: How RUNTESTS Searches for Test Cases
% When you call |runtests| with no input arguments:
%
%   >> runtests
%
% it automatically searches for all the test cases in the current directory.  It
% looks for test cases in three types of M-files:
%
% 1. An M-file function whose name begins or ends with "test" or "Test" and that does
% not return an output argument.  Such a function is considered to be a single
% test case.
%
% 2. An M-file function whose name begins or ends with "test" or "Test" and that returns
% an output argument that is a test suite.  Such a function is considered to contain
% subfunction-style test cases.  Each subfunction whose name begins or ends with "test"
% or "Test" is a test case.
%
% 3. An M-file that defines a subclass of TestCase.  Each method beginning or ending with
% "test" or "Test" is a test case.
%
% |runtests| uses the |TestSuite| static methods |fromName| and |fromPwd| to
% automatically construct the test suites.
%
% Here are a couple of examples.
%
% |TestSuite.fromName| takes an M-file name, determines what
% kind of test file it is, and returns a cell array of test case objects.

cd examples_general
test_suite_1 = TestSuite.fromName('testSetupExample')

%%
% |TestSuite.fromPwd| returns a test suite based on all the test files in the
% current directory.

test_suite_2 = TestSuite.fromPwd()

%%
% <../index.html Back to MATLAB xUnit Test Framework>

%%
% Copyright 2008-2010 The MathWorks, Inc.

author: ricedsp


text:
function [updateField, updateCertainty] = ...
    morphon3d(moving, movingCertainty, fixed, fixedCertainty, ...
    scale, varargin)
% MORPHON3D Estimates a displacement field between moving and fixed using the Morphon
%
% [updateField, updateCertainty] = morphon3d(...
%     moving, movingCertainty, fixed, fixedCertainty, scale)
%
% INPUT ARGUMENTS
% moving                - Moving image
% movingCertainty       - Certainty mask of moving image
% fixed                 - Fixed image
% fixedCertainty        - Certainty mask of fixed image
% scale                 - Current scale
%
% OPTIONAL INPUT ARGUMENTS
% 'solution'            - Solution to apply when estimating the
%                         displacement field (3)
% 'applyCertainty'      - Apply certainty when smoothing tensor elements
%                         and equation system elements (true)
% 'tensorSigma'         - Sigma of Gaussian to smooth tensor (2.5)
% 'elementsSigma'       - Sigma of Gaussian to smooth equation system (2.5)
%
% OUTPUT ARGUMENTS
% updateField 			- Estimated displacement
% updateCertainty 		- Certainty of displacement

% Copyright (c) 2011 Daniel Forsberg
% danne.forsberg@outlook.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 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/>.

% Setup default parameters
solution = 3;
applyCertainty = true;
tensorSigma = 2.5;
elementsSigma = 2.5;

for k=1:2:length(varargin),         % overwrites default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

%***************************************************************
% Load quadrature filters
%***************************************************************

load quadratureFiltersForMorphonRegistration3D

numberOfFilterDirections = 6;

%***************************************************************
% Determine image dimension
%***************************************************************

sizeImage = size(fixed);

%***************************************************************
% Initialize some variables
%***************************************************************

dk = zeros([sizeImage numberOfFilterDirections]);
ck = zeros([sizeImage numberOfFilterDirections]);
t11 = zeros(sizeImage);
t12 = zeros(sizeImage);
t13 = zeros(sizeImage);
t22 = zeros(sizeImage);
t23 = zeros(sizeImage);
t33 = zeros(sizeImage);
updateField = cell(3,1);

%***************************************************************
% Compute ck, dk and T
%***************************************************************

for k = 1 : numberOfFilterDirections
    % Filter moving and fixed
    qa = imfilter(moving, f{k}, 'same', 'conv');
    qb = imfilter(fixed, f{k}, 'same', 'conv');

    % Estimate dk and ck
    qq = qa.*conj(qb);
    Aqq = abs(qq);
    dk(:,:,:,k) = atan2(imag(qq),real(qq));
    ck(:,:,:,k) = Aqq;
    ck(:,:,:,k) = ck(:,:,:,k).^0.5.*(cos(dk(:,:,:,k)/2).^2);

    % Estimate T (Ta and Tb)
    Aqb = abs(qb);
    t11 = t11 + Aqb*m11{k};
    t12 = t12 + Aqb*m12{k};
    t13 = t13 + Aqb*m13{k};
    t22 = t22 + Aqb*m22{k};
    t23 = t23 + Aqb*m23{k};
    t33 = t33 + Aqb*m33{k};
end

% Add certainty masks
if ~isempty(fixedCertainty)
    ck = bsxfun(@times,ck,fixedCertainty);
end
if ~isempty(movingCertainty)
    ck = bsxfun(@times,ck,movingCertainty);
end

%***************************************************************
% Average and normalize Ta and Tb
%***************************************************************

if applyCertainty
    Tcert = sqrt(t11.^2 + 2*t12.^2 + 2*t13.^2 + t22.^2 + 2*t23.^2 + t33.^2);
    t11 = normgauss_smoothing(t11, Tcert, tensorSigma);
    t12 = normgauss_smoothing(t12, Tcert, tensorSigma);
    t13 = normgauss_smoothing(t13, Tcert, tensorSigma);
    t22 = normgauss_smoothing(t22, Tcert, tensorSigma);
    t23 = normgauss_smoothing(t23, Tcert, tensorSigma);
    t33 = normgauss_smoothing(t33, Tcert, tensorSigma);
else
    t11 = gauss_smoothing(t11, tensorSigma);
    t12 = gauss_smoothing(t12, tensorSigma);
    t13 = gauss_smoothing(t13, tensorSigma);
    t22 = gauss_smoothing(t22, tensorSigma);
    t23 = gauss_smoothing(t23, tensorSigma);
    t33 = gauss_smoothing(t33, tensorSigma);
end

tNorm = (sqrt(t11.^2 + 2*t12.^2 + 2*t13.^2 + t22.^2 + 2*t23.^2 + t33.^2) + eps);
t11 = t11./max(tNorm(:));
t12 = t12./max(tNorm(:));
t13 = t13./max(tNorm(:));
t22 = t22./max(tNorm(:));
t23 = t23./max(tNorm(:));
t33 = t33./max(tNorm(:));

%***************************************************************
% Solutions based upon minimizing a least square problem,
% i.e. runParameters.solution = 0-5
%***************************************************************

if (solution >= 0 && solution <= 5)

    [a11, a12, a13, a22, a23, a33, b1, b2, b3] = ...
        estimate_displacement_based_upon_least_square_solution3d(...
            solution, 1, ...
            sizeImage, ...
            dk, ck, ...
            filterDirection, ...
            numberOfFilterDirections, ...
            t11, t12, t13, t22, t23, t33);

    %***********************************************************
    % Reqularize axx and bx
    %***********************************************************

    a11 = gauss_smoothing(a11, elementsSigma);
    a12 = gauss_smoothing(a12, elementsSigma);
    a13 = gauss_smoothing(a13, elementsSigma);
    a22 = gauss_smoothing(a22, elementsSigma);
    a23 = gauss_smoothing(a23, elementsSigma);
    a33 = gauss_smoothing(a33, elementsSigma);

    b1 = gauss_smoothing(b1, elementsSigma);
    b2 = gauss_smoothing(b2, elementsSigma);
    b3 = gauss_smoothing(b3, elementsSigma);

    %***********************************************************
    % Solve equation system and compute incremental displacement and
    % incremental certainty
    %***********************************************************

    updateCertainty = a11 + a22 + a33;

    norm = 1./(a11.*a22.*a33 - a11.*a23.*a23 - a12.*a12.*a33 + a12.*a23.*a13 + a13.*a12.*a23 - a13.*a22.*a13 + eps);

    updateField{2} = norm.*((b3.*(a12.*a23 - a13.*a22)) - (b2.*(a12.*a33 - a13.*a23)) + (b1.*(a22.*a33 - a23.*a23)));
    updateField{1} = norm.*((b2.*(a11.*a33 - a13.*a13)) - (b3.*(a11.*a23 - a13.*a12)) - (b1.*(a12.*a33 - a23.*a13)));
    updateField{3} = norm.*((b3.*(a11.*a22 - a12.*a12)) - (b2.*(a11.*a23 - a12.*a13)) + (b1.*(a12.*a23 - a22.*a13)));

end

%***************************************************************
% Step size compensation
%***************************************************************

% Scale independent ampfactor
displacementAmplification = 2;

for k = 1 : 3
    updateField{k} = ...
        -displacementAmplification*updateField{k};
end

% Fix to gain more power on fine detail scales
updateCertainty = updateCertainty.*2^-scale;

author: fordanic

repo: image-registration


text:
function [data,header] = ipcc2mat(fname,r)
% PURPOSE: Import IPCC climate observation data into Matlab (3D) matrix
% -------------------------------------------------------------------
% USAGE: [data,header] = ipcc2mat(fname,r)
% where: fname = filename of ipcc climate observation file in ascii .dat
%                format as provided at http://www.ipcc-data.org
%        r = optional flag to rotate& flip data matrix for pretty printing in
%            mat format (1) or not (0). Default is 1 (rotate),
% Output:
%        data = 3D matrix in format (rows/cols/month) usually (360,720,12)
%
% Note:
%
% Example: [data,header]=ipcc2mat('ctmp8190.dat',1);
%          figure, imagesc(data)
%
% Felix Hebeler, Geography Dept., University Zurich, May 2007.

if ~exist('r','var')
    r=1;
end
nodata = -9999;

%% open file
fid = fopen(fname);
if fid==-1
  error('File not found or permission denied');
end

%% read header ifno and data
h1= fgetl(fid); % this should be text
h2= fgetl(fid); % this the info
data = [];
data = [data; fscanf(fid, '%f')];
fclose(fid);

%% make a nice header
strs = {'grd_sz','xmin','ymin','xmax','ymax','n_cols','n_rows','n_months'};
header={};
for i = 1:length(strs)
    [t,h2]=strtok(h2);
    eval(['header.',strs{i},'=str2num(t);'])
end

%% replace NaN and reshape data
data(data==nodata)=NaN;
data=reshape(data,header.n_cols,header.n_rows,header.n_months);

%% rotate (only works for 2D matrices)
if r==1
    t=nan(header.n_rows,header.n_cols,header.n_months);
    for d = 1:size(data,3)
        t(:,:,d)=fliplr(rot90(data(:,:,d),3));
    end
end
data=t;
clear t;

author: Sable

repo: mcbench-benchmarks


text:
function [kappa] = kappaindex(X,G,N)
%  USAGE kappa=kappaindex(X,G,N)
%
%  Computing Kappa Index for Clustering Agreement by Hossein Mobahi
%
% X is a vector of length M (number of data samples).
%   Each entry of X is associated with the cluster index for that sample
% G is a vector in the same format of X belonging to ground truth indecies
%   or a secondary clustering
% N is the number of clusters
%
% This measure only works when X and G have the same number of clusters
% In addition, the set of indecies must be integers from one through
% the number of clusters.
% However, they might be equivalent through a permutation of indecies
%
% Example:
% >> X=[1 1 2 2 3 3]; G=[2 2 1 1 3 3]; kappaindex(X,G,3)
% ans = 1
% >> X=[1 1 2 2 4 4 3]; G=[4 2 2 2 1 1 3]; kappaindex(X,G,4)
% ans = 0.8056
%
%

% Check all permutations of cluster indecies
% And pick the one with maximal agreement

P=perms(1:N);

maxn=-inf;
for i=1:size(P,1)
 n=0;
 for j=1:length(X)
  if G(j)==P(i,X(j))
   n=n+1;
  end
 end
 if n>maxn
  maxn=n;
  maxi=i;
 end
end

% Build matrix
A=zeros(N,N);
for j=1:length(X)
  x=P(maxi,X(j));
  g=G(j);
  A(x,g)=A(x,g)+1;
end

% Compute the kappa index
t=0;
for i=1:N
 t=t+sum(A(i,:))*sum(A(:,i));
end
kappa=(length(X)*trace(A)-t)/(length(X)*length(X)-t);


text:
function R= VBA_spm_autocorr(y)
% computes simple autocorrelation function of signal y
% function R = spm_autocorr(y)
% Note: the autocorrelation is computed from lag 0 to lag t (length of the
% times series), and back. The symmetry of the autocorrelation function can
% be eyeballed using the following command:
% figure,plot(fftshift(R)')
% This way, the 0-lag correlation (always 1) is in the centre of the axes.
% IN:
%   - y: nXt time series matrix
% OUT:
%   - R: nXt autocorrelation matrix

[n,t] = size(y);
R = zeros(n,t*2);
% standardize y
my = VBA_nanmean(y,2);
sy = sqrt(VBA_nanvar(y,2));
sy(sy==0) = 1; % correct for identically constant time series
y = y - repmat(my,1,t);
y = diag(1./sy)*y;
% compute auto-correlation using FFT
hasNan = any(any(isnan(y)));
if hasNan
     disp('= Warning: Autocorrelation with missing data could yield unreliable results');
end
for i=1:n
    if hasNan
        fr = nanfft(y(i,:),t*2);
    else
        fr = fft(y(i,:),t*2);
    end
    tfr = fr';
    S = fr(:).*tfr(:);
    R(i,:) = ifft(S)'/t;

end

end

author: MBB-team

repo: VBA-toolbox


text:
function geometry_test172 ( )

%*****************************************************************************80
%
%% TEST172 tests QUAT_CONJ, QUAT_INV, QUAT_MUL, QUAT_NORM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
  q1 = [ 2.0, 3.0, 4.0, 5.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST172\n' );
  fprintf ( 1, '  QUAT_CONJ conjugates a quaternion;\n' );
  fprintf ( 1, '  QUAT_INV inverts a quaternion;\n' );
  fprintf ( 1, '  QUAT_MUL multiplies quaternions.\n' );
  fprintf ( 1, '  QUAT_NORM computes the norm.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Q1 =            %10f  %10f  %10f  %10f\n', q1(1:4) );

  fprintf ( 1, '  Norm ( Q1 ) = %10f\n', quat_norm ( q1 ) );

  q2 = quat_conj ( q1 );

  fprintf ( 1, '  Q2 = conj(Q1) = %10f  %10f  %10f  %10f\n', q2(1:4) );

  q3 = quat_mul ( q1, q2 );

  fprintf ( 1, '  Q3 = Q1*Q2 =    %10f  %10f  %10f  %10f\n', q3(1:4) );

  q2 = quat_inv ( q1 );

  fprintf ( 1, '  Q2 = inv(Q1) =  %10f  %10f  %10f  %10f\n', q2(1:4) );

  q3 = quat_mul ( q1, q2 );

  fprintf ( 1, '  Q3 = Q1*Q2 =    %10f  %10f  %10f  %10f\n', q3(1:4) );

  return
end


text:
function varargout = spm_bsplinc(varargin)
% Return volume of B-spline coefficients
% FORMAT c = spm_bsplinc(V,d)
% 	V - a structure of image volume information (from spm_vol)
% 	d(1:3) - degree of B-spline (from 0 to 7) along different dimensions
%       d(4:6) - 1/0 to indicate wrapping along the dimensions
% 	c - returned volume of B-spline coefficients
%
% This function deconvolves B-splines from volume V, returning
% coefficients, c.  These coefficients are then passed to spm_bsplins
% in order to sample the data using B-spline interpolation.
%
%_______________________________________________________________________
%
% References:
%	M. Unser, A. Aldroubi and M. Eden.
%	"B-Spline Signal Processing: Part I-Theory,"
%	IEEE Transactions on Signal Processing 41(2):821-832 (1993).
%
%	M. Unser, A. Aldroubi and M. Eden.
%	"B-Spline Signal Processing: Part II-Efficient Design and
%	Applications,"
%	IEEE Transactions on Signal Processing 41(2):834-848 (1993).
%
%	M. Unser.
%	"Splines: A Perfect Fit for Signal and Image Processing,"
%	IEEE Signal Processing Magazine, 16(6):22-38 (1999)
%
%	P. Thevenaz and T. Blu and M. Unser.
%	"Interpolation Revisited"
%	IEEE Transactions on Medical Imaging 19(7):739-758 (2000).
%_______________________________________________________________________
% @(#)spm_bsplinc.m	2.3 John Ashburner 02/07/30

%-This is merely the help file for the compiled routine
error('spm_bsplinc.c not compiled.');


text:
function m = ema(x,N)
L = length(x);
m = zeros(L,1);
w = 2/(N+1);
m(1) = x(1);
for i=2:L
    m(i) = w*(x(i)-m(i-1)) + m(i-1);
end


text:
function [M]=IpaS(d, a)
%A special bidiagonal matrix in the QTT-format
%   [M]=IPAS(D, A)
%   Generates I+a*S_{-1} matrix in the QTT-format:
%   1 0 0 0
%   a 1 0 0
%   0 a 1 0
%   0 0 a 1
% Convenient for Crank-Nicolson and time gradient matrices
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

M = cell(d,1);
M{1} = zeros(2,2,0);
M{1}(:,:,1)=[1,0;a,1];
if (d>1)
    M{1}(:,:,2)=[0,a;0,0];
end;
for i=2:d-1
    M{i}=zeros(2,2,2,2);
    M{i}(:,:,1,1)=eye(2);
    M{i}(:,:,2,1)=[0,0;1,0];
    M{i}(:,:,2,2)=[0,1;0,0];
end;
if (d>1)
    M{d} = zeros(2,2,2);
    M{d}(:,:,1)=eye(2);
    M{d}(:,:,2)=[0,0;1,0];
end;
M=tt_matrix(M); %Bydlocode
end

author: oseledets

repo: TT-Toolbox


text:
classdef (Abstract) audio < matlab.mixin.Copyable
%AUDIO Abstract superclass providing audio-related properties and methods.
%
%   The IOSR.DSP.AUDIO superclass is an abstract class that defines
%   audio-related properties and methods. As an abstract class, it cannot
%   be instantiated.
%
%   IOSR.DSP.AUDIO is a subclass of MATLAB.MIXIN.COPYABLE, which is a
%   subclass of the HANDLE class. Hence subclasses of AUDIO are handle
%   classes.
%
%   IOSR.DSP.AUDIO properties:
%       filename        - Name of the corresponding audio file
%       fs              - Sampling frequency (Hz)
%       signal          - The sampled data (read-only)
%
%   IOSR.DSP.AUDIO methods:
%       sound           - Replay the audio signal
%     Static methods:
%       up_down_mix     - Remix a source to a designated channel count
%       normalize       - Normalize an audio signal
%       spat            - Spatialise a sound using a SOFA file
%
%   See also IOSR.BSS.MIXTURE, IOSR.BSS.SOURCE, HANDLE,
%            MATLAB.MIXIN.COPYABLE.

%   Copyright 2016 University of Surrey.

    properties (AbortSet)
        filename            % Name of the audio file
        fs = 16000          % Sampling frequency (Hz)
    end

    properties (Abstract, Dependent, SetAccess = protected)
        signal              % The sampled data (read-only)
    end

    properties (Access = protected)
        rendered = false    % Flag indicates whether files match object
    end

    methods % public methods

        function obj = audio(filename,fs)
        %AUDIO Create the audio object

            if nargin > 0
                obj.filename = filename;
            end
            if nargin > 1
                obj.fs = fs;
            end

        end

        function sound(obj)
        %SOUND Replay the audio signal
        %
        %   IOSR.DSP.AUDIO.SOUND() replays the audio signal.

            obj.replay(obj.signal)

        end

        % validate properties

        % validate filename
        function set.filename(obj,val)
            assert(ischar(val) || isempty(val), 'iosr:audio:invalidFile', 'FILENAME must be a char array or an empty array')
            obj.filename = val;
            obj.property_changed('filename',val);
        end

        % validate fs
        function set.fs(obj,val)
            assert(isscalar(val), 'iosr:audio:invalidFs', 'FS must be a scalar')
            assert(val > 0, 'iosr:audio:invalidFs', 'FS must be greater than 0')
            obj.fs = val;
            obj.property_changed('fs',val);
        end

    end

    methods (Access = protected)

        function property_changed(obj,name,val) %#ok<INUSD>
        %PROPERTY_CHANGED handle property changes

            obj.rendered = false;
        end

        function replay(obj,x)
        %REPLAY Replay any audio signal

            sound(obj.normalize(x),obj.fs) % replay

        end

        function set_properties(obj,vgin,propNames)
        %SET_PROPERTIES Parse varargin input to set object properties
        %
        %   IOSR.DSP.AUDIO.SET_PROPERTIES(VGIN,PROPNAMES) searches through
        %   the cell array VGIN (a cell array of parameter/value pairs) for
        %   properties specified in the cell array PROPNAMES. Any matching
        %   properties are written to OBJ.

            % count arguments
            nArgs = length(vgin);
            if round(nArgs/2)~=nArgs/2
               error('iosr:audio:nameValuePair',[class(obj) ' needs propertyName/propertyValue pairs'])
            end
            % overwrite defults
            for pair = reshape(vgin,2,[]) % pair is {propName;propValue}
               IX = strcmpi(pair{1},propNames); % find match parameter names
               if any(IX)
                  % do the overwrite
                  obj.(propNames{IX}) = pair{2};
               else
                  error('iosr:audio:unknownOption','%s is not a recognized parameter name',pair{1})
               end
            end

        end

        function cpObj = copyElement(obj)
        %COPYELEMENT Overload copy method with additional functionality

            % Make a shallow copy of all properties
            cpObj = copyElement@matlab.mixin.Copyable(obj);

        end

    end

    methods (Static)

        function [y,gain] = normalize(x)
        %NORMALIZE Normalize an audio signal
        %
        %   Y = IOSR.DSP.AUDIO.NORMALIZE(X) normalises X such that
        %   ABS(Y(:))<1.
        %
        %   [Y,GAIN] = IOSR.DSP.AUDIO.NORMALIZE(X) returns the scalar gain
        %   GAIN applied to X in order to normalize it.

            gain = (1-((2*16)/(2^16)))/max(abs(x(:)));
            y = x.*gain;

        end

        function y = spat(SOFAfile,x,azimuth,elevation,fs)
        %SPAT Spatialise a sound using a SOFA file
        %
        %   IOSR.DSP.AUDIO.SPAT(SOFAFILE,X,AZIMUTH,ELEVATION) spatialises
        %   the signal X using spatialisation data contained in the SOFA
        %   file SOFAFILE at the corresponding AZIMUTH and ELEVATION. The
        %   sampling frequency is that of the SOFA data. X should be a
        %   column vector; if it is not, it will be down-mixed to a column
        %   vector.
        %
        %   IOSR.DSP.AUDIO.SPAT(...,FS) uses the sampling frequency FS for
        %   the spatialisation. If the SOFA data are not at the sampling
        %   rate FS, the data are resampled.

            assert(size(x,1)>1, 'iosr:audio:invalidX', 'X should be a column vector. Matrices will be downmixed to a column vector');

            if size(x,2)>1 % mix to mono
                x = audio.up_down_mix(x,1);
            end
            SOFAobj = SOFAload(SOFAfile); % load SOFA object
            if ~exist('fs','var')==1 % set default FS
                fs = SOFAobj.Data.SamplingRate;
            end

            % convert navigational coordinates to spherical coordinates
            if azimuth<0
                [azimuth,elevation] = nav2sph(azimuth,elevation);
            end

            % find index of corresponding IR is SOFA obj
            dist = (SOFAobj.SourcePosition(:,1)-azimuth).^2 + ...
                (SOFAobj.SourcePosition(:,2)-elevation).^2;
            [~,idx] = min(dist);

            % return IR
            IR = squeeze(SOFAobj.Data.IR(idx,:,:))';

            % resample?
            if SOFAobj.Data.SamplingRate~=fs
                disp(['Resampling audio to ' num2str(fs) ' Hz.'])
                IR = resample(IR,fs,SOFAobj.Data.SamplingRate);
            end

            % convolve
            y = zeros(length(x)+size(IR,1)-1,size(IR,2));
            for c = 1:size(IR,2)
                y(:,c) = iosr.dsp.convFft(x,IR(:,c));
            end

            % crop
            y = y(1:length(x),:);

        end

        function y = up_down_mix(x,N)
        %UP_DOWN_MIX Remix a source to a designated channel count
        %
        %   Y = IOSR.DSP.AUDIO.UP_DOWN_MIX(X,N) mixes the signal X to have
        %   the specified channel count N. Data for each channel in X
        %   should be stored down its columns. For input channel count M,
        %   the up-/down-mix coefficients are calculated as
        %   ONES(M,N)./(M*N).

            M = size(x,2);
            y = x*(ones(M,N)./(M*N));

        end

    end

    methods (Static, Access = protected)

        function ensure_path(filename)
        %ENSURE_PATH Ensure path exists for specified filename

            % get path
            filepath = fileparts(filename);
            if ~isempty(filepath)
                % ensure path exits
                if exist(filepath,'dir')~=7
                    mkdir(filepath)
                end
            end

        end

    end

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
function view = adjustMRImage(view);
% view = adjustMRImage(view);
%
% A button-down function for the MR image
% on a view window. This is set to act when
% the user clicks on the image, and does
% different things depending upon which button
% was pressed. It does the following:
%
% Left click ->
% Middle click -> Adjust brightness, contrast
% Right click ->
% Double click ->
%
% ras 11/04
fig = view.ui.windowHandle;
axs = gcbo;

% figure out which button was pressed
button = get(fig, 'SelectionType');
if strcmp(button,'open')   % double click
   button = 4;
elseif strcmp(button,'normal')  % left
   button = 1;
elseif strcmp(button,'extend')  % middle / shift-left
   button = 2;
elseif strcmp(button,'alt')     % right / ctrl-left
   button = 3;
else
    % weird selection -- just quit
    return
end

switch button
    case 1,
%         % a test
%         buttonisdown = 1;
%         p1 = get(gca,'CurrentPoint');
%         ax = axis;
%         tic
%         while buttonisdown
%             p2 = get(gca,'CurrentPoint');
%             dx = p2(1,1) - p1(1,1);
%             dy = p2(1,2) - p1(1,2);
%             view = adjBrightness(view,dx/ax(4));
%             tmp = get(gcf,'SelectionType');
%             buttonisdown = (toc<1);
%             p1 = p2;
%         end
    case 2,
        % adjust brightness/contrast of image, while button down
        view = adjBrightness(view,0.1);
    case 3,
        % adjust brightness/contrast of image, while button down
        view = adjBrightness(view,-0.1);
    case 4,
        if isequal(view.viewType,'Flat')
            % set params for thresholding curvature

        else
            % a crude manual setting of brightness
            p1 = get(gca,'CurrentPoint');
            ax = axis;
            delta = 2*p1(1,2)/ax(4) - 1;
            view = adjBrightness(view,delta);
        end
end

view = refreshScreen(view);

% reset the call to this functions
cbstr = sprintf('%s = adjustMRImage(%s);',view.name,view.name);
set(gca,'ButtonDownFcn',cbstr);
tmp = get(gca,'Children');
set(tmp,'ButtonDownFcn',cbstr);

return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function view = adjBrightness(view,delta);
% alter the brightness of the image
cmap = view.ui.anatMode.cmap;
cmap = brighten(cmap,delta);
numGrays = view.ui.mapMode.numGrays;
view.ui.anatMode.cmap = cmap;
view.ui.ampMode.cmap(1:numGrays,:) = cmap;
view.ui.phMode.cmap(1:numGrays,:) = cmap;
view.ui.coMode.cmap(1:numGrays,:) = cmap;
view.ui.mapMode.cmap(1:numGrays,:) = cmap;
return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function adjContrast(delta);
% alter the contrast of the image
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function H = topple_contactHorizontal(th,P)

% AUTOMATICALLY GENERATED  --  DO NOT EDIT
% Computes horizontal component of contact force for a stick toppling from rest.

m = P.m;
g = P.g;
L = P.L;
I = P.I;

H = -(L.^2.*g.*m.^2.*(3.*sin(2.*th) - 4.*sin(th)))./(2.*(I + L.^2.*m));

end

author: MatthewPeterKelly

repo: dscTutorials


text:
% OP_UDOTN_VDOTN: assemble the matrix M = [m(i,j)], m(i,j) = (mu u_j * n, v_i * n), with n the exterior normal vector.
%
%   mat = op_udotn_vdotn (spu, spv, msh, coeff);
%   [rows, cols, values] = op_udotn_vdotn (spu, spv, msh, coeff);
%
% INPUT:
%
%  spu:   structure representing the space of trial functions (see sp_vector/sp_eval_boundary_side)
%  spv:   structure representing the space of test functions  (see sp_vector/sp_eval_boundary_side)
%  msh:   structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_eval_boundary_side)
%  coeff: physical parameter
%
% OUTPUT:
%
%   mat:    assembled matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2009, 2010 Carlo de Falco, Rafael Vazquez
% Copyright (C) 2011, 2014, 2017 Rafael Vazquez
%
%    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/>.

function varargout = op_udotn_vdotn (spu, spv, msh, coeff)

  shpu = reshape (spu.shape_functions, spu.ncomp, msh.nqn, spu.nsh_max, msh.nel);
  shpv = reshape (spv.shape_functions, spv.ncomp, msh.nqn, spv.nsh_max, msh.nel);

  rows = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights .* coeff;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:, iel)))
      shpu_iel = reshape (shpu(:, :, :, iel), spu.ncomp, msh.nqn, 1, spu.nsh_max);
      shpv_iel = reshape (shpv(:, :, :, iel), spv.ncomp, msh.nqn, spv.nsh_max, 1);

      shpu_dot_n = sum (bsxfun (@times, shpu_iel, msh.normal(:,:,iel)), 1);
      shpv_dot_n = sum (bsxfun (@times, shpv_iel, msh.normal(:,:,iel)), 1);

      jacdet_iel = reshape (jacdet_weights(:,iel), [1, msh.nqn]);
      shpu_times_jw = bsxfun (@times, jacdet_iel, shpu_dot_n);
      tmp1 = bsxfun (@times, shpu_times_jw, shpv_dot_n);
      elementary_values = reshape (sum (tmp1, 2), spv.nsh_max, spu.nsh_max);

      [rows_loc, cols_loc] = ndgrid (spv.connectivity(:,iel), spu.connectivity(:,iel));
      indices = rows_loc & cols_loc;
      rows(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = rows_loc(indices);
      cols(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = cols_loc(indices);
      values(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = elementary_values(indices);
      ncounter = ncounter + spu.nsh(iel)*spv.nsh(iel);

    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_udotn_vdotn: singular map in element number %d', iel)
    end
  end

  if (nargout == 1 || nargout == 0)
    varargout{1} = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), spv.ndof, spu.ndof);
  elseif (nargout == 3)
    varargout{1} = rows(1:ncounter);
    varargout{2} = cols(1:ncounter);
    varargout{3} = values(1:ncounter);
  else
    error ('op_udotn_vdotn: wrong number of output arguments')
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function pred = nnrpred2(Xtest,X,y,y_exp,K,dist_type,pret_type)

% prediction of new samples with knn regression model
%
% pred = nnrpred(Xtest,X,y,K,dist_type,pret_type)
%
% ------------ INPUT ---------------------------------------------------
% Xtest:        dataset to be predicted [n_test x p] n objects, p variables
% X:            training data matrix (n x p)
% y:            training y vector (n x 1)
% K:            number of neighbors
% dist_type:    'euclidean' Euclidean distance
%               'mahalanobis' Mahalanobis distance
%               'cityblock' City Block metric
%               'minkowski' Minkowski metric
% pret_type:    'cent' cenering
%               'scal' variance scaling
%               'auto' for autoscaling (centering + variance scaling)
%               'rang' range scaling (0-1)
%               'fp'   fingerprints
%
% ------------ OUTPUT --------------------------------------------------
% pred is a structure conyaining
% y_pred            predicted y vector [n_test x 1]
% y_pred_weighted   predicted wighted y vector [n_test x 1]
% neighbors         list of k neighbors for each predicted sample [n_test x k]
%
% version 2.0 - may 2012
% Davide Ballabio
% Milano Chemometrics and QSAR Research Group
% www.disat.unimib.it/chm

[n,p] = size(Xtest);
[X_scal_train,param] = data_pretreatment(X,pret_type);
X_scal = test_pretreatment(Xtest,param);
Xd = [X_scal;X_scal_train];
D = nnr_calc_dist(X_scal_train,X_scal,dist_type,pret_type);
neighbors = zeros(n,K);
yc=zeros(1,n);
yc_weighted=zeros(1,n);
w=zeros(n,K);
y_exp(find(isnan(y_exp)))=y(find(isnan(y_exp)));
for i=1:n
    D_in = D(i,:);
    [d_tmp,n_tmp] = sort(D_in);
    neighbors(i,:) = n_tmp(1:K);
    d_neighbors = d_tmp(1:K);
    if d_neighbors(1)<1e-5 %&& d_neighbors(2)>d_neighbors(1)
        d_neighbors(1)=0;
    end
    if d_neighbors(1)==0 %&& isnan(y_exp(neighbors(i,1)))%&& d_neighbors(2)~=0
        yc(i) = mean(y(neighbors(i,:)));
        yc_weighted(i)=y(neighbors(i,1));
        w(i,1)=1;
    else
        %y(find(~isnan(y_exp)))=y_exp(find(~isnan(y_exp)));
        [yc(i),yc_weighted(i),w(i,:)] = nnrcalcy2(y_exp(neighbors(i,:)),d_neighbors,K);
    end

    dc(i,:)=d_neighbors;
end

pred.neighbors  = neighbors;
pred.y_pred = yc';
pred.y_pred_weighted = yc_weighted';
pred.D=D;
pred.dc=dc;
pred.w=w;

author: kmansouri

repo: OPERA

save_path: github-repos/MATLAB/kmansouri-OPERA