text:
function h = max( f, g, dim )
%MAX   Maximum value of a SEPARABLEAPPROX in one direction.
%   MAX(f) returns a chebfun representing the maximum of the SEPARABLEAPPROX along the
%   y direction, i.e, MAX(f) = @(x) max( f ( x, : ) )
%
%   MAX(f, [], dim) returns a CHEBFUN representing the maximum of f along the
%   DIM direction. If DIM = 1 is along the y-direction and DIM = 2 is along the
%   x-direction.
%
%   WARNING: This function is not always accurate to the expected precision.
%
%   For the global maximum use MAX2.
%
% See also MIN, MAX2, MIN2, MINANDMAX2.

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

% Empty check:
if ( isempty( f ) )
    error('CHEBFUN:SEPARABLEAPPROX:max:input', 'SEPARABLEAPPROX is empty');
end

% Default to max of one separableApprox:
if ( nargin < 2 )
    g = [];
end

% Default to maximum along the y direction:
if ( nargin < 3 )
    dim = 1;
end

% Do not allow max(F, G):
if ( nargin > 1 && ~isempty( g ) )
    error('CHEBFUN:SEPARABLEAPPROX:max:twoSeparableApproxInputs', ...
        'Unable to maximize two SEPARABLEAPPROX objects.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We have no idea how to achieve this in an efficient way. This
% is an attempt to return a result, but typically it won't be accurate to
% more than 4-5 digits.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dom = f.domain;
n = 2049;
if ( dim == 1 )
    vals = sample(f, n, n);
    h = chebfun( max( vals ).', dom(1:2), 'splitting', 'on' );
    h = simplify( h.', [], 'globaltol' );
elseif ( dim == 2 )
    vals = sample(f, n, n);
    h = chebfun( max( vals, [], 2 ), dom(3:4), 'splitting', 'on' );
    h = simplify( h, [], 'globaltol' );
elseif ( dim == 0 )
    error('CHEBFUN:SEPARABLEAPPROX:max:dim', ...
        'Dimension argument must be a positive integer scalar within indexing range.')
else
   % return the SEPARABLEAPPROX. This is analogous to that MAX() command in
   % MATLAB.
   h = f;
end

end


text:
function c = correlation_besselj ( n, rho, rho0 )

%*****************************************************************************80
%
%% CORRELATION_BESSELJ evaluates the Bessel J correlation function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Petter Abrahamsen,
%    A Review of Gaussian Random Fields and Correlation Functions,
%    Norwegian Computing Center, 1997.
%
%  Parameters:
%
%    Input, integer N, the number of arguments.
%
%    Input, real RHO(N,1), the arguments.
%
%    Input, real RHO0, the correlation length.
%
%    Output, real C(N,1), the correlations.
%
  rho = rho ( : );

  rhohat = abs ( rho ) / rho0;

  c = besselj ( 0, rhohat );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% This function select utterance names from ID cell array that also appears
% in key cell array. It also select the corresponding features in data
% arrya.
% The function assumes key is already sorted.
% [ID_selected, data_selected] = selectData(key, ID, data, FollowKeyIdx)
%
function [ID_selected, data_selected] = selectData(key, ID, data, FollowKeyIdx)
if nargin<4
    FollowKeyIdx = 0;
end

if FollowKeyIdx==0
    ID_selected = {};
    data_selected = {};
end
for i=1:length(ID)
    found = searchInCellBinary(key, ID{i}, 1, length(key));
    if found > 0
        if FollowKeyIdx
            ID_selected{found} = ID{i};
            data_selected{found} = data{i};
        else
            ID_selected{end+1} = ID{i};
            data_selected{end+1} = data{i};
        end
    end
end

author: singaxiong

repo: SignalGraph


text:
function [dphi,dlambda,h] = togeod(a,finv,X,Y,Z)
%TOGEOD   Subroutine to calculate geodetic coordinates
%         latitude, longitude, height given Cartesian
%         coordinates X,Y,Z, and reference ellipsoid
%         values semi-major axis (a) and the inverse
%         of flattening (finv).

%  The units of linear parameters X,Y,Z,a must all agree (m,km,mi,ft,..etc)
%  The output units of angular quantities will be in decimal degrees
%  (15.5 degrees not 15 deg 30 min).  The output units of h will be the
%  same as the units of X,Y,Z,a.

%  Copyright (C) 1987 C. Goad, Columbus, Ohio
%  Reprinted with permission of author, 1996
%  Fortran code translated into MATLAB
%  Kai Borre 03-30-96
% Changed according to Matlab ver. 6.5.1, 9 February 2004

h = 0;
tolsq = 1.e-10;
maxit = 1000;
% compute radians-to-degree factor
rtd = 180/pi;
% compute square of eccentricity
if finv < 1.e-20
   esq = 0;
else
   esq = (2-1/finv)/finv;
end
oneesq = 1-esq;
% first guess
% P is distance from spin axix
P = sqrt(X^2+Y^2);
% direct calculation of longitude
if P > 1.e-20
   dlambda = atan2(Y,X)*rtd;
else
   dlambda = 0;
end;
if (dlambda < 0)
   dlambda = dlambda + 360;
end
% r is distance from origin (0,0,0)
r = sqrt(P^2+Z^2);
if r > 1.e-20
   sinphi = Z/r;
else
   sinphi = 0;
end
dphi = asin(sinphi);
% initial value of height  =  distance from origin minus
% approximate distance from origin to surface of ellipsoid
if r < 1.e-20
   h = 0;
   return;
end;
h = r-a*(1-sinphi*sinphi/finv);
% iterate
for i = 1:maxit
   sinphi = sin(dphi);
   cosphi = cos(dphi);
   % compute radius of curvature in prime vertical direction
   N_phi = a/sqrt(1-esq*sinphi*sinphi);
   % compute residuals in P and Z
   dP = P - (N_phi + h) * cosphi;
   dZ = Z - (N_phi*oneesq + h) * sinphi;
   % update height and latitude
   h = h+(sinphi*dZ+cosphi*dP);
   dphi = dphi+(cosphi*dZ-sinphi*dP)/(N_phi + h);
   % test for convergence
   if (dP*dP + dZ*dZ < tolsq)
      break;
   end

   % Not Converged--Warn user
   if i == maxit
      fprintf([' Problem in TOGEOD, did not converge in %2.0f',...
            ' iterations\n'],i)
   end;
end;
dphi = dphi*rtd;
%%%%%%%% end togeod.m  %%%%%%%%%%%%%%%%%%%%%%

author: yandld


text:
function eventLog = ma_executeNBodyCoast(coastEvent, initialState, eventNum, celBodyData)
%ma_executeNBodyCoast Summary of this function goes here
%   Detailed explanation goes here

	soiSkipIds = coastEvent.soiSkipIds;
    refBody = coastEvent.refBody;
    massLoss = coastEvent.massloss;
    forceModel = coastEvent.forceModel;
	maxPropTime = coastEvent.maxPropTime;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Get global coast ODE solver event functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    events = ma_getGlobalEvents(initialState, celBodyData);

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Propagate through the n revs, if applicable
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    eventLog = [];
    revs = coastEvent.revs;
    if(revs > 0)
        rVect = initialState(2:4);
        vVect = initialState(5:7);

        bodyID = initialState(8);
        bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);

        [sma, ecc, ~, ~, ~, tru] = getKeplerFromState(rVect, vVect, bodyInfo.gm);

        if(ecc < 1)
            revsInitstate = initialState;
            for(i=1:revs) %#ok<*NO4LP>
                period = computePeriod(sma, bodyInfo.gm);
                eventLogCoastRev1 = ma_executeCoast_nBody_goto_dt(period/2, revsInitstate, eventNum, forceModel, true, soiSkipIds, massLoss, maxPropTime, events, true, celBodyData);
                eventLogCoastRev2 = ma_executeCoast_nBody_goto_tru(tru, eventLogCoastRev1(end,:), eventNum, forceModel, true, soiSkipIds, refBody, massLoss, maxPropTime, events, celBodyData);

                eventLog = [eventLog;eventLogCoastRev1; eventLogCoastRev2]; %#ok<AGROW>
                revsInitstate = eventLogCoastRev2(end,:);
            end

            coastINIState = eventLogCoastRev2(end,:);
        else
            coastINIState = initialState;
            warnStr = ['Could not perform ',num2str(revs),' revolutions during coast, initial orbit is not elliptical.'];
            addToExecutionWarnings(warnStr, eventNum, bodyID, celBodyData);
        end
    else
        coastINIState = initialState;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Do the coast according to its sub-type
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    type = coastEvent.coastType;
    value = coastEvent.coastToValue;
    switch type
        case 'goto_ut'
            eventLogCoast = ma_executeCoast_nBody_goto_ut(value, coastINIState, eventNum, forceModel, true, soiSkipIds, massLoss, maxPropTime, events, celBodyData);

        case 'goto_dt'
            eventLogCoast = ma_executeCoast_nBody_goto_dt(value, coastINIState, eventNum, forceModel, true, soiSkipIds, massLoss, maxPropTime, events, true, celBodyData);

        case 'goto_tru'
            eventLogCoast = ma_executeCoast_nBody_goto_tru(value, coastINIState, eventNum, forceModel, true, soiSkipIds, refBody, massLoss, maxPropTime, events, celBodyData);

        case 'goto_apo'
            eventLogCoast = ma_executeCoast_nBody_goto_tru(pi, coastINIState, eventNum, forceModel, true, soiSkipIds, refBody, massLoss, maxPropTime, events, celBodyData);

        case 'goto_peri'
            eventLogCoast = ma_executeCoast_nBody_goto_tru(0, coastINIState, eventNum, forceModel, true, soiSkipIds, refBody, massLoss, maxPropTime, events, celBodyData);

        case 'goto_asc_node'
            eventLogCoast = ma_executeCoast_nBody_goto_node('asc', coastINIState, eventNum, forceModel, true, soiSkipIds, refBody, massLoss, maxPropTime, events, celBodyData);

        case 'goto_desc_node'
            eventLogCoast = ma_executeCoast_nBody_goto_node('desc', coastINIState, eventNum, forceModel, true, soiSkipIds, refBody, massLoss, maxPropTime, events, celBodyData);

        case 'goto_soi_trans'
            eventLogCoast = ma_executeCoast_nBody_goto_soi_trans(coastINIState, eventNum, forceModel, soiSkipIds, massLoss, maxPropTime, events, celBodyData);

        otherwise
            error(['Did not recongize coast of type ', type]);
    end

    eventLog = [eventLog; eventLogCoast];
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  BPSO and VPSO source codes version 1.0                           %
%                                                                   %
%  Developed in MATLAB R2011b(7.13)                                 %
%                                                                   %
%  Author and programmer: Seyedali Mirjalili                        %
%                                                                   %
%         e-Mail: ali.mirjalili@gmail.com                           %
%                 seyedali.mirjalili@griffithuni.edu.au             %
%                                                                   %
%       Homepage: http://www.alimirjalili.com                       %
%                                                                   %
%   Main paper: S. Mirjalili and A. Lewis, "S-shaped versus         %
%               V-shaped transfer functions for binary Particle     %
%               Swarm Optimization," Swarm and Evolutionary         %
%               Computation, vol. 9, pp. 1-14, 2013.                %
%                                                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [gBest,gBestScore,ConvergenceCurve]=BPSO(noP,Max_iteration,BPSO_num,CostFunction,noV)

%Initial Parameters for PSO
w=2;              %Inirtia weight
wMax=0.9;         %Max inirtia weight
wMin=0.4;         %Min inirtia weight
c1=2;
c2=2;
Vmax=6;

Velocity=zeros(noP,noV);%Velocity vector
Position=zeros(noP,noV);%Position vector

%////////Cognitive component/////////
pBestScore=zeros(noP);
pBest=zeros(noP,noV);
%////////////////////////////////////

%////////Social component///////////
gBestScore=inf;
gBest=zeros(1,noV);
%///////////////////////////////////

ConvergenceCurve=zeros(1,Max_iteration); %Convergence vector

%Initialization
for i=1:size(Position,1) % For each particle
    for j=1:size(Position,2) % For each variable
        if rand<=0.5
            Position(i,j)=0;
        else
            Position(i,j)=1;
        end
    end
end
for i=1:noP
    pBestScore(i)=inf;
end

for l=1:Max_iteration

    %Calculate cost for each particle
    for i=1:size(Position,1)
        fitness=CostFunction(Position(i,:));

        if(pBestScore(i)>fitness)
            pBestScore(i)=fitness;
            pBest(i,:)=Position(i,:);
        end
        if(gBestScore>fitness)
            gBestScore=fitness;
            gBest=Position(i,:);
        end
    end

    %update the W of PSO
    w=wMax-l*((wMax-wMin)/Max_iteration);
    %Update the Velocity and Position of particles
    for i=1:size(Position,1)
        for j=1:size(Position,2)
            %Equation (1)
            Velocity(i,j)=w*Velocity(i,j)+c1*rand()*(pBest(i,j)-Position(i,j))+c2*rand()*(gBest(j)-Position(i,j));

            if(Velocity(i,j)>Vmax)
                Velocity(i,j)=Vmax;
            end
            if(Velocity(i,j)<-Vmax)
                Velocity(i,j)=-Vmax;
            end

            if BPSO_num==1
                s=1/(1+exp(-2*Velocity(i,j))); %S1 transfer function
            end
            if BPSO_num==2
                s=1/(1+exp(-Velocity(i,j)));   %S2 transfer function
            end
            if BPSO_num==3
                s=1/(1+exp(-Velocity(i,j)/2)); %S3 transfer function
            end
            if BPSO_num==4
               s=1/(1+exp(-Velocity(i,j)/3));  %S4 transfer function
            end

            if BPSO_num<=4 %S-shaped transfer functions
                if rand<s % Equation (4) and (8)
                    Position(i,j)=1;
                else
                    Position(i,j)=0;
                end
            end

            if BPSO_num==5
                s=abs(erf(((sqrt(pi)/2)*Velocity(i,j)))); %V1 transfer function
            end
            if BPSO_num==6
                s=abs(tanh(Velocity(i,j))); %V2 transfer function
            end
            if BPSO_num==7
                s=abs(Velocity(i,j)/sqrt((1+Velocity(i,j)^2))); %V3 transfer function
            end
            if BPSO_num==8
                s=abs((2/pi)*atan((pi/2)*Velocity(i,j))); %V4 transfer function (VPSO)
            end

            if BPSO_num>4 && BPSO_num<=8 %V-shaped transfer functions
                if rand<s %Equation (10)
                    Position(i,j)=~Position(i,j);
                end
            end
        end
    end
    ConvergenceCurve(l)=gBestScore;
end
end

author: Sable

repo: mcbench-benchmarks


text:
function Population = Truncation(Population,N)
% Limit the size of final popualtion in RVEA*

%------------------------------- 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".
%--------------------------------------------------------------------------

    Choose = true(1,length(Population));
    Cosine = 1 - pdist2(Population.objs,Population.objs,'cosine');
    Cosine(logical(eye(length(Cosine)))) = 0;
    while sum(Choose) > N
        Remain   = find(Choose);
        Temp     = sort(-Cosine(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Choose(Remain(Rank(1))) = false;
    end
    Population = Population(Choose);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Gcone_twv_test.m
% test the Gcone object for "thread within view" version by Brian Gonzalez
% Copyright 2015-07-20, Jeff Fessler, University of Michigan

pn = jf_protected_names;

systypes = 'sf2,v:all';
systypes = {'sf2', systypes}
nsys = numel(systypes);

f.is_zxy = true;
if f.is_zxy
	to_zxy = @(x) permute(x, [3 1 2]); to_xyz = @(z) permute(z, [2 3 1]);
else
	to_zxy = @(x) x;
end

% small systems for basic tests
if 1 && (0 || ~isvar('Ac')), printm 'setup small'
	f.down = 16;

	dfs_list = [0 inf inf]; % arc flat parallel
	dsd_list = [949.075 949.075 inf]; % arc flat parallel

	f.dy = 0.7; % stress test SF with non-square
	igs = image_geom('nx', 512, 'ny', 480, 'nz', 416, ...
		'dx', 1, 'dy', f.dy, 'dz', 0.5, ...
		'offset_x', 2.9, 'offset_y', 3.5, ...
		'offset_z', -3.4, ...
		'mask', 'all-but-edge', ... % trick: for hct2
		'down', f.down);

	for kk = 1:numel(dfs_list)

		cgs = ct_geom('ge2', 'nt', 320, ...
			'source_z0', -20, 'pitch', 0.5, ... % test helix
			'dfs', dfs_list(kk), ... % arc or flat
			'dsd', dsd_list(kk), ... % fan or parallel beam
			'down', f.down);

		for ii=1:nsys
			systype = systypes{ii};
			printm('testing type %s dfs=%g dsd=%g', ...
				systype, cgs.dfs, cgs.dsd)

			args = {cgs, igs, 'type', systype};
			args = {args{:}, 'zxy', f.is_zxy};

			Ac = Gcone(args{:});

		%	im(Ac * to_zxy(single(igs.mask)))
		%	im(to_xyz(Ac' * cgs.ones)), prompt

			tester_tomo2(Ac, to_zxy(igs.mask)) % paces
			test_adjoint(Ac, 'big', 1, 'tol', 5e-5)
		end % systype
	end % dfs
end % small

if ~isvar('x0'), printm 'x0 big'
	f.down = 4;
	igb = image_geom('nx', 512, 'ny', 480, 'nz', 416, ...
		'dx', 1, 'dz', 0.5, ...
		'dy', 0.7, ... % stress test
		'offset_x', 12.9, 'offset_y', 3.5, 'offset_z', -3.4, ...
		'down', f.down);
	ell = [3*igb.dx 5*igb.dx -2*igb.dz ...
		igb.dx*igb.nx/3 igb.dy*igb.ny/4 igb.zfov/4 ...
		0 0 10];
	x0 = ellipsoid_im(igb, ell, 'oversample', 2);
end

% big systems for accuracy tests
if 0 || ~isvar('Ab'), printm 'setup big'
	cgb = ct_geom('ge1', 'nt', 320, ...
		'source_z0', -40, 'pitch', 0.5, ... % test helix
		'down', f.down);
%		'dfs', inf, ... % flat detector
%		'dsd', inf, 'dfs', inf, ... % parallel beam

	for ii=1:nsys
		systype = systypes{ii};
		Ab{ii} = Gcone(cgb, igb, 'type', systype, 'zxy', f.is_zxy);
	end
end

if ~isvar('ya'), printm 'analytical projections'
	ya = ellipsoid_proj(cgb, ell, 'oversample', 2);
%	im clf, im(ya)
%prompt
end

if 0 || ~isvar('yb'), printm 'discrete projections'
	nrmse = @(x,y) norm(y(:)-x(:)) / norm(x(:)) * 100;
	x00 = to_zxy(x0);
	for ii=1:nsys
		cpu etic
		yb{ii} = Ab{ii} * x00;
		f.time(ii) = cpu('etoc');
		printm('%s: time %g nrmse %g %%', ...
			systypes{ii}, f.time(ii), nrmse(ya, yb{ii}))
	end
end

if 0
	im_toggle(ya(:,end/2,:), yb{1}(:,end/2,:), ...
		yb{1}(:,end/2,:) - ya(:,end/2,:))
	im(ya(:,end/2,:) - yb{1}(:,end/2,:))
	im(yb{1}(:,end/2,:) - yb{2}(:,end/2,:))
end

if 1, printm 'look at error in worst views'
	im('pl', 2, nsys+1)
	for ii=1:numel(systypes)
		err = yb{ii} - ya;
		tmp = reshape(err, [], cgb.na);
		tmp = sqrt(mean(tmp.^2)); % error in each view
		ia = imax(tmp); % worst view
		im(ii, err(:,:,ia)), cbar h
		titlef('%s ia=%d', systypes{ii}, ia)
		im('subplot', ii+1+nsys)
		if im
			plot(1:cgb.na, tmp), axis tight
			titlef('%s RMS', systypes{ii})
			xlabelf('view')
		end
	end
	if 1
		err = yb{2} - yb{1};
		tmp = reshape(err, [], cgb.na);
		tmp = sqrt(mean(tmp.^2)); % error in each view
		ia = imax(tmp); % worst view
		im(nsys+1, err(:,:,ia)), cbar h
		titlef('%s vs %s', systypes{1}, systypes{2})
	end
return
end

if 1, printm 'projection profiles'
	it = cgb.nt;
	it = round(cgb.nt/2); % middle
	it = it + [-2 0 2];
	ia = imin(abs(cgb.ad - 45));
%	ia = ceil(ia/2);
	pro = @(y) col(y(:,it,ia));
	arg = [];
	for ii=1:numel(systypes)
		arg = [arg pro(yb{ii})];
	end
	if im
		clf, plot([arg pro(ya)])
		text(10, 200, sprintf('ia=%d', ia))
		text(10, 400, sprintf('ang=%g', cgb.ad(ia)))
		legend(systypes{:}, 'true')
		axisy(0, 1.2 * max(ya(:)))
		grid
	end
return
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function nlp = imposeNLPConstraint(obj, nlp, ep, nzy, load_path)
    % impose virtual constraints as NLP constraints in the trajectory
    % optimization problem 'nlp' of the dynamical system
    %
    % For each (vector) virtual constraints, we will enforce the following
    % constraints at the first node:
    % {y, y', ..., y^N} % where N = RelativeDegree - 1
    % and this constraint:
    % y^(N+1) + ep_N y^N + ... ep_1 y = 0
    % at all nodes.
    %
    % Parameters:
    % nlp: the trajectory optimization NLP @type TrajectoryOptimization
    % ep: the coefficients of the derivatives @type rowvec
    % nzy: the derivatives that does not requires to be zero at the first
    % node @type rowvec
    %
    % @note For example, ep = [kp, kd] for typical vector relative degree 2
    % holonomic virtual constraints.
    %
    % @note For example, nzy = [0] for a relative degree 1 nonholonomic
    % virtual constraint (velocity outputs) means we do not require y to
    % be zero at the first node. Another example, nzy = [0,0] for a
    % relative degree 2 virtual constraints, we do not require both y and
    % y' to be zero at the first node, on the other hand, nzy = [1,0] for a
    % relative degree 2 virtual constraints, we only enforce y being zero
    % at the first node.

    if nargin < 5
        load_path = [];
    end

    % local variables for speed
    rel_deg = obj.RelativeDegree;
    is_holonomic = obj.Holonomic;
    is_state_based = strcmp(obj.PhaseType, 'StateBased');
    model = obj.Model;
    ya = obj.ActualFuncs;
    yd = obj.DesiredFuncs;

    assert(~isempty(ya) || ~isempty(yd),...
        'The virtual constraint functions are not compiled yet. Please run compile method first.');

    %% validate input arguments
    validateattributes(nlp, {'NonlinearProgram'},...
        {'scalar'},'VirtualConstraint.imposeNLPConstraint',...
        'nlp');
    nlpOptions = nlp.Options;
    n_node     = nlp.NumNode;
    %%
    if strcmp(nlpOptions.CollocationScheme,'PseudoSpectral')
        error('The PseudoSpectral is incompatible with virtual constraints.');
    end

    if nargin > 2
        validateattributes(ep, {'double'},...
            {'vector','numel',rel_deg,'positive','real'},...
            'VirtualConstraint.imposeNLPConstraint','ep');
    else
        warning('The coefficient vector (ep) of outputs are not defined. Use 1 for all derivatives.');
        ep = ones(1,rel_deg);
        for l= 1:rel_deg
            ep(l) = nchoosek(rel_deg,l-1)*10^(rel_deg - l + 1);
        end
    end
    if nargin > 3
        validateattributes(nzy, {'double'},...
            {'vector','numel',rel_deg,'binary'},...
            'VirtualConstraint.imposeNLPConstraint','nzy');
    else
        nzy = ones(1,rel_deg);

    end

    % the name suffix of functions
    name = [obj.Name '_' model.Name];

    % desired output parameters and its name in the var table of the NLP
    if isa(obj.OutputParams,'SymVariable')
        a_var = {SymVariable(tomatrix(obj.OutputParams(:)))};
        a_name = obj.OutputParamName;
    else
        a_var = {};
        a_name = {};
    end

    % phase variable parameters and its name in the var table of the NLP
    if isa(obj.PhaseParams,'SymVariable')
        p_var = {SymVariable(tomatrix(obj.PhaseParams(:)))};
        p_name = obj.PhaseParamName;
    else
        p_var = {};
        p_name = {};
    end

    % offset variable parameters and its name in the var table of the NLP
    if isa(obj.OffsetParams,'SymVariable')
        c_var = {SymVariable(tomatrix(obj.OffsetParams(:)))};
        c_name = obj.OffsetParamName;
    else
        c_var = {};
        c_name = {};
    end

    % states and their name representation in the var table of the NLP
    switch model.Type
        case 'FirstOrder'
            x_var = {model.States.x};
            dx_var = {model.States.dx};
            x_name = {'x'};
            dx_name = {'dx'};

            dX = model.States.dx;
        case 'SecondOrder'
            x_var = {model.States.x, model.States.dx};
            dx_var = {model.States.ddx};
            x_name = {'x','dx'};
            dx_name = {'ddx'};

            dX = [model.States.dx; model.States.ddx];
    end

    if is_holonomic
        q_var = {model.States.x};
        q_name = {'x'};
    else
        q_var = {model.States.x,model.States.dx};
        q_name = {'x','dx'};
    end

    % if the desired outputs are time-based, get the function for
    % converting the horizon time T to the actual node time t
    if is_state_based
        t_var = {};
        t_name = {};
        aux_var = {};
    else
        t = SymVariable('t');
        k = SymVariable('k');
        T  = SymVariable('t',[2,1]);
        nNode = SymVariable('nNode');
        tsubs = T(1) + ((k-1)./(nNode-1)).*(T(2)-T(1));

        if ~isnan(nlpOptions.ConstantTimeHorizon)
            t_var = {};
            t_name = {};
            aux_var = {T,k,nNode};
        else
            t_var = {T};
            t_name = 'T';
            aux_var = {k,nNode};
        end
    end

    y_fun = cell(rel_deg+1,1);

    %% y(x,a,p) = 0
    if nzy(1)==1
        if isempty(load_path)
            y = ya{1} - yd{1};
            if ~is_state_based
                %% Time-based outputs, need to incoorporates the time variable
                y = subs(y,t,tsubs);
            end

            % holonomic virtual constraints
            y_fun{1} = SymFunction(['y_' name], y, [t_var, q_var, a_var, p_var, c_var], aux_var);
        else
            % holonomic virtual constraints
            y_fun{1} = SymFunction(['y_' name], [], [t_var, q_var, a_var, p_var, c_var], aux_var);
            y_fun{1} = load(y_fun{1},load_path);
        end
        if ~isempty(aux_var)
            switch length(aux_var)
                case 2
                    % k, nNode
                    aux_data = {1,n_node};

                case 3
                    % T, k, nNode
                    aux_data = {nlpOptions.ConstantTimeHorizon, 1,n_node};

            end
        else
            aux_data = [];
        end
        % add constraint at the first node
        nlp = addNodeConstraint(nlp, y_fun{1}, [t_name, q_name ,a_name, p_name, c_name], 'first',...
            0, 0, 'Nonlinear', aux_data);

    end

    %% higher order derivatives
    if rel_deg > 1

        for i=2:rel_deg
            if nzy(i) == 1
                %% state-based output, no need to use time variable
                if isempty(load_path)
                    dy = ya{i} - yd{i};
                    if ~is_state_based
                        %% Time-based outputs, need to incoorporates the time variable
                        dy = subs(dy,t,tsubs);
                    end
                    y_fun{i} = SymFunction(['d' num2str(i-1) 'y_' name], dy, [t_var, x_var, a_var, p_var], aux_var);
                else
                    y_fun{i} = SymFunction(['d' num2str(i-1) 'y_' name], [], [t_var, x_var, a_var, p_var], aux_var);
                    y_fun{i} = load(y_fun{i}, load_path);
                end
                % add constraint at the first node

                nlp = addNodeConstraint(nlp, y_fun{i}, [t_name, x_name ,a_name, p_name], 'first',...
                    0, 0, 'Nonlinear', aux_data);
            end
        end
    end

    %% the highest order derivatives imposed at all nodes (feedback linearization)
    node_list = 1:1:n_node;
    dim = obj.Dimension;
    vars   = nlp.OptVarTable;
    ceq_err_bound = nlpOptions.EqualityConstraintBoundary;

    y_dynamics(1,n_node) = NlpFunction();

    if isempty(load_path)
        % state-based output, no need to use time variable
        if is_state_based
            ddy = (ya{rel_deg+1} - yd{rel_deg+1})*dX;
        else
            ddy = ya{rel_deg+1}*dX - yd{rel_deg+1};
        end

        %         ep_s = SymVariable('k',[rel_deg,1]);

        %         for j=1:1:rel_deg
        %             ddy = ddy + ep_s(j)*(ya{j} - yd{j});
        %         end
        % Time-based outputs, need to incoorporates the time variable
        if ~is_state_based
            ddy = subs(ddy,t,tsubs);
        end

        %         ddy_fun = SymFunction(['d' num2str(rel_deg) 'y_' name], ddy, [t_var, x_var, dx_var, a_var, p_var, c_var], [aux_var,{ep_s}]);
        ddy_fun = SymFunction(['d' num2str(rel_deg) 'y_' name], ddy, [t_var, x_var, dx_var, a_var, p_var, c_var], [aux_var]);
    else
        %         ep_s = SymVariable('k',[rel_deg,1]);

        %         ddy_fun = SymFunction(['d' num2str(rel_deg) 'y_' name], [], [t_var, x_var, dx_var, a_var, p_var, c_var], [aux_var,{ep_s}]);
        ddy_fun = SymFunction(['d' num2str(rel_deg) 'y_' name], [], [t_var, x_var, dx_var, a_var, p_var, c_var], [aux_var]);
        ddy_fun = load(ddy_fun, load_path);

    end

    for i=node_list
        idx = node_list(i);
        if nlpOptions.DistributeParameters
            param_index = idx;
        else
            param_index = 1;
        end

        if ~isempty(a_name)
            a_deps = vars.(a_name)(param_index);
        else
            a_deps = {};
        end
        if ~isempty(p_name)
            p_deps = vars.(p_name)(param_index);
        else
            p_deps = {};
        end
        if ~isempty(c_name)
            c_deps = vars.(c_name)(param_index);
        else
            c_deps = {};
        end

        if nlpOptions.DistributeTimeVariable
            time_index = idx;
        else
            time_index = 1;
        end

        if ~isempty(t_name)
            t_deps = vars.(t_name)(time_index);
        else
            t_deps = {};
        end
        if ~isempty(aux_var)
            switch length(aux_var)
                case 2
                    % k, nNode
                    aux_data = {idx,n_node};

                case 3
                    % T, k, nNode
                    aux_data = {nlpOptions.ConstantTimeHorizon, idx,n_node};

            end
        else
            aux_data = [];
        end
        x_deps = cellfun(@(x)vars.(x)(idx),x_name,'UniformOutput',false);
        dx_deps = cellfun(@(x)vars.(x)(idx),dx_name,'UniformOutput',false);

        y_dynamics(i) = NlpFunction('Name',[obj.Name '_output_dynamics'],...
            'Dimension',dim,'SymFun',ddy_fun,'lb',-ceq_err_bound,...
            'ub',ceq_err_bound,'Type','Nonlinear',...
            'DepVariables',[t_deps, x_deps{:}, dx_deps{:}, a_deps, p_deps, c_deps]',...
            'AuxData', {aux_data});
        %             'AuxData', {[aux_data,{ep}]});

    end

    obj.OutputFuncs = [y_fun;{ddy_fun}];
    %     obj.OutputFuncsName_ = cellfun(@(f)f.Name, obj.OutputFuncs,'UniformOutput',false);

    % add output dynamics at all nodes
    % add dynamical equation constraints
    nlp = addConstraint(nlp,[obj.Name '_output_dynamics'],'all',y_dynamics);

end

author: ayonga

repo: frost-dev

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


text:
function scaled_regions = map_boxes_to_regions_mult_scales(region_params, bboxes, image_size, scales)
% map_boxes_to_regions_mult_scales: given a set of bounding boxes, it
% creates for each of them a region (acording to the region_params provided)
% and it maps each of those regions from the original image with size
% image_size to the appropriate image scale from set of scales given from
% the parameter scales.
%
% 1) region_params: (type struct) the region pooling parameters. Some of
% its fields are:
%   a) sz_conv_standard: (scalar value) the last convolution size
%   b) step_standard: (scalar value) is the stride in pixels of the output
%   of the last convolutional layer of the network where the regions are
%   going to be projected (for VGG16 is equal to 16).
% 2) boxes: a N x 4 array with the bounding box coordinates in the form of
% [x0,y0,x1,y1] (where (x0,y0) is the top-left corner and (x1,y1) the
% bottom left corner)
% 3) image_size: 2 x 1 or 1 x 2 array with the size of the original image
% 4) scales: NS x 1 or 1 x NS array with the image scales that are used. NS
% is the number of images. The i-th value is the size in pixels of the
% smallest dimension of the image in the i-th scale.
%
% OUTPUT:
% 1) scaled_regions: is a N x 9 array with the N output regions. Each region is
% represented by 9 values [scale_id, xo0, yo0, xo1, yo1, xi0, yi0, xi1, yi1]
% that correspond to its outer rectangle [xo0, yo0, xo1, yo1] and its inner
% rectangle [xi0, yi0, xi1, yi1]. scale_id corresponds to the scale id to
% which each region is mapped.
%
% This file is part of the code that implements the following paper:
% Title      : "LocNet: Improving Localization Accuracy for Object Detection"
% Authors    : Spyros Gidaris, Nikos Komodakis
% Institution: Universite Paris Est, Ecole des Ponts ParisTech
% ArXiv link : http://arxiv.org/abs/1511.07763
% code       : https://github.com/gidariss/LocNet
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2016 Spyros Gidaris
%
% Title     : "LocNet: Improving Localization Accuracy for Object Detection"
% ArXiv link: http://arxiv.org/abs/1511.07763
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------
regions = create_regions_from_boxes( region_params, bboxes );
[ scaled_regions ] = map_regions_to_mult_scales(region_params, regions, image_size, scales );
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function y = mean_nan(x)

count = sum(~isnan(x));
x(find(isnan(x)))=0;
y = sum(x)./count;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypemethod  @@sym {} max (@var{a})
%% @deftypemethodx @@sym {} max (@var{a}, @var{b})
%% @deftypemethodx @@sym {} max (@var{a}, [], @var{dim})
%% @deftypemethodx @@sym {[@var{r}, @var{I}] =} max (@dots{})
%% Return maximum value of a symbolic vector or vectors.
%%
%% Example:
%% @example
%% @group
%% max(sym(1), sym(2))
%%   @result{} (sym) 2
%% max([1 2*sym(pi) 6])
%%   @result{} (sym) 2⋅π
%% [M, I] = max([1 2*sym(pi) 6])
%%   @result{} M = (sym) 2⋅π
%%   @result{} I = 2
%% @end group
%% @end example
%%
%% @seealso{@@sym/min}
%% @end deftypemethod

function [z, I] = max(A, B, dim)

  if (nargout <= 1)
    if (nargin == 1)
      if (isvector(A))
        z = pycall_sympy__ ('return Max(*_ins[0])', A);
      else
        z = max(A, [], 1);
      end
    elseif (nargin == 2)
      z = elementwise_op ('Max', sym(A), sym(B));
    elseif (nargin == 3)
      assert (isempty (B))
      assert (logical(dim == 1) || logical(dim == 2))

      cmd = { '(A, dim) = _ins'
              'if not A.is_Matrix:'
              '    A = sp.Matrix([A])'
              'if dim == 0:'
              '    if A.rows == 0:'
              '        return A'
              '    return Matrix([[Max(*A.col(i)) for i in range(0, A.cols)]])'
              'elif dim == 1:'
              '    if A.cols == 0:'
              '        return A'
              '    return Matrix([Max(*A.row(i)) for i in range(0, A.rows)])' };
      z = pycall_sympy__ (cmd, A, dim - 1);
    else
      print_usage ();
    end
    return
  end

  % dealing with the index (2nd output) is complicated, defer to min
  if (nargin == 1)
    [z, I] = min(-A);
    z = -z;
  elseif (nargin == 3)
    [z, I] = min(-A, -B, dim);
    z = -z;
  else
    print_usage ();
  end

end

%% many other tests are in @sym/min

%!test
%! % simple
%! assert (isequal (max([sym(10) sym(11)]), sym(11)))

%!test
%! syms x y
%! assert (isequal (children (max (x, y)), [x y]))

repo: octsympy


text:
function [d,F,G]=principalAngVec(A,B)
%%PRINCANGPRINCVEC Given an mXp matrix and an mXq matrix with p>=q, compute
%           the principal angles and vectors between the matrices. This can
%           be used in canonical correlation analysis to find combinations
%           of variables with maximum correlation between each other.
%           Specifically, the principle angles are defined recursively by
%           d(k)=cos(theta_k)=F(:,k)'*G(:,k)=
%               max_{F(:,k)} max_{G(:,k)} F(:,k)'*G(:,k)
%           where in the optimization F(:,k)  is a unit vector in the space
%           spanned by A that is orthogonal to all F(:,i) with i<k, and
%           G(:,k) is similarly a unit vector in the space spanned by B
%           that is orthogonal to all G(,i) with i<k.
%
%INPUTS: A An mXp matrix with independent columns.
%        B An mXq matrix with independent columns q<=p.
%
%OUTPUTS: d A qX1 vector of the cosine values of the principle angles (Use
%           acos to get the angles). The angles all satisfy 0<=theta<=pi/2.
%         F mXq The first set of principal vectors.
%         G mXq The second set of principal vectors.
%
%This function implements Algorithm 6.4.3 in Chapter 6.4 of [1].
%
%EXAMPLE:
%This function can be used to quickly tell if two matrices are identicle
%except for a reordering and scaling of the columns. If they are, then the
%principal angles will all be zero (the cosines will all be 1). We
%demonstrate this here:
% A=[8,0,1;
%    9,2,9;
%    1,5,9;
%    9,9,4;
%    6,9,8];
% B=-4*A(:,[3,1,2]);%Scale and permute the columns.
% [d,F,G]=principalAngVec(A,B)
%One will see that the three d values are all 1 within finite precision
%bounds.
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

q=size(B,2);

if(q>size(A,2))
    error('The number of columns of B must be <= the number of columns of A.')
end

[QA,~]=qr(A,0);
[QB,~]=qr(B,0);
[Y,D,Z]=svd(QA'*QB);
F=QA*Y(:,1:q);
G=QB*Z(:,1:q);
d=diag(D);

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function vbrfa2010_testbed_experiment(D, common_nu, common_tau, maxiter)

testbed_experiment_vbrfa(D, common_nu, common_tau, maxiter, 'robustness', ...
                         'independent-t');
testbed_experiment_vbrfa(D, common_nu, common_tau, maxiter, 'robustness', ...
                         'multivariate-t');
testbed_experiment_vbrfa(D, common_nu, common_tau, maxiter, 'robustness', ...
                         'none');

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [ filters ] = func_fbcsp( Dat, varargin )
%FUNC_FBCSP Summary of this function goes here
%   Detailed explanation goes here
opt=opt_cellToStruct(varargin{:});

if isstruct(Dat)
    if isfield(Dat, 'fs')
        opt.fs=Dat.fs;
    elseif sifield(opt, 'fs')

    else
        error('OpenBMI: "fs" is missing')
    end
else
    if ~isfield(Dat, 'fs')
        error('OpenBMI: "fs" is missing')
    end
end

if iscell(Dat)  % cell type binary classes
    if length(Dat)==2
        C1=prep_selectClass(Dat, {'class', opt.classes{1}});
        C2=prep_selectClass(Dat, {'class', opt.classes{2}});
    else
        error('OpenBMI: The classes should be binary');
    end
else
    if isfield(opt, 'classes')
        if length(opt.classes) ~=2
            error('OpenBMI: The classes should be binary');
        else
            C1=prep_selectClass(Dat, {'class', opt.classes{1}});
            C2=prep_selectClass(Dat, {'class', opt.classes{2}});
        end
    end
end

%% Specral filtering
[nDat, nTri, nCH] = size( C1.x ); % same size of [dat, trial, ch] for each class
X1 = reshape( C1.x, [nDat*nTri,nCH ] );
X2 = reshape( C2.x, [nDat*nTri,nCH ] );

x_flt=cell(2, length(opt.Filters)) % binary,
for i=1:length(opt.Filters)
    x_flt{1, i}= prep_filter(X1,{'frequency', opt.Filters(i,:); 'fs', opt.fs});
    x_flt{2, i}= prep_filter(X2,{'frequency', opt.Filters(i,:); 'fs', opt.fs});
end

%% Spatial filtering
CSPFilter = cell( 1, length(opt.Filters) );
for i=1:length(opt.Filters)
    X1 = x_flt{1, i};
    X2 = x_flt{2, i};

    %% Here is basic CSP filtering.
    S1 = cov(X1(:,:));
    S2 = cov(X2(:,:));
    [W,D] = eig(S1, S1+S2);
    %%
    CSPFilter{i}.W = W( :, [1:opt.numCSPPatterns, end - opt.numCSPPatterns+1:end] );
    Dd = diag(D);
    CSPFilter{i}.D = Dd([1:opt.numCSPPatterns, end - opt.numCSPPatterns+1:end]);
end

%% Feature extraction

features=cell(2, length(opt.Filters))
for i=1:length(opt.Filters)
    for j=1:2 % for two classes
        dat=func_projection(x_flt{j,i}, CSPFilter{i}.W);
        dat=reshape( dat, [nDat,nTri,size(dat,2)]);
        features{j,i} = squeeze( log(var(dat, 0, 1)) ); % log-variance feature
    end
end

%% Mutual information
kernelWidth=1;
for i=1:length(opt.Filters)
    miValue(i) = proc_mutual_information( features{1, i}, features{2, i}, kernelWidth );
end

end

author: PatternRecognition

repo: OpenBMI


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

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

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

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

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

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

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function ui = mrViewOverlay(ui);
%
%  ui = mrViewOverlay(ui);
%
% Compute overlays for each image,  combine with
% colormap and threshold settings to produce a
% set of True-Color images to display in
% ui.display.images.
%
% ras 07/05.
if ~exist('ui', 'var') | isempty(ui),  ui = mrViewGet; end
if ishandle(ui),     ui = get(ui, 'UserData');         end

%%%% get needed params
s = ui.settings.space;
b = ui.settings.bounds;
slices = ui.display.slices;
oris = ui.display.oris;
method = ui.settings.mapInterp;

% preference for how to combine multiple overlays:
% 1) average two colors: [1 0 0] + [.2 0 1] => [.6 0 .5];
% 2) add and saturate: [1 0 0] + [.2 0 1] => [1 0 1];
% 3) opaque: [1 0 0] + [.2 0 1] => [.2 0 1]; % last one on top
if ispref('VISTA', 'alphaMethod')
    alphaMethod = getpref('VISTA', 'alphaMethod');
else
    alphaMethod = 1; % average
end

%%%%% get a list of which overlays to show (ignoring ones that are
%%%%% hidden
isHidden = [ui.overlays.hide];
overlayList = find( ~isHidden );
if isempty(overlayList),  return; end

%%%%% loop across overlays,  computing / displaying each in turn
for i = 1:length(ui.display.images)

    imgSz = size(ui.display.images{i});

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % compute each overlay image,  mask of where to show it %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    overlay = []; mask = logical([]);
    for o = overlayList
        ii = find(overlayList==o);
        m = ui.overlays(o).mapNum;

        %%%%% get map volume for this timepoint
        vol = ui.maps(m).data;
        if ndims(vol)>3
            t = ui.overlays(o).subVol;
            vol = vol(:,:,:,t);
        end

        % get the coords for this slice in the base data space
        imgCoords = mrViewGet(ui, 'DataCoords', i);

        % restrict coords to those which pass the threshold
        % of this overlay
        [coords, imageMask] = mrViewRestrict(ui, imgCoords, o, method);

        % keep track of the pixels which pass threshold
        mask(:,:,ii) = logical(reshape(imageMask, imgSz));

        % xform coords into map data coordinates
        C = inv(ui.maps(m).baseXform) * [coords; ones(1, size(coords, 2))];

        % get the values for the overlay map
        sz = size(vol);
        if(method(1)=='n')
            vals = myCinterp3(vol,[sz(1) sz(2)], sz(3), round(C([2,1,3],:)'), 0.0);
        elseif(method(1)=='l')
            vals = myCinterp3(vol,[sz(1) sz(2)], sz(3), C([2,1,3],:)', 0.0);
        else
            vals = interp3(vol, C(2,:), C(1,:), C(3,:), method);
        end

        % store vals in overlay matrix
        overlaySlice = zeros(imgSz);
        overlaySlice(imageMask) = vals;
        overlay(:,:,ii) = overlaySlice;

        % compute clip vals for overlay
        if isequal(ui.overlays(o).clim, 'auto'),  clim = [min(vals) max(vals)];
        else,                                     clim = ui.overlays(o).clim;
        end

        % update colorbar's color limits, label, units
        ui.overlays(o).cbar.clim = clim;
		ui.overlays(o).cbar.label = ui.maps(m).name;
		ui.overlays(o).cbar.units = ui.maps(m).dataUnits;
    end

    % set NaN values to 0
    mask(isnan(overlay)) = false;
    overlay(isnan(overlay)) = 0;

    % compute # of overlays present @ each pixel (for transparency)
    overlaysPerPixel = sum(double(mask), 3);

    % find indices of pixels w/ any overlays
    ok = find(overlaysPerPixel>0);

    % initialize those parts of image to be black:

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plug in overlays,  taking into account overlap %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % initialize underlay image: black out parts where
    % overlays will be,  make truecolor
    img = rescale(ui.display.images{i}, ui.settings.clim, [1 256]);
    img(ok) = 0; img = ind2rgb(img, ui.settings.cmap);

    for o = overlayList
        ii = find(overlayList==o);
        nColors = size(ui.overlays(o).cbar.cmap, 1);

        % get pixel mask for this overlay
        M = mask(:,:,ii);

        if any(M(:)>0)
            % use overlay data to form an index into the color map:
            vals = overlay(:,:,ii);
            vals = vals(M);
            I = rescale(vals, ui.overlays(o).cbar.clim, [1 nColors]);

            % manually plug values in to each color channel of image
            % (use the alphaMethod pref to determine how to combine
            % multiple overlays)
            for ch = 1:3
                tmp = img(:,:,ch);
                switch alphaMethod
                    case 1,   % standard alpha -- average colors
                        W = 1 ./ overlaysPerPixel(M>0); % weights
                        tmp(M) = tmp(M) + ui.overlays(o).cbar.cmap(I,ch) .* W;
                    case 2,   % add channels and saturate
                        tmp(M) = tmp(M) + ui.overlays(o).cbar.cmap(I,ch);
                        tmp(M) = min(tmp(M), 1);
                    case 3,   % opaque overlays
                        tmp(M) = ui.overlays(o).cbar.cmap(I,ch);
                end
                img(:,:,ch) = tmp;
            end
        end
    end

    ui.display.images{i} = img;
    clear overlay mask
end

return

% % attempt to use alpha map for display
%         % set so that other overlays can peek through
%         mask = uint8(255 .* (1/length(ui.overlays)) .* mask);
%         % combine color map and overlay values to make
%         % truecolor image
%         img = ind2rgb(img, ui.overlays(o).cmap);
%        % display image
%         otherDims = setdiff(1:3, ui.display.oris(i));
%         xx = b(otherDims(2), 1):b(otherDims(2), 2);
%         yy = b(otherDims(1), 1):b(otherDims(1), 2);
%         h = image(xx, yy, img, 'AlphaData', mask, ...
%             'AlphaDataMapping', 'none', ...
%             'Parent', ui.display.axes(i));
%         ui.overlays(o).imgHandle(i) = h;

% old comments for above:
% Compute and display overlay images in a mrViewer UI.
%
% NOTE: The method used to display overlay images involves
% using matlab alpha channels,  which makes it more difficult
% to export these images outside matlab. Moreover,  to overlay
% multiple maps which each have different color maps,  ok've found
% it necessary to convert each image to a true color image.
% Because of the increased memory requirements -- 3 color planes
% plus an alpha mask per overlay image -- ok've decided to
% combine the computation and display steps into one function.
%
% If greater use is made of 'hidden' views,  which don't use
% the mrViewer UI but use the structs,  then this should probably
% be rewritten to keep them separate. As well,  it would be probably
% be better to do away with the use of alpha channels,  and just
% produce a single truecolor image for underlay + all overlays.
% Failing that,  it's probably simpler and more elegant to use
% alpha channels.
%


text:
function [ ishow, list_num, list, nregion ] = voxels_region_3d ( maxlist, ...
  nx, ny, nz, ishow )

%*****************************************************************************80
%
%% VOXELS_REGION_3D arranges a set of voxels into contiguous regions in 3D.
%
%  Discussion:
%
%    On input, the ISHOW array contains zero and nonzero values.  The nonzero
%    values are taken to be active voxels.  On output, the zero voxels remain
%    zero, and all the active voxels have been assigned a value which now
%    indicates membership in a region, or group of contiguous voxels.
%
%    On output, the array LIST contains information about the regions.
%    The last used element of LIST is LIST_NUM.
%
%    The number of elements in region NREGION is NELEM = LIST(LIST_NUM).
%    The (I,J,K) indices of the last element in this region are in
%    LIST(LIST_NUM-3) through LIST(LIST_NUM-1), and the first element is
%    listed in LIST(LIST_NUM-3*NELEM), LIST(LIST_NUM-3*NELEM+1),
%    LIST(LIST_NUM-3*NELEM+2).
%
%    The number of elements in NREGION-1 is listed in LIST(LIST_NUM-3*NELEM-1),
%    and the (I,J,K) indices of the these elements are listed there.
%
%  Picture:
%
%    Input:
%
%      0  2  0  0 17  0  3
%      0  0  3  0  1  0  4
%      1  0  4  8  8  0  7
%      3  0  6 45  0  0  0
%      3 17  0  5  9  2  5
%
%    Output:
%
%      0  1  0  0  2  0  3
%      0  0  2  0  2  0  3
%      4  0  2  2  2  0  3
%      4  0  2  2  0  0  0
%      4  4  0  2  2  2  2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer MAXLIST, the maximum length of the array used to
%    list the elements of the regions.
%
%    Input, integer NX, NY, NZ, the number of voxels in the X, Y and
%    Z directions.
%
%    Input, integer ISHOW(NX,NY,NZ), the only significance to
%    the entries is whether they are zero or nonzero.
%
%    Output, integer ISHOW(NX,NY,NZ), the nonzero entries have now been revalued
%    so that contiguous entries have the same value, indicating a grouping into
%    a region.
%
%    Output, integer LIST_NUM, the number of entries of LIST that were used.
%    However, if MAXLIST < LIST_NUM, then there was not enough space in
%    LIST to store the data properly, and LIST should not be used,
%    although the data in ISHOW should be correct.
%
%    Output, integer LIST(MAXLIST), contains, in stack form, a list
%    of the indices of the elements in each region.
%
%    Output, integer NREGION, the number of regions discovered.
%

%
%  Reset all nonzero entries of ISHOW to -1.
%
  for i = 1 : nx
    for j = 1 : ny
      for k = 1 : nz
        if ( ishow(i,j,k) ~= 0 )
          ishow(i,j,k) = -1;
        end
      end
    end
  end
%
%  Start the number of items in the region list at 0.
%
  list_num = 0;
%
%  Start the number of regions at 0.
%
  nregion = 0;
%
%  The stack begins empty.
%
  nstack = 0;
%
%  Search for an unused "ON" voxel from which we can "grow" a new region.
%
  for i = 1 : nx
    for j = 1 : ny
      for k = 1 : nz
%
%  We found a voxel that is "ON", and does not belong to any region.
%
        if ( ishow(i,j,k) == -1 )
%
%  Increase the number of regions.
%
          nregion = nregion + 1;
%
%  Add this voxel to the region.
%
          ishow(i,j,k) = nregion;
%
%  Add this voxel to the stack.
%
          stack(nstack+1) = i;
          stack(nstack+2) = j;
          stack(nstack+3) = k;

          stack(nstack+4) = 1;

          nstack = nstack + 4;
%
%  Add this voxel to the description of the region.
%
          nelements = 1;

          if ( list_num + 3 <= maxlist )
            list(list_num+1) = i;
            list(list_num+2) = j;
            list(list_num+3) = k;
          end

          list_num = list_num + 3;

          while ( 1 )
%
%  Find all neighbors of BASE that are "ON" but unused.
%  Mark them as belonging to this region, and stack their indices.
%
            ibase = stack(nstack-3);
            jbase = stack(nstack-2);
            kbase = stack(nstack-1);

            ilo = max ( ibase-1, 1 );
            ihi = min ( ibase+1, nx );
            jlo = max ( jbase-1, 1 );
            jhi = min ( jbase+1, ny );
            klo = max ( kbase-1, 1 );
            khi = min ( kbase+1, nz );

            nabes = 0;

            for i2 = ilo : ihi
              for j2 = jlo : jhi
                for k2 = klo : khi
%
%  We found a neighbor to our current search point, which is "ON" and unused.
%
                  if ( ishow(i2,j2,k2) == -1 )
%
%  Increase the number of neighbors.
%
                    nabes = nabes + 1;
%
%  Mark the neighbor as belonging to the region.
%
                    ishow(i2,j2,k2) = nregion;
%
%  Add the neighbor to the stack.
%
                    stack(nstack+1) = i2;
                    stack(nstack+2) = j2;
                    stack(nstack+3) = k2;

                    nstack = nstack + 3;
%
%  Add the neighbor to the description of the region.
%
                    nelements = nelements + 1;

                    if ( list_num + 3 <= maxlist )
                      list(list_num+1) = i2;
                      list(list_num+2) = j2;
                      list(list_num+3) = k2;
                    end

                    list_num = list_num + 3;

                  end
                end
              end
            end
%
%  If any new neighbors were found, take the last one as the basis
%  for a deeper search.
%
            if ( 0 < nabes )
              stack(nstack+1) = nabes;
              nstack = nstack + 1;
              continue
            end
%
%  If the current search point had no new neighbors, drop it from the stack.
%
            ncan = stack(nstack) - 1;
            nstack = nstack - 3;
            stack(nstack) = ncan;
%
%  If there are still any unused candidates at this level, take the
%  last one as the basis for a deeper search.
%
            if ( 0 < stack(nstack) )
              continue
            end
%
%  If there are no more unused candidates at this level, then we need
%  to back up a level in the stack.  If there are any candidates at
%  that earlier level, then we can still do more searching.
%
            nstack = nstack - 1;

            if ( nstack <= 0 )
              break
            end

          end
%
%  If we have exhausted the stack, we have completed this region.
%  Tag the number of elements to the end of the region description list.
%
          list_num = list_num + 1;
          if ( list_num <= maxlist )
            list(list_num) = nelements;
          end

        end

      end
    end
  end
%
%  Print some warnings.
%
  if ( maxlist < list_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'VOXELS_REGION - Warning!\n' );
    fprintf ( 1, '  MAXLIST was too small to list the regions.\n' );
    fprintf ( 1, '  Do not try to use the LIST array!\n' );
    fprintf ( 1, '  The ISHOW data is OK, however.\n' );
  end

  return
end


text:
% SP_TO_VTK: Export to VTK format for plotting.
%
%  sp_to_vtk (u, space, geometry, npts, filename, fieldname, [options], [lambda_lame, mu_lame])
%  sp_to_vtk (u, space, geometry, pts, filename, fieldname, [options], [lambda_lame, mu_lame])
%
% INPUT:
%
%     u:           vector of dof weights
%     space:       object representing the space of discrete functions (see sp_vector)
%     geometry:    geometry structure (see geo_load)
%     npts:        number of points along each parametric direction where to evaluate
%     pts:         cell array with the coordinates along each parametric direction of the points where to evaluate
%     filename:    name of the output file.
%     fieldnames:  how to name the saved variables in the vtk file
%     options:     cell array with the fields to plot
%                   accepted options are 'value' (default), 'gradient',
%                   and for vectors also 'curl', 'divergence', 'stress'
%     lambda_lame: function handle to the first Lame coefficient (only needed to compute 'stress')
%     mu_lame:     function handle for the second Lame coefficient (only needed to compute 'stress')
%
% OUTPUT:
%
%    none
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2012, 2015 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function sp_to_vtk (u, space, geometry, npts, filename, fieldname, varargin)

  [eu, F] = sp_eval (u, space, geometry, npts, varargin{:});

  msh_to_vtk (F, eu, filename, fieldname);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function post_3d ( p, t )

%*****************************************************************************80
%
%% POST_3D is a postprocessor for the tetrahedronal mesh.
%
%  Copyright:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Modified:
%
%    31 July 2009
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, P(N,3), the coordinates of nodes.
%
%    Input, integer T(NT,4), the indices of nodes that make up the tetrahedrons.
%
  [ node_num, junk ] = size ( p );
  [ tetra_num, junk ] = size ( t );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  %d nodes were used.\n', node_num );
  fprintf ( 1, '  %d tetrahedrons were created.\n', tetra_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  QUALITY MEASURES:\n');
  fprintf ( 1, '  A value of 1.00 is best, and 0.00 is the worst.\n' );

  q = simp_qual_3d ( p, t, 1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Measure Min         Mean        Max         Var\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  #1  %10f  %10f  %10f  %10f\n', ...
    min ( q ), mean ( q ), max ( q ), var ( q ) );

  bad = find ( q == min ( q ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The tetrahedron of minimum quality is number %d\n', bad(1) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Node        X             Y             Z\n' );
  fprintf ( 1, '\n' );

  for j = 1 : 4
    n = t(bad(1),j);
    fprintf ( 1, '  %4d  %12f  %12f  %12f\n',  n, p(n,1:3) );
  end

  return
end


text:
function [ChannelLabels, Time, Frequency, Data, Info] = read_besa_tfc(FILENAME)

% READ_BESA_TFC imports data from a BESA *.tfc file
%
% Use as
%   [DataType, ConditionName, Channels, Time, Frequency, Data] = read_besa_tfc(FILENAME)
%
% This reads data from the BESA Time-Frequency-Coherence output data file
% FILENAME and returns the following data:
%   ConditionName: name of analyzed condition
%   ChannelLabels: character array of channel labels
%   Time: array of sampled time instants
%   Frequency: array of sampled frequencies
%   Data: 3D data matrix with indices (channel,time,frequency)
%   Info: Struct containing additional information:
%       DataType: type of the exported data
%       ConditionName: name of analyzed condition
%       NumbeOfTrials: Number of trials on which the data is based
%       StatisticsCorrection: Type of statistics correction for multiple testing
%       EvokedSignalSubtraction: Type of evoked signal subtraction

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

fp = fopen(FILENAME);

VersionNumber = fscanf(fp,'VersionNumber=%s ',1);
DataType = fscanf(fp,'DataType=%s ',1);
ConditionName = fscanf(fp,'ConditionName=%s ',1);
try
    if ConditionName == 'Condition'
        ConditionName = [ConditionName,' ',fscanf(fp,'%s ',1)];
    end
catch
end
NumberTrials = fscanf(fp,'NumberTrials=%i ',1);
NumberTimeSamples = fscanf(fp,'NumberTimeSamples=%i ');
TimeStartInMS = fscanf(fp,'TimeStartInMS=%f ',1);
TimeIntervalInMS = fscanf(fp,'IntervalInMS=%f ',1);
NumberFrequencies = fscanf(fp,'NumberFrequencies=%i ');
FreqStartInHZ = fscanf(fp,'FreqStartInHz=%f ',1);
FreqIntervalInHZ = fscanf(fp,'FreqIntervalInHz=%f ',1);
NumberChannels = fscanf(fp,'NumberChannels=%i ');

% New file versions (BESA 5.0.8 and higher) include more information in the .tfc file header; skip that
vers=0;                 % new tfc file format
try
    StatisticsCorrection = fscanf(fp,'StatisticsCorrection=%s ',1);
    EvokedSignalSubtraction = fscanf(fp,'EvokedSignalSubtraction=%s',1);
catch
    vers=1;             % old tfc file format
end

% Handle possible future extensions of the tfc file header
i=1;
while i<1000
    a = fscanf(fp,'%c',1);
    if strcmp(a,sprintf('\n'))
        i=1000;
    end
    i=i+1;
end

% Generate return values
% FIXME the following statement does not work for MATLAB 7.2 on XP (see mail from Stephan Bickel)
Time = [TimeStartInMS:TimeIntervalInMS:(NumberTimeSamples-1)*TimeIntervalInMS+TimeStartInMS];
Frequency = [FreqStartInHZ:FreqIntervalInHZ:(NumberFrequencies-1)*FreqIntervalInHZ+FreqStartInHZ];
if vers == 1
    Info = struct('DataType',{DataType},'ConditionName',{ConditionName},'NumberOfTrials',{NumberTrials});
else
    Info = struct('DataType',{DataType},'ConditionName',{ConditionName},'NumberOfTrials',{NumberTrials},...
        'StatisticsCorrection',{StatisticsCorrection},'EvokedSignalSubtraction',{EvokedSignalSubtraction});
end

if ~contains(DataType,'COH')
    for Channel=1:NumberChannels
        ChannelLabels(Channel) = cellstr(fscanf(fp,'%s ',1));
    end
else
    for Channel=1:NumberChannels
        ChannelLabels(Channel) = cellstr(fscanf(fp,'%s ',3));
    end
end

ChannelLabels=char(ChannelLabels);

Data = zeros(NumberChannels,NumberTimeSamples,NumberFrequencies);
for Channel=1:NumberChannels
    Data(Channel,:,:) = fscanf(fp,'%f',[NumberTimeSamples,NumberFrequencies]);
end
fclose(fp);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Z = tenfun(fun,varargin)
%TENFUN Apply a function to each element in a tensor.
%
%   TENFUN(F,X,...) applies the function specified by the function
%   handle F to the given arguments.  Either both arguments
%   must be tensors, or one is a tensor and the other is a scalar/MDA.
%
%   Examples
%   Z = tenfun(@(x)(x+1),X) %<-- increase every element by one
%   Z = tenfun(@eq,X,1) %<-- logical comparison of X with scalar
%   Z = tenfun(@plus,X,Y) %<-- adds the two tensors X and Y.
%   Z = tenfun(@max,X,Y,Z) %<-- max over all elements in X,Y,Z
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if nargin < 2
    %  error('TENFUN requires at least two input arguments')
    error('Not enough input arguments.');
end

if ~isa(fun, 'function_handle')
    error('First argument must be a function handle.');
end

%% Case I: NARGIN == 2 (one function and one tensor)
if (nargin == 2) && isa(varargin{1},'tensor')
    Z = varargin{1};
    Z.data(:) = fun(Z.data(:));
    return;
end

%% Determine if function is binary
% Note that we swap the arguments below if 2nd argument is sparse, but need
% to take special measures for those functions that don't commute.
binfuns = {@plus,@minus,@eq,@ge,@gt,@le,@lt,@ne,@and,@or,@xor, ...
           @power,@times,@ldivide,@rdivide};
isbinary = false;
for i = 1 : numel(binfuns)
    if isequal(fun,binfuns{i})
        isbinary = true;
        break;
    end
end

%% Case II: NARGIN == 3 and function is binary
if (nargin == 3) && isbinary

    X = varargin{1};
    Y = varargin{2};

    % Case IIa: X is a scalar/MDA and Y is a tensor
    if isnumeric(X) && isa(Y,'tensor')
        Z = Y;
        Z.data = fun(X,Y.data);
        return;
    end

    % Case IIb: X is a tensor and Y is a scalar/MDA
    if isa(X,'tensor') && isnumeric(Y)
        Z = X;
        Z.data = fun(X.data,Y);
        return;
    end

    % Case IIc: X and Y are both tensors
    if isa(X,'tensor') && isa(Y,'tensor')
        if ~(isequal(size(X),size(Y)))
            error('Tensor size mismatch.')
        end
        data = fun(X.data,Y.data);
        Z = tensor(data,size(X));
        return;
    end

    % Case IId: Either X or Y is a sptensor
    if isa(X,'tensor') && isa(Y,'sptensor')
        if isequal(fun,@lt)
            Z = gt(Y,X);
        elseif isequal(fun,@le)
            Z = ge(Y,X);
        elseif isequal(fun,@gt)
            Z = lt(Y,X);
        elseif isequal(fun,@ge)
            Z = le(Y,X);
        elseif isequal(fun,@minus)
            Z = plus(-Y,X);
        elseif isequal(fun,@power)
            error('Cannot do array power with a sparse and dense tensor');
        elseif isequal(fun,@ldivide)
            error('Cannot do ldivide with a sparse and dense tensor');
        elseif isequal(fun,@rdivide)
            error('Cannot do rdivide with a sparse and dense tensor');
        else
            Z = fun(Y,X);
        end
        return;
    end

    % Case IIe: Either X or Y is not a tensor
    error('For a binary function, arguments must be either two tensors or one tensor and a scalar/MDA.');
end

%% Case III: one function and the rest are same-sized tensors

X = varargin;
n = numel(X);
sz = size(X{1});
m = prod(sz);
Y = zeros(m,n);
for j = 1:n
    Y(:,j) = X{j}.data(:);
end
data = zeros(m,1);
for i = 1:m
    data(i) = fun(Y(i,:));
end
Z = tensor(data,sz);
return;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function refreshdisp(str,prevstr,iteration)

if ~exist('iteration','var')
    iteration=2;
end

if iteration==1
    fprintf(str)
else
    fprintf(char(8*ones(1,length(prevstr))));
    fprintf(str);
end

author: flatironinstitute

repo: NoRMCorre


text:
function [ferns,hsPr] = fernsClfRemoveTrainingData( hs, remove_data_idx, ferns )
% Compute cross-validation error of random ferns classifier
%
% See "Fast Keypoint Recognition in Ten Lines of Code" by Mustafa Ozuysal,
% Pascal Fua and Vincent Lepetit, CVPR07.
%
% Dimensions:
%  M - number ferns
%  S - fern depth
%  F - number features
%  N - number input vectors
%  H - number classes
%
% USAGE
%  [ferns,hsPr] = fernsClfTrain( data, hs, [varargin] )
%
% INPUTS
%  data     - [NxF] N length F feature vectors
%  hs       - [Nx1] target output labels in [1,H]
%  varargin - additional params (struct or name/value pairs)
%   .S        - [10] fern depth (ferns are exponential in S)
%   .M        - [50] number of ferns to train
%   .thrr     - [0 1] range for randomly generated thresholds
%   .bayes    - [1] if true combine probs using bayes assumption
%   .ferns    - [] if given reuse previous ferns (recompute pFern)
%   .cv_sets  - [] if given, N x 1 vector containing the index of each
%   .cv_nsets - [10] if cv_sets not specified, then the data is split into
%               cv_nsets randomly chosen sets.
%
% OUTPUTS
%  fracwrong_cv  - [1x1] fraction of training examples predicted
%               incorrectly in hsPr_cv
%  hsPr_cv    - [Nx1] predicted cross-validation-based classification of
%               each training example
%  probs_cv   - [NxH] predicted cross-validation-based output label
%               probabilities for each training example
%  ferns      - learned fern model from all training data with the
%               following fields
%   .fids     - [MxS] feature ids for each fern for each depth
%   .thrs     - [MxS] threshold corresponding to each fid
%   .pFern    - [2^SxHxM] learned log probs at fern leaves
%   .bayes    - if true combine probs using bayes assumption
%   .inds     - [NxM] cached indices for original training data
%   .H        - number classes
%
% EXAMPLE
%  N=5000; H=5; d=2; [xs0,hs0,xs1,hs1]=demoGenData(N,N,H,d,1,1);
%  fernPrm=struct('S',4,'M',50,'thrr',[-1 1],'bayes',1);
%  tic, [ferns,hsPr0]=fernsClfTrain(xs0,hs0,fernPrm); toc
%  tic, hsPr1 = fernsClfApply( xs1, ferns ); toc
%  e0=mean(hsPr0~=hs0); e1=mean(hsPr1~=hs1);
%  fprintf('errors trn=%f tst=%f\n',e0,e1); figure(1);
%  subplot(2,2,1); visualizeData(xs0,2,hs0);
%  subplot(2,2,2); visualizeData(xs0,2,hsPr0);
%  subplot(2,2,3); visualizeData(xs1,2,hs1);
%  subplot(2,2,4); visualizeData(xs1,2,hsPr1);
%
% See also fernsClfTrain, fernsClfApply, fernsInds
%
% Kristin Branson 2011-06-09, based on fernsClfTrain from:
%
% Piotr's Image&Video Toolbox      Version 2.50
% Copyright 2010 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Lesser GPL [see external/lgpl.txt]

% get additional parameters and check dimensions
%dfs={'S',10,'M',50,'thrr',[0 1],'bayes',1,'ferns',[],'cv_sets',[],'cv_nsets',10};
%[S,M,thrr,bayes,ferns,cv_sets,cv_nsets]=getPrmDflt(varargin,dfs,1);
H=max(hs);
[M,S] = size(ferns.fids);
assert(all(hs>0)); assert(S<=20);

% grab inds for data we are removing
inds_remove = ferns.inds(remove_data_idx,:);

% remove inds from ferns
ferns.inds(remove_data_idx,:) = [];

% update counts
pFern = ferns.counts;
edges = 1:2^S;
for m = 1:M,
  for h = 1:H,
    pFern(:,h,m) = pFern(:,h,m) - reshape(histc(inds_remove(hs(remove_data_idx)==h,m),edges),[2^S,1]);
  end
end
ferns.counts = pFern;

% update pFern
if( ferns.bayes<=0 )
  norm = 1./sum(pFern,2);
  % KB: use bsxfun instead of loop
  pFern = bsxfun(@times,pFern,norm);
  %for h=1:H, pFern(:,h,:)=pFern(:,h,:).*norm; end
else
  norm = 1./sum(pFern,1);
  % KB: use bsxfun instead of loop
  pFern = bsxfun(@times,pFern,norm);
  %for s=1:2^S, pFern(s,:,:)=pFern(s,:,:).*norm; end
  pFern=log(pFern);
end

% store pFern and compute output values
ferns.pFern=pFern; clear pFern;
if(nargout==2), hsPr=fernsClfApply([],ferns,ferns.inds); end

end

author: kristinbranson

repo: JAABA


text:
function X = tendiag(v,sz)
%TENDIAG Creates a tensor with v on the diagonal.
%
%   TENDIAG(V) creates a tensor with N dimensions, each of size N, where N
%   is the number of elements of V. The elements of V are placed on the
%   superdiagonal.
%
%   TENDIAG(V,SZ) is the same as above but creates a tensor of size SZ. If
%   SZ is not big enough, the tensor will be enlarged to accommodate the
%   elements of V on the superdiagonal.
%
%   Examples
%   X = tendiag([0.1 0.22 0.333]) %<-- creates a 3x3x3 tensor
%
%   See also TENSOR, SPTENDIAG.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

% Make sure v is a column vector
v = reshape(v,[numel(v) 1]);

N = numel(v);
if ~exist('sz','var')
    sz = repmat(N,1,N);
end

X = tenzeros(sz);
subs = repmat((1:N)', 1, length(sz));
X(subs) = v;