text:
classdef CMatrixTest < matlab.unittest.TestCase
   properties (TestParameter)
      type = {@ddouble}
      lhsSize = {[0, 0], [1, 1], [30, 20], [30, 1], [1 20]}
      rhsSize = {[0, 0], [1, 1], [30, 20], [30, 1], [1 20]}
      lhsMode = struct('dense', 0, 'sparse', 1);
      rhsMode = struct('dense', 0, 'sparse', 1);
   end

   methods (Static)
      function [okay, A1, B1] = generateMatrices(lhsSize, rhsSize, lhsMode, rhsMode)
         Am = lhsSize(1); An = lhsSize(2); Bm = rhsSize(1); Bn = rhsSize(2);
         A1 = randn(Am, An); B1 = randn(Bm, Bn);
         okay = false;

         if lhsMode == 1, A1 = sparse(A1); end
         if rhsMode == 1, B1 = sparse(B1); end

         % check compatibility
         m = Am; n = An;
         if (Bm ~= 1), m = Bm; end
         if (Bn ~= 1), n = Bn; end
         if (Am ~= 1 && Am ~= m), return; end
         if (Bm ~= 1 && Bm ~= m), return; end
         if (An ~= 1 && An ~= n), return; end
         if (Bn ~= 1 && Bn ~= n), return; end

         okay = true;
      end
   end

   methods (Test)
      function comparisons(testCase, type, lhsSize, rhsSize, lhsMode, rhsMode)
         [okay, A1, B1] = CMatrixTest.generateMatrices(lhsSize, rhsSize, lhsMode, rhsMode);
         if (~okay), return; end

         % lt, gt, le, ge
         A2 = type(A1); B2 = type(B1);
         testCase.verifyEqual(A2 < B2, A1 < B1)
         testCase.verifyEqual(A2 > B2, A1 > B1)
         testCase.verifyEqual(A2 <= B2, A1 <= B1)
         testCase.verifyEqual(A2 >= B2, A1 >= B1)

         % ne, eq, and, or, not
         A1 = round(A1); B1 = round(B1);
         A2 = type(A1); B2 = type(B1);
         testCase.verifyEqual(A2 ~= B2, A1 ~= B1)
         testCase.verifyEqual(A2 == B2, A1 == B1)
         testCase.verifyEqual(A2 & B2, A1 & B1)
         testCase.verifyEqual(A2 | B2, A1 | B1)
         testCase.verifyEqual(~A2, ~A1)
      end

      function arithmetic(testCase, type, lhsSize, rhsSize, lhsMode, rhsMode)
         [okay, A1, B1] = CMatrixTest.generateMatrices(lhsSize, rhsSize, lhsMode, rhsMode);
         Ceps = double(eps(type(1))) + eps;
         if (~okay), return; end

         % plus, minus, times, rdivide, ldivide
         A2 = type(A1); B2 = type(B1);
         testCase.verifyEqual(double(A2 + B2), A1 + B1, 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(A2 - B2), A1 - B1, 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(A2 .* B2), A1 .* B1, 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(A2 ./ B2), A1 ./ B1, 'AbsTol', Ceps*1e4, 'RelTol', Ceps*1e4)
         testCase.verifyEqual(double(A2 .\ B2), A1 .\ B1, 'AbsTol', Ceps*1e4, 'RelTol', Ceps*1e4)

         % max, min
         testCase.verifyEqual(double(max(A2,B2)), max(A1,B1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(min(A2,B2)), min(A1,B1), 'AbsTol', Ceps*1e4)

         % horzcat, vertcat
         if size(A2,1) == size(B2,1)
            testCase.verifyEqual(double([A2 B2]), [A1 B1], 'AbsTol', Ceps*1e4)
            testCase.verifyEqual(double([A2 B2 A2 B2]), [A1 B1 A1 B1], 'AbsTol', Ceps*1e4)
         end

         if size(A2,2) == size(B2,2)
            testCase.verifyEqual(double([A2;B2]), [A1;B1], 'AbsTol', Ceps*1e4)
            testCase.verifyEqual(double([A2;B2;A2;B2]), [A1;B1;A1;B1], 'AbsTol', Ceps*1e4)
         end

         % times, rdivide, ldivide
         if size(A1,2) == size(B2',1)
            testCase.verifyEqual(double(A2 * B2'), A1 * B1', 'AbsTol', Ceps*1e6)
         end

         if size(A1,1) == size(B2,1)
            testCase.verifyEqual(double(A2 \ B2), A1 \ B1, 'AbsTol', Ceps*1e6)
         end

         if size(A1,2) == size(B2,2)
            testCase.verifyEqual(double(A2 / B2), A1 / B1, 'AbsTol', Ceps*1e6)
         end
      end

      function unaryOperations(testCase, type, lhsSize, lhsMode)
         A1 = randn(lhsSize);
         if lhsMode == 1, A1 = sparse(A1); end
         Ceps = double(eps(type(1))) + eps;

         % size, length, numel
         A2 = type(A1);
         testCase.verifyEqual(size(A2), size(A1))
         testCase.verifyEqual(size(A2, 1), size(A1, 1))
         testCase.verifyEqual(size(A2, 2), size(A1, 2))
         [a1, b1] = size(A1);
         [a2, b2] = size(A2);
         testCase.verifyEqual(a2, a1)
         testCase.verifyEqual(b2, b1)
         testCase.verifyEqual(length(A2), length(A1))
         testCase.verifyEqual(numel(A2), numel(A1))

         % isscalar, isvector, ismatrix, isempty, isrow, iscolumn, issymmetric
         testCase.verifyEqual(isscalar(A2), isscalar(A1))
         testCase.verifyEqual(isvector(A2), isvector(A1))
         testCase.verifyEqual(ismatrix(A2), ismatrix(A1))
         testCase.verifyEqual(isempty(A2), isempty(A1))
         testCase.verifyEqual(isrow(A2), isrow(A1))
         testCase.verifyEqual(iscolumn(A2), iscolumn(A1))
         testCase.verifyEqual(issymmetric(A2), issymmetric(A1))

         if (size(A1,1) == size(A1,2))
            A1 = A1 + A1'; A2 = type(A1);
            testCase.verifyEqual(issymmetric(A2), issymmetric(A1))
         end

         % uminus, uplus, abs, transpose, ctranspose, sqrt, all, any, nnz, diag
         testCase.verifyEqual(double(-A2), -A1, 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(+A2), +A1, 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(abs(A2)), abs(A1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(transpose(A2)), transpose(A1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(ctranspose(A2)), ctranspose(A1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(sqrt(abs(A2)+1)), sqrt(abs(A1)+1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(all(A2), all(A1))
         testCase.verifyEqual(all(A2,1), all(A1,1))
         testCase.verifyEqual(all(A2,2), all(A1,2))
         testCase.verifyEqual(any(A2), any(A1))
         testCase.verifyEqual(any(A2,1), any(A1,1))
         testCase.verifyEqual(any(A2,2), any(A1,2))
         testCase.verifyEqual(nnz(A2), nnz(A1))
         testCase.verifyEqual(double(diag(A2)), diag(A1), 'AbsTol', Ceps*1e4)

         % sum, prod, max, min, norm
         testCase.verifyEqual(double(sum(A2)), sum(A1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(sum(A2,1)), sum(A1,1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(sum(A2,2)), sum(A1,2), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(sum(A2,'all')), sum(A1,'all'), 'AbsTol', Ceps*1e4)

         testCase.verifyEqual(double(prod(A2)), prod(A1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(prod(A2,1)), prod(A1,1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(prod(A2,2)), prod(A1,2), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(prod(A2,'all')), prod(A1,'all'), 'AbsTol', Ceps*1e4)

         testCase.verifyEqual(double(max(A2)), max(A1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(max(A2,[],1)), max(A1,[],1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(max(A2,[],2)), max(A1,[],2), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(max(A2,[],'all')), max(A1,[],'all'), 'AbsTol', Ceps*1e4)

         testCase.verifyEqual(double(min(A2)), min(A1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(min(A2,[],1)), min(A1,[],1), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(min(A2,[],2)), min(A1,[],2), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(min(A2,[],'all')), min(A1,[],'all'), 'AbsTol', Ceps*1e4)

         if isvector(A1)
            testCase.verifyEqual(double(norm(A2)), norm(A1), 'AbsTol', Ceps*1e4)
         end

         % find
         testCase.verifyEqual(find(A2), find(A1))
         [i2,j2] = find(A2); [i1,j1] = find(A1);
         testCase.verifyEqual(i2, i1)
         testCase.verifyEqual(j2, j1)
         [i2,j2,v2] = find(A2); [i1,j1,v1] = find(A1);
         testCase.verifyEqual(i2, i1)
         testCase.verifyEqual(j2, j1)
         testCase.verifyEqual(double(v2), v1, 'AbsTol', Ceps*1e4)

         % chol
         H1 = A1 * A1' + speye(size(A1,1));
         testCase.verifyEqual(double(chol(type(H1))), chol(H1), 'AbsTol', Ceps*1e6)
      end

      function externalFunc(testCase, type, lhsSize)
         Am = lhsSize(1); An = lhsSize(2);
         Ceps = double(eps(type(1))) + eps;
         typename = class(type(1.0));

         % ones, zeros, eye, rand, randn, randi, sparse
         A1 = ones(Am, An, typename); A2 = ones(Am, An);
         testCase.verifyTrue(all(A1 == A2, 'all'))
         A1 = zeros(Am, An, typename); A2 = zeros(Am, An);
         testCase.verifyTrue(all(A1 == A2, 'all'))
         A1 = eye(Am, An, typename); A2 = eye(Am, An);
         testCase.verifyTrue(all(A1 == A2, 'all'))
         rng(1); A1 = rand(Am, An, typename); rng(1); A2 = rand(Am, An);
         testCase.verifyTrue(all(A1 == A2, 'all'))
         rng(1); A1 = randn(Am, An, typename); rng(1); A2 = randn(Am, An);
         testCase.verifyTrue(all(A1 == A2, 'all'))
         rng(1); A1 = randi(10, Am, An, typename); rng(1); A2 = randi(10, Am, An);
         testCase.verifyTrue(all(A1 == A2, 'all'))
         A1 = full(sprand(Am, An, 0.5)); A2 = type(A1);
         testCase.verifyEqual(double(sparse(A2)), sparse(A1), 'AbsTol', 1e4*Ceps);
         testCase.verifyEqual(double(full(A2)), full(A1), 'AbsTol', 1e4*Ceps);
         A1 = (sprand(Am, An, 0.5)); A2 = type(A1);
         testCase.verifyEqual(double(sparse(A2)), sparse(A1), 'AbsTol', 1e4*Ceps);
         testCase.verifyEqual(double(full(A2)), full(A1), 'AbsTol', 1e4*Ceps);
         Carr = type(1:Am);
         testCase.verifyEqual(double(sparse(1:Am, 1:Am, Carr)), sparse(1:Am, 1:Am, 1:Am));
         testCase.verifyEqual(double(sparse(1:Am, 1:Am, Carr, Am, Am)), sparse(1:Am, 1:Am, 1:Am, Am, Am));
         testCase.verifyEqual(double(sparse(1:Am, 1:Am, Carr, Am, Am, Am)), sparse(1:Am, 1:Am, 1:Am, Am, Am, Am));

         % .x
         b = A2.x;
         A2.x = b;
      end

      function otherTests(testCase, type, lhsMode)
         A1 = sprandn(5, 6, 0.6);
         if lhsMode == 0, A1 = full(A1); end
         A2 = type(A1);
         Ceps = double(eps(type(1))) + eps;

         % subsref
         testCase.verifyEqual(double(A2(1:3,:)), A1(1:3,:), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(A2(:,2:4)), A1(:,2:4), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(A2(1:3,2:4)), A1(1:3,2:4), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(A2(1:3,2:end)), A1(1:3,2:end), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(A2(5:end)), A1(5:end), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(A2(:)), A1(:), 'AbsTol', Ceps*1e4)
         testCase.verifyEqual(double(A2(2:2:8)), A1(2:2:8), 'AbsTol', Ceps*1e4)

         % subsasgn
         A1(1:3,:) = 5; A2(1:3,:) = 5;
         testCase.verifyEqual(double(A2),A1, 'AbsTol', Ceps*1e4)
         A1(5:end) = 5; A2(5:end) = 5;
         testCase.verifyEqual(double(A2),A1, 'AbsTol', Ceps*1e4)
         A1(:,2:4) = []; A2(:,2:4) = [];
         testCase.verifyEqual(double(A2),A1, 'AbsTol', Ceps*1e4)

         disp(A2);

         % dissect, amd, symbfact
         A1 = sprandn(6, 6, 0.6);
         A1 = A1 * A1';
         A2 = type(A1);
         testCase.verifyEqual(dissect(A2),dissect(A1))
         testCase.verifyEqual(amd(A2),amd(A1))
         testCase.verifyEqual(symbfact(A2),symbfact(A1))
         [r21, r22] = etree(A2);
         [r11, r12] = etree(A1);
         testCase.verifyEqual(r21,r11)
         testCase.verifyEqual(r22,r12)

         % chol
         try
            A = randn(5,5);
            A = A * A';
            A = sparse(A) * NaN;
            Z = chol(type(A));
            Z'\randn(5,1);
         end

         A = sprandn(30,30,0.3);
         H = type(A * A');
         R = chol(H);
         D = H - R' * R;
         testCase.verifyLessThan(double(norm(D(:))), Ceps*1e4)

         % solves for triangular matrix
         A = type(tril(sprandn(30,30,0.3) + speye(30)));
         b = randn(30,1);
         x = A\b;
         testCase.verifyLessThan(double(norm(A*x-b)), Ceps*1e4)

         A = type(triu(sprandn(30,30,0.3) + speye(30)));
         b = randn(30,1);
         x = A\b;
         testCase.verifyLessThan(double(norm(A*x-b)), Ceps*1e4)

         A = type(tril(randn(30,30) + eye(30)));
         b = randn(30,1);
         x = A\b;
         testCase.verifyLessThan(double(norm(A*x-b)), Ceps*1e4)

         A = type(triu(randn(30,30) + eye(30)));
         b = randn(30,1);
         x = A\b;
         testCase.verifyLessThan(double(norm(A*x-b)), Ceps*1e4)
      end

      function cholTests(testCase)
         load('..\..\Problem\LPnetlib\lp_80bau3b.mat')

         A = Problem.A;
         p = colamd(A');
         A = A(p,:);
         w = rand(size(A,2),1);
         R = chol((A*(w.*A')));
         x = rand(size(A,1),1);
         z = R\(R'\x);
         d = full(diag(R));

         o = AdaptiveChol(A);
         o.factorize(diag(sparse(w)));
         z2 = o.solve(x);
         testCase.verifyEqual(double(z), z2, 'AbsTol', eps*1e4)

         d2 = o.diagonal();
         testCase.verifyEqual(double(d), d2, 'AbsTol', eps*1e4)

         ls = o.leverageScore(100);
         testCase.verifyTrue(all(ls<1.5));

         try
            o.factorize(diag(sparse(rand(12,1))));
         end
      end
   end
end

author: opencobra

repo: cobratoolbox


text:
%LAMBDA2RGB RGB chromaticity coordinates
%
% RGB = LAMBDA2RG(LAMBDA) is the rg-chromaticity coordinate (1x2) for
% illumination at the specific wavelength LAMBDA [m]. If LAMBDA is a
% vector (Nx1), then P (Nx2) is a vector whose elements are the chromaticity
% coordinates at the corresponding elements of LAMBDA.
%
% RGB = LAMBDA2RG(LAMBDA, E) is the rg-chromaticity coordinate (1x2) for an
% illumination spectrum E (Nx1) defined at corresponding wavelengths
% LAMBDA (Nx1).
%
% References::
%  - Robotics, Vision & Control, Section 10.2,
%    P. Corke, Springer 2011.
%
% See also CMFRGB, LAMBDA2XY.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
function [r,g] = lambda2rg(lambda, e)
    if nargin == 1,
        RGB = cmfrgb(lambda);
    elseif nargin == 2,
        RGB = cmfrgb(lambda, e);
    end
    cc = tristim2cc(RGB);

    if nargout == 1
        r = cc;
    elseif nargout == 2
        r = cc(:,1);
        g = cc(:,2);
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
function x = cg_ge ( n, a, b, x )

%*****************************************************************************80
%
%% CG_GE uses the conjugate gradient method for a general storage matrix.
%
%  Discussion:
%
%    The linear system has the form A*x=b, where A is a positive-definite
%    symmetric matrix, stored as a full storage matrix.
%
%    The method is designed to reach the solution to the linear system
%      A * x = b
%    after N computational steps.  However, roundoff may introduce
%    unacceptably large errors for some problems.  In such a case,
%    calling the routine a second time, using the current solution estimate
%    as the new starting guess, should result in improved results.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Frank Beckman,
%    The Solution of Linear Equations by the Conjugate Gradient Method,
%    in Mathematical Methods for Digital Computers,
%    edited by John Ralston, Herbert Wilf,
%    Wiley, 1967,
%    ISBN: 0471706892,
%    LC: QA76.5.R3.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N,N), the matrix.
%
%    Input, real B(N), the right hand side vector.
%
%    Input/output, real X(N).
%    On input, an estimate for the solution, which may be 0.
%    On output, the approximate solution vector.
%

%
%  Initialize
%    AP = A * x,
%    R  = b - A * x,
%    P  = b - A * x.
%
  ap(1:n,1) = a * x;

  r(1:n,1) = b(1:n,1) - ap(1:n,1);
  p(1:n,1) = b(1:n,1) - ap(1:n,1);
%
%  Do the N steps of the conjugate gradient method.
%
  for it = 1 : n
%
%  Compute the matrix*vector product AP = A*P.
%
    ap(1:n,1) = a * p;
%
%  Compute the dot products
%    PAP = P*AP,
%    PR  = P*R
%  Set
%    ALPHA = PR / PAP.
%
    pap = p(1:n,1)' * ap(1:n,1);
    pr =  p(1:n,1)' * r(1:n,1);

    if ( pap == 0.0 )
      return;
    end

    alpha = pr / pap;
%
%  Set
%    X = X + ALPHA * P
%    R = R - ALPHA * AP.
%
    x(1:n,1) = x(1:n,1) + alpha * p(1:n,1);
    r(1:n,1) = r(1:n,1) - alpha * ap(1:n,1);
%
%  Compute the vector dot product
%    RAP = R*AP
%  Set
%    BETA = - RAP / PAP.
%
    rap = r(1:n,1)' * ap(1:n,1);

    beta = - rap / pap;
%
%  Update the perturbation vector
%    P = R + BETA * P.
%
    p(1:n,1) = r(1:n,1) + beta * p(1:n,1);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function tifd=friedman(tfr,hat,t,method,trace);
%FRIEDMAN Instantaneous frequency density.
%	TIFD = FRIEDMAN(TFR,HAT,T,METHOD,TRACE) computes the
%	time-instantaneous frequency density (defined by Friedman [1])
%	of a reassigned time-frequency representation.
%
%	TFR   : time-frequency representation, (N,M) matrix.
%	HAT   : complex matrix of the reassignment vectors.
%	T     : the time instant(s)	(default : (1:M)).
%	METHOD: chosen representation	(default : 'tfrrsp').
%	TRACE : if nonzero, the progression of the algorithm is shown
%					(default : 0).
%	TIFD  : time instantaneous-frequency density. When called without
%	        output arguments, FRIEDMAN runs TFRQVIEW.
%
%	WARNING : TIFD is not an energy distribution, but an estimated
%	-------        probability distribution !
%
%	Example :
%	 sig=fmlin(128,0.1,0.4); h=tftb_window(47,'Kaiser');
%	 t=1:2:127; [tfr,rtfr,hat]=tfrrpwv(sig,t,128,h);
%	 friedman(tfr,hat,t,'tfrrpwv',1);
%
%	See also : RIDGES.

%	F. Auger, August 1994, Decembre 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%
%	[1] : D. H. Friedman, "Instantaneous Frequency vs Time : An
%	      Interpretation of the Phase Structure of Speech", Proc. IEEE
%	      ICASSP, pp. 29.10.1-4, Tampa, 1985.

if (nargin<3),
 error('At least 2 parameters required');
end;

[tfrrow,tfrcol]=size(tfr);
[hatrow,hatcol]=size(hat);

if (nargin==2),
 t=1:tfrcol; method='tfrrsp'; trace=0;
elseif (nargin==3),
 method='tfrrsp'; trace=0;
elseif (nargin==4),
 trace=0;
end;

[trow,tcol] = size(t);
if (trow~=1),
 error('T must only have one row');
elseif (tfrrow~=hatrow)|(tfrcol~=hatcol),
 error('tfr and hat must have the same size');
end;

tifd=zeros(tfrrow,tfrcol);
bins=0.5+(0:tfrrow-1);
threshold=sum(sum(tfr))*0.5/(tfrrow*tfrcol);

if trace, fprintf ('\nFriedman distribution: \n'); end;

for j=1:tfrcol,
 if trace, disprog(j,tfrcol,10); end;
 indices=find(tfr(:,j)>threshold);
 if (length(indices)>=1),
  [occurences,trash]=hist(real(hat(indices,j)),bins);
  tifd(:,j)=occurences';
 end;
end;
tifd=tifd/sum(sum(tifd));

method=upper(method);
if nargout==0,
 tfrqview(tifd,[],t,method);
end

author: HeLiangHIT


text:
function net = GNGUpdate(oSignals,net)
% Update GNG

%--------------------------------------------------------------------------
% Copyright Yiping Liu
% Please contact {yiping0liu@gmail.com} if you have any problem.
%--------------------------------------------------------------------------

    %% Parameters
    [nSig,D]=size(oSignals);
    maxAge = net.maxAge;
    lambda = net.lambda;
    epsilon_a = net.epsilon_a;
    epsilon_nb = net.epsilon_nb;
    alpha = net.alpha;
    delta = net.delta;
    maxHP = net.maxHP;
    maxIter = net.maxIter;
    Node = net.Node;
    Err  = net.Err;
    edge = net.edge;
    age  = net.age;
    hp  = net.hp;

    %% Randamize Input Signal
    ran = randperm(nSig);
    oSignals = oSignals(ran,:);

    %% Input Signal Normalization [0-1]
    oMax = max(oSignals,[],1);
    oMin = min(oSignals,[],1);
    oRange = oMax-oMin;
    Signals = (oSignals - repmat(oMin,nSig,1))./repmat(oRange,nSig,1);

    %% Update GNG
    for nitr = 1:maxIter
        % Step 0: Initialization. Start with two neural units (nodes) selected from input data
        if size(Node,1) <= 2
            Ni = 2;
            Xmin = min(Signals,[],1);
            Xmax = max(Signals,[],1);
            for i = 1:Ni
                Node(i,:) = unifrnd(Xmin, Xmax);
            end
            Err = [0; 0];
            edge = zeros(2,2);
            age  = zeros(2,2);
            hp = ones(1,2).*maxHP;
        end

        for numSig = 1:nSig

            % Step 1: Input one signal
            pattern = Signals(numSig,:);

            % Step 2: Find the two nearest nodes ra and rb to new signal
            d = pdist2(pattern, Node);
            [~, SortOrder] = sort(d);
            ra = SortOrder(1);
            rb = SortOrder(2);

            % Step 2.5: change HP
            hp = hp-1;
            hp(ra) = maxHP;
            hp(rb) = hp(rb)+1;

            % Steps 3: Increment the age of all edges emanating from ra
            age(ra, :) = age(ra, :) + 1;
            age(:, ra) = age(:, ra) + 1;

            % Step 4: Add the squared distance to a local error counter variable
            Err(ra) = Err(ra) + d(ra)^2;

            % Step 5: Move ra and its topological neighbors towards singal
            Node(ra,:) = Node(ra,:) + epsilon_a * (pattern - Node(ra,:));
            for j = find(edge(ra,:)==1)
                Node(j,:) = Node(j,:) + epsilon_nb * (pattern - Node(j,:));% for Neighbor nodes which are connecting to ra.
            end

            % Step 6:
            % If ra and rb are connected by an edge, set the age of this edge to zero.
            % If such an edge does not exist, create it.
            edge(ra,rb) = 1;
            edge(rb,ra) = 1;
            age(ra,rb) = 0;
            age(rb,ra) = 0;

            % Step 7(1):
            % Remove edges from node if age>maxAge.
            edge(age>maxAge) = 0;

            % Step 7(2):
            % Remove dead node and their edges.
            DeadNodes = (hp<=0);
            edge(DeadNodes, :) = [];
            edge(:, DeadNodes) = [];
            age(DeadNodes, :) = [];
            age(:, DeadNodes) = [];
            Node(DeadNodes, :) = [];
            Err(DeadNodes) = [];
            hp(DeadNodes) = [];

            % Step 8: Node Insertion Procedure.
            % rnew: new node
            % r1max: node which has maximum accumulated error
            % r2max: neighbor node of r1max
            if mod(numSig, lambda) == 0 && net.maxNode > size(Node,1)
                [~, r1max] = max(Err);
                [~, r2max] = max(edge(:,r1max).*Err);
                rnew = size(Node,1) + 1;
                Node(rnew,:) = (Node(r1max,:) + Node(r2max,:))/2;
                edge(r1max,r2max) = 0;
                edge(r2max,r1max) = 0;
                edge(r1max,rnew) = 1;
                edge(rnew,r1max) = 1;
                edge(rnew,r2max) = 1;
                edge(r2max,rnew) = 1;
                age(rnew,:) = 0;
                age(:,rnew) = 0;
                Err(r1max) = alpha * Err(r1max);
                Err(r2max) = alpha * Err(r2max);
                Err(rnew) = Err(r1max);
                hp(rnew) = maxHP;
            end

            % Step 9: Decrease the error of all units.
            Err = delta*Err;

        end
    end

    %% Output net
    net.Node = Node;
    net.Err = Err;
    net.edge = edge;
    net.age = age;
    net.hp = hp;

    %% Expansion (Algorithm 3)
    % Accociate signal to its closest node
    Distance = pdist2(Node,Signals);
    [~,pi] = min(Distance,[],1);
    % Node Label based on edge (detect sub-network)
    connection = graph(edge ~= 0);
    NetLabel = conncomp(connection);
    % Data Label based on Node
    DataLable = NetLabel(pi);
    % Expand each sub-network
    for i=1:max(NetLabel)
        subData = find(DataLable == i);
        subNet = find(NetLabel == i);
        if length(subData)<=1 || length(subNet)<=1
            continue;
        end
        DataMax = max(Signals(subData,:),[],1);
        DataMin = min(Signals(subData,:),[],1);
        NetMax = max(Node(subNet,:),[],1);
        NetMin = min(Node(subNet,:),[],1);
        DataRange = DataMax-DataMin;
        NetRange = NetMax-NetMin;
        Ratio = DataRange./NetRange;
        for k = 1:D
            for j = subNet
                Node(j,k)= (Node(j,k)- NetMin(k)).*Ratio(k)+DataMin(k);
            end
        end
    end
    net.NodeS = Node;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef nme_gen_acc_nln < mp.nme_gen_acc & mp.nme_wrapper_ac_nln

%   MATPOWER
%   Copyright (c) 2019-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

    properties
        nme_class = @mp.nme_gen_acc;
    end

    methods
        function obj = nme_gen_acc_nln()
            obj@mp.nme_gen_acc();
            obj.nme_wrapper_ac_nln_init();
        end

        function build_params(obj, nm, dm)
            build_params@mp.nme_gen_acc(obj, nm, dm);
            obj.build_nln_params(nm, dm);
        end

        function nk = count(obj, dm)
            obj.nme.count(dm);
            nk = count@mp.nme_gen_acc(obj, dm);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
clear all; close all;
I=imread('coins.png');
I=imnoise(I, 'gaussian', 0, 0.01);
I=im2double(I);
M=2*size(I,1);
N=2*size(I,2);
u=-M/2:(M/2-1);
v=-N/2:(N/2-1);
[U,V]=meshgrid(u, v);
D=sqrt(U.^2+V.^2);
D0=50;
W=30;
H=double(or(D<(D0-W/2), D>D0+W/2));
J=fftshift(fft2(I, size(H, 1), size(H, 2)));
K=J.*H;
L=ifft2(ifftshift(K));
L=L(1:size(I,1), 1:size(I, 2));
figure;
subplot(121);
imshow(I);
subplot(122);
imshow(L);

author: HuangCongQing


text:
function [nrm] = normals(pnt, tri, opt)

% NORMALS compute the surface normals of a triangular mesh
% for each triangle or for each vertex
%
% Use as
%   [nrm] = normals(pnt, tri, opt)
% where opt is either 'vertex' (default) or 'triangle'.
%
% See also PCNORMALS, PROJECTTRI

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

if nargin<3
  opt='vertex';
elseif (opt(1)=='v' || opt(1)=='V')
  opt='vertex';
elseif (opt(1)=='t' || opt(1)=='T')
  opt='triangle';
else
  ft_error('invalid optional argument');
end

npnt = size(pnt,1);
ntri = size(tri,1);

% shift to center
pnt(:,1) = pnt(:,1)-mean(pnt(:,1),1);
pnt(:,2) = pnt(:,2)-mean(pnt(:,2),1);
pnt(:,3) = pnt(:,3)-mean(pnt(:,3),1);

% compute triangle normals
% nrm_tri = zeros(ntri, 3);
% for i=1:ntri
%   v2 = pnt(tri(i,2),:) - pnt(tri(i,1),:);
%   v3 = pnt(tri(i,3),:) - pnt(tri(i,1),:);
%   nrm_tri(i,:) = cross(v2, v3);
% end

% vectorized version of the previous part
v2 = pnt(tri(:,2),:) - pnt(tri(:,1),:);
v3 = pnt(tri(:,3),:) - pnt(tri(:,1),:);
nrm_tri = cross(v2, v3);

if strcmp(opt, 'vertex')
  % compute vertex normals
  nrm_pnt = zeros(npnt, 3);
  for i=1:ntri
    nrm_pnt(tri(i,1),:) = nrm_pnt(tri(i,1),:) + nrm_tri(i,:);
    nrm_pnt(tri(i,2),:) = nrm_pnt(tri(i,2),:) + nrm_tri(i,:);
    nrm_pnt(tri(i,3),:) = nrm_pnt(tri(i,3),:) + nrm_tri(i,:);
  end
  % normalise the direction vectors to have length one
  nrm = nrm_pnt ./ (sqrt(sum(nrm_pnt.^2, 2)) * ones(1,3));
else
  % normalise the direction vectors to have length one
  nrm = nrm_tri ./ (sqrt(sum(nrm_tri.^2, 2)) * ones(1,3));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fast cross product to replace the MATLAB standard version
function [c] = cross(a,b)
c = [a(:,2).*b(:,3)-a(:,3).*b(:,2) a(:,3).*b(:,1)-a(:,1).*b(:,3) a(:,1).*b(:,2)-a(:,2).*b(:,1)];

author: spm


text:
function RegMeanSquareAffine(handles, metric)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC stateS;
indexS = planC{end};

    [originF, spacingF] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationBaseDataset));

    [originM, spacingM] = getScanOriginSpacing(planC{indexS.scan}(stateS.imageRegistrationMovDataset));

    minstep = str2double(get(handles.affine_minstep, 'string'));
    maxstep = str2double(get(handles.affine_maxstep, 'string'));
    iternum = str2double(get(handles.affine_iternum, 'string'));
    tscale = str2double(get(handles.affine_tscale, 'string'));

    output = cell(1, 18);

    FImg = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).scanArray;
    MImg = planC{indexS.scan}(stateS.imageRegistrationMovDataset).scanArray;

    dimF = size(FImg);
    dimM = size(MImg);
    %generate cliped base dataset
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_baseTrans'));
    if ~isempty(clipBox)
        FImg = FImg(clipBox(2):clipBox(4), clipBox(1):clipBox(3), :);
        originF(1) = originF(1) + spacingF(1)*double(clipBox(1)-1);
        originF(2) = originF(2) + spacingF(2)*double(uint16(dimF(1)-clipBox(4)));
    end
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_baseSag'));
    if ~isempty(clipBox)
        FImg = FImg(:, :, clipBox(2):clipBox(4));
        originF(3) = originF(3) + spacingF(3)*double(clipBox(2)-1);
    end

    %generate cliped move dataset
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_movTrans'));
    if ~isempty(clipBox)
        MImg = MImg(clipBox(2):clipBox(4), clipBox(1):clipBox(3), :);
        originM(1) = originM(1) + spacingM(1)*double(clipBox(1)-1);
        originM(2) = originM(2) + spacingM(2)*double(uint16(dimM(1)-clipBox(4)));
    end
    clipBox = uint16(getappdata(stateS.handle.CERRSliceViewer, 'clipBox_movSag'));
    if ~isempty(clipBox)
        MImg = MImg(:, :, clipBox(2):clipBox(4));
        originM(3) = originM(3) + spacingM(3)*double(clipBox(2)-1);
    end

%downsample the input datasets
    dsampled = 0;
    if (get(handles.dsampleCheck,'Value') == get(handles.dsampleCheck,'Max'))

        tic;
        xyScale =  2; zScale = 2;
        spacingF = [spacingF(1)*xyScale spacingF(2)*xyScale spacingF(3)*zScale];
        spacingM = [spacingM(1)*xyScale spacingM(2)*xyScale spacingM(3)*zScale];
        disp('downSampling ...');
%         FImg = imdesample3d(FImg,xyScale,zScale);
%         MImg = imdesample3d(MImg,xyScale,zScale);
%         FImg=GPReduce2D(FImg,1);
%         MImg=GPReduce2(MImg,1);
        FImg = FImg(1:xyScale:end, 1:xyScale:end, 1:zScale:end);
        MImg = MImg(1:xyScale:end, 1:xyScale:end, 1:zScale:end);

        toc;
        dsampled = 1;
    end

    switch metric
        case 'mean squares'
            metricIndex = 0;
        case 'normalized correlation'
            metricIndex = 1;
    end

    tic;

%flip the source datasets on X for itk coordinate system
    fdim = 1;
    FImg = flipdim(FImg, fdim);
    MImg = flipdim(MImg, fdim);

% transform initializing modes 0:MomentON 1:GeometryOn 2:initail transform On
    initMode = get(handles.InitTrans, 'value');

%prepare the initial transform matrix
    transMB = planC{indexS.scan}(stateS.imageRegistrationBaseDataset).transM;
    transMM = planC{indexS.scan}(stateS.imageRegistrationMovDataset).transM;
    if isempty(transMM), transMM = eye(4); end;
    if isempty(transMB), transMB = eye(4); end;
    transM = inv(transMB) * transMM;

%flip the moving dataset to match the initial transM
    Tf = eye(4);
    flipM = 0;
    if (get(handles.flipX, 'value'))
        MImg = flipdim(MImg, 2);
        Tf = [-1 0 0 2*centerM(1); 0 1 0 0; 0 0 1 0; 0 0 0 1];
        flipM = 2;
    end
    if (get(handles.flipY, 'value'))
        MImg = flipdim(MImg, 1);
        Tf = [1 0 0 0; 0 -1 0 2*centerM(2); 0 0 1 0; 0 0 0 1];
        flipM = 1;
    end
    if (get(handles.flipZ, 'value'))
        MImg = flipdim(MImg, 3);
        Tf = [1 0 0 0; 0 1 0 0; 0 0 -1 2*centerM(3); 0 0 0 1];
        flipM = 3;
    end
    transM = transM * inv(Tf);
    rotM = transM(1:3, 1:3); transV = transM(1:3, 4);
    rotM(1,2) = -rotM(1,2); rotM(2,1) = -rotM(2,1); %rotM = rotM';
    transV(3) = -transV(3);%transV = -transV;

%call registration method
    try
        [im, Rotation, Offset] = Affine3D(int16(FImg), originF, spacingF, ...
                                      int16(MImg), originM, spacingM, ...
                                      minstep, maxstep, iternum, tscale, metricIndex, rotM, transV,initMode);
    catch
        [im, Rotation, Offset] = Affine3D_64(int16(FImg), originF, spacingF, ...
                                      int16(MImg), originM, spacingM, ...
                                      minstep, maxstep, iternum, tscale, metricIndex, rotM, transV,initMode);
    end

    toc;

    im = flipdim(im, fdim);
    %if (flipM), im = flipdim(im, flipM); end;

    output{1} = ['Angle (radians)   = ' num2str(Offset(4))];
    output{2} = ['Angle (degrees)   = ' num2str(Offset(5))];
    output{3} = ['Translation X = ' num2str(Offset(1))];
    output{4} = ['Translation Y = ' num2str(Offset(2))];
    output{5} = ['Translation Z = ' num2str(Offset(3))];
    output{6} = ['Iterations = ' num2str(Offset(6))];
    output{7} = ['Metric Value (mean square) = ' num2str(Offset(7))];
    output{8} = ['Rotation Center X = ' num2str(Offset(9))];
    output{9} = ['Rotation Center Y = ' num2str(Offset(10))];
    output{10} = ['Rotation Center Z = ' num2str(Offset(11))];
    output{11} = ['Offset X = ' num2str(Offset(12))];
    output{12} = ['Offset Y = ' num2str(Offset(13))];
    output{13} = ['Offset Z = ' num2str(Offset(14))];
    output{14} = ['Rotation Center X = ' num2str(Offset(15))];
    output{15} = ['Rotation Center Y = ' num2str(Offset(16))];
    output{16} = ['Rotation Center Z = ' num2str(Offset(17))];

    set(handles.OutputList, 'string', output);

%update the transM;
    Tv = [Offset(1) Offset(2) Offset(3)];
    Cv = [Offset(15) Offset(16) Offset(17)];
    Cv(3) = -Cv(3);

    rot = reshape(Rotation, 3,3);
    offset = [Offset(12) Offset(13) Offset(14)]';

    TM = eye(4);
    TM(:,4) = [offset; 1];

    RM = eye(4);
    RM(1:3, 1:3) = rot;

    newTransform = inv(TM*RM);
    newTransform =  transMB * newTransform * Tf;

    scanSetM = stateS.imageRegistrationMovDataset;
    scanSetF = stateS.imageRegistrationBaseDataset;
    planC{indexS.(stateS.imageRegistrationMovDatasetType)}(scanSetM).transM = newTransform;

% save the resampled dataset
    if (get(handles.saveCheck, 'value'))
        if (~dsampled)
            planC{indexS.scan}(end+1) = planC{indexS.scan}(scanSetF);
            planC{indexS.scan}(end).scanArray = im;
        else
            %imReg = resample(scanSetF, scanSetM); % need to up-sample im ????
            planC{indexS.scan}(end+1) = planC{indexS.scan}(scanSetF);
            planC{indexS.scan}(end).scanArray = imReg;
        end
        controlFrame('refresh');
    end

    sliceCallBack('refresh');

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function alpha = p10_alpha ( )

%*****************************************************************************80
%
%% P10_ALPHA returns ALPHA for problem 10.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real ALPHA, the value of ALPHA.
%
  alpha = 0.0;

  return
end


text:
function add_nuisance_to_SPMcfg(Xn)
% Adds a matrix Xn to the end of covariates of interest in
% xX structure in SPMcfg.mat
%
% :Usage:
% ::
%
%     function add_nuisance_to_SPMcfg(Xn)
%
% :Inputs:
%
%   **oXn:**
%        should contain ALL nuisance covariates and intercepts
%        as in output of tor_get_physio.m
%
% This function is automatically run by tor_get_physio.m
%
% ..
%    Tor Wager
%..

    if exist([pwd filesep 'SPMcfg.mat']) == 2
        load SPMcfg
    else
        warning(['No SPMcfg.mat file found in current directory:' pwd])
    end

    xX.X = [xX.X(:,xX.iC) Xn];

    XNamesB = xX.Xnames(xX.iB);

    for i = max(xX.iC)+1:size(xX.X,2)-length(XNamesB), xX.Xnames{i} = ['Nuisance ' num2str(i)];, end

    xX.Xnames(size(xX.X,2)-length(XNamesB)+1:size(xX.X,2)) = XNamesB;
    xX.iB = max(xX.iC)+1:size(xX.X,2);
    F_iX0.iX0 = xX.iB;
    save SPMcfg F_iX0 SPMid Sess VY xGX xM xX xsDes

    fprintf(1,'\n SPMcfg xX modified and saved.\n')

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function determ = hankel_n_determinant ( n )

%*****************************************************************************80
%
%% HANKEL_N_DETERMINANT returns the determinant of the HANKEL_N matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real DETERM, the determinant.
%
  determ = r8_mop ( floor ( n / 2 ) ) * n ^ n;

  return
end


text:
function old = applytofig(h,varargin)
%APPLYTOFIG  Apply EXPORTFIG options to figure.
%   OLD = APPLYTOFIG(H) applies the default EXPORTFIG options to
%   the figure with handle H and returns the state-difference
%   structure in OLD.
%   OLD = APPLYTOFIG(H, OPTIONS) applies OPTIONS to H as described
%   in EXPORTFIG.
%   OLD = APPLYTOFIG(...,PARAM1,VAL1,PARAM2,VAL2,...) applies the
%   specified parameter-value pairs to H as described in EXPORTFIG.
%
%   Use RESTOREFIG(H,OLD) to restore the figure to its state before
%   applyfig was called.
%
%   See also EXPORTFIG, PREVIEWFIG, RESTOREFIG.

%  Copyright 2000-2009 The MathWorks, Inc

temp = tempname;
old = exportfig(h,temp, varargin{:},'applystyle',1);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [cvx_w, ATest, dSigPredict, dSigTest, rows, R2] = t_mrdArcuateXvalidate(dSig,A,rows,ndir)
%
% function [cvx_w, ATest, dSigTest, rows, R2] = t_mrdArcuateXvalidate(dSig,A,rows,ndir)
%
% Cross-validate fiber predictions.
%
% dSig = the vector of diffusion measurements
% A    = the fiber prediction matrix
% rows = which rows to use
%
% Example:
%    t_mrdArcuateXvalidate(dSig,A,rows,ndir)
%
% See also:  t_mrdTensors, t_mrdViewFibers, dwiLoad, dtiGet/Set,
%            t_mrdFiberPredictions, t_mrdArcuatePRedictions
%
% (c) Stanford VISTA Team

% The user must either pass in the rows to fit and hold out or the number
% of directions so that we can randomly hold out rows
if exist('rows','var') && ~isempty(rows);
    rows = logical(rows);
elseif ~exist('rows','var') || isempty(rows);
    % This is a randomly chosen direction to hold out for each voxel
    outVols = ceil(rand(length(dSig)./ndir,1).*ndir);
    % Now we will make a vector that has a 1 for each row of dsig to fit
    % and a 0 for each row to hold out for cross validation
    rows = [];
    for ii = 1:length(outVols)
        tmp = ones(ndir,1);
        tmp(outVols(ii))=0;
        nextrow = length(rows)+1;
        rows(nextrow:nextrow+ndir-1) = tmp;
    end
   rows = logical(rows);
end

% This is the data we will hold out for cross validation
tmp        = full(A);   % we're not sure if we have to make it a full before indexing
dSigTest   = dSig(~rows);
ATest      = tmp(~rows,:);

% now run the CVX code to solve the L1-minimization problem:
ATrain    = tmp(rows,:); clear tmp
n         = size(ATrain,2);
dSigTrain = dSig(rows);
fFraction = 0.2; % fraction of the weights over which weights are nromalized

l = 0;                 % Lower and upper bounds on the weights
u = 1;
cvx_solver sedumi;     % sdpt3
cvx_precision('low')  % We can handle low precision during testing.

cvx_begin              % start te cvx environment
   variable cvx_w(n)   % set the variable we are looking to fit in the cvx environment
   minimize(norm(ATrain * cvx_w - dSigTrain,1)) % minimize using L1 norm
   subject to
     norm(cvx_w(1:n),1) <= fFraction*n;
     cvx_w >= l;
     cvx_w <= u;
cvx_end

% compute the predicted signal
dSigPredict = ATest*cvx_w;

% Amount of deviation in the data explained by the model
% relative to the total variance in the data.
R2 = 100 * (1-sum((dSigPredict - dSigTest).^2) / sum((dSigTest-mean(dSigTest)).^2));
%R2 = corr(dSigPredict,dSigTest)^2;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef FEMPlotter < handle

    properties (Access = private)
        dim
        mesh
        displacement
        dispCoords
    end

    methods (Access = public)

        function obj = FEMPlotter(cParams)
            obj.init(cParams);
        end

        function plot(obj)
            ndimf = obj.dim.ndimField;
            switch ndimf
                case 2
                    obj.plotFem2D();
                case 3
                    obj.plotFem3D();
            end
        end

    end

    methods (Access = private)

        function init(obj, cParams)
            obj.dim          = cParams.dim;
            obj.mesh         = cParams.mesh;
            obj.displacement = cParams.displacement;
            obj.dispCoords   = obj.computeDisplacedCoords();
        end

        function dispCoords = computeDisplacedCoords(obj)
            ndimf = obj.dim.ndimField;
            ndof = obj.dim.ndof;
            coords = zeros(ndof,1);
            for i = 1:ndimf
                dofs = i:ndimf:ndof;
                coor = obj.mesh.coord(:,i);
                coords(dofs) = coor;
            end
            delta = obj.displacement;
            dispCoords = coords + delta;
        end

        function plotFem2D(obj)
            obj.plotNodes2D();
            obj.plotDisplacement2D();
        end

        function plotFem3D(obj)
            obj.plotNodes3D();
            obj.plotDisplacement3D();
        end

        function plotNodes2D(obj)
            Tn = obj.mesh.connec;
            x  = obj.mesh.coord(:,1);
            y  = obj.mesh.coord(:,2);
            figure()
            hold on
            colormap jet;
            plot(x(Tn)',y(Tn)','--','linewidth',0.5);
        end

        function plotDisplacement2D(obj)
            ndimf  = obj.dim.ndimField;
            ndof   = obj.dim.ndof;
            coords = obj.dispCoords;
            x = coords(1:ndimf:ndof);
            y = coords(2:ndimf:ndof);
            Tn   = obj.mesh.connec;
            plot(x(Tn)',y(Tn)','-k','linewidth',0.5);
        end

        function plotNodes3D(obj)
            Tn = obj.mesh.connec;
            x  = obj.mesh.coord(:,1);
            y  = obj.mesh.coord(:,2);
            z  = obj.mesh.coord(:,3);
            figure()
            hold on
            colormap jet;
            plot3(x(Tn)',y(Tn)',z(Tn)','--','linewidth',0.5);
        end

        function plotDisplacement3D(obj)
            ndimf  = obj.dim.ndimf;
            ndof   = obj.dim.ndof;
            coords = obj.dispCoords;
            x = coords(1:ndimf:ndof);
            y = coords(2:ndimf:ndof);
            z = coords(3:ndimf:ndof);
            Tn   = obj.mesh.connec;
            plot3(x(Tn)',y(Tn)',z(Tn)','-k','linewidth',0.5);
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


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

%%
clear; close all; clc;

%% Syntax
% |[V]=squircle(r,f,n);|

%% Description
% The squircle is a shape intermediate between a circle and a square. This
% function outputs the coordinates V for a squircle using the input radius
% r, the weighting parameter f, and the number of points n.
%
% See also: https://en.wikipedia.org/wiki/Squircle
%
% This implementation is based on:
% https://arxiv.org/vc/arxiv/papers/1604/1604.02174v1.pdf

%% Examples
%

%% Example 1: Constructing a squircle
n=200; %Number of points on curve
r=2; %Squircle radius

%%

f=0.5;
[V]=squircle(r,f,n);

%%

t=linspace(0,2*pi,n)';
Vc=r*[cos(t) sin(t) zeros(size(t))]; %Circle coordinates
Vs=r*[-1 -1 0; 1 -1 0; 1 1 0; -1 1 0; -1 -1 0];

%%
% Visualization

cFigure; hold on
hp1=plotV(Vc,'r-' ,'LineWidth',2);
hp2=plotV(Vs,'b-' ,'LineWidth',2);
hp3=plotV(V ,'k.-','LineWidth',5);

axis equal; axis tight; grid on; box on;
set(gca,'FontSize',25);
axis(r*[-1.25 1.25 -1.25 1.25]);
legend([hp1 hp2 hp3],{'Circle','Square',['Squircle ','f=',num2str(f)]},'Location','NorthEastOutSide')
drawnow;

%% Example 2: Animate effect of varying f

n=104; %Number of points on curve
r=1; %Squircle radius

nSteps=25; %Number of animation steps
f=linspace(0,1,nSteps);

%%

hf=cFigure; hold on
ht=gtitle(['f=',sprintf('%.3f',f(1))],25);
c=spectral(nSteps);
for q=1:1:nSteps
    [V]=squircle(r,f(q),n);
    hp1=plotV(V,'k-','LineWidth',2);
    hp1.Color=c(q,:);
end

hp2=plotV(V,'k.-','LineWidth',4,'MarkerSize',25);
axis equal; axis tight; grid on; box on;
colormap(c); colorbar; caxis([0 1]);
set(gca,'FontSize',25);
drawnow;

%%
%Populate the animation structure

%Create the time vector
animStruct.Time=linspace(0,1,nSteps);

for q=1:1:nSteps
    [V]=squircle(r,f(q),n);

    %Set entries in animation structure
    animStruct.Handles{q}=[hp2 hp2 ht]; %Handles of objects to animate
    animStruct.Props{q}={'XData','YData','String'}; %Properties of objects to animate
    animStruct.Set{q}={V(:,1),V(:,2),['f=',sprintf('%.3f',f(q))]}; %Property values for to set in order to animate
end

anim8(hf,animStruct);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function varargout = arg_guidialog(func,varargin)
% Create an input dialog that displays input fields for a Function and Parameters.
% Parameters = arg_guidialog(Function, Options...)
%
% The Parameters that are passed to the function can be used to override some of its defaults. The
% function must declare its arguments via arg_define. In addition, only a Subset of the function's
% specified arguments can be displayed.
%
% If this function is called with no arguments from within a function that uses arg_define it brings
% up a dialog for the function's defined arguments, with values assigned based on the function's
% inputs (contents of in varargin) and returns an argument struct if the dialog was confirmed with a
% click on OK, and an empty value otherwise. The calling function can exit on empty, or otherwise
% resume with entered parameters.
%
% In:
%   Function : the function for which to display arguments
%
%   Options... : optional name-value pairs; possible names are:
%                 'Parameters' : cell array of parameters to the Function to override some of its
%                                defaults.
%
%                 'Subset' : Cell array of argument names to which the dialog shall be restricted;
%                            these arguments may contain . notation to index into arg_sub and the
%                            selected branch(es) of arg_subswitch/arg_subtoggle specifiers. Empty
%                            cells show up in the dialog as empty rows.
%
%                 'Title' : title of the dialog (by default: functionname())
%
%                 'Invoke' : whether to invoke the function directly; in this case, the output
%                            arguments are those of the function (default: true, unless called in
%                            the form g = arg_guidialog; from within some function)
%
%                 'ShowGuru' : whether to show parameters marked as guru-level. (default: according
%                              to env_startup setting)
%
% Out:
%   Parameters : If Invoke was set to true, this is either the results of the function call or empty
%                (if the dialog was cancelled). If Invoke was set to false, this is either a struct
%                that is a valid input to the Function, or an empty value if cancel was clicked.
%
% Examples:
%   % bring up a configuration dialog for the given function
%   settings = arg_guidialog(@myfunction)
%
%   % bring up a config dialog with some pre-specified defaults
%   settings = arg_guidialog(@myfunction,'Parameters',{4,20,'test'})
%
%   % bring up a config dialog which displays only a subset of the function's arguments (in a particular order)
%   settings = arg_guidialog(@myfunction,'Subset',{'blah','xyz',[],'flag'})
%
%   % bring up a dialog, and invoke the function with the selected settings after the user clicks OK
%   settings = arg_guidialog(@myfunction,'Invoke',true)
%
%   % from within a function
%   args = arg_define(varargin, ...);
%   if nargin < 1
%       % optionally allow the user to enter arguments via a GUI
%       args = arg_guidialog;
%       % user clicked cancel?
%       if isempty(args), return; end
%   end
%
% See also:
%   arg_guidialog_ex, arg_guipanel, arg_define
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-10-24

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

global tracking;

if ~exist('func','var')
    % called with no arguments, from inside a function: open function dialog
    func = hlp_getcaller;
    varargin = {'Parameters',evalin('caller','varargin'),'Invoke',nargout==0};
end

if ischar(func)
    func = str2func(func); end

% parse arguments...
hlp_varargin2struct(varargin, ...
    {'params','Parameters','parameters'},{}, ...
    {'subset','Subset'},{}, ...
    {'dialogtitle','title','Title'}, [char(func) '()'], ...
    {'buttons','Buttons'},[], ...
    {'invoke','Invoke'},true, ...
    {'show_guru','ShowGuru'},tracking.gui.show_guru);
oldparams = params;

% obtain the argument specification for the function
rawspec = arg_report('rich', func, params); %#ok<*NODEF>

% extract a list of sub arguments...
if ~isempty(subset) && subset{1}==-1
    % user specified a set of items to *exclude*
    % convert subset to setdiff(all-arguments,subset)
    allnames = fieldnames(arg_tovals(rawspec));
    subset(1) = [];
    subset = allnames(~ismember(allnames,[subset 'arg_direct']));
end
[spec,subset] = obtain_items(rawspec,subset,'',show_guru);

% trim unnecessary blank margins at the beginning and end
while isempty(spec{1}) && isempty(subset{1})
    spec = spec(2:end);
    subset = subset(2:end);
end
while isempty(spec{end}) && isempty(subset{end})
    spec = spec(1:end-1);
    subset = subset(1:end-1);
end

% create an inputgui() dialog...
geometry = repmat({[0.6 0.35]},1,length(spec)+length(buttons)/2);
geomvert = ones(1,length(spec)+length(buttons)/2);
% turn the spec into a UI list...
uilist = {};
for k = 1:length(spec)

    s = spec{k};
    if isempty(s)
        uilist(end+1:end+2) = {{} {}};
    else
        if isempty(s.help)
            error(['Cannot display the argument ' subset{k} ' because it contains no description.']);
        else

            tag = subset{k};

            uilist{end+1} = {'Style','text', 'string',s.help{1}, 'fontweight','bold'};
            % depending on the type, we introduce different types of input widgets here...
            if iscell(s.range) && strcmp(s.type,'char')
                % string popup menu
                uilist{end+1} = {'Style','popupmenu', 'string',s.range,'value',find(strcmp(s.value,s.range)),'tag',tag};
            elseif strcmp(s.type,'logical')
                if iscell(s.range) % length(s.range)>1
                    % multiselect
                    uilist{end+1} = {'Style','listbox', 'string',s.range, 'value',find(ismember(s.range,s.value)),'tag',tag,'min',1,'max',100000};
                    geomvert(k) = min(3.5,length(s.range));
                else
                    % checkbox
                    uilist{end+1} = {'Style','checkbox', 'string','(set)', 'value',double(s.value),'tag',tag};
                end
            elseif strcmp(s.type,'char')
                % string edit
                uilist{end+1} = {'Style','edit', 'string', s.value,'tag',tag};
            else
                % expression edit
                if isinteger(s.value)
                    s.value = double(s.value); end
                uilist{end+1} = {'Style','edit', 'string', hlp_tostring(s.value),'tag',tag};
            end
            % append the tooltip string
            if length(s.help) > 1
                uilist{end} = [uilist{end} 'tooltipstring', regexprep(s.help{2},'\.\s+(?=[A-Z])','.\n')]; end
        end
    end

    if ~isempty(buttons) && k==buttons{1}
        % render a command button
        uilist(end+1:end+2) = {{} buttons{2}};
        buttons(1:2) = [];
    end
end

% invoke the GUI, obtaining a list of output values...
helptopic = char(func);
try
    if helptopic(1) == '@'
        fn = functions(func);
        tmp = struct2cell(fn.workspace{1});
        helptopic = class(tmp{1});
    end
catch
    disp('Cannot deduce help topic.');
end

[outs,dummy,okpressed] = inputgui('geometry',geometry, 'uilist',uilist,'helpcom',['env_doc ' helptopic], 'title',dialogtitle,'geomvert',geomvert); %#ok<ASGLU>

if ~isempty(okpressed)

    % remove blanks from the spec
    spec = spec(~cellfun('isempty',spec));
    subset = subset(~cellfun('isempty',subset));
    % turn the raw specification into a parameter struct (a non-direct one, since we will mess with
    % it)
    params = arg_tovals(rawspec);

    % for each parameter produced by the GUI...
    for k = 1:length(outs)
        s = spec{k}; % current specifier
        v = outs{k}; % current raw value
        % do type conversion according to spec
        if iscell(s.range) && strcmp(s.type,'char')
            v = s.range{v};
        elseif strcmp(s.type,'expression')
            v = eval(v);
        elseif strcmp(s.type,'logical')
            if length(s.range)>1
                v = s.range(v);
            else
                v = logical(v);
            end
        elseif strcmp(s.type,'char')
            % noting to do
        elseif ~isempty(v) && ischar(v)
            v = eval(v);
        end
        % assign the converted value to params struct...
        params = assign(params,subset{k},v);

    end

    % now send the result through the function to check for errors and obtain a values structure...
    params = arg_report('rich',func,{params});
    params = arg_tovals(params);

    % invoke the function, if so desired
    if ischar(func)
        func = str2func(func); end
    if invoke
        [varargout{1:nargout(func)}] = func(oldparams{:},params);
    else
        varargout = {params};
    end
else
    varargout = {[]};
end

% obtain a cell array of spec entries by name from the given specification
function [items,ids] = obtain_items(rawspec,requested,prefix,show_guru)
if ~exist('prefix','var')
    prefix = ''; end
items = {};
ids = {};
% determine what subset of (possibly nested) items is requested
if isempty(requested)
    % look for a special argument/property arg_dialogsel, which defines the standard dialog
    % representation for the given specification
    dialog_sel = find(cellfun(@(x)any(strcmp(x,'arg_dialogsel')),{rawspec.names}));
    if ~isempty(dialog_sel)
        requested = rawspec(dialog_sel).value; end
end
if isempty(requested)
    % empty means that all items are requested
    for k=1:length(rawspec)
        items{k} = rawspec(k);
        ids{k} = [prefix rawspec(k).names{1}];
    end
else
    % otherwise we need to obtain those items
    for k=1:length(requested)
        if ~isempty(requested{k})
            try
                [items{k},subid] = obtain(rawspec,requested{k});
            catch
                error(['The specified identifier (' prefix requested{k} ') could not be found in the function''s declared arguments.']);
            end
            ids{k} = [prefix subid];
        end
    end
end
% splice items that have children (recursively) into this list
for k = length(items):-1:1
    % special case: switch arguments are not spliced, but instead the argument that defines the
    % option popupmenu will be retained
    if ~isempty(items{k}) && ~isempty(items{k}.children) && (~iscellstr(items{k}.range) || isempty(requested))
        [subitems, subids] = obtain_items(items{k}.children,{},[ids{k} '.'],show_guru);
        if ~isempty(subitems)
            % and introduce blank rows around them
            items = [items(1:k-1) {{}} subitems {{}} items(k+1:end)];
            ids = [ids(1:k-1) {{}} subids {{}} ids(k+1:end)];
        end
    end
end

% remove items that cannot be displayed
retain = cellfun(@(x)isempty(x)||x.displayable&&(show_guru||~x.guru),items);
items = items(retain);
ids = ids(retain);

% remove double blank rows
empties = cellfun('isempty',ids);
items(empties(1:end-1) & empties(2:end)) = [];
ids(empties(1:end-1) & empties(2:end)) = [];

% obtain a spec entry by name from the given specification
function [item,id] = obtain(rawspec,identifier,prefix)
if ~exist('prefix','var')
    prefix = ''; end
% parse the . notation
dot = find(identifier=='.',1);
if ~isempty(dot)
    [head,rest] = deal(identifier(1:dot-1), identifier(dot+1:end));
else
    head = identifier;
    rest = [];
end
% search for the identifier at this level
names = {rawspec.names};
for k=1:length(names)
    if any(strcmp(names{k},head))
        % found a match!
        if isempty(rest)
            % return it
            item = rawspec(k);
            id = [prefix names{k}{1}];
        else
            % obtain the rest of the identifier
            [item,id] = obtain(rawspec(k).children,rest,[prefix names{k}{1} '.']);
        end
        return;
    end
end
error(['The given identifier (' head ') was not found among the function''s declared arguments.']);

% assign a field with dot notation in a struct
function s = assign(s,id,v)
% parse the . notation
dot = find(id=='.',1);
if ~isempty(dot)
    [head,rest] = deal(id(1:dot-1), id(dot+1:end));
    if ~isfield(s,head)
        s.(head) = struct(); end
    s.(head) = assign(s.(head),rest,v);
else
    s.(id) = v;
end

repo: REST


text:
function B=stdUnifRectBeamPattern(uv,numRows,numCols)
%%STDUNIFRECTBEAMPATTERN Obtain the value of the beampattern for a
%       rectangular array consisting of uniformly weighted isotropic
%       elements arranged in a rectangular grid and spaced a half-
%       wavelength apart. As defined in Chapter 2.2 of [1], the beampattern
%       is the frequency-wavenumber response function evaluated versus
%       direction. It describes the complex gain of the array to an input
%       plane wave. The values returned by this function are real.
%
%INPUTS: uv A 2XN set of N directions offset from the pointing direction of
%           the array in terms of direction cosines. sum(uv.^2,1) should
%           all be less than or equal to 1 for results to be valid. No
%           warnings or errors are given for values outside this range.
%   numRows The number of rows of element in the array.
%   numCols The number of columns of elements in the array.
%
%OUTPUTS: B The beampattern evaluated at the points in uv.
%
%The beam pattern of a rectangular array of uniformly weighted isotropic
%antennas is given in terms of sine functions in Equation 4.26 of [1]. The
%beampattern plays a roal in determining the gain in the radar range
%equation, which itself plays a role in determining the detection
%probability of a target.
%
%EXAMPLE:
%Here, we will plot the beam pattern magnitude power in decibels over the
%viewable region.
% numRows=40;
% numCols=20;
% numPoints=100;
% uVals=linspace(-1,1,numPoints);
% [u,v]=meshgrid(uVals,uVals);
% uv=[u(:).';v(:).'];
% B=stdUnifRectBeamPattern(uv,numRows,numCols);
% %Set the values of the response outside of the valid region to zero.
% sel=sum(uv.^2,1)>1;
% B(sel)=0;
% B=reshape(B,numPoints,numPoints);
%
% figure(1)
% clf
% surf(u,v,20*log10(abs(B)),'EdgeColor','None')
% axis([-1 1 -1 1 -50 0])
% caxis([-50 0])
% colormap(jet(256))
% view(20,45)
% light()
% h1=xlabel('u');
% h2=ylabel('v');
% h3=zlabel('Array Response');
% title('Array Power Response in Decibels')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h3,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] H. L. Van Trees, Optimum Array Processing. New York:
%    Wiley-Interscience, 2002.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=numRows;
N=numCols;

rowVals=(1/M)*sin(pi*(M/2)*uv(1,:))./sin((pi/2)*uv(1,:));
rowVals(~isfinite(rowVals))=1;%Remove singularity

colVals=(1/N)*sin(pi*(N/2)*uv(2,:))./sin((pi/2)*uv(2,:));
colVals(~isfinite(colVals))=1;%Remove singularity

B=rowVals.*colVals;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ fn ] = getBatchWrapper_Remat_rgbdiff(opts, numThreads, trainopts)

fn = @(imdb,batch, moreopts) getBatch(imdb,batch,opts,numThreads, trainopts, moreopts) ;
end

% -------------------------------------------------------------------------
function [im,labels] = getBatch(imdb, batch, opts, numThreads, trainopts, moreopts)
% -------------------------------------------------------------------------
opts.nFrames = imdb.images.nFrames(batch);

images = cell(1,numel(opts.nFrames));
labels = imdb.images.label(batch) ;
for k = 1:numel(opts.nFrames)
  images{k}  = strcat([imdb.images.name{batch(k)} filesep], images{k} ) ;
%   images{k}  = strcat([imdb.images.name{batch(k)}(1:end-4) filesep], images{k} ) ;
end
if ~isempty(moreopts) || nargin < 6
  for f = fieldnames(moreopts)'
    f = char(f) ;
    trainopts.(f) = moreopts.(f);
  end
end

opts.nFramesPerVid = trainopts.nFramesPerVid;
opts.numAugments = trainopts.numAugments;
opts.frameSample = trainopts.frameSample;
opts.augmentation =  trainopts.augmentation;

if isfield(trainopts, 'temporalStride')
  opts.temporalStride =  trainopts.temporalStride(randi(numel(trainopts.temporalStride))) ; % shuffle;
end
if isfield(trainopts, 'nFrameStack')
  opts.nFrameStack = trainopts.nFrameStack;
else
  opts.nFrameStack = trainopts.nFramesPerVid;
end
if isfield(trainopts, 'keepFramesDim')
  keepFramesDim = trainopts.keepFramesDim;
else
  keepFramesDim = false;
end
if isfield(trainopts, 'cheapResize')
  opts.cheapResize =  trainopts.cheapResize;
end

if isfield(trainopts, 'prefetch')
    opts.prefetch = trainopts.prefetch;
end
if isfield(trainopts, 'frameList')
    opts.frameList = trainopts.frameList;
end

im = cnn_ucf101_get_im_diff_batch(images, opts, ...
                            'numThreads', numThreads, ...
                          'flowDir', imdb.flowDir, ...
                          'imageDir', imdb.imageDir) ;

if iscell(im), return; end

if nargout == 1 % to work with dagnn & simplenn code
  if ndims(im) > 4 && ~keepFramesDim
    sz = size(im);
    nFrames = sz(5:end);
    im = permute(im, [1 2 3 5 4]);
    im = reshape(im, sz(1), sz(2), sz(3), []);
  end
  im = {'input', im(:,:,end-2:end,:,:), 'label', labels, 'input_flow', im(:,:,1:end-3,:,:)} ;
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function rule = subsrule(obj, name)
    % Return the string subsititution rule for the symbolic variable vector
    %
    % Return values:
    % rule: the string subsititution rules @type SymExpression

    assert(isempty(regexp(name, '\W', 'once')),...
        'SymExpression:invalidSymbol', ...
        'Invalid symbol name, can NOT contain special characters.');

    assert(isempty(regexp(name, '_', 'once')),...
        'SymExpression:invalidSymbol', ...
        'Invalid symbol name, can NOT contain ''_''.');

    if isscalar(obj)
        str = eval_math(['{' obj.s, '->HoldForm@Global`', name,'[[0]]}'],false);
        rule = SymExpression(str);
    elseif isvectorform(obj)
        siz = numel(obj);
        fstr= obj.s;
        str = eval_math(['((' fstr '[[#+1]]-> HoldForm@Global`', name '[[#]]&)/@(Range[', num2str(siz),']-1))']);
        rule = SymExpression(str,false);
    elseif ismatrix(obj)
        siz = numel(obj);
        fsym= tovector(obj);
        fstr = fsym.s;
        str = eval_math(['((' fstr '[[#+1]]-> HoldForm@Global`', name '[[#]]&)/@(Range[', num2str(siz),']-1))']);
        rule = SymExpression(str,false);
    else
        error('Unsupported format.');
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [obj_out, varargout] = image_math(obj1, varargin)
% Perform simple mathematical and boolean operations on image objects
%
% :Usage:
% ::
%
%    obj_out = image_math(obj1, [optional inputs, e.g., a 2nd object, keywords])
%
% For objects: Type methods(object_name) for a list of special commands
%              Type help object_name.method_name for help on specific
%              methods.
%
% :Inputs:
%
%   **obj1:**
%        An image_vector object
%
% :Optional Inputs:
%
%   **obj2:**
%        An additional image_vector object
%   **{'add', 'plus'}:**
%        Keyword to perform image-wise addition of images in obj1
%                              and obj2.  Assumes these are paired/matched objects.
%   **{'subtract', 'minus'}:**
%        Keyword to perform image-wise subtraction of images
%                              in obj1 and obj2
%   **{'cat', 'concatenate'}:**
%        Concatenate obj1 and obj2 image-wise.  Requires same
%                              number of voxels in both image sets.  Returns effects
%                              codes of 1, -1 in obj_out.Y.
%   **{'power'}:**
%        Keyword to raise data to power element-wise; obj.dat = obj.dat.^b;
%                              Followed by exponent to apply (b)
%
%   **{'rmssd'}:**
%        Keyword to perform root-mean-square successive differences.
%        Meaningful if image series is a timeseries.
%
%        [obj, rmssd, wh_outliers] = image_math(obj, 'rmssd');
%        Does not affect object.
%
% :Outputs:
%
%   **obj_out:**
%        The result - an image_vector object
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2015  Tor Wager
%
%     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/>.
% ..

% ..
% DEFAULTS AND INPUTS
% ..

keyword = '';  % initalize optional variables to default values here.
% keyword. should enter one...
obj2 = [];
my_exponent = [];

% optional inputs with default values
% -----------------------------------

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case {'add', 'plus'}, keyword = 'plus'; varargin{i} = [];
            case {'subtract', 'minus'}, keyword = 'minus'; varargin{i} = [];
            case {'cat', 'concatenate'}, keyword = 'cat'; varargin{i} = [];

            case 'power', keyword = 'power'; varargin{i+1} = my_exponent;
                %case 'basistype', basistype = varargin{i+1}; varargin{i+1} = [];

            case 'rmssd', keyword = 'rmssd'; varargin{i} = [];

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end

    elseif isa(varargin{i}, 'image_vector')

        obj2 = varargin{i}; varargin{i} = [];

    end
end

% -------------------------------------------------------------------------
% DATA CHECKS
% -------------------------------------------------------------------------

switch keyword

    case ''
        disp('Nothing to do.')
        return

    case 'power'
        %Check if image_vector object
        if ~isa(obj1,'image_vector')
            error('Input Data is not an image_vector object')
        end

    case {'plus', 'minus', 'cat'}
        n1 = size(obj1.dat, 2);
        n2 = size(obj2.dat, 2);

        y1 = ones(n1, 1);
        y2 = ones(n2, 1);

        %Check if image_vector object
        if ~isa(obj1,'image_vector') || ~isa(obj2,'image_vector')
            error('Input Data is not an image_vector object')
        end

        %Check number of rows
        if size(obj1.dat,1) ~= size(obj2.dat,1)
            % Voxel list mismatch
            % Try to fix

            obj1 = replace_empty(obj1);
            obj2 = replace_empty(obj2);

            if size(obj1.dat,1) ~= size(obj2.dat,1)
                % If still different, give a warning and resample.

                warning('Image_math: Objects being operated on do not appear to be in the same space. Resampling 2nd object, but proceed with caution!!');

                obj2 = resample_space(obj2, obj1);
            end

            %error('number of voxels is different between objects.')
        end

        if strcmp(keyword, 'cat')
            % we are done

        else

            % The rest is for plus/minus, which assume paired images

            %Check number of columns
            if n1 ~= n2
                error('Sizes of objects do not match.');
            end

        end

    case 'rmssd'

    otherwise
        warning(['Unknown keyword:' keyword]);
        return

end % keyword ; data checks

% -------------------------------------------------------------------------
% RUN OPERATION
% -------------------------------------------------------------------------

switch keyword

    case ''
        disp('Nothing to do.')
        return

    case 'plus'

        obj_out = obj1;
        obj_out.dat = obj_out.dat + obj2.dat;
        obj_out.history{end+1} = 'Image-wise addition operation by image_math';
        obj_out.dat_descrip = cell(1, 3);
        obj_out.dat_descrip{1} = 'Names of images added in next cells, 1st set plus 2nd';
        obj_out.dat_descrip{2} = obj1.fullpath;
        obj_out.dat_descrip{3} = obj2.fullpath;
        obj_out.image_names = [];
        obj_out.fullpath = [];
        obj_out.files_exist = false;

    case 'minus'

        obj_out = obj1;
        obj_out.dat = obj_out.dat - obj2.dat;
        obj_out.history{end+1} = 'Image-wise subtraction operation by image_math';
        obj_out.dat_descrip = cell(1, 3);
        obj_out.dat_descrip{1} = 'Names of images subtracted in next cells, 1st set minus 2nd';
        obj_out.dat_descrip{2} = obj1.fullpath;
        obj_out.dat_descrip{3} = obj2.fullpath;
        obj_out.image_names = [];
        obj_out.fullpath = [];
        obj_out.files_exist = false;

    case 'cat'

        obj_out = obj1;
        obj_out.dat = [obj_out.dat obj2.dat];
        obj_out.history{end+1} = 'Image-wise concatenation operation by image_math';
        obj_out.image_names = char(obj1.image_names, obj2.image_names);
        obj_out.fullpath = char(obj1.fullpath, obj2.fullpath);
        if isfield(obj_out, 'Y')
            obj_out.Y = [y1; -y2];
            obj_out.Y_descrip = 'Effects codes for image set A (1) and B (-1) added by image_math';
        end
    case 'power'
        obj_out = obj1;
        obj_out.dat = obj_out.dat .^ my_exponent;
        obj_out.history{end+1} = sprintf('Raised to %3.0f element-wise by image_math', my_exponent);

    case 'rmssd'

        obj_out = obj1;

        sdiffs = diff(obj1.dat')';
        sdiffs = [mean(sdiffs, 2) sdiffs]; % keep in image order
        rmssd = ( mean(sdiffs .^ 2) ) .^ .5; % rmssd - root mean square successive diffs

        % avoid first time point being very different and influencing distribution and plots.
        rmssd(1) = median(rmssd);

        % z-scores of rmssd
        wh_outliers = scale(rmssd) > 3;

        varargout{1} = rmssd;
        varargout{2} = wh_outliers;

    otherwise
        warning(['Unknown keyword:' keyword]);
        return

end % switch keyword; operation

end % main function


text:
function u = ic_test01 ( x_num, x, t )

%*****************************************************************************80
%
%% IC_TEST01 sets the initial condition for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, the number of nodes.
%
%    Input, real X(X_NUM,1), the node coordinates.
%
%    Input, real T, the time.
%
%    Output, real U(X_NUM,1), the initial value of U.
%
  u = zeros ( x_num, 1 );

  u(1:x_num,1) = 50.0;

  return
end


text:
%--- help for prepare_transition_routine ---
%
%  INTERNAL FUNCTION
%
%
%    Other functions named prepare_transition_routine
%
%       dsge/prepare_transition_routine
%       generic/prepare_transition_routine
%

author: jmaih


text:
function  Theta0 = InitialPara_random_MarkovDep_NN(X,K)
if nargin<2
    K=150;
end
[P,N]=size(X);

%--------Initialize Parameter D,S,Z,Pi,gamma_epsi,gamma_s-----------------
D=zeros(P,K);
for k=1:K
    D(:,k)=randn(P,1)*sqrt(1/P);
end

S=zeros(K,N);
for n=1:N
    S(:,n)=randn(K,1);
end

Z = zeros(K,1);
%Z = ones(K,1);
Delta = ones(K,1);
Tao = ones(K,1);
Pi = 0.5*ones(K,1);
gamma_epsi = 1e3*ones(1,N);

%sampe sparse component
gamma_s = 1;
S2 = randn(P,N);
Z2 = zeros(P,N);
Pi2 = 0.5*ones(P,N);

Theta0.D = D;
Theta0.S = S;
Theta0.Z = Z;
Theta0.Delta = Delta;
Theta0.Tao = Tao;
Theta0.Pi = Pi;
Theta0.gamma_epsi = gamma_epsi;
Theta0.S2 = S2;
Theta0.Z2 = Z2;
Theta0.gamma_s = gamma_s;
Theta0.Pi2 = Pi2;
end

author: andrewssobral

repo: lrslibrary


text:
function y2 = y_astronomical_to_common ( y )

%*****************************************************************************80
%
%% Y_ASTRONOMICAL_TO_COMMON converts an Astronomical year to a Common year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the astronomical year.
%
%    Output, integer Y2, the Common year.
%
  if ( y <= 0 )
    y2 = y - 1;
  else
    y2 = y;
  end

  return
end

author: johannesgerer


text:
function out = DSLInit(startPos, endPos, map, scalling)

%% generate shape pattern
    radius = 10/scalling;
    mat = zeros(radius*2,radius*2);
    dista = @(a,b) sqrt(a*a+b*b);

    shapePattern = [];

    for x = 1:1:2*radius
        for y =1:1:2*radius
            if dista(x-radius-0.5,y-radius-0.5) > radius -1 && dista(x-radius-0.5,y-radius-0.5) < radius
                mat(x,y) = 1;
                shapePattern = [shapePattern; floor(x-radius) floor(y-radius) 0 0];
            end
        end
    end

    neighbours = [
        0 1 0 0;
        -1 0 0 0;
        0 -1 0 0;
        1 1 0 0;
        -1 -1 0 0;
        1 0 0 0;
        1 -1 0 0;
        -1 1 0 0;
    ]';

%% prepare all data

		out.map = map;
    out.startPos = [startPos(2:-1:1); 0; 0];
    out.endPos = [endPos(2:-1:1); 0; 0];
    startPos = out.startPos;

    out.sacalling = scalling;
    out.pattern = shapePattern';
    out.ucc = neighbours;
    out.height = ceil(length(out.map(:,1)));
    out.width = ceil(length(out.map(1,:)));
    out.graph = zeros(out.height, out.width,5);
    out.graph(:,:,1:2) = inf;
    out.graph(:,:,3) = false;
    out.graph(:,:,5) = -1;
    out.kM = 0;
    SQRT2 = sqrt(2)-1;
    out.comparator = DStarcmp;
    out.stack = java.util.PriorityQueue(180247, out.comparator);

    setRhs(out.endPos, 0);
    setQ(out.endPos);
    out.endPos(3:4) = [heur(out.endPos); 0];
    add(out.stack, out.endPos);

    %-----------------------------------------------------------
    function setQ(s)
         out.graph(s(1), s(2), 3) = true;
    end
    function setg(s, val)
        out.graph(s(1), s(2),1) = val;
    end
    function setRhs(s, val)
        out.graph(s(1), s(2),2) = val;
    end
    function out = heur(s)
        k = abs(startPos - s);
        out = SQRT2*min(k(1:2)) + max(k(1:2));
    end

end

author: LazyFalcon


text:
function [pvec, pstruct] = tapas_hgf_ar1_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pvec    = NaN(1,length(ptrans));
pstruct = struct;

l = r.c_prc.n_levels;

pvec(1:l)         = ptrans(1:l);                           % mu_0
pstruct.mu_0      = pvec(1:l);
pvec(l+1:2*l)     = exp(ptrans(l+1:2*l));                  % sa_0
pstruct.sa_0      = pvec(l+1:2*l);
pvec(2*l+1:3*l)   = tapas_sgm(ptrans(2*l+1:3*l),1);        % phi
pstruct.phi       = pvec(2*l+1:3*l);
pvec(3*l+1:4*l)   = ptrans(3*l+1:4*l);                     % m
pstruct.m         = pvec(3*l+1:4*l);
pvec(4*l+1:5*l-1) = exp(ptrans(4*l+1:5*l-1));              % ka
pstruct.ka        = pvec(4*l+1:5*l-1);
pvec(5*l:6*l-1)   = ptrans(5*l:6*l-1);                     % om
pstruct.om        = pvec(5*l:6*l-1);
pvec(6*l)         = exp(ptrans(6*l));                      % al
pstruct.al        = pvec(6*l);

return;

author: translationalneuromodeling

repo: tapas


text:
% The example shows how to transforms between different conventions
%  of Euler angles. In particular from ZYZ to ZYX
%
%
%
%
function test_transform_euler_angles()
% initial
convention1 = 'ZYZ'
alpha1 = pi/4;
beta1 = pi/4;
gamma1 = pi/4;
euler1 = [alpha1 beta1 gamma1]

R1 = Rot(alpha1, 'z')*Rot(beta1, 'y')*Rot(gamma1, 'z')

convention2 = 'ZYX'
[euler2_1, euler2_2] = rot2euler(R1,convention2)

% check that yields the same rotation matrix
R2_1 = Rot(euler2_1(1), 'z')*Rot(euler2_1(2), 'y')*Rot(euler2_1(3), 'x')
R2_2 = Rot(euler2_2(1), 'z')*Rot(euler2_2(2), 'y')*Rot(euler2_2(3), 'x')

norm(R1-R2_1)
norm(R1-R2_2)

repo: ARTE


text:
function [CTmaterial materialNames] = generateMaterialMap(CTdensity, materialMap);

% JC Feb 2007
% Assign material numbers based on the density (g/cm3).
% Use to show material numbers/names in CERR GUI.

% Usage:
% [CTmaterial] = generateMaterialMap(materialMap, planC, scanNumber);

% output:
    % CTmaterial := same size as the uniformized CT scan, containing
    % material numbers.
    % materialNames := defined material names: eg. 'Air'    'Lung'
    % 'Water'    'Bone'

% input:
% >> materialMap(1)
% ans =
%           name: 'Air'
%     minDensity: 0
%     maxDensity: 0.0500

CTmaterial = zeros(size(CTdensity));
materialNames = cell(size(materialMap));

for i = 1 : length(materialMap),
    CTmaterial(materialMap(i).minDensity <= CTdensity & CTdensity < materialMap(i).maxDensity) = i;
    materialNames{i} = getfield(materialMap, {1,i}, 'name');
end

% If CTmaterial still has "0" elements, it means that the "materialMap"
% doesn't cover the whole CT values.
x = ismember(unique(CTmaterial), [1:length(materialMap)]);
if ismember(0,x)
warning('Input materialMap does NOT cover the whole CT scan values. Please check the CT values for voxels with "0" material number')
end

return;


text:
function [Tmin,Tsec] = sec2min(T)

% Transforms the given number of seconds into minutes and residual seconds

Tmin = floor(T/60);
Tsec = round((T-Tmin*60)*10)/10;

author: InverseTampere

repo: TreeQSM


text:
function demoImMaxGeodesicPath(varargin)
%DEMOIMGEODESICPATH  Demo for function imGeodesicPath
%
%   output = demoImGeodesicPath(input)
%
%   Example
%   demoImGeodesicPath
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-02-22,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

%% Initialisation de l'image a traiter

% lecture image
img = imread('rice.png');

% segmentation basique
img2 = img - imopen(img, ones(30, 30));
bin =img2 > 50;

% etiquetage
lbl = bwlabel(bin, 4);

figure;
imshow(label2rgb(lbl, 'jet', 'w'));

%% Calcule des chemins geodesiques maximaux

% affiche image
imshow(img);
hold on;

% pour chaque particule, calcule et affiche le chemin
tic;
for i=1:max(lbl(:))
    path = imMaxGeodesicPath(lbl == i);
    plot(path(:,1), path(:,2), 'color', 'r', 'linewidth', 2);
end

% affiche le temps ecoule
t = toc;
disp(sprintf('Temps ecoule : %f', t)); %#ok<DSPS>

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function success=testConvertClassifierWithScoreFeature()

testDataDirName='/groups/branson/bransonlab/projects/JAABA/test_data';
%testDataDirName='/Users/taylora/jaaba/test_data';

projectFileName=fullfile(testDataDirName,'metafooing_project.mat');
classifierFileName=fullfile(testDataDirName,'metafooing_labels_and_classifier.mat');
gtExpDirNames={};
jabFileName='/tmp/metafooing_labels_and_classifier.jab';

nameOfExpInMetafooingClassifier='GMR_71G01_AE_01_TrpA_Rig2Plate14BowlD_20110707T154929';

% Make sure the fooing score are _not_ already in the perframe dir
scorePerframeFileNameAbs= ...
  sprintf('%s/%s/perframe/scores_fooing.mat', ...
          testDataDirName, ...
          nameOfExpInMetafooingClassifier);
if exist(scorePerframeFileNameAbs,'file') ,
  cmd=sprintf('rm -f "%s"',scorePerframeFileNameAbs);
  system(cmd);
end

% Also deleted the converted .jab file if it exsits
if exist(jabFileName,'file') ,
  cmd=sprintf('rm -f "%s"',jabFileName);
  system(cmd);
end

% Convert the old-style classifier
gum=GrandlyUnifyModel();
gum.setProjectFileName(projectFileName);
gum.setClassifierFileName(classifierFileName);
for i=1:length(gtExpDirNames)
  gum.addGTExpDirName(gtExpDirNames{i})
end
gum.convert(jabFileName);
gum=[];  %#ok

% Open the converted file
gtMode=false;
data=JLabelData('setstatusfn',@(str)(fprintf('%s\n',str)), ...
                'clearstatusfn',@()(nop()));
data.openJabFile(jabFileName,gtMode);

% Display the JLabelData object
data

% Do a few spot-checks
labels=data.labels;
if length(labels(1).t0s{1}) ~= 2 ,
  data.closeJabFile();
  error('.jab file has the wrong number of labels');
end
if abs(labels(1).timestamp{1}(2)-735364.459128824) > 1e-6 ,
  data.closeJabFile();
  error('the one timestamp I checked is wrong');
end

% Clean-up
data.closeJabFile();
success=true;

end

author: kristinbranson

repo: JAABA


text:
% PWD2PATH()
%
% add current working directory (pwd) to path.

P = path;
path(pwd,P);

repo: SelfExSR


text:
function display(ipa)
% INPLACE/DISPLAY Display an inplace array from the command line.
%
% Example:
%    ipa = inplace(ones(5));
%    ipa

disp([inputname(1), ' = ']);
disp(ipa.get_a());


text:
function [sys,x0,str,ts]=Plant1_fhan(t,x,u,flag)

switch flag,
    case 0,
        [sys,x0,str,ts]=mdlInitializeSizes;
    case 1,
        sys=mdlDerivatives(t,x,u);
    case 3,
        sys=mdlOutputs(x);
    case {2,4,9},
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes
    sizes=simsizes;
    sizes.NumContStates=2;
    sizes.NumDiscStates=0;
    sizes.NumOutputs=1;
    sizes.NumInputs=1;
    sizes.DirFeedthrough=1;
    sizes.NumSampleTimes=1;
    sys=simsizes(sizes);
    x0=[0;0];
    str=[];
    ts=[0 0];
function sys=mdlDerivatives(t,x,u)
    sys(1)=x(2);
    sys(2)=cos(0.6*t)*x(1)+cos(0.7*t)*x(2)+u;
function sys=mdlOutputs(x)
    sys=x(1);

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function demo()
%
% demonstration file for NMFLibrary.
%
% This file illustrates how to use this library.
% This demonstrates Frobenius-norm based multiplicative updates (MU) algorithm and
% hierarchical alternative least squares (Hierarchical ALS) algorithm.
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on Apr. 05, 2017

    clc;
    clear;
    close all;

    %% generate synthetic data of (mxn) matrix
    m = 500;
    n = 100;
    V = rand(m,n);

    %% Initialize of rank to be factorized
    rank = 5;

    %% perform factroization
    options.verbose = 1;
    % MU
    options.alg = 'mu';
    [w_mu, infos_mu] = fro_mu_nmf(V, rank, options);
    % Hierarchical ALS
    options.alg = 'hals';
    [w_hals, infos_hals] = als_nmf(V, rank, options);
    % Accelerated Hierarchical ALS
    options.alg = 'acc_hals';
    [w_acchals, infos_acchals] = als_nmf(V, rank, options);

    %% plot
    display_graph('epoch','cost', {'Fro-MU', 'HALS', 'Acc-HALS'}, {w_mu, w_hals, w_acchals}, {infos_mu, infos_hals, infos_acchals});
    display_graph('time','cost', {'Fro-MU', 'HALS', 'Acc-HALS'}, {w_mu, w_hals, w_acchals}, {infos_mu, infos_hals, infos_acchals});

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function deep_attention_test()
close all;
clc;
% clear mex;
clear is_valid_handle; % to clear init_key
run(fullfile(fileparts(mfilename('fullpath')), 'startup'));
opts.vocal = true;

%% -------------------- CONFIG --------------------
opts.dataset{1} = 'pascal';
%opts.dataset{2} = 'DUT';
%opts.dataset{3} = 'MIT300';
%opts.dataset{4} = 'MIT1003';
%opts.dataset{5} = 'Toronto';
opts.caffe_version          = 'caffe_faster_rcnn';
opts.gpu_id                 = auto_select_gpu;
active_caffe_mex(opts.gpu_id, opts.caffe_version);
opts.use_gpu                = true;

attention_model_dir              = fullfile(pwd, 'models'); %% VGG
attention_model                  = load_model(attention_model_dir);

attention_net = caffe.Net(attention_model.attention_net_def, 'test');
attention_net.copy_from(attention_model.attention_net);
%attention_net.params('final_attention_pred',1).get_data()
% if opts.use_gpu
%     attention_model.image_means = gpuArray(attention_model.image_means);
% end

% set gpu/cpu
if opts.use_gpu
    caffe.set_mode_gpu();
else
    caffe.set_mode_cpu();
end

for kk = 1:length(opts.dataset)
    opts.attention_dir = fullfile(pwd, 'datasets', 'result', opts.dataset{kk});
    mkdir_if_missing(opts.attention_dir);

    files=dir(['datasets/' opts.dataset{kk} '/*.jpg']);
    files=struct2cell(files)' ;
    image_names=files(:,1);

    for imnum = 1:length(image_names)
        imnum
        image = imread(['datasets/' opts.dataset{kk} '/' image_names{imnum}]);
        [w,h,c] = size(image);
%       if opts.use_gpu
%           im = gpuArray(image);
%       end
        [~, ~, ~, final_attentionmap] = fixationmap_detect(attention_model.conf, attention_net, image);
        final_attentionmap  = imresize(final_attentionmap,[w,h]);
        imwrite(final_attentionmap,[opts.attention_dir '/'  image_names{imnum}(1:end-4) '.jpg']);
    end
end
caffe.reset_all();
clear mex;

end

function model = load_model(model_dir)
    ld             = load(fullfile(model_dir, 'model'));
    model = ld.attention_model;
    clear ld;
    %% load attention model
    model.attention_net_def ...
                                = fullfile(model_dir, model.attention_net_def);
    model.attention_net ...
                                = fullfile(model_dir, model.attention_net);
end

author: wenguanwang

repo: deepattention


text:
%
% y = ndetrend(x,dim,type,breaks)
%
% NDETREND removes a trend from an N-dimensional array <x> along a given
% dimension <dim>.
%
%       x: array to be detrended
%     dim: dimension along which to detrend, default is first nonsingleton
%    type: linear detrend (1) or demean (0), default is linear.
%  breaks: breakpoint indices for a piecewise linear trend
%
% See also DETREND, MEAN

%   Based on DETREND.M
%   Copyright 1984-2004 The MathWorks, Inc.
%   Modified by Bill Winter May 2006
function x = ndetrend(x,dim,o,b)
siz = size(x);                                      % array size
if nargin < 2, dim = find(siz > 1,1); end           % default: nonsingleton
N = siz(dim);                                       % dimension length
a = ones(N,1);                                      % constant
if nargin < 3 || o                                  % default: linear
    if nargin < 4, b = []; end                      % default: no breaks
    b = unique([1;b(:)]);                           % breaks unique
    b(b < 1 | b > N-1) = [];                        % breaks within array
    l = length(b);                                  % number linear pieces
    a = [zeros(N,l) a];                             % preallocate linear
    M = N - b;                                      % length of linear piece
    for k = 1:l, a(1+b(k):N,k) = (1:M(k))'/M(k); end% linear pieces
end
if length(siz) > 2
    if dim ~= 1                                     % permute 1 with dim
        per = 1:length(siz);
        per([1 dim]) = [dim 1];
        siz([1 dim]) = siz([dim 1]);
        x = builtin('permute',x,per);
    end
    for k = 1:prod(siz(3:end)), x(:,:,k) = x(:,:,k) - a*(a\x(:,:,k));end
    if dim ~= 1, x = builtin('permute',x,per);end   % depermute 1 with dim
elseif dim == 2, x = x - (x/a')*a';                 % right-regress
elseif dim == 1, x = x - a*(a\x);                   % left-regress
end

author: Sable

repo: mcbench-benchmarks


text:
function [ x, w ] = ccfi_1 ( n, ell )

%*****************************************************************************80
%
%% CCFI_1 returns a Boyd quadrature rule for the Laguerre integral.
%
%  Discussion:
%
%    The Laguerre integral I(f) is:
%      I(f) = integral ( 0 <= x < +oo ) f(x) dx
%    and the quadrature rule which approximates I(f) is
%      Q(f) = sum ( 1 <= i <= n ) w(i) * f(x(i))
%
%    The parameter ELL controls the mapping between [-1,+1] and [0,+oo),
%    and can be initially chosen to be 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John Boyd,
%    Exponentially convergent Fourier-Chebyshev quadrature schemes on
%    bounded and infinite intervals,
%    Journal of Scientific Computing,
%    Volume 2, Number 2, 1987, pages 99-109.
%
%  Parameters:
%
%    Input, integer N, the number of points to use in the rule.
%
%    Input, real ELL, the mapping parameter.
%
%    Output, real X(N), W(N), the points and weights for the quadrature rule.
%
  t = ( pi * ( 1 : n ) / ( n + 1 ) )';
  x = ell * ( cot ( 0.5 * t ) ) .^ 2;

  w = zeros ( n, 1 );
  for i = 1 : n
    for j = 1 : n
      w(i) = w(i) + sin ( j * t(i) ) * ( 1.0 - cos ( j * pi ) ) / j;
    end
  end
  w = w * 2.0 * ell .* sin ( t ) ./ ( 1.0 - cos ( t ) ) .^ 2 * 2.0 / ( n + 1 );

  return
end

repo: jburkardt-m