text:
function quality=meshquality(node,elem)
%
% quality=meshquality(node,elem)
%
% compute the Joe-Liu mesh quality measure of a tetrahedral mesh
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2011/02/26
%
% input:
%    node:  node coordinates of the mesh (nn x 3)
%    elem:  element table of a tetrahedral mesh (ne x 4)
%
% output:
%    quality: a vector of the same length as size(elem,1), with
%           each element being the Joe-Liu mesh quality metric (0-1) of
%           the corresponding element. A value close to 1 represents
%           higher mesh quality (1 means equilateral tetrahedron);
%           a value close to 0 means nearly degenerated element.
%
% reference:
%    A. Liu, B. Joe, Relationship between tetrahedron shape measures,
%                    BIT 34 (2) (1994) 268-287.
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(size(elem,2)>4)
    elem=elem(:,1:4);
end
enum=size(elem,1);
vol=elemvolume(node,elem);
edges=meshedge(elem);
ed=node(edges(:,1),:)-node(edges(:,2),:);
ed=sum((ed.*ed)');
ed=sum(reshape(ed,[enum length(ed)/enum])')';

quality=12*((3*vol).^(2/3))./ed;

repo: REST


text:
function x=nsst_rec2(dst,shear_f,lpfilt)
% This function performs the inverse (local) nonsubsampled shearlet transform as given
% in G. Easley, D. Labate and W. Lim, "Sparse Directional Image Representations
% using the Discrete Shearlet Transform", Appl. Comput. Harmon. Anal. 25 pp.
% 25-46, (2008).
%
% Input:
%
% dst          - the nonsubsampled shearlet coefficients
%
% shear_f      - the cell array containing the shearing filters
%
% lpfilt       - the filter to be used for the Laplacian
%                Pyramid/ATrous decomposition using the codes
%                written by Arthur L. Cunha
%
% Output
%
% x         - the reconstructed image
%
% Code contributors: Glenn R. Easley, Demetrio Labate, and Wang-Q Lim.
% Copyright 2011 by Glenn R. Easley. All Rights Reserved.
%

level=length(dst)-1;
y{1}=dst{1};
for i=1:level,
      l=size(dst{i+1},3);
      for k=1:l,
          dst{i+1}(:,:,k)=conv2p(conj(shear_f{i}(:,:,k)),dst{i+1}(:,:,k));
      end
      y{i+1} = real(sum(dst{i+1},3));
end

x=real(atrousrec(y,lpfilt));

author: thomaskuestner


text:
function [atm] = Pa2atm(Pa)
% Convert units of pressure from pascals to atmospheres.
% Chad A Greene 2012
atm = Pa/101325;

author: Sable

repo: mcbench-benchmarks


text:
function item = update_deps(item, varargin)

% function item = update_deps(item, varargin)
% This function will run cfg_dep/update_deps in all leaf (cfg_entry,
% cfg_menu, cfg_files) nodes of a configuration tree and update their
% dependency information (mod_job_ids) if necessary.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: update_deps.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok

if ~isempty(item.val) && isa(item.val{1}, 'cfg_dep')
    item.val{1} = update_deps(item.val{1},varargin{:});
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%
%  Another demo on the Jacobian for the Sawyer robot.
% JACOBIAN DEMO

% Copyright (C) 2017, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE 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 ARTE.  If not, see <http://www.gnu.odrg/licenses/>.
close all;

fprintf('\nThe demo shows how to compute the end effectors speed as a function of the joint speeds and viceversa')
%
%robot=load_robot('RETHINK','SAWYER');

%First compute the linear and angular speeds of the end effector given the
%joint speeds qd = [2 2 2 2 2 2 2] rad/s.
% Compute at joint position q = [0.1 0.1 0.1 0.1 0.1 0.1] rad
q = [-0.15 0.0 0.0 0 0.0 0.0 0.0]';
%qd_1 = [2 2 2 2 2 2 2]';

qd_1 = [0 1 0 0 0 0 0 ]';

drawrobot3d(robot, q)

% The result V is  [Vx Vy Vz Wx Wy Wz] [m/s m/s m/s rad/s rad/s rad/s]
% direct jacobian Here
V_1 = compute_end_velocity(robot, q, qd_1)

% Now, at the same position q, we would like to compute the join speeds qd
% that will bring the end effector to the speed V. The result should match
% the values in qd as defined earlier
qd_2 = compute_joint_velocity(robot, q, V_1)

V_2 = compute_end_velocity(robot, q, qd_2)

%cambia q de manera que qd_1 y qd_2 sean diferentes!

%SINGULARITIES
% The abb irb140 shows a singularity point at q = [0 0 0 0 0 0]. This means
% that the manipulator Jacobian J cannot be inverted (det(J)=0). If we
% repeat the former instructions for q = [0 0 0 0 0 0] we get a Warning
% indicating that J is badly conditioned and an unaccurate result.
% (Uncomment the following 4 lines to test this case)
%q = [0 0 0 0 0 0];
%qd = [1 1 1 1 1 1];
%V = compute_end_velocity(robot, q, qd)
%qd = compute_joint_velocity(robot, q, V)

repo: ARTE


text:
% RES = corrDn(IM, FILT, EDGES, STEP, START, STOP)
%
% Compute correlation of matrices IM with FILT, followed by
% downsampling.  These arguments should be 1D or 2D matrices, and IM
% must be larger (in both dimensions) than FILT.  The origin of filt
% is assumed to be floor(size(filt)/2)+1.
%
% EDGES is a string determining boundary handling:
%    'circular' - Circular convolution
%    'reflect1' - Reflect about the edge pixels
%    'reflect2' - Reflect, doubling the edge pixels
%    'repeat'   - Repeat the edge pixels
%    'zero'     - Assume values of zero outside image boundary
%    'extend'   - Reflect and invert
%    'dont-compute' - Zero output when filter overhangs input boundaries
%
% Downsampling factors are determined by STEP (optional, default=[1 1]),
% which should be a 2-vector [y,x].
%
% The window over which the convolution occurs is specfied by START
% (optional, default=[1,1], and STOP (optional, default=size(IM)).
%
% NOTE: this operation corresponds to multiplication of a signal
% vector by a matrix whose rows contain copies of the FILT shifted by
% multiples of STEP.  See upConv.m for the operation corresponding to
% the transpose of this matrix.

% Eero Simoncelli, 6/96, revised 2/97.

function res = corrDn(im, filt, edges, step, start, stop)

%% NOTE: THIS CODE IS NOT ACTUALLY USED! (MEX FILE IS CALLED INSTEAD)

fprintf(1,'Warning: You should compile the MEX code for "corrDn", found in the MEX subdirectory.  It is MUCH faster.\n');

%------------------------------------------------------------
%% OPTIONAL ARGS:

if (exist('edges') == 1)
  if (strcmp(edges,'reflect1') ~= 1)
    warning('Using REFLECT1 edge-handling (use MEX code for other options).');
  end
end

if (exist('step') ~= 1)
	step = [1,1];
end

if (exist('start') ~= 1)
	start = [1,1];
end

if (exist('stop') ~= 1)
	stop = size(im);
end

%------------------------------------------------------------

% Reverse order of taps in filt, to do correlation instead of convolution
filt = filt(size(filt,1):-1:1,size(filt,2):-1:1);

tmp = rconv2(im,filt);
res = tmp(start(1):step(1):stop(1),start(2):step(2):stop(2));

author: ricedsp


text:
function z = comp_fftanalytic(s)
%COMP_FFTANALYTIC Compute analytic representation
%
%   Usage: z = comp_fftanalytic(s);
%
%   `comp_fftanalytic(s)` computes the analytic representation of s.
%   The analytic representation is computed through the FFT of f.
%

% AUTHOR: Jordy van Velthoven

Ls = size(s,1);
H = floor(Ls/2);

z = fft(s);
z(2:Ls-H,:) = 2*z(2:Ls-H,:);
z(H+2:Ls,:) = 0;
z = ifft(z);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function deri = deri_inherit(deri_cost)
    global config mem;
    deri = deri_cost .* config.DERI_NONLINEARITY(mem.output) / config.batch_size;
end

author: jimmy-ren


text:
function [projectedData, mapping] = ml_projection(X,options)
%ML_PROJECTION Does a proection of the data.
%
%   input -----------------------------------------------------------------
%
%       o X : (N x D), a data set with N samples each being of dimension D.
%
%       o options : structure containing the options for the projection
%                   method in question.
%
%               options.method_name  = 'PCA' : this options should be given
%               options.nbDimensions =  3    : the number of dimensions to
%                                             keep when doing the
%                                             projection.
%
%   output ----------------------------------------------------------------
%
%      o projectedData : data after the projection, keeping only the number
%      of dimensions given in the options
%      o mapping : a structure containing the main parameter, results of
%      the algorithm (eigen values ...)
%
%
%

if ~isfield(options,'method_name')
    error('field of method_name of structure options was not defined!');
end

if ~isfield(options,'nbDimensions')
    warning('options.nbDimensions not specified!, setting default value, options.nbDimensions=2');
    options.nbDimensions = 2;
end

if ~isfield(options,'neighbors') && (strcmp(options.method_name,'Isomap') || strcmp(options.method_name,'LLE') || strcmp(options.method_name,'Laplacian'))
    warning('options.neighbors not specified!, setting default value, options.neighbors=7');
    options.neighbors = 7;
end

nbDimensions        = options.nbDimensions;

switch options.method_name

    case {'PCA', 'MDS'}
        [projectedData, mapping] = compute_mapping(X, options.method_name, nbDimensions);

    case 'KPCA'

        [kernel,kpar]                = sanitise_kernel_input(options);
        norm_K = true;

        if isfield(options,'norm_K')
            if (options.norm_K == false)
            display('Kernel (Gram) matrix will NOT be normalized!');
            norm_K = false;
            end
        end

        if length(kpar) == 1
            [projectedData, mapping] = compute_mapping(X, 'KPCA',nbDimensions,norm_K,kernel,kpar(1));
        else
            [projectedData, mapping] = compute_mapping(X, 'KPCA',nbDimensions,norm_K,kernel,kpar(1),kpar(2));
        end

    case {'Isomap'}
	neighbors = options.neighbors;
    [projectedData, mapping] = compute_mapping(X, options.method_name, nbDimensions, neighbors);

    case {'Laplacian'}
	neighbors = options.neighbors;
    sigma = options.sigma;
    [projectedData, mapping] = compute_mapping(X, options.method_name, nbDimensions, neighbors, sigma);

    otherwise
        error(['no such method [' options.method_name '] handled!']);

end

end

author: epfl-lasa


text:
function y = huber_pos( x, M, t )

%HUBER_POS   Monotonic huber penalty function.
%   For a vector X, HUBER_POS(X) computes the monotonic Huber-style function
%
%                      0     if 0>=X
%       HUBER_POS(X) = X^2   if 0<=X<=1
%                      2*X-1 if    X>=1
%
%   HUBER_POS(X,M) is the monotonic Huber-style penalty function of
%   halfwidth M, M.^2.*HUBER_POS(X./M). M must be real and positive.
%
%   HUBER_POS(X,M,T) computes the monotonic Huber-style penalty function
%   with halfwidth M and concomitant scale T:
%
%       HUBER_POS(X,M,T) = T.*HUBER_POS(X./T,M) if T > 0
%                          +Inf                 if T <= 0
%
%   See the help file for HUBER for information about this usage.
%
%   For matrices and N-D arrays, the penalty function is applied to each
%   element of X independently. M and T must be compatible with X in the same
%   sense as .*: one must be a scalar, or they must have identical size.
%
%   Disciplined convex programming information:
%       HUBER_POS is jointly convex in X and T. It is nondecreasing in X and
%       nonincreasing in T. Therefore, when used in CVX specifications, X
%       must be convex and T must be concave (or affine). Both must be real.

%
% Check arguments
%

error( nargchk( 1, 3, nargin ) );
if ~isreal( x ),
    error( 'First argument must be real.' );
end
if nargin < 2,
    M = 1;
elseif ~isreal( M ) || any( M( : ) <= 0 ),
    error( 'Second argument must be real and positive.' );
end
if nargin < 3,
    t = 1;
elseif ~isreal( t ),
    error( 'Third argument must be real.' );
end
sz = cvx_size_check( x, M, t );
if isempty( sz ),
    error( 'Sizes are incompatible.' );
end

%
% Compute result
%

y = max( x, 0 );
z = min( y, M );
y = t .* z .* ( 2 * y - z );
q = t <= 0;
if nnz( q ),
    if length( t ) == 1,
        y = Inf * ones( sy );
    else
        y( q ) = Inf;
    end
end

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


text:
function gT = rbfwhiteKernGradX(kern, t1, t2)

% RBFWHITEKERNGRADX Gradient of RBF-WHITE kernel with respect to a point t.
% FORMAT
% DESC computes the gradient of the RBF-WHITE kernel with respect to the
% input positions.
% ARG kern : kernel structure for which gradients are being computed.
% ARG t1 : locations against which gradients are being computed.
% RETURN gT : the returned gradients. The gradients are returned in
% a matrix which is numData x numInputs x numData. Where numData is
% the number of data points and numInputs is the number of input
% dimensions in t1 (currently always one).
%
% FORMAT
% DESC computes the gradident of the RBF-WHITE kernel with respect to the
% input positions where both the row positions and column positions are
% provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG t1 : row locations against which gradients are being computed.
% ARG t2 : column locations against which gradients are being computed.
% RETURN gT : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in t1, numData2 is the number of data
% points in t2 and numInputs is the number of input dimensions in t1
% (currently always one).
%
% SEEALSO rbfwhiteKernParamInit, kernGradX, rbfwhiteKernDiagGradX
%
% COPYRIGHT : David Luengo, 2009

% KERN

error('rbfwhiteKernGradX not yet implemented.')

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%  Figure 10.73      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% Fig. 10.73
% Data for RTP Demo 3-3-99
% Data provided by Dr. Gwen van der Linden
% Data is from System Identification Studies
InputFlux=[3.460064464376177e-1 1.177299050104922e-1 2.838023866104041e-2;
   3.880303397347619e-11 8.024902450324316e-2 1.807231516460469e-2;
   8.004191616976514e-9 2.721604310757543e-3 3.171348842079633e-2];
M_inv=diag([1.000040130716728 5.557442686788876 13.63821806414694]);
Radiation=[5.47621193859299e-2 -8.570695054070524e-3 -8.296135532988507e-4...
      -4.536181077856052e-2;
   -8.570695054070524e-3 8.570946319867835e-3 -1.621311365067015e-7...
      -8.913466080455817e-8;
   -8.296135532988507e-4 -1.621311365067015e-7 8.299854517643017e-4...
      -2.097673289443245e-7];
Conduction=[3.559939609150268e-7 -1.113667477845243e-7 -1.976161155515125e-7...
      -4.701109757899004e-8;
   -1.113667477845243e-7 1.160207476868843e-2 -2.502736022145532e-3...
      -9.099227379795117e-3;
   -1.976161155515125e-7 -2.502736022145532e-3 6.37364815665867e-3...
      -3.870714518397587e-3];
ScaleTemp=diag([0.01 0.01 0.01 0.01]);

clf;
%
sim('fig10_72')
%plot(tout,r,'-');
%hold on;
%plot(tout,y,'--');
%xlabel('Time (sec)');
%ylabel('Temperature (K)');
%hold off;
%pause;
ii=240:876;
plot(tout(ii),r(ii),'-',tout(ii),y(ii,2),'--');
legend('r','y');
grid on;
xlabel('Time (sec)');
ylabel('Temperature (K)');
title('Fig. 10.73(a) Temperature tracking response');
pause;
hold off;
ii=240:876;
plot(tout(ii),u(ii),'-');
xlabel('Time (sec)');
ylabel('Lamp voltage (V)');
grid on;
legend('u');
title('Fig. 10.73(b) Control effort');
%grid
nicegrid


text:
function value = sdot ( n, dx, incx, dy, incy )

%*****************************************************************************80
%
%% SDOT forms the dot product of two vectors.
%
%  Discussion:
%
%    This routine uses unrolled loops for increments equal to one.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vectors.
%
%    Input, real DX(*), the first vector.
%
%    Input, integer INCX, the increment between successive entries in DX.
%
%    Input, real DY(*), the second vector.
%
%    Input, integer INCY, the increment between successive entries in DY.
%
%    Output, real VALUE, the sum of the product of the
%    corresponding entries of DX and DY.
%
  value = 0.0;

  if ( n <= 0 )
    return
  end
%
%  Code for unequal increments or equal increments
%  not equal to 1.
%
  if ( incx ~= 1 | incy ~= 1 )

    if ( 0 <= incx )
      ix = 1;
    else
      ix = ( - n + 1 ) * incx + 1;
    end

    if ( 0 <= incy )
      iy = 1;
    else
      iy = ( - n + 1 ) * incy + 1;
    end

    for i = 1 : n
      value = value + dx(ix) * dy(iy);
      ix = ix + incx;
      iy = iy + incy;
    end
%
%  Code for both increments equal to 1.
%
  else

    m = mod ( n, 5 );

    for i = 1 : m
      value = value + dx(i) * dy(i);
    end

    for i = m+1 : 5 : n

      value = value + dx(i  ) * dy(i  ) ...
                    + dx(i+1) * dy(i+1) ...
                    + dx(i+2) * dy(i+2) ...
                    + dx(i+3) * dy(i+3) ...
                    + dx(i+4) * dy(i+4);
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ y, next ] = p11_equil ( neqn, next )

%*****************************************************************************80
%
%% P11_EQUIL returns equilibrium solutions of problem p11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%
  if ( next == 0 )
    next = 1;
    y(1:neqn,1) = 0.0;
  else
    next = 0;
    y = [];
  end

  return
end


text:
function [ fea, out ] = ex_convdiff3( varargin )
%EX_CONVDIFF3 1D Time dependent convection and diffusion equation example.
%
%   [ FEA, OUT ] = EX_CONVDIFF3( VARARGIN ) 1D time dependendt convection and diffusion equation on
%   a line with an infinately oscillating exact solution. Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       w           scalar {1.5*pi}        Simulation parameter
%       U           scalar {2}             Simulation parameter
%       a           scalar {1}             Convection velocity
%       nu          scalar {0.1}           Diffusion coefficient
%       hmax        scalar {1/25}          Max grid cell size
%       dt          scalar {0.02}          Time step size
%       ischeme     scalar {3}             Time stepping scheme
%       sfun        string {sflag2}        Shape function
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
  'w',        1.5*pi; ...
  'U',        2; ...
  'a',        1; ...
  'nu',       0.1; ...
  'hmax',     1/25; ...
  'dt'        0.02; ...
  'ischeme'   3; ...
  'sfun',     'sflag1'; ...
  'iplot',    1; ...
  'tol',      3e-2; ...
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

w  = opt.w;
U  = opt.U;
a  = opt.a;
nu = opt.nu;
l1 = ['(',num2str(a),'+sqrt(',num2str(a^2),'+',num2str(4*nu*w),'*i))/',num2str(2*nu)];
l2 = ['(',num2str(a),'-sqrt(',num2str(a^2),'+',num2str(4*nu*w),'*i))/',num2str(2*nu)];
refsol  = ['(exp(',l1,'*x)-exp(',l2,'*x))/(exp(',l1,')-exp(',l2,'))*',num2str(U),'*exp(i*',num2str(w),'*t)'];
refsol0 = ['(exp(',l1,'*x)-exp(',l2,'*x))/(exp(',l1,')-exp(',l2,'))*',num2str(U)];

% Grid generation.
fea.grid = linegrid( 1/opt.hmax, 0, 1 );

% Problem definition.
fea.sdim  = { 'x' };
fea = addphys( fea, @convectiondiffusion );
fea.phys.cd.sfun = { opt.sfun };
fea.phys.cd.eqn.coef{2,4} = { opt.nu };
fea.phys.cd.eqn.coef{3,4} = { opt.a  };
fea = parsephys(fea);

% Parse and solve problem.
fea = parseprob( fea );
fea.bdr.d{1} = 0;
fea.bdr.d{2} = [num2str(U),'*cos(',num2str(w),'*t)'];
fea.bdr.n    = cell(1,2);
x = fea.grid.p';
if( strcmp( opt.sfun,'sflag2' ) )
  x = [ x; (x(2:end)+x(1:end-1))/2 ];
end
init = real( eval( refsol0 ) );
[fea.sol.u,tlist] = solvetime( fea, 'fid', fid, 'init', init, 'ischeme', opt.ischeme, 'tstep', opt.dt, 'tmax', 1 );

% Postprocessing.
if( opt.iplot>0 )
  figure;
  if( opt.iplot>1 )
    i_sol_list = 1:numel(tlist);
  else
    i_sol_list = numel(tlist);
  end
  [~,ix] = sort( x );
  for i_sol=i_sol_list
    t = tlist(i_sol);
    clf
    postplot( fea, 'surfexpr', 'c', 'solnum', i_sol );
    hold on
    u_r = real( eval( refsol ) );
    plot( sort(x), u_r(ix), 'r--' );
    title( ['Solution at time ',num2str(t)])
    xlabel( 'x' )
    drawnow
  end
end

% Error checking.
for i_sol=1:numel(tlist)
  u_i = fea.sol.u(:,i_sol);
  t   = tlist(i_sol);
  u_r = real( eval( refsol ) );
  errnm(i_sol) = norm( u_i - u_r )/norm( u_r );
end
out.err  = errnm;
out.pass = all( errnm<opt.tol );

if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function pde = TorusTimeInitial3(mum,mup,sigm,sigp,epsm,epsp,omega,x0,y0,z0,r1,r2,a,b,intPt)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,...
    'f1',@f1,'f2',@f2,'f3',@f3,...
    'fm1',@fm1,'fm2',@fm2,'fm3',@fm3,...
    'fp1',@fp1,'fp2',@fp2,'fp3',@fp3,...
    'Mu',@Mu,'Mum',@Mum,'Mup',@Mup,'one',@one,...
    'Epslon',@Epslon,'Epslonm',@Epslonm,'Epslonp',@Epslonp,...
    'Sig',@Sig,'Sigm',@Sigm,'Sigp',@Sigp,...
    'E1',@E1,'E2',@E2,'E3',@E3,...
    'Et1',@Et1,'Et2',@Et2,'Et3',@Et3);

pde.mum = mum;
pde.mup = mup;
pde.sigm = sigm;
pde.sigp = sigp;
pde.epsm = epsm;
pde.epsp = epsp;
%% interface function
    function u = intf(x,y,z)
        u = (sqrt((x-x0).^2+(y-y0).^2)-r2).^2 + (z-z0).^2 - r1^2;
    end

%% exact solution

%% right hand side function
    function u = f1(x,y,z,t)
        u = fm1(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = fp1(x(id),y(id),z(id),t);
    end
    function u = f2(x,y,z,t)
        u = fm2(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = fp2(x(id),y(id),z(id),t);
    end
    function u = f3(x,y,z,t)
        u = fm3(x,y,z,t);
        id = intf(x,y,z) > 0;
        u(id) = fp3(x(id),y(id),z(id),t);
    end

    function u = fm1(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fm2(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fm3(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fp1(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fp2(x,y,z,t)
        u = zeros(size(x));
    end
    function u = fp3(x,y,z,t)
        u = zeros(size(x));
    end

%%%%% intial condition
    function u = E1(x,y,z,t)
        u = exp(-b*(a*(x-intPt)-omega*t).^2);
%         id = pulregion(x,y,z) > 0;
%         u(id) = zeros(size(x(id)));
    end
    function u = E2(x,y,z,t)
        u = -exp(-b*(a*(x-intPt)-omega*t).^2);
    end
    function u = E3(x,y,z,t)
        u = zeros(size(x));
    end
    function u = Et1(x,y,z,t)
        u = exp(-b*(a*(x-intPt)-omega*t).^2)*2*omega*b.*(a*(x-intPt)-omega*t);
    end
    function u = Et2(x,y,z,t)
        u = -exp(-b*(a*(x-intPt)-omega*t).^2)*2*omega*b.*(a*(x-intPt)-omega*t);
    end
    function u = Et3(x,y,z,t)
        u = zeros(size(x));
    end
%     function u = pulregion(x,y,z)
%        u = (x-intPt).^2/(0.15)^2 - 1;
%     end
%% Diffusion coefficient function
    function u = Mu(x,y,z)
        u = Mum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Mup(x(id),y(id),z(id));
    end
    function u = Mum(x,y,z)
        u = mum^(-1)*ones(size(x));
    end
    function u = Mup(x,y,z)
        u = mup^(-1)*ones(size(x));
    end

%% Mass coefficient function
    function u = Epslon(x,y,z)
        u = Epslonm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Epslonp(x(id),y(id),z(id));
    end
    function u = Epslonm(x,y,z)
        u = epsm*ones(size(x));
    end
    function u = Epslonp(x,y,z)
        u = epsp*ones(size(x));
    end

    function u = Sig(x,y,z)
        u = Sigm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Sigp(x(id),y(id),z(id));
    end
    function u = Sigm(x,y,z)
        u = sigm*ones(size(x));
    end
    function u = Sigp(x,y,z)
        u = sigp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end

repo: ifem


text:
%this script demos the usage of evaluation routines
% the result file 'demo.val.pred.txt' on validation data is evaluated
% against the ground truth

fprintf('CLASSIFICATION WITH LOCALIZATION TASK\n');

meta_file = '../data/meta_clsloc.mat';
%pred_file='demo.val.pred.loc.txt';
%pred_file='/data/vision/torralba/deeplearning/visualization/prediction_imagenetValSet_imagenetCNNaveSumDeep.txt'
pred_file='/data/vision/torralba/deeplearning/visualization/prediction_imagenetValSet_googlenet_imagenet.txt'

ground_truth_file='../data/ILSVRC2014_clsloc_validation_ground_truth.txt';
blacklist_file='../data/ILSVRC2014_clsloc_validation_blacklist.txt';
num_predictions_per_image=5;
optional_cache_file = 'cache_groundtruth.mat';

fprintf('pred_file: %s\n', pred_file);
fprintf('ground_truth_file: %s\n', ground_truth_file);
fprintf('blacklist_file: %s\n', blacklist_file);

if isempty(optional_cache_file)
    fprintf(['NOTE: you can specify a cache filename and the ground ' ...
             'truth data will be automatically cached to save loading time ' ...
             'in the future\n']);
end

% num_val_files = -1;
% while num_val_files ~= 50000
%     if num_val_files ~= -1
%         fprintf('That does not seem to be the correct directory. Please try again\n');
%     end
%     %ground_truth_dir=input('Please enter the path to the Validation bounding box annotations directory: ', 's');
     ground_truth_dir='/data/vision/torralba/deeplearning/imagenet_toolkit/val_bbox';
%     %fprintf('ground_truth_dir: %s\n', ground_truth_dir);
%     val_files = dir(sprintf('%s/*.xml',ground_truth_dir));
%     num_val_files = numel(val_files);
% end

error_cls = zeros(num_predictions_per_image,1);
error_loc = zeros(num_predictions_per_image,1);

for i=1:num_predictions_per_image
    [error_cls(i) error_loc(i)] = eval_clsloc(pred_file,ground_truth_file,ground_truth_dir,...
                                              meta_file,i, blacklist_file,optional_cache_file);
end

disp('# guesses vs clsloc error vs cls-only error');
disp([(1:num_predictions_per_image)',error_loc,error_cls]);

%% Result: alexFullConv network
% # guesses vs clsloc error vs cls-only error
%     1.0000    0.7461    0.5077
%     2.0000    0.6859    0.3841
%     3.0000    0.6569    0.3209
%     4.0000    0.6401    0.2835
%     5.0000    0.6270    0.2556

%% Result: googleNet_imagenet
% # guesses vs clsloc error vs cls-only error
%     1.0000    0.6691    0.3450
%     2.0000    0.6158    0.2294
%     3.0000    0.5937    0.1786
%     4.0000    0.5813    0.1503
%     5.0000    0.5734    0.1314

author: zhoubolei

repo: CAM

save_path: github-repos/MATLAB/zhoubolei-CAM


text:
function [ xpoly, ypoly, zpoly, avg_resid ] = sv2poly( xpos, ypos, zpos, ...
    times, xvel, yvel, zvel )
%SV2POLY Compute best fit polynomials from 3D state vectors
%
% This function determines a polynomial order that achieves good precision
% without being unreasonably high or overfitting.  One would generally only
% use this function on state vectors from an orbital system, where the
% positions/velocity would reasonably follow a smooth path.
%
% "Optimal" polynomial order depends on both the length of the collect
% around its orbit and the noise in the state vector data.
%
% Written by: Wade Schwartzkopf, NGA/R
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Check arguments
if ~isequal(size(xpos), size(ypos), size(zpos), size(times)) || ...
        (exist('zvel','var') && ...
        ~isequal(size(xvel), size(yvel), size(zvel), size(times)))
    error('SV2POLY:INVALID_INPUT_ARGS','All state vector components must be same size/shape.');
end
% Vectorize inputs to assure column vectors
times = times(:);
xpos = xpos(:);
ypos = ypos(:);
zpos = zpos(:);
if exist('zvel', 'var')
    xvel = xvel(:);
    yvel = yvel(:);
    zvel = zvel(:);
end

% Because of the scale of the values used in orbital state vectors,
% MATLAB's polyfit oftens detects them as badly conditioned, but results
% work out OK anyway for reasonable polynomail orders. We will assure in
% this function that polynomial order is reasonable.
old_state = warning('off','MATLAB:polyfit:RepeatedPointsOrRescale');

%% Method 1.  Naive approach.
% Hardcode polynomial order.  This is usually sufficient.  5th order
% generally works well for most spaceborne SAR collects of typical length
% (< 200 s), often putting evaluation of the polynomials well within 1 mm
% of error.
% polyorder = min(5, numel(times) - 1);

if exist('zvel', 'var')
    %% Method 2.  Cross check with velocity.
    % One could also find the order of polynomial that most accurately
    % describes this position, but use velocity (which is presumably
    % measured independently) as cross-validation so that the data is not
    % being overfit.
    for polyorder = 1:(numel(times)-1)
        xpoly = polyfit(times, xpos, polyorder);
        ypoly = polyfit(times, ypos, polyorder);
        zpoly = polyfit(times, zpos, polyorder);
        error_list(polyorder) = mean(sqrt(sum(([xvel, yvel, zvel] - ...
            [polyval(polyder(xpoly), times), ...
            polyval(polyder(ypoly), times), ...
            polyval(polyder(zpoly), times)]).^2, 2)));
    end
    % Increase order only as long as it results in a "significant" decrease
    % (half) in the error of the velocity.
    polyorder = find(error_list(1:(end-1))<(2*error_list(2:end)), 1, 'first');
    if isempty(polyorder), polyorder = numel(error_list); end  % All orders improve quality
else
    %% Method 3. "Significant" improvement.
    % One could also use the lowest order of polynomial that significantly
    % reduces the error of the position fit.
    for polyorder = 1:(numel(times)-2)  % polyorder of n-1 will always be exact match
        % mu term used for error computation in loop since we are
        % potentially computing high order polynomials, generally well
        % beyond what is required for a good fit, and these will likely be
        % badly conditioned. For our final fit, after the polynomial order
        % is determined, we will assume order and fit is reasonable-- even
        % if MATLAB would have thrown a warning.
        [xpoly, ~, xmu] = polyfit(times, xpos, polyorder);
        [ypoly, ~, ymu] = polyfit(times, ypos, polyorder);
        [zpoly, ~, zmu] = polyfit(times, zpos, polyorder);
        error_list(polyorder) = mean(sqrt(sum(([xpos, ypos, zpos] - ...
            [polyval(xpoly, times, [], xmu), ...
            polyval(ypoly, times, [], ymu), ...
            polyval(zpoly, times, [], zmu)]).^2, 2)));
    end
    % Increasing polynomial order by one must result in error being cut in
    % half in order to be considered "significant".
    polyorder = find(error_list(1:(end-1))<(2*error_list(2:end)), 1, 'first');
    if isempty(polyorder), polyorder = numel(error_list); end  % All orders improve quality
end

% Once optimal polynomial order is determined, do actual polynomial fit
xpoly  = polyfit(times, xpos, polyorder);
ypoly  = polyfit(times, ypos, polyorder);
zpoly  = polyfit(times, zpos, polyorder);

warning(old_state);

% Compute final residuals if requested
if nargout>3
    avg_resid = mean(sqrt(sum(([xpos, ypos, zpos] - ...
        [polyval(xpoly, times), ...
        polyval(ypoly, times), ...
        polyval(zpoly, times)]).^2, 2)));
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function run_experiment_baseline_conv5(imdb_video)
%% Experiment entry point

	opts.gpus = 1;

	if nargin < 1
	    imdb_video = [];
	end

    opts.join.method = 'xcorr';
    opts.branch.conf.num_out = [96 256 384 384 32];
    opts.branch.conf.num_in = [ 3 48 256 192 192];
    opts.branch.conf.conv_stride = [ 2 1 1 1 1];
    opts.branch.conf.pool_stride = [ 2 1];

    opts.exemplarSize = 127;
   	opts.train.numEpochs = 100;

	experiment(imdb_video, opts);

end

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function [ t, wts ] = sgqf ( nt, aj, bj, zemu )

%*****************************************************************************80
%
%% SGQF computes knots and weights of a Gauss Quadrature formula.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with simple knots from the Jacobi matrix and the zero-th
%    moment of the weight function, using the Golub-Welsch technique.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real AJ(NT), the diagonal of the Jacobi matrix.
%
%    Input, real BJ(NT), the subdiagonal of the Jacobi
%    matrix, in entries 1 through NT-1.  On output, BJ has been overwritten.
%
%    Input, real ZEMU, the zero-th moment of the weight function.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Exit if the zero-th moment is not positive.
%
  if ( zemu <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SGQF - Fatal error!\n' );
    fprintf ( 1, '  ZEMU <= 0.\n' );
    error ( 'SGQF - Fatal error!' );
  end
%
%  Set up vectors for IMTQLX.
%
  wts = zeros ( nt, 1 );

  wts(1) = sqrt ( zemu );
  wts(2:nt) = 0.0;
%
%  Diagonalize the Jacobi matrix.
%
  [ t, wts ] = imtqlx ( nt, aj, bj, wts );

  wts(1:nt) = wts(1:nt).^2;

  return
end


text:
function [lambda_vec, error_train, error_val] = ...
    validationCurve(X, y, Xval, yval)
%VALIDATIONCURVE Generate the train and validation errors needed to
%plot a validation curve that we can use to select lambda
%   [lambda_vec, error_train, error_val] = ...
%       VALIDATIONCURVE(X, y, Xval, yval) returns the train
%       and validation errors (in error_train, error_val)
%       for different values of lambda. You are given the training set (X,
%       y) and validation set (Xval, yval).
%

% Selected values of lambda (you should not change this)
lambda_vec = [0 0.001 0.003 0.01 0.03 0.1 0.3 1 3 10]';

% You need to return these variables correctly.
error_train = zeros(length(lambda_vec), 1);
error_val = zeros(length(lambda_vec), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return training errors in
%               error_train and the validation errors in error_val. The
%               vector lambda_vec contains the different lambda parameters
%               to use for each calculation of the errors, i.e,
%               error_train(i), and error_val(i) should give
%               you the errors obtained after training with
%               lambda = lambda_vec(i)
%
% Note: You can loop over lambda_vec with the following:
%
%       for i = 1:length(lambda_vec)
%           lambda = lambda_vec(i);
%           % Compute train / val errors when training linear
%           % regression with regularization parameter lambda
%           % You should store the result in error_train(i)
%           % and error_val(i)
%           ....
%
%       end
%
%

for i = 1:length(lambda_vec),
  lambda = lambda_vec(i);
  theta = trainLinearReg(X, y, lambda);
  error_train(i) = linearRegCostFunction(X, y, theta, 0);
  error_val(i) = linearRegCostFunction(Xval, yval, theta, 0);
endfor

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

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
function Matrices = removeExplorationConstraints(Matrices);
candidates = find((~any(Matrices.G,2)) & (sum(Matrices.E | Matrices.E,2) == 1));
if ~isempty(candidates)
    Matrices.bndA = -Matrices.E(candidates,:);
    Matrices.bndb = Matrices.W(candidates,:);
    Matrices.G(candidates,:) = [];
    Matrices.E(candidates,:) = [];
    Matrices.W(candidates,:) = [];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [sig_e, bay,model] = bay_errorbar(model,Xt, type, nb, bay)
% Compute the error bars for a one dimensional regression problem
%
% >> sig_e = bay_errorbar({X,Y,'function',gam,sig2}, Xt)
% >> sig_e = bay_errorbar(model, Xt)
%
% The computation takes into account the estimated noise variance
% and the uncertainty of the model parameters, estimated by
% Bayesian inference. sig_e is the estimated standard deviation of
% the error bars of the points Xt. A plot is obtained by replacing
% Xt by the string 'figure'.
%
%
% Full syntax
%
%     1. Using the functional interface:
%
% >> sig_e = bay_errorbar({X,Y,'function',gam,sig2,kernel,preprocess}, Xt)
% >> sig_e = bay_errorbar({X,Y,'function',gam,sig2,kernel,preprocess}, Xt, type)
% >> sig_e = bay_errorbar({X,Y,'function',gam,sig2,kernel,preprocess}, Xt, type, nb)
% >> sig_e = bay_errorbar({X,Y,'function',gam,sig2,kernel,preprocess}, 'figure')
% >> sig_e = bay_errorbar({X,Y,'function',gam,sig2,kernel,preprocess}, 'figure', type)
% >> sig_e = bay_errorbar({X,Y,'function',gam,sig2,kernel,preprocess}, 'figure', type, nb)
%
%       Outputs
%         sig_e         : Nt x 1 vector with the [$ \sigma^2$] errorbands of the test data
%       Inputs
%         X             : N x d matrix with the inputs of the training data
%         Y             : N x 1 vector with the inputs of the training data
%         type          : 'function estimation' ('f')
%         gam           : Regularization parameter
%         sig2          : Kernel parameter
%         kernel(*)     : Kernel type (by default 'RBF_kernel')
%         preprocess(*) : 'preprocess'(*) or 'original'
%         Xt            : Nt x d matrix with the inputs of the test data
%         type(*)       : 'svd'(*), 'eig', 'eigs' or 'eign'
%         nb(*)         : Number of eigenvalues/eigenvectors used in the eigenvalue decomposition approximation
%
%     2. Using the object oriented interface:
%
% >> [sig_e, bay, model] = bay_errorbar(model, Xt)
% >> [sig_e, bay, model] = bay_errorbar(model, Xt,       type)
% >> [sig_e, bay, model] = bay_errorbar(model, Xt,       type, nb)
% >> [sig_e, bay, model] = bay_errorbar(model, 'figure')
% >> [sig_e, bay, model] = bay_errorbar(model, 'figure', type)
% >> [sig_e, bay, model] = bay_errorbar(model, 'figure', type, nb)
%
%       Outputs
%         sig_e     : Nt x 1 vector with the [$ \sigma^2$] errorbands of the test data
%         model(*)  : Object oriented representation of the LS-SVM model
%         bay(*)    : Object oriented representation of the results of the Bayesian inference
%       Inputs
%         model     : Object oriented representation of the LS-SVM model
%         Xt        : Nt x d matrix with the inputs of the test data
%         type(*)   : 'svd'(*), 'eig', 'eigs' or 'eign'
%         nb(*)     : Number of eigenvalues/eigenvectors used in the eigenvalue decomposition approximation
%
% See also:
%   bay_lssvm, bay_optimize, bay_modoutClass, plotlssvm

% Copyright (c) 2002,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.ac.be/sista/lssvmlab

if iscell(model), model = initlssvm(model{:}); end

if model.type(1)~='f',
  error(['confidence bounds only for function estimation. For' ...
	 ' classification, use ''bay_modoutClass(...)'' instead;']);
end

eval('type;','type=''svd'';');
eval('nb;','nb=model.nb_data;');
if ~(strcmpi(type,'svd') | strcmpi(type,'eig') | strcmpi(type,'eigs') | strcmpi(type,'eign')),
  error('Eigenvalue decomposition via ''svd'', ''eig'', ''eigs'' or ''eign''...');
end

if strcmpi(type,'eign')
  warning('The resulting errorbars are most probably not very usefull...');
end

if ~isstr(Xt),

  eval('[sig_e, bay] = bay_confb(model,Xt,type,nb,bay);',...
       '[sig_e, bay] = bay_confb(model,Xt,type,nb);');

else

  grid = 50;
  [X,Y] = postlssvm(model,model.xtrain,model.ytrain);
  eval('[sig_e, bay] = bay_confb(model,X,type,nb,bay);',...
       '[sig_e, bay] = bay_confb(model,X,type,nb);');

  % plot the curve including confidence bound
  sige = sqrt(sig_e);
  Yt = simlssvm(model,X);

  figure;
  hold on;
  title(['LS-SVM_{\gamma=' num2str(model.gam(1)) ', \sigma^2=' num2str(model.kernel_pars(1)) ...
	 '}^{' model.kernel_type(1:3) '} and its 95% (2\sigma) error bands']);

  if model.x_dim==1,
    xlabel('X');
    ylabel('Y');
    [~,si] = sort(X);
    plot(X(si),Yt(si),'k'); hold on;
    plot(X(si),Yt(si)+2.*sige(si),'-.r');
    plot(X(si),Yt(si)-2.*sige(si),':r');
    plot(X(si),Y(si),'k*'); hold off;
  else
    xlabel('time');
    ylabel('Y');
    plot(Yt,'k'); hold on;
    plot(Yt+2.*sige,'-.r');
    plot(Yt-2.*sige,':r');
    plot(Y,'k*'); hold off;
  end

end

function [sig_e, bay] = bay_confb(model,X,type,nb,bay)
% see formula's thesis TvG blz 126

nD = size(X,1);
%tol = .0001;

%
% calculate the eigenvalues
%
eval('bay;','[c1,c2,c3,bay] = bay_lssvm(model,1,type,nb);');
omega = kernel_matrix(model.xtrain(model.selector,1:model.x_dim), ...
		      model.kernel_type, model.kernel_pars);
oo = ones(1,model.nb_data)*omega;

% kernel values of  X
theta = kernel_matrix(model.xtrain(model.selector, 1:model.x_dim), ...
		      model.kernel_type, model.kernel_pars, X);
for i=1:nD,
  kxx(i,1) = feval(model.kernel_type, X(i,:),X(i,:), model.kernel_pars);
end

Zc = eye(model.nb_data) - ones(model.nb_data)./model.nb_data;

Hd = (Zc*bay.Rscores);
Hd = Hd*diag(1./bay.mu - (bay.mu+ bay.zeta*bay.eigvals).^-1)*Hd';

% forall x
for i=1:nD,
  term1(i,1) = bay.zeta^-1 + kxx(i)/bay.mu - theta(:,i)'*Hd*theta(:,i);
  term2(i,1) = 2/model.nb_data*sum(theta(:,i)'*Hd*omega) - 2/bay.mu/model.nb_data* sum(theta(:,i));
end

% once
term3 = 1/(bay.zeta*model.nb_data) ...
	+ 1/(bay.mu*model.nb_data^2)* sum(oo)  ...
	-1/(model.nb_data^2)* oo*Hd*oo';

sig_e = term1+term2+term3;

author: peterhcharlton

repo: RRest


text:
function [ type, typeCost ] = assignVanishingType( lines, vp, tol, area )
%ASSIGNVANISHINGTYPE Summary of this function goes here
%   Detailed explanation goes here
if nargin<=3
    area = 10;
end

numLine = size(lines,1);
numVP = size(vp, 1);
typeCost = zeros(numLine, numVP);
% perpendicular
for vid = 1:numVP
    cosint = dot( lines(:,1:3), repmat( vp(vid,:), [numLine 1]), 2);
    typeCost(:,vid) = asin(abs(cosint));
end
% infinity
for vid = 1:numVP
    valid = true(numLine,1);
    for i = 1:numLine
        us = lines(i,5);
        ue = lines(i,6);
        u = [us;ue]*2*pi-pi;
        v = computeUVN(lines(i,1:3), u, lines(i,4));
        xyz = uv2xyzN([u v], lines(i,4));
        x = linspace(xyz(1,1),xyz(2,1),100);
        y = linspace(xyz(1,2),xyz(2,2),100);
        z = linspace(xyz(1,3),xyz(2,3),100);
        xyz = [x' y' z'];
        xyz = xyz ./ repmat(sqrt(sum(xyz.^2,2)),[1 3]);
        ang = acos( abs(dot(xyz, repmat(vp(vid,:), [100 1]), 2)));
        valid(i) = ~any(ang<area*pi/180);
    end
    typeCost(~valid,vid) = 100;
end

[I, type] = min(typeCost,[],2);
type(I>tol) = numVP+1;

end

author: zouchuhang

repo: LayoutNet


text:
function b = rhs(disc, f)
%RHS   Discretize the right-hand side of a linear system for VALSDISCRETIZATION.
%   B = RHS(DISC, F) returns a discrete version, B,  of the function (or
%   chebmatrix) F, as defined by the discretization DISC.
%
% See also MATRIX, REDUCE.

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

% Developers note: This method works as follows.
%   The original function(s) are discretized, then reduced in dimension the same
%   way as the operator. The constraints are prepended to the top of the vector.

xOut = equationPoints(disc);
for k = 1:numel(f.blocks)
    if ( ~isnumeric(f.blocks{k}) )
        f.blocks{k} = feval(f.blocks{k}, xOut);
    end
end
b = cell2mat(f.blocks);

% Developer note:
%   The continuity conditions go above the constraints. See getConstraints().

% Prepend the values of the constraints and continuity conditions.
L = disc.source;

if ( ~isempty(L.constraint) )
    b = [ L.constraint.values ; b ];
end
if ( ~isempty(L.continuity) )
    b = [ L.continuity.values ; b ];
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% This script provides a tutorial application of the Gaussian-Copula Mutual
% Information (GCMI) estimator with a 2 category discrete stimulus (face
% vs noise images) event-related EEG data set.

% This script uses cell mode - cells are delimited by %% lines and can
% be run with:
% ctrl-enter (windows, linux) or cmd-enter (mac)
% or "Run Section" button from the toolbar
% or right click -> Evaluate currect section

% You need the gcmi/matlab directory on your matlab path.
% If you are running from the gcmi project tree you can add this with the
% following command:
% addpath('../matlab')

% Questions / comments : robince@gmail.com

%% Download and load data

% this will attempt to download the eeg data from the internet (~210MB)
% alternatively you can download the file eeg_face_vs_noise.mat manually
% and place it in the same directory as this script
fname = 'eeg_face_vs_noise.mat';
data_url = ['https://www.robince.net/data/gcmi/data/' fname];

% checks for data in current working directory
if ~exist(fname)
    disp(sprintf('Downloading %s ...', data_url))
    websave(fname,data_url);
    disp('Done.')
end

load(fname);
% csddat : [channels x time x trials]
% permute to trials first axis for looping
csddat = permute(csddat, [3 1 2]);
[Ntrl, Nch, Nt] = size(csddat);

% stim : stimulus class on each trial (0, 1)

%% Calculate GCMI at a specific sensor and time point

% this is the strongest effect
chi = 41;
ti = 168;

% for a t-test we contrast the data from the two classes
faceeeg = csddat(stim==0,chi,ti);
noiseeeg = csddat(stim==1,chi,ti);
[h,p,ci,stats] = ttest2(faceeeg, noiseeeg, 'VarType', 'unequal');
t = stats.tstat

% for MI we use the stimulus labels for each trial
% this works in the same way if there are more than 2 classes
% Reminder: stim must take values 0 or 1
I = gcmi_model_cd(csddat(:,chi,ti), stim, 2)

%% Calculate GCMI across all sensors and time point

% following the commonly used mass-univeriate approach we repeat
% the calculation above for each sensor and time point
% we could call gcmi_cd as above inside the loop, but here we first
% normalise the data separately so we can reuse it for the permutation
% testing

% Gaussian-copula normalisation, works along first axis, applied to each
% other dimension independently
ceeg = copnorm(csddat);
Ieeg = zeros(Nch,Nt);
for ti=1:Nt
    for chi=1:Nch
        % as the data has been copula-normalised we can use the
        % Gaussian parametric estimator (whatever the distribution was
        % originally)
        Ieeg(chi,ti) = mi_model_gd(ceeg(:,chi,ti), stim, 2, true, true);
    end
end

figure
imagesc(time,[],Ieeg)
xlabel('Time')
ylabel('Channels')

%% Permutation test for signifiance with the method of maximum statistics
Nperm = 100;

% to reduce computation time, we consider a single time point here - the
% same method should normally be applied over all sensors and time points
ti = 168;
I = Ieeg(:,ti);
Iperm = zeros(Nch,Nperm);
for pi=1:Nperm
    idx = randperm(Ntrl);
    % randomly permute stimulus labels
    pstim = stim(idx);
    % repeat mass-univariate MI calculation
    for chi=1:Nch
        Iperm(chi,pi) = mi_model_gd(ceeg(:,chi,ti), pstim, 2, true, true);
    end
end

% method of maximum statistics
% maximum values across permutations
Imax = max(Iperm,[],1);
thresh = prctile(Imax, 99);
Isig = I>thresh;

% if you have eeglab we can plot the topologies
if exist('topoplot')
    figure
    subplot(1,2,1)
    topoplot(I, chanlocs, 'maplimits', [0 max(I)]);
    title(sprintf('MI, t = %d ms', time(ti)))
    colorbar

    subplot(1,2,2)
    topoplot(Isig, chanlocs, 'maplimits', [0 1])
    title(sprintf('MI, p<0.01, t = %d ms', time(ti)));
    colorbar

    colormap parula
end


text:
function [cst X A B C] = test05(m, n)
% function [cst X A B C] = test05(m, n)
% All intputs are optional.
%
% Test case for the Euclidean space (classical nonlinear optimization).
% We generate three matrices A, B and C and try to solve the Sylvester
% equation AX + XB = C by minimizing the residual ||AX+XB-C||^2.
%

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    if ~exist('m', 'var') || isempty(m)
        m = 3;
    end
    if ~exist('n', 'var') || isempty(n)
        n = 10;
    end

    % Create the problem structure
    problem.M = euclideanfactory(m, n);

    % Define the problem data
    A = randn(m);
    B = randn(n);
    C = randn(m, n);

    % Define the problem cost function
    problem.cost = @(X) .5*norm(A*X+X*B-C, 'fro')^2;
    problem.grad = @(X) A'*(A*X+X*B-C) + (A*X+X*B-C)*B';
    problem.hess = @(X, H) A'*(A*H+H*B) + (A*H+H*B)*B';

    % Check differnetials consistency.
    checkgradient(problem);
    checkhessian(problem);

    % Solve
    [X cst] = trustregions(problem);

end

author: NicolasBoumal

repo: manopt


text:
function out = true( in )
%TRUE for initialization of nested cell arrays
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

out = cellfun(@(x) true(x), in.data{2,1}, 'UniformOutput', false);
out = reconFlatCellMatrix(out, in.data{1,1});
out = TRAFO(out,in.meta);

end

author: thomaskuestner


text:
% LFMatlabPathSetup - convenience function to add the light field toolbox to matlab's path
%
% It may be convenient to add this to your startup.m file.
%
% User guide: LFToolbox.pdf

% Copyright (c) 2013-2020 Donald G. Dansereau

function LFMatlabPathSetup

if (~isdeployed)
	% Find the path to this script, and use it as the base path
	LFToolboxPath = fileparts(mfilename('fullpath'));

	fprintf('Adding paths for LF Toolbox ');
	addpath( fullfile(LFToolboxPath) );
	addpath( fullfile(LFToolboxPath, 'SupportFunctions') );
	addpath( fullfile(LFToolboxPath, 'SupportFunctions', 'CameraCal') );

	fprintf('%s, done.\n', LFToolboxVersion);
end

repo: LFToolbox


text:
function []=save_optimal_network(meth_Net,BrainNet,label,result_dir,varargin)
% Generate one of the output, i.e., the group averaged network for each group generated by the optimized parameters.
switch meth_Net
    case {'PC','tHOFC','aHOFC'}
        opt_BrainNet=BrainNet;
    case {'SR','SLR','SGR','GSR','WSR','WSGR','SSGSR','dHOFC'}
        Acc_para=varargin{1};
        [~,B]=max(Acc_para);
        opt_BrainNet=BrainNet{B(1)};
end

label_negative=find(label==-1);
label_positive=find(label==1);
BrainNet_negative_mean=mean(opt_BrainNet(:,:,label_negative),3);
BrainNet_positive_mean=mean(opt_BrainNet(:,:,label_positive),3);

%figure;
figure('visible','off');
subplot(1,2,1);
imagesc(BrainNet_negative_mean);
colormap jet
colorbar
axis square
xlabel('ROI');
ylabel('ROI');
title('label = -1');

subplot(1,2,2);
imagesc(BrainNet_positive_mean);
colormap jet
colorbar
axis square
xlabel('ROI');
ylabel('ROI');
title('label = 1');
print(gcf,'-r1000','-dtiff',char(strcat(result_dir,'/Mean_optimal_network.tiff')));
save (char(strcat(result_dir,'/Mean_optimal_negativeLabel_network.mat')),'BrainNet_negative_mean');
save (char(strcat(result_dir,'/Mean_optimal_positiveLabel_network.mat')),'BrainNet_positive_mean');

author: zzstefan

repo: BrainNetClass


text:
classdef testAmplificatorTensorForRectangleInclusion < ...
         testAmplificatorTensorForInclusions

    properties (Access = protected)
       % testName = 'AmplificatorTensorForRectangleInclusion';
    end

    methods (Access = protected)

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function x = TT_to_TTeMPS( tt )
    %TT_to_TTeMPS Convert from TT Toolbox format.
    %   A = TT_to_TTeMPS( tt ) takes the tt_tensor object tt created using the
    %   TT Toolbox 2.x from Oseledets et al. and converts it into a TTeMPS tensor.
    %   Toolbox needs to be installed, of course.
    %
    %   See also TTeMPS_to_TT.

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

    cores = {};
    ps = tt.ps;
    for i = 1:tt.d
        cores{i} = reshape( tt.core( ps(i):ps(i+1)-1 ), ...
                                [tt.r(i), tt.n(i), tt.r(i+1)] );
    end

    x = TTeMPS( cores );

end


text:
function model = multimodelCreate(inputDim, outputDim, varargin)

% MULTIMODELCREATE Create a MULTIMODEL model.
% The MULTIMODEL is a way of performing multi-task learning by sharing
% model parameters across a range of models. The default (simple)
% assumption is that the data is conditionally independent given the
% parameters, i.e. the log likelihood is the sum of the log likelihood of
% the models.
%
% SEEALSO : modelCreate
%
% FORMAT
% DESC creates a multi-task learning wrapper
% model structure given an options structure.
% ARG inputDim : the input dimension of the model.
% ARG outputDim : the output dimension of the model.
% ARG options : an options structure that determines the form of the model.
% RETURN model : the model structure with the default parameters placed in.
%
% SEEALSO : multimodelOptions, multimodelParamInit, modelCreate
%
% COPYRIGHT : Neil D. Lawrence, 2007
%
% MODIFICATIONS : Mauricio A. Alvarez, 2009

% MLTOOLS

options = varargin{end};
model.numModels = options.numModels;
model.type = 'multimodel';
model.compType = options.type;
model.inputDim = inputDim;
model.outputDim = outputDim;

if numel(outputDim) == 1
    % Indices of parameters to be trained separately for each model.
    model.separateIndices = options.separate;
    model.numSep = length(model.separateIndices);
    for i = 1:model.numModels
        varargput = cell(1, length(varargin)-1);
        for j = 1:length(varargput)
            varargput{j} = varargin{j}{i};
        end
        % MAURICIO : temporarily changed to allow different options for each model
        if ~iscell(options.compOptions)
            model.comp{i} = modelCreate(model.compType, inputDim, outputDim, ...
                varargput{:}, options.compOptions);
        else
            model.comp{i} = modelCreate(model.compType, inputDim, outputDim, ...
                varargput{:}, options.compOptions{i});
        end
    end
    if isfield(model.comp{i}, 'numParams');
        model.numParams = model.comp{1}.numParams;
    else
        model.numParams = length(modelExtractParam(model.comp{1}));
    end
    model.sharedIndices = 1:model.numParams;
    model.sharedIndices(model.separateIndices) = [];
    model.numParams = model.numParams + (model.numModels-1)*model.numSep;
else
    model.separateIndices = options.separate;
    model.numSep = length(model.separateIndices);
    for i = 1:model.numModels
        varargput = cell(1, length(varargin)-1);
        for j = 1:length(varargput)
            varargput{j} = varargin{j}{i};
        end
        fprintf('Creating model number: %d\n', i)
        model.comp{i} = modelCreate(model.compType, inputDim, outputDim(i), ...
            varargput{:}, options.compOptions{i});

    end
    if isfield(model.comp{1}, 'nParams');
        model.numParams = 0;
        for i =1:model.numModels,
            model.numParams = model.numParams + model.comp{i}.nParams;
        end
    else
        model.numParams = length(modelExtractParam(model.comp{1}));
    end
    model.numParams = model.numParams + (model.numModels-1)*model.numSep;
end

if isfield(options, 'optimiser') && ~isempty(options.optimiser)
    model.optimiser = options.optimiser;
end


text:
function test_ft_spike_rate()

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_spike_rate

clear
spikesPerTrial = 1;
nTrials = 100;
shapePar = 2;
scalePar = 3;
data = [];
time       = linspace(0,1,1000);
data.trial(1:nTrials) = {zeros(1,length(time))};
data.time(1:nTrials) = {time};
for iTrial = 1:nTrials
  for iUnit = 1:4
      data.trial{iTrial}(iUnit,:) = rand(1,1000)<iTrial/100;
      if iUnit==4
        data.trial{iTrial}(iUnit,:) = rand(1,1000);
      end
  end
end
data.fsample = 1000;
data.hdr = [];
data.cfg.trl = [];
data.label{1} = 'chan1';
data.label{end+1} = 'chan2';
data.label{end+1} = 'chan3';
data.label{end+1} = 'chan4';

% show that the psth works also with the poisson format
cfg.spikechannel = 1:3;
spike = ft_checkdata(data,'datatype', 'spike', 'feedback', 'yes');

%%
cfgRate = [];
cfgRate.outputunit = 'spikecount';
cfgRate.keeptrials = 'yes';
cfgRate.vartriallen = 'no';
cfgRate.latency = [-10 20];
RateCnt1 = ft_spike_rate(cfgRate,spike);

%%
cfgRate = [];
cfgRate.outputunit = 'spikecount';
cfgRate.keeptrials = 'no';
cfgRate.trials = 1:2:30;
RateCnt1 = ft_spike_rate(cfgRate,spike);

cfgRate.trials = 32:2:60;
RateCnt2 = ft_spike_rate(cfgRate,spike);
cfgRate.trials = 61:2:99;
RateCnt3 = ft_spike_rate(cfgRate,spike);

cfg.stimuli = [0 pi/8 pi/4];
cfg.method = 'orientation';
stat = ft_spike_rate_orituning(cfg,RateCnt1,RateCnt2,RateCnt3);
%%

cfgRate = [];
cfgRate.outputunit = 'spikecount';
cfgRate.keeptrials = 'yes';
cfgRate.trials = 1:5:30;
RateCnt1 = ft_spike_rate(cfgRate,spike);
disp('expect about 10 times the number of spikes as the trial number')
[RateCnt1.cfg.trials(:) RateCnt1.trial]

author: fieldtrip

repo: fieldtrip


text:
% SYNTAX:
% [yAvg, yAvgStd, tHRF, nTrials] = hmrS_SessAvg_Nirs(yAvgRuns, yAvgStdRuns, ySum2Runs, tHRFRuns, mlActRuns, nTrialsRuns)
%
% UI NAME:
% Run_Average
%
% DESCRIPTION:
% Calculate avearge HRF of all runs for one subject.
%
% INPUTS:
% yAvgRuns:
% yAvgStdRuns:
% tHRFRuns:
% mlActRuns:
% nTrialsRuns:
% trange: defines the range for the block average
% thresh: Threshold for excluding channels if it's data deviates too much
%         from mean
%
% OUTPUTS:
% yavg: the averaged results
% yAvgStd: the standard deviation across trials
% tHRF: the time vector
% nTrials: the number of trials averaged for each condition across all runs
%
% USAGE OPTIONS:
% Run_Average_on_Concentration_Data:  [dcAvg, dcAvgStd, tHRF, nTrials]    = hmrS_SessAvg_Nirs(dcAvgRuns, dcAvgStdRuns, dcSum2Runs, tHRFRuns, mlActRuns, nTrialsRuns)
% Run_Average_on_Delta_OD_Data:       [dodAvg, dodAvgStd, tHRF, nTrials]  = hmrS_SessAvg_Nirs(dodAvgRuns, dodAvgStdRuns, dodSum2Runs, tHRFRuns, mlActRuns, nTrialsRuns)
%

function [yAvg, yAvgStd, tHRF, nTrials] = hmrS_SessAvg_Nirs(yAvgRuns, yAvgStdRuns, ySum2Runs, tHRFRuns, mlActRuns, nTrialsRuns)

yAvg = [];
yAvgStd = [];
tHRF = [];
nTrials_tot = [];
grp1=[];

for iRun = 1:length(yAvgRuns)

    yAvg      = yAvgRuns{iRun};
    yAvgStd   = yAvgStdRuns{iRun};
    ySum2     = ySum2Runs{iRun};
    tHRF      = tHRFRuns{iRun};
    nTrials   = nTrialsRuns{iRun};
    mlAct     = mlActRuns{iRun};

    if isempty(yAvg)
        break;
    end

    nCond = size(nTrials,2);

    if ndims(yAvg) == (3-(nCond<2))

        if isempty(mlAct)
            mlAct = 1:size(yAvg,2);
        end

        % grab tHRF to make common for group average
        if iRun==1
            grp1 = zeros(size(yAvg,1), size(yAvg,2), nCond);
            grp1Sum2 = zeros(size(yAvg,1), size(yAvg,2), nCond);
            nTrials_tot = zeros(size(yAvg,2), nCond);
        end

        for iC = 1:nCond
            nT = nTrials(iC);
            if nT>0
                if iRun==1
                    grp1(:,mlAct,iC) = yAvg(:,mlAct,iC) * nT;
                    grp1Sum2(:,mlAct,iC) = ySum2(:,mlAct,iC);
                    nTrials_tot(mlAct,iC) = nT;
                else
                    for iCh=1:length(mlAct) %size(yAvg,2)
                        grp1(:,mlAct(iCh),iC) = grp1(:,mlAct(iCh),iC) + interp1(tHRF,yAvg(:,mlAct(iCh),iC),tHRF(:)) * nT;
                        grp1Sum2(:,mlAct(iCh),iC) = grp1Sum2(:,mlAct(iCh),iC) + interp1(tHRF,ySum2(:,mlAct(iCh),iC),tHRF(:));
                        nTrials_tot(mlAct(iCh),iC) = nTrials_tot(mlAct(iCh),iC) + nT;
                    end
                end
            end
        end
        yAvg    = [];
        yAvgStd = [];
        if ~isempty(grp1)
            for iC = 1:size(grp1,3)
                for iCh = 1:size(grp1,2)
                    yAvg(:,iCh,iC) = grp1(:,iCh,iC) / nTrials_tot(iCh,iC);

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

    elseif ndims(yAvg) == (4-(nCond<2))

        if isempty(mlAct)
            mlAct = 1:size(yAvg,3);
        end

        % grab tHRF to make common for group average
        if iRun==1
            grp1 = zeros(size(yAvg,1), size(yAvg,2), size(yAvg,3), nCond);
            grp1Sum2 = zeros(size(yAvg,1), size(yAvg,2), size(yAvg,3), nCond);
            nTrials_tot = zeros(size(yAvg,3), nCond);
        end

        for iC = 1:1:nCond
            nT = nTrials(iC);
            if nT>0
                if iRun==1
                    grp1(:,:,mlAct,iC) = yAvg(:,:,mlAct,iC) * nT;
                    grp1Sum2(:,:,mlAct,iC) = ySum2(:,:,mlAct,iC);
                    nTrials_tot(mlAct,iC) = nT;
                else
                    for iCh=1:length(mlAct) %size(yAvg,3)
                        for iHb=1:size(yAvg,2)
                            grp1(:,iHb,mlAct(iCh),iC) = grp1(:,iHb,mlAct(iCh),iC) + interp1(tHRF,yAvg(:,iHb,mlAct(iCh),iC),tHRF(:)) * nT;
                            grp1Sum2(:,iHb,mlAct(iCh),iC) = grp1Sum2(:,iHb,mlAct(iCh),iC) + interp1(tHRF,ySum2(:,iHb,mlAct(iCh),iC),tHRF(:));
                        end
                        nTrials_tot(mlAct(iCh),iC) = nTrials_tot(mlAct(iCh),iC) + nT;
                    end
                end
            end
        end
        yAvg    = [];
        yAvgStd = [];
        if ~isempty(grp1)
            for iC = 1:size(grp1,4)
                for iCh = 1:size(grp1,3)
                    yAvg(:,:,iCh,iC) = grp1(:,:,iCh,iC) / nTrials_tot(iCh,iC);

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

author: BUNPC


text:
function test04 ( m, n )

%*****************************************************************************80
%
%% TEST04 estimates integrals in M dimensions, using N points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of evaluation points.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  P00_Q returns the integral of F over [0,1]^m.\n' );
  fprintf ( 1, '  Here, we use spatial dimension M = %d\n', m );
  fprintf ( 1, '  The number of sample points is N = %d\n', n );

  seed = 123456789;
  [ x, seed ] = r8mat_uniform_01 ( m, n, seed );

  prob_num = p00_prob_num ( );

  for prob = 1 : prob_num

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Problem %d\n', prob );

    [ c, seed ] = p00_c ( prob, m, seed );
    r8vec_print ( m, c, '  C parameters:' );

    [ w, seed ] = p00_w ( prob, m, seed );
    r8vec_print ( m, w, '  W parameters:' );

    fprintf ( 1, '\n' );
    fprintf ( 1, '      Exact Integral     Q\n' );
    fprintf ( 1, '\n' );

    q1 = p00_q ( prob, m, c, w );

    f = p00_f ( prob, m, c, w, n, x );
    q2 = sum ( f(1:n) ) / n;

    fprintf ( 1, '  %14g  %14g\n', q1, q2 );

  end

  return
end


text:
%ZSSD Sum of squared differences
%
% M = ZSSD(I1, I2) is the zero-mean sum of squared differences between the
% two equally sized image patches I1 and I2.  The result M is a scalar that
% indicates image similarity, a value of 0 indicates identical pixel patterns
% and is increasingly positive as image dissimilarity increases.
%
% Notes::
% - The ZSSD similarity measure is invariant to changes in image brightness
%   offset.
%
% See also SDD, SAD, NCC, ISIMILARITY.

% 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 m = zssd(w1, w2)

	w1 = w1 - mean(w1(:));
	w2 = w2 - mean(w2(:));

	m = (w1-w2).^2;
    m = sum(m(:));

author: petercorke

repo: machinevision-toolbox-matlab


text:
%MAKESHAPEIMAGES  One-line description here, please.
%
%   output = makeShapeImages(input)
%
%   Example
%   makeShapeImages
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2021-07-09,    using Matlab 9.8.0.1323502 (R2020a)
% Copyright 2021 INRAE.

% read data
img = imread('cameraman.tif');

% flipH
imgFlip = imFlip(img, 1);
imwrite(imgFlip, 'cameraman_flipH.png', 'png');

% rotate 90
imgRot90 = imRotate90(img);
imwrite(imgRot90, 'cameraman_rot90.png', 'png');

% crop border
imgCropBorder = imCropBorder(img, [20 40 30 50]);
imwrite(imgCropBorder, 'cameraman_cropBorder.png', 'png');

% add border
imgAddBorder = imAddBorder(img, [20 40 30 50]);
imwrite(imgAddBorder, 'cameraman_addBorder.png', 'png');

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [adfreq, n, ts, fn, ad] = plx_ad(filename, channel)
% plx_ad(filename, channel): read all a/d data for the specified channel from a .plx or .pl2 file
%              returns raw a/d values
%
% [adfreq, n, ts, fn, ad] = plx_ad(filename, channel)
% [adfreq, n, ts, fn, ad] = plx_ad(filename, 0)
% [adfreq, n, ts, fn, ad] = plx_ad(filename, 'FP01')
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   channel - 0-based channel number or channel name
%
%           a/d data come in fragments. Each fragment has a timestamp
%           and a number of a/d data points. The timestamp corresponds to
%           the time of recording of the first a/d value in this fragment.
%           All the data values stored in the vector ad.
%
% OUTPUT:
%   adfreq - digitization frequency for this channel
%   n - total number of data points
%   ts - array of fragment timestamps (one timestamp per fragment, in seconds)
%   fn - number of data points in each fragment
%   ad - array of raw a/d values

adfreq = 0;
n = 0;
ts = -1;
fn = -1;
ad = -1;

if nargin ~= 2
    error 'expected 2 input arguments';
end

[ filename, isPl2 ] = internalPL2ResolveFilenamePlx( filename );
if isPl2 == 1
    % return data from .pl2 file
    pl2Channel = channel;
    % PL2 uses 1-based channel numbers
    if ischar(channel) == 0
        pl2Channel = channel+1;
    end
    pl2ad = PL2Ad(filename, pl2Channel);
    if numel(pl2ad.Values) > 0
        adfreq = pl2ad.ADFreq;
        n = numel(pl2ad.Values);
        ts = pl2ad.FragTs;
        fn = pl2ad.FragCounts;
        ad = pl2ad.Values;
        % convert millivolts to raw a/d values
        pl2 = PL2GetFileIndex(filename);
        channelNumber = internalPL2ResolveChannel(pl2.AnalogChannels, pl2Channel);
        ad = round(ad/(pl2.AnalogChannels{channelNumber}.CoeffToConvertToUnits * 1000));
    end
    return;
end

% return plx data
channelNumber = plx_ad_resolve_channel(filename, channel);
if channelNumber == -1
    fprintf('\n plx_ad: no header for the specified A/D channel.');
    fprintf('\n    use plx_ad_info(filename) to print the list of a/d channels\n');
    return
end
[adfreq, n, ts, fn, ad] = mexPlex(2, filename, channelNumber);
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Function to difference each column in a matrix
%
% Description
%  Function to difference each column in a matrix
%
%  This version is a slightly later one than that described in the
%  above published 2013 CSL paper [2]. The algorithm here rather than using binary
%  decision trees using artificial neural networks and combines the
%  features used in the CSL paper with those proposed in Ishi et al.
%  (2008). This updated version has been submitted to CSL for a special
%  issue on glottal source processing on April 14th 2013. It will have
%  reference [1].
%
% Inputs
%  mat  : [samples] [NxM] Feature matrix
%
% Outputs
%  delta_mat : [samples] [NxM] Differenced feature matrix
%
% Example
%  Please see the HOWTO_glottalsource.m example file.
%
% References
%  [1] Drugman, T., Kane, J., Gobl, C., `Automatic Analysis of Creaky
%       Excitation Patterns', Submitted to Computer Speech and
%       Language.
%  [2] Kane, J., Drugman, T., Gobl, C., (2013) `Improved automatic
%       detection of creak', Computer Speech and Language 27(4), pp.
%       1028-1047.
%  [3] Drugman, T., Kane, J., Gobl, C., (2012) `Resonator-based creaky
%       voice detection', Interspeech 2012, Portland, Oregon, USA.
%
% Copyright (c) 2013 University of Mons, FNRS & 2013 Trinity College Dublin
%
% License
%  This code is a part of the GLOAT toolbox with the following
%  licence:
%  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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Authors
%  Thomas Drugman <thomas.drugman@umons.ac.be> & John Kane <kanejo@tcd.ie>

function delta_mat = get_delta_mat(mat)

delta_mat=zeros(size(mat));
N=size(mat,2);

for n=1:N
    delta_mat(2:end,n)=diff(mat(:,n));
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function MRSIStruct = setAffineMatrix(MRSIStruct, newAffineMatrix)
    MRSIStruct.affineMatrix = newAffineMatrix;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
addpath(genpath('./solver/SDPNAL+v1.0/'), path);
run('./solvers/SDPNAL+v1.0/startup.m');
addpath(genpath('./utils'));
addpath('./utils/kdtree/toolbox/');

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
% envproj() - plot envelopes of projections of selected ICA component
%             projections against envelope of the original data
%
% Usage:    >> [envdata] = envproj(data,weights,compnums);
%           >> [envdata] = envproj(data,weights,compnums, ...
%                                  title,limits,chanlist,compnames,colors);
% Inputs:
%   data      = runica() input data (chans,frames) <- best one epoch only!
%   weights   = unmixing weight matrix (runica() weights*sphere)
%   compnums  = list of component numbers to project and plot
%
% Optional inputs:
%   title     = 'fairly short plot title' (in quotes) (0 -> none)
%   limits    = [xmin xmax ymin ymax] (x's in msec)
%                       (0 or both y's 0 -> [min max])
%   chanlist  = list of data channels to use for max/min (0 -> all)
%   compnames = file of component labels (exactly 5 chars per line,
%               trailing '.' = space) (default|0 -> compnums)
%   colors    = file of color codes, 3 chars per line  ('.' = space)
%                       (0 -> default color order (black/white first))
%
% Output:
%   envdata   = [2,(frames*length(compnums)+1)] envelopes of data, comps.
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 1/1997
%
% See also: envtopo()

% Copyright (C) 01-23-97 Scott Makeig, SCCN/INC/UCSD, scott@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

% 03-19-97 use datamean instead of frames/baseframes, diag(cov()) for var() -sm
% 04-03-97 changed name to envproj() -sm
% 05-20-97 used sum-of-squares for var instead of diag() to allow long data
%          read MAXPLOTDATACHANS from icadefs.m -sm
% 06-07-97 changed order of args to conform to runica -sm
% 06-18-97 read MAXENVPLOTCHANS instead of MAXPLOTDATACHANS -sm
% 07-23-97 dropped datamean from args; mean is distributed among components -sm
% 10-11-97 range of max/min for default ylimits over chanlist only! -sm
% 11-05-97 disallowed white lines unless axis color is white -sm & ch
% 11-13-97 rm'd errorcode variable to improve limits handling -sm
% 11-18-97 added legend to plot; varied line types -sm
% 11-30-97 added f=figure below to keep previous plot from being altered -sm
% 12-01-97 added thicker lines for data envelope -sm
% 12-08-97 debugged thicker lines for data envelope -sm
% 12-10-97 changed f=figure below to f=gcf to avoid matlab bug (gcf confusion) -sm
% 12-10-97 implemented compnames arg -sm
% 12-13-97 took out compnames=compnames'; bug -sm
% 02-06-98 added 'show-largest' feature to legend when numcomps > 7 -sm
% 02-09-98 fixed numcomps bug -sm
% 02-12-98 test for outofbounds compnums early -sm
% 04-28-98 made to work for rectangular weights -sm
% 06-05-98 made NEG_UP optional -sm
% 02-22-99 added FILL mode default when numcomps==1 -sm
% 03-14-99 made envelope finding code more efficient -sm
% 12-19-00 updated icaproj() args -sm
% 01-25-02 reformated help & license, added links -ad

function [envdata] = envproj(data,weights,compnums,titl,limits,chanlist,compnames,colors)

FILL = 1; % 1;  % 1 = fill in component envelope unless ncomps>1
DATA_LINEWIDTH = 1.5; % 2
COMP_LINEWIDTH = 1.5; % 1
MAXENVPLOTCHANS = 256;

if nargin < 3,
    help envproj
    fprintf('envproj(): must have at least three arguments.\n');
    return
end

icadefs        % load ENVCOLORS & MAXENVPLOTCHANS default variables
MAX_LEG = 7;   % maximum number of component line types to label in legend
NEG_UP = 0;    % 1 = plot negative up

%
%%%%%%%%%%%% Substitute for omitted arguments %%%%%%%%%%%%%%%%%%%%%%%%
%
if nargin < 8,
    colors = 'envproj.col';
elseif colors(1)==0,
    colors = 'envproj.col';
end

[chans,frames] = size(data);

if nargin < 7
   compnames = 0;
end

if nargin < 6
    chanlist = 0;
end
if nargin < 5,
    limits = 0;
end
if nargin < 4,
    titl = 0;
end
%
%%%%%%%%%%%%%%%%%%%%%%%% Test data size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[wr,wc]       = size(weights);
%
%%%%%%%%%%%%%%%% Substitute for zero arguments %%%%%%%%%%%%%%%%%%%%%%%%%
%
if chanlist == 0,
    chanlist = [1:chans];
end;
if compnums == 0,
    compnums = [1:wr];
end;
if size(compnums,1)>1,        % handle column of compnums !
    compnums = compnums';
end;
numcomps = length(compnums);
if numcomps > MAXENVPLOTCHANS,
    fprintf(...
'envproj(): cannot plot more than %d channels of data at once.\n',...
                  MAXENVPLOTCHANS);
    return
end;

if max(compnums)>wr
 fprintf('\nenvproj(): Component index %d out of bounds (1:%d).\n',...
                                             max(compnums),wr);
 return
end

if min(compnums)<1,
 fprintf('\nenvproj(): Component index %d out of bounds (1:%d).\n',...
                                             min(compnums),wr);
 return
end

if FILL>0 & numcomps == 1
   FILL = 1;
else
   FILL = 0;
end
%
%%%%%%%%%%%%%%%%%%%% Read and adjust limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  if limits==0,
    xmin=0;xmax=1;
    ymin=min(min(data(chanlist,:)));
    ymax=max(max(data(chanlist,:)));
  else
    if length(limits)~=4,
        fprintf( ...
 'envproj(): limits should be 0 or an array [xmin xmax ymin ymax].\n');
        return
      end;
        if limits(1,1) == 0 & limits(1,2) ==0,
            xmin=0;
            xmax=0;
        else
            xmin = limits(1,1);
            xmax = limits(1,2);
        end;
         if limits(1,3) == 0 & limits(1,4) ==0,
            ymin=0;
            ymax=0;
        else
            ymin = limits(1,3);
            ymax = limits(1,4);
        end;
  end;

  if xmax == 0 & xmin == 0,
    x = (0:1:frames-1);
    xmin = 0;
    xmax = frames-1;
  else
    dx = (xmax-xmin)/(frames-1);
    x=xmin*ones(1,frames)+dx*(0:frames-1); % construct x-values
  end;

  if ymax == 0 & ymin == 0,
    ymax=max(max(data(chanlist,:)));
    ymin=min(min(data(chanlist,:)));
  end

%
%%%%%%%%%%%%%%%%%%%% Read the color names %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    if ~isstr(colors)
        fprintf('envproj(): color file name must be a string.\n');
        return
    end
    cid = fopen(colors,'r');
    if cid <3,
        fprintf('envproj(): cannot open file %s.\n',colors);
        return
    else
        colors = fscanf(cid,'%s',[3 MAXENVPLOTCHANS]);
        colors = colors';
        [r c] = size(colors);
        for i=1:r
            for j=1:c
                if colors(i,j)=='.',
                    colors(i,j)=' ';
                end;
            end;
        end;
    end;
[rr cc] = size(colors);

%
%%%%%%%%%%%% Compute projected data for single components %%%%%%%%%%%%%
%
projdata = [data(chanlist,:)];
envdata = [min(projdata); max(projdata)]; % begin with envelope of data
fprintf('envproj(): Projecting component(s) ');
maxvars = zeros(1,length(compnums));
n=1;
for c=compnums,            % for each component
  fprintf('%d ',c);
  % [icaproj] = icaproj(data,weights,compindex); % new arg order 12/00
  proj = icaproj(data,weights,c); % let offsets be
                                         % distributed among comps.
  [val,i] = max(sum(proj.*proj)); % find max variance
  maxvars(n) = val;
  proj = proj(chanlist,:);
  projdata = [projdata proj];
  envtmp = [ min(proj) ; max(proj)];
  envdata = [envdata envtmp];
                 % append envelope of projected data sets onto envdata
  % Note: size(envdata) = [length(chanlist)  frames*(numcomps+1)]
  n = n+1;
end;
fprintf('\n');
%
%%%%%%%%%%%%%%%%%%%%%%%% Make the plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if ~isstr(titl)
  if titl==0,
    titl = ' ';
  else
    fprintf('envproj(): titl unrecognized.\n');
    return
  end
end

set(gcf,'Color',BACKCOLOR); % set the background color
set(gca,'Color','none'); % set the axis color to the background color
if NEG_UP
  set(gca,'Ydir','reverse');
end
n = 1;
H = zeros(1,numcomps+1);

for c=1:numcomps+1;
  %
  % Plot envelopes
  %
  if n <= 6
    if n == 1              % thick lines, store H for legend()
      H(n)=plot(x,envdata(1,(c-1)*frames+1:c*frames)',...
                                     colors(c),'Linewidth',DATA_LINEWIDTH);
      hold on;
           plot(x,envdata(2,(c-1)*frames+1:c*frames)',...
                                     colors(c),'Linewidth',DATA_LINEWIDTH);
    else
       if FILL
            fillx = [x x(frames:-1:1)];
            filly = [envdata(1,(c-1)*frames+1:c*frames) ...
                            envdata(2,c*frames:-1:(c-1)*frames+1)];
            i = find(isnan(filly)==1);
            if ~isempty(i)
                fprintf('Replacing %d NaNs in envelope with 0s.\n',length(i)/2);
                filly(i)=0;
            end
            H(n) = fill(fillx,filly,colors(c));
       else
            H(n)= plot(x,envdata(1,(c-1)*frames+1:c*frames)',...
                                     colors(c),'Linewidth',COMP_LINEWIDTH);
             hold on;
            plot(x,envdata(2,(c-1)*frames+1:c*frames)',...
                                     colors(c),'Linewidth',COMP_LINEWIDTH);
       end
    end
  else
    H(n)= plot(x,envdata(1,(c-1)*frames+1:c*frames)',...
                                     colors(c),'LineStyle',':','LineWidth',DATA_LINEWIDTH);
    hold on;
          plot(x,envdata(2,(c-1)*frames+1:c*frames)',...
                                     colors(c),'LineStyle',':','LineWidth',DATA_LINEWIDTH);
  end
  n = n+1;
end
set(gca,'Color','none'); % set the axis color to the background color

l=xlabel('Time (ms)');      % xaxis label
set(l,'FontSize',14);

l=ylabel('Potential (uV)'); % yaxis label
set(l,'FontSize',14);

if xmax > 0  & xmin < 0
  plot([0 0],[ymin ymax],'k','linewidth',1.5); % plot vertical line at time 0
end

if ymax ~= 0 |ymin~= 0,
    axis([min(x),max(x),ymin,ymax]);
else
    axis([min(x),max(x),min(envdata(1,:)'),max(envdata(2,:)')]);
end
legtext = [ 'Data ' ];
if compnames(1) ~= 0
    pid = fopen(compnames,'r');
    if pid < 3
       fprintf('envproj(): file %s not found.\n',compnames);
       return
    end
    compnames = fscanf(pid,'%s',[5 numcomps]);
    compnames = [ ['Data.']' compnames];
    compnames = compnames';
    if size(compnames,1) ~= numcomps+1
      fprintf('envproj(): no. of compnames must equal no. of compnums.\n');
      return
    end
    for c=1:5 % remove padding with .'s
      for r = 1:size(compnames,1)
         if compnames(r,c) == '.'
           compnames(r,c) = ' ';
         end
      end
    end
    legtext = compnames;
else
  legtext = ['Data'];
  for c = compnums
      legtext = strvcat(legtext,int2str(c));
  end
end

if numcomps<MAX_LEG+1
  if FILL <= 0  % no legend in FILL mode
       l=legend(H,legtext,0);      % plot legend
  end
else
  fprintf('Finding largest components for legend...\n');
  [maxvars,i] = sort(maxvars);
  maxvars = maxvars(length(maxvars):-1:1); % sort large-to-small
  i= i(length(i):-1:1);
  H = [H(1) H(i+1)];
  H = H(1:MAX_LEG+1);
  legtext = [legtext(1,:);legtext(i+1,:)];
  legtext = legtext(1:MAX_LEG+1,:);
  legend(H,legtext,0); % MAX_LEG largest-peaked comps
end

set(gca,'FontSize',12);
fprintf('If desired, use mouse to move legend or >> delete(legend).\n');

%
%%%%%%%%%%%%% Compute percentage of variance accounted for %%%%%%%%%%%%
%
sumdata = zeros(length(chanlist),frames);
for e=2:numcomps+1,              % sum the component activities
    sumdata = sumdata + projdata(:,(e-1)*frames+1:e*frames);
end

sigssqr = sum(data(chanlist,:).*data(chanlist,:))/(length(chanlist)-1);
dif    = data(chanlist,:)-sumdata;
difssqr = sum(dif.*dif)/(length(chanlist)-1);
pvaf   = round(100.0*(1.0-difssqr/sigssqr)); % % variance accounted for
rtitl = ['(' int2str(pvaf) '%)'];
%titl = [titl '   ' rtitl];

t=title(titl);              % plot title
set(t,'FontSize',14);

author: PatternRecognition

repo: OpenBMI


text:
%Thanks for Dr. Ma Zheng's X-ray diffraction Datas.
%You can test ma.mat, ma2.mat yourself.
%Enjoy!!!
load ma3.mat;
plot(A,B,'-');
grid on;
hold on;
t=fpeak(A,B,30,[23,90,700,inf]);
plot(t(:,1),t(:,2),'o');
title('\fontsize{24}Perfect!');
hold off;

author: Sable

repo: mcbench-benchmarks


text:
function [sig,params] = prod_sigmoid_logdist(w,data1,data2,ndx1,ndx2,ddim)
%
% Algorithm: sig = distribute(ndx1,sigmoid(
%                   log(
%                   sum(bsxfun(@minus,M*data_1,c).^2,1)
%                  )))
%                   *
%                  distribute(ndx2,sigmoid(
%                   log(
%                   sum(bsxfun(@minus,M*data_2,c).^2,1)
%                  )))
%
%
% Inputs:
%  w: is vec([M,c]), where M is ddim-by-D and c is ddim-by-1
%     Use w=[] to let output sld be an MV2DF function handle.
%
%  data_1: D-by-T1 matrix
%  data_2: D-by-T2 matrix
%  ndx1,ndx2: indices of size 1 by T to distribute T1 and T2 segs over T
%             trials
%
%  ddim: the first dimension of the M matrix
%
% Outputs:
%   sig: function handle (if w=[]), or numeric
%   params.get_w0(ssat): returns w0 for optimization initialization,
%                        0<ssat<1 is required average sigmoid output.
%   params.tail: is tail of parameter w, which is not consumed by this
%                function.

if nargin==0
    test_this();
    return;
end

datadim = size(data1,1);
assert(datadim==size(data2,1),'data1 and data2 must have same number of rows');
assert(length(ndx1)==length(ndx2));
assert(max(ndx1)<=size(data1,2));
assert(max(ndx2)<=size(data2,2));

wsz = ddim*(datadim+1);
[whead,wtail] = splitvec_fh(wsz,w);
params.get_w0 = @(ssat) init_w0(ssat);
params.tail = wtail;

sqd1 = square_distance_mv2df([],data1,ddim); %Don't put whead in here,
sqd2 = square_distance_mv2df([],data2,ddim); %or here. Will cause whead to be called twice.

sig1 = one_over_one_plus_w_mv2df(sqd1);
sig2 = one_over_one_plus_w_mv2df(sqd2);

% distribute over trials
distrib = duplicator_fh(ndx1,size(data1,2));
sig1 = distrib(sig1);
distrib = duplicator_fh(ndx2,size(data2,2));
sig2 = distrib(sig2);

sigh = dottimes_of_functions([],sig1,sig2);
sig = sigh(whead);

    function w0 = init_w0(ssat)
    W0 = randn(ddim,datadim+1); %subspace projector
    W0(:,end) = 0; % centroid from which distances are computed

    d0 = (1-ssat)/ssat;
    s = sigh(W0(:));
    d = (1-s)./s;
    W0 = sqrt(d0/median(d))*W0;
    w0 = W0(:);
    end

end

function test_this()

K = 5;
N1 = 10;
N2 = 2*N1;
data1 = randn(K,N1);
data2 = randn(K,N2);
ndx1 = [1:N1,1:N1];
ndx2 = 1:N2;

ddim = 3;
ssat = 0.01;
[sys,params] = prod_sigmoid_logdist([],data1,data2,ndx1,ndx2,ddim);

w0 = params.get_w0(ssat);
test_MV2DF(sys,w0);

sig = prod_sigmoid_logdist(w0,data1,data2,ndx1,ndx2,ddim),

end

author: nesl


text:
function  ACFTA  = readACFTA(fid)
%
% <<<<<<<<<< CLASSIFICATION: UNCLASSIFIED >>>>>>>>>>
%
% READACFTA Create structure for metadata fields.
%   ACFTA = READACFTA(FID) returns a structure of fields and their
%   associated values.

% Written by Matt Donath, NGA, matthew.b.donath@nga.ic.gov
% References:
%  > MIL-STD-2500A, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.0
%  > MIL-STD-2500C, NATIONAL IMAGERY TRANSMISSION FORMAT VERSION 2.1
%  > STDI-0001,     NATIONAL SUPPORT DATA EXTENSIONS (SDE)(VERSION 1.3/CN2)
%                   FOR THE NATIONAL IMAGERY TRANSMISSION FORMAT (NITF)
%  > STDI-0002,     THE COMPENDIUM OF CONTROLLED EXTENSIONS (CE) FOR THE
%                   NATIONAL IMAGERY TRANSMISSION FORMAT (NITF) VERSION 2.1
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%%  Create the metadata structure for ACFTA.

ACFTA.CETAG = fread(fid,6,'uint8=>char')';
ACFTA.CEL = str2double(fread(fid,5,'uint8=>char')');

ACFTA.AC_MSN_ID = fread(fid,10,'uint8=>char')';
ACFTA.SCTYPE = fread(fid,1,'uint8=>char')';
ACFTA.SCNUM = fread(fid,4,'uint8=>char')';
ACFTA.SENSOR_ID = fread(fid,3,'uint8=>char')';
ACFTA.PATCH_TOT = str2double(fread(fid,4,'uint8=>char')');
ACFTA.MTI_TOT = str2double(fread(fid,3,'uint8=>char')');
ACFTA.PDATE = fread(fid,7,'uint8=>char')';
ACFTA.IMHOSTNO = str2double(fread(fid,3,'uint8=>char')');
ACFTA.IMREQID = str2double(fread(fid,5,'uint8=>char')');
ACFTA.SCENE_SOURCE = fread(fid,1,'uint8=>char')';
ACFTA.MPLAN = fread(fid,2,'uint8=>char')';
ACFTA.ENTLOC = fread(fid,21,'uint8=>char')';
ACFTA.ENTELV = str2double(fread(fid,6,'uint8=>char')');
ACFTA.EXITLOC = fread(fid,21,'uint8=>char')';
ACFTA.EXITELV = str2double(fread(fid,6,'uint8=>char')');
ACFTA.TMAP = str2double(fread(fid,7,'uint8=>char')');
ACFTA.RCS = str2double(fread(fid,3,'uint8=>char')');
ACFTA.ROW_SPACING = str2double(fread(fid,7,'uint8=>char')');
ACFTA.COL_SPACING = str2double(fread(fid,7,'uint8=>char')');
ACFTA.SENSERIAL = fread(fid,4,'uint8=>char')';
ACFTA.ABSWVER = fread(fid,7,'uint8=>char')';

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
% by Tolga Birdal
% A sample application and a function for solving the maximum inscribed
% circle problem.
% Unlike my other submission "Maximum Inscribed Circle using Distance
% Transform", this algorithm is subpixel accurate. It operates only on the
% polygon and not the image points. Therefore, if the polygon is given in
% sub-pixels, the result will be accurate.
% I use an O(n log(n)) algorithm as follows:
% Construct the Voronoi Diagram of the polygon.
% For Voronoi nodes which are inside the polygon:
%       Find the node with the maximum distance to edges in P. This node is
%       the centre of the maximum inscribed circle.
%
% For more details on the problem itself please checkout my previous
% submission as mentioned above.
%
% To speed things up:
% Replace "inpolygon" function by Bruno Lunog's faster implementation:
% "2D polygon interior detection" :
% http://www.mathworks.com/matlabcentral/fileexchange/27840-2d-polygon-inte
% rior-detection
% Copyright (c) 2011, Tolga Birdal <http://www.tbirdal.me>

function []=test_max_circle2()
close all;
figure,

% Read the boundary segmented and the original image. Original image is
% only needed for visualization
Iorg=imread('hand.jpg');
I=imread('hand_contour.png');
imshow(Iorg);

% obtain the boundary
[y,x]=find(I>0);
contour = bwtraceboundary(logical(I), [y(1),x(1)], 'E', 8);

% get the circle
[cx,cy,r]=find_inner_circle(contour(:,2),contour(:,1));

% plot
theta = [linspace(0,2*pi) 0];
hold on, plot(x, y,'g.', 'MarkerSize',1);
hold on, plot(cos(theta)*r+cx,sin(theta)*r+cy,'r', 'LineWidth', 1);

end

% given a polygon [x,y] find the inner circle [cx,cy,r]
function [cx,cy,r]=find_inner_circle(x,y)

% make a voronoi diagram
[vx,vy]=voronoi(x,y);

% find voronoi nodes inside the polygon [x,y]
Vx=vx(:);
Vy=vy(:);
% Here, you could also use a faster version of inpolygon
IN=inpolygon(Vx,Vy, x,y);
ind=find(IN==1);
Vx=Vx(ind);
Vy=Vy(ind);

% maximize the distance of each voronoi node to the closest node on the
% polygon.
minDist=0;
minDistInd=-1;
for i=1:length(Vx)
    dx=(Vx(i)-x);
    dy=(Vy(i)-y);
    r=min(dx.*dx+dy.*dy);
    if (r>minDist)
        minDist=r;
        minDistInd=i;
    end
end

% take the center and radius
cx=Vx(minDistInd);
cy=Vy(minDistInd);
r=sqrt(minDist);

end


text:
function k = findrows(A, b)
%FINDROWS Find indices of a given row within a matrix.
%
%   FINDROWS(A, B) returns a column vector with the indices of the rows
%   in the matrix A that are identical to the row vector B.  If no rows
%   in A are identical to B, an empty vector is returned.
%
%   The methods uses a for-loop, but it uses less memory and is in many
%   cases a lot faster than the vectorized methods
%
%      find( all( A == repmat(b, size(A, 1), 1), 2 ) )
%      find( all( A == b(ones(size(A, 1), 1),:), 2 ) )
%
%   See also FIND, FINDCOLS.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:51:15 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

if isempty(A)
    k = [];
    return
end

k = find( A(:,1) == b(1));
top = size(A, 2);
for j = 2:top
    k = k(A(k,j) == b(j));
    if isempty(k)
        return
    end
end

author: yalmip

repo: YALMIP


text:
% This script plots the Gaussian sensing matrices generated by
% ex_generate_sample_sensing_matrices.m script.
% Please run the dependencies before running this script.

%% Initialization code
clear all; close all; clc;
% Let us load sample Rademacher sensing matrices
load('bin/gaussian_sensing_matrices.mat');

[numMatrices, numSets] = size(gaussianSensingMatrices);

% We pickup the first set
phis = gaussianSensingMatrices(:, 1);
spx.graphics.figure.full_screen;
for i=1:numMatrices
    % We pickup the matrix
    phi  = phis{i};
    % We get its size
    [M, N]  = size(phi);
    % We draw the matrix
    subplot(111);
    imagesc(phi) ;
    colormap(gray);
    axis image;
    xlabel(sprintf('\\Phi [%dx%d]', M, N));
    pause;
end

author: indigits

repo: sparse-plex


text:
function polyhedra(varargin)
%POLYHEDRA Index of classical polyhedral meshes
%
%   Polyhedra are specific meshes, with additional assumptions:
%   * the set of faces is assumed to enclose a single 3D domain
%   * each face has a neighbor face for each edge
%   * some functions also assume that normals of all faces point outwards
%
%   Most polyhedron creation functions follow the patterns:
%   * [V, F] = createXXX();     % returns vertex and face arrays
%   * [V, E, F] = createXXX();  % returns also edge array
%   * M = createXXX();          % return a data structure with 'vertices',
%                               % 'edges' and 'faces' fields.
%
%   Example
%   % create a soccer ball mesh and display it
%   [n, f] = createSoccerBall;
%   drawMesh(n, f, 'faceColor', 'g', 'linewidth', 2);
%   axis equal;
%
%   See also
%   meshes3d
%   createCube, createCubeOctahedron, createIcosahedron, createOctahedron
%   createRhombododecahedron, createTetrahedron, createTetrakaidecahedron
%   createDodecahedron, createSoccerBall, createMengerSponge
%   steinerPolytope, minConvexHull
%   polyhedronNormalAngle, polyhedronMeanBreadth
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2008-10-13,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% HISTORY

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


%   getBox:  Given the current state, returns a number from 1 to 162
%             designating the region of the state space encompassing the current state.
%             Returns a value of -1 if a failure state is encountered.

function box = getBox5(theta,thetaDot,x,xDot)
theta = rad2deg(theta);
thetaDot = rad2deg(thetaDot);
if (x < -2.4 || x > 2.4  || theta < -12 || theta > 12)
    box = 379;
else

if(theta>=-12 && theta<=-10)
    thetaBucket = 1;
elseif(theta>-10 && theta<=-8)
    thetaBucket = 2;
elseif(theta>-8 && theta<=-6)
    thetaBucket = 3;
elseif(theta>-6 && theta<=-4)
    thetaBucket = 4;
elseif(theta>-4 && theta<=-2)
    thetaBucket = 5;
elseif(theta>-2 && theta<=-1)
    thetaBucket = 6;
elseif(theta>-1 && theta<=0)
    thetaBucket = 7;
elseif(theta>0 && theta<=1)
    thetaBucket = 8;
elseif(theta>1 && theta<=2)
    thetaBucket = 9;
elseif(theta>2 && theta<=4)
    thetaBucket = 10;
elseif(theta>4 && theta<=6)
    thetaBucket = 11;
elseif(theta>6 && theta<=8)
    thetaBucket = 12;
elseif(theta>8 && theta<=10)
    thetaBucket = 13;
elseif(theta>10 && theta<=12)
    thetaBucket = 14;
end

if (x<-0.8&&x>=-2.4)
	xBucket = 1;
elseif (x<=0.8&&x>=-0.8)
	xBucket = 2;
elseif (x<=2.4&&x>0.8)
	xBucket = 3;
end

if (xDot<-0.5)
	xDotBucket = 1;
elseif (xDot>=-0.5&&xDot<=0.5)
	xDotBucket = 2;
else
	xDotBucket = 3;
end

if (thetaDot<-50)
	thetaDotBucket = 1;
elseif (thetaDot>=-50&&thetaDot<=50)
	thetaDotBucket = 2;
else
	thetaDotBucket = 3;
end

box = sub2ind([14,3,3,3],thetaBucket, thetaDotBucket,xBucket,xDotBucket);
end
return;


text:
function r = plus(p,q)
% this function adds the x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'CellVector')&&isa(q, 'CellVector'))
    r=p;
    r.xvalue = p.xvalue+q.xvalue;
    r.yvalue = p.yvalue+q.yvalue;
    r.zvalue = p.zvalue+q.zvalue;
elseif isa(p, 'CellVector')
    r=p;
    r.xvalue = p.xvalue+q;
    r.yvalue = p.yvalue+q;
    r.zvalue = p.zvalue+q;
else
    r=q;
    r.xvalue = p+q.xvalue;
    r.yvalue = p+q.yvalue;
    r.zvalue = p+q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool