text:
%% ModalSrc
% Concrete subclass of <Source.html |Source|> representing an electric dipole
% distribution that generates a mode of a waveguide.

%%% Description
% |ModalSrc| is used to store the transverse _J_-field distribution on a plane
% to generate a mode of a waveguide.  The plane should be orthogonal to the axis
% of the waveguide.  The mode is calculated by an external mode solver.  To
% assist the mode solver to calculate the mode, users need to provide an
% estimate of the effective refractive constant of the mode to the constructor
% of |ModalSrc|.

%%% Construction
%  src = ModalSrc(normal_axis, intercept, opts)
%  src = ModalSrc(normal_axis, intercept, opts, KA)
%
% *Input Arguments*
%
% * |normal_axis|: direction normal to the plane over which the electric dipoles
% distribute.  It is also the axis of the waveguide.  It should be one of
% |Axis.x|, |Axis.y|, |Axis.z|.
% * |intercept|: location of the plane in the |normal_axis| direction.
% * |opts|: options structure that controls the behavior of the mode solver.
% See below for more details.
% * |KA|: integral of the norm of the surface current density.  It indicates the
% strength of the mode amplitude.  If unassigned, the deault value |KA = 1| is
% used.
%
% *Options Structure*
%
% |opts| should have a string property |opts.clue| that indicates the type of
% the clue that the mode solver utilizes to find an appropriate waveguide mode.
% The clue can be either the order number of the waveguide mode, or the
% effective refractive index of the waveguide mode.
%
% * |opts.clue = 'order'| commands the mode solver to find the |n|th-order mode,
% where |n| is the order number you should provide as |opts.order = n|.  For
% example, |opts| with |opts.clue = 'order'| and |opts.order = 1| commands the
% mode solver to find the fundamental mode.
% * |opts.clue = 'guess'| commands the mode solver to find the mode with the
% effective index close to a "guess".  The guess of the effective index should
% be provided as |opts.neff = neff_guess|.  For example, |opts| with |opts.clue
% = 'guess'| and |opts.neff = 2.0 - 0.01i| commands the mode solver to find the
% mode whose effective index is close to |2.0 - 0.01i|.
%
% In general, |opts.clue = 'guess'| finds the mode much faster than |opts.clue =
% 'order'|.

%%% Note
% In the finite-difference grid, |ModalSrc| excites dipoles at the _E_-field
% points.  This poses a condition on |intercept| argument in the constructor:
% |intercept| should be at a dual grid point in the |normal_axis| direction.
% Therefore, make sure that |intercept| does not overlap with the locations of
% the vertices of <Shape.html |Shape|> in the |normal_axis| direction; otherwise
% dynamic grid generation in <moxwell_run.html |maxwell_run|> will fail.

%%% Example
%   % Create an instance of ModalSrc.
%   modeopts.clue = 'order';
%   modeopts.order = 1;
%   src =  ModalSrc(Axis.y, -1000, modeopts);  % y = -1000 should not be primary grid point
%
%   % Use the constructed src in maxwell_run().
%   [E, H] = maxwell_run({INITIAL ARGUMENTS}, 'SRCJ', src);

%%% See Also
% <Plane.html |PlaneSrc|>, <TFSFPlaneSrc.html |TFSFPlaneSrc|>, <ModalSrc.html
% |ModalSrc|>, <maxwell_run.html |maxwell_run|>

classdef ModalSrc < Source

	properties (SetAccess = immutable)
		normal_axis  % plane normal axis: one of Axis.x, Axis.y, Axis.z
		intercept  % intercept between plane and normal axis
		KA  % surface integral of surface current density K; for z-normal plane, volume integral of (|Jx|+|Jy|)
		opts  % options to control method to calculate mode
	end

	properties (SetAccess = private)
		grid2d  % instance of Grid2d
		osc  % instance of Oscillation
		E2d  % {Ex2d Ey2d Ez2d}: cell array of Scalar2d for E on this plane
		H2d  % {Hx2d Hy2d Hz2d}: cell array of Scalar2d for H on this plane
		JMh  % 2D array J or M in horizontal direction on this plane: e.g., Jx for normal == z
		JMv  % 2D array J or M in vertical direction on this plane: e.g., Jy for normal == z
		neff  % effective n
	end

	% Properties for dispersion relation and propagation length
	properties (Dependent, SetAccess = immutable)
		ispreped  % true if this ModalSrc is prepared; false otherwise
		beta_r  % real part of complex wavevector
		Lp  % propagation length
	end

	methods
		function this = ModalSrc(normal_axis, intercept, opts, KA)
			chkarg(istypesizeof(normal_axis, 'Axis'), ...
				'"normal_axis" should be instance of Axis.');
			chkarg(istypesizeof(intercept, 'real'), '"intercept" should be real.');

			if nargin < 3  % no opts
				opts.clue = 'order';
				opts.order = 1;  % fundamental mode
			end
			chkarg(istypesizeof(opts, 'struct'), '"opts" should be structure.');

			chkarg(isequal(opts.clue, 'guess') || isequal(opts.clue, 'order'), ...
				'"opts.clue" should be ''guess'' or ''order'' (string).');
			if isequal(opts.clue, 'guess')
				chkarg(isfield(opts, 'neff') && istypesizeof(opts.neff, 'complex'), ...
					'"opts.neff" should be complex.');
				% if opts.clue == 'guess', opts can have additional field
				% opts.H2d.
			else  % opts.clue == 'order'
				if ~isfield(opts, 'order')
					opts.order = 1;  % fundamental mode
				end
				chkarg(istypesizeof(opts.order, 'int') && opts.order > 0, ...
					'"opts.order" should be positive integer.');
			end

			if nargin < 4  % no KA
				KA = 1.0;
			end
			chkarg(istypesizeof(KA, 'real'), '"KA" should be real.');

			lgrid = cell(1, Axis.count);
			laltgrid = cell(1, Axis.count);
			lgrid{normal_axis} = intercept;
			plane = Plane(normal_axis, intercept);
			this = this@Source(lgrid, laltgrid, plane);

			this.normal_axis = normal_axis;
			this.intercept = intercept;
			this.opts = opts;
			this.KA = KA;
			this.JMh = [];
			this.JMv = [];
		end

		function beta_r = get.beta_r(this)
			beta = 2*pi*this.neff / this.osc.in_L0();
			beta_r = real(beta);
		end

		function Lp = get.Lp(this)
			beta = 2*pi*this.neff / this.osc.in_L0();
			Lp = -1/imag(beta);
		end

		function ispreped = get.ispreped(this)
			ispreped = ~isempty(this.JMh) && ~isempty(this.JMh);
		end

		function setEH(this, neff, osc, E_cell, H_cell, ge, grid3d)
			chkarg(istypesizeof(neff, 'complex'), '"neff" should be complex.');
			this.neff = neff;

			chkarg(istypesizeof(osc, 'Oscillation'), '"osc" should be instance of Oscillation.');
			this.osc = osc;

			chkarg(istypesizeof(grid3d, 'Grid3d'), '"grid3d" should be instance of Grid3d.');
			this.grid2d = Grid2d(grid3d, this.normal_axis);

			chkarg(istypesizeof(ge, 'GT'), '"ge" should be instance of GT.');

			Nh = this.grid2d.N(Dir.h);
			Nv = this.grid2d.N(Dir.v);
			h = this.grid2d.axis(Dir.h);
			v = this.grid2d.axis(Dir.v);

			assert(istypesizeof(E_cell, 'complexcell', [1 Axis.count], [Nh Nv]), ...
				'"E_cell" should be length-%d row cell array whose each element is %d-by-%d array with complex elements.', Axis.count, Nh, Nv);
			assert(istypesizeof(H_cell, 'complexcell', [1 Axis.count], [Nh Nv]), ...
				'"H_cell" should be length-%d row cell array whole each element is %d-by-%d array with complex elements.', Axis.count, Nh, Nv);

			if this.gt == ge  % source is J
				this.JMh = H_cell{v};
				this.JMv = -H_cell{h};
			else  % source is M
				this.JMh = E_cell{v};
				this.JMv = -E_cell{h};
			end

			this.E2d = cell(1, Axis.count);
			this.H2d = cell(1, Axis.count);
			for w = Axis.elems
				this.E2d{w} = array2scalar(E_cell{w}, this.grid2d, ge, FT.e, w, osc, PhysQ.E, this.intercept);
				this.H2d{w} = array2scalar(H_cell{w}, this.grid2d, ge, FT.h, w, osc, PhysQ.H, this.intercept);
			end
		end

		function [index_cell, JMw_patch] = generate_kernel(this, w_axis, grid3d)
			assert(~isempty(this.JMh) && ~isempty(this.JMv), '"JMh" and "JMv" are not set in this ModalSrc.');
			if w_axis == this.normal_axis
				JMw_patch = [];
				index_cell = cell(1, Axis.count);
			else
				g2d = Grid2d(grid3d, this.normal_axis);
				assert(isequal(g2d, this.grid2d), ...
					'%s-normal cross section of "grid3d" is different from the one set with JMh and JMv.', char(this.normal_axis));

				h = this.grid2d.axis(Dir.h);
				v = this.grid2d.axis(Dir.v);
				n = this.normal_axis;

				g = this.gt;
				ind_n = ind_for_loc(this.intercept, n, g, grid3d);

				% Set index_cell.
				index_cell = {':', ':', ':'};
				index_cell{n} = ind_n;

				% Set Jw_patch.
				dn = grid3d.dl{n,g}(ind_n);
				dha = grid3d.dl{h, alter(this.gt)};
				dhg = grid3d.dl{h, this.gt};
				dva = grid3d.dl{v, alter(this.gt)};
				dvg = grid3d.dl{v, this.gt};

				dVh = dn .* (dha.' * dvg);
				dVv = dn .* (dhg.' * dva);

				KA_curr = abs(this.JMh) .* dVh + abs(this.JMv) .* dVv;
				KA_curr = sum(KA_curr(:));  % KA_curr is real
				JMhv = [this.JMh(:); this.JMv(:)];
				[~, i_pf] = max(abs(JMhv));
				phasefactor = JMhv(i_pf)/abs(JMhv(i_pf));
				KA_curr = KA_curr * phasefactor;  % KA_curr is complex
				norm_factor = this.KA/KA_curr;  % normalization factor

				if w_axis == h
					JMw_patch = norm_factor .* this.JMh;
				else
					assert(w_axis == v);
					JMw_patch = norm_factor .* this.JMv;
				end

% 				if h > v  % h == Axis.z and v == Axis.x if this.normal_axis == Axis.y
% 					Jw_patch = permute(Jw_patch, int([Dir.v, Dir.h]));
% 				end
				JMw_patch = ipermute(JMw_patch, int([h v n]));
			end
		end
	end

end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function x = chebyshev1 ( n )

%*****************************************************************************80
%
%% CHEBYSHEV1 returns the Type 1 Chebyshev points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2018
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), the points.
%
  x(1:n) = cos ( ( 2.0 * (0:n-1) + 1.0 ) * pi / 2 / n );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = mrdivide(X,Y)
%MRDIVIDE (overloaded)

if (isa(Y,'sdpvar'))
    if Y.dim(1)*Y.dim(2) == 1
        y = X*Y^-1;
        return
    else
        error('Division of matrix variables not possible (maybe you meant ./).')
    end
end

try

    % Quick exit on scalar division
    if X.dim(1)*X.dim(2) == 1 & length(Y)==1
        if isinf(Y)
            y = 0;
        else
            y = X;
            y.basis = y.basis/Y;
        end
        return
    end

    lmi_variables = getvariables(X);
    nv = length(lmi_variables);
    y  = X;
    n = X.dim(1);
    m = X.dim(2);
    if (n==1) & (m==1) % SPECIAL CODE FOR FREAKY MATLAB BUG
        y.basis = sparse(reshape(reshape(full(X.basis(:,1)),n,m)/Y,n*m,1));
        for i = 1:nv
            temp = reshape(full(X.basis(:,i+1)),n,m)/Y;
            y.basis(:,i+1) = sparse(temp(:));
        end;
    else
        y.basis = reshape(reshape(X.basis(:,1),n,m)/Y,n*m,1);
        for i = 1:nv
            temp = reshape(X.basis(:,i+1),n,m)/Y;
            y.basis(:,i+1) = temp(:);
        end;
    end

    y.dim(1) = size(temp,1);
    y.dim(2) = size(temp,2);
catch
    error(lasterr);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function d = p03_fd ( p )

%*****************************************************************************80
%
%% P03_FD is a signed distance function for problem 3.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    06 February 2006
%
%  Parameters:
%
%    Input, real P, one or more points.
%
%    Output, real D, the signed distance of each point to the boundary of the region.
%
  d1 = drectangle ( p, -1.0, 1.0, -1.0, 1.0 );
  d2 = dcircle ( p, 0.0, 0.0, 0.4 );

  d = ddiff ( d1, d2 );

  return
end


text:
function [uniqueSparse, grayGraph] = Mesh2Connections(mesh)
%
%  [uniqueSparse,grayGraph] = Mesh2Connections(mesh)
%
% AUTHOR:  Maher Khoury
% DATE:    08.19.99
% PURPOSE:
%  Takes in a Mesh structure and computes a connections matrix.
%
%  The Mesh structure has data in a striplist and triangles.
%  These are OpenGL constructs.  A striplist is organized so that
%  the first three nodes defines triangle, the 2-4 define another
%  triangle, 3-5, and so forth.
%
%
% INPUT:
%   Mesh
%
% SEE ALSO:  MrReadMrM to understand more about the Mesh structure.
%
% 09.03.99  BW.
%    Rounding the mesh early caused problems.  So, I leave the
%    vertices as floats.  The rounding only happens now when the
%    data are converted to nodes in Connections2Graph.
%

%    vertices = round(mesh.vertices); -- MK rounded.  Why?  I
%    removed this rounding from here.
%  There is also an issue of when the vertices should have 1
% added to them.  This should probably happen in Connections2Graph?
%
vertices = mesh.vertices;

%% Get number of nodes
%
numNodes = size(vertices,1);

%% Initialize nodeConnection sparse matrix
%
nodeConnections = sparse(numNodes,numNodes);

%% Build unique list.  nodeIdx maps the unique node entries back
%  to the redundant node list.  uniqueIdx are the nodes in
%  vertices that are instances of unique node.
%    uniqueVertices = vertices(nodeIdx)
%    vertices    = uniqueVertices(uniqueIdx)
%
[uniqueVertices, nodeIdx, uniqueIdx] = unique(vertices,'rows');

%% Fill the sparse matrix with stripList
%
h = mrvWaitbar(0,'Reading Edges from Stripped List');
for k=1:size(mesh.stripList,1),
  % Start at the 3rd vertex on the list, and build your triangles
  % looking backwards.
  %
   vertex = mesh.stripList(k,1)+3;
   nodeConnections(vertex,vertex-1)=1;
   nodeConnections(vertex,vertex-2)=1;
   nodeConnections(vertex-2,vertex-1)=1;
   nodeConnections(vertex-1,vertex)=1;
   nodeConnections(vertex-2,vertex)=1;
   nodeConnections(vertex-1,vertex-2)=1;

   % Go on to the next vertex if you are between the current one
   % and the next strip list location.
   %
   while vertex < mesh.stripList(k,1)+mesh.stripList(k,2),
      vertex = vertex + 1;
      nodeConnections(vertex,vertex-1)=1;
      nodeConnections(vertex,vertex-2)=1;
      nodeConnections(vertex-1,vertex)=1;
      nodeConnections(vertex-2,vertex)=1;

      if ~mod(vertex,100)
	mrvWaitbar(vertex/(mesh.triangleOffset-1));
      end
   end
end

close(h);

%% Check statement
%
if ~((vertex+1)==mesh.triangleOffset)
   disp('Problem !!');
   return;
end

ptr = mesh.triangleOffset;

%% Fill the sparse matrix with triangle
%
h = mrvWaitbar(0,'Reading Edges from Triangles');
while ptr<size(vertices,1),
   nodeConnections(ptr,ptr+1)=1;nodeConnections(ptr,ptr+2)=1;nodeConnections(ptr+1,ptr+2)=1;
   nodeConnections(ptr+1,ptr)=1;nodeConnections(ptr+2,ptr)=1;nodeConnections(ptr+2,ptr+1)=1;
   ptr = ptr+3;
   if ~mod(ptr,100)
     mrvWaitbar(ptr/numNodes);
   end
end
close(h);

%% So now we have an unconnected set of strips and triangles in the nodeConnections structure
%  We go through the unique list, find all the strips and triangles the nodes in that list
%  belong to and we merge them.  We then copy the new connection lines back in the structure

%    BW.  Nothing is connected to itself up to here.  But,
%    something in here causes same nodes to be connected to
%    themselves.  Let's clear up this logic and see why.
%
tmp = nodeConnections;
h = mrvWaitbar(0,'Removing redundant vertices from connections...');
for k=1:length(uniqueVertices),

  % For each unique node, find rows of nodeConnections
  % that are really the same and refer to that unique node.
  %
  sameNodes = find(uniqueIdx==k)  ;

  % Now, we are going to collapse the larger matrix by copying
  % all of the entries from all of the same nodes into each other.
  % I think this can produce some diagonal entries. -- BW

  % This sums the nodeConnections along the first (row) dimension
  %
  mergedNodes = sum(nodeConnections(sameNodes,:),1);

  % repmat belows up the mergedNodes to the same size or
  % nodeConnections(sameNodes,:)
  % The two are or'd together to make the entries all 0s and 1s
  %
  nodeConnections(sameNodes,:) = or(nodeConnections(sameNodes,:),...
      repmat(mergedNodes,size(nodeConnections(sameNodes,:),1),1)) ;

  % And then this is done again but for the corresponding columns
  %
  nodeConnections(:,sameNodes') = or(nodeConnections(:,sameNodes'),...
      nodeConnections(sameNodes,:)');

  %
  if ~mod(k,100)
    mrvWaitbar(k/length(nodeIdx));
  end
end
close(h);

%% Compacting the connections matrix
%
uniqueSparse = nodeConnections(sort(nodeIdx),sort(nodeIdx)');

% Did the user want a gray graph, too?
%
if nargout == 1,
   return;
elseif nargout == 2
  [grayGraph.nodes, grayGraph.edges] = Connections2Graph(uniqueSparse,uniqueVertices);
  grayGraph.dimdist = mesh.parameters.structural.voxSize;
else
  error('Mesh2Connections:  Wrong number of output arguments.');
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function r = join(rs, joinall)
%JOIN Join similar structures of arrays to one structure of arrays
%
%   JOIN can be used to combine statistics collected from several
%   independent simulations which collect samples to structure
%   of arrays.
%
%   R = JOIN(RS) returns structure of arrays R, RS is array or
%   cell-array of similar structures of arrays. Arrays of R are
%   contanation of arrays in RS. Arrays having first dimension 1
%   are not concataned, but copied from the first element of RS.
%
%   R = JOIN(RS, JOINALL) if JOINALL is 1, arrays having first
%   dimension 1 are concataned.
%
%   See also
%     THIN

% Copyright (C) 2000-2006  Aki Vehtari
%
% This software is distributed under the GNU General Public
% Licence (version 3 or later); please refer to the file
% Licence.txt, included with the software, for details.

if nargin <2
  joinall=false;
end

[m,n]=size(rs);

r=[];
if (m==1 & n==1)
  r=rs;
else
  if iscell(rs)
    rs=[rs{:}];
  end
  if isstruct(rs)
    r=rs(1);
    names = fieldnames(r);
    for i1=1:size(names,1)
      switch names{i1}
       case 'rstate'
        continue
       case 'type'
        r.type=rs(1).(names{1});
        continue
      end
      tmp1=getfield(r,names{i1});
      if iscell(tmp1)
	%tmp2=eval(['cat(1,rs(:).' names{i1} ')']);
	tmp2=cat(1,rs(:).(names{i1}));
	for i2=1:length(tmp1(:))
	  %eval(['r.' names{i1} '{i2}=cat(1,tmp2{:,i2});']);
	  r.(names{i1}){i2}=cat(1,tmp2{:,i2});
	end
      elseif length(tmp1) > 1 || joinall
	%eval(['r.' names{i1} '=cat(1,rs.' names{i1} ');']);
	r.(names{i1})=cat(1,rs.(names{i1}));
      end
    end
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [uniqueEventTypes, boundaryEvents, hardEvents, hardFrames] = ...
                                            getBoundaryEvents(events)
% Returns the event types and boundary events for EEGLAB events structure
% (as extracted from EEG.event).
   uniqueEventTypes = [];
   boundaryEvents = [];
   hardEvents = [];
   hardFrames = [];
   if isempty(events)
       return;
   end
   types = cellfun(@num2str, {events.type}, 'UniformOutput', false);
   uniqueEventTypes = unique(types);
   boundaryMask = strcmpi(types, 'boundary');
   if sum(boundaryMask > 0) && isfield(events, 'duration')
        boundaryEvents = (1:length(events));
        boundaryEvents = boundaryEvents(boundaryMask);
        boundaryDurations = cell2mat({events(boundaryMask).duration}');
        hard = isnan(boundaryDurations);
        if sum(hard) > 0
            hardEvents = boundaryEvents(hard);
            hardFrames = round(cell2mat({events(hardEvents).latency}));
        end
   end
end

author: VisLab

repo: EEG-Clean-Tools


text:
% This script test under 10-fold Cross Validation, divided by ID, how dtfv works
function prob = test_dtfv_CV(part,feaMap, method)
addpath('../../../toolbox/libsvm/matlab');
%part = 1;
%fea_path = '../dataset_trial/FVs';
trainfile = ['../Scripts_by_ID/trainVideo',num2str(part),'.txt'];
testfile = ['../Scripts_by_ID/testVideo',num2str(part),'.txt'];

fid = fopen(trainfile);
C = textscan(fid, '%s');
fclose(fid);
fea_types = {'mbhx','mbhy'};
fea_dim = 24576*2;%+24576+27648;

num_v = length(C{1})/2;
train_fea = zeros(num_v, fea_dim);
train_lab = zeros(num_v, 1);
for i = 1:num_v
    [pathstr,name,ext] = fileparts(C{1}{2*i-1});
    videoname = [name,ext];
    s = findstr(name, 'lie');
    lab = isempty(s);
    fea = [];
    for j = 1: length(fea_types)
        %fea_file = fullfile(fea_path,[videoname,'.',fea_types{j},'.fv.txt']);
        %disp(fea_file)
        %fea = [fea, importdata(fea_file)];
        fea = [fea, feaMap([videoname,'.',fea_types{j},'.fv.txt'])];
    end
    train_fea(i,:) = fea;
    train_lab(i) = lab;
end

%% test phase
fid = fopen(testfile);
C = textscan(fid, '%s');
fclose(fid);
num_v = length(C{1})/2;
test_fea = zeros(num_v, fea_dim);
test_lab = zeros(num_v, 1);
for i = 1:num_v
    [pathstr,name,ext] = fileparts(C{1}{2*i-1});
    videoname = [name,ext];
    s = findstr(name, 'lie');
    lab = isempty(s);
    fea = [];
    for j = 1: length(fea_types)
        %fea_file = fullfile(fea_path,[videoname,'.',fea_types{j},'.fv.txt']);
        %fea = [fea, importdata(fea_file)];
        fea = [fea, feaMap([videoname,'.',fea_types{j},'.fv.txt'])];
    end
    test_fea(i,:) = fea;
    test_lab(i) = lab;
end

switch(method)
case 'NN'
    net = feedforwardnet(10);
    net.trainFcn = 'trainscg';
    net = configure(net, train_fea', train_lab');
    net = train(net, train_fea', train_lab');
    prob = net(test_fea');
case 'tree'
    tc = fitctree(train_fea, train_lab);
    [label,score,node,cnum] = predict(tc, test_fea);
    prob = score(:,1);
case 'randforest'
    BaggedEnsemble = TreeBagger(50,train_fea,train_lab,'OOBPred','On');
    [label,scores] = predict(BaggedEnsemble, test_fea);
    prob = scores(:,1);
case 'bayes'
    flag = bitand(var(train_fea(train_lab==1,:))>1e-10,var(train_fea(train_lab==0,:))>1e-10); %clear 0 variance features
    O1 = fitNaiveBayes(train_fea(:,flag), train_lab);
    C1 = posterior(O1, test_fea(:,flag));
    prob = C1(:,1);
case 'log'
    B = glmfit(train_fea, [train_lab ones(size(train_lab,1),1)], 'binomial', 'link', 'logit')
    Z = repmat(B(1), size(test_lab,1),1) + test_fea*B(2:end);
    prob = 1 ./ (1 + exp(-Z));
    prob = 1-prob;
case 'boost'
    ens = fitensemble(train_fea,train_lab,'AdaBoostM1',100,'Tree')
    [~, prob] = predict(ens,test_fea)
    prob = prob(:,1);
case 'linearsvm'
    %model = svmtrain(train_lab, train_fea, '-t 0 -b 1 -q')
    %fprintf('Finished training.\n');
    %[pred, acc, prob] = svmpredict(test_lab,test_fea, model, '-b 1 -q');
    %prob = prob(:,2);

    model = svmtrain(train_lab, train_fea, '-t 0 -q')
    [pred, acc, prob] = svmpredict(test_lab,test_fea, model, ' -q');
    lie_id = find(prob<0); % find score <0 sample
    if ~isempty(lie_id) % if have such sample
        if pred(lie_id(1)) == 0 % if pred of the sample is deceptive, reverse the score
            isign = -1;
        else
            isign = 1;
        end
    else % if all sample scores >0
        if pred(1) == 1 % if pred is truthful, then reverse all scores
            isign = -1;
        else
            isign = 1;
        end
    end

    prob = isign*prob;
case 'kernelsvm'
    %model = svmtrain(train_lab, tmptrain_fea, '-t 0 -q -b 1');
    %fprintf('Finished training.\n');
    %[pred, acc, prob] = svmpredict(test_lab, tmptest_fea, model, '-q -b 1');
    %prob = prob(:,2);
    model = svmtrain(train_lab, train_fea, '-t 1 -c 1 -g 1 -q');
    [pred, acc, prob] = svmpredict(test_lab, test_fea, model, '-q');

    lie_id = find(prob<0);
    if ~isempty(lie_id)
        if pred(lie_id(1)) == 0
            isign = -1;
        else
            isign = 1;
        end
    else
        if pred(1) == 1
            isign = -1;
        else
            isign = 1;
        end
    end
    prob = isign*prob;
end

author: Doubaibai

repo: DARE

save_path: github-repos/MATLAB/Doubaibai-DARE


text:
function newxy=screenxyabs(monitor,xy)
% newxy = screenxyabs(monitor,xy)  SCREEN X-Y ABSolute coordinates.
%   Calculates absolute screen coordinates from the relative monitor screen
%   coordinates. This function simplifies the use of multiple monitors
%   in Matlab through the use of relative coordinates for each monitor.
%   See also SCREENXYMON.M for calculation of relative coordinates.
%   "newxy" is a vector of xy pairs (pixels) calculated for each input xy pair.
%   "monitor" is a number associated with each number (Default=1).
%       Monitor numbering is defined as the row number of the corresponding
%       monitor in the root property 'MonitorPositions'.
%   "xy" is a vector of xy coordinate pairs (pixels) relative to the monitor's
%       origin (1,1) which is the coordinate of the monitor's lower left corner.
%   Examples:
%       y=screenxyabs(2,[1,1]);  y=screenxyabs(2,[1,1,1024,768]);

%   Copyright 2006 Mirtech, Inc.
%   created 09/25/2006  by Mirko Hrovat on Matlab Ver. 7.2
%   Mirtech, Inc.       email: mhrovat@email.com
%   Uses getmondim.m

switch nargin
    case 2
    case 1
        xy=[];
    case 0
        monitor=[];
        xy=[];
    otherwise
        error ('  Too many arguments!')
end  % switch
if isempty(monitor),    monitor=1;  end
if isempty(xy),         xy=[1,1];   end
xypairs=fix(numel(xy)/2);
xyformat=size(xy);
xy=reshape(xy,[2,xypairs]);     % make each row correspond to a different xy pair
mpos=getmondim(monitor);        % get monitor's origin
if isempty(mpos),
    error ('  Monitor %g does not exist',monitor)
else
    newxy=reshape(xy+repmat([mpos(1)-1;mpos(2)-1],[1,xypairs]),xyformat);
end

author: Sable

repo: mcbench-benchmarks


text:
function f = vlmop2(x)

dim = 2;

transl = 1./sqrt(dim);
part1 = (x(:,1) - transl).^2 + (x(:,2) - transl).^2;
part2 = (x(:,1) + transl).^2 + (x(:,2) + transl).^2;

f(:,1) = 1 - exp( -part1 );
f(:,2) = 1 - exp( -part2 );

author: Eric-Bradford

repo: TS-EMO


text:
function legendre_polynomial_test ( )

%*****************************************************************************80
%
%% LEGENDRE_POLYNOMIAL_TEST tests the LEGENDRE_POLYNOMIAL library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2014
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LEGENDRE_POLYNOMIAL_TEST:\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the LEGENDRE_POLYNOMIAL library.\n' );

  legendre_polynomial_test01 ( );
  legendre_polynomial_test015 ( );
  legendre_polynomial_test016 ( );
  legendre_polynomial_test02 ( );
  legendre_polynomial_test03 ( );
  legendre_polynomial_test04 ( );

  p = 5;
  b = 0.0;
  legendre_polynomial_test05 ( p, b );

  p = 5;
  b = 1.0;
  legendre_polynomial_test05 ( p, b );

  p = 5;
  e = 0;
  legendre_polynomial_test06 ( p, e );

  p = 5;
  e = 1;
  legendre_polynomial_test06 ( p, e );

  legendre_polynomial_test07 ( );
  legendre_polynomial_test08 ( );
  legendre_polynomial_test09 ( );
  legendre_polynomial_test095 ( );

  p = 5;
  legendre_polynomial_test10 ( p );

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

  return
end


text:
% test of shortest path using Dijkstra/A*

name = 'erdos';
name = 'rand_network';
name = 'rand';
name = 'square';
name = 'sanfrancisco';
name = 'northamerica';
name = 'sanjoaquin';
name = 'california';
name = 'oldenburg';
name = 'rand_clusters_triangulation';
name = 'randn_triangulation';
name = 'nefertiti.off';
name = 'rand_triangulation';

options.proba = 0.1;
options.alpha = 0.3;
options.beta = 0.3;

n = Inf;
switch lower(name)
    case 'square'
        n = 50^2;
    case {'rand_triangulation','randn_triangulation', 'rand_clusters_triangulation'}
        n = 5000;
    case 'rand'
        n = 100;
    case 'erdos'
        n = 100;
    case 'rand_network'
        n = 200;

end

disp('Loading file');
[A,vertex] = load_graph(name,n,options);

n = min(n,size(A,1));
A = A(1:n,1:n);
vertex = vertex(:,1:n);
A = max(A,A');

% select two points
clf;
hold on;
gplot(A,vertex');
axis tight;
axis square;
axis off;
start_point = ginput(1)';
plot( start_point(1), start_point(2), 'ro'  );
end_point = ginput(1)';
plot( end_point(1), end_point(2), 'bx'  );
hold off;

% find nearest node
d = compute_distance_to_points(vertex,start_point);
[tmp,start_point] = min(d);
d = compute_distance_to_points(vertex,end_point);
[tmp,end_point] = min(d);

% compute heuristic
v = vertex';
p = v(end_point,:);
H = sqrt( (p(1)-v(:,1)).^2 + (p(2)-v(:,2)).^2 );

% perform front propagation without heuristic
clear options;
options.end_points = end_point;
options.H = [];
disp('Performing classical Dijkstra');
[D,S] = perform_dijkstra(A, start_point, options);

% perform front propagation with euclidean heuristic
clear options;
options.end_points = end_point;
options.H = H;
disp('Performing classical A*');
[D1,S1] = perform_dijkstra(A, start_point, options);

% extract paths
path  = perform_dijkstra_path_extraction(A,D,end_point);
path1 = perform_dijkstra_path_extraction(A,D1,end_point);

% plot the paths
options.point_size = 5;
options.graph_style = 'k';
options.far_point_style = '';
clf;
subplot(1,2,1);
plot_dijkstra(A, vertex, S, path, start_point,end_point, options );
axis square;
title('Classical Dijkstra');
subplot(1,2,2);
plot_dijkstra(A, vertex, S1, path1, start_point,end_point, options );
axis square;
title('Algorihtm A^*');

author: gpeyre

repo: matlab-toolboxes


text:
function [ obj ] = mh_invert( obj )
% Metropolized Gibbs sampling on collapsed HUGE model.
%
% This is a protected method of the tapas_Huge class. It cannot be called
% from outside the class.
%
%

obj = obj.mh_init( );
% number of iterations and burn-in
nIt = obj.options.nvp.numberofiterations;
if isempty(nIt)
    nIt = 2e5;
end
nBi = obj.options.nvp.burnin;
if isempty(nBi) || nBi >= nIt
    nBi = fix(nIt/2);
end

% inverse chain temperature
invTemp = obj.options.nvp.inversetemperature;

% reserve memory
q_nk = zeros(obj.N, obj.K);
obj.trace = struct();
obj.trace.smp = repmat(obj.aux.sample, nIt, 1);
obj.trace.lpr = repmat(obj.aux.lpr, nIt, 1);
psrf = repmat(obj.aux.sample, fix(nIt/obj.const.nPsrf) + 1, 1);

%% ===== MAIN LOOP =====
for iIt = 1:nIt

    % subject level
    for iMhDcm = 1:obj.options.mh.nSteps.dcm
        % sample DCM (parameters)
        obj = mh_sample_dcm( obj, invTemp );
        % sample lambda (hyperparameters)
        obj = mh_sample_noise( obj, invTemp );
    end

    % group level
    if mod(iIt, obj.options.mh.nSteps.knKm) == 0
        k = mod(iIt/obj.options.mh.nSteps.knKm - 1,obj.K) + 1;
        obj = mh_sample_kmhop( obj, k );
    else
        % sample weights pi
        if obj.K > 1
            for iMh = 1:obj.options.mh.nSteps.weights
                obj = mh_sample_weights( obj );
            end
        end
        % sample cluster mu and Sigma
        for iMh = 1:obj.options.mh.nSteps.clusters
            obj = mh_sample_cluster( obj );
        end
    end

    % save sample
    obj.trace.smp(iIt) = obj.aux.sample;
    obj.trace.lpr(iIt) = obj.aux.lpr;

    % adapt proposal step size
    if (iIt <= nBi/3) && (mod(iIt, obj.const.mhAdapt(1)) == 0)
        obj = mh_adapt(obj, iIt);
    end

    % accumulate cluster assigment estimate
    if iIt > nBi
        tmp = bsxfun(@times, exp(obj.aux.rho), obj.aux.sample.pi);
        q_nk = q_nk + bsxfun(@rdivide, tmp, sum(tmp, 2));
    end
    % convergence monitoring
    if mod(iIt, obj.const.nPsrf) == 0
        iPsrf = iIt/obj.const.nPsrf;
        psrf(iPsrf) = mh_psrf(obj.trace.smp(1:iIt), 4);
    end
    if obj.options.nvp.verbose && mod(iIt, 100) == 1
        fprintf('Iteration %u\n', iIt);
    end

% -------------------------------------
end%         END MAIN LOOP
% -------------------------------------
%%          Post Processing
% estimates for assignment probability
q_nk = q_nk/(nIt - nBi);

% acceptance ratio
ratio = struct();
obj.aux.nProp.sp(1) = obj.aux.nProp.sp(1)*obj.N;
for parameter = {'pi', 'mu', 'kappa', 'theta', 'lambda', 'sp'}
    ratio.(parameter{1}) = obj.aux.nAccept.(parameter{1})./...
        obj.aux.nProp.(parameter{1});
end
% posterior mean and quantiles
postMean    = struct();
postVar     = struct();
postQuant   = struct();
for parameter = {'pi', 'mu', 'kappa', 'theta_c', 'theta_h', 'lambda'}
    tmp = reshape([obj.trace.smp(nBi+1:end).(parameter{1})], ...
        [size(obj.aux.sample.(parameter{1})), nIt - nBi]);
    postMean.(parameter{1})  = mean(tmp, 3);
    postVar.(parameter{1})   = var(tmp, 0, 3);
    postQuant.(parameter{1}) = quantile(tmp, obj.options.quantiles, 3);
end
% cumulative probability levels for quantiles
postQuant.levels = obj.options.quantiles;

% calculate PSRF post burn-in
psrf(end) = mh_psrf(obj.trace.smp(nBi + 1:end), 4);

% collect posterior summaries
obj.posterior = struct('nIt', nIt, 'nBi', nBi, 'q_nk', q_nk, ...
    'ratio', ratio, 'psrf', psrf, 'mean', postMean, 'variance', postVar, ...
    'quantile', postQuant, 'lvarBold', obj.aux.lvarBold, ...
    'nrv', exp(-postMean.lambda));

% thin MC chain
if ~isempty(obj.options.nTrace)
    % keep only nTrace samples from post-burn-in phase
    nTrace  = min(nIt - nBi, obj.options.nTrace);
    nThin   = fix((nIt - nBi)/nTrace);
    % select samples uniformly
    obj.trace.smp = obj.trace.smp(end-nThin*(nTrace - 1):nThin:end);
    obj.trace.lpr = obj.trace.lpr(end-nThin*(nTrace - 1):nThin:end);
end

end

%---------------------------------
%            SAMPLING
%---------------------------------
%% SAMPLING: weights (pi)
function [ obj ] = mh_sample_weights( obj )

% propose (in unconstrained space)
prop_piu = obj.aux.sample.pi_u + ...
    randn(1,obj.K-1)*obj.aux.step.pi;
% transform to unit simplex
prop_pis = 1./(1 + exp(log(obj.K-1:-1:1) - prop_piu));
prop_pic = cumprod(1-prop_pis);
prop_pi = [prop_pis(1),-diff(prop_pic),prop_pic(end)];

% evaluate joint (in unconstrained space)
% log-prior on pi
prop_lpr = log(prop_pi)*(obj.prior.alpha_0 - 1) + ...
   sum(log(prop_pis)) + sum(log(1-prop_pis)) + ...
   sum(log(prop_pic(1:end-1)));
prop_lpr = max(prop_lpr, -realmax);
% log-conditional of theta_c given pi
prop_lcd = log(exp(obj.aux.rho)*prop_pi') + obj.aux.rho_max;

% accept/reject
obj.aux.nProp.pi = obj.aux.nProp.pi + 1;
a = exp(prop_lpr - obj.aux.lpr.pi...
        + sum(prop_lcd - obj.aux.lpr.theta_c));
%         a = exp(prop_lpr - obj.aux.lpr.pi);
if ~isnan(a) && ~isinf(a) && rand()<a
    obj.aux.sample.pi_u = prop_piu;
    obj.aux.sample.pi = prop_pi;
    obj.aux.lpr.pi = prop_lpr;
    obj.aux.lpr.theta_c = prop_lcd;
    obj.aux.nAccept.pi = obj.aux.nAccept.pi + 1;
end

end

%% SAMPLING: cluster parameters (mu, kappa = - log(sigma^2))
function [ obj ] = mh_sample_cluster( obj )

obj.aux.nProp.mu = obj.aux.nProp.mu + 1;
obj.aux.nProp.kappa = obj.aux.nProp.mu;

for k = 1:obj.K

    prop_kappa = obj.aux.sample.kappa(k,:);
    % mean
    tmp = randn(1,obj.idx.P_c)*obj.aux.step.mu(k);
    prop_mu = obj.aux.sample.mu(k,:) + tmp*...
        obj.aux.transform.mu(:,:,k);
    dlq = 0; % delta log-proposal density

    % evaluate log-conditional
    % prior
    prop_dmu_k = prop_mu - obj.prior.m_0;
    prop_lpr_mu = -prop_dmu_k*obj.prior.T_0*prop_dmu_k'/2;

    % log-conditional theta_c given mu and kappa
    prop_dtheta_c = bsxfun(@minus, obj.aux.sample.theta_c, prop_mu);
    tmp = obj.aux.l2pi - .5*prop_dtheta_c.^2*exp(prop_kappa') ...
        + .5*sum(prop_kappa);
    if obj.K > 1
        prop_rho = bsxfun(@plus, obj.aux.rho, obj.aux.rho_max);
        prop_rho(:,k) = tmp;
        prop_rho_max = max(prop_rho, [], 2);
        prop_rho = bsxfun(@minus, prop_rho, prop_rho_max);
        prop_lcd = log(exp(prop_rho)*obj.aux.sample.pi') ...
            + prop_rho_max; % log-sum-exp
    else
        prop_rho_max = tmp;
        prop_lcd = prop_rho_max;
        prop_rho = zeros(obj.N, 1);
    end

    % accept/reject
    a = exp(sum(prop_lcd - obj.aux.lpr.theta_c) ...
          + prop_lpr_mu - obj.aux.lpr.mu(k) + dlq);
%             a = exp(prop_lpr_mu - obj.aux.lpr.mu(k));
    if ~isnan(a) && ~isinf(a) && rand()<a
        obj.aux.nAccept.mu(k) = obj.aux.nAccept.mu(k) + 1;
        obj.aux.sample.mu(k,:) = prop_mu;
        obj.aux.rho = prop_rho;
        obj.aux.rho_max = prop_rho_max;
        obj.aux.lpr.theta_c = prop_lcd;
        obj.aux.lpr.mu(k) = prop_lpr_mu;
    end

    % precision
    prop_kappa = prop_kappa + randn(1,obj.idx.P_c)*obj.aux.step.kappa(k);

    % evaluate log-conditional
    % prior
    prop_dkappa = prop_kappa - obj.prior.s_0;
    prop_lpr_kappa = -.5*prop_dkappa.^2*obj.prior.nu_0;
    % log-conditional theta_c given mu and kappa
    prop_dtheta_c = bsxfun(@minus, obj.aux.sample.theta_c, ...
        obj.aux.sample.mu(k,:));
    tmp = obj.aux.l2pi - .5*prop_dtheta_c.^2*exp(prop_kappa') ...
        + .5*sum(prop_kappa);
    if obj.K > 1
        prop_rho = bsxfun(@plus, obj.aux.rho, obj.aux.rho_max);
        prop_rho(:,k) = tmp;
        prop_rho_max = max(prop_rho, [], 2);
        prop_rho = bsxfun(@minus, prop_rho, prop_rho_max);
        prop_lcd = log(exp(prop_rho)*obj.aux.sample.pi') ...
            + prop_rho_max;
    else
        prop_rho_max = tmp;
        prop_lcd = prop_rho_max;
        prop_rho = zeros(obj.N, 1);
    end

    % accept/reject
    a = exp(sum(prop_lcd - obj.aux.lpr.theta_c) ...
          + prop_lpr_kappa - obj.aux.lpr.kappa(k));
%             a = exp(prop_lpr_kappa - obj.aux.lpr.kappa(k));
    if ~isnan(a) && ~isinf(a) && rand()<a
        obj.aux.nAccept.kappa(k) = obj.aux.nAccept.kappa(k) + 1;
        obj.aux.sample.kappa(k,:) = prop_kappa;
        obj.aux.rho = prop_rho;
        obj.aux.rho_max = prop_rho_max;
        obj.aux.lpr.theta_c = prop_lcd;
        obj.aux.lpr.kappa(k) = prop_lpr_kappa;
    end

end

end

%% SAMPLING: DCM parameter (theta)
function [ obj ] = mh_sample_dcm( obj, invTemp )
nProp = obj.aux.nProp.theta + obj.aux.nProp.sp(1) + 1;
bGmm = mod(nProp, obj.options.mh.nSteps.knGmm) == 0;
if bGmm
    obj.aux.nProp.sp(1) = obj.aux.nProp.sp(1) + 1;
else
    obj.aux.nProp.theta = obj.aux.nProp.theta + 1;
end
prec = -.5*exp(obj.aux.sample.kappa);

for n = 1:obj.N

    % propose
    if bGmm % sample from GMM
        k = randsample(obj.K, 1, true, obj.aux.sample.pi);
        prop_theta_c = obj.aux.sample.mu(k,:) + ...
            randn(1, obj.idx.P_c).*exp(-.5*obj.aux.sample.kappa(k,:));
        prop_theta_h = obj.prior.mu_h + randn(1,obj.idx.P_h)./...
            sqrt(diag(obj.prior.Pi_h))'; % assume diagonal precision
    else % sample from Gaussian kernel
%             prop_theta_c = obj.aux.sample.theta_c(n,:) + ...
%                 randn(1, obj.idx.P_c).*obj.aux.step.theta_c(n,:);
        tmp = (randn(1, obj.idx.P_c + obj.idx.P_h).*obj.aux.step.theta(n,:))* ...
            obj.aux.transform.theta;
        prop_theta_c = obj.aux.sample.theta_c(n,:) + tmp(1:obj.idx.P_c);
        prop_theta_h = obj.aux.sample.theta_h(n,:) + tmp(obj.idx.P_c+1:end);
    end

    % evaluate joint
    % log-prior
    prop_dtheta_c = bsxfun(@minus, prop_theta_c, obj.aux.sample.mu);
    prop_rho = obj.aux.l2pi + sum(prec.*prop_dtheta_c.^2, 2)' ...
        + 0.5*sum(obj.aux.sample.kappa, 2)';
    prop_rho_max = max(prop_rho);
    prop_rho = prop_rho - prop_rho_max;
    prop_lpr_c = log(exp(prop_rho)*obj.aux.sample.pi(:)) + prop_rho_max;
    prop_dtheta_h = prop_theta_h - obj.prior.mu_h;
    prop_lpr_h = -.5*prop_dtheta_h*obj.prior.Pi_h*prop_dtheta_h';
    % log-likelihood
    if invTemp
        prop_epsilon = obj.bold_gen( [prop_theta_c, prop_theta_h], ...
            obj.data(n), obj.inputs(n), obj.options.hemo, obj.R, ...
            obj.L, obj.idx );
        tmp = obj.aux.sample.lambda(n,:) - obj.aux.lvarBold(n,:);
        prop_llh = -.5*sum(prop_epsilon.^2*exp(tmp)') ...
           +.5*obj.aux.q_r(n)*sum(tmp);
    else % skip calculation of BOLD signal if inverse temp is zero
        prop_llh = 0;
        prop_epsilon = 0;
    end

    % accept/reject
    a = exp(~bGmm*(prop_lpr_c - obj.aux.lpr.theta_c(n) ...
        + prop_lpr_h - obj.aux.lpr.theta_h(n)) ...
        + invTemp*(prop_llh - obj.aux.lpr.llh(n)));

    if ~isnan(a) && ~isinf(a) && rand()<a
        if bGmm
            obj.aux.nAccept.sp(1) = obj.aux.nAccept.sp(1) + 1;
        else
            obj.aux.nAccept.theta(n) = obj.aux.nAccept.theta(n) + 1;
        end
        obj.aux.sample.theta_c(n,:) = prop_theta_c;
        obj.aux.sample.theta_h(n,:) = prop_theta_h;
        obj.aux.lpr.theta_c(n) = prop_lpr_c;
        obj.aux.lpr.theta_h(n) = prop_lpr_h;
        obj.aux.rho(n,:) = prop_rho;
        obj.aux.rho_max(n) = prop_rho_max;
        obj.aux.epsilon{n} = prop_epsilon;
        obj.aux.lpr.llh(n) = prop_llh;
    end

end
end

%% SAMPLING: Signal-to-noise (lambda)
function [ obj ] = mh_sample_noise( obj, invTemp )

obj.aux.nProp.lambda = obj.aux.nProp.lambda + 1;
for n = 1:obj.N

    % propose in unconstrained space
    prop_lambda = obj.aux.sample.lambda(n,:) + ...
        randn(1, obj.R).*obj.aux.step.lambda(n);

    % evaluate log-conditional
    % prior
    prop_dlambda = prop_lambda - obj.prior.lambda_0;
    prop_lpr = -.5*prop_dlambda.^2*obj.prior.omega_0;
    % likelihood
    tmp = prop_lambda - obj.aux.lvarBold(n,:);
    prop_llh = -.5*sum(obj.aux.epsilon{n}.^2*exp(tmp)') ...
        +.5*obj.aux.q_r(n)*sum(tmp);

    % accept/reject
    a = exp(prop_lpr - obj.aux.lpr.lambda(n) ...
            + invTemp*(prop_llh - obj.aux.lpr.llh(n)));

    if ~isnan(a) && ~isinf(a) && rand()<a
        obj.aux.nAccept.lambda(n) = obj.aux.nAccept.lambda(n) + 1;
        obj.aux.sample.lambda(n,:) = prop_lambda;
        obj.aux.lpr.lambda(n) = prop_lpr;
        obj.aux.lpr.llh(n) = prop_llh;
    end
end
end

%% SAMPLING: k-means-based mode hopping
% A special proposal for clustering part of HUGE model (pi, mu and kappa)
function [ obj ] = mh_sample_kmhop(obj, k)
% track acceptance rate
obj.aux.nProp.sp(2) = obj.aux.nProp.sp(2) + 1;

% --- do kmeans for k = 1,...,K, and define q(...|k) ---
% set up prior
tmpm = exp(obj.prior.s_0+1./2./obj.prior.nu_0');
tmpv = (exp(1./obj.prior.nu_0')-1).*exp(2*obj.prior.s_0+1./obj.prior.nu_0');
b0 = tmpm./tmpv;
a0 = tmpm.^2./tmpv;

% get data
X = obj.aux.sample.theta_c;

% reserve memory
qmm = repmat(obj.prior.m_0,obj.K,1,obj.K);
qms = repmat(diag(obj.prior.T_0)',obj.K,1,obj.K);
qkm = repmat(obj.prior.s_0,obj.K,1,obj.K);
qks = repmat(obj.prior.nu_0',obj.K,1,obj.K);
qpm = zeros(obj.K,obj.K);
for k1 = 1:obj.K
    % do kmeans
    idx = kmeans(X,k1,'Replicates',100);
    % assemble parameters for q
    for k2 = 1:obj.K
        nk = nnz(idx(:)==k2);
        if nk > 0
            mk = mean(X(idx(:)==k2,:),1);
            ak = a0 + (nk + 1)/2;
            bk = b0 + 0.5.*sum(bsxfun(@minus,X(idx(:)==k2,:),mk).^2,1);
            tmpm = ak./bk;
            tmpv = ak./bk.^2;
            % q(kappa)
            qks(k2,:,k1) = 1./log(tmpv./tmpm.^2 + 1);
            qkm(k2,:,k1) = log(tmpm) - 1./2./qks(k2,:,k1);
            % q(mu)
            qms(k2,:,k1) = diag(obj.prior.T_0)' + nk.*exp(obj.prior.s_0);
            qmm(k2,:,k1) = (diag(obj.prior.T_0)'.*obj.prior.m_0 + ...
                nk.*exp(obj.prior.s_0).*mk)./qms(k2,:,k1);
        end
        qpm(k1,k2) = nk + obj.prior.alpha_0(k2);
    end
end

% --- generate proposal for k in 1,...,K, and eval p ---
k1 = k;

prop_lpr_mu = zeros(1,obj.K);
prop_lpr_kappa = zeros(1,obj.K);
prop_rho = obj.aux.rho;

prop_pi = zeros(1,obj.K);
prop_kappa = zeros(obj.K,obj.idx.P_c);
prop_mu = zeros(obj.K,obj.idx.P_c);
for k2 = 1:obj.K
    % draw kappa*
    prop_kappa(k2,:) = qkm(k2,:,k1) + randn(1,obj.idx.P_c)./sqrt(qks(k2,:,k1));
    % draw mu*
    prop_mu(k2,:) = qmm(k2,:,k1) + randn(1,obj.idx.P_c)./sqrt(qms(k2,:,k1));
    % draw pi*
    prop_pi(1,k2) = gamrnd(qpm(k1,k2),1);

    % eval log p(mu*)
    prop_dmu_k = prop_mu(k2,:) - obj.prior.m_0;
    prop_lpr_mu(k2) = -prop_dmu_k*obj.prior.T_0*prop_dmu_k'/2;
    % eval log p(kappa*)
    prop_dkappa = prop_kappa(k2,:) - obj.prior.s_0;
    prop_lpr_kappa(k2) = -.5*prop_dkappa.^2*obj.prior.nu_0;

    % log-conditional theta_c given mu and kappa
    prop_dtheta_c = bsxfun(@minus, obj.aux.sample.theta_c, prop_mu(k2,:));
    prop_rho(:,k2) = obj.aux.l2pi - .5*prop_dtheta_c.^2*exp(prop_kappa(k2,:)') ...
        + .5*sum(prop_kappa(k2,:));

end
prop_pi = prop_pi./sum(prop_pi);

% eval log p(pi*)
prop_pis = [prop_pi(1),prop_pi(2:end-1)./(1 - cumsum(prop_pi(1:end-2)))];
prop_pic = cumprod(1-prop_pis);
prop_piu = tapas_huge_logit(prop_pis) - log(1./(obj.K-1:-1:1));
% log-prior on pi
prop_lpr_pi = log(prop_pi)*(obj.prior.alpha_0 - 1) + ...
   sum(log(prop_pis)) + sum(log(1-prop_pis)) + ...
   sum(log(prop_pic(1:end-1)));
prop_lpr_pi = max(prop_lpr_pi, -realmax);

% eval log p(theta|pi*,mu*,kappa*)
prop_rho_max = max(prop_rho, [], 2);
prop_rho = bsxfun(@minus, prop_rho, prop_rho_max);
prop_lcd = log(exp(prop_rho)*prop_pi') + prop_rho_max; % log-sum-exp

% --- iterate over all permutations and eval q ---
% eval log q(pi*,kappa*,mu*)
[prop_lq_pi,prop_lq_kappa,prop_lq_mu] = eval_lq_perm(prop_pi,...
    prop_kappa,prop_mu,qpm,qkm,qks,qmm,qms,obj.K);

% eval log q(pi,kappa,mu)
[smp_lq_pi,smp_lq_kappa,smp_lq_mu] = eval_lq_perm(obj.aux.sample.pi,...
    obj.aux.sample.kappa,obj.aux.sample.mu,qpm,qkm,qks,qmm,qms,obj.K);

m_prop_lq_mu = logmeanexp(prop_lq_mu);
m_prop_lq_kappa = logmeanexp(prop_lq_kappa);
m_prop_lq_pi = logmeanexp(prop_lq_pi);

m_smp_lq_mu = logmeanexp(smp_lq_mu);
m_smp_lq_kappa = logmeanexp(smp_lq_kappa);
m_smp_lq_pi = logmeanexp(smp_lq_pi);

% --- eval MH acceptance ratio ---
a = [sum(prop_lcd - obj.aux.lpr.theta_c) ...
    ;sum(prop_lpr_kappa - obj.aux.lpr.kappa) ...
    ;sum(prop_lpr_mu - obj.aux.lpr.mu) ...
    ;prop_lpr_pi - obj.aux.lpr.pi ...
    ;sum(m_smp_lq_kappa - m_prop_lq_kappa) ...
    ;sum(m_smp_lq_mu - m_prop_lq_mu) ...
    ;m_smp_lq_pi - m_prop_lq_pi ...
];
% figure(1);clf;stem(a);
a = min(1,exp(sum(a)));

% accept/reject
if ~isnan(a) && ~isinf(a) && rand()<a
    obj.aux.nAccept.sp(2) = obj.aux.nAccept.sp(2) + 1;
    obj.aux.sample.pi_u = prop_piu;
    obj.aux.sample.pi = prop_pi;
    obj.aux.sample.mu = prop_mu;
    obj.aux.sample.kappa = prop_kappa;
    obj.aux.lpr.pi = prop_lpr_pi;
    obj.aux.lpr.mu = prop_lpr_mu;
    obj.aux.lpr.kappa = prop_lpr_kappa;
    obj.aux.lpr.theta_c = prop_lcd;
    obj.aux.rho = prop_rho;
    obj.aux.rho_max = prop_rho_max;
end

end

% evaluate log q(pi*,kappa*,mu*) for all permutations
function [lq_pi,lq_kappa,lq_mu] = eval_lq_perm(p_pi,p_kappa,p_mu,qpm,qkm,qks,qmm,qms,K)
pidx = perms(1:K);

lq_mu = zeros(size(pidx));
lq_kappa = zeros(size(pidx));
lq_pi = zeros(size(pidx));

x = p_pi;
z = [x(1),x(2:end-1)./(1 - cumsum(x(1:end-2)))];
c = cumprod(1-z);

for k1 = 1:K
    m_pi = zeros(K,K);
    m_mu = zeros(size(m_pi));
    m_kappa = zeros(size(m_pi));
    % mesh log q(pi)
    tmp = gammaln(sum(qpm(k1,:))) - sum(gammaln(qpm(k1,:))) + ... % const
    sum(log(z)) + sum(log(1-z)) + ... % jacobian
    sum(log(c(1:end-1))); % jacobian
    m_pi(:,:) = (qpm(k1,:)' - 1)*log(x(:)') + tmp/K; % x

    % mesh log q(mu,kappa)
    for k2 = 1:K
        for k3 = 1:K
            % log p(mu*)
            m_mu(k2,k3) = -.5*(p_mu(k3,:) - qmm(k2,:,k1)).^2*qms(k2,:,k1)' -.5*log(2*pi) +.5*sum(log(qms(k2,:,k1)));
            % log p(kappa*)
            m_kappa(k2,k3) = -.5*(p_kappa(k3,:) - qkm(k2,:,k1)).^2*qks(k2,:,k1)' -.5*log(2*pi) +.5*sum(log(qks(k2,:,k1)));
        end
    end

    linInd = sub2ind([K,K], repmat(1:K,size(pidx,1),1), pidx);
    lq_mu(:,k1) = sum(m_mu(linInd),2);
    lq_kappa(:,k1) = sum(m_kappa(linInd),2);
    lq_pi(:,k1) = sum(m_pi(linInd),2);

end
end

% log(mean(exp(log_probability)))
function lmep = logmeanexp(lpr)
tmp = max(lpr(:));
lmep = log(mean(exp(lpr(:) - tmp))) + tmp;
end

%---------------------------------
%              MISC
%---------------------------------
%% adapt proposal distribution
function [ obj ] = mh_adapt(obj, iIt)

% adapt step size for:
for parameter = {'pi', 'mu', 'kappa', 'theta', 'lambda'}
    % current acceptance rate
    ratio = obj.aux.nAccept.(parameter{1})./obj.aux.nProp.(parameter{1});
    % correction factor
    tmp = ratio - obj.const.mhRate;
    tmp = exp(sign(tmp).*abs(tmp).^3.*obj.const.mhReg);
    % adapt step size
    obj.aux.step.(parameter{1}) = obj.aux.step.(parameter{1}).*tmp;
    % reset counter
    obj.aux.nProp.(parameter{1}) = 0;
    obj.aux.nAccept.(parameter{1}) = ...
        zeros(size(obj.aux.nAccept.(parameter{1})));

end

step = fix(iIt/obj.const.mhTrans);
idx = iIt:-step:1;

% match proposal distribution to posterior covariance for ...
% ... cluster means
mu = reshape([obj.trace.smp(idx).mu], obj.K, obj.idx.P_c, []);
for k = 1:obj.K
    % calculate SVD on empirical covariance of posterior samples
    tmp = cov(permute(mu(k,:,:), [3 2 1]));
    [rotation, scales] = svd(tmp);
    scales = sqrt(diag(scales));
    % limit smallest proposal step size to 5% of maximum step size
    scales(scales < max(scales)*.05) = max(scales)*.05;
    if max(scales) > 0
        obj.aux.transform.mu(:,:,k) = bsxfun(@times, rotation, scales')';
        tmp = sum(log(scales));
        obj.aux.step.mu(k) = obj.aux.step.mu(k)*exp(...
            (obj.aux.logdet.mu(k) - tmp)/obj.idx.P_c);
        obj.aux.logdet.mu(k) = tmp;
    end
end

% ... subject means
theta = [reshape([obj.trace.smp(idx).theta_c], obj.N, obj.idx.P_c, []), ...
         reshape([obj.trace.smp(idx).theta_h], obj.N, obj.idx.P_h, [])];
pooled = zeros(obj.N, size(theta, 3), obj.idx.P_c + obj.idx.P_h);
for n = 1:obj.N
    tmp = permute(theta(n,:,:), [1 3 2]);
    pooled(n, :, :) = bsxfun(@minus, tmp, mean(tmp, 2));
end
pooled = reshape(pooled, [], obj.idx.P_c + obj.idx.P_h);
% calculate SVD on empirical covariance of posterior samples
tmp = cov(pooled);

[rotation, scales] = svd(tmp);
scales = sqrt(diag(scales));
% limit smallest proposal step size to 5% of maximum step size
scales(scales < max(scales)*.05) = max(scales)*.05;
if max(scales) > 0
    obj.aux.transform.theta = bsxfun(@times, rotation, scales')';
    tmp = sum(log(scales));
    tmp = exp((obj.aux.logdet.theta - tmp)/(obj.idx.P_c + obj.idx.P_h));
    for n = 1:obj.N
        obj.aux.step.theta(n) = obj.aux.step.theta(n)*tmp;
    end
    obj.aux.logdet.theta = sum(log(scales));
end

end

%% potential scale reduction factor
function [ psrf ] = mh_psrf( trace, nChains )

nIt = length(trace);
psrf = struct();
for parameter = fieldnames(trace)'
    smpSize = [size(trace(1).(parameter{1})), nIt];
    try
        tmp = permute(reshape([trace(:).(parameter{1})], smpSize), [3 1 2]);

        psrf.(parameter{1}) = tapas_huge_psrf( tmp, nChains );
    catch
        warning('TAPAS:HUGE:convergence', [ 'Potential scale reduction ' ...
            'factor could not be calculated for %s.'], parameter{1});
        psrf.(parameter{1}) = NaN(smpSize);
    end
end

end

author: translationalneuromodeling

repo: tapas


text:
function openface=volface(t)
%
% openface=volface(t)
%
% find the surface patches of a volume
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2009/10/13
%
% input:
%      t: input, volumetric element list, dimension (ne,4)
%
% output:
%      openface: list of faces of the specified volume
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

openface=surfedge(t);

repo: REST


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

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

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

all_preds = all_theta * X';
[max_vals, max_ndxs] = max(all_preds);
p = max_ndxs';

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

end

repo: machine-learning


text:
classdef PCA < handle
    %PCA  Principal Component Analysis class
    %
    % The class is used to compute a special basis for a set of vectors.
    % The basis will consist of eigenvectors of the covariance matrix
    % computed from the input set of vectors. The class PCA can also
    % transform vectors to/from the new coordinate space defined by the
    % basis. Usually, in this new coordinate system, each vector from the
    % original set (and any linear combination of such vectors) can be
    % quite accurately approximated by taking its first few components,
    % corresponding to the eigenvectors of the largest eigenvalues of the
    % covariance matrix. Geometrically it means that you compute a
    % projection of the vector to a subspace formed by a few eigenvectors
    % corresponding to the dominant eigenvalues of the covariance matrix.
    % And usually such a projection is very close to the original vector.
    % So, you can represent the original vector from a high-dimensional
    % space with a much shorter vector consisting of the projected vector's
    % coordinates in the subspace. Such a transformation is also known as
    % Karhunen-Loeve Transform, or KLT. See
    % [PCA](https://en.wikipedia.org/wiki/Principal_component_analysis).
    %
    % ## Example
    % The following shows a quick example of how to reduce dimensionality
    % of samples from 10 to 3.
    %
    %     Xtrain = randn(100,10);
    %     Xtest  = randn(100,10);
    %     pca = cv.PCA(Xtrain, 'MaxComponents',3);
    %     Y = pca.project(Xtest);
    %     Xapprox = pca.backProject(Y);
    %
    % The class also implements the save/load pattern to regular MAT-files,
    % so we can do the following:
    %
    %     pca = cv.PCA(randn(100,5));
    %     save out.mat pca
    %     clear pca
    %
    %     load out.mat
    %     disp(pca)
    %
    % See also: cv.PCA.PCA, cv.SVD
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    % cached properties of their C++ couterpart
    properties (Access = private, Hidden)
        p_eigenvectors
        p_eigenvalues
        p_mean
    end

    % publicly exposed properties (backed by the cached ones)
    properties (Dependent)
        % Eigenvectors of the covariation matrix.
        eigenvectors
        % Eigenvalues of the covariation matrix.
        eigenvalues
        % Mean value subtracted before the projection and added after
        % the back projection.
        mean
    end

    methods
        function this = PCA(varargin)
            %PCA  PCA constructors
            %
            %     pca = cv.PCA()
            %     pca = cv.PCA(data, 'OptionName', optionValue, ...)
            %     pca = cv.PCA(S)
            %
            % ## Input
            % * __data__ input samples stored as matrix rows or matrix columns
            % * __S__ A scalar struct to be imported
            %
            % ## Options
            % Same set of options as the cv.PCA.compute method.
            %
            % The constructor creates an empty PCA object without an argument.
            % Use cv.PCA.compute to analyze data samples. The second form
            % accepts data samples with the same arguments as cv.PCA.compute.
            % The third form create a new PCA object from a scalar struct with
            % 'eigenvectors', 'eigenvalues', and 'mean' fields.
            %
            % See also: cv.PCA, cv.PCA.compute
            %
            this.id = PCA_(0,'new');
            if nargin > 0
                if nargin==1 && isstruct(varargin{1}) && isscalar(varargin{1})
                    S = varargin{1};
                    this.eigenvectors = S.eigenvectors;
                    this.eigenvalues = S.eigenvalues;
                    this.mean = S.mean;
                else
                    this.compute(varargin{:});
                end
            end
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     obj.delete()
            %
            % See also: cv.PCA
            %
            if isempty(this.id), return; end
            PCA_(this.id, 'delete');
        end

        function read(this, fname_or_str, varargin)
            %READ  Read PCA objects from file
            %
            %     obj.read(filename)
            %     obj.read(str, 'FromString',true)
            %     obj.read(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __filename__ Name of the file to read.
            % * __str__ String containing serialized object you want to load.
            %
            % ## Options
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized object.
            %   default false
            %
            % Loads `eigenvalues`, `eigenvectors` and `mean` from specified
            % storage.
            %
            % See also: cv.PCA.write
            %
            PCA_(this.id, 'read', fname_or_str, varargin{:});
        end

        function varargout = write(this, filename)
            %WRITE  Write PCA objects to file
            %
            %     obj.write(filename)
            %     str = obj.write(filename)
            %
            % ## Input
            % * __filename__ Name of the file to write to.
            %
            % ## Output
            % * __str__ optional output. If requested, the object is persisted
            %   to a string in memory instead of writing to disk.
            %
            % Writes `eigenvalues`, `eigenvectors` and `mean` to the specified
            % file.
            %
            % See also: cv.PCA.read
            %
            [varargout{1:nargout}] = PCA_(this.id, 'write', filename);
        end

        function compute(this, data, varargin)
            %COMPUTE  Performs Principal Component Analysis on the supplied dataset
            %
            %     pca.compute(data)
            %     pca.compute(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __data__ input samples stored as the matrix rows or as the
            %   matrix columns
            %
            % ## Options
            % * __DataAs__ Data layout option. Default 'Row'. One of:
            %   * __Row__ indicates that the input samples are stored as
            %     matrix rows.
            %   * __Col__ indicates that the input samples are stored as
            %     matrix columns.
            % * __MaxComponents__ Maximum number of components that PCA should
            %   retain. default 0 (all the components are retained).
            % * __RetainedVariance__ Percentage of variance that PCA should
            %   retain. Using this parameter will let the PCA decided how many
            %   components to retain but it will always keep at least 2.
            %   default 1.0 (all the components are retained).
            % * __Mean__ Optional mean value. By default, the mean is computed
            %   from the data. Not set by default.
            %
            % **Note**: `RetainedVariance` and `MaxComponents` are mutually
            % exclusive options, and shoudn't be used together.
            %
            % The method performs PCA of the supplied dataset. It is safe
            % to reuse the same PCA structure for multiple datasets. That
            % is, if the structure has been previously used with another
            % dataset, the existing internal data is reclaimed and the new
            % `eigenvalues`, `eigenvectors` and `mean` are allocated and
            % computed.
            %
            % The computed `eigenvalues` are sorted from the largest to the
            % smallest and the corresponding `eigenvectors` are stored as
            % cv.PCA.eigenvectors rows.
            %
            % See also: cv.PCA.PCA, cv.PCA.project, cv.PCA.backProject
            %
            PCA_(this.id, 'compute', data, varargin{:});

            % invalidate cached properties
            this.p_eigenvectors = [];
            this.p_eigenvalues = [];
            this.p_mean = [];
        end

        function Y = project(this, X)
            %PROJECT  Projects vector(s) to the principal component subspace
            %
            %     Y = pca.project(X)
            %
            % ## Input
            % * __X__ input vector(s); must have the same dimensionality and
            %   the same layout as the input data used at PCA phase, that is,
            %   if 'Row' was specified, then `size(X,2)==size(data,2)` (vector
            %   dimensionality) and `size(X,1)` is the number of vectors to
            %   project, and the same is true for the 'Col' case.
            %
            % ## Output
            % * __Y__ output vectors (PC coefficients); in case of 'Col', the
            %   output matrix has as many columns as the number of input
            %   vectors, this means that `size(Y,2)==size(X,2)` and the number
            %   of rows match the number of principal components (for example,
            %   `MaxComponents` parameter passed to the constructor).
            %
            % The method project one or more vectors to the principal
            % component subspace, where each vector projection is
            % represented by coefficients in the principal component basis.
            %
            % See also: cv.PCA, cv.PCA.backProject
            %
            Y = PCA_(this.id, 'project', X);
        end

        function X = backProject(this, Y)
            %BACKPROJECT  Reconstructs vectors from their PC projections
            %
            %     X = pca.backProject(Y)
            %
            % ## Input
            % * __Y__ coordinates of the vectors in the principal component
            %   subspace, the layout and size are the same as of
            %   cv.PCA.project output vectors.
            %
            % ## Output
            % * __X__ reconstructed vectors; the layout and size are the same
            %   as of cv.PCA.project input vectors.
            %
            % The method is the inverse operation to cv.PCA.project. It
            % takes PC coordinates of projected vectors and reconstruct the
            % original vectors. Unless all the principal components have
            % been retained, the reconstructed vectors are different from
            % the originals. But typically, the difference is small if the
            % number of components is large enough (but still much smaller
            % than the original vector dimensionality). As a result, PCA is
            % used.
            %
            % See also: cv.PCA.compute, cv.PCA.project
            %
            X = PCA_(this.id, 'backProject', Y);
        end
    end

    %% Getters/Setters for dependent properties
    methods
        function value = get.eigenvectors(this)
            if isempty(this.p_eigenvectors)
                this.p_eigenvectors = PCA_(this.id, 'get', 'eigenvectors');
            end
            value = this.p_eigenvectors;
        end
        function set.eigenvectors(this, value)
            PCA_(this.id, 'set', 'eigenvectors', value);
            this.p_eigenvectors = value;
        end

        function value = get.eigenvalues(this)
            if isempty(this.p_eigenvalues)
                this.p_eigenvalues = PCA_(this.id, 'get', 'eigenvalues');
            end
            value = this.p_eigenvalues;
        end
        function set.eigenvalues(this, value)
            PCA_(this.id, 'set', 'eigenvalues', value);
            this.p_eigenvalues = value;
        end

        function value = get.mean(this)
            if isempty(this.p_mean)
                this.p_mean = PCA_(this.id, 'get', 'mean');
            end
            value = this.p_mean;
        end
        function set.mean(this, value)
            PCA_(this.id, 'set', 'mean', value);
            this.p_mean = value;
        end
    end

    methods (Hidden)
        function S = struct(this)
            %STRUCT  Converts to a struct array
            %
            %     S = struct(obj)
            %
            % ## Output
            % * __S__ output struct array
            %
            % See also: cv.PCA
            %
            S = struct('eigenvectors',{this.eigenvectors},...
                       'eigenvalues', {this.eigenvalues},...
                       'mean',        {this.mean});
        end

        function S = saveobj(this)
            %SAVEOBJ  Serialization before save
            %
            %     S = obj.saveobj()
            %
            % ## Output
            % * __S__ output struct.
            %
            % Called by the `save` function, when saving object to a MAT-file.
            %
            % See also: cv.PCA.loadobj, saveobj, loadobj
            %
            S = struct(this);
        end
    end

    methods (Static, Hidden)
        function this = loadobj(S)
            %LOADOBJ  Deserialization after load
            %
            %     obj = loadobj(S)
            %
            % ## Input
            % * __S__ input struct.
            %
            % Called by the `load` function, when loading object from a
            % MAT-file.
            %
            % See also: cv.PCA.PCA, cv.PCA.saveobj, saveobj, loadobj
            %
            this = cv.PCA(S);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function varargout=lpcstable(varargin)
%
% For calling details please see v_lpcstable.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcstable.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcstable.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcstable(varargin{:});
else
    v_lpcstable(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
% Here the training data is adapted from UCI ML repository, 'zoo' data

dtreeCPD=tree_CPD;

% load data
fname = fullfile(BNT_HOME, 'examples', 'static', 'uci_data', 'zoo', 'zoo1.data')
data=load(fname);
data=data';

data=transform_data_into_bnt_format(data, []);

% learn decision tree from data
ns=2*ones(1,17);
ns(13)=6;
ns(17)=7;
dtreeCPD1=learn_params(dtreeCPD,1:17,data,ns,[],'stop_cases',5); % a node with less than 5 cases will not be splitted

% evaluate on data
[score,outputs]=evaluate_tree_performance(dtreeCPD1,1:17,data,ns,[]);
fprintf('Accuracy in old training data %6.3f\n',score);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function shmDataSimulated(x0, springConstant, noiseFactor, outFileName)
% Generate simulated pendulum data

%   Copyright 2008-2009 The MathWorks, Inc.
%   $Revision: 35 $    $Date: 2009-05-29 15:27:34 +0100 (Fri, 29 May 2009) $

if nargin<4
    outFileName = 'shmData.xls';
end

mass = 1;  % kg

% Determine the natural frequency of the pendulum
omega0 = sqrt(springConstant/mass);

% Simulate the pendulum motion
res = shmSimulation(x0,springConstant,mass);

% Add noise to position and velocity
positionNoise = randn(size(res.Position))*noiseFactor*x0;
velocityNoise = randn(size(res.Position))*noiseFactor*x0*omega0;

% Include this noise into the position and velocity
res.Position = res.Position + positionNoise;
res.Velocity = res.Velocity + velocityNoise;

% Format the data for writing to the Excel spreadsheet.
colHeadings = {'Time', 'Position', 'Velocity'};
dataArray = [res.Time res.Position res.Velocity];

% Create the cell array to write to the Excel file by combining the column
% headings and the numeric simulation results.  To do this we need to
% convert the array of numeric data to a cell array using num2cell.  Use
% cell2mat to convert a cell array of numbers back to a matrix.
dataArray = [colHeadings; num2cell(dataArray)];
% Save the data
xlswrite(outFileName,dataArray);


text:
function [verts, faces, normals] = read_off(filename)

% read_off - read data from OFF file.
%
%   [verts, faces, normals] = read_off(filename);
%
%   'verts' is a 'nb.vert x 3' array specifying the position of the vertices.
%   'faces' is a 'nb.face x 3' array specifying the connectivity of the mesh.
%   'normals' is a 'nb.normal * 3' array specifying the normal of each vertex.
%   Copyright (c) 2003 Gabriel Peyr?
%   Changed by (2009) JJCAO

fid = fopen(filename,'r');
if( fid==-1 )
    warning('Can''t open the file.');
    verts = [];faces=[];normals=[];
    return;
end

fgetl(fid);
str = fgetl(fid);
tmp = sscanf(str,'%d %d');
nverts = tmp(1);
nfaces = tmp(2);
str = fgetl(fid);
tmp = sscanf(str,'%f %f %f');

frewind(fid);
fgetl(fid); fgetl(fid);
if length(tmp) < 4 % only x y z
    [A,cnt] = fscanf(fid,'%f %f %f', 3*nverts);
    if cnt~=3*nverts
        warning('Problem in reading vertices.');
    end
    A = reshape(A, 3, cnt/3);
    verts = A';
    normals = [];
else % x y z nx ny nz
    [A,cnt] = fscanf(fid,'%f %f %f %f %f %f', 6*nverts);
    if cnt~=6*nverts
        warning('Problem in reading vertices.');
    end

    A = reshape(A, 6, cnt/6);
    A = A';
    verts = A(:,1:3);
    normals = A(:,4:6);
end

if nfaces > 0 % read Face 1  1088 480 1022
    [A,cnt] = fscanf(fid,'%d %d %d %d\n', 4*nfaces);
    if cnt~=4*nfaces
        warning('Problem in reading faces.');
    end
    A = reshape(A, 4, cnt/4);
    faces = A(2:4,:)+1;
    faces = faces';
else
    faces = [];
end

fclose(fid);
return;

author: taiya

repo: cloudcontr

save_path: github-repos/MATLAB/taiya-cloudcontr


text:
function [idx, theta, rho] = PolarCoordinatesGrid(y, r, nDivs)
%   PolarCoordinatesGrid - Calculate index of grid cell and
%   polar coordinates for a given objective vector
%
%   [idx, theta, rho] = PolarCoordinatesGrid(y, r, nDivs)
%
%   Input:
%   y - individual's objectives
%   r - reference point
%   nDivs - number of divisions in each right angle
%
%   Output:
%   idx - index of grid cell where individual resides
%   theta - polar angles
%   rho - radius

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Roman Denysiuk

%% calculate polar coordinates

% determine dimensionality
m = numel(r);

% calculate radius
rho = norm(y-r);

% calculate angles
theta = zeros(m-1,1);

for i = 0:m-2

    u = (y(m-i) - r(m-i))/rho;

    for j = 0:i-1
        u = u/cos(theta(j+1));
    end

    u = min( max(u, -1), 1); % ensure feasibility

    theta(i+1) = asin(u);

    theta(i+1) = min( max(theta(i+1), eps), pi/2-eps);
end

%% calculate grid index based on polar coordinates

% calculate grid coordinates (integer values)
G = ceil(2*nDivs*theta/pi);

% calculate index of pop member using grid coordinates
% formula: idx=G(1)*ndiv^0+(G(2)-1)*ndiv^1+(G(3)-1)*ndiv^2 ...

idx = G(1);
for i = 2:numel(G)
    idx = idx + (G(i)-1)*power(nDivs, i-1);
end

end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [bus, gen, branch, f, success, info, et, g, jac, xr, pimul] = ...
    uopf(varargin)
%UOPF  Solves combined unit decommitment / optimal power flow.
%   [RESULTS, SUCCESS] = UOPF(MPC, MPOPT)
%
%   Returns either a RESULTS struct and an optional SUCCESS flag, or individual
%   data matrices, the objective function value and a SUCCESS flag. In the
%   latter case, there are additional optional return values. See Examples
%   below for the possible calling syntax options.
%
%   Examples:
%       Output argument options:
%
%       results = uopf(...)
%       [results, success] = uopf(...)
%       [bus, gen, branch, f, success] = uopf(...)
%       [bus, gen, branch, f, success, info, et, g, jac, xr, pimul] = uopf(...)
%
%       Input arguments options:
%
%       uopf(mpc)
%       uopf(mpc, mpopt)
%       uopf(mpc, userfcn, mpopt)
%       uopf(mpc, A, l, u)
%       uopf(mpc, A, l, u, mpopt)
%       uopf(mpc, A, l, u, mpopt, N, fparm, H, Cw)
%       uopf(mpc, A, l, u, mpopt, N, fparm, H, Cw, z0, zl, zu)
%
%       uopf(baseMVA, bus, gen, branch, areas, gencost)
%       uopf(baseMVA, bus, gen, branch, areas, gencost, mpopt)
%       uopf(baseMVA, bus, gen, branch, areas, gencost, userfcn, mpopt)
%       uopf(baseMVA, bus, gen, branch, areas, gencost, A, l, u)
%       uopf(baseMVA, bus, gen, branch, areas, gencost, A, l, u, mpopt)
%       uopf(baseMVA, bus, gen, branch, areas, gencost, A, l, u, ...
%                                   mpopt, N, fparm, H, Cw)
%       uopf(baseMVA, bus, gen, branch, areas, gencost, A, l, u, ...
%                                   mpopt, N, fparm, H, Cw, z0, zl, zu)
%
%   See OPF for more information on input and output arguments.
%
%   Solves a combined unit decommitment and optimal power flow for a single
%   time period. Uses an algorithm similar to dynamic programming. It proceeds
%   through a sequence of stages, where stage N has N generators shut down,
%   starting with N=0. In each stage, it forms a list of candidates (gens at
%   their Pmin limits) and computes the cost with each one of them shut down.
%   It selects the least cost case as the starting point for the next stage,
%   continuing until there are no more candidates to be shut down or no
%   more improvement can be gained by shutting something down.
%   If MPOPT.verbose (see MPOPTION) is true, it prints progress
%   info, if it is > 1 it prints the output of each individual opf.
%
%   See also OPF, RUNUOPF.

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

%%----- initialization -----
t0 = tic;       %% start timer

%% process input arguments
[mpc, mpopt] = opf_args(varargin{:});

%% options
if mpopt.verbose    %% turn down verbosity one level for calls to opf
    mpopt = mpoption(mpopt, 'verbose', mpopt.verbose-1);
end

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[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;

%%-----  do combined unit commitment/optimal power flow  -----

%% check for sum(Pmin) > total load, decommit as necessary
on   = find( mpc.gen(:, GEN_STATUS) > 0 & ~isload(mpc.gen) );   %% gens in service
onld = find( mpc.gen(:, GEN_STATUS) > 0 &  isload(mpc.gen) );   %% disp loads in serv
load_capacity = sum(mpc.bus(:, PD)) - sum(mpc.gen(onld, PMIN)); %% total load capacity
Pmin = mpc.gen(on, PMIN);
while sum(Pmin) > load_capacity
    %% shut down most expensive unit
    avgPmincost = totcost(mpc.gencost(on, :), Pmin) ./ Pmin;
    [junk, i] = fairmax(avgPmincost);   %% pick one with max avg cost at Pmin
    i = on(i);                          %% convert to generator index

    if mpopt.verbose
        fprintf('Shutting down generator %d so all Pmin limits can be satisfied.\n', i);
    end

    %% set generation to zero
    mpc.gen(i, [ PG QG GEN_STATUS ]) = 0;

    %% update minimum gen capacity
    on  = find( mpc.gen(:, GEN_STATUS) > 0 & ~isload(mpc.gen) );   %% gens in service
    Pmin = mpc.gen(on, PMIN);
end
if ~any(mpc.gen(:, GEN_STATUS) > 0)     %% don't bother to run anything if
    success = 0;                        %% everything has been shut down
    results0 = mpc;
    results0.success = success;
    results0.f = NaN;
    results0.et = 0;
    if mpopt.verbose
        fprintf('Infeasible problem, Pmin limits cannot be satisfied without shutting down all generators.\n');
    end
else
    %% run initial opf
    [results, success] = opf(mpc, mpopt);

    %% best case so far
    results1 = results;

    %% best case for this stage (ie. with n gens shut down, n=0,1,2 ...)
    results0 = results1;
    mpc.bus = results0.bus;     %% use these V as starting point for OPF

    while 1
        %% get candidates for shutdown
        candidates = find(results0.gen(:, MU_PMIN) > 0 & results0.gen(:, PMIN) > 0);
        if isempty(candidates)
            break;
        end
        done = 1;   %% do not check for further decommitment unless we
                    %%  see something better during this stage
        for i = 1:length(candidates)
            k = candidates(i);
            %% start with best for this stage
            mpc.gen = results0.gen;

            %% shut down gen k
            mpc.gen(k, [ PG QG GEN_STATUS ]) = 0;

            %% run opf
            if any(mpc.gen(:, GEN_STATUS) > 0)
                [results, success] = opf(mpc, mpopt);
            else
                success = 0;
            end

            %% something better?
            if success && results.f < results1.f
                results1 = results;
                k1 = k;
                done = 0;   %% make sure we check for further decommitment
            end
        end

        if done
            %% decommits at this stage did not help, so let's quit
            break;
        else
            %% shutting something else down helps, so let's keep going
            if mpopt.verbose
                fprintf('Shutting down generator %d.\n', k1);
            end

            results0 = results1;
            mpc.bus = results0.bus;     %% use these V as starting point for OPF
        end
    end
end

%% compute elapsed time
et = toc(t0);

%% finish preparing output
if nargout > 0
  success = results0.success;
  if nargout <= 2
    results0.et = et;
    bus = results0;
    gen = success;
  else
    [bus, gen, branch, f, info, xr, pimul] = deal(results0.bus, results0.gen, ...
                    results0.branch, results0.f, results0.raw.info, ...
                    results0.raw.xr, results0.raw.pimul);
    if isfield(results0, 'g')
      g = results0.g;
    end
    if isfield(results0, 'dg')
      jac = results0.dg;
    end
  end
elseif results0.success
  results0.et = et;
  printpf(results0, 1, mpopt);
end

function [val, idx] = fairmax(x)
%FAIRMAX    Same as built-in MAX, except breaks ties randomly.
%   [VAL, IDX] = FAIRMAX(X) takes a vector as an argument and returns
%   the same output as the built-in function MAX with two output
%   parameters, except that where the maximum value occurs at more
%   than one position in the  vector, the index is chosen randomly
%   from these positions as opposed to just choosing the first occurance.
%
%   See also MAX.

val = max(x);               %% find max value
i   = find(x == val);       %% find all positions where this occurs
n   = length(i);            %% number of occurences
idx = i( fix(n*rand)+1 );   %% select index randomly among occurances

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function ccii_1_tabulate ( )

%*****************************************************************************80
%
%% CCII_1_TABULATE tabulates CCII_1 quadrature rules for the Hermite integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CCII_1_TABULATE\n' );
  fprintf ( 1, '  Tabulate CCII_1 quadrature rules for the Hermite integral.\n' );
  fprintf ( 1, '  Density function rho(x) = 1.\n' );
  fprintf ( 1, '  Region: -oo < x < +oo.\n' );
  fprintf ( 1, '  Exactness: NONE.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Print the first 5 rules:\n' );
  fprintf ( 1, '   I          X(I)            W(I)\n' );
  fprintf ( 1, '\n' );

  ell = 1.0;

  for n = 1 : 5

    [ x, w ] = ccii_1 ( n, ell );

    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %2d  %14.6f  %14.6f\n', i, x(i), w(i));
    end
    fprintf ( 1, ' Sum                  %14.6f  %14.6f\n', sum ( w(1:n) ) );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Print the first 5 nested rules:\n' );
  fprintf ( 1, '   I          X(I)            W(I)\n' );
  fprintf ( 1, '\n' );

  ell = 1.0;

  for j = 1 : 5

    n = 2^j - 1;

    [ x, w ] = ccii_1 ( n, ell );

    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %2d  %14.6f  %14.6f\n', i, x(i), w(i) );
    end
    fprintf ( 1, ' Sum                  %14.6f\n', sum ( w(1:n) ) );

  end

  return
end

repo: jburkardt-m


text:
function hermite_test04 ( )

%*****************************************************************************80
%
%% HERMITE_TEST04 interpolates the Runge function using equally spaced data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_TEST04\n' );
  fprintf ( 1, '  HERMITE computes the Hermite interpolant to data.\n' );
  fprintf ( 1, '  Here, f(x) is the Runge function\n' );
  fprintf ( 1, '  and the data is evaluated at equally spaced points.\n' );
  fprintf ( 1, '  As N increases, the maximum error grows.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N     Max | F(X) - H(F(X)) |\n' );
  fprintf ( 1, '\n' );

  for n = 3 : 2 : 15

    nd = 2 * ( n - 1 ) + 1;
    ns = 10 * ( n - 1 ) + 1;
    xlo = -5.0;
    xhi = +5.0;
    x = linspace ( xlo, xhi, n );

    y(1:n) = 1.0 ./ ( 1.0 + x(1:n).^2 );
    yp(1:n) = - 2.0 * x(1:n) ./ ( 1.0 + x(1:n).^2 ).^2;

    [ xd, yd, xdp, ydp ] = hermite_interpolant ( n, x, y, yp );
%
%  Compare exact and interpolant at sample points.
%
    xs = linspace ( xlo, xhi, ns );

    ys = dif_vals ( nd, xd, yd, ns, xs );

    max_dif = 0.0;
    for i = 1 : ns
      xt = xs(i);
      yt = 1.0 / ( 1.0 + xt * xt );
      max_dif = max ( max_dif, abs ( ys(i) - yt ) );
    end

    fprintf ( 1, '  %4d  %14g\n', n, max_dif );

  end

  return
end


text:
function S = randskew(n, N)
% Generates random skew symmetric matrices with normal entries.
%
% function S = randskew(n)
% function S = randskew(n, N)
%
% S is an n-by-n-by-N array where each slice S(:, :, i) for i = 1..N is a
% random skew-symmetric matrix with upper triangular entries distributed
% independently following a normal distribution (Gaussian, zero mean, unit
% variance).
%
% By default, N = 1.
%
% See also: randrot randsym randskewh

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Sept. 25, 2012.
% Contributors:
% Change log:
%       June 19, 2019 (NB):
%           Now handles the case n = 1 properly.

    if nargin < 2
        N = 1;
    end

    if n == 1
        S = zeros(1, 1, N);
        return;
    end

    % Subindices of the (strictly) upper triangular entries of an n-by-n
    % matrix
    [I, J] = find(triu(ones(n), 1));

    K = repmat(1:N, n*(n-1)/2, 1);

    % Indices of the strictly upper triangular entries of all N slices of
    % an n-by-n-by-N matrix
    L = sub2ind([n n N], repmat(I, N, 1), repmat(J, N, 1), K(:));

    % Allocate memory for N random skew matrices of size n-by-n and
    % populate each upper triangular entry with a random number following a
    % normal distribution and copy them with opposite sign on the
    % corresponding lower triangular side.
    S = zeros(n, n, N);
    S(L) = randn(size(L));
    S = S - multitransp(S);

end

author: NicolasBoumal

repo: manopt


text:
function varargout=writehtk(varargin)
%
% For calling details please see v_writehtk.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_writehtk.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: writehtk.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_writehtk(varargin{:});
else
    v_writehtk(varargin{:});
end


text:
classdef NodalFieldPlotter < handle

    properties (Access = private)
        mesh
        field
    end

    methods (Access = public)

        function obj = NodalFieldPlotter(cParams)
            obj.init(cParams)
        end

        function plot(obj)
            x = obj.mesh.coord(:,1);
            y = obj.mesh.coord(:,2);
            z = obj.field;
            %figure()
            trisurf(obj.mesh.connec,x,y,z)
            view(0,90)
            colorbar
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh = cParams.mesh;
            obj.field = cParams.field;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function FVToolStartUp()
%
% SYNOPSIS:
%   FVToolStartUp()
%
% PARAMETERS:
%   No perameter
%
% RETURNS:
%   None
%
% EXAMPLE:
%   n.a.
%
% SEE ALSO:
%   PVTinitialize, FVTdemo

%{
Copyright (c) 2012-2021 Ali Akbar Eftekhari
All rights reserved.

Redistribution and use in source and binary forms, with or
without modification, are permitted provided that the following
conditions are met:

    *   Redistributions of source code must retain the above copyright notice,
        this list of conditions and the following disclaimer.
    *   Redistributions in binary form must reproduce the above
        copyright notice, this list of conditions and the following
        disclaimer in the documentation and/or other materials provided
        with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%}
try
p = mfilename('fullpath');
file_name = mfilename;
current_path = p(1:end-1-length(file_name));
addpath([current_path '/Boundary']);
addpath([current_path '/Calculus']);
addpath([current_path '/Classes']);
addpath([current_path '/Discretization']);
addpath([current_path '/MeshGeneration']);
addpath([current_path '/Solvers']);
addpath([current_path '/Utilities']);
addpath([current_path '/Visualization']);
addpath([current_path '/Examples']);
addpath([current_path '/Physics']);
addpath([current_path '/Tests']);

try
    addpath([current_path '/PhysicalProperties']);
    addpath([current_path '/FieldGeology']);
catch
    disp(['Some of the physical functions are not available in this copy.' ...
    ' It does not affect the functionality of the FVMtool']);
end

% Check for other solvers
    % check for AGMG availability
    cd(current_path);
    cd('Solvers');
    if exist('AGMG_3.2', 'dir') == 7
        addpath([pwd '/AGMG_3.2']);
        disp('AGMG 3.2 linear solver is available.');
    elseif exist('AGMG_3.0', 'dir') == 7
        addpath([pwd '/AGMG_3.0']);
        disp('AGMG 3.0 linear solver is available.');
    else
        disp('AGMG 3.x linear solver is NOT available (Not necessary).');
    end
    % check for Factorize availability
    if exist('Factorize', 'dir') == 7
        addpath([pwd '/Factorize']);
        disp('Factorize is available.');
    end
    cd(current_path);
% end of check for other solvers

% Check for the PVT package
if exist('PVTtoolbox', 'dir') == 7
    addpath([current_path '/PVTtoolbox']);
    if exist('PVTinitialize', 'file') == 2
        cd('PVTtoolbox');
        PVTinitialize();
        cd(current_path);
        disp('PVTtoolbox has started successfully.');
    else
        disp('PVTtoolbox is found but cannot be initialized.');
    end
else
    disp('PVTtoolbox is NOT available (Not necessary).');
end
disp('FiniteVolumeToolbox has started successfully.');
catch err
    error('An error occured while tryng to start the FiniteVolumeToolbox.')
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function obj = D_func_hand_eye_new(in1,in2,in3,in4)
coef_f0_q_sym1 = in1(:,1);
coef_f0_q_sym2 = in1(:,2);
coef_f0_q_sym3 = in1(:,3);
coef_f0_q_sym4 = in1(:,4);
coef_f0_q_sym5 = in1(:,5);
coef_f0_q_sym6 = in1(:,6);
coef_f0_q_sym7 = in1(:,7);
coef_f0_q_sym8 = in1(:,8);
coef_f0_q_sym9 = in1(:,9);
coef_f0_q_sym10 = in1(:,10);
coef_f1_q_sym1 = in2(:,1);
coef_f0_q_sym12 = in1(:,12);
coef_f1_q_sym2 = in2(:,2);
coef_f0_q_sym13 = in1(:,13);
coef_f1_q_sym3 = in2(:,3);
coef_f0_q_sym14 = in1(:,14);
coef_f1_q_sym4 = in2(:,4);
coef_f0_q_sym15 = in1(:,15);
coef_f1_q_sym5 = in2(:,5);
coef_f0_q_sym16 = in1(:,16);
coef_f1_q_sym6 = in2(:,6);
coef_f0_q_sym17 = in1(:,17);
coef_f1_q_sym7 = in2(:,7);
coef_f1_q_sym8 = in2(:,8);
coef_f0_q_sym19 = in1(:,19);
coef_f1_q_sym9 = in2(:,9);
coef_f0_q_sym20 = in1(:,20);
coef_f0_q_sym21 = in1(:,21);
coef_f2_q_sym1 = in3(:,1);
coef_f2_q_sym2 = in3(:,2);
coef_f0_q_sym23 = in1(:,23);
coef_f2_q_sym3 = in3(:,3);
coef_f2_q_sym4 = in3(:,4);
coef_f2_q_sym5 = in3(:,5);
coef_f2_q_sym6 = in3(:,6);
coef_f2_q_sym7 = in3(:,7);
coef_f2_q_sym8 = in3(:,8);
coef_f2_q_sym9 = in3(:,9);
coef_f3_q_sym1 = in4(:,1);
coef_f3_q_sym2 = in4(:,2);
coef_f3_q_sym3 = in4(:,3);
coef_f3_q_sym4 = in4(:,4);
coef_f3_q_sym5 = in4(:,5);
coef_f3_q_sym6 = in4(:,6);
coef_f3_q_sym7 = in4(:,7);
coef_f3_q_sym8 = in4(:,8);
coef_f3_q_sym9 = in4(:,9);
coef_f1_q_sym10 = in2(:,10);
coef_f1_q_sym12 = in2(:,12);
coef_f1_q_sym13 = in2(:,13);
coef_f1_q_sym14 = in2(:,14);
coef_f1_q_sym15 = in2(:,15);
coef_f1_q_sym16 = in2(:,16);
coef_f1_q_sym17 = in2(:,17);
coef_f1_q_sym19 = in2(:,19);
coef_f1_q_sym20 = in2(:,20);
coef_f1_q_sym21 = in2(:,21);
coef_f1_q_sym23 = in2(:,23);
coef_f2_q_sym10 = in3(:,10);
coef_f2_q_sym12 = in3(:,12);
coef_f2_q_sym13 = in3(:,13);
coef_f2_q_sym14 = in3(:,14);
coef_f2_q_sym15 = in3(:,15);
coef_f2_q_sym16 = in3(:,16);
coef_f2_q_sym17 = in3(:,17);
coef_f2_q_sym19 = in3(:,19);
coef_f2_q_sym20 = in3(:,20);
coef_f2_q_sym21 = in3(:,21);
coef_f2_q_sym23 = in3(:,23);
coef_f3_q_sym10 = in4(:,10);
coef_f3_q_sym12 = in4(:,12);
coef_f3_q_sym13 = in4(:,13);
coef_f3_q_sym14 = in4(:,14);
coef_f3_q_sym15 = in4(:,15);
coef_f3_q_sym16 = in4(:,16);
coef_f3_q_sym17 = in4(:,17);
coef_f3_q_sym19 = in4(:,19);
coef_f3_q_sym20 = in4(:,20);
coef_f3_q_sym21 = in4(:,21);
coef_f3_q_sym23 = in4(:,23);
t2 = coef_f1_q_sym1.*2.0;
t3 = coef_f2_q_sym1.*2.0;
t4 = coef_f3_q_sym1.*2.0;
t5 = -coef_f0_q_sym2;
t6 = -coef_f0_q_sym3;
t7 = -coef_f0_q_sym4;
t8 = -coef_f1_q_sym1;
t10 = -coef_f2_q_sym1;
t12 = -coef_f3_q_sym1;
t9 = -t2;
t11 = -t3;
t13 = -t4;
obj = reshape([coef_f0_q_sym1-coef_f1_q_sym2,-coef_f2_q_sym2,-coef_f3_q_sym2,-coef_f1_q_sym3,coef_f0_q_sym1-coef_f2_q_sym3,-coef_f3_q_sym3,-coef_f1_q_sym4,-coef_f2_q_sym4,coef_f0_q_sym1-coef_f3_q_sym4,coef_f0_q_sym5-coef_f1_q_sym12,-coef_f2_q_sym12,-coef_f3_q_sym12,coef_f0_q_sym6-coef_f1_q_sym13,coef_f0_q_sym5-coef_f2_q_sym13,-coef_f3_q_sym13,coef_f0_q_sym7-coef_f1_q_sym14,-coef_f2_q_sym14,coef_f0_q_sym5-coef_f3_q_sym14,coef_f0_q_sym8-coef_f1_q_sym15,coef_f0_q_sym6-coef_f2_q_sym15,-coef_f3_q_sym15,coef_f0_q_sym9-coef_f1_q_sym16,coef_f0_q_sym7-coef_f2_q_sym16,coef_f0_q_sym6-coef_f3_q_sym16,coef_f0_q_sym10-coef_f1_q_sym17,-coef_f2_q_sym17,coef_f0_q_sym7-coef_f3_q_sym17,-coef_f1_q_sym19,coef_f0_q_sym8-coef_f2_q_sym19,-coef_f3_q_sym19,-coef_f1_q_sym20,coef_f0_q_sym9-coef_f2_q_sym20,coef_f0_q_sym8-coef_f3_q_sym20,-coef_f1_q_sym21,coef_f0_q_sym10-coef_f2_q_sym21,coef_f0_q_sym9-coef_f3_q_sym21,-coef_f1_q_sym23,-coef_f2_q_sym23,coef_f0_q_sym10-coef_f3_q_sym23,coef_f0_q_sym12+coef_f1_q_sym5+t5+t8,coef_f2_q_sym5+t10,coef_f3_q_sym5+t12,coef_f0_q_sym13+coef_f1_q_sym6+t6,coef_f0_q_sym12+coef_f2_q_sym6+t5,coef_f3_q_sym6,coef_f0_q_sym14+coef_f1_q_sym7+t7,coef_f2_q_sym7,coef_f0_q_sym12+coef_f3_q_sym7+t5,coef_f0_q_sym15+coef_f1_q_sym5+coef_f1_q_sym8+t5+t9,coef_f0_q_sym13+coef_f2_q_sym5+coef_f2_q_sym8+t6+t11,coef_f3_q_sym5+coef_f3_q_sym8+t13,coef_f0_q_sym16+coef_f1_q_sym9,coef_f0_q_sym14+coef_f2_q_sym9+t7,coef_f0_q_sym13+coef_f3_q_sym9+t6,coef_f1_q_sym5+coef_f0_q_sym17+coef_f1_q_sym10+t5+t9,coef_f2_q_sym5+coef_f2_q_sym10+t11,coef_f0_q_sym14+coef_f3_q_sym5+coef_f3_q_sym10+t7+t13,coef_f1_q_sym6+coef_f0_q_sym19+t6,coef_f0_q_sym15+coef_f2_q_sym6+t5,coef_f3_q_sym6,coef_f1_q_sym7+coef_f0_q_sym20+t7,coef_f0_q_sym16+coef_f2_q_sym7,coef_f0_q_sym15+coef_f3_q_sym7+t5,coef_f1_q_sym6+coef_f0_q_sym21+t6,coef_f0_q_sym17+coef_f2_q_sym6+t5,coef_f0_q_sym16+coef_f3_q_sym6,coef_f1_q_sym7+coef_f0_q_sym23+t7,coef_f2_q_sym7,coef_f0_q_sym17+coef_f3_q_sym7+t5,coef_f1_q_sym8+t8,coef_f0_q_sym19+coef_f2_q_sym8+t6+t10,coef_f3_q_sym8+t12,coef_f1_q_sym9,coef_f0_q_sym20+coef_f2_q_sym9+t7,coef_f0_q_sym19+coef_f3_q_sym9+t6,coef_f1_q_sym8+coef_f1_q_sym10+t9,coef_f0_q_sym21+coef_f2_q_sym8+coef_f2_q_sym10+t6+t11,coef_f0_q_sym20+coef_f3_q_sym8+coef_f3_q_sym10+t7+t13,coef_f1_q_sym9,coef_f0_q_sym23+coef_f2_q_sym9+t7,coef_f0_q_sym21+coef_f3_q_sym9+t6,coef_f1_q_sym10+t8,coef_f2_q_sym10+t10,coef_f0_q_sym23+coef_f3_q_sym10+t7+t12],[3,28]);

author: zarathustr

repo: LibQPEP


text:
% EPOCH = DATENUM2EPOCH(TIME) translates Matlab numeric
% date format into Unix epoch date format.

% Author: Michael West. Modified Glenn Thompson.
% $Date$
% $Revision$

function epoch = datenum2epoch(time)

epoch = zeros(size(time));

if admin.antelope_exists()
    for n = 1:numel(time)
        epoch(n) = str2epoch(datestr(time(n),'mm/dd/yyyy HH:MM:SS.FFF'));
    end
else
    % Does not account for leap seconds.
    for n = 1:numel(time)
        epoch(n) = (time - datenum(1970,1,1,0,0,0))*86400;
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function [ n_data, x, bip ] = airy_bi_prime_values ( n_data )

%*****************************************************************************80
%
%% AIRY_BI_PRIME_VALUES returns some values of the Airy function Bi'(x).
%
%  Discussion:
%
%    The Airy functions Ai(X) and Bi(X) are a pair of linearly independent
%    solutions of the differential equation:
%
%      W'' - X * W = 0
%
%    In Mathematica, the function can be evaluated by:
%
%      AiryBiPrime[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real BIP, the derivative of the Airy BI function.
%
  n_max = 11;

  bip_vec = [ ...
     0.4482883573538264E+00, ...
     0.4515126311496465E+00, ...
     0.4617892843621509E+00, ...
     0.4800490287524480E+00, ...
     0.5072816760506224E+00, ...
     0.5445725641405923E+00, ...
     0.5931444786342857E+00, ...
     0.6544059191721400E+00, ...
     0.7300069016152518E+00, ...
     0.8219038903072090E+00, ...
     0.9324359333927756E+00 ];

  x_vec = [ ...
     0.0E+00, ...
     0.1E+00, ...
     0.2E+00, ...
     0.3E+00, ...
     0.4E+00, ...
     0.5E+00, ...
     0.6E+00, ...
     0.7E+00, ...
     0.8E+00, ...
     0.9E+00, ...
     1.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    bip = 0.0;
  else
    x = x_vec(n_data);
    bip = bip_vec(n_data);
  end

  return
end

author: johannesgerer


text:
function combo_test11 ( )

%*****************************************************************************80
%
%% COMBO_TEST11 tests KNAPSACK_REORDER and KNAPSACK_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  mass_limit = 26.0;
  p = [ 24.0, 13.0, 23.0, 15.0, 16.0 ]';
  w = [ 12.0,  7.0, 11.0,  8.0,  9.0 ]';

  fprintf ( 1, ' \n' );
  fprintf ( 1, 'COMBO_TEST11\n' );
  fprintf ( 1, '  KNAPSACK_REORDER reorders the knapsack data.\n' );
  fprintf ( 1, '  KNAPSACK_01 solves the 0/1 knapsack problem.\n' );

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Object, Profit, Mass, "Profit Density"\n' );
  fprintf ( 1, ' \n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %7.3f  %7.3f  %7.3f\n', i, p(i), w(i), p(i) / w(i) );
  end

  [ p, w ] = knapsack_reorder ( n, p, w );

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  After reordering by Profit Density:\n' );
  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Object, Profit, Mass, "Profit Density"\n' );
  fprintf ( 1, ' \n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %7.3f  %7.3f  %7.3f\n', i, p(i), w(i), p(i) / w(i) );
  end

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Total mass restriction is %f\n', mass_limit );

  [ x, mass, profit ] = knapsack_01 ( n, mass_limit, p, w );

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Object, Density, Choice, Profit, Mass\n' );
  fprintf ( 1, ' \n' );
  for i = 1 : n
    fprintf ( 1, '%6d  %7.3f  %7.3f  %7.3f  %7.3f\n', ...
      i, p(i) / w(i), x(i),  x(i) * p(i), x(i) * w(i) );
  end

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Total:                  %7.3f  %7.3f\n', profit, mass );

  return
end