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

author: jmaih


text:
function c=ref_dftiv(f)
%REF_DFT  Reference Discrete Fourier Transform Type IV
%   Usage:  c=ref_dftiv(f);
%
%   This is highly experimental!

L=size(f,1);
W=size(f,2);

% Create weights.
w=sqrt(1/L);

% Create transform matrix.
F=zeros(L);

for m=0:L-1
  for n=0:L-1
    F(m+1,n+1)=w*exp(2*pi*i*(m+.5)*(n+.5)/L);
  end;
end;

% Compute coefficients.
c=F'*f;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [x,m,e] = RunningAverage(X,Y,varargin)

%RunningAverage - Compute running linear or angular average.
%
% Computes the running average of y=f(x). Variable y can be linear or circular
% (use radians). The error bars are standard errors of the mean for linear
% data, or 95% confidence intervals for circular data.
%
%  USAGE
%
%    [x,m,e] = RunningAverage(x,y,<options>)
%
%    x              x variable (e.g., time)
%    y              values at x
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'window'      averaging window size (default (max-min)/10)
%     'overlap'     overlap between successive windows (default = 0.8*window)
%     'limits'      x limits to use instead of min and max
%     'type'        either 'linear' or 'circular' (default 'linear')
%    =========================================================================
%
%  OUTPUT
%
%    x              new x variable
%    m              running average
%    e              sem for linear variables, otherwise 95% confidence intervals

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% 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.

% Default values
type = 'linear';
limits = [min(X) max(X)];
nBins = 10;
window = 0;
overlap = [];

if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help RunningAverage'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help RunningAverage'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'type',
			type = varargin{i+1};
			if ~isstring_FMAT(type,'linear','circular'),
				error('Incorrect value for property ''type'' (type ''help RunningAverage'' for details).');
			end

		case 'window',
			window = varargin{i+1};
			if ~isdscalar(window,'>0'),
				error('Incorrect value for property ''window'' (type ''help RunningAverage'' for details).');
			end

		case 'overlap',
			overlap = varargin{i+1};
			if ~isdscalar(overlap,'>=0'),
				error('Incorrect value for property ''overlap'' (type ''help RunningAverage'' for details).');
			end

		case 'limits',
			limits = varargin{i+1};
			if ~isdvector(limits,'#2','<'),
				error('Incorrect value for property ''limits'' (type ''help RunningAverage'' for details).');
			end

		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help RunningAverage'' for details).']);

	end
end

x0 = limits(1);
x1 = limits(2);
if window == 0,
	window = (x1-x0)/nBins;
end
if isempty(overlap),
	overlap = 0.8*window;
end

% Loop through data
i = 1;
xi = x0+window/2;
while xi+window/2 <= x1,
	x(i,1) = xi;
	ok = InIntervals(X,xi+[-0.5 0.5]*window);
	if sum(ok) == 0,
		m(i,1) = NaN;
		e(i,:) = [NaN NaN];
	elseif strcmp(type,'circular'),
		[M,C] = CircularConfidenceIntervals(Y(ok));
		m(i,1) = M;
		e(i,:) = C;
	else
		m(i,1) = nanmean(Y(ok));
		n = sum(ok);
		s = nanstd(Y(ok))/sqrt(n);
		e(i,1) = m(i)-s;
		e(i,2) = m(i)+s;
	end
	xi = xi + window - overlap;
	i = i + 1;
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% READEGI - read EGI Simple Binary datafile (versions 2,3,4,5,6,7).
%	      Return header info, EEG data, and any event data.
% Usage:
%   >> [head, TrialData, EventData, CatIndex] = readegi(filename, dataChunks, forceversion)
%
% Required Input:
%   filename = EGI data filename
%
% Optional Input:
%   dataChunks = vector containing desired frame numbers(for unsegmented
%                datafiles) or segments (for segmented files). If this
%                input is empty or is not provided then all data will be
%                returned.
%   forceversion = integer from 2 to 7 which overrides the EGI version read from the
%                file header.  This has been occasionally necessary in cases where
%                the file format was incorrectly indicated in the header.
% Outputs:
%   head = struct containing header info (see READEGIHDR )
%   TrialData = EEG channel data
%   EventData = event codes
%   CatIndex  = segment category indices
%
% Author: Cooper Roddey, SCCN, 13 Nov 2002
%
% Note: code derived from C source code written by
%       Tom Renner at EGI, Inc.
%
% See also: READEGIHDR

% Copyright (C) 2002 , Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function  [head, TrialData, EventData, SegmentCatIndex] = readegi(filename, dataChunks,forceversion)

if nargin <1 || nargin > 3,
    help readegi;
    return;
end

if nargout < 2 || nargout > 4,
	error('2 to 4 output args required');
end

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

if ~isvector(dataChunks)
     error('dataChunks must either be empty or a vector');
end

[fid,message] = fopen(filename,'rb','b');
if (fid < 0),
   error(message);
end

% get our header structure
fprintf('Importing binary EGI data file ...\n');
if exist('forceversion')
   head = readegihdr(fid,forceversion);
else
   head = readegihdr(fid);
end

% do we have segmented data?
segmented = 0;
switch(head.version),
    case {2,4,6}
       segmented = 0;
    case {3,5,7}
       segmented = 1;
end

% each type of event has a dedicated "channel"
FrameVals = head.nchan + head.eventtypes;

if (segmented)

    desiredSegments = dataChunks;
    if isempty(desiredSegments),
       desiredSegments = [1:head.segments];
    end

    nsegs = length(desiredSegments);

    TrialData = zeros(FrameVals,head.segsamps*nsegs);
    readexpected = FrameVals*head.segsamps*nsegs;

else

    desiredFrames = dataChunks;
    if isempty(desiredFrames),
       desiredFrames = [1:head.samples];
    end

    nframes = length(desiredFrames);

    TrialData = zeros(FrameVals,nframes);
    readexpected = FrameVals*nframes;

end

% get datatype from version number
switch(head.version)
   case {2,3}
      datatype = 'integer*2';
   case {4,5}
      datatype = 'float32';
   case {6,7}
      datatype = 'float64';
   otherwise,
      error('Unknown data format');
end

% read in epoch data
readtotal = 0;
j=0;
SegmentCatIndex = [];
if (segmented),

    for i=1:head.segments,
	segcatind = fread(fid,1,'integer*2');
	segtime = fread(fid,1,'integer*4');

	[tdata, count] = ...
	   fread(fid,[FrameVals,head.segsamps],datatype);

       % check if this segment is one of our desired ones
        if ismember(i,desiredSegments)

           j=j+1;
           SegmentCatIndex(j) = segcatind;
           SegmentStartTime(j) = segtime;

           TrialData(:,[1+(j-1)*head.segsamps:j*head.segsamps]) = tdata;

           readtotal = readtotal + count;
        end

	if (j >= nsegs), break; end

    end

else
    % read unsegmented data

    % if dataChunks is empty, read all frames
    if isempty(dataChunks)

        [TrialData, readtotal] = fread(fid, [FrameVals,head.samples],datatype);

    else   % grab only the desiredFrames
           % This could take a while...

        for i=1:head.samples,

          [tdata, count] = fread(fid, [FrameVals,1],datatype);

          % check if this segment is a keeper
          if ismember(i,desiredFrames),
              j=j+1;
              TrialData(:,j) = tdata;
              readtotal = readtotal + count;
          end

          if (j >= nframes), break; end

        end
    end
end

fclose(fid);

if ~isequal(readtotal, readexpected)
     error('Number of values read not equal to the number expected.');
end

EventData = [];
if (head.eventtypes > 0),
	EventData = TrialData(head.nchan+1:end,:);
	TrialData = TrialData(1:head.nchan,:);
end

% convert from A/D units to microvolts
if ( head.bits ~= 0 && head.range ~= 0 )
       TrialData = (head.range/(2^head.bits))*TrialData;
end

% convert event codes to char
% ---------------------------
head.eventcode = char(head.eventcode);

%--------------------------- ISVECTOR ---------------
function retval = isvector(V)

   s = size(V);
   retval = ( length(s) < 3 ) & ( min(s) <= 1 );

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

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% CREATE_SRC Create a source of files & objects
%
% Usage
%    src = CREATE_SRC(directory, objects_fun)
%
% Input
%    directory (char): The directory in which the files are found.
%    objects_fun (function handle): Given a filename, objects_fun returns its
%       constituent objects and their respective classes.
%
% Output
%    src (struct): The source corresponding to all data files (.jpg, .wav,
%       .au) contained in directory and their objects as defined by
%       objects_fun.
%
% See also
%    PREPARE_DATABASE

function src = create_src(directory,objects_fun)
	if nargin < 1
		error('Must specify directory!');
	end

	if nargin < 2
		objects_fun = @default_objects_fun;
	end

	files = find_files(directory);

	if isempty(files)
		error('No data files found in the specified directory!');
	end

	objects = [];
	classes = [];

	for ind = 1:length(files)

		[file_objects,file_classes] = objects_fun(files{ind});

		if ~isempty(file_objects)
			file_objects = arrayfun(@(s)(setfield(s,'ind',ind)),file_objects);
		end

		objects = [objects file_objects];
		classes = [classes file_classes];
	end

	file_inds = [objects.ind];

	[file_inds,I,J] = unique(file_inds);

	J = num2cell(J);

	[objects.ind] = J{:};

	files = files(file_inds);

	fields = fieldnames(objects)';
	fields(strmatch('ind',fields)) = [];
	fields(strmatch('u1',fields)) = [];
	fields(strmatch('u2',fields)) = [];

	objects = orderfields(objects,[{'ind','u1','u2'} fields]);

	[classes,temp,obj_class] = unique(classes);

	obj_class = num2cell(obj_class);

	[objects.class] = obj_class{:};

	src.classes = classes;
	src.files = files;
	src.objects = objects;
end

function files = find_files(directory)
	extensions = {'au','wav','jpg','png'};

	dir_list = dir(directory);

	files = {};

	for k = 1:length(dir_list)
		name = dir_list(k).name;

		% Hidden file or current/upper directory? Skip.
		if name(1) == '.'
			continue;
		end

		% Depending on file type, recurse or add audio file.
		if dir_list(k).isdir
			files = [files find_files(fullfile(directory,name))];
		else
			found = 0;
			for l = 1:length(extensions)
				len = length(extensions{l});
				if length(name) > len+1 && ...
				   strcmpi(name(end-len:end),['.' extensions{l}])
					found = l;
					break;
				end
			end
			if found > 0
				files = [files fullfile(directory,name)];
			end
		end
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function fixed_num = p09_fixed_num ( )

%*****************************************************************************80
%
%% P09_FIXED_NUM returns the number of fixed points in problem 09.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer FIXED_NUM, the number of fixed points.
%
  fixed_num = 16;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function alpha = sphericalAngle(p1, p2, p3)
%SPHERICALANGLE Compute angle between points on the sphere.
%
%   ALPHA = sphericalAngle(P1, P2, P3)
%   Computes angle (P1, P2, P3), i.e. the angle, measured at point P2,
%   between the direction (P2, P1) and the direction (P2, P3).
%   The result is given in radians, between 0 and 2*PI.
%
%   Points are given either as [x y z] (there will be normalized to lie on
%   the unit sphere), or as [phi theta], with phi being the longitude in [0
%   2*PI] and theta being the elevation on horizontal [-pi/2 pi/2].
%
%
%   NOTE:
%   this is an 'oriented' version of the angle computation, that is, the
%   result of sphericalAngle(P1, P2, P3) equals
%   2*pi-sphericalAngle(P3,P2,P1). To have the more classical relation
%   (with results given betwen 0 and PI), it suffices to take the minimum
%   of angle and 2*pi-angle.
%
%   Examples
%     % Use inputs as cartesian coordinates
%     p1 = [0 1 0];
%     p2 = [1 0 0];
%     p3 = [0 0 1];
%     alpha = sphericalAngle(p1, p2, p3)
%     alpha =
%         1.5708
%
%     % Use inputs as spherical coordinates
%     sph1 = [.1 0];
%     sph2 = [0 0];
%     sph3 = [0 .1];
%     alphas = sphericalAngle(sph1, sph2, sph3)
%     alphas =
%         1.5708
%
%
%   See also:
%   geom3d, angles3d, spheres, sph2cart

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005.
%

%   HISTORY
%   23-05-2006 fix bug for points with angle from center > pi/2
%   05-06-2013 fix bug for points given as spherical coordinates, better
%       support for multiple inputs

% test if points are given as matlab spherical coordinates
if size(p1, 2) == 2
    [x, y, z] = sph2cart(p1(:,1), p1(:,2), ones(size(p1,1), 1));
    p1 = [x y z];
    [x, y, z] = sph2cart(p2(:,1), p2(:,2), ones(size(p2,1), 1));
    p2 = [x y z];
    [x, y, z] = sph2cart(p3(:,1), p3(:,2), ones(size(p3,1), 1));
    p3 = [x y z];
end

% normalize points
p1  = normalizeVector3d(p1);
p2  = normalizeVector3d(p2);
p3  = normalizeVector3d(p3);

% create the plane tangent to the unit sphere and containing central point
plane = createPlane(p2, p2);

% project the two other points on the plane
pp1 = planePosition(projPointOnPlane(p1, plane), plane);
pp3 = planePosition(projPointOnPlane(p3, plane), plane);

% compute angle on the tangent plane
pp2 = zeros(max(size(pp1, 1), size(pp3,1)), 2);
alpha = angle3Points(pp1, pp2, pp3);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ errors ] = test_sparsify(  )
    gsp_reset_seed(0);

errors = 0;
errors = errors + test1();

errors = errors + test2();

end

function [ errors ] = test1( )

errors = 0;
try
    epsilon = 0.4;
    param.distribute = 1;
    param.Nc = 20;
    G = gsp_sensor(256,param);
    G2 = gsp_graph_sparsify(G,epsilon);
    figure(100);
    gsp_plot_graph(G);
    title('Original graph')
    figure(101);
    gsp_plot_graph(G2);
    title('Sparsified graph')
    close(100);
    close(101);
   fprintf('SPARSIFY: test 1 ok\n');
catch
    errors = errors +1;
    warning('SPARSIFY: test 1 error')
end

end

function [ errors ] = test2( )

errors = 0;
    N = 100;
    epsilon = 0.6;
    G = gsp_sensor(N);
    gsp_reset_seed(0);
    G2 = gsp_graph_sparsify(G,epsilon);
    gsp_reset_seed(0);
    L = gsp_graph_sparsify_old(G.L,epsilon);

if sum(sum(abs(G2.L-L)))<1e-10

   fprintf('SPARSIFY: test  ok\n');
else
    errors = errors+1;
    badness = sum(sum(abs(G2.L-L)))
    warning('SPARSIFY: test 2 error')
end

end

repo: gspbox


text:
function [ a, b] = p03_interval ( )

%*****************************************************************************80
%
%% P03_INTERVAL returns a starting interval for optimization for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, two points defining an interval in which
%    the local minimizer should be sought.
%
  a = -2.0;
  b = +2.0;

  return
end


text:
function g = graf(matriks)

n=length(matriks);
theta = 2*pi./n;
angle = theta*(0:(n-1));
radius = n* ones(1,length(angle));
textRadius =  n*( ones(1,length(angle)) + 0.2);
x = radius .* cos(angle);
y = radius .* sin(angle);
textX = textRadius .* cos(angle);
textY =  textRadius .* sin(angle);
fHandle = figure;
figure(fHandle);
dummy = 1:1:n;
edges = combntns(dummy,2);

mat1 = zeros(n,n);
mat2 = zeros(n,n);
for i=1:1:n
    for j=1:1:n
        if matriks(i,j)==1
            mat1(i,j)=1;             % edge yang positif
        end
        if matriks(i,j)==-1
            mat2(i,j)=1;              % edge yang negatif
        end
    end
end

for m = 1:n
    h = rectangle('Position',[x(m)-0.5, y(m)-0.5, 1, 1], 'Curvature', [1,1], 'FaceColor','k');
end

grid off;
axis off;

for j=1:n
    text(textX(j), textY(j), int2str(j),'Color','black');
    nodeUpdated(j).output = [];
    indices = find(mat1(j,:));
    for k=1:length(indices)
        multiplicity(k) = matriks(j,indices(k));
    end

    for m=1:length(indices)
        nodeUpdated(j).output = [nodeUpdated(j).output repmat(indices(m),1,multiplicity(m))];
        line([x(j) x(indices(m))] , [y(j) y(indices(m))],'LineWidth',3);
    end
end

for j=1:n
    nodeUpdated(j).output = [];
    indices = find(mat2(j,:));
    for k=1:length(indices)
        multiplicity(k) = matriks(j,indices(k));
    end

    for m=1:length(indices)
        nodeUpdated(j).output = [nodeUpdated(j).output repmat(indices(m),1,multiplicity(m))];
        line([x(j) x(indices(m))] , [y(j) y(indices(m))],'Color','k','LineStyle',':','LineWidth',3);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function e2 = fem2d_l2_error_serene ( nx, ny, x, y, u, exact )

%*****************************************************************************80
%
%% FEM2D_L2_ERROR_SERENE: L2 error norm of a finite element solution.
%
%  Discussion:
%
%    The finite element method has been used, over a rectangle,
%    involving a grid of NXxNY nodes, with serendipity elements used
%    for the basis.
%
%    The finite element coefficients have been computed, and a formula for the
%    exact solution is known.
%
%    This function estimates E2, the L2 norm of the error:
%
%      E2 = Integral ( X, Y ) ( U(X,Y) - EXACT(X,Y) )^2 dX dY
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, the number of nodes in the X and Y directions.
%
%    Input, real X(NX), Y(NY), the grid coordinates.
%
%    Input, real U(*), the finite element coefficients.
%
%    Input, function EQ = EXACT(X,Y), returns the value of the exact
%    solution at the point (X,Y).
%
%    Output, real E2, the estimated L2 norm of the error.
%
  e2 = 0.0;
%
%  Quadrature definitions.
%
  quad_num = 3;
  abscissa(1) = -0.774596669241483377035853079956;
  abscissa(2) = 0.000000000000000000000000000000;
  abscissa(3) = 0.774596669241483377035853079956;
  weight(1) = 0.555555555555555555555555555556;
  weight(2) = 0.888888888888888888888888888889;
  weight(3) = 0.555555555555555555555555555556;
%
%  Compute the matrix entries by integrating over each element.
%
  ex_num = ( nx - 1 ) / 2;
  ey_num = ( ny - 1 ) / 2;

  for ey = 1 : ey_num

    s = 2 * ey - 1;
    mm = 2 * ey;
    n = 2 * ey + 1;

    ys = y(s);
    ym = y(mm);
    yn = y(n);

    yy(1) = y(n);
    yy(2) = y(n);
    yy(3) = y(n);
    yy(4) = y(mm);
    yy(5) = y(s);
    yy(6) = y(s);
    yy(7) = y(s);
    yy(8) = y(mm);

    for ex = 1 : ex_num

      w = 2 * ex - 1;
      cc = 2 * ex;
      e = 2 * ex + 1;

      xe = x(e);
      xc = x(cc);
      xw = x(w);

      xx(1) = x(e);
      xx(2) = x(cc);
      xx(3) = x(w);
      xx(4) = x(w);
      xx(5) = x(w);
      xx(6) = x(cc);
      xx(7) = x(e);
      xx(8) = x(e);
%
%  Node indices
%
%  3  2  1  wn  cn  en
%  4     8  wm      em
%  5  6  7  ws  cs  es
%
      node(1) = ( 3 * ey     ) * ey_num + 2 * ey + 2 * ex + 1;
      node(2) = ( 3 * ey     ) * ey_num + 2 * ey + 2 * ex;
      node(3) = ( 3 * ey     ) * ey_num + 2 * ey + 2 * ex - 1;
      node(4) = ( 3 * ey - 1 ) * ey_num + 2 * ey +     ex - 1;
      node(5) = ( 3 * ey - 3 ) * ey_num + 2 * ey + 2 * ex - 3;
      node(6) = ( 3 * ey - 3 ) * ey_num + 2 * ey + 2 * ex - 2;
      node(7) = ( 3 * ey - 3 ) * ey_num + 2 * ey + 2 * ex - 1;
      node(8) = ( 3 * ey - 1 ) * ey_num + 2 * ey +     ex;

      for qx = 1 : quad_num

        xq = ( ( 1.0 - abscissa(qx) ) * x(e)   ...
             + ( 1.0 + abscissa(qx) ) * x(w) ) ...
               / 2.0;

        for qy = 1 : quad_num

          yq = ( ( 1.0 - abscissa(qy) ) * y(n)   ...
               + ( 1.0 + abscissa(qy) ) * y(s) ) ...
                 / 2.0;

          wq = weight(qx) * ( x(e) - x(w) ) / 2.0 ...
             * weight(qy) * ( y(n) - y(s) ) / 2.0;

          v = basis_serene ( xq, yq, xw, ys, xe, yn, xx, yy );

          uq = 0.0;
          for k = 1 : 8
            uq = uq + u(node(k)) * v(k);
          end

          eq = exact ( xq, yq );
          e2 = e2 + wq * ( uq - eq )^2;

        end
      end
    end
  end

  e2 = sqrt ( e2 );

  return
end


text:
%%%%%%%%%%%%%%%
%%
if true,
    clear all;close all;
    varList = {'CR_dtr','numcell_full','CellXYZ','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn'};
else
    % or without 'CR_dtr', keeping it from previous step:
    clearvars -except 'CR_dtr'; % 'CR_z_dtr';
    varList = {'numcell_full','CellXYZ','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn'};
end

%%
M_dir = GetFishDirectories();
M_stimset = [1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2];

save_dir = GetCurrentDataDir();

dshift = 2; % differential shift between fictive and fluo, manually chosen

%% MANUAL
for i_fish = 1,
    disp(['i_fish = ', num2str(i_fish)]);
    tic
    %% load data
    load(fullfile(M_dir{i_fish},['Fish' num2str(i_fish) '_direct_load_full_v6.mat']),varList{:}); % '_direct_load_nodiscard.mat'
%     load(fullfile(datadir,'frame_turn.mat'),'frame_turn');

    if size(frame_turn,1)<size(frame_turn,2),
        frame_turn = frame_turn';
    end

    %% Manual occasional frame corrections: ONLY RUN ONCE!!!!!!!!!

    if false,
        % add 1 frame at start: ------------- for which fish again??
        CR_dtr = horzcat(CR_dtr(:,1),CR_dtr);
        frame_turn = vertcat(frame_turn(1,:),frame_turn);
    end

    % correction of an error in Fish #1
    if i_fish == 1,
        CR_dtr = horzcat(CR_dtr(:,1:20),CR_dtr);
%         CR_z_dtr = horzcat(CR_z_dtr(:,1:20),CR_z_dtr);
    end

    %% index processing
    if M_stimset(i_fish)==1, % fish 1-7
        M_period = {480,160,160,320,480,140,150}; %,{120,150,360}};
        periods = M_period{i_fish};
        period = periods;

        nrep = floor(size(CR_dtr,2)/period)-1;

        shift = period; % (fixed 9/23/15; before = period+dshift)

        IX_all = 1+shift:period*nrep+shift;
        CellResp = CR_dtr(:,IX_all);
%         CellRespZ = CR_z_dtr(:,IX_all); % z-scored

        nCells = size(CR_dtr,1);
        CellRespAvr = mean(reshape(CR_dtr(:,IX_all+dshift),nCells,period,[]),3);
%         CellRespAvrZ = mean(reshape(CR_z_dtr(:,IX_all+dshift),nCells,period,[]),3);

        if i_fish<6,
            stimrangenames = {'16 permut: B/W/phototaxis*2'};
        else % i_fish = 6 or 7
            stimrangenames = {'phototaxis'};
        end
        %         stimrangenames = {'rep average','all reps','rep #1','rep #2','last rep'};
        ix_avr = 1:period;
        ix_all = 1:period*nrep;
        tlists = {ix_avr, ix_all};
        for i = [1,2,nrep],
            ix = 1+period*(i-1):period*i;
            tlists = [tlists, ix];
        end

        stim_full = frame_turn(IX_all,17)';
        stim_full = round(stim_full);

        stimAvr = stim_full(1:period);

    else % multiple protocols
        % process raw stimulus code
        [stimset,stim_full_raw] = StimulusKeyPreprocessing(frame_turn,i_fish); % StimulusKeyPreprocessing(frame_turn,i_fish,'isplotting') to show plot

        %% variables to save later in struct
        periods = [stimset.period];
        stimrangenames = {stimset.name};

        %% find time-lists = list of time-frames numbers for a certain stimulus
        nTypes = length(stimset);
        tlists_raw = cell(1,nTypes+1); % time-lists, i.e. frame indices
        IX_all = [];
        for i_type = 1:nTypes,
            nSets = length(stimset(i_type).starts);
            IX = [];
            for i_set = 1:nSets,
                IX = horzcat(IX,stimset(i_type).starts(i_set):stimset(i_type).stops(i_set));
            end
            tlists_raw{i_type} = IX;
            IX_all = horzcat(IX_all,IX);
        end
        tlists_raw{nTypes+1} = IX_all;
%         tlists_raw{nTypes+2} = sort(IX_all);

        % This is the main data to store and load to GUI
        CellResp = CR_dtr(:,IX_all);
%         CellRespZ = CR_z_dtr(:,IX_all);
        stim_full = stim_full_raw(IX_all);

        %% get tlists corresponding to IX_all (tlists_raw ~ 1:totalframe#)
        tlists = cell(1,nTypes+1); % time-lists, i.e. frame indices
        for i_type = 1:nTypes+1,
            [~,IX] = ismember(tlists_raw{i_type},IX_all);
            tlists{i_type} = IX(find(IX));
        end

        %% find average
        numcell = size(CR_dtr,1);
        shift = -dshift; % circshift fluo left by 2 frames

        CellRespAvr = [];
        CellRespAvrZ = [];
        stimAvr = [];
        for i = 1:nTypes,
            M = circshift(CellResp(:,tlists{i}),shift,2);
            CellRespAvr = horzcat(CellRespAvr,mean(reshape(M,numcell,periods(i),[]),3));
%             CellRespAvrZ = horzcat(CellRespAvrZ,mean(reshape(M,numcell,periods(i),[]),3));
            % stim from stim_full
            m = stim_full(tlists{i});
            stimAvr = horzcat(stimAvr,m(1:periods(i)));
        end

    end

    %% prepare fictive data
    rows = [7,8,9,13,14];
    F = frame_turn(:,rows)';

    F(2,:) = -F(2,:);
    Fictive = F(:,IX_all);

    % find averages
    if M_stimset(i_fish)==1, % fish 1-7
        FictiveAvr = mean(reshape(Fictive,length(rows),period,[]),3);
    else
        FictiveAvr = [];
        for i = 1:nTypes,
            avr = mean(reshape(F(:,tlists_raw{i}),length(rows),periods(i),[]),3);
            FictiveAvr = horzcat(FictiveAvr,avr);
        end
    end

    % normalizations
    for i = 1:3,
        m = FictiveAvr(i,:);
        FictiveAvr(i,:) = (m-min(m))/(max(m)-min(m));
        m = Fictive(i,:);
        Fictive(i,:) = (m-min(m))/(max(m)-min(m));
    end
    m = FictiveAvr(4:5,:);
    FictiveAvr(4:5,:) =  (m-min(min(m)))/(max(max(m))-min(min(m)));
    m = Fictive(4:5,:);
    Fictive(4:5,:) =  (m-min(min(m)))/(max(max(m))-min(min(m)));

    %% compile CONST
    CellRespAvrZ = zscore(CellRespAvr')';

    data = [];
    names = {'numcell_full','CellXYZ','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','periods','shift','dshift',...
        'CellRespAvr','CellRespAvrZ','Fictive','FictiveAvr','stim_full','stimAvr',... %'CellResp',
        'tlists','stimrangenames'};
    if M_stimset(i_fish) > 1,
        names = [names,{'tlists_raw','stimset'}];
    end

    for i = 1:length(names), % use loop to save variables into fields of CONST
        eval(['data.',names{i},'=', names{i},';']);
    end

    %% and save

    %%% new method with partitioning of main data
    newfishdir = fullfile(save_dir,['Data_F' num2str(i_fish) '_full.mat']);
    dimCR = size(CellResp);
    save(newfishdir,'data','dimCR','-v6');
    % custom function:
    SaveFileInPartsAppendv6(newfishdir,CellResp,'CellResp');
%     SaveFileInPartsAppendv6(newfishdir,CellRespZ,'CellRespZ');

    toc; beep

end

%% compute z-score
% disp('compute z-score')
% cell_resp_z_full = zscore(cell_resp_full')';

author: xiuyechen

repo: FishExplorer


text:
function [ a, more ] = tree_rb_lex_next ( n, a, more )

%*****************************************************************************80
%
%% TREE_RB_LEX_NEXT generates rooted binary trees in lexicographic order.
%
%  Discussion:
%
%    The information definining the tree of N nodes is stored in a vector
%    of 0's and 1's, in preorder traversal form.  Essentially, the
%    shape of the tree is traced out with a pencil that starts at the root,
%    and stops at the very last null leaf.  The first time that a (non-null)
%    node is encountered, a 1 is added to the vector, and the left
%    descendant of the node is visited next.  When the path returns from
%    the first descendant, the second descendant is visited.  When then path
%    returns again, the path passes back up from the node to its parent.
%    A null leaf is encountered only once, and causes a zero to be added to
%    the vector, and the path goes back up to the parent node.
%
%    The lexicographic order is used to order the vectors of 1's and 0's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2013
%
%  Reference:
%
%    Frank Ruskey,
%    Combinatorial Generation,
%    To appear.
%
%  Parameters:
%
%    Input, integer N, the number of nodes in the rooted binary
%    tree.  N should be odd.
%
%    Input/output, integer A(N), the preorder traversal form for
%    the previous/next rooted binary tree.
%
%    Output, logical MORE, is TRUE if the next rooted binary tree was
%    returned on this call, or FALSE if there are no more rooted binary
%    trees, and the output of the previous call was the last one.
%
  if ( ~ more )
    a(1:2:n-2) = 1;
    a(2:2:n-1) = 0;
    a(n) = 0;
    more = 1;
    return
  end
%
%  Find the last 1 in A.
%
  k = n;
  while ( a(k) == 0 )
    k = k - 1;
  end
  q = n - k - 1;
%
%  Find the last 0 preceding the last 1 in A.
%  If there is none, then we are done, because 11...1100..00
%  is the final element.
%
  while ( 1 )

    if ( k == 1 )
      more = 0;
      return
    end

    if ( a(k) == 0 )
      break
    end

    k = k - 1;

  end

  p = n - k - q - 1;
  a(k) = 1;
  a(k+1:n-2*p+1) = 0;
  a(n-2*p+2:2:n-2) = 1;
  a(n-2*p+3:2:n-1) = 0;
  a(n) = 0;

  return
end


text:
function test_issue1067

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_prepare_layout ft_prepare_mesh ft_prepare_sourcemodel ft_prepare_headmodel ft_prepare_neighbours

% For some fields there are multiple ways to obtain them. The desired use in the
% high-level functions is to call the function that is the highest in the hierarchy.
%
% ft_prepare_xxx -> ft_fetch_xxx -> ft_read_xxx/loadvar

%%

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/issue1067'));

%% layout

cfg = [];
cfg.layout = 'CTF151.lay';
cfg.layout = ft_prepare_layout(cfg); % using FT_PREPARE_XXX
assert(isstruct(cfg.layout));
cfg.layout = ft_prepare_layout(cfg); % using FT_PREPARE_XXX again

cfg = [];
cfg.layout = 'CTF151.mat';
cfg.layout = ft_prepare_layout(cfg); % using FT_PREPARE_XXX
assert(isstruct(cfg.layout));
cfg.layout = ft_prepare_layout(cfg); % using FT_PREPARE_XXX again

%% headshape

cfg = [];
cfg.headshape = 'Subject01.shape';
cfg.headshape = ft_prepare_mesh(cfg); % using FT_PREPARE_XXX
assert(isstruct(cfg.headshape));
cfg.headshape = ft_prepare_mesh(cfg); % using FT_PREPARE_XXX again

cfg = [];
cfg.headshape = 'Subject01_headshape.mat';
cfg.headshape = ft_prepare_mesh(cfg); % using FT_PREPARE_XXX
assert(isstruct(cfg.headshape));
cfg.headshape = ft_prepare_mesh(cfg); % using FT_PREPARE_XXX again

%% sourcemodel

cfg = [];
cfg.sourcemodel = 'Subject01_sourcemodel.mat';
cfg.sourcemodel = ft_prepare_sourcemodel(cfg); % using FT_PREPARE_XXX
assert(isstruct(cfg.sourcemodel));
cfg.sourcemodel = ft_prepare_sourcemodel(cfg); % using FT_PREPARE_XXX again

%% headmodel

cfg = [];
cfg.headmodel = 'Subject01.hdm';
cfg.headmodel = ft_prepare_headmodel(cfg); % using FT_PREPARE_XXX
assert(isstruct(cfg.headmodel));
cfg.headmodel = ft_prepare_headmodel(cfg); % using FT_PREPARE_XXX again

cfg = [];
cfg.headmodel = 'Subject01_headmodel.mat';
cfg.headmodel = ft_prepare_headmodel(cfg); % using FT_PREPARE_XXX
assert(isstruct(cfg.headmodel));
cfg.headmodel = ft_prepare_headmodel(cfg); % using FT_PREPARE_XXX again

%% neighbours

cfg = [];
cfg.neighbours = 'ctf151_neighb.mat';
cfg.neighbours = ft_prepare_neighbours(cfg); % using FT_PREPARE_XXX
assert(isstruct(cfg.neighbours));
cfg.neighbours = ft_prepare_neighbours(cfg); % using FT_PREPARE_XXX again

%% grad

cfg = [];
cfg.grad = 'Subject01_grad.mat';
cfg.grad = ft_fetch_sens(cfg); % using FT_FETCH_XXX
assert(isstruct(cfg.grad));
cfg.grad = ft_fetch_sens(cfg); % using FT_FETCH_XXX again

%% elec

cfg = [];
cfg.elec = 'Subject01_elec.mat';
cfg.elec = ft_fetch_sens(cfg); % using FT_FETCH_XXX
assert(isstruct(cfg.elec));
cfg.elec = ft_fetch_sens(cfg); % using FT_FETCH_XXX again

%% opto

cfg = [];
cfg.opto = 'opto.mat';
cfg.opto = ft_fetch_sens(cfg); % using FT_FETCH_XXX
assert(isstruct(cfg.opto));
cfg.opto = ft_fetch_sens(cfg); % using FT_FETCH_XXX again

%% event

cfg = [];
cfg.event = 'Subject01_event.mat';
cfg.event = ft_read_event(cfg.event); % using FT_READ_XXX
assert(isstruct(cfg.event));

%% mri

cfg = [];
cfg.mri = 'Subject01.mri';
cfg.mri = ft_read_mri(cfg.mri);  % using FT_READ_XXX
assert(isstruct(cfg.mri));

cfg = [];
cfg.mri = 'Subject01_mri.mat';
cfg.mri = ft_read_mri(cfg.mri);  % using FT_READ_XXX
assert(isstruct(cfg.mri));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Degree-preserving random rewiring.
% Every rewiring decreases the assortativity (pearson coefficient).
%
% Note 1: There are rare cases of neutral rewiring (pearson coefficient stays the same within numerical error).
% Note 2: Assume unweighted undirected graph.
%
% INPUTS: edge list, el and number of rewirings, k (integer)
% OUTPUTS: rewired edge list
%
% Other routines used: degrees.m, edgeL2adj.m
% GB: last updated, Sep 27 2012

function el = rewireDisassort(el,k)

[deg,~,~]=degrees(edgeL2adj(el));

rew=0;

while rew<k

    % pick two random edges
    ind = randi(length(el),1,2);
    edge1=el(ind(1),:); edge2=el(ind(2),:);

    if length(intersect(edge1(1:2),edge2(1:2)))>0; continue; end % the two edges cannot overlap

    nodes=[edge1(1) edge1(2) edge2(1) edge2(2)];
    [~,Y]=sort(deg(nodes));

    % connect nodes(Y(1))-nodes(Y(4)) and nodes(Y(2))-nodes(Y(3))
    if ismember([nodes(Y(1)),nodes(Y(4)),1],el,'rows') || ismember([nodes(Y(2)),nodes(Y(3)),1],el,'rows'); continue; end

    el(ind(1),:)=[nodes(Y(1)),nodes(Y(4)),1];
    el(ind(2),:)=[nodes(Y(2)),nodes(Y(3)),1];

    [~,inds1] = ismember([edge1(2),edge1(1),1],el,'rows');
    el(inds1,:)=[nodes(Y(4)),nodes(Y(1)),1];

    [~,inds2] = ismember([edge2(2),edge2(1),1],el,'rows');
    el(inds2,:)=[nodes(Y(3)),nodes(Y(2)),1];

    rew=rew+1;

end

author: aeolianine

repo: octave-networks-toolbox


text:
function [a,count] = nanplus(a,b,count,alpha)

if nargin < 4
  alpha = 1;
elseif isnan(alpha)
  return;
end
ind = ~isnan(b);

a(ind) = a(ind) + alpha * b(ind);

if nargin > 2
  count = count + alpha * ind;
end

author: mtex-toolbox

repo: mtex

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


text:
function [jac,invjac,phi,dphidx,dphidy] = deriv(s,t,xl,yl)
%deriv  evaluates derivatives of bilinear shape functions
%   [jac,invjac,phi,dphidx,dphidy] = deriv(s,t,xl,yl);
%   input
%          s         reference element x coordinate
%          t         reference element y coordinate
%          xl        physical element x vertex coordinates
%          yl        physical element y vertex coordinates
%   output
%          jac       elementwise jacobian (evaluated at (s,t))
%          invjac    elementwise inverse of jacobian
%          phi       elementwise shape functions
%          dphidx    x derivatives of phi
%          dphidy    y derivatives of phi
%
%   IFISS function: DJS; 4 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
      nel=length(xl(:,1));
      zero_v = zeros(nel,1);
      one_v = ones(nel,1);
%
% evaluate shape functions
     [phi_e,dphids,dphidt] = shape(s,t);
%
      dxds = zero_v;
      dxdt = zero_v;
      dyds = zero_v;
      dydt = zero_v;
	   jac = zero_v;
    invjac = zero_v;
%
      for ivtx = 1:4
         dxds(:) = dxds(:) + xl(:,ivtx) .* one_v*dphids(ivtx);
         dxdt(:) = dxdt(:) + xl(:,ivtx) .* one_v*dphidt(ivtx);
         dyds(:) = dyds(:) + yl(:,ivtx) .* one_v*dphids(ivtx);
         dydt(:) = dydt(:) + yl(:,ivtx) .* one_v*dphidt(ivtx);
      end
%
      jac(:) = dxds(:).*dydt(:) - dxdt(:).*dyds(:);
%
% check element Jacobian
      if any(jac < 1e-9)
      fprintf('Bad element warning ...\n')
         if any(jac <= 0.0)
         error('singular Jacobian ... Aborted ...')
         end
      end
      invjac(:) = one_v ./ jac(:);
%
      for ivtx = 1:4
         phi(:,ivtx) = phi_e(ivtx)*one_v;
		 dphidx(:,ivtx) = ( dphids(:,ivtx).*dydt(:) - dphidt(:,ivtx).*dyds(:));
         dphidy(:,ivtx) = (-dphids(:,ivtx).*dxdt(:) + dphidt(:,ivtx).*dxds(:));
      end
      return


text:
function linpack_s_test30 ( )

%*****************************************************************************80
%
%% TEST30 tests STRDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  lda = n;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST30\n' );
  fprintf ( 1, '  For a triangular matrix,\n' );
  fprintf ( 1, '  STRDI computes the determinant or inverse.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Lower triangular matrix A.
%
  [ a, seed ] = r4mat_uniform_01 ( n, n, seed );

  for i = 1 : n
    for j = i+1 : n
      a(i,j) = 0.0;
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Lower triangular matrix A:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %14f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  job = 110;

  [ a, det, info ] = strdi ( a, lda, n, job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The determinant = %f * 10 ^ %f\n ', det(1), det(2) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The inverse matrix:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %14f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end

%
%  Upper triangular matrix A.
%
  [ a, seed ] = r4mat_uniform_01 ( n, n, seed );

  for i = 1 : n
    for j = 1 : i - 1
      a(i,j) = 0.0;
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Upper triangular matrix A:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %14f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  job = 111;

  [ a, det, info ] = strdi ( a, lda, n, job );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The determinant = %f * 10 ^ %f\n ', det(1), det(2) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The inverse matrix:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %14f', a(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function [fMshift, fProbability, fAICc, mProblikelihood, bH, bH2] = calc_loglikelihood_dM2(mCat1, mCat2)
% function [fMshift, fProbability, fAICc, mProblikelihood, bH] = calc_loglikelihood_dM2(mCat1, mCat2);
% ----------------------------------------------------------------------------------------------------
% Calculate log-likelihood estimation of magnitude shift dM between to periods; use power-law
% approximation to first time period for comparison
%
% Incoming variable
% mCat1 : EQ catalog period 1
% mCat2 : EQ catalog period 2 (Observed catalog)
%
% Outgoing variable
% fProbability    : log-likelihood probabilty
% fMshift         : Magnitude shift with the lowest  max. lieklihood score
% fAICc           : Corrected Akaike Information Criterion
% mProblikelihood : Solution matrix shift and likelihood score
% bH              : Result of KSTEST2 hypothesis test at 0.05 significance level
% bH2             : Result of KSTEST2 hypothesis test at 0.05 significance level on EMR-model
%
% Author: J. Woessner, woessner@seismo.ifg,.ethz.ch
% last update: 02.02.04

% Initialize
mProblikelihood = [];
vfProbability = [];
vMshift = [];
vMc = [];
fBinning = 0.1;

% Determine exact time period
fPeriod1 = max(mCat1(:,3)) - min(mCat1(:,3));
fPeriod2 = max(mCat2(:,3)) - min(mCat2(:,3));

mCat1Mod = mCat1;

for fMshift = -0.4:0.1:0.4
    % Apply shift
    mCat1Mod(:,6) = mCat1(:,6)+fMshift;
    % Initialize values
    fMinMag = floor(min([min(mCat1Mod(:,6)) min(mCat2(:,6))]));
    fMaxMag = ceil(max([max(mCat1Mod(:,6)) max(mCat2(:,6))]));
    %% Calculate model for best fitting Mc
     try
        [mResult, fMls, fMc, fMu, fSigma, mDataPred, vPredBest, fBvalue] = calc_McCdfnormal(mCat1Mod, fBinning);
        vPredFMD = mDataPred(:,1)'./fPeriod1;
        vMags = mDataPred(:,2)';
        [vModFMD,vBin1] = hist(mCat1Mod(:,6),fMinMag:0.1:fMaxMag);
        [vObsFMD,vBin2] = hist(mCat2(:,6),min(vMags):0.1:max(vMags));
%             figure_w_normalized_uicontrolunits(40)
%             plot(vBin2, vObsFMD./fPeriod2,'bo')
%             hold on
%             plot(vBin1, vModFMD./fPeriod1,'go')
%             plot(vMags,vPredFMD,'r*')
%             hold off
%             drawnow
        % Time normalization and round due to Poisson distribution calculation in calc_log10poisspdf
        vObsFMD = ceil(vObsFMD./fPeriod2);
        % Calculate the likelihoods for both models
        vProb_ = calc_log10poisspdf2(vObsFMD',vPredFMD');
        % Sum the probabilities
        fProbability = (-1) * sum(vProb_);
        vfProbability = [vfProbability; fProbability];
        vMshift = [vMshift; fMshift];
        vMc = [vMc; fMc];
    catch
        vfProbability = [vfProbability; NaN];
        vMshift = [vMshift; NaN];
        vMc = [vMc; NaN];
    end; % of try-catch
end

try
    %%% Find the minimum loglikelihood score: if the minimum score is obtained several times, calculate MEAN
    %%% of the magnitude shift
    vdMloglikeli = [vfProbability vMshift];
    vSel = (vdMloglikeli == nan(vdMloglikeli(:,1)));
    vdMloglikeli = vdMloglikeli(vSel,:);
    if length(vdMloglikeli(:,1)) > 1
        fProbability = min(vdMloglikeli(:,1));
        fMshift = roundn(mean(vdMloglikeli(:,2)),-1);
    else
        fProbability = vdMloglikeli(:,1);
        fMshift = vdMloglikeli(:,2);
    end
    % KS-Test
    [bH,fPval,fKsstat] = kstest2(roundn(mCat2(:,6),-1),roundn(mCat1(:,6)+fMshift,-1),0.05,0);
    % Solution matrix
    mProblikelihood = [vMc vMshift vfProbability];
    nDegFree = 1; % Magnitude shift is the degree of freedom
    n_samples = length(mCat1(:,6))+length(mCat2(:,6));
    %% Corrected Akaike Information Criterion (AICc)
    fAICc = -2*(-fProbability)+2*nDegFree+2*nDegFree*(nDegFree+1)/(n_samples-nDegFree-1);

    % KSTest on Model
    vMag1 = [];
    mCat1Mod(:,6) = mCat1(:,6)+fMshift;
    [mResult, fMls, fMc, fMu, fSigma, mDataPred, vPredBest, fBvalue] = calc_McCdfnormal(mCat1Mod, fBinning);
    mPredFMD = mDataPred;
    %mPredFMD(:,1) = ceil(mDataPred(:,1)./fPeriod1);
    vSel = (mPredFMD(:,2) ~= 0); % Remove bins with zero frequency of zero events
    mData2 = mPredFMD(vSel,:);
    for nCnt=1:length(mData2(:,1))
        fM = repmat(mData2(nCnt,2),mData2(nCnt,1),1);
        vMag1 = [vMag1; fM];
    end
    [bH2,fPval2,fKsstat2] = kstest2(roundn(mCat2(:,6),-1),roundn(vMag1,-1),0.05,0);

catch
    mProblikelihood = [NaN NaN NaN];
    fAICc = NaN;
    fMshift = NaN;
    fProbability = NaN;
    bH = NaN;
    bH2 = NaN;
end; % of try-catch

author: CelsoReyes


text:
%% 2D Advection-Diffusion D2Q4
% by Manuel Diaz
clc; clear; %close all;

%% Domain Grid
L = 100;    % Length
H = 100;    % Height
n = 100;    % Nodes in x
m = 100;    % Nodes in y
dx = L/n; dy = H/m; dt = 1;
x = 1:dx:L; y = 1:dy:H;

%% Initial Condition
f1   = zeros(m,n);
f2   = zeros(m,n);
f3   = zeros(m,n);
f4   = zeros(m,n);
feq1 = zeros(m,n);
feq2 = zeros(m,n);
feq3 = zeros(m,n);
feq4 = zeros(m,n);
rho  = zeros(m,n); % initial value of the dependent variable rho = T(x,t)

%% Constants
u     = 0.1;
v     = -0.2; % <- interesting result if we change to -v!
ck    = dx/dt;
csq   = (dx^2)/(dt^2);
alpha = 0.25;
omega = 1/(2*alpha/(dt*csq)+0.5);
w     = [1/4,1/4,1/4,1/4];
cx    = [  1, -1,  0,  0];
cy    = [  0,  0,  1, -1];
link  = [  1,  2,  3,  4];
tEnd  = 400; % time steps

%% Movie Parameters
tPlot = 10; frames= tEnd/tPlot; count = 0;
figure(1)
colordef white %black

%% Boundary Conditions
twall = 1;

%% Main Loop
for cycle = 2:tEnd;
    % Collision process
    rho = f1 + f2 + f3 + f4;

    feq1 = 0.25 * rho *(1+2*u/ck);
    feq2 = 0.25 * rho *(1-2*u/ck);
    feq3 = 0.25 * rho *(1+2*v/ck);
    feq4 = 0.25 * rho *(1-2*v/ck);

    f1 = w(1) * rho;
    f2 = w(2) * rho;
    f3 = w(3) * rho;
    f4 = w(4) * rho;

    f1 = omega * feq1 + (1-omega) * f1;
    f2 = omega * feq2 + (1-omega) * f2;
    f3 = omega * feq3 + (1-omega) * f3;
    f4 = omega * feq4 + (1-omega) * f4;

    % Streaming process
    f1 = stream2d( f1 , [cy(1),cx(1)]);
    f2 = stream2d( f2 , [cy(2),cx(2)]);
    f3 = stream2d( f3 , [cy(3),cx(3)]);
    f4 = stream2d( f4 , [cy(4),cx(4)]);

    % Boundary conditions
    f1(:,1) = twall*w(1) + twall*w(2) - f2(:,1); %Dirichlet BC
    f3(:,1) = twall*w(3) + twall*w(4) - f4(:,1); %Dirichlet BC

    f1(m,:) = 0;        %Dirichlet BC
    f2(m,:) = 0;        %Dirichlet BC
    f3(m,:) = 0;        %Dirichlet BC
    f4(m,:) = 0;        %Dirichlet BC

    f1(:,n) = 0;        %Dirichlet BC
    f2(:,n) = 0;        %Dirichlet BC
    f3(:,n) = 0;        %Dirichlet BC
    f4(:,n) = 0;        %Dirichlet BC

    f1(1,:) = f1(2,:);  %Neumann BC
    f2(1,:) = f2(2,:);  %Neumann BC
    f3(1,:) = f3(2,:);  %Neumann BC
    f4(1,:) = f4(2,:);  %Neumann BC

    % Visualization
    if mod(cycle,tPlot) == 1
        count = count + 1;
        contourf(rho)
        colormap hot
        colorbar('location','southoutside')
        axis tight
        M(count) = getframe;
    end

    % Animated gif file
    if mod(cycle,tPlot) == 1
        F = getframe;
        if count == 1
            [im,map] = rgb2ind(F.cdata,256,'nodither');
            im(1,1,1,tEnd/tPlot) = 0;
        end
        im(:,:,1,count) = rgb2ind(F.cdata,map,'nodither');
    end
end

%% Make Movie
movie(M,1,10); % movie(M,n,fps)

%% Export to Gif
imwrite(im,map,'advec_diff_d2q4.gif','DelayTime',0,'LoopCount',3)

repo: Aero-matlab


text:
function bicycle_lock ( )

%*****************************************************************************80
%
%% BICYCLE_LOCK finds the combination on a typical bicycle lock.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
  timestamp ( )
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BICYCLE_LOCK\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A bicycle combination lock consists of 3 dials,\n' );
  fprintf ( 1, '  each having 10 symbols, 0 through 9.\n' );
  fprintf ( 1, '  We seek to determine the combination C.\n' );
%
%  Set the combination randomly.
%  We can think of the combination as a number between 0 and 999.
%
  rng ( 'shuffle' );
  c = randi ( [ 0, 999 ], 1, 1 );
%
%  To find the combination, simply generate every number between 0 and 999,
%  and then try it.
%
  for a = 0 : 999

    if ( a == c )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  The combination is %d!\n', c );
      break
    end

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BICYCLE_LOCK\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
%% SQUARESTOKE Stokes equations on the unit square
%
%   SQUARESTOKE computes P2-P1 approximations of the Stokes equations in
%   the unit square on a sequence of meshes obtained by uniform refinement.
%   It plots the approximation error (pressue in L2 norm and velocity in H1
%   norm) vs the number of dof. Other types of FEM approximation can be
%   computed similarly.
%
% See also StokesP2P1, collidingflow, squarePoisson
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all;
clear variables;
%% Set up
maxIt = 4;
N = zeros(maxIt,1);
h = zeros(maxIt,1);
erru = zeros(maxIt,1);
errp = zeros(maxIt,1);

%% Generate initial mesh
[node,elem] = squaremesh([0 1 0 1], 0.25);
bdFlag = setboundary(node,elem,'Dirichlet');

%% PDE and options
pde = Stokesdata2;
% pde = Stokesdata3;

%% Finite Element Method
for k = 1:maxIt
    % refine mesh
    [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
    % solve the equation
    [soln,eqn] = StokesP2P1(node,elem,bdFlag,pde);
%     [soln,eqn] = StokesP2P0(node,elem,bdFlag,pde);
%     [soln,eqn] = StokesisoP2P1(node,elem,bdFlag,pde);
%     [soln,eqn] = StokesisoP2P0(node,elem,bdFlag,pde);
%     [soln,eqn] = StokesCRP0(node,elem,bdFlag,pde);
%     [soln,eqn] = StokesP1bP1(node,elem,bdFlag,pde);
    uh = soln.u;
    ph = soln.p;
    N(k) = length(uh)+length(ph);
    h(k) = 1./(sqrt(size(node,1))-1);
    if N(k) < 2e3 % show solution for small size
        figure(1);  showresult(node,elem,ph);
    end
    % compute error
    uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
    erru(k) = sqrt((uh-uI(:))'*eqn.A*(uh-uI(:)));
    errp(k) = getL2error(node,elem,pde.exactp,ph);
end

%% Plot convergence rates
figure(2);
showrateh2(h,erru,1,'-*','|u_I-u_h|_1',...
           h,errp,1,'m-+','|| p-p_h||');

fprintf('\n');
disp('Table: Error')
colname = {'#Dof','h','|u_I-u_h|_1','||p-p_h||'};
disptable(colname,N,[],h,'%0.3e',erru,'%0.5e',errp,'%0.5e');
% figure;
% set(gcf,'Units','normal');
% set(gcf,'Position',[0.25,0.25,0.55,0.4]);

repo: ifem


text:
function [ rectangles ] = pruneRectangle( rotImg, rectangles )
%PRUNERECTANGLE Prune rectangles based on gradients
%   Rectangles should present stronger gradients on edges

rotImg = imresize(rotImg, 0.25);
[H,W,~] = size(rotImg);
% img = zeros(H,W);

[FX, FY] = gradient(double(rgb2gray(rotImg)));
F = sqrt(FX.^2 + FY.^2);
[FD, ~, ~] = dt2(F, 0.1, 0, 0.1, 0 );

for vid = 1:6
    numrect = size(rectangles(vid).highPixelBox,1);
    valid = false(numrect, 1);
    rects = rectangles(vid).xyzBox;
    for i = 1:numrect
        rect = [rects(i,1:3); rects(i,4:6); rects(i,7:9); rects(i,10:12)];
        lines = lineFromTwoPoint(rect([1 2 3 4],:), rect([2 3 4 1],:));
        panoReport = paintParameterLine( lines, W, H);
        G = FD(panoReport(:)>0);
        SG = sort(G, 'descend');
        N = round(length(G)/2);

%         valid(i) = sum(SG(1:N))/N;
        if sum(SG(1:N))/N > 20
            valid(i) = true;
        end
    end

    rectangles(vid).highPixelBox = rectangles(vid).highPixelBox(valid,:);
    rectangles(vid).count = sum(valid);
    rectangles(vid).xyzBox = rectangles(vid).xyzBox(valid,:);
    rectangles(vid).score = rectangles(vid).score(valid);
end

end

author: zouchuhang

repo: LayoutNet


text:
function R = RPYtoRot_ZXY(phi,theta,psi)
%RPYtoRot_ZXY Converts roll, pitch, yaw to a body-to-world Rotation matrix
%   The rotation matrix in this function is world to body [bRw] you will
%   need to transpose this matrix to get the body to world [wRb] such that
%   [wP] = [wRb] * [bP], where [bP] is a point in the body frame and [wP]
%   is a point in the world frame
%   written by Daniel Mellinger
%

R = [cos(psi)*cos(theta) - sin(phi)*sin(psi)*sin(theta), ...
     cos(theta)*sin(psi) + cos(psi)*sin(phi)*sin(theta), ...
     -cos(phi)*sin(theta); ...
     -cos(phi)*sin(psi),...
     cos(phi)*cos(psi), ...
     sin(phi);...
     cos(psi)*sin(theta) + cos(theta)*sin(phi)*sin(psi),...
     sin(psi)*sin(theta) - cos(psi)*cos(theta)*sin(phi),...
     cos(phi)*cos(theta)];

end

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
classdef nnl2norm < nntest
  methods (Test)

    function basic(test)
      x = test.randn([5 5 3 3]) ;
      y = vl_nnl2norm(x) ;

      % check derivatives with numerical approximation
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nnl2norm(x, dzdy) ;
      test.der(@(x) vl_nnl2norm(x), x, dzdy, dzdx, 1e-3*test.range) ;
    end

  end
end

repo: MFT


text:
function correct_rate = test_speaker_identification
    % Test speaker identification with TIDIGITS
    USE_CACHE = false;

    USE_GMM = true;
    USE_BINARYSVM = false;
    USE_MULTISVM = false;
    USE_DTW = false;

    db = gendb('tidigits_gyro');
    db = filterdb(db, 'age', '1-Adults');

    if USE_CACHE
        load('speaker_identification_features', 'features', 'labels');
    else
        [features, labels] = ...
            get_speaker_features_and_labels(db, 'speaker');
        savesave('speaker_identification_features', 'features', 'labels');
    end

    db = filterdb(db, 'type', 'WOMAN');
%     db = filterdb(db, 'device', '0094e779d7d1841f');
    db = filterdb(db, 'device', '00a697fa469633a5');
    speakers = get_speakers(db);

    % choose random speakers
    SPEAKERS_NUM = 5;
    TOTAL_NUM = length(speakers);
    speaker_ids = randperm(TOTAL_NUM);
    speaker_ids = sort(speaker_ids(1:SPEAKERS_NUM));
    % set mask
    mask = false(size(labels));
    for i = 1:length(speaker_ids)
       mask(labels == speaker_ids(i)) = true;
    end
    % filter speakers
    features = features(:, mask);
    labels = labels(mask);

    % Choose train and test sets
    cp = cvpartition(labels, 'holdout', 0.1);
    train_labels = labels(cp.training);
    test_labels = labels(cp.test);

    if USE_GMM || USE_MULTISVM || USE_BINARYSVM
        train_features = cell2mat(features(cp.training));
        test_features = cell2mat(features(cp.test));
    end

    if USE_GMM
        NUM_OF_GAUSSIANS = 8;
        NUM_OF_ITERATIONS = 20;
        [mu_train, sigma_train, c_train] = ...
            GMM.gmm_training(train_features, train_labels, NUM_OF_GAUSSIANS, ...
                             NUM_OF_ITERATIONS);
        class = GMM.gmm_classification(test_features, mu_train, sigma_train, c_train);
    end

    if USE_MULTISVM
        class = multisvm(train_features', train_labels, test_features', ...
            'tolkkt', 1e-2, 'kktviolationlevel', 0.1);
    end

    if USE_BINARYSVM
        s = svmtrain(train_features', train_labels, ...
            'tolkkt', 1e-2, 'kktviolationlevel', 0.1);
        class = svmclassify(s, test_features');
    end

    if USE_DTW
        train_features = features(cp.training);
        test_features = features(cp.test);

        class = zeros(size(test_labels));
        display('Classifying...');
        progressbar;
        for i = 1:length(test_labels)
            sample = test_features{i};
            class(i) = dtw_classify_sample(sample, train_features, train_labels);
            progressbar(i/length(test_labels));
        end
    end

    if ~USE_DTW
        u = unique(train_labels);
        class = u(class);
    end

    correct = class == test_labels;
    correct_rate = sum(correct)/length(correct);
    display(correct_rate);
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function [xn, yn, zn] = findnearest (x, y, z, x0, y0, N)
d = (x0 - x).^2 + (y0 - y).^2;
[md, i] = min (d);
xn1 = x(i); yn1 = y(i); zn1 = z(i);
xn2 = []; yn2 = []; zn2 = [];
if N > 1
    x = [x(1:i-1); x(i+1:end)];
    y = [y(1:i-1); y(i+1:end)];
    z = [z(1:i-1); z(i+1:end)];
    [xn2, yn2, zn2] = findnearest (x, y, z, x0, y0, N-1);
end
xn = [xn1; xn2];
yn = [yn1; yn2];
zn = [zn1; zn2];


text:
function v = rotate(v,q,varargin)
% rotate vector3d by rotation or orientation
%
% Syntax
%   v = rotate(v,20*degree) % rotation about the z-axis
%   rot = rotation.byEuler(10*degree,20*degree,30*degree)
%   v = rotate(v,rot)
%
% Description
%  Either |v| or |rot| are single elements or both have the same size. The
%  ouptut |v| will have the same size as the biger of both input arrays.
%
% Input
%  v - @vector3d
%  q - @quaternion
%
% Output
%  r - q * v
%

if isnumeric(q), q = axis2quat(zvector,q);end

if ~isa(q,'rotation')
  [a,b,c,d] = double(q);
  i = [];
else
  [a,b,c,d,i] = double(q);
end
[x,y,z] = double(v);

n = b.^2 + c.^2 + d.^2;
s = 2*(x.*b + y.*c + z.*d);

a_2 = 2*a;
a_n  = a.^2 - n;

v.x = a_2.*(c.* z - y.*d) + s.*b + a_n.*x;
v.y = a_2.*(d.* x - z.*b) + s.*c + a_n.*y;
v.z = a_2.*(b.* y - x.*c) + s.*d + a_n.*z;

if ~isempty(i)
  if numel(i)>1
    i = logical(i);
    v.x(i) = -v.x(i);
    v.y(i) = -v.y(i);
    v.z(i) = -v.z(i);
  elseif i
    v.x = -v.x;
    v.y = -v.y;
    v.z = -v.z;
  end
end

if isa(q,'orientation')

  if isa(q.SS,'crystalSymmetry')
    v = Miller(v,q.SS);
    v.dispStyle = MillerConvention(v.dispStyle);
    v.dispStyle = make4Digit(v.dispStyle,q.SS);
  else % convert to vector3d
    v = vector3d(v);
  end

end


text:
%tests checkStoichiometricConsistency and minCardinalityConservationRelaxationVector
if ~exist('movel','var')
    load 'Recon3DModel_301.mat'
end
printLevel=0;

model = findSExRxnInd(model,[],printLevel-1);

N = model.S(:,model.SIntRxnBool);

%Recon3DModel_301 is stoichiometrically consistent, so check should be
%positive
[isConsistent, m, model] = checkStoichiometricConsistency(model, printLevel);
assert(isConsistent==1)

[mlt,nlt]=size(N');
feasTol = getCobraSolverParams('LP', 'feasTol');
param.eta=feasTol*100;
param.checkConsistency=0;
param.epsilon=1e-4;
param.nonRelaxBool=false(mlt,1);
param.checkFeasibility = 0;
param.printLevel=printLevel;

%Recon3DModel_301 is stoichiometrically consistent, so no relaxations
%should be needed
[relaxRxnBool, solutionRelax] = minCardinalityConservationRelaxationVector(N,param);
assert(nnz(relaxRxnBool)==0);

[relaxRxnBool, solutionRelax] = minCardinalityConservationRelaxationVector(N);
assert(nnz(relaxRxnBool)==0);

if solutionRelax.stat==1
    %conserved if relaxation is below epsilon
    relaxRxnBool=abs(solutionRelax.x)>=param.eta;
    if printLevel>1
        fprintf('%g%s\n',norm(N(:,~relaxRxnBool)'*solutionRelax.z),' = ||N''*z|| (should be zero)')
    end
    if printLevel>1
        fprintf('%s\n',[int2str(nnz(relaxRxnBool)) '/' int2str(length(relaxRxnBool)) ' reactions relaxed.'])
    end
else
    disp(solutionRelax)
    error('solve for minimum cardinality of conservation relaxation vector failed')
end

assert(isConsistent & nnz(relaxRxnBool)==0);

author: opencobra

repo: cobratoolbox


text: freeFemSettings = SettingsFreeFemMeshGenerator();

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [err,p,f] = GetPath (s, allmc)
%
%   [err,PathString,FileString] = GetPath (s, allmc)
%
%Purpose:
%   Breaks the string s into a path and file part
%
%
%Input Parameters:
%   s is a string like How/Didley/Doo
%   allmc : 1 --> search for path using both types of file-separators / and \
%           0 --> search for path using filesep's output (default)
%
%Output Parameters:
%
%if nargout = 1
%  returns PathString
%
%Otherwise:
%   err : 0 No Problem
%       : 1 Mucho Problems
%
%   PathString is a string like How/Didley/
%   FileString is a string like Doo
%
%More Info :
%
%
%
%
%     Author : Ziad Saad
%     Date : Thu Apr 23 11:41:06 CDT 1998

%Define the function name for easy referencing
FuncName = 'GetPath';

if (nargin == 1),
   allmc = 0;
end

%initailize return variables
err = 1;
p = '.';
f = '';

N = length(s);
if (N == 0),
	err = ErrEval(FuncName,'Err_Emtpy string');
   if (nargout == 1) err = ''; end
   return;
end

if (allmc == 0),
   [i] = findstr(s,filesep);
else
   [i] = union(findstr(s,'/'), findstr(s,'\'));
end
if (isempty(i)),
	p = ['.' filesep];
	f = s;
elseif (max(i)==N);
	p = s;
	f = '';
else
	p = s(1:max(i));
	f = s(max(i)+1:N);
end

if (nargout < 2) err = p; else err = 0; end

return;


text:
function output = memtic(action, counter)

% MEMTIC start a MATLAB memory recorder
%
%   MEMTIC and MEMTOC functions work together to measure memory usage.
%   MEMTIC, by itself, saves the current memory footprint that MEMTOC
%   uses later to measure the memory that was used between the two.
%
% Use as
%   MEMTIC
%   MEMTOC
% to print the estimated memory use on screen, or
%   MEMTIC
%   M = MEMTOC
% to return the estimated memory (in bytes) in variable M, or
%   C = MEMTIC
%   M = MEMTOC(C)
% to specifically estimate the memory use between a well-defined tic/toc pair.
%
% Note that MATLAB uses internal memory allocation, garbage collection, shallow
% copies of variables, and virtual memory. Due to the advanced handling of
% memory for its variables, it is not easy and in certain cases not possible to
% make a reliable and reproducible estimate based on the memory information
% provided by the operating system.
%
% Example: measure the memory increase due to allocating a lot of memory.
% Doing a "clear x" following the allocation and priot to MEMTOC does not
% affect the memory that is reported.
%
%   memtic
%   n = 125; x = cell(1,n);
%   for i=1:n
%     x{i} = randn(1000,1000); % 8kB per item
%     disp(i);
%   end
%   whos x
%   memtoc
%
% See also TIC, TOC

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

persistent state

if nargin<1
  action = 'tic';
end

% the memtic/memtoc functions make use of a low-level mex file that interacts directly with the operating system
% do not fail if the mex file does not exist
if isempty(strfind(which('memprofile'), mexext))
  switch action
    case 'tic'
      if nargout
        output = nan;
      end
    case 'toc'
      if nargout
        output = nan;
      end
    otherwise
      ft_error('invalid input argument #1');
  end % switch
  return
end % if mex file does not exist

% get the current usage, this will start memprofile if needed
memstat = memprofile('info');

switch action
  case 'tic'
    if nargin>1
      ft_error('the counter cannot be specified as imput argument');
    end

    counter = length(state)+1;
    state(counter).running = true;
    state(counter).time    = memstat(end).time;
    state(counter).mem     = memstat(end).mem;

    if nargout
      output = counter;
    end

  case 'toc'
    if nargin<2
      % take the latest
      counter = length(state);
    elseif counter<1 || counter>numel(state)
      ft_error('invalid counter');
    end

    if counter==0
      memused = 0;
    else
      sel = [memstat.time]>=state(counter).time;
      memused = max([memstat(sel).mem]);       % determine the maximum
      memused = memused - state(counter).mem;  % subtract the initial value
    end

    if nargout
      output = memused;
    end

  otherwise
    ft_error('invalid input argument #1');
end % switch


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%  	System response to sinusoidal inputs

t=0:.1:30;
x1=3*cos(t+pi/3);
plot(t,x1);
legend('input signal x(t)')
ylim([-4 4]);

figure
w0=1;
Hw0=(3*(j*w0)^2+4j*w0+2)/(-w0^2+j*w0+3)
magn=abs(Hw0)
phas=angle(Hw0)
y1=3*abs(Hw0)*cos(t+pi/3+angle(Hw0));
plot(t,y1);
legend('system response')
ylim([-7 7]);

figure
num=[3 4 2];
den=[1 1 3];
yls=lsim(num,den,x1,t);
plot(t,yls);
legend('system response by lsim')

figure
plot(t,x1,t,y1,'o')
legend('x(t)','y(t)')
ylim([-7 9]);


text:
function merged_sample = merge_samples(sample1, sample2, w1, w2, sample_merge_type)
% Merge sample1 and sample2 using weights w1 and w2. The type of merging is
% decided by sample_merge_type.
% The sample_merge_type can be
% 1) Merge: The output is the weighted sum of the input samples
% 2) Replace: The output is the first sample. ie w2 is assumed to be 0

% Normalise the weights so that they sum to one
alpha1 = w1/(w1+w2);
alpha2 = 1 - alpha1;

% Build the merged sample
if strcmpi(sample_merge_type, 'replace')
    merged_sample = sample1;
elseif strcmpi(sample_merge_type, 'merge')
    num_feature_blocks = numel(sample1);

    merged_sample = cell(1, 1, num_feature_blocks);
    for k = 1:num_feature_blocks
        merged_sample{k} = alpha1*sample1{k} + alpha2*sample2{k};
    end
else
    error('Invalid sample merge type');
end

author: martin-danelljan

repo: ECO


text:
function check13Color(img)
%
%
%        check13Color(img)
%
%
%        Input:
%           -img: an image to be tested if it is an one or three color
%                 channels image
%
%     Copyright (C) 2013  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

col = size(img, 3);
if(col ~= 3 && col ~= 1)
    error('The image has to be an RGB or luminance image.');
end

author: banterle


text:
function [ ns, xyz ] = sphere_cubed_points_face ( n, i1, j1, k1, i2, ...
  j2, k2, ns, xyz )

%*****************************************************************************80
%
%% SPHERE_CUBED_POINTS_FACE: points on one face of a cubed sphere grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of sections into which each face of
%    the cube is to be divided.
%
%    Input, integer I1, J1, K1, I2, J2, K2, the logical indices, between 0
%    and N, of two corners of the face grid.  It is guaranteed that I1 <= I2,
%    J1 <= J2, and K1 <= K2.
%
%    Input, integer NS, the number of points.
%
%    Input, real XYZ(3,NS), distinct points on the unit sphere
%    generated by a cubed sphere grid.
%
%    Output, integer NS, the number of points.
%
%    Output, real XYZ(3,NS), distinct points on the unit sphere
%    generated by a cubed sphere grid.
%
  for i = i1 : i2

    if ( i1 < i2 )
      xc = tan ( ( 2 * i - n ) * 0.25 * pi / n );
    elseif ( i1 == 0 )
      xc = -1.0;
    elseif ( i1 == n )
      xc = +1.0;
    else
      xc = 0.0;
    end

    for j = j1 : j2

      if ( j1 < j2 )
        yc = tan ( ( 2 * j - n ) * 0.25 * pi / n );
      elseif ( j1 == 0 )
        yc = -1.0;
      elseif ( j1 == n )
        yc = +1.0;
      else
        yc = 0.0;
      end

      for k = k1 : k2

        if ( k1 < k2 )
          zc = tan ( ( 2 * k - n ) * 0.25 * pi / n );
        elseif ( k1 == 0 )
          zc = -1.0;
        elseif ( k1 == n )
          zc = +1.0;
        else
          zc = 0.0;
        end

        xyzn = sqrt ( xc^2 + yc^2 + zc^2 );

        ns = ns + 1;
        xyz(1,ns) = xc / xyzn;
        xyz(2,ns) = yc / xyzn;
        xyz(3,ns) = zc / xyzn;

      end
    end
  end

  return
end

repo: jburkardt-m


text:
function volume = sphere_volume_nd ( dim_num, r )

%*****************************************************************************80
%
%% SPHERE_VOLUME_ND computes the volume of an implicit sphere in ND.
%
%  Discussion:
%
%    An implicit sphere in ND satisfies the equation:
%
%      sum ( ( X(1:N) - CENTER(1:N) )**2 ) = R**2
%
%    where R is the radius and CENTER is the center.
%
%    Results for the first few values of N are:
%
%    DIM_NUM  Volume
%    -     -----------------------
%    2                PI    * R**2
%    3     (4/3)    * PI    * R**3
%    4     (1/2)    * PI**2 * R**4
%    5     (8/15)   * PI**2 * R**5
%    6     (1/6)    * PI**3 * R**6
%    7     (16/105) * PI**3 * R**7
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the space.
%
%    Input, real R, the radius of the sphere.
%
%    Output, real VOLUME, the volume of the sphere.
%
  volume = r^dim_num * sphere_unit_volume_nd ( dim_num );

  return
end


%% SPHP-garden
% %{
imfolder = 'images\SPHP-garden';
im_n = 2;
imfile = cell(im_n,1);
imfile{1} = [imfolder '\' 'garden_01.jpg'];
imfile{2} = [imfolder '\' 'garden_02.jpg'];

im = cell(im_n,1);
for ii = 1:im_n
    im{ii} = imread(imfile{ii});
end

edge_list = [1,2];

imsize = zeros(im_n,3);

for ii = 1:im_n
    imsize(ii,:) = size(im{ii});
    if imsize(ii,1) > 720
        scale = 720/size(im{ii}, 1);
        im{ii} = imresize(im{ii}, scale);

        imsize(ii,:) = size(im{ii});
    end
end

mosaic = REW_mosaic( im, edge_list, 1, [], 0, imfolder );
%}


text:
function mu = circularmean(points, weights, angleidx)

% mu = circularmean(points, weights, angleidx)

weights = weights(:)'/sum(weights(:));
mu = sum(bsxfun(@times, points, weights),2);
for i=angleidx(:)'
    sinx = sum(sin(points(i,:)).*weights);
    cosx = sum(cos(points(i,:)).*weights);
    mu(i) = atan2(sinx, cosx);
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
 function F_enhance= RMSHEg_gray( F )
% clear all;
% close all;
% clc;
x0=F;
%x0=imread('om.bmp');
%x0=imread('11.bmp');
 %x0=imread('12.bmp');
%y0=rgb2ycbcr(x0);
luma=x0;%(:,:,1);
[m,n]=size(luma);
% cb=y0(:,:,2);
% cr=y0(:,:,3);
k=mean(mean(luma));
XM=k;
r=round(k);
l=length(luma(luma<=r));
         listindex=find(luma<=r);
 k1=reshape(luma(listindex),l,1);%obtained below fist mean values put it into an arry
 k2=mean(k1);%find mean of mean in first half
 k3=round(k2);%round to nearest integer to find pdf and cdf till first quarter
  sum=0;
 l1=length(luma(luma<=k3));%find all  values in luma falls below fist mean
 listindex1=find(luma<=k3);%all  values in luma falls below fist mean given in their locations
 k4=reshape(luma(listindex1),l1,1);%arrange all fist mean values in a vector
     xpdf=hist(k4,[0:k3]);%pdf from 0:r
      xpdf=xpdf/l1;%normalized pdf to get nk/n,l=sum of xpdf,total no of pixels form 0 to r.
     % plot(xpdf);
%      plot(imhist(luma));
%       xlabel('gray levels up to mean');
%       ylabel('pdf up to mean');
%       title('histogram for half an image up to mean');
      sk=xpdf*triu(ones(k3+1));
%       figure(2);
%       plot(sk);
%       xlabel('gray levels upto 1st mean');
%       ylabel('cdf upto mean');
%       title('cdf for half  of an image up to ist mean');
      alpha=0.6;
      for l2=0:k3
              list1=find(k4==l2);%find value in an vector i.e converted from matrix
              %list(list1)=alpha*sk(l2+1)*(k3+1)+(1-alpha)*(k3+1);
              list(list1)=sk(l2+1)*(k3+1);%map dont disturb to get bhe as
              %it is 13/3/2011
              ert(l2+1)=sk(l2+1)*(k3+1);
      end
          p=zeros(m,n);
      p(listindex1)= list;
%      figure(3);
%      imshow(p);
%      xlabel('gray levels up to first mean');
%       ylabel('luma component equilized image up to first mean');
%      title('processed luma image up to first mean');
     k=mean(mean(luma));
r=round(k);
l=length(luma(luma<=r));
        listindex=find(luma<=r);
 k1=reshape(luma(listindex),l,1);%obtained below fist mean values put it into an arry
 k2=mean(k1);%find mean of mean in first half
 k3=round(k2);%round to nearest integer to find pdf and cdf till first quarter
%   sum=0;
   b=k3;
%  l2=length(luma(luma>k3));%find all  values in luma falls below fist mean
  listindex2=find((luma>k3)&(luma<=r));%all  values in luma falls below fist mean given in their locations
  l2=length(listindex2);
  k5=reshape(luma(listindex2),l2,1);%arrange all fist mean values in a vector
      x2pdf=hist(k5,[k3+1:r]);%pdf from 0:r
       x2pdf=x2pdf/l2;%normalized pdf to get nk/n,l=sum of xpdf,total no of pixels form 0 to r.
%        figure(4);
%        plot(x2pdf);
%        xlabel('gray levels 2nd mean');
%        ylabel('pdf of 2nd mean');
%        title('histogram for 2nd mean');
       sk2=x2pdf*triu(ones(r-k3));
%        figure(5);
%        plot(sk2);
%        xlabel('gray levels of 2nd mean');
%        ylabel('cdf upto mean');
%        title('cdf for half  of an image of 2nd mean');
       k2u=1;
       for l3=k3+1:r
               list2=find(k5==l3);%find value in an vector i.e converted from matrix
               %list1(list2)=alpha*(sk2(k2u)*(r-k3)+(k3+1))+(1-alpha)*(k3+1);
               list1(list2)=(k3+1)+(sk2(k2u))*(r-k3);%map dont disturb to
               %get BHE 13/3/2011
               ert(l3)=(k3+1)+(sk2(k2u))*(r-k3);
               k2u=k2u+1;
       end
           p1=zeros(m,n);
       p1(listindex2)= list1;
%       figure(6);
%       imshow(p1);
%      xlabel('gray levels up to first 2nd mean');
%       ylabel('luma component equilized image 2nd mean');
%      title('processed luma image up to 2nd mean');
 %lupper30=length(luma(luma>r);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 lupper30=length(luma(luma>r));
%for i=0:r
    %if(luma(luma<=r))
        listindexupper3=find(luma>r);
   % end
%end
 k1upper=reshape(luma(listindexupper3),lupper30,1);
 mean3=mean(k1upper);
  r3=round(mean3);
 %length30=length((luma<=r3)&(luma>r));
  listindexupper30=find((luma>r)&(luma<=r3));
  length30=length(listindexupper30);
  k30upper=reshape(luma(listindexupper30),length30,1);

  %length30=length((luma<=r3)&(luma>r));
%  sum=0;
      xpdfupper30=hist(k30upper,[r+1:r3]);%pdf from r+1:r3
      xpdfupper30=xpdfupper30/length30;%normalized pdf to get nk/n,l=sum of xpdf,total no of pixels form r+1 to 255.
%       figure(7);
%      plot(xpdfupper30);
%      xlabel('gray levels   3rd mean');
%      ylabel('pdf of 3rd mean');
%      title('histogram for upper half an image  3rd mean');
      skupper30=xpdfupper30*triu(ones(r3-r));
%       figure(8);
%       plot(skupper30);
%       xlabel('gray levels after mean');
%       ylabel('cdf after mean');
%       title('cdf for upper half  of an image after mean');
      k3u=1;
      for k3upper=(r+1):r3
              list1upper30=find(k30upper==k3upper);%find value in an vector i.e converted from matrix
              listnew(list1upper30)=(r+1)+skupper30(k3u)*(r3-r);%map dont
              %disturg to get original heq 14/3/2011
              %listnew(list1upper30)=alpha*(skupper30(k3u)*(r3-r)+(r+1))+(1-alpha)*(r+1);
              ert(k3upper)=(r+1)+skupper30(k3u)*(r3-r);
         k3u=k3u+1;
     end

          p2=zeros(m,n);

     p2(listindexupper30)= listnew;
%     figure(9);
%      imshow(p2);
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     listindexupper40=find(luma>r3);
     lupper40=length(listindexupper40);
   % end
%end
 k1upper40=reshape(luma(listindexupper40),lupper40,1);
 %sum=0;
 %for i=0:r
     xpdfupper40=hist(k1upper40,[r3+1:255]);%pdf from r+1:255
     xpdfupper40=xpdfupper40/lupper40;%normalized pdf to get nk/n,l=sum of xpdf,total no of pixels form r+1 to 255.
%      figure(10);
%      plot(xpdfupper40);
%      xlabel('gray levels after 4mean');
%      ylabel('pdf after 4mean');
%      title('histogram for upper half an image after 4mean');
     skupper40=xpdfupper40*triu(ones(255-r3));
%      figure(11);
%      plot(skupper40);
%      xlabel('gray levels after 4mean');
%      ylabel('cdf after 4mean');
%      title('cdf for upper half  of an image after 4mean');
     k4u=1;
     for k4upper=(r3+1):255
         %if(xpdfupper(k2upper)>r)
             list1upper40=find(k1upper40==k4upper);%find value in an vector i.e converted from matrix
             %for k2u=1:58
             %listnew4(list1upper40)=alpha*(198+skupper40(k4u)*(255-r3))+(1-alpha)*(r3+1);
             listnew4(list1upper40)=(r3+1)+skupper40(k4u)*(255-r3);%map
             %dont disturb to get original bbhe 14/3/2011
         %end
         ert(k4upper)=(r3+1)+skupper40(k4u)*(255-r3);
         k4u=k4u+1;
     end

%      for i=0:l-1
          p3=zeros(m,n);
%          if (p(listindex))
%              p(:)=list;
%          end
%      end

     p3(listindexupper40)= listnew4;
%     figure(12);
%     imshow(p3);
%     xlabel('gray levels after 4mean');
%      ylabel('luma component equilized image after 4mean');
%      title('processed luma image after 4mean');
     om=p+p1+p2+p3;
      F_enhance=om;
 end
     %ommmmm=p1+p;
%     figure(13);
% %     imshow(ommmmm);
%      colormap('gray');
%       xlabel('gray level');
%      ylabel('combined lower and upper half luma component equilized image');
%      title('combined luma image');
%     figure(8);
%     imshow(uint8(om));
%     for j1=0:255
%         count=0;
%         for i1=0:m*n-1
%             if om(i1+1)==j1
%                 count=count+1;
%             end
%         end
%             prob(j1+1)=count/m*n;
%     end
%     figure(16);
%     plot(prob);
%
%    for j2=0:255
%         count1=0;
%         for i2=0:m*n-1
%             if luma(i2+1)==j2
%                 count1=count1+1;
%             end
%         end
%             prob2(j2+1)=count1/m*n;
%     end
%     figure(17);
%     plot(prob2);
%     xlabel('gray levels after mean');
%      ylabel('luma component equilized image after mean');
%      title('processed luma image after mean');
%     ommmmm=p1+p;
%     figure(7);
%      colormap('gray');
%       xlabel('gray level');
%      ylabel('combined lower and upper half luma component equilized image');
%      title('combined luma image');
%     image(ommmmm);
% % % %      cat1=cat(3,om,cb,cr);
% % % %     figure(14);
% % % %     imshow(cat1);
% % % %     xlabel('gray level(ycbcr)');
% % % %      ylabel('combined lower and upper half luma,cromablue,croma red component equilized image');
% % % %      title('luma croma b and r color processed image');
% % % %     catconversion=ycbcr2rgb(cat1);
% % % %     figure(15);
% % % %     imshow(catconversion);
%     xlabel('gray level(rgb)');
%      ylabel('combined lower and upper half  RGB component equilized image');
%      title('converted from ycbcr2rgb color(RGB) processed image');
%      YM=mean(mean(om));
%      AMBE=abs(XM-YM);
%      disp('Absolute mean brightness error=');
%      disp(AMBE);
%      MSE = MeanSquareError(luma, om);
%     % figure(16);
% disp('Mean Square Error = ');
% disp(MSE);
% PSNR = PeakSignaltoNoiseRatio(luma, om);
% %figure(17);
% disp('Peak Signal to Noise Ratio = ');
% disp(PSNR);
% E=entropy(uint8(om));
% disp('Entropy=');
% disp(E);
%      figure(16);
%      plot(imhist(om));

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function file_datenum = filedatenum(fname)
% <cpp> Purpose: Returns the modification date (datenum) of a file
% Usage: file_datenum = filedatenum('myfile.dat')
% Author: Charles Plum, Varolii Corp, 21 August 2009
%*********************************************************************

% Extract the file path from the argument fname
dir_name = fname(1:max(strfind (fname,'\')));

%If a path was given, then extract the file name from the fname arg
if (dir_name)
    file_name = fname(max(strfind (fname,'\'))+1:end);
else % assume the local directory is the path
    dir_name = pwd;
    file_name = fname;
end

% Not elegant but it works ...
d = dir(dir_name);
for ifile=1:length(d)
    if (strcmp(file_name,char(d(ifile).name)))
      file_datenum = d(ifile).datenum;
      return;
    end
end
file_datenum = 0; % Return zero if no file found
return

repo: mcbench-benchmarks


text:
function [cmap]=buildcmap(colors)
% [cmap]=buildcmap(colors)
%
% This function can be used to build your own custom colormaps. Imagine if
% you want to display rainfall distribution map. You want a colormap which
% ideally brings rainfall in mind, which is not achiveved by colormaps such
% as winter, cool or jet and such. A gradient of white to blue will do the
% task, but you might also use a more complex gradient (such as
% white+blue+red or colors='wbr'). This function can be use to build any
% colormap using main colors rgbcmyk. In image processing, w (white) can be
% used as the first color so that in the output, the background (usually
% with 0 values) appears white. In the example of rainfall map, 'wb' will
% produce a rainfall density map where the background (if its DN values are
% 0) will appear as white.
%
% Inputs:
%  colors: string (char) of color codes, any sequence of rgbcmywk
%  representing different colors (such as 'b' for blue) is acceptable. If a
%  gradient of white to blue is needed, colors would be 'wb'; a rainbow of
%  white+blue+red+green would be 'wbrg'.
%
% Example:
%  [cmap]=buildcmap('wygbr');
% %try the output cmap:
% im=imread('cameraman.tif');
% imshow(im), colorbar
% colormap(cmap) %will use the output colormap
%
% First version: 14 Feb. 2013
% sohrabinia.m@gmail.com
%--------------------------------------------------------------------------

if nargin<1
    colors='wrgbcmyk';
end

if ~ischar(colors)
    error(['Error! colors must be a variable of type char with '...
        'color-names, such as ''r'', ''g'', etc., '...
        'type ''help buildcmap'' for more info']);
end

ncolors=length(colors)-1;

bins=round(255/ncolors);
% diff1=255-bins*ncolors;

vec=zeros(300,3);

switch colors(1)
    case 'w'
        vec(1,:)=1;
    case 'r'
        vec(1,:)=[1 0 0];
    case 'g'
        vec(1,:)=[0 1 0];
    case 'b'
        vec(1,:)=[0 0 1];
    case 'c'
        vec(1,:)=[0 1 1];
    case 'm'
        vec(1,:)=[1 0 1];
    case 'y'
        vec(1,:)=[1 1 0];
    case 'k'
        vec(1,:)=[0 0 0];
end

for i=1:ncolors
 beG=(i-1)*bins+1;
 enD=i*bins+1; %beG,enD
 switch colors(i+1)
     case 'w'
         vec(beG:enD,1)=linspace(vec(beG,1),1,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),1,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),1,bins+1)';%colors(i+1),beG,enD,
     case 'r'
         vec(beG:enD,1)=linspace(vec(beG,1),1,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),0,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),0,bins+1)';%colors(i+1),beG,enD
     case 'g'
         vec(beG:enD,1)=linspace(vec(beG,1),0,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),1,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),0,bins+1)';%colors(i+1),beG,enD
     case 'b'
         vec(beG:enD,1)=linspace(vec(beG,1),0,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),0,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),1,bins+1)';%colors(i+1),beG,enD
     case 'c'
         vec(beG:enD,1)=linspace(vec(beG,1),0,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),1,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),1,bins+1)';%colors(i+1),beG,enD
     case 'm'
         vec(beG:enD,1)=linspace(vec(beG,1),1,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),0,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),1,bins+1)';
     case 'y'
         vec(beG:enD,1)=linspace(vec(beG,1),1,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),1,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),0,bins+1)';
     case 'k'
         vec(beG:enD,1)=linspace(vec(beG,1),0,bins+1)';
         vec(beG:enD,2)=linspace(vec(beG,2),0,bins+1)';
         vec(beG:enD,3)=linspace(vec(beG,3),0,bins+1)';
 end
end
cmap=vec(1:bins*ncolors,:);
end %end of buildcmap


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

% RBFINFWHITEKERNGRADX Gradient of RBF-WHITE kernel (with integration limits
% between minus infinity and infinity) 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 gradient 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
% and t2 (currently always one).
%
% SEEALSO rbfinfwhiteKernParamInit, kernGradX, rbfinfwhiteKernDiagGradX
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin < 3
  t2 = t1;
end
if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end

gT = zeros(size(t1, 1), 1, size(t2, 1));

% Parameters of the kernel required in the computation
variance = kern.variance;
invWidth = kern.inverseWidth;

for i = size(t1, 1)
    gT(i, 1, :) = - ( variance * invWidth * (t1(i)-t2) * sqrt(invWidth/pi) / 4) ...
        .* exp(- invWidth * ((t1(i)-t2).^2) / 4);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% [INPUT]
% pods = A vector of floats [0,1] of length n representing the probabilities of default.
% g = A boolean n^2-by-n matrix representing the posterior density orthants.
% p = A vector of floats [0,1] of length n^2 representing the posterior density probabilities.
%
% [OUTPUT]
% jpod = A float [0,1] representing the Joint Probability of Default.
% fsi = A float [1,n] representing the Financial Stability Index.
% pce = A float [0,1] representing the Probability of Cascade Effects.
% dide = A float n-by-n matrix [0,1] representing the Distress Dependency.
% si = A row vector of floats [0,1] of length n representing the Systemic Importance.
% sv = A row vector of floats [0,1] of length n representing the Systemic Vulnerability.
% cojpods = A row vector of floats [0,1] of length n representing the Conditional Joint Probabilities of Default.

function [jpod,fsi,pce,dide,si,sv,cojpods] = cross_entropy_metrics(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('pods',@(x)validateattributes(x,{'double'},{'real' 'finite' '>=' 0 '<=' 1 'vector' 'nonempty'}));
        ip.addRequired('g',@(x)validateattributes(x,{'double'},{'real' 'finite' 'binary' '2d' 'nonempty'}));
        ip.addRequired('p',@(x)validateattributes(x,{'double'},{'real' 'finite' '>=' 0 '<=' 1 'vector' 'nonempty'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    [pods,g,p] = validate_input(ipr.pods,ipr.g,ipr.p);

    nargoutchk(7,7);

    [jpod,fsi,pce,dide,si,sv,cojpods] = cross_entropy_metrics_internal(pods,g,p);

end

function [jpod,fsi,pce,dide,si,sv,cojpods] = cross_entropy_metrics_internal(pods,g,p)

    n = numel(pods);
    g_refs = sum(g,2);

    jpod = p(g_refs == n,:);
    fsi = min(max(sum(pods,'omitnan') / (1 - p(g_refs == 0,:)),1),n);
    pce = sum(p(g_refs >= 2,:)) / sum(p(g_refs >= 1,:));

    dide = eye(n);

    for i = 1:n
        for j = 1:n
            if (isnan(pods(j)))
                dide(i,j) = NaN;
            elseif (i ~= j)
                dide(i,j) = p((g_refs == 2) & (g(:,i) == 1) & (g(:,j) == 1),:) / pods(j);
            end
        end
    end

    dide_pods = ((dide - eye(n)) .* repmat(pods,1,n));
    si = sum(dide_pods,2);
    sv = sum(dide_pods,1).';

    jpods = ones(n,1) .* jpod;
    cojpods = (jpods ./ pods).';

end

function [pods,g,p] = validate_input(pods,g,p)

    pods = pods(:);
    n = numel(pods);

    if (n < 2)
        error('The value of ''pods'' is invalid. Expected input to be a vector containing at least 2 elements.');
    end

    k = n^2;

    [kg,ng] = size(g);

    if ((kg ~= k) || (ng ~= n))
        error(['The value of ''g'' is invalid. Expected input to be a matrix of size ' num2str(k) 'x' num2str(n) '.']);
    end

    kp = numel(p);

    if (kp ~= k)
        error(['The value of ''p'' is invalid. Expected input to be a vector containing ' num2str(k) ' elements.']);
    end

    p = p(:);

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
%% Example: Proton Treatment Plan with Manipulated CT values
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2017 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
% In this example we will show
% (i) how to load patient data into matRad
% (ii) how to setup a proton dose calculation
% (iii) how to inversely optimize the pencil beam intensities directly from command window in MATLAB.
% (iv) how to re-optimize a treatment plan
% (v) how to manipulate the CT cube by adding noise to the cube
% (vi) how to recalculate the dose considering the manipulated CT cube and the previously optimized pencil beam intensities
% (vii) how to compare the two results

%% Patient Data Import
% Let's begin with a clear Matlab environment and import the prostate
% patient into your workspace.

matRad_rc; %If this throws an error, run it from the parent directory first to set the paths

load('PROSTATE.mat');

%% Treatment Plan
% The next step is to define your treatment plan labeled as 'pln'. This
% structure requires input from the treatment planner and defines
% the most important cornerstones of your treatment plan.

pln.radiationMode           = 'protons';
pln.machine                 = 'generic_MCsquare';
pln.numOfFractions          = 30;
pln.propOpt.bioOptimization = 'const_RBExD';
pln.propStf.gantryAngles    = [90 270];
pln.propStf.couchAngles     = [0 0];
pln.propStf.bixelWidth      = 3;
pln.propStf.numOfBeams      = numel(pln.propStf.gantryAngles);
pln.propStf.isoCenter       = ones(pln.propStf.numOfBeams,1) * matRad_getIsoCenter(cst,ct,0);
pln.propOpt.runDAO          = 0;
pln.propOpt.runSequencing   = 0;

% dose calculation settings
pln.propDoseCalc.doseGrid.resolution.x = 3; % [mm]
pln.propDoseCalc.doseGrid.resolution.y = 3; % [mm]
pln.propDoseCalc.doseGrid.resolution.z = 3; % [mm]

%% Generate Beam Geometry STF
stf = matRad_generateStf(ct,cst,pln);

%% Dose Calculation
dij = matRad_calcParticleDoseMC(ct,stf,pln,cst);

%% Inverse Optimization for IMPT
resultGUI = matRad_fluenceOptimization(dij,cst,pln);

%% Calculate quality indicators
[dvh,qi]       = matRad_indicatorWrapper(cst,pln,resultGUI);
ixRectum       = 1;
display(qi(ixRectum).D_5);

%%
% Let's change the optimization parameter of the rectum in such a way that it
% will be better spared. We increase the penalty and lower the threshold
% of the squared overdose objective function. Afterwards we re-optimize
% the treatment plan and evaluate dose statistics one more time.

objective = cst{ixRectum,6}{1}; %This gives a struct
objective = matRad_DoseOptimizationFunction.createInstanceFromStruct(objective); %Now we turn it into a class
objective = objective.setDoseParameters(40); %We can simply call this function to change the/all dose parameter(s)
cst{ixRectum,6}{1} = struct(objective); % We put it back as struct

cst{ixRectum,6}{1}.parameters{1} = 40;
cst{ixRectum,6}{1}.penalty = 500;
resultGUI               = matRad_fluenceOptimization(dij,cst,pln);
[dvh2,qi2]              = matRad_indicatorWrapper(cst,pln,resultGUI);
display(qi2(ixRectum).D_5);

%% Plot the Resulting Dose Slice
% Let's plot the transversal iso-center dose slice
slice = round(pln.propStf.isoCenter(1,3)./ct.resolution.z);
figure
imagesc(resultGUI.RBExDose(:,:,slice)),colorbar, colormap(jet)

%%
% Now let's simulate a range undershoot by scaling the relative stopping power cube by 3.5% percent
ct_manip         = ct;
ct_manip.cubeHU{1} = 1.035*ct_manip.cubeHU{1};

%% Recalculate Plan with MC square
% Let's use the existing optimized pencil beam weights and recalculate the RBE weighted dose
resultGUI_noise = matRad_calcDoseDirectMC(ct_manip,stf,pln,cst,resultGUI.w);

%%  Visual Comparison of results
% Let's compare the new recalculation against the optimization result.
plane      = 3;
doseWindow = [0 max([resultGUI.RBExDose(:); resultGUI_noise.RBExDose(:)])];

figure,title('original plan')
matRad_plotSliceWrapper(gca,ct,cst,1,resultGUI.RBExDose,plane,slice,[],0.75,colorcube,[],doseWindow,[]);
figure,title('manipulated plan')
matRad_plotSliceWrapper(gca,ct_manip,cst,1,resultGUI_noise.RBExDose,plane,slice,[],0.75,colorcube,[],doseWindow,[]);

% Let's plot single profiles along the beam direction
ixProfileY = round(pln.propStf.isoCenter(1,1)./ct.resolution.x);

profileOrginal = resultGUI.RBExDose(:,ixProfileY,slice);
profileNoise   = resultGUI_noise.RBExDose(:,ixProfileY,slice);

figure,plot(profileOrginal,'LineWidth',2),grid on,hold on,
       plot(profileNoise,'LineWidth',2),legend({'original profile','noise profile'}),
       xlabel('mm'),ylabel('Gy(RBE)'),title('profile plot')

%% Quantitative Comparison of results
% Compare the two dose cubes using a gamma-index analysis.

doseDifference   = 2;
distToAgreement  = 2;
n                = 1;

[gammaCube,gammaPassRateCell] = matRad_gammaIndex(...
    resultGUI_noise.RBExDose,resultGUI.RBExDose,...
    [ct.resolution.x, ct.resolution.y, ct.resolution.z],...
    [doseDifference distToAgreement],slice,n,'global',cst);

repo: matRad


text:
function parse_data_pl2(filename)

% Requires Matlab Offline SDK (Plexon Inc) in the Matlab searchpath
% Based in code written by Azriel Ghadooshahy, Jan 2016

%Gather relevant info
pl2 = PL2GetFileIndex(filename);
[unused, fname, ext] = fileparts(filename);

%List of possible different continuous sources
%query source-channels for whether or not recording was enabled
source_names_cont = cellfun(@(x) x.SourceName,pl2.AnalogChannels,'uniformoutput',0); %it's possible that they haven't the same size
cont_record_chs = cellfun(@(x) x.RecordingEnabled,pl2.AnalogChannels,'uniformoutput',1);

%Load PL2 recorded continuous data
source_string = 'WB'; %could also be 'SPKC' or 'FP'
inds_in_source = ismember(source_names_cont,source_string);
inds_in_source = find(inds_in_source & cont_record_chs);
sourcenums = cellfun(@(x) x.Channel,pl2.AnalogChannels(inds_in_source),'uniformoutput',1);
%cdat = cell(1,numel(sourcenums));
channels = sourcenums;
fprintf('%d channels in the file\n', numel(sourcenums) )

for n = 1:numel(sourcenums)
    num_in_source = sourcenums(n);
    handler = PL2AdBySource(filename,source_string,num_in_source);
    lts = handler.FragCounts;
    sr = handler.ADFreq;
    fout = fopen([fname '_' num2str(num_in_source) '.pl2ch'],'w');
    fwrite(fout,handler.Values,'double');
    fclose(fout);
    fprintf('Channel %d out. \n', num_in_source);
end

save('pl2_meta_data','sr','lts','channels')

author: csn-le


text:
function [Pa] = inH2O2Pa(inH2O)
% Convert pressure from inches of water column at 4 degrees to pascals
% Chad Greene 2012
Pa = inH2O*249.089;

author: Sable


text:
function output = ...
  SynthesizeSpeech(spectrum_envelope, source)
% Synthesize speech from harmonic information
%
% output = ...
%  SynthesizeSpeech(spectrum_envelope, source)
%
% Arguments
% spectrum_envelope: structure with following fields
%   harmonic_power_dB: Harmonics level in dB
%   sampling_frequency: sampling frequency 9Hz)
% source: structure with following fields
%   f0: fundamental frequency (Hz)
%   aperiodicity_matrix: residual level in dB (Has to be calibrated)
%   vuv : voiced/unvoiced indicator of each frame: 1: voiced
%
% Return value
%   output : synthesized signal

% Copyright 2016 Google Inc. All Rights Reserved
% Author: hidekik@google.com (Hideki Kawahara)

% check input and prepare information for synthesis
narginchk(2, 2);
tx = source.frame_time;
f0 = source.f0;
fs = source.sampling_frequency;
harmonic_level = spectrum_envelope.harmonic_power_dB;
harmonic_deviation = source.aperiodicity_matrix;
unvoiced_frame = double(source.vuv < 0.5);
[event_index, event_locations, f0i] = CalculateEventLocations(tx, fs, f0);
unvoiced_sample = ...
  interp1(tx, double(unvoiced_frame), event_locations, 'linear', 'extrap');
harmonics_level_at_event = ...
  interp1(tx, harmonic_level', event_locations, 'linear', 'extrap')';
harmonics_deviation_at_event = ...
  StabilizeHarmonicDeviations(tx, harmonic_deviation, event_locations);
output = ...
  SynthesizeByEventBasedMethod(event_index, event_locations,...
  harmonics_level_at_event, ...
  harmonics_deviation_at_event, fs, f0i, ...
  unvoiced_sample);
end

author: google


text:
function slide34

	a = [0 0 0;
		3 0 0;
		3 1 0;
		0 1 0;
		0 0 1;
		3 0 1;
		3 1 1;
		0 1 1];

	b = [1 2 6 5;
		2 3 7 6;
		3 4 8 7;
		4 1 5 8;
		1 2 3 4;
		5 6 7 8];

	p1 = patch('faces',b,...
		'vertices',a,...
		'facecolor',[.5 .5 .5],...
		'edgecolor',[1,1,1],...
		'facealpha',0.5);

	view(3)
	axis([-3 7 -3 5 -3 5])
	grid on

end