text:
function status = besa_save2Swf(file_path, file_name, data_matrix, ...
    time_samples, source_labels, data_scale_factor, time_scale_factor)

% BESA_SAVE2SWF writes a data matrix into ASCII-source waveform file format
%
% Parameters:
%     [file_path]
%         Full path to the folder where the file should be saved.
%
%     [file_name]
%         The name of the file where the output should be written.
%
%     [data_matrix]
%         A 2D data matrix to be saved. It should have the size
%         [NumberOfChannels x NumberOfTimeSamples].
%
%     [time_samples]
%         An array with the time samples corresponding to the data. The
%         size of the array should be [1 x NumberOfTimeSamples].
%
%     [source_labels]
%         An array containing the channel labels for the current data.
%
%     [data_scale_factor]
%         A constant used for converting the data values from Tesla into fT
%         (10e15) or from Volt into uV (10e6). If the data is already in fT
%         or uV the just set it to 1.
%
%     [time_scale_factor]
%         A constant used for converting the time samples from seconds into
%         milliseconds (1000). If the time samples are already in
%         milliseconds just set the constant to 1.
%
% Return:
%     [status]
%         The status of the writing process. Still unused.
%

% Copyright (C) 2013, BESA GmbH
%
% File name: besa_save2Swf.m
%
% This file is part of MATLAB2BESA.
%
%    MATLAB2BESA 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.
%
%    MATLAB2BESA 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 MATLAB2BESA. If not, see <http://www.gnu.org/licenses/>.
%
% Author: Todor Jordanov
% Created: 2012-08-27

n_sources = size(source_labels, 2); %size(data_matrix, 1);

% Check if data is vectorized or multiplex.
isMul = 1;  % boolean variable containing the info if the file is
            %vectorized or multiplexed.
tmp1 = size(data_matrix);
if(tmp1(1) == n_sources)

    n_tbins = tmp1(2);
    isMul = 0;

else

    n_tbins = tmp1(1);

end

% Scale the time.
start_time = time_scale_factor * time_samples(1);
time_step = time_scale_factor * (time_samples(2) - time_samples(1));

% Scale the data.
data_matrix = data_matrix .* data_scale_factor;

% Change to the directory where the data should be saved.
cd(file_path)

% Open file for writing.
fid = fopen(file_name, 'w');

% Write the first line of the header.
fprintf(fid, ...
    'Npts= %d   TSB= %f DI= %f SB= 1.000 SC= 200.0 Nchan= %d\n', ...
    n_tbins, start_time, time_step, n_sources);

% Write the data matrix to the file.
for i=1:n_sources

    fprintf(fid, '%s: ', source_labels{i});

    for j=1:n_tbins

        if(isMul == 0)

            fprintf(fid, '%f ', data_matrix(i, j));

        else

            fprintf(fid, '%f ', data_matrix(j, i));

        end

    end

    fprintf(fid, '\n');

end

fclose(fid);

status = 1;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ h, hmin, hmax ] = p12_stepsize ( option )

%*****************************************************************************80
%
%% P12_STEPSIZE returns step sizes for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, real H, HMIN, HMAX, suggested values for the
%    initial step, the minimum step, and the maximum step.
%
  h =     2.000;
  hmin =  0.001;
  hmax = 10.000;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function X = Inv_CoarseCurveCoeff(C);
% Inv_CoarseCurveCoeff: Reconstruct the low-frequency subband from
%                           the coarsest scale curvelet coefficients
%  Usage:
%    X = Inv_CoarseCurveCoeff(C);
%  Inputs:
%    C   m by m matrix, m = 2*2^L
%  Outputs:
%    X   m by m matrix
% See Also
%   CoaseCurveCoeff, Inv_Curvelet02Xform
%
% By Emmanuel Candes, 2003-2004

X = fft2_mid0(C)/sqrt(prod(size(C)));

author: thomaskuestner


text:
%QUATERNION	constructor for quaternion objects
%
% 	QUATERNION([s v1 v2 v3])	from 4 elements
% 	QUATERNION(v, theta)		from vector plus angle
% 	QUATERNION(R)			from a 3x3 or 4x4 matrix
% 	QUATERNION(q)			from another quaternion

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function q = Quaternion(a1, a2)

	if nargin == 0,
		q.s = 1;
		q.v = [0 0 0];
		q = class (q, 'Quaternion');

	elseif nargin == 1
		if isa(a1, 'Quaternion')
			q = a1;
			q = class(q, 'Quaternion');
		elseif isreal (a1) && size(a1) == 1
			q.s = a1(1);
			q.v = [0,0,0];
			q = class(q, 'Quaternion');
		elseif isreal (a1) && all (size (a1) == [1 3]) # Quaternion (vector part)
			q.s = 0;
			q.v = a1(1:3);
			q = class(q, 'Quaternion');
		elseif all(size(a1) == [3 3])
			q = Quaternion( tr2q(a1) );
		elseif all(size(a1) == [4 4])
			q = Quaternion( tr2q(a1(1:3,1:3)) );
		elseif all(size(a1) == [1 4])
			q.s = a1(1);
			q.v = a1(2:4);
			q = class(q, 'Quaternion');
		else
			error('unknown dimension of input');
		end
	elseif nargin == 2
		if  isscalar(a1) && isvector(a2)
			q.s = cos(a1/2);
			q.v = (sin(a1/2)*unit(a2(:)'));
			q = class(q, 'Quaternion');
		end
	end
endfunction
%TR2Q	Convert homogeneous transform to a unit-quaternion
%
%	Q = tr2q(T)
%
%	Return a unit quaternion corresponding to the rotational part of the
%	homogeneous transform T.
%
%	See also Q2TR

% Ryan Steidnl based on Robotics Toolbox for MATLAB (v6 and v9)
%

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function q = tr2q(t)
	qs = sqrt(trace(t)+1)/2.0;
	kx = t(3,2) - t(2,3);	% Oz - Ay
	ky = t(1,3) - t(3,1);	% Ax - Nz
	kz = t(2,1) - t(1,2);	% Ny - Ox

	if (t(1,1) >= t(2,2)) & (t(1,1) >= t(3,3))
		kx1 = t(1,1) - t(2,2) - t(3,3) + 1;	% Nx - Oy - Az + 1
		ky1 = t(2,1) + t(1,2);			% Ny + Ox
		kz1 = t(3,1) + t(1,3);			% Nz + Ax
		add = (kx >= 0);
	elseif (t(2,2) >= t(3,3))
		kx1 = t(2,1) + t(1,2);			% Ny + Ox
		ky1 = t(2,2) - t(1,1) - t(3,3) + 1;	% Oy - Nx - Az + 1
		kz1 = t(3,2) + t(2,3);			% Oz + Ay
		add = (ky >= 0);
	else
		kx1 = t(3,1) + t(1,3);			% Nz + Ax
		ky1 = t(3,2) + t(2,3);			% Oz + Ay
		kz1 = t(3,3) - t(1,1) - t(2,2) + 1;	% Az - Nx - Oy + 1
		add = (kz >= 0);
	end

	if add
		kx = kx + kx1;
		ky = ky + ky1;
		kz = kz + kz1;
	else
		kx = kx - kx1;
		ky = ky - ky1;
		kz = kz - kz1;
	end
	nm = norm([kx ky kz]);
	if nm == 0,
		q = Quaternion([1 0 0 0]);
	else
		s = sqrt(1 - qs^2) / nm;
		qv = s*[kx ky kz];

		q = Quaternion([qs qv]);

	end
endfunction

author: petercorke

repo: robotics-toolbox-matlab


text:
function P=sympoly(varargin)
% sympoly: creator for sympoly objects
% usage: P=sympoly
% usage: P=sympoly(scalar_numeric_variable)
% usage: P=sympoly(array_numeric_variable)
% usage: P=sympoly('variablename')
% usage: sympoly varname1 varname2 varname3 ...
%
% Arguments:
%  varname - character string that represents a valid Matlab
%          variable name.
%
%          Example variable name: 'x'
%
%          or
%
%          any scalar constant or array of constants
%
% Examples of use:
%  P=sympoly       -->  creates a constant monomial, P(x) = 0
%
%  P=sympoly('x')  -->  creates a linear monomial, P(x) = x
%
%  P=sympoly(1)    -->  creates a constant monomial == 1, P(x) = 1
%
%  P=sympoly(hilb(3))  -->  creates a 3x3 matrix of constant sympoly
%                       variables, in this case a 3x3 Hilbert matrix.
%
%  sympoly a b x y -->  creates sympoly variables with those names in
%                       the caller workspace. This mode is equivalent
%                       to 4 calls:
%                        sympoly('a')
%                        sympoly('b')
%                        sympoly('x')
%                        sympoly('y')
%
%  Q=sympoly(P)    -->  copies an existing sympoly

% How was sympoly called? There are 5 options.

% sympoly('x') with an output argument --> create a linear
% sympoly variable, then return it in P.
if (nargin == 0)
  % this mode creates a scalar sympoly, as a constant == 0
  P = sympoly(0);

elseif (nargin==1)
  % Its only a single input. What kind?
  inp = varargin{1};

  if isa(inp,'sympoly')
    % Its already a sympoly. no-op.
    P = inp;
  elseif isstruct(inp) && isfield(inp,'Var') && isfield(inp,'Exponent') && isfield(inp,'Coefficient')
    % it is a struct, that wants desperately to be a sympoly
      % define this as a sympoly
      P=class(inp,'sympoly');
  elseif ischar(inp)
    % A string that contains the name of a new sympoly

    if nargout == 0
      % assign a variable with this name in the caller
      % workspace
      str = [inp,'=sympoly(''',inp,''');'];

      % evalin will put it in place in the caller workspace
      evalin('caller',str)
    else
      % Create a sympoly, to be returned as an output
      P.Var = {inp};
      P.Exponent = 1;
      P.Coefficient = 1;

      % define this as a sympoly
      P=class(P,'sympoly');

    end

  elseif isnumeric(inp) || isa(inp,'logical')
    % A numeric scalar or array

    % initialize
    P.Var = {''};
    P.Exponent = 0;
    P.Coefficient = 0;

    % replicate
    P = repmat(P,size(inp));

    % and stuff the coefficients
    for i = 1:numel(inp)
      P(i).Coefficient = inp(i);
    end

    % define this as a sympoly
    P=class(P,'sympoly');

  end

else
  % nargin is > 1. nargout must = 0
  if nargout>0
    error 'Mulitple sympolys can only be created as: "sympoly x y z w"'
  end

  % just loop over the inputs
  for i = 1:nargin
    % A list of several variables to create

    % assign variables with this name in the caller
    % workspace
    inp = varargin{i};
    str = [inp,'=sympoly(''',inp,''');'];

    % evalin will put it in place in the caller workspace
    evalin('caller',str)

  end
end

% make sure that the appropriate sympoly method is
% used whenever one of the arguments is a sympoly
% and any numeric type as the other argument.
superiorto('double','single','int8','uint8','int16', ...
  'uint16','int32','uint32','int64','uint64','logical')

author: Sable

repo: mcbench-benchmarks


text:
function [mu, B] = clg_Mstep_simple(w, Y, YY, YTY, X, XX, XY)
% CLG_MSTEP_SIMPLE Same as CLG_MSTEP, but doesn;t estimate Sigma,  so is slightly faster
% function [mu, B] = clg_Mstep_simple(w, Y, YY, YTY, X, XX, XY)
%
% See clg_Mstep for details.
% Unlike clg_Mstep, there are no optional arguments, which are slow to process
% if this function is inside a tight loop.

[Ysz Q] = size(Y);

if isempty(X) % no regression
  %B = [];
  B2 = zeros(Ysz, 1, Q);
  for i=1:Q
    B(:,:,i) = B2(:,1:0,i); % make an empty array of size Ysz x 0 x Q
  end
  [mu, Sigma] = mixgauss_Mstep(w, Y, YY, YTY);
  return;
end

N = sum(w);
%YY = YY + cov_prior; % regularize the scatter matrix

% Set any zero weights to one before dividing
% This is valid because w(i)=0 => Y(:,i)=0, etc
w = w + (w==0);

Xsz = size(X,1);
% Append 1 to X to get Z
ZZ = zeros(Xsz+1, Xsz+1, Q);
ZY = zeros(Xsz+1, Ysz, Q);
for i=1:Q
  ZZ(:,:,i) = [XX(:,:,i)  X(:,i);
	       X(:,i)'    w(i)];
  ZY(:,:,i) = [XY(:,:,i);
	       Y(:,i)'];
end

mu = zeros(Ysz, Q);
B = zeros(Ysz, Xsz, Q);
for i=1:Q
  % eqn 9
  if rcond(ZZ(:,:,i)) < 1e-10
    sprintf('clg_Mstep warning: ZZ(:,:,%d) is ill-conditioned', i);
    %probably because there are too few cases for a high-dimensional input
    ZZ(:,:,i) = ZZ(:,:,i) + 1e-5*eye(Xsz+1);
  end
  %A = ZY(:,:,i)' * inv(ZZ(:,:,i));
  A = (ZZ(:,:,i) \ ZY(:,:,i))';
  B(:,:,i) = A(:, 1:Xsz);
  mu(:,i) = A(:, Xsz+1);
end


text:
function B = biharm(f)
%BIHARM   Biharmonic operator applied to a CHEBFUN3.
%   B = BIHARM(F) returns a CHEBFUN3 object B representing the biharmonic
%   operator applied to a CHEBFUN3 object F.
%
%   This is shorthand for BIHARMONIC(F).
%
% See also CHEBFUN3/BIHARMONIC.

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

% Call Biharmonic:
B = biharmonic(f);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%GET_BOUNDARY_LAYERS Find boundary layers.
%
%   [ BL ] = GET_BOUNDARY_LAYERS( FEA, IND_BL, HL, NL, AR, ARG )
%   Find boundary layers in FEA struct and return the struct BL
%   with the following fields for connected wall boundaries
%
%       boundaries - list of boundaries in layer
%       svertices  - list of single vertices in layer
%       sedges     - list of single edges in layer (only 3D)
%       gmsh       - gmsh boundary field
%
%   The optional array IND_BL is used to override automatic detection
%   of wall boundaries (boundary condition selection 1 for
%   Navier-Stokes, swirl flow, Brinkman equations, and fluid-structure
%   interation physics modes). If IND_BL is empty layers for all
%   boundaries will be returned.
%
%   The optional numerical arguments HL, NL (default 5), and AR
%   (default 1.2) specifies the layer height, number of cells in
%   layer, and aspect ratios for gmsh output.
%
%   ARG is an optional struct to give parameters for Gmsh boundary
%   layer generation with the fields hmax, hmaxb, and thmax (default 5*pi/3).

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
classdef TestEMDL1
    %TestEMDL1

    methods (Static)
        function test_1
            sig1 = randn(20,1);
            sig2 = randn(20,1);
            d = cv.EMDL1(sig1, sig2);
            validateattributes(d, {'numeric'}, {'scalar', 'real'});
        end

        function test_2
            sig1 = randn(20,1);
            d = cv.EMDL1(sig1, sig1);
            validateattributes(d, {'numeric'}, {'scalar', 'real'});
        end

        function test_error_argnum
            try
                cv.EMDL1();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function tree = struct2xml(s,rootname)
%STRUCT2XML Convert structure to an XML tree object and specify its root name
%   TREE = STRUCT2XML(S,ROOTNAME) converts the structure S into an XML
%   representation TREE (an XMLTree object) with ROOTNAME as the root tag, if
%   provided. Only conventional objects (char, double, sparse, *int*) are
%   accepted in S's fields.
%
%   Example
%      report = struct('name','John','marks',...
%                      struct('maths',17,'physics',12));
%      tree = struct2xml(report);
%      save(tree,'report.xml');
%
%   See also XMLTREE.

%   Copyright 2003 Guillaume Flandin.
%   $Revision: 1.1 $  $Date: 2003/03/11 21:26: $

error(nargchk(1,2,nargin));
if ~isstruct(s)
	error('[STRUCT2XML] Input argument must be a struct.');
end
if nargin == 1,
	if ~isempty(inputname(1)),
		rootname=inputname(1);
	else
		rootname='root';
	end
end
% Create an empty XML tree
tree = xmltree;

% Root element is the input argument name
tree = set(tree,root(tree),'name',rootname);

% Recursively walk inside the structure
tree = sub_struct2xml(tree,s,root(tree));

%=======================================================================
function tree = sub_struct2xml(tree,s,uid)

	if isstruct(s)
		for k=1:length(s)
			names = fieldnames(s(k));
			for i=1:length(names)
				if iscell(getfield(s(k),names{i}))
					for j=1:length(getfield(s(k),names{i}))
						[tree, uid2] = add(tree,uid,'element',names{i});
						tree = sub_struct2xml(tree,getfield(s(k),names{i},{j}),uid2);
					end
				else
					[tree, uid2] = add(tree,uid,'element',names{i});
					tree = sub_struct2xml(tree,getfield(s(k),names{i}),uid2);
				end
			end
		end
	else
		switch class(s)
			case 'char'
				tree = add(tree,uid,'chardata',s); %need to handle char arrays...
			case 'cell'
				% if a cell is present here, it comes from: getfield(s(k),names{i},{j})
				tree = sub_struct2xml(tree,s{1},uid);
			case {'double','sparse'}
				tree = add(tree,uid,'chardata',sub_num2str(s));
			case {'int8','uint8','int16','uint16','int32','uint32'}
				% need a specific function because num2str only works with double
			otherwise
				error(sprintf('[STRUCT2XML] Cannot convert from %s to char.',class(s)));
		end
	end

%=======================================================================
function s = sub_num2str(n)  % to be improved for ND arrays (code from Ph. Ciuciu)
	[N,P] = size(n);

	if N>1 | P>1
		s=['['];
		w = ones(1,P);w(P)=0;
		v = ones(N,1);v(N)=0;
		for k=1:N
			for i=1:P
				s = [s num2str(n(k,i)) repmat(',',1,w(i))];
			end
			s = [s repmat(';',1,v(k))];
		end
		s = [s ']'];
	else
		s = num2str(n);
	end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Local Regression and Likelihood, Figure 11.1
%
% Variable degree fit, uses module 'vord'. Note there is
% no `lower' degree here; it defaults to 0.
%
% Author: Catherine Loader

load ethanol;

figure('Name','fig11_1a: Variable degree fit');
fit = locfit(E,NOx,'deg',3,'nn',0.3,'module','vord');
lfplot(fit);

figure('Name','fig11_1b: Variable degree fit');
x = fit.fit_points.evaluation_points';
z = predict(fit,'fitp','what','deg');
plot(x,z,'o');
xlabel('Fitting Point');
ylabel('Degree');


text:
function bbci_trigger(varargin)
%BBCI_TRIGGER - Sending triggers to the EEG (or whatever)
%
%Synopsis:
%  bbci_trigger(VALUE)
%  bbci_trigger(BBCI, VALUE)
%
%Arguments:
%  VALUE - Value that is sent as a trigger. Most trigger functions will
%          accept only integers from 1 to 255 as VALUE.
%  BBCI - BBCI variable. This function uses the field BBCI.trigger.
%
%Description:
%  This function calls one of a selection of specific trigger functions.
%  If the BBCI variable is specified, BBCI.trigger.fcn is called with
%    BBCI.trigger.param{:} as additional parameters.
%  Otherwise, the function BTB.Acq.TriggerFcn is called (global variable BTB)
%  with BTB.Acq.TriggerParam{:} as addiotional parameters.

if nargin>1 && isstruct(varargin{1}),
  bbci= varargin{1};
  bbci.trigger.fcn(varargin{2}, bbci.trigger.param{:});
elseif nargin==1,
  global BTB
  BTB.Acq.TriggerFcn(varargin{1}, BTB.Acq.TriggerParam{:});
end

author: bbci


text:
function [cost]=costf2(x)
global impvol; global strike; global T; global F0; global r;

for i=1:length(T)
%mok(i)=max(blsimpv(F0, strike(i), r,T(i),HestonCall(F0,strike(i),r,T(i),x(1),x(2),x(3),x(4),x(5),0), 3),0);
%cost(i)=impvol(i)-mok(i);
cost(i)=blsprice(F0,strike(i),r,T(i),impvol(i))-HestonCall(F0,strike(i),r,T(i),x(1),x(2),x(3),x(4),x(5),0);
end


text:
function energy = computeEnergyBinaryPairwise( unaryTerms, pairwiseTerms, labels )
%computeEnergyBinaryPairwise computes the value of the energy with unary and pairwise potentials

%% check input
numLabels = 2;

if ~isnumeric(unaryTerms) || ~ismatrix(unaryTerms) || size(unaryTerms, 2) ~= 2
    error('Incorrect format for unaryTerms, has to be numNodes x 2')
end
numNodes = size(unaryTerms, 1);

if ~isnumeric(pairwiseTerms) || ~ismatrix(pairwiseTerms) || size(pairwiseTerms, 2) ~= 6
    error('Incorrect format for pairwiseTerms, has to be numEdges x 6')
end
numEdges = size(pairwiseTerms, 1);

if ~isnumeric(labels) || ~isvector(labels) || length(labels) ~= numNodes
    error('Incorrect format for labels, has to be numNodes x 1')
end
labels = labels(:);
if any(labels > numLabels) || any(labels < 1)
    error('Incorrect values for labels, has to be an integer from 1 to numLabels')
end

%% computation
energy = sum( unaryTerms((1 : numNodes)' + numNodes * (labels - 1) ) );

label1 = labels(pairwiseTerms(:, 1));
label2 = labels(pairwiseTerms(:, 2));
labelMap = [ 3, 4; 5, 6];
jointLabelMap = labelMap( label1 + 2 * (label2 - 1) );

energy = energy + sum( pairwiseTerms( (1 : numEdges)' + numEdges * (jointLabelMap - 1) ) );

end

author: aosokin


text:
function [outexpdir] = ConvertLarvaeSpecies2Trx(inexpdir,rootoutputdir,varargin)

%% set parameters

% default parameters
trxfilestr = 'trx.mat';
moviefilestr = 'movie.wmv';
makesoftlink = true;
intrxfilestr = 'trx.mat';
inmovieformat = 'wmv';
arena_radius_mm = 60;
swapxy = false;

[trxfilestr,moviefilestr,makesoftlink,intrxfilestr,inmovieformat,arena_radius_mm,swapxy] = ...
  myparse(varargin,...
  'trxfilestr',trxfilestr,...
  'moviefilestr',moviefilestr,...
  'makesoftlink',makesoftlink,...
  'intrxfilestr',intrxfilestr,...
  'inmovieformat',inmovieformat,...
  'arena_radius_mm',arena_radius_mm,...
  'swapxy',swapxy);

%% get paths to input data

[rootdatadir,inexperiment_name] = fileparts(inexpdir); %#ok<ASGLU>

files = dir(fullfile(inexpdir,['*.',inmovieformat]));
if isempty(files),
  error('No %s movie found in %s',inmovieformat,inexpdir);
end
if numel(files) > 1,
  warning('More than one %s movie found in %s, choosing %s',inmovieformat,inexpdir,files(1).name);
  files = files(1);
end
inmoviefile = fullfile(inexpdir,files(1).name);

intrxfile = fullfile(inexpdir,intrxfilestr);
if ~exist(intrxfile,'file'),
  error('Could not find trx file %s',intrxfile);
end

%% create output directory

% make sure the root output directory exists
if ~exist(rootoutputdir,'dir'),
  mkdir(rootoutputdir);
end

% create the experiment directory
outexperiment_name = inexperiment_name;
outexpdir = fullfile(rootoutputdir,outexperiment_name);
if ~exist(outexpdir,'dir'),
  mkdir(outexpdir);
end

%% load in trx data
%
% load(intrxfile,'trx');
% nlarvae = numel(trx); %#ok<NODEF>

%% get arena location
%
% arena = full(trx(1).arena);
% stats = regionprops(arena,'centroid','majoraxislength');
% arena_center_x = stats.Centroid(1);
% arena_center_y = stats.Centroid(2);
% arena_radius = stats.MajorAxisLength/2;

%% swap x and y if nec
%
% if swapxy,
%   for i = 1:nlarvae,
%     tmp = trx(i).x;
%     trx(i).x = trx(i).y; %#ok<*AGROW>
%     trx(i).y = tmp;
%     isdata = ~cellfun(@isempty,trx(i).spine);
%     trx(i).spine(isdata) = cellfun(@(x) x(:,[2,1]),trx(i).spine(isdata),'UniformOutput',false);
%     trx(i).contour = cellfun(@(x) x(:,[2,1]),trx(i).contour,'UniformOutput',false);
%   end
% end

%% remove extra stuff at the end of fields

% fns = {'x','y','a','b','theta','area'};
% for i = 1:numel(trx),
%   isdiscrepency = false;
%   for j = 1:numel(fns),
%     fn = fns{j};
%     if numel(trx(i).(fn)) ~= trx(i).nframes,
%       if ~isdiscrepency,
%         fprintf('Warning: for larva %d, size of fields and nframes = %d, endframe=%d - firstframe=%d -> %d do not match: ',...
%           i,trx(i).nframes,trx(i).lastframe,trx(i).firstframe,trx(i).lastframe-trx(i).firstframe+1);
%         isdiscrepency = true;
%       end
%       fprintf('|%s|=%d ',fn,numel(trx(i).(fn)));
%       trx(i).(fn) = trx(i).(fn)(1:trx(i).nframes); %#ok<AGROW>
%     end
%   end
%   if isdiscrepency,
%     fprintf('\n');
%   end
% end

%% convert to mm

% pxpermm = arena_radius / arena_radius_mm;
%
% for i = 1:nlarvae,
%   trx(i).x_mm = (trx(i).x - arena_center_x) / pxpermm;
%   trx(i).y_mm = (trx(i).y - arena_center_y) / pxpermm;
%   trx(i).a_mm = trx(i).a / pxpermm;
%   trx(i).b_mm = trx(i).b / pxpermm;
%   trx(i).theta_mm = trx(i).theta;
%   trx(i).pxpermm = pxpermm;
%   % convert spines to matrix
%   if i == 1,
%     nspinepts = max(cellfun(@(x) size(x,1),trx(i).spine));
%   end
%   trx(i).xspine = nan(nspinepts,trx(i).nframes);
%   trx(i).yspine = nan(nspinepts,trx(i).nframes);
%   isdata = ~cellfun(@isempty,trx(i).spine);
%   trx(i).xspine(:,isdata) = cell2mat(cellfun(@(x) x(:,1),trx(i).spine(isdata),'UniformOutput',false));
%   trx(i).yspine(:,isdata) = cell2mat(cellfun(@(x) x(:,2),trx(i).spine(isdata),'UniformOutput',false));
%   trx(i).xspine_mm = (trx(i).xspine - arena_center_x) / pxpermm;
%   trx(i).yspine_mm = (trx(i).yspine - arena_center_y) / pxpermm;
%
%   trx(i).xcontour = cellfun(@(x) x(:,1),trx(i).contour,'UniformOutput',false);
%   trx(i).ycontour = cellfun(@(x) x(:,2),trx(i).contour,'UniformOutput',false);
%   trx(i).xcontour_mm = cellfun(@(x) (x - arena_center_x)/pxpermm,trx(i).xcontour,'UniformOutput',false);
%   trx(i).ycontour_mm = cellfun(@(x) (x - arena_center_y)/pxpermm,trx(i).ycontour,'UniformOutput',false);
%
%   trx(i).area_mm = trx(i).area / pxpermm^2;
%
% end
%
% %% create dt, compress arena, rename lastframe as endframe
%
% for i = 1:nlarvae,
%   trx(i).dt = repmat(1/trx(i).fps,[1,trx(i).nframes-1]);
%   trx(i).arena = struct('arena_center_mm_x',0,'arena_center_mm_y',0,'arena_radius_mm',arena_radius_mm,...
%     'arena_center_px_x',arena_center_x,'arena_center_px_y',arena_center_y,'arena_radius_px',arena_radius);
%   trx(i).endframe = trx(i).lastframe;
% end
%
% %% remove redundant fields
%
% trx = rmfield(trx,{'spine','contour','lastframe'}); %#ok<NASGU>
%
% %% create the trx file
%
% outtrxfile = fullfile(outexpdir,trxfilestr);
% save(outtrxfile,'trx');

%% copy over the movie file

outmoviefile = fullfile(outexpdir,moviefilestr);
if isunix && makesoftlink,
  if exist(outmoviefile,'file'),
    delete(outmoviefile);
  end
  cmd = sprintf('ln -s %s %s',inmoviefile,outmoviefile);
  unix(cmd);
else
  [success,msg] = copyfile(inmoviefile,outmoviefile);
  if ~success,
    error('Error copying file %s to %s: %s',inmoviefile,outmoviefile,msg);
  end
end

author: kristinbranson

repo: JAABA


text:
function [ xtab, weight ] = legendre_compute ( norder )

%*****************************************************************************80
%
%% LEGENDRE_COMPUTE computes a Gauss-Legendre quadrature rule.
%
%  Discussion:
%
%    The integration interval is [ -1, 1 ].
%
%    The weight function is w(x) = 1.0.
%
%    The integral to approximate:
%
%      Integral ( -1 <= X <= 1 ) F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= NORDER ) WEIGHT(I) * F ( XTAB(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2005
%
%  Author:
%
%    Original FORTRAN77 version by Philip Davis, Philip Rabinowitz
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%  Parameters:
%
%    Input, integer NORDER, the order of the rule.
%    NORDER must be greater than 0.
%
%    Output, real XTAB(NORDER), the abscissas of the rule.
%
%    Output, real WEIGHT(NORDER), the weights of the rule.
%    The weights are positive, symmetric, and should sum to 2.
%
  if ( norder < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEGENDRE_COMPUTE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of NORDER = %d\n', norder );
    error ( 'LEGENDRE_COM - Fatal error!' );
  end

  e1 = norder * ( norder + 1 );

  m = floor ( ( norder + 1 ) / 2 );

  for i = 1 : floor ( ( norder + 1 ) / 2 )

    mp1mi = m + 1 - i;

    t = ( 4 * i - 1 ) * pi / ( 4 * norder + 2 );

    x0 = cos(t) * ( 1.0 - ( 1.0 - 1.0 / ( norder ) ) / ( 8 * norder * norder ) );

    pkm1 = 1.0;
    pk = x0;

    for k = 2 : norder
      pkp1 = 2.0 * x0 * pk - pkm1 - ( x0 * pk - pkm1 ) / k;
      pkm1 = pk;
      pk = pkp1;
    end

    d1 = norder * ( pkm1 - x0 * pk );

    dpn = d1 / ( 1.0 - x0 * x0 );

    d2pn = ( 2.0 * x0 * dpn - e1 * pk ) / ( 1.0 - x0 * x0 );

    d3pn = ( 4.0 * x0 * d2pn + ( 2.0 - e1 ) * dpn ) / ( 1.0 - x0 * x0 );

    d4pn = ( 6.0 * x0 * d3pn + ( 6.0 - e1 ) * d2pn ) / ( 1.0 - x0 * x0 );

    u = pk / dpn;
    v = d2pn / dpn;
%
%  Initial approximation H:
%
    h = - u * ( 1.0 + 0.5 * u * ( v + u * ( v * v - d3pn / ( 3.0 * dpn ) ) ) );
%
%  Refine H using one step of Newton's method:
%
    p = pk + h * ( dpn + 0.5 * h * ( d2pn + h / 3.0 ...
      * ( d3pn + 0.25 * h * d4pn ) ) );

    dp = dpn + h * ( d2pn + 0.5 * h * ( d3pn + h * d4pn / 3.0 ) );

    h = h - p / dp;

    xtemp = x0 + h;

    xtab(mp1mi) = xtemp;

    fx = d1 - h * e1 * ( pk + 0.5 * h * ( dpn + h / 3.0 ...
      * ( d2pn + 0.25 * h * ( d3pn + 0.2 * h * d4pn ) ) ) );

    weight(mp1mi) = 2.0 * ( 1.0 - xtemp * xtemp ) / fx / fx;

  end

  if ( mod ( norder, 2 ) == 1 )
    xtab(1) = 0.0;
  end
%
%  Shift the data up.
%
  nmove = floor ( ( norder + 1 ) / 2 );
  ncopy = norder - nmove;

  for i = 1 : nmove
    iback = norder + 1 - i;
    xtab(iback) = xtab(iback-ncopy);
    weight(iback) = weight(iback-ncopy);
  end
%
%  Reflect values for the negative abscissas.
%
  for i = 1 : norder - nmove
    xtab(i) = - xtab(norder+1-i);
    weight(i) = weight(norder+1-i);
  end

  return
end


text:
function z = piecewise_eval(x,breakpoints,funs)
% PIECEWISE_EVAL: evaluates a piecewise function of x
% usage: y = PIECEWISE_EVAL(x,breakpoints,funs)
%
% arguments (input)
%  x    - vector or array of points to evaluate though the function
%
%  breakpoints - list of n breakpoints, -inf and +inf are implicitly
%         the first and last breakpoints. A function with only two
%         pieces has only one explicit breakpoint. In the event that
%         you want to define a function with breakpoints [a,b,c],
%         and only two functions, but you do not care what happens
%         for x < a or x > b, then you should specify only the
%         breakpoint b. Alternatively, one could specify all 3
%         breaks, and force the function to return NaN above and
%         below those limits.
%
%         x(i) will be identified as falling in interval (j) if
%         break(j) <= x(i) < break(j+1)
%
%  funs - cell array containing n+1 functions as scalar constants,
%         strings, anonymous functions, inline functions, or any
%         simple matlab function name.
%
%         Note: use .*, ./, .^ where appropriate in the function
%
%         These functions need not be differentiable or even
%         continuous across the breaks.
%
% arguments (output)
%  z    - evaluated function, result is same shape as x
%
% Example usage:
%  For       x < -5, y = 2
%  For -5 <= x < 0,  y = sin(x)
%  For  0 <= x < 2,  y = x.^2
%  For  2 <= x < 3,  y = 6
%  For  3 <= x,      y = inf
%
%  y = piecewise_eval(-10:10,[-5 0 2 3],{2,'sin(x)','x.^2',6,inf})

n=length(breakpoints);
% there must be n+1 funs for n breaks
if length(funs)~=(n+1)
  error 'funs and breakpoints are incompatible in size'
end

if any(diff(breakpoints)<=0)
  error 'Breakpoints must be both distinct and increasing'
end

% ensure the functions are feval-able
for i=1:(n+1)
  if ischar(funs{i})
    % A string. Make it a function
    f=inline(funs{i});
    funs{i} = f;
  elseif isa(funs{i},'function_handle') || isa(funs{i},'inline')
    % a function handle or an inline. do nothing.
  elseif isnumeric(funs{i}) | isnan(funs{i}) | isinf(funs{i})
    % A scalar value was supplied, may be NaN or inf.
    % Make it a function.
    funs{i}=@(x) funs{i};
  else
    % It must be something that feval can handle
    % directly, so leave funs{i} alone.
  end
end

% initialize as nans
z=nan(size(x));

% below the first break
k=(x<breakpoints(1));
z(k)=feval(funs{1},x(k));

left = k;
for i=2:n
  k=(~left) & (x<breakpoints(i));
  if any(k)
    z(k)=feval(funs{i},x(k));
    left = k | left;
  end
end

% over the top
k=(x>=breakpoints(end));
z(k)=feval(funs{end},x(k));


text:
clear; clc;

addpath(genpath('../libs/exportFig'));
addpath(genpath('../libs/fun4MeanShift'));
addpath(genpath('../libs/layerExt'));
addpath(genpath('../libs/layerExt'));
path_to_matconvnet = '../libs/matconvnet-1.0-beta23_modifiedDagnn';
run(fullfile(path_to_matconvnet, 'matlab', 'vl_setupnn'));
addpath(genpath(fullfile('dependencies', 'matconvnet','examples')));
%% read matconvnet model
load('imdb_complete_on_server.mat');
imdb.meta.meanvalue = reshape(imdb.meta.meanvalue,[1 1 3]);

% set GPU
gpuId = 1;
gpuDevice(gpuId);

modelName = 'basemodel_resnet101_div4_loopNum0.mat';
netbasemodel = load( fullfile('./basemodel', modelName) );
netbasemodel = dagnn.DagNN.loadobj(netbasemodel);
netbasemodel.move('gpu') ;

netbasemodel.removeLayer('obj_MM');
netbasemodel.removeLayer('obj_reg');
netbasemodel.removeLayer('output_cosSim');
netbasemodel.removeLayer('output_l2norm');

%netbasemodel.mode = 'test' ;
netbasemodel.mode = 'normal' ;
net.conserveMemory = 1;
sName = 'res9_conv';

lName = 'output_l2norm';
netbasemodel.addLayer(lName, L2normalization(), sName, lName) ;
sName = lName;
outputIdx_l2norm = netbasemodel.layers(netbasemodel.getLayerIndex('output_l2norm')).outputIndexes;
netbasemodel.vars(outputIdx_l2norm).precious = 1;

loopNum = 16; % let it loop and visualize the intermediate results
predInstanceMaskMat_loop0 = cell(1,3);
predInstanceMaskMat_loops = cell(loopNum,3);

sName_l2norm = lName;
for loopIdx = 1:loopNum
    [netbasemodel, sName, sName_l2norm] = addOneLoop_forMeanShiftGrouping(netbasemodel, sName_l2norm, loopIdx);

    keepLayerName = sprintf('loop%d_meanshift_S_is_XX', loopIdx);
    netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).block.analyzeGradient = false;

    keepLayerName = sprintf('loop%d_meanshift_Y_l2norm', loopIdx);
    netbasemodel.vars(netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).outputIndexes).precious = 1;

    rmLayerName = sprintf('loop%d_meanshift_cosSim', loopIdx);
    netbasemodel.removeLayer(rmLayerName); % remove layer
end
%% visualize single image
prefix4saving = './results/'; % all visual results are saved here
if ~isdir(prefix4saving)
    mkdir(prefix4saving);
end
path_to_image = './images/';
imgList = dir(fullfile(path_to_image, '*jpg'));
for idx = 1:length(imgList)
    fprintf('%d %s\n', idx, imgList(idx).name)
    cur_path_to_image = fullfile(path_to_image, imgList(idx).name);
    imgOrg = imread(cur_path_to_image);

    [~,curFileName,curFileExt] = fileparts(cur_path_to_image);
    cur_path_to_annot = fullfile(path_to_image, [curFileName,'.mat']);
    cur_segMap = load(cur_path_to_annot);
    cur_instMap = cur_segMap.gtMat.instMap;
    cur_segMap = cur_segMap.gtMat.segMap;
    gtOrg = cur_instMap;

    sz = size(gtOrg);
    reSZ = round(sz/8)*8;

    img = imresize(imgOrg, reSZ);
    im = bsxfun(@minus, single(img), imdb.meta.meanvalue) ;
    inputs = {'data', gpuArray(single(im))};

    netbasemodel.eval(inputs) ;
    %% gather result
    feaMap = gather(netbasemodel.vars(outputIdx_l2norm).value);
    feaMapSize = size(feaMap);
    predInstanceMaskMat_loop0{1} = reshape(feaMap, [], size(feaMap,3));
    rng(777); randProj = randn(feaMapSize(3), 3);
    predInstanceMaskMat_loop0{1} = predInstanceMaskMat_loop0{1}*randProj;
    predInstanceMaskMat_loop0{1} = reshape( predInstanceMaskMat_loop0{1}, [feaMapSize(1), feaMapSize(2), 3] );
    predInstanceMaskMat_loop0{1} = rescaleFeaMap(predInstanceMaskMat_loop0{1});
    predInstanceMaskMat_loop0{1} = imresize(predInstanceMaskMat_loop0{1}, [sz(1),sz(2)]);

    rng(77); randProj = randn(feaMapSize(3), 3);
    predInstanceMaskMat_loop0{2} = reshape(feaMap, [], size(feaMap,3));
    predInstanceMaskMat_loop0{2} = predInstanceMaskMat_loop0{2}*randProj;
    predInstanceMaskMat_loop0{2} = reshape( predInstanceMaskMat_loop0{2}, [feaMapSize(1), feaMapSize(2), 3] );
    predInstanceMaskMat_loop0{2} = rescaleFeaMap(predInstanceMaskMat_loop0{2});
    predInstanceMaskMat_loop0{2} = imresize(predInstanceMaskMat_loop0{2}, [sz(1),sz(2)]);

    rng(7); randProj = randn(feaMapSize(3), 3);
    predInstanceMaskMat_loop0{3} = reshape(feaMap, [], size(feaMap,3));
    predInstanceMaskMat_loop0{3} = predInstanceMaskMat_loop0{3}*randProj;
    predInstanceMaskMat_loop0{3} = reshape( predInstanceMaskMat_loop0{3}, [feaMapSize(1), feaMapSize(2), 3] );
    predInstanceMaskMat_loop0{3} = rescaleFeaMap(predInstanceMaskMat_loop0{3});
    predInstanceMaskMat_loop0{3} = imresize(predInstanceMaskMat_loop0{3}, [sz(1),sz(2)]);

    for loopidx = 1:loopNum
        keepLayerName = sprintf('loop%d_meanshift_Y_l2norm', loopIdx);
        keepLayerName = netbasemodel.layers(netbasemodel.getLayerIndex(keepLayerName)).outputIndexes;
        feaMap = gather(netbasemodel.vars(keepLayerName).value);
        feaMapSize = size(feaMap);
        predInstanceMaskMat_loops{loopidx,1} = reshape(feaMap, [], size(feaMap,3));
        rng(777); randProj = randn(feaMapSize(3), 3);
        predInstanceMaskMat_loops{loopidx,1} = predInstanceMaskMat_loops{loopidx,1}*randProj;
        predInstanceMaskMat_loops{loopidx,1} = reshape( predInstanceMaskMat_loops{loopidx,1}, [feaMapSize(1), feaMapSize(2), 3] );
        predInstanceMaskMat_loops{loopidx,1} = rescaleFeaMap(predInstanceMaskMat_loops{1});
        predInstanceMaskMat_loops{loopidx,1} = imresize(predInstanceMaskMat_loops{loopidx,1}, [sz(1),sz(2)]);

        rng(77); randProj = randn(feaMapSize(3), 3);
        predInstanceMaskMat_loops{loopidx,2} = reshape(feaMap, [], size(feaMap,3));
        predInstanceMaskMat_loops{loopidx,2} = predInstanceMaskMat_loops{loopidx,2}*randProj;
        predInstanceMaskMat_loops{loopidx,2} = reshape( predInstanceMaskMat_loops{loopidx,2}, [feaMapSize(1), feaMapSize(2), 3] );
        predInstanceMaskMat_loops{loopidx,2} = rescaleFeaMap(predInstanceMaskMat_loops{loopidx,2});
        predInstanceMaskMat_loops{loopidx,2} = imresize(predInstanceMaskMat_loops{loopidx,2}, [sz(1),sz(2)]);

        rng(7); randProj = randn(feaMapSize(3), 3);
        predInstanceMaskMat_loops{loopidx,3} = reshape(feaMap, [], size(feaMap,3));
        predInstanceMaskMat_loops{loopidx,3} = predInstanceMaskMat_loops{loopidx,3}*randProj;
        predInstanceMaskMat_loops{loopidx,3} = reshape( predInstanceMaskMat_loops{loopidx,3}, [feaMapSize(1), feaMapSize(2), 3] );
        predInstanceMaskMat_loops{loopidx,3} = rescaleFeaMap(predInstanceMaskMat_loops{loopidx,3});
        predInstanceMaskMat_loops{loopidx,3} = imresize(predInstanceMaskMat_loops{loopidx,3}, [sz(1),sz(2)]);
    end

    %% show figures
    imgFig = figure(1);
    set(imgFig, 'Position', [100 100 1500 800]) % [1 1 width height]
    subplot(2,3,1);
    imshow(uint8(img)); title('original image');
    subplot(2,3,2);
    imagesc(gtOrg); axis off image; title('instance mask');
    subplot(2,3,4);
    imagesc(predInstanceMaskMat_loop0{1}); axis off image; title('Loop-0 randProj 3-dim')
    subplot(2,3,5);
    imagesc(predInstanceMaskMat_loop0{2}); axis off image; title('Loop-0 randProj 3-dim')
    subplot(2,3,6);
    imagesc(predInstanceMaskMat_loop0{3}); axis off image; title('Loop-0 randProj 3-dim')
    export_fig(sprintf('%s/id%d_summary.jpg', prefix4saving, idx));

    imgFig3 = figure(2);
    set(imgFig3, 'Position', [100 100 1000 1000]) % [1 1 width height]
    winH = 4;
    winW = 3;
    curWinIdx = 1;
    for loopidx = 1:5:loopNum
        subplot(winH,winW,curWinIdx); curWinIdx = curWinIdx + 1;
        imagesc(predInstanceMaskMat_loops{loopidx,1}); axis off image; title(sprintf('Loop-%d randProj (3-dim)',loopidx))
        subplot(winH,winW,curWinIdx); curWinIdx = curWinIdx + 1;
        imagesc(predInstanceMaskMat_loops{loopidx,2}); axis off image; title(sprintf('Loop-%d randProj (3-dim)',loopidx))
        subplot(winH,winW,curWinIdx); curWinIdx = curWinIdx + 1;
        imagesc(predInstanceMaskMat_loops{loopidx,3}); axis off image; title(sprintf('Loop-%d randProj (3-dim)',loopidx))
    end
    export_fig(sprintf('%s/id%d_summaryLoops.jpg', prefix4saving, idx));
end
%% leaving blank

author: aimerykong


text:
function [PErrorOutput, PErrorOutputmean, PErrorOutputvar] = PError(N,Ptrue,Pestimate)

% This function computes the proportion erros per pixel per endmember
%       given true and estimated proportion values.
%
% SYNTAX : [PErrorOutput, PErrorOutputmean, PErrorOutputvar] = PError(N,Ptrue,Pestimate)
%
% INPUTS
%     N:             Total number of data points. N=N1*N2 in hyperspectral
%                    image data.
%     Ptrue:         True proportion values. N*M, where M is the number of
%                    endmembers.
%     Pestimate:     Estimated proportion values. N*M.
% OUTPUTS
%     PErrorOutput:       PError per pixel
%     PErrorOutputmean:   Mean of PErrorOutput, PError per pixel per
%                         endmember. (This is the result in the Tables in the paper)
%     PErrorOutputvar:    Variance of PErrorOutput across endmembers.
%
% This product is Copyright (c) 2014 University of Missouri
% All rights reserved.
%
% 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.
%   3. Neither the name of the University nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF MISSOURI 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 UNIVERSITY 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.

  M = size(Ptrue,2);
  t =   (Ptrue- Pestimate).^2;
  PErrorOutput = sum(sum(t)) / N; % PError per pixel

PErrorOutputmean = PErrorOutput/M; %Mean of PErrorOutput, PError per pixel per endmember
PErrorOutputvar= var(PErrorOutput); %Variance of PErrorOutput

%%Plot histogram across iterations for one N. FOR DEMO PURPOSES.
% figure(100);
% plot(PErrorOutput);
% xlabel('N');ylabel('PError');
% title('Plot of PErrorOutput');
% figure(200);
% hist(PErrorOutput);
% xlabel('PError');ylabel('# of Data Points');
% title('Histogram of PErrorOutput');

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function menu = MakeNumberedMenu(names) % e.g. name = {Cluster.name} (note {})
menu = [{'(choose)'},names];
for j=2:length(menu),menu(j)={[num2str(j-1) ': ' menu{j}]};end
end

author: xiuyechen

repo: FishExplorer


text:
function [B,elapse] = ITQ_compress(A, model)
%   This is a wrapper function of ITQ testing.
%
%	Usage:
%	[B,elapse] = ITQ_compress(A, model)
%
%	      A: Rows of vectors of data points. Each row is sample point
%     model: The model generated by ITQ_learn.
%
%	      B: The binary code of the input data A. Each row is sample point
%    elapse: The coding time (testing time).
%
%
%
%   version 2.0 --Nov/2016
%   version 1.0 --Jan/2013
%
%   Written by  Yue Lin (linyue29@gmail.com)
%               Deng Cai (dengcai AT gmail DOT com)
%

tmp_T = tic;

Ym = A * model.pc;
Ym = Ym * model.R;
B = (Ym > 0);

elapse = toc(tmp_T);
end

author: ZJULearning

repo: MatlabFunc


text:
function title = p30_title ( )

%*****************************************************************************80
%
%% P30_TITLE returns a title for problem 30.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Branin RCOS Function';

  return
end


text:
% Check that online inference gives same results as filtering for various algorithms

N = 3;
Q = 2;
ss = N*2;

rand('state', 0);
randn('state', 0);

obs_size = 1;
discrete_obs = 0;
bnet = mk_chmm(N, Q, obs_size, discrete_obs);
ns = bnet.node_sizes_slice;

engine = {};
engine{end+1} = hmm_inf_engine(bnet);
E = length(engine);

onodes = (1:N)+N;

T = 4;
ev = cell(ss,T);
ev(onodes,:) = num2cell(randn(N, T));

filter = 1;
loglik2 = zeros(1,E);
for e=1:E
  [engine2{e}, loglik2(e)] = enter_evidence(engine{e}, ev, 'filter', filter);
end

loglik = zeros(1,E);
marg1 = cell(E,N,T);
for e=1:E
  ll = zeros(1,T);
  engine{e} = dbn_init_bel(engine{e});
  for t=1:T
    [engine{e}, ll(t)] = dbn_update_bel(engine{e}, ev(:,t), t);
    for i=1:N
      marg1{e,i,t} = dbn_marginal_from_bel(engine{e}, i);
    end
  end
  loglik1(e) = sum(ll);
end

assert(approxeq(loglik1, loglik2))

a = zeros(E,N,T);
for e=1:E
  for t=1:T
    for i=1:N
      marg2{e,i,t} = marginal_nodes(engine2{e}, i, t);
      a(e,i,t) = (approxeq(marg2{e,i,t}.T(:), marg1{e,i,t}.T(:)));
    end
  end
end

assert(all(a(:)==1))

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function ind = sub2ind2( siz, sub )
% Improved version of sub2ind.
%
% Almost the same as sub2ind, except always returns only a single output
% that contains all the subscript locations.  Also handles multiple linear
% subscripts at the same time more conveniently then matlab's version. See
% help for sub2ind for more info.
%
% USAGE
%  ind = sub2ind2( siz, sub )
%
% INPUTS
%  siz     - size of array into which sub is an index
%  sub     - sub(i,:) is the ith set of subscripts into the array.
%
% OUTPUTS
%  ind     - linear index (or vector of indicies) into given array
%
% EXAMPLE
%  ind = sub2ind2( [10,10], [10 2] )      % 20
%  ind = sub2ind2( [10,10], [9 2; 10 2] ) % 19, 20
%
% See also SUB2IND, IND2SUB2, SUBSTOARRAY
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if(isempty(sub)); ind=[]; return; end;
n = length(siz);

% error check (commented out to speed up substantially)
if( size(sub,2)~=n )
  error('Incorrect dimension for sub');
  % for i = 1:n; if( any( sub(:,i)<1 ) || any( sub(:,i)>siz(i) ) )
  %   error('subscript out of range'); end;
end

k = [1 cumprod(siz(1:end-1))];
ind = 1;
for i = 1:n; ind = ind + (sub(:,i)-1)*k(i); end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function planC = fixStructureZVals(planC)
%Compares zVals stored in structure contours to CT zValues.
%Correct by replacing the structure's zVals with those of
%the CT slice it was contoured on.
%
%We have encountered some archives where the two do not match,
%which causes problems with generating 3D datasets.
%
% JRA 10/24/03
%
% TODO: Vectorize.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};
ctZVal = [planC{indexS.scan}.scanInfo.zValue];

for structNum=1:length(planC{indexS.structures})
    errorFound = 0;
    errorDiff = [];
    contour = planC{indexS.structures}(structNum).contour;
    for sliceNum=1:length(contour)
        for pointsSet=1:length(contour(sliceNum).segments)
            if ~isempty(contour(sliceNum).segments(pointsSet).points)
                structureZVals = contour(sliceNum).segments(pointsSet).points(:,3);

                uniqStrZvals = unique(structureZVals);
                if length(uniqStrZvals) > 1
                    for uniqInd = 1:length(uniqStrZvals)
                        indZvals = find(structureZVals==uniqStrZvals(uniqInd));
                        numInds(uniqInd) = length(indZvals);
                    end
                    [jnk,maxInd] = max(numInds);

                    uniqStructVal = uniqStrZvals(maxInd);
                    CERRStatusString(['Structure "',planC{indexS.structures}(structNum).structureName,'" has inconsistent contours defined at z = ', num2str(ctZVal(ind)), ' slice location.'])
                    CERRStatusString(['Using the contour points defined at z = ',num2str(uniqStructVal)])
                else
                    uniqStructVal = uniqStrZvals;
                end

                %ind = findnearest(ctZVal, unique(structureZVals));
                ind = findnearest(ctZVal, uniqStructVal);

                wrongZVals = structureZVals ~= ctZVal(ind);

                if find(wrongZVals)
                    errorFound = 1;
                    errorDiff = [errorDiff abs(structureZVals' - ctZVal(ind))];
                    contour(sliceNum).segments(pointsSet).points(find(wrongZVals),3) = ctZVal(ind);
                end
            end
        end
    end
    %Now check for bad rasterSegment zValues, that may have been created in
    %old plans from uncorrected contour data.
    %     for sliceNum=1:length(contour)
    [rasterSegs, planC, isError] = getRasterSegments(structNum, planC);
    %rasterSegs = planC{indexS.structures}(structNum).rasterSegments;
    if ~isempty(rasterSegs)
        if rasterSegs(:,1) ~= [planC{indexS.scan}.scanInfo(rasterSegs(:,6)).zValue]';
            rasterSegs(:,1) = [planC{indexS.scan}.scanInfo(rasterSegs(:,6)).zValue]';
        end
    end
    planC{indexS.structures}(structNum).rasterSegments = rasterSegs;
    %     end

    if errorFound
        CERRStatusString(['Correcting inconsistent zValues in ' planC{indexS.structures}(structNum).structureName '. Median correction is ' num2str(median(errorDiff)) ' cm.']);
        planC{indexS.structures}(structNum).contour = contour;
    end
end
CERRStatusString('');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef test_net < matlab.unittest.TestCase

  properties
    num_output
    model_file
    net
  end

  methods (Static)
    function model_file = simple_net_file(num_output)
      model_file = tempname();
      fid = fopen(model_file, 'w');
      fprintf(fid, [ ...
        'name: "testnet" force_backward: true\n' ...
        'layer { type: "DummyData" name: "data" top: "data" top: "label"\n' ...
        'dummy_data_param { num: 5 channels: 2 height: 3 width: 4\n' ...
        '    num: 5 channels: 1 height: 1 width: 1\n' ...
        '    data_filler { type: "gaussian" std: 1 }\n' ...
        '    data_filler { type: "constant" } } }\n' ...
        'layer { type: "Convolution" name: "conv" bottom: "data" top: "conv"\n' ...
        '  convolution_param { num_output: 11 kernel_size: 2 pad: 3\n' ...
        '    weight_filler { type: "gaussian" std: 1 }\n' ...
        '    bias_filler { type: "constant" value: 2 } }\n' ...
        '    param { decay_mult: 1 } param { decay_mult: 0 }\n' ...
        '    }\n' ...
        'layer { type: "InnerProduct" name: "ip" bottom: "conv" top: "ip"\n' ...
        '  inner_product_param { num_output: ' num2str(num_output) ...
        '    weight_filler { type: "gaussian" std: 2.5 }\n' ...
        '    bias_filler { type: "constant" value: -3 } } }\n' ...
        'layer { type: "SoftmaxWithLoss" name: "loss" bottom: "ip" bottom: "label"\n' ...
        '  top: "loss" }' ]);
      fclose(fid);
    end
  end
  methods
    function self = test_net()
      self.num_output = 13;
      self.model_file = caffe.test.test_net.simple_net_file(self.num_output);
      self.net = caffe.Net(self.model_file, 'train');
      % also make sure get_solver runs
      caffe.get_net(self.model_file, 'train');

      % fill in valid labels
      self.net.blobs('label').set_data(randi( ...
        self.num_output - 1, self.net.blobs('label').shape));

      delete(self.model_file);
    end
  end
  methods (Test)
    function self = test_blob(self)
      self.net.blobs('data').set_data(10 * ones(self.net.blobs('data').shape));
      self.verifyEqual(self.net.blobs('data').get_data(), ...
        10 * ones(self.net.blobs('data').shape, 'single'));
      self.net.blobs('data').set_diff(-2 * ones(self.net.blobs('data').shape));
      self.verifyEqual(self.net.blobs('data').get_diff(), ...
        -2 * ones(self.net.blobs('data').shape, 'single'));
      original_shape = self.net.blobs('data').shape;
      self.net.blobs('data').reshape([6 5 4 3 2 1]);
      self.verifyEqual(self.net.blobs('data').shape, [6 5 4 3 2 1]);
      self.net.blobs('data').reshape(original_shape);
      self.net.reshape();
    end
    function self = test_layer(self)
      self.verifyEqual(self.net.params('conv', 1).shape, [2 2 2 11]);
      self.verifyEqual(self.net.layers('conv').params(2).shape, 11);
      self.verifyEqual(self.net.layers('conv').type(), 'Convolution');
    end
    function test_forward_backward(self)
      self.net.forward_prefilled();
      self.net.backward_prefilled();
    end
    function test_inputs_outputs(self)
      self.verifyEqual(self.net.inputs, cell(0, 1))
      self.verifyEqual(self.net.outputs, {'loss'});
    end
    function test_save_and_read(self)
      weights_file = tempname();
      self.net.save(weights_file);
      model_file2 = caffe.test.test_net.simple_net_file(self.num_output);
      net2 = caffe.Net(model_file2, 'train');
      net2.copy_from(weights_file);
      net3 = caffe.Net(model_file2, weights_file, 'train');
      delete(model_file2);
      delete(weights_file);
      for l = 1:length(self.net.layer_vec)
        for i = 1:length(self.net.layer_vec(l).params)
          self.verifyEqual(self.net.layer_vec(l).params(i).get_data(), ...
            net2.layer_vec(l).params(i).get_data());
          self.verifyEqual(self.net.layer_vec(l).params(i).get_data(), ...
            net3.layer_vec(l).params(i).get_data());
        end
      end
    end
  end
end

repo: deepImageAestheticsAnalysis


text:
% qqdiagram() - Empirical quantile-quantile diagram.
%
% Description:
%               The quantiles (percentiles) of the input distribution Y are plotted (Y-axis)
%               against the corresponding quantiles of the input distribution X.
%               If only X is given, the corresponding quantiles are plotted (Y-axis)
%               against the quantiles of a Gaussian distribution ('Normal plot').
%               Two black dots indicate the lower and upper quartiles.
%               If the data in X and Y belong the same distribution the plot will be linear.
%               In this case,the red and black reference lines (.-.-.-.-) will overlap.
%               This will be true also if the data in X and Y belong to two distributions with
%               the same shape, one distribution being rescaled and shifted with respect to the
%               other.
%               If only X is given, a line is plotted to indicate the mean of X, and a segment
%               is plotted to indicate the standard deviation of X. If the data in X are normally
%               distributed, the red and black reference lines (.-.-.-.-) will overlap.
%
% Usage:
%   >>  ah  =  qqdiagram( x, y, pk );
%
% Inputs:
%   x       - vector of observations
%
% Optional inputs:
%   y       - second vector of observation to compare the first to
%   pk      - the empirical quantiles will be estimated at the values in pk [0..1]
%
% Author: Luca Finelli, CNL / Salk Institute - SCCN, 20 August 2002
%
% Reference: Stahel W., Statistische Datenanalyse, Vieweg, Braunschweig/Wiesbaden, 1995
%
% See also:
%   quantile(), signalstat(), eeglab()

% Copyright (C) 2002 Luca Finelli, Salk/SCCN, La Jolla, CA
%
% Reference: Stahel, W. Statistische Datenanalyse, Vieweg, Braunschweig/Wiesbaden 1995
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function qqdiagram( x , y, pk )

if nargin < 1
	help qqdiagram;
	return;
end;

if (nargin == 3 & (any(pk > 1) | any(pk < 0)))
    error('qqdiagram(): elements in pk must be between 0 and 1');
end

if nargin==1
	y=x;
	nn=max(1000,10*length(y))+1;
	x=randn(1,nn);
end

if nargin < 3
	nx=sum(~isnan(x));
	ny=sum(~isnan(y));
   	k=min(nx,ny);
    pk=((1:k) - 0.5) ./ k;  % values to estimate the empirical quantiles at
else
    k=length(pk);
end

if nx==k
    xx=sort(x(~isnan(x)));
else
    xx=quantile(x(~isnan(x)),pk);
end

if ny==k
    yy=sort(y(~isnan(y)));
else
    yy=quantile(y(~isnan(y)),pk);
end

% QQ diagram
plot(xx,yy,'+')
hold on

% x-axis range
maxx=max(xx);
minx=min(xx);
rangex=maxx-minx;
xmin=minx-rangex/50;
xmax=maxx+rangex/50;

% Quartiles
xqrt1=quantile(x,0.25); xqrt3=quantile(x,0.75);
yqrt1=quantile(y,0.25); yqrt3=quantile(y,0.75);

plot([xqrt1 xqrt3],[yqrt1 yqrt3],'k-','LineWidth',2); % IQR range

% Drawing the line
sigma=(yqrt3-yqrt1)/(xqrt3-xqrt1);
cy=(yqrt1 + yqrt3)/2;

if nargin ==1
    maxy=max(y);
    miny=min(y);
    rangey=maxy-miny;
	ymin=miny-rangey/50;
	ymax=maxy+rangey/50;

	plot([(miny-cy)/sigma (maxy-cy)/sigma],[miny maxy],'r-.') % the line
    % For normally distributed data, the slope of the plot line
    % is equal to the ratio of the standard deviation of the distributions
	plot([0 (maxy-mean(y))/std(y)],[mean(y) maxy],'k-.') % the ideal line

	xlim=get(gca,'XLim');
	plot([1 1],[ymin  mean(y)+std(y)],'k--')
	plot([1 1],[mean(y)  mean(y)+std(y)],'k-','LineWidth',2)
        % textx = 1.0;
        % texty = mean(y)+3.0*rangey/50.0;
	% text(double(textx), double(texty),' St. Dev.','horizontalalignment','center')
    set(gca,'xtick',get(gca,'xtick'));  % show that vertical line is at 1 sd
	plot([0 0],[ymin  mean(y)],'k--')
	plot(xlim,[mean(y) mean(y)],'k--')
	% text(double(xlim(1)), double(mean(y)+rangey/50),'Mean X')
	plot([xqrt1  xqrt3],[yqrt1 yqrt3],'k.','MarkerSize',10)
	set(gca,'XLim',[xmin xmax],'YLim',[ymin ymax])
	xlabel('Standard Normal Quantiles')
	ylabel('X Quantiles')
else
    cx=(xqrt1 + xqrt3)/2;
    maxy=cy+sigma*(max(x)-cx);
	miny=cy-sigma*(cx-min(x));

	plot([min(x) max(x)],[miny maxy],'r-.'); % the line
    xlabel('X Quantiles');
    ylabel('Y Quantiles');
end

author: PatternRecognition

repo: OpenBMI


text:
function d = rows(X)
warning('This function is deprecated')

d = size(X,1);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%% QUAD_JOB_ITHACA divides the quadrature problem into tasks.
%
%  Discussion:
%
%    Each task is given a portion of the interval [0,1] over which
%    to integrate.
%
%    Each task returns the estimated integral over its subinterval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUAD_JOB_ITHACA\n' );
  fprintf ( 1, '  Set up and execute a job of independent tasks.\n' );
  fprintf ( 1, '  The computation of an integral over [0,1] is divided\n' );
  fprintf ( 1, '  into tasks integrating from A to B.\n' );
  fprintf ( 1, '  The integral estimate is the sum of the task results.\n' );
%
%  Define the job.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define the job:\n' );

  job = createJob ( ...
    'configuration', 'ithaca_2009b', ...
    'FileDependencies', { 'quad_task.m' } );

  task_num = 4;
  n = 100001;
  n_per_task = 1 + ( n - 1 ) / task_num;
%
%  Define the tasks that constitute the job.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define %d tasks:\n', task_num );

  for task = 1 : task_num
    a = ( task - 1 ) / task_num;
    b =   task       / task_num;
    task_id = createTask ( job, @quad_task, 1, { n_per_task, a, b } );
  end
%
%  Submit the job.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Submit the job:\n' );

  submit ( job );
%
%  Wait for results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wait for job to finish:\n' );

  timestamp ( );
  wait ( job );
  timestamp ( );
%
%  The following commands will capture and print error messages
%  from the tasks, if something bad occurs.
%
  errmsgs = get ( job.Tasks, {'ErrorMessage'} );
  nonempty = ~cellfun ( @isempty, errmsgs );
  celldisp ( errmsgs(nonempty) );
%
%  Retrieve the output from all tasks.
%
  results = getAllOutputArguments ( job );
%
%  The output is a cell array.
%  We much prefer a MATLAB matrix!
%
  results = cell2mat ( results );
%
%  Sum the results to get the estimate.
%
  quad = sum ( results(:,1) );
  exact = pi;
  err = abs ( quad - exact );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Exact integral = %f\n', exact );
  fprintf ( 1, '  Estimate =       %f\n', quad );
  fprintf ( 1, '  Error =          %e\n', err );
%
%  Clean up by freeing memory associated with the job.
%
  destroy ( job );
%
%  Terminate;
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUAD_JOB_ITHACA:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

author: johannesgerer

repo: jburkardt-m


text:
function [varargout] = size2(x, dim, cellflag)

% SIZE2 Size of cells in cell-array.
%
% Use as:
%   [varargout] = size2(x, dim, cellflag)
%
% Input arguments:
%   x   = cell-array
%   dim = dimension argument, default []
%   cellflag = true (default) or false. If false, size2 behaves as size,
%     otherwise the size of each of the cells is returned.

if nargin<3, cellflag = false; end
if nargin<2, dim      = [];    end

if nargout>1 && ~isempty(dim)
  error('Too many output arguments.');
end

if ischar(cellflag) && strcmp(cellflag, 'cell')
  cellflag = true;
end

if ~cellflag
  % traditional functionality
  numdim = ndims(x);

  if isempty(dim)
    siz    = zeros(1,numdim);
    for k = 1:numdim
      switch k
        case 1
          siz(k) = numel(x(:,1,1,1,1,1,1,1,1,1));
        case 2
          siz(k) = numel(x(1,:,1,1,1,1,1,1,1,1));
        case 3
          siz(k) = numel(x(1,1,:,1,1,1,1,1,1,1));
        case 4
          siz(k) = numel(x(1,1,1,:,1,1,1,1,1,1));
        case 5
          siz(k) = numel(x(1,1,1,1,:,1,1,1,1,1));
        case 6
          siz(k) = numel(x(1,1,1,1,1,:,1,1,1,1));
        case 7
          siz(k) = numel(x(1,1,1,1,1,1,:,1,1,1));
        case 8
          siz(k) = numel(x(1,1,1,1,1,1,1,:,1,1));
        case 9
          siz(k) = numel(x(1,1,1,1,1,1,1,1,:,1));
        case 10
          siz(k) = numel(x(1,1,1,1,1,1,1,1,1,:));
        otherwise
      end
    end
  else
    switch dim
      case 1
        siz = numel(x(:,1,1,1,1,1,1,1,1,1));
      case 2
        siz = numel(x(1,:,1,1,1,1,1,1,1,1));
      case 3
        siz = numel(x(1,1,:,1,1,1,1,1,1,1));
      case 4
        siz = numel(x(1,1,1,:,1,1,1,1,1,1));
      case 5
        siz = numel(x(1,1,1,1,:,1,1,1,1,1));
      case 6
        siz = numel(x(1,1,1,1,1,:,1,1,1,1));
      case 7
        siz = numel(x(1,1,1,1,1,1,:,1,1,1));
      case 8
        siz = numel(x(1,1,1,1,1,1,1,:,1,1));
      case 9
        siz = numel(x(1,1,1,1,1,1,1,1,:,1));
      case 10
        siz = numel(x(1,1,1,1,1,1,1,1,1,:));
      otherwise
    end
  end

  if nargout<=1
    varargout{1} = siz;
  elseif nargout>=numdim
    for k = 1:numdim,           varargout{k} = siz(k); end
    for k = (numdim+1):nargout, varargout{k} = 1;      end
  elseif nargout<numdim
    for k = 1:(nargout-1),      varargout{k} = siz(k); end
    varargout{nargout} = prod(siz((k+1):end));
  end
  return;
else
  % operate on the individual cells
%   siz    = ones(numel(x), 10);
%   numdim = zeros(numel(x), 1);
%   for k = 1:numel(x)
%     numdim(k) = ndims(x{k});
%     siz(k,1:numdim(k)) = size(x{k});
%   end
  numdim = cellfun(@ndims, x);
  siz    = cellfun(@size,  x, 'uniformoutput', false);
  if all(numdim==numdim(1))
    siz = reshape(cat(2,siz{:}),numdim(1),[])';
  end

  if isempty(dim)
    if nargout<=1
      varargout{1} = siz(:,1:max(numdim));
    elseif nargout>=max(numdim)
      for k = 1:numdim,           varargout{k} = reshape(siz(:,k), size(x)); end
      for k = (numdim+1):nargout, varargout{k} = ones(size(x));              end
    end
  else
    varargout{1} = reshape(siz(:,dim), size(x));
  end
  return;
end

author: fieldtrip

repo: fieldtrip


text:
function [x,t,intersects]=constrainedLSLineSpher(a,b,alphaVal)
%%CONSTRAINEDLSLINESPHER This function solves the real optimization problem
%       arg min _{x,t} norm(a+b*t-x) such that norm(x)=alphaVal
%       where t is a scalar and x is a vector. That is, it finds the
%       closest point on the unit sphere to a line in an arbitrary number
%       of dimensions. If the line does not intersect the sphere, there
%       will be 1 solution. If the line intersects the sphere, there can be
%       1 or two solutions. All solutions and inputs are real.
%
%INPUTS: a A real xDimX1 vector.
%        b A real xDimX1 vector.
% alphaVal The norm of the solution. If this is omittd or ane mpty matrix
%          is passed, the default of 1 is used.
%
%OUTPUTS: x An xDimX1 or xDimX2 set of solutions.
%         t The t value(s) associated with the solution(s) in x. If b=0,
%           then t will be 0, though its value does not matter.
% intersects If the line intersects the ellipsoid this is true. Otherwise,
%           this is false.
%
%We want to solve arg min _{x,t} norm(a+b*t-x)^2 such that
%norm(x)^2=alphaVal. First, we consider two special cases. The first is if
%b is all zero, in which case t does not matter and we are finding the unit
%vector x such that solves min _{x} norm(a-x). The solution is just
%x=sqrt(alphaVal)*(a/norm(a)). Note that if a=all zeros, then a bunch of
%NaNs will be returned for x --every point on the sphere is the "closest".
%
%The second special case is if the line intersects the unit sphere. We do
%not know a priori if the line does intersect the sphere, so we solve the
%problem and if the solutions are complex, then it does not intersect the
%sphere. In this case, we just substitute the line equation into the
%constraint:
%norm(a+b*t)^2=alphaVal
%This is a quadratic equation and can be solved using the quadratic
%formula.
%
%Finally, if the solutions to the above are complex, we solve the general
%problem. In this instance, the Lagrangian for equality constrained
%optimization is
%L=norm(a+b*t-x)^2+lambda*(norm(x)^2-alphaVal)
%where lambda is the Lagrgian parameter. Taking the gradient with respect
%to x and the derivative with respect to t and setting it equal to 0 leads
%to the following linear system of equations:
%[-b',                   norm(b)^2;
%  (1+lambda)*eye(xDim),        -b] *[x;t] == [-a'*b;a]
%The matrix on the left can be explicitly inverted as
%1/(lambda*norm(b)^2)*[b, (1/(1+lambda)*(b*b'+lambda*norm(b)^2*eye(xDim));
%                       1+lambda,  b'];
%This leads to solutions for t and x as
%t=-a'b/norm(b)^2
%x=(1/(1+lambda))*(a-(1/norm(b)^2)*(a'*b)*b)
%Substituting x into the constraint leads to the solutions for lambda:
%lambda=-1+norm(a-(1/norm(b)^2)*(a'*b)*b)
%lambda=-1-norm(a-(1/norm(b)^2)*(a'*b)*b)
%The solution for lambda to us is the one that minimizes norm(a+b*t-x) for
%the given t.
%
%EXAMPLE:
%In this example, the unit circle is plotted along with some lines. The
%points on the circle closest to the lines are plotted in red and black
%circles are used to mark the location on each line that is closest to the
%circle. Dashed lines connect the closest points on the lines to the
%closest points on the circle when the lines do not intersect the circle.
%Generate a unit circle to plot.
% figure(1)
% clf
% hold on
% %Draw a circle.
% drawEllipse([0;0],eye(2,2),1,'linewidth',2)
% axis([-2,2,-2,2])
% axis square
%
% t=linspace(-10,10,500);
% %Consider 4 lines.
% a=[1,0,0, -2;
%    2,0,-2, -1];
% b=[4,1,-2,3;
%    3,1,-1,-1];
% numLines=size(a,2);
%
% for k=1:numLines
%     xy=bsxfun(@plus,a(:,k),bsxfun(@times,b(:,k),t));
%     plot(xy(1,:),xy(2,:),'linewidth',2)
%     [x,ts]=constrainedLSLineSpher(a(:,k),b(:,k));
%
%     scatter(x(1,:),x(2,:),200,'.r')
%     xySol=bsxfun(@plus,a(:,k),bsxfun(@times,b(:,k),ts));
%     scatter(xySol(1,:),xySol(2,:),100,'ok')
%     if(length(ts)==1)
%         %If the line either grazes the unit circle or does not intercept
%         %the unit circle, connect the point on the line to the nearest
%         %point on the circle.
%         plot([x(1),xySol(1)],[x(2),xySol(2)],'--k')
%     end
% end
%
%October 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

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

%Check for the special case of b being 0. In this instance,
%x=sqrt(alphaVal)*a/norm(a).
if(all(b==0))
    normA=norm(a);
    intersects=normA==1;

    x=sqrt(alphaVal)*(a/normA);
    t=0;%t does not matter.
    return;
end

xDim=size(a,1);

%Deal with the case where the line intersects the unit sphere. This means
%solve for the intersection points and see if they are real. If they are
%complex, then the line does not intersect the unit sphere.
aMag2=dot(a,a);
bMag2=dot(b,b);
ab=dot(a,b);
t=roots([bMag2,2*ab,aMag2-alphaVal]).';
if(~isempty(t)&&isreal(t))
    %The line intersects the sphere of radius alphaVal
    numSol=length(t);
    x=zeros(xDim,numSol);

    for curSol=1:numSol
       x(:,curSol)=a+b*t(curSol);
       %Force normalization (deal with any finite precision issues).
       x(:,curSol)=sqrt(alphaVal)*(x(:,curSol)/norm(x(:,curSol)));
    end
    intersects=true;
    return;
end

%In this instance, the line does not intersect the sphere.
intersects=false;
t=-ab/bMag2;

magVal=norm(a-b*(ab/bMag2))/sqrt(alphaVal);
lambda=zeros(2,1);
lambda(1)=-1+magVal;
lambda(2)=-1-magVal;

%For the return values.
x=zeros(xDim,2);
for k=1:2
    x(:,k)=(1/(1+lambda(k)))*a-(1/(bMag2*(1+lambda(k))))*ab*b;
end

linePt=a+b*t;

if(norm(linePt-x(:,1))<norm(linePt-x(:,2)))
    x=x(:,1);
else
    x=x(:,2);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function t_most_30b_1_1_17(quiet)
%T_MOST_30B_1_1_17  Tests for MOST.

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

if nargin < 1
    quiet = 0;
end

n_tests = 34;

t_begin(n_tests, quiet);

casename = 't_case30_most';
fudging = struct( ...       %% paramters for fudging reserve contract for sopf2
    'fudge',    0.05, ...   %% initial value (MW)
    'step',     0.01, ...   %% if necessary, increase by this amount and retry (MW)
    'lim',      0.1);       %% upper limit (MW), give up if no convergence
                            %% with fudge equal to this limit

%% options
algs.dc     = {'DEFAULT'};  %% opf.dc.solver sequence to try for c3sopf (DC run)
algs.ac     = {'DEFAULT'};  %% opf.ac.solver sequence to try for c3sopf (AC run)
mpopt = mpoption('verbose', 0, 'out.all', 0);
mpopt = mpoption(mpopt, 'opf.violation', 5e-7, 'mips.comptol', 5e-8);
mpopt = mpoption(mpopt, 'sopf.force_Pc_eq_P0', 0);  %% don't constrain contracted == base case dispatch
if have_feature('linprog')
    if have_feature('linprog_ds')
        mpopt = mpoption(mpopt, 'linprog.Algorithm', 'dual-simplex');
    else
        mpopt = mpoption(mpopt, 'linprog.Algorithm', 'simplex');
    end
end
mpoptac = mpoption(mpopt, 'model', 'AC');
mpoptdc = mpoption(mpopt, 'model', 'DC');
mpopt = mpoption(mpopt, 'most.solver', algs.dc{1});

%% turn off warnings
s7 = warning('query', 'MATLAB:nearlySingularMatrix');
s6 = warning('query', 'MATLAB:nearlySingularMatrixUMFPACK');
warning('off', 'MATLAB:nearlySingularMatrix');
warning('off', 'MATLAB:nearlySingularMatrixUMFPACK');

%% define named indices into data matrices
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
    CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
    CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
    CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
    CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
    CT_MODCOST_X] = idx_ct;

%% reserve and delta offers
xgd_table.colnames = {
    'PositiveActiveReservePrice', ...
            'PositiveActiveReserveQuantity', ...
                    'NegativeActiveReservePrice', ...
                            'NegativeActiveReserveQuantity', ...
                                    'PositiveActiveDeltaPrice', ...
                                            'NegativeActiveDeltaPrice', ...
};
xgd_table.data = [
    10.1    15      10.0    15      0.1     0.0;
    10.3    30      10.2    30      0.3     0.2;
    10.5    20      10.4    20      0.5     0.4;
    10.7    25      10.6    25      0.7     0.6;
    20.1    25      20.0    25      60.1    60.0;
    20.3    15      20.2    15      15.1    15.0;
    20.5    30      20.4    30      60.3    60.2;
    20.7    15      20.6    15      15.3    15.2;
    30.1    15      30.0    15      60.3    60.4;
    30.3    30      30.2    30      30.1    30.0;
    30.5    25      30.4    25      60.7    60.6;
    30.7    30      30.6    30      30.3    30.2;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
    0.001   50      0.002   50      0       0;
];

%% contingency table
% label probty  type        row column      chgtype newvalue
contab = [
    1   0.002   CT_TBRCH    1   BR_STATUS   CT_REP  0;      %% line 1-2
    2   0.002   CT_TBRCH    2   BR_STATUS   CT_REP  0;      %% line 1-3, all power from gen1 flows via gen2
    3   0.002   CT_TBRCH    3   BR_STATUS   CT_REP  0;      %% line 2-4, a path to loads @ buses 7 & 8
    4   0.002   CT_TBRCH    5   BR_STATUS   CT_REP  0;      %% line 2-5, a path to loads @ buses 7 & 8
    5   0.002   CT_TBRCH    6   BR_STATUS   CT_REP  0;      %% line 2-6, a path to loads @ buses 7 & 8
    6   0.002   CT_TBRCH    36  BR_STATUS   CT_REP  0;      %% line 28-27, tie line between areas 1 & 3
    7   0.002   CT_TBRCH    15  BR_STATUS   CT_REP  0;      %% line 4-12, tie line between areas 1 & 2
    8   0.002   CT_TBRCH    12  BR_STATUS   CT_REP  0;      %% line 6-10, tie line between areas 1 & 3
    9   0.002   CT_TBRCH    14  BR_STATUS   CT_REP  0;      %% line 9-10, tie line between areas 1 & 3
    10  0.002   CT_TGEN     1   GEN_STATUS  CT_REP  0;      %% gen 1 at bus 1
    11  0.002   CT_TGEN     2   GEN_STATUS  CT_REP  0;      %% gen 2 at bus 2
    12  0.002   CT_TGEN     3   GEN_STATUS  CT_REP  0;      %% gen 3 at bus 22
    13  0.002   CT_TGEN     4   GEN_STATUS  CT_REP  0;      %% gen 4 at bus 27
    14  0.002   CT_TGEN     5   GEN_STATUS  CT_REP  0;      %% gen 5 at bus 23
    15  0.002   CT_TGEN     6   GEN_STATUS  CT_REP  0;      %% gen 6 at bus 13
    20  0.010   CT_TGEN     0   PMIN        CT_REL  1.1;    %% 10% load increase
    20  0.010   CT_TGEN     0   QMIN        CT_REL  1.1;
    21  0.010   CT_TGEN     0   PMIN        CT_REL  0.9;    %% 10% load decrease
    21  0.010   CT_TGEN     0   QMIN        CT_REL  0.9;
];
clist = unique(contab(:, CT_LABEL));
nc = length(clist);

%% load the case
mpc = loadcase(casename);
gbus = mpc.gen(:, GEN_BUS);

%%-----  get c3sopf results  -----
rdc = c3sopf_retry(algs.dc, mpc, xgd_table.data, contab, mpoptdc);
% rac = c3sopf_retry(algs.ac, mpc, xgd_table.data, contab, mpoptac);
% save t_most2_soln rdc rac -v6
% s = load('t_most2_soln');
s.rdc = rdc;
% s.rac = rac;

%%-----  set up data for DC run (most)  -----
ng = size(mpc.gen, 1);      %% number of gens
xgd = loadxgendata(xgd_table, mpc);
md = loadmd(mpc, [], xgd, [], contab);

%%-----  do DC run (most)  -----
r = most(md, mpopt);

%%-----  test the results  -----
t = 'success1';
t_ok(s.rdc.opf_results.success, t);
t = 'success2';
t_is(r.QP.exitflag, 1, 12, t);

t = 'f';
t_is(r.results.f, s.rdc.opf_results.f, 4, t);

t = 'Pg : base';
t_is(r.flow(1,1,1).mpc.gen(:, PG), s.rdc.base.gen(:, PG), 5, t);
t = 'Pg : cont ';
for k = 1:nc
    t_is(r.flow(1,1,k+1).mpc.gen(:, PG), s.rdc.cont(k).gen(:, PG), 5, sprintf('%s %d', t, k));
end

% t = 'gen : base';
% t_is(r.flow(1,1,1).mpc.gen(:,1:MU_PMIN), s.rdc.base.gen(:,1:MU_PMIN), 3, t);
% t = 'gen : cont ';
% for k = 1:nc
%     t_is(r.flow(1,1,k+1).mpc.gen(:,1:MU_PMIN), s.rdc.cont(k).gen(:,1:MU_PMIN), 3, sprintf('%s %d', t, k));
% end

t = 'energy prices';
t_is(r.results.GenPrices, s.rdc.energy.prc.sum_bus_lam_p(gbus), 6, t);

t = 'Pc';
t_is(r.results.Pc, s.rdc.energy.Pc, 4, t);

t = 'Gmin';
t_is(r.results.Pc - r.results.Rpm, s.rdc.energy.Gmin, 4, t);

t = 'Gmax';
t_is(r.results.Pc + r.results.Rpp, s.rdc.energy.Gmax, 4, t);

t = 'upward contingency reserve quantities';
t_is(r.results.Rpp, s.rdc.reserve.qty.Rp_pos, 4, t);

t = 'downward contingency reserve quantities';
t_is(r.results.Rpm, s.rdc.reserve.qty.Rp_neg, 4, t);

t = 'upward contingency reserve prices';
t_is(r.results.RppPrices, s.rdc.reserve.prc.Rp_pos, 6, t);

t = 'downward contingency reserve prices';
t_is(r.results.RpmPrices, s.rdc.reserve.prc.Rp_neg, 6, t);

t = 'contingency physical ramp price';
[vv, ll] = r.om.get_idx();
Ramp_P_max = zeros(ng, nc);
sum_muPmax = zeros(ng, 1);
sum_muPmin = zeros(ng, 1);
for k = 1:nc+1
    ii = find(r.flow(1,1,k).mpc.gen(:, GEN_STATUS) > 0);
    if k > 1
        Ramp_P_max(ii,k-1) = (r.QP.lambda.mu_u(ll.i1.rampcont(1,1,k):ll.iN.rampcont(1,1,k)) - r.QP.lambda.mu_l(ll.i1.rampcont(1,1,k):ll.iN.rampcont(1,1,k))) / mpc.baseMVA;
    end
    sum_muPmax(ii) = sum_muPmax(ii) + r.flow(1,1,k).mpc.gen(ii, MU_PMAX);
    sum_muPmin(ii) = sum_muPmin(ii) + r.flow(1,1,k).mpc.gen(ii, MU_PMIN);
end
t_is(Ramp_P_max, s.rdc.energy.mu.Ramp_P_max, 2, t);

t = 'sum_muPmax';
t_is(sum_muPmax, s.rdc.energy.sum_muPmax, 1, t);

t = 'sum_muPmin';
t_is(sum_muPmin, s.rdc.energy.sum_muPmin, 0.9, t);

t = 'Rpmax_pos';
Rpmax_pos = (r.QP.lambda.upper(vv.i1.Rpp(1):vv.iN.Rpp(1)) - r.QP.lambda.lower(vv.i1.Rpp(1):vv.iN.Rpp(1))) / mpc.baseMVA;
t_is(Rpmax_pos, s.rdc.reserve.mu.Rpmax_pos, 6, t);

t = 'Rpmax_neg';
Rpmax_neg = (r.QP.lambda.upper(vv.i1.Rpm(1):vv.iN.Rpm(1)) - r.QP.lambda.lower(vv.i1.Rpm(1):vv.iN.Rpm(1))) / mpc.baseMVA;
t_is(Rpmax_neg, s.rdc.reserve.mu.Rpmax_neg, 6, t);

% g1 = s.rdc.base.gen(:, PG);
% g2 = r.flow(1,1,1).mpc.gen(:, PG);
% for k = 1:nc
%     g1 = [ g1 s.rdc.cont(k).gen(:, PG) ];
%     g2 = [ g2 r.flow(1,1,k+1).mpc.gen(:, PG) ];
% end
% [m,n] = size(g1);
% for j = 1:n
%     fprintf('\n');
%     for i = 1:m
%         fprintf('%9.2f  %9.2f\n', g1(i,j), g2(i,j));
%     end
% end

%%-----  do AC run (most)  -----
%mostac;

%% turn warnings back on
warning(s7.state, 'MATLAB:nearlySingularMatrix');
warning(s6.state, 'MATLAB:nearlySingularMatrixUMFPACK');

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% pop_rejkurt() - rejection of artifact in a dataset using kurtosis
%                 of activity (i.e. to detect peaky distribution of
%                 activity).
%
% Usage:
%   >> pop_rejkurt( INEEG, typerej) % pop-up interative window mode
%   >> [OUTEEG, locthresh, globthresh, nrej] = ...
%		= pop_rejkurt( INEEG, typerej, elec_comp, ...
%                   locthresh, globthresh, superpose, reject, vistype);
%
% Graphical interface:
%   "Electrode" - [edit box] electrodes or components (number) to take into
%                 consideration for rejection. Same as the 'elec_comp'
%                 parameter from the command line.
%   "Single-channel limit" - [edit box] pkurtosis limit in terms of
%                 standard-dev. Same as 'locthresh' command line
%                 parameter.
%   "All-channel limit" - [edit box] kurtosis limit in terms of
%                 standard-dev (all channel regrouped). Same as
%                 'globthresh' command line parameter.
%   "Display with previous rejection" - [edit box] can be either YES or
%                 NO. This edit box corresponds to the command line input
%                 option 'superpose'.
%   "Reject marked trials" - [edit box] can be either YES or NO. This edit
%                 box corresponds to the command line input option 'reject'.
%   "visualization type" - [edit box] can be either REJECTRIALS or EEGPLOT.
%                 This edit box corresponds to the command line input
%                 option 'vistype'.
%
% Inputs:
%   INEEG      - input dataset
%   typerej    - type of rejection (0 = independent components; 1 = eeg
%              data). Default is 1. For independent components, before
%              thresholding, the activity is normalized for each
%              component.
%   elec_comp  - [e1 e2 ...] electrodes (number) to take into
%              consideration for rejection
%   locthresh  - activity kurtosis limit in terms of standard-dev.
%   globthresh - global limit (for all channel). Same unit as above.
%   superpose  - 0=do not superpose pre-labelling with previous
%              pre-labelling (stored in the dataset). 1=consider both
%              pre-labelling (using different colors). Default is 0.
%   reject     - 0=do not reject labelled trials (but still store the
%              labels. 1=reject labelled trials. Default is 1.
%   vistype    - visualization type. 0 calls rejstatepoch() and 1 calls eegplot()
%              default is 0.
%
% Outputs:
%   OUTEEG     - output dataset with updated kurtosis array
%   locthresh  - electrodes probability of activity thresholds in terms
%              of standard-dev.
%   globthresh - global threshold (where all electrode activity are
%              regrouped).
%   nrej       - number of rejected sweeps
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: rejkurt(), rejstatepoch(), pop_rejepoch(), eegplot(), eeglab()

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

% 01-25-02 reformated help & license -ad
% 03-07-02 added srate argument to eegplot call -ad
% 03-08-02 add eeglab options -ad

function [EEG, locthresh, globthresh, nrej, com] = pop_rejkurt( EEG, icacomp, elecrange, ...
                       		locthresh, globthresh, superpose, reject, vistype, topcommand);
com = '';
if nargin < 1
   help pop_rejkurt;
   return;
end;
if nargin < 2
   icacomp = 1;
end;
if exist('reject') ~= 1
    reject = 1;
end;
if icacomp == 0
	if isempty( EEG.icasphere )
	    ButtonName=questdlg( 'Do you want to run ICA now ?', ...
                         'Confirmation', 'NO', 'YES', 'YES');
    	switch ButtonName,
        	case 'NO', disp('Operation cancelled'); return;
        	case 'YES', [ EEG com ] = pop_runica(EEG);
    	end % switch
	end;
end;

if nargin < 3

	% which set to save
	% -----------------
	promptstr   = { [ fastif(icacomp, 'Electrode', 'Component') ' (number(s); Ex: 2 4 5):' ], ...
					[ fastif(icacomp, 'Single-channel', 'Single-component') ' limit(s) (std. dev(s).: Ex: 2 2 2.5):'], ...
					[ fastif(icacomp, 'All-channel', 'All-component') ' limit (std. dev(s).: Ex: 2.1 2 2):'], ...
               		'Display with previously marked rejections? (YES or NO)', ...
         			'Reject marked trial(s)? (YES or NO)', ...
         			'Visualization mode (REJECTTRIALS or EEGPLOT)' };
	inistr      = { fastif(icacomp, ['1:' int2str(EEG.nbchan)], ['1:' int2str(size(EEG.icaweights,1))])...
					fastif(icacomp, '3', '5'),  ...
					fastif(icacomp, '3', '5'), ...
               		'YES', ...
            		'NO', ...
            		'REJECTTRIALS' };

	result       = inputdlg2( promptstr, fastif(~icacomp, 'Trial rejection using comp. kurtosis -- pop_rejkurt()', 'Trial rejection using data kurtosis -- pop_rejkurt()'), 1,  inistr, 'pop_rejkurt');
	size_result  = size( result );
	if size_result(1) == 0 return; end;
	elecrange    = result{1};
	locthresh    = result{2};
	globthresh   = result{3};
	switch lower(result{4}), case 'yes', superpose=1; otherwise, superpose=0; end;
	switch lower(result{5}), case 'yes', reject=1; otherwise, reject=0; end;
	switch lower(result{6}), case 'rejecttrials', vistype=0; otherwise, vistype=1; end;
end;

if ~exist('vistype') vistype = 0; end;
if ~exist('reject') reject = 0; end;
if ~exist('superpose') superpose = 1; end;

if isstr(elecrange) % convert arguments if they are in text format
	calldisp = 1;
	elecrange = eval( [ '[' elecrange ']' ]  );
	locthresh = eval( [ '[' locthresh ']' ]  );
	globthresh = eval( [ '[' globthresh ']' ]  );
else
	calldisp = 0;
end;

if isempty(elecrange)
	error('No electrode selectionned');
end;

% compute the joint probability
% -----------------------------
if icacomp == 1
	fprintf('Computing kurtosis for channels...\n');
    if isempty(EEG.stats.kurtE )
		[ EEG.stats.kurtE rejE ] = rejkurt( EEG.data, locthresh, EEG.stats.kurtE, 1);
	end;
	[ tmp rejEtmp ] = rejkurt( EEG.data(elecrange, :,:), locthresh, EEG.stats.kurtE(elecrange, :), 1);
    rejE    = zeros(EEG.nbchan, size(rejEtmp,2));
	rejE(elecrange,:) = rejEtmp;

	fprintf('Computing all-channel kurtosis...\n');
	tmpdata = permute(EEG.data, [3 1 2]);
	tmpdata = reshape(tmpdata, size(tmpdata,1), size(tmpdata,2)*size(tmpdata,3));
	[ EEG.stats.kurt rej ] = rejkurt( tmpdata, globthresh, EEG.stats.kurt, 1);
else
	fprintf('Computing joint probability for components...\n');
    % test if ICA was computed
    % ------------------------
    icaacttmp = eeg_getica(EEG);
    if isempty(EEG.stats.icakurtE )
		[ EEG.stats.icakurtE rejE ] = rejkurt( icaacttmp, locthresh, EEG.stats.icakurtE, 1);
	end;
	[ tmp rejEtmp ] = rejkurt( icaacttmp(elecrange, :,:), locthresh, EEG.stats.icakurtE(elecrange, :), 1);
	rejE    = zeros(size(icaacttmp,1), size(rejEtmp,2));
	rejE(elecrange,:) = rejEtmp;

	fprintf('Computing global joint probability...\n');
	tmpdata = permute(icaacttmp, [3 1 2]);
	tmpdata = reshape(tmpdata, size(tmpdata,1), size(tmpdata,2)*size(tmpdata,3));
	[ EEG.stats.icakurt rej] = rejkurt( tmpdata, globthresh, EEG.stats.icakurt, 1);
end;
rej = rej' | max(rejE, [], 1);
fprintf('%d/%d trials marked for rejection\n', sum(rej), EEG.trials);

if calldisp
	if vistype == 1 % EEGPLOT -------------------------
	    if icacomp == 1 macrorej  = 'EEG.reject.rejkurt';
	        			macrorejE = 'EEG.reject.rejkurtE';
	    else			macrorej  = 'EEG.reject.icarejkurt';
	        			macrorejE = 'EEG.reject.icarejkurtE';
	    end;

		colrej = EEG.reject.rejkurtcol;
		eeg_rejmacro; % script macro for generating command and old rejection arrays

	    if icacomp == 1
	        eegplot( EEG.data(elecrange,:,:), 'srate', ...
		      EEG.srate, 'limits', [EEG.xmin EEG.xmax]*1000 , 'command', command, eegplotoptions{:});
	    else
	        eegplot( icaacttmp(elecrange,:,:), 'srate', ...
		      EEG.srate, 'limits', [EEG.xmin EEG.xmax]*1000 , 'command', command, eegplotoptions{:});
	    end;
    else % REJECTRIALS -------------------------
	  	if icacomp	== 1
			[ rej, rejE, n, locthresh, globthresh] = ...
				rejstatepoch( EEG.data, EEG.stats.kurtE(elecrange,:), 'global', 'on', 'rejglob', EEG.stats.kurt, ...
						'threshold', locthresh, 'thresholdg', globthresh, 'normalize', 'off'  );
		else
			[ rej, rejE, n, locthresh, globthresh] = ...
				rejstatepoch( icaacttmp, EEG.stats.icakurtE(elecrange,:), 'global', 'on', 'rejglob', EEG.stats.icakurt, ...
						'threshold', locthresh, 'thresholdg', globthresh, 'normalize', 'off' );
		end;
		nrej = n;
	end;
else
	% compute rejection locally
	rejtmp = max(rejE(elecrange,:),[],1);
	rej = rejtmp | rej;
	nrej =  sum(rej);
	fprintf('%d trials marked for rejection\n', nrej);
end;
if ~isempty(rej)
	if icacomp	== 1
		EEG.reject.rejkurt = rej;
		EEG.reject.rejkurtE = rejE;
	else
		EEG.reject.icarejkurt = rej;
		EEG.reject.icarejkurtE = rejE;
	end;
    if reject
        EEG = pop_rejepoch(EEG, rej, 0);
    end;
end;
nrej = sum(rej);

com = [ com sprintf('%s = pop_rejkurt(%s,%s);', inputname(1), ...
		inputname(1), vararg2str({icacomp,elecrange,locthresh,globthresh,superpose,reject})) ];
if nargin < 3 & nargout == 2
	locthresh = com;
end;

return;


text:
%CODEGENERATOR.GENGRAVLOAD Generate code for gravitational load
%
%  G = cGen.gengravload() is a symbolic vector (1xN) of joint load
%  forces/torques due to gravity.
%
% Notes::
% - Side effects of execution depends on the cGen flags:
%   - saveresult: the symbolic expressions are saved to
%     disk in the directory specified by cGen.sympath
%   - genmfun: ready to use m-functions are generated and
%     provided via a subclass of SerialLink stored in cGen.robjpath
%   - genslblock: a Simulink block is generated and stored in a
%     robot specific block library cGen.slib in the directory
%     cGen.basepath
%   - genccode: generates C-functions and -headers in the directory
%     specified by the ccodepath property of the CodeGenerator object.
%   - mex: generates robot specific MEX-functions as replacement for the
%     m-functions mentioned above. Access is provided by the SerialLink
%     subclass. The MEX files rely on the C code generated before.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator, CodeGenerator.geninvdyn, CodeGenerator.genfdyn.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function [G] = gengravload(CGen)

%% Derivation of symbolic expressions
CGen.logmsg([datestr(now),'\tDeriving gravitational load vector']);

q = CGen.rob.gencoords;
tmpRob = CGen.rob.nofriction;
G = tmpRob.gravload(q);

CGen.logmsg('\t%s\n',' done!');

%% Save symbolic expressions
if CGen.saveresult
    CGen.logmsg([datestr(now),'\tSaving symbolic expression for gravitational load']);

    CGen.savesym(G,'gravload','gravload.mat');

    CGen.logmsg('\t%s\n',' done!');
end

% M-Functions
if CGen.genmfun
    CGen.genmfungravload;
end

% Embedded Matlab Function Simulink blocks
if CGen.genslblock
    CGen.genslblockgravload;
end

%% C-Code
if CGen.genccode
    CGen.genccodegravload;
end

%% MEX
if CGen.genmex
    CGen.genmexgravload;
end

end

author: petercorke

repo: robotics-toolbox-matlab


text:
function F = transpose( F )
% .' transpose of a DISKFUN
%    F = TRANSPOSE(F) is the same as F.'. The orientation of the DISKFUNV
%    is transposed.
%
%    See also DISKFUNV/CTRANSPOSE

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

F.isTransposed = ~F.isTransposed;

end

author: chebfun

repo: chebfun