text:
 function pot = potential_func(type, delta, param)
%function pot = potential_func(type, delta, param)
%|
%| Define roughness penalty potential functions (using function handles).
%|
%| The penalty will have the form
%|	R(x) = sum_k w_k * potential_k([Cx]_k, delta_k)
%| where w_k is provided elsewhere, not here!
%|
%| in
%|	type		quad, huber, hyper2, hyper3, cauchy, lange1, lange3, ...
%|			recommended: 'hyper3'
%|	delta		scalar or image-sized array;
%|			"cutoff" parameter for edge-preserving regularization
%|	param		optional additional parameter(s) for some choices
%| out
%|	pot.delta	local copy
%|	anonymous functions:
%|	pot.potk(pot, C*x)	potential function value
%|	pot.wpot(pot, C*x)	potential 'weights' (aka half-quad. curvatures)
%|	pot.dpot(pot, C*x)	potential derivative
%|
% Copyright 2004-5-18, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if streq(type, 'test') potential_func_test, return, end

pot.delta = [];
pot.param = [];
if isvar('delta')
	pot.delta = delta(:);
	clear delta
end
if isvar('param')
	pot.param = param(:);
	clear param
end

% trick: huber2 is just a huber with delta / 2
% so that weighting function drops to 1/2 at delta, like hyper3 etc.
if streq(type, 'huber2')
	type = 'huber';
	pot.delta = pot.delta / 2;
end

% trick: hyper3 is just a hyperbola with delta scaled by sqrt(3)
% to approximately "match" the properties of 'cauchy' (old erroneous 'hyper')
if streq(type, 'hyper3')
	type = 'hyper2';
	pot.delta = pot.delta / sqrt(3);
end

switch type

%
% quadratic potential function
%
case 'quad'
	potk = '(abs(t).^2)/2';
	wpot = 'ones(size(t))';
	dpot = 't';

%
% huber potential function
%
case 'huber'
	potk = 'huber_pot(t, pot.delta)';
	wpot = 'huber_wpot(t, pot.delta)';
	dpot = 'huber_dpot(t, pot.delta)';

%
% cauchy penalty: d^2 / 2 * log(1 + (t/d)^2)
% Not convex!
%
case 'cauchy'
	potk = 'pot.delta.^2 / 2 .* log(1 + abs(t ./ pot.delta).^2)';
	wpot = '1 ./ (1 + abs(t ./ pot.delta).^2)';
	dpot = 't ./ (1 + abs(t ./ pot.delta).^2)';

%
% Geman&McClure penalty: d^2 / 2 * (t/d)^2 / (1 + (t/d)^2)
% Not convex!
%
case 'geman&mcclure'
	potk = 'pot.delta.^2 / 2 .* (t/pot.delta).^2 ./ (1 + abs(t ./ pot.delta).^2)';
	wpot = '1 ./ (1 + abs(t ./ pot.delta).^2).^2';
	dpot = 't ./ (1 + abs(t ./ pot.delta).^2).^2';

%
% hyperbola penalty: d^2 * [ sqrt(1 + (t/d)^2) - 1 ]
%
case 'hyper2'
	potk = 'pot.delta.^2 .* (sqrt(1 + abs(t ./ pot.delta).^2) - 1)';
	wpot = '1 ./ sqrt(1 + abs(t ./ pot.delta).^2)';
	dpot = 't ./ sqrt(1 + abs(t ./ pot.delta).^2)';

case 'hyper'
	error 'use "cauchy" or "hyper3" not "hyper" now'

%
% Lange1 penalty
%
case 'lange1'
	potk = 't.^2 / 2 ./ (1+abs(t./pot.delta))';
	wpot = '(1 + abs(t ./ pot.delta) / 2) ./ (1 + abs(t ./ pot.delta)).^2';
	dpot = ['t .* (' wpot ')'];

%
% Lange3 penalty
%
case 'lange3'
	potk = 'pot.delta.^2 .* (abs(t./pot.delta) - log(1+abs(t./pot.delta)))';
	wpot = '1 ./ (1 + abs(t ./ pot.delta))';
	dpot = 't ./ (1 + abs(t ./ pot.delta))';

%
% li98cfs
%
case 'li98cfs'
	% f = @(x) atan(x) / x - 0.5; fsolve(f, 2.3)
	pot.delta = pot.delta / 2.3311;
	potk = 'pot.delta.^2 .* ((t ./ pot.delta) .* atan(t ./ pot.delta) - 0.5 * log(1 + (t ./ pot.delta).^2))';
	wpot = 'ir_li98cfs_wpot(t, pot.delta)';
	dpot = ['t .* ' wpot];

%
% qgg2: q-generalized gaussian for p=2, due to Thibault, Sauer, Bouman
% q = "param", same as lange1 when q=1
%
case 'qgg2'
	potk = 't.^2 / 2 ./ (1+abs(t./pot.delta).^(2-pot.param))';
	wpot = ['(1 + abs(t ./ pot.delta).^(2-pot.param) * pot.param / 2) ' ...
		' ./ (1 + abs(t ./ pot.delta).^(2-pot.param)).^2'];
	dpot = ['t .* (' wpot ')'];

%
% genhub : "generalized" Huber, same as huber when p=2 and q=1
% p = param(1), q = param(2)
%
case 'genhub'
	potk = ['0.5 * abs(t) .^ p .* (abs(t) <= d) + ' ...
		'0.5 * (p ./ q .* d .^ (p-q) .* abs(t) .^ q + (1 - p ./ q) .* d .^ p) .* (abs(t) > d)'];
	reps = {'p', 'zp', 'q', 'zq', 'd', 'zd', ...
		'zp', 'pot.param(1)', 'zq', 'pot.param(2)', 'zd', 'pot.delta'};
	potk = strreps(potk, reps{:});
	wpot = ['p / 2 .* (abs(t) .^ (p-2)) .* (abs(t) <= d) + ' ...
		'p / 2 .* (d .^ (p-q)) .* ((abs(t) + (abs(t)==0)) .^ (q-2)) .* (abs(t) > d)']; % trick to avoid 0^-:
	wpot = strreps(wpot, reps{:});
	dpot = ['t .* (' wpot ')'];

%
% stevenson:94:dpr
% p = param(1), q = param(2), same as huber when p=2 and q=1 ???
%
case 'stevenson94dpr'
	potk = ['0.5 * abs(t) .^ p .* (abs(t) <= d) + ' ...
		'0.5 * ( (p .* d .^ (p-1) .* abs(t) - p .* d .^ p + (1 ./ q) .^ (1 ./ (q-1)) ) .^ q + d .^ p - (1 ./ q) .^ (q ./ (q-1)) ) .* (abs(t) > d)'];
	reps = {'p', 'zp', 'q', 'zq', 'd', 'zd', ...
		'zp', 'pot.param(1)', 'zq', 'pot.param(2)', 'zd', 'pot.delta'};
	potk = strreps(potk, reps{:});
	wpot = 'ones(size(t))'; % fix: fake for now
	dpot = ['t .* (' wpot ')'];

otherwise
	try % trick: try new strum version
		p = potential_fun(type, pot.delta, pot.param);
		pot.potk = @(dum, t) p.potk(t);
		pot.wpot = @(dum, t) p.wpot(t);
		pot.dpot = @(dum, t) p.dpot(t);
	catch
		fail('Unknown potential "%s"', type)
	end
end

if ~isfield(pot, 'potk')
	pot.potk = ir_str2func(['@(pot,t) ' potk]);
	pot.wpot = ir_str2func(['@(pot,t) ' wpot]);
	pot.dpot = ir_str2func(['@(pot,t) ' dpot]);
end

% test routine
% examine potential functions after rescaling.
function potential_func_test

delta = 10; tmax = 4;
plist = {'quad', 'huber2', 'hyper3', 'lange1', 'lange3', ...
	'cauchy', 'qgg2', 'gf1'};
%plist = {'quad', 'li98cfs', 'hyper3', 'huber2'}; % show li98cfs roughly hyper3
%plist = {'huber', 'genhub', 'quad'};%, 'stevenson94dpr'};
%plist = {'hyper3', 'qgg2'}; delta = 20; tmax = 10;
%plist = {'huber', 'hyper3', 'qgg2'}; delta = 10; tmax = 8;
%plist = {'qgg2', 'gf1'};
t = tmax * linspace(-delta, delta, 401)';
for ii=1:length(plist)
	type = plist{ii};
	leg{ii} = [plist{ii} ' \delta = ' num2str(delta)];
	if streq(type, 'qgg2')
		param = 1.2;
		leg{ii} = [leg{ii} ' q = ' num2str(param)];
	elseif streq(type, 'gf1')
		param = [0.0558 1.6395];
	elseif streq(type, 'genhub')
		param = [1.9 1.1];
		leg{ii} = [leg{ii} sprintf('p=%g q=%g', param(1), param(2))];
	elseif streq(type, 'stevenson94dpr')
		param = [2 2.01];
		leg{ii} = [leg{ii} sprintf('p=%g q=%g', param(1), param(2))];
	else
		param = [];
	end
	pot = potential_func(type, delta, param);
	pp(:,ii) = pot.potk(pot, t);
	pw(:,ii) = pot.wpot(pot, t);
	pd(:,ii) = pot.dpot(pot, t);
end

if im
	clf
	subplot(411), plot(t, pp), title 'potk'
	axis tight
	axisy([-0.0 2.5] * delta^2)
	legend(leg, 'location', 'north')
	subplot(412), plot(t, pw), title 'wpot'
	axisy(0, 1.1)
	subplot(413), plot(t, pd), title 'dpot'
	axisy([-1 1] * 1.1 * delta)
	% check derivatives
	subplot(414)
	plot(t, pd)
	title 'dpot check'
	hold on
	d = diffc(pp) / (t(2)-t(1));
	plot(t(1:end-1), d(1:end-1,:), '--')
	hold off
	axisy([-1 1] * 1.1 * delta)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function y = wprctile(x, p, w)
% WPRCTILE Percentiles of a weighted sample.
%
%    Description
%      Y = PRCTILE(X, P, W) returns percentiles of the values in X. P
%       is a scalar or a vector of percent values. W is a vector of
%       unnormalized weights for samples. Length of W has to be same
%       as length of X. X need to be a a vector. Y is the same size as
%       P, and Y(i) contains the P(i)-th percentile.
%
%      Example
%       y = prctile(x,50,w); % the median of x given sample weights w
%
%    See also wmean, prctile

% BUGS: Accepts only vector valued X

% Copyright (c) 2000-2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

x=sort(x);
p=p./100;
y=p;
ww=cumsum(w);ww=ww./ww(end);
for j=1:length(p)
  wi=min(find(ww>=p(j)));
  if wi==1
    y(j)=x(1);
  else
    w1=ww(wi-1);x1=x(wi-1);
    y(j)=x1+(x(wi)-x1).*(p(j)-w1)./(ww(wi)-w1);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% GEOMETRIC COLLISION AVOIDANCE AGENT (agent_2D_vectorSharing.m) %%%%%%%%%
% This programs contains a generic agent object with the 2017 3D geometric
% avoidance alogorithm applied to conduct course corrections.

% Author: James A. Douthwaite

classdef agent_2D_vectorSharing < agent_2D & agent_vectorSharing
    %% ////////////////////// MAIN CLASS METHODS //////////////////////////
    methods
        % Constructor
        function [this] = agent_2D_vectorSharing(varargin)

            % Call the super class
            this@agent_2D(varargin);

            % //////////////////// SENSOR PARAMETERS //////////////////////
            [this.SENSORS] = this.GetDefaultSensorParameters();     % Default sensing
%             [this.SENSORS] = this.GetCustomSensorParameters();       % Experimental sensing
            % /////////////////////////////////////////////////////////////

            % //////////////// Check for user overrides ///////////////////
            this = this.ApplyUserOverrides(varargin); % Recursive overrides
            % /////////////////////////////////////////////////////////////
        end
        % Setup - X = [x y psi dx dy dpsi]
        function [this] = setup(this,localXYZVelocity,localXYZrotations)
            % This function calculates the intial state for a generic
            % object.
            % The default state vector:
            % [x y psi dx dy dpsi]
            % INITIALISE THE 2D STATE VECTOR WITH CONCANTINATED VELOCITIES
            [this] = this.setup_2DVelocities(localXYZVelocity,localXYZrotations);
        end
        % Main
        function [this] = main(this,ENV,varargin)
            % This function is designed to house a generic agent process
            % cycle that results in an acceleration vector in the global axis.
            % INPUTS:
            % varargin - Cell array of inputs
            % >dt      - The timestep
            % >objects - The detectable objects cell array of structures
            % OUTPUTS:
            % obj      - The updated project

            % INPUT HANDLING
            dt = ENV.dt;

            % PLOT AGENT FIGURE
            visualiseProblem = 0;
            visualiseAgent = 1;
            if this.objectID == visualiseAgent && visualiseProblem == 1
                overHandle = figure(100+this.objectID);
                hold on; grid on;
                axis equal;
                xlabel('x_{m}'); ylabel('y_{m}'); zlabel('z_{m}');
            end

            % //////////// CHECK FOR NEW INFORMATION UPDATE ///////////////
            [this,obstacleSet,agentSet] = this.GetAgentUpdate(ENV,varargin{1});       % IDEAL INFORMATION UPDATE
            % /////////////////////////////////////////////////////////////

            % /////////////////// WAYPOINT TRACKING ///////////////////////
            % Design the current desired trajectory from the waypoint.
            desiredHeadingVector = this.GetTargetHeading();
            desiredVelocity = desiredHeadingVector*this.v_nominal;

            % ////////////////// OBSTACLE AVOIDANCE ///////////////////////
            % Modify the desired velocity with the augmented avoidance velocity.
            avoidanceSet = [obstacleSet,agentSet];
            algorithm_start = tic; algorithm_indicator = 0;
            if ~isempty(avoidanceSet)
                algorithm_indicator = 1;
                % GET THE UPDATED DESIRED VELOCITY
                 [desiredHeadingVector,desiredSpeed] = this.GetAvoidanceCorrection(desiredVelocity,visualiseProblem);
                 desiredVelocity = desiredHeadingVector*desiredSpeed;
            end
            algorithm_dt = toc(algorithm_start);                           % Stop timing the algorithm

            desiredVelocity

            % ///////////////////// CONTROLLER ////////////////////////////
            this = this.Controller(dt,desiredVelocity);

            % ////////////// RECORD THE AGENT-SIDE DATA ///////////////////
            this = this.writeAgentData(ENV,algorithm_indicator,algorithm_dt);      % Record when the algorithm is ran
            this.DATA.inputNames = {'$v_x$ (m/s)','$v_y$ (m/s)','$\dot{\psi}$ (rad/s)'};
            this.DATA.inputs(1:length(this.DATA.inputNames),ENV.currentStep) = this.localState(4:6);         % Record the control inputs

            % // DISPLAY CONFLICT RESOLUTION
            if this.objectID == visualiseAgent && visualiseProblem == 1
                this = this.GetAnimationFrame(overHandle,ENV,'resolutionZone.gif');
                close(overHandle);
            end
        end
    end
    %% /////////////////////// AUXILLARY METHODS //////////////////////////
    methods (Access = public)
        % GET THE AVOIDANCE CORRECTION
        function [heading,speed] = GetAvoidanceCorrection(this,desiredVelocity,visualiseProblem)
            % This function calculates the collision avoidance velocity in
            % light of the current obstacles

            % Check we aren't stopping
            [heading,speed] = this.nullVelocityCheck(desiredVelocity);
            if speed == 0
                return
            end

            % AGENT KNOWLEDGE
            [p_i,v_i,r_i] = this.GetAgentMeasurements(); % Its own position, velocity and radius

            % Define the obstacle list
            obstacleIDs = [this.MEMORY([this.MEMORY.type] ~= OMAS_objectType.waypoint).objectID];

            % MOVE THROUGH THE PRIORITISED OBSTACLE SET
            optimalSet = [];
            for j = 1:numel(obstacleIDs)
                % Fetch the obstacle trajectories
                p_j = this.GetLastMeasurementByID(obstacleIDs(j),'position');
                % Neighbour conditions
                neighbourConditionA = j < this.maxNeighbours;            % Maximum number of neighbours
                neighbourConditionB = norm(p_j) < this.neighbourDist;       % [CONFIRMED]
                if ~neighbourConditionA || ~neighbourConditionB
                    continue
                end

                % Get further obstacle information
                v_j = this.GetLastMeasurementByID(obstacleIDs(j),'velocity');
                r_j = this.GetLastMeasurementByID(obstacleIDs(j),'radius');

                % OBSTACLE KNOWLEDGE
                p_j = p_j + p_i;
                v_j = v_j + v_i;                                           % Convert relative parameters to absolute

                % COMPUTE THE VECTOR SHARING PROBLEM
                [Voptimal] = this.Define2DVectorSharingVelocity(...
                    desiredVelocity,...
                    p_i,v_i,r_i,...
                    p_j,v_j,r_j,...
                    visualiseProblem);
                optimalSet = horzcat(optimalSet,[Voptimal;abs(norm(desiredVelocity - Voptimal))]);
            end

            % INTERPRETING THE AVOIDANCE VELOCITY SET
            if ~isempty(optimalSet)
                inputDim = 2;
                % THE CLOSEST OBSTACLE
                %                 avoidanceVelocity = optimalSet(1:inputDim,1);
                % FIND THE MINIMUM MAGNITUDE DEVIATION FROM THE DESIRED
                [~,minIndex] = min(optimalSet((inputDim+1),:),[],2);       % Return the index of the smallest vector
                avoidanceVelocity = optimalSet(1:inputDim,minIndex);
            else
                % NOTHING TO AVOID, CHOOSE OPTIMAL VELOCITY
                avoidanceVelocity = desiredVelocity;
            end

            % CHECK VELOCITY IS PERMISSIBLE
            if any(isnan(avoidanceVelocity))
%                 heading = [1;0];
%                 speed = 0;
                speed   = norm(desiredVelocity);
                heading = desiredVelocity/speed;    % Retain forward direction
            else
                [heading,speed] = this.nullVelocityCheck(avoidanceVelocity);
            end
        end
        % DEFINE THE VECTOR SHARING AVOIDANCE PROBLEM
        function [U_a] = Define2DVectorSharingVelocity(obj,desiredVelocity,p_a,v_a,r_a,p_b,v_b,r_b,visualiseProblem)
            % This function calculates the avoidance vectors based on the
            % principle of vector sharing.
            % INPUTS:
            % desiredVelocity - The true optimal vector
            % p_a - The agent's absolute position
            % v_a - The agent's absolute velocity
            % r_a - The agent's radius
            % p_b - The obstacle's absolute position
            % v_b - The obstacle's absolute velocity
            % r_b - The obstacle's radius
            % visualiseProblem - Plot flag
            % OUTPUTS:
            % U_a - The optimal vector heading, scaled by the desired
            %       velocity.

            % Generate the 3D inputs
            p_a = [p_a;0]; v_a = [v_a;0];
            p_b = [p_b;0]; v_b = [v_b;0];
            % Pass to the 3D function
            [U_a] = obj.Define3DVectorSharingVelocity(...
                desiredVelocity,...
                p_a,v_a,r_a,...
                p_b,v_b,r_b,...
                visualiseProblem);
            % Reform the inputs for 2D application
            U_a = U_a(1:2,1);
        end
    end
    % The function uses the same sensor characteristics from the parent
    % class "agent_vectorSharing.m".
end
% AGENT STATE VECTOR [x;y;phi;xdot;ydot;phidot]

repo: OpenMAS


text:
function patches = img2patch_batch(img, patchSize, stride)

a = patchSize(1); b = patchSize(2); c = stride;

if stride > 1
	cut = mod(size(img), patchSize);
	img = img(1:end-cut(1), 1:end-cut(2));
end

nPatches = prod(size(img)/c);

% Extended image from the original image
img = [ img img(:, 1:b-c); img(1:a-c, :) img(1:a-c, 1:b-c) ];
imgSize = size(img);

% Linear indices for the 1st patch
stride = [1 imgSize(1)];
limit = (patchSize-1).*stride;
ind = 1;
for dim = 1:numel(imgSize)
	ind = bsxfun (@plus, ind(:), 0:stride(dim):limit(dim));
end

% Linear indices for all patches
slides = imgSize - patchSize;
limit = slides.*stride;
stride = stride*c;
for dim = 1:numel(imgSize)
	ind = bsxfun( @plus, ind(:), 0:stride(dim):limit(dim) );
end
patches = reshape(img(ind(:)), prod(patchSize), nPatches);

end % img2patch_batch


text:
classdef Projector < handle

    properties (Access = protected)
        mesh
        connec
    end

    methods (Static, Access = public)

        function obj = create(cParams)
            p = ProjectorFactory();
            obj = p.create(cParams);
        end

    end

    methods (Static, Access = protected)

        function ord = determineQuadratureOrder(fun)
            switch fun.fType
                case 'L2'
                    ord = 'QUADRATIC';
                case 'FE'
%                     ord = 'LINEAR';
                    ord = 'QUADRATIC'; % needed to project P1 to P1D
                case 'GAUSSPOINTS'
                    ord = fun.quadrature.order;
            end
        end

    end

    methods (Access = protected)

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

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [bnet, LL, engine] = learn_params_em(engine, evidence, max_iter, thresh)
% LEARN_PARAMS_EM Set the parameters of each adjustable node to their ML/MAP values using batch EM.
% [bnet, LLtrace, engine] = learn_params_em(engine, data, max_iter, thresh)
%
% data{i,l} is the value of node i in case l, or [] if hidden.
%   Suppose you have L training cases in an O*L array, D, where O is the num observed
%   scalar nodes, and N is the total num nodes.
%   Then you can create 'data' as follows, where onodes is the index of the observable nodes:
%      data = cell(N, L);
%      data(onodes,:) = num2cell(D);
%   Of course it is possible for different sets of nodes to be observed in each case.
%
% We return the modified bnet and engine.
% To see the learned parameters for node i, use the construct
%   s = struct(bnet.CPD{i}); % violate object privacy
% LLtrace is the learning curve: the vector of log-likelihood scores at each iteration.
%
% max_iter specifies the maximum number of iterations. Default: 10.
%
% thresh specifies the thresold for stopping EM. Default: 1e-3.
% We stop when |f(t) - f(t-1)| / avg < threshold,
% where avg = (|f(t)| + |f(t-1)|)/2 and f is log lik.

if nargin < 3, max_iter = 10; end
if nargin < 4, thresh = 1e-3; end

verbose = 1;

loglik = 0;
previous_loglik = -inf;
converged = 0;
num_iter = 1;
LL = [];

while ~converged & (num_iter <= max_iter)
  [engine, loglik] = EM_step(engine, evidence);
  if verbose, fprintf('EM iteration %d, ll = %8.4f\n', num_iter, loglik); end
  num_iter = num_iter + 1;
  converged = em_converged(loglik, previous_loglik, thresh);
  previous_loglik = loglik;
  LL = [LL loglik];
end
if verbose, fprintf('\n'); end

bnet = bnet_from_engine(engine);

%%%%%%%%%

function [engine, loglik] = EM_step(engine, cases)

bnet = bnet_from_engine(engine); % engine contains the old params that are used for the E step
CPDs = bnet.CPD; % these are the new params that get maximized
num_CPDs = length(CPDs);
adjustable = zeros(1,num_CPDs);
for e=1:num_CPDs
  adjustable(e) = adjustable_CPD(CPDs{e});
end
adj = find(adjustable);
n = length(bnet.dag);

for e=adj(:)'
  CPDs{e} = reset_ess(CPDs{e});
end

loglik = 0;
ncases = size(cases, 2);
for l=1:ncases
  evidence = cases(:,l);
  [engine, ll] = enter_evidence(engine, evidence);
  loglik = loglik + ll;
  hidden_bitv = zeros(1,n);
  hidden_bitv(isemptycell(evidence))=1;
  for i=1:n
    e = bnet.equiv_class(i);
    if adjustable(e)
      fmarg = marginal_family(engine, i);
      CPDs{e} = update_ess(CPDs{e}, fmarg, evidence, bnet.node_sizes, bnet.cnodes, hidden_bitv);
    end
  end
end

for e=adj(:)'
  CPDs{e} = maximize_params(CPDs{e});
end

engine = update_engine(engine, CPDs);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


function dataout=modeDR_AEC_BSS(nummics, numrefs, datain)
%
% Perform dr, aec, and bss.
% nummics:              no. of mic channels
% numrefs:              no. of reference channels
% datain:               input data
% dataout:              output data
%

%% perform stft
addpath('stft2');
% fft size
fftsize=512;
stftshift=fftsize/2;

M=nummics;
R=numrefs;
N=M;

Xtf=cell(M+R, 1);
for m=1:M+R
    Xtf{m}=stft(datain(:, m), stftshift, fftsize, false);
end
[K, T]=size(Xtf{1});

Ytf=cell(M, 1);
for m=1:M
    Ytf{m}=zeros(K, T);
end

%% params
% dr filter length
DR_FLEN=5;
% forgetting factor for dr and aec
DRAEC_FORGET=0.999;
% forgetting factor for bss
BF_FORGET=0.999;
% the shape parameter of the source prior
GAMMA=0.2;
%
% used to keep stable
%
VAR_BIAS=0.01;
STABLE_EPS=1e-3;
DRAEC_DIAGLOAD=1e-6;
BF_DIAGLOAD=1e-6;

%% space for dr
% current mic data backup
Miccurrent=zeros(K, M);
% delayed mic data
drfsize=M*DR_FLEN;
Micdelay=zeros(K, drfsize);

% mic and delay correlation
Cmd=cell(K, 1);
for k=1:K
    Cmd{k}=zeros(M, drfsize);
end

% delay auto correlation
Cdd=cell(K, 1);
for k=1:K
    Cdd{k}=zeros(drfsize, drfsize);
end

% reverb path
Reverbpath=cell(K, 1);
for k=1:K
    Reverbpath{k}=zeros(M, drfsize);
end

%% space for aec
% mic-reference correlation
Cmr=cell(K, 1);
for k=1:K
    Cmr{k}=zeros(M, R);
end

% reference auto correlation
Crr=cell(K, 1);
for k=1:K
    Crr{k}=zeros(R, R);
end

% echo path
Echopath=cell(K, 1);
for k=1:K
    Echopath{k}=zeros(M, R);
end

%% space for bss
% the weighted correlation matrices
C1=cell(K, 1);
C2=cell(K, 1);
for k=1:K
    C1{k}=STABLE_EPS*eye(M, M);
    C2{k}=STABLE_EPS*eye(M, M);
end

% demixing matrices
Demix=cell(K, 1);
for k=1:K
    Demix{k}=eye(N, M);
end

%% perform iteration
for tau=1:T
    %% perform dr
    % direct and early reverberation
    Early=zeros(K, M);

    %
    % shift in new data
    %
    % shift old mic data back
    Micdelay=circshift(Micdelay, M, 2);
    % shift in delayed mic data
    Micdelay(:, 1:M)=Miccurrent;

    % mic data backup
    for m=1:M
        Miccurrent(:, m)=Xtf{m}(:, tau);
    end

    for k=1:K
        % calculate late reverberation
        ref=Micdelay(k, :).';
        late=Reverbpath{k}*ref;

        % direct and early reverberation
        mic=Miccurrent(k, :).';
        early=mic-late;
        % output data
        Early(k, :)=early.';

        %
        % calculate nonlinearity
        %
        xsq=abs(mic).^2;
        ysq=abs(early).^2;

        phi=0;
        for m=1:M
            if ysq(m)<=xsq(m)
                phi=phi+ysq(m);
            else
                phi=phi+xsq(m);
            end
        end

        phi=(1-DRAEC_FORGET)*(phi+VAR_BIAS)^((GAMMA-2)/2);

        % update mic ref correlation
        Cmd{k}=DRAEC_FORGET*Cmd{k}+phi*(mic*ref');

        % update ref auto-correlation
        Cdd{k}=DRAEC_FORGET*Cdd{k}+phi*(ref*ref');

        % update reverb path
        Reverbpath{k}=Cmd{k}/(Cdd{k}+DRAEC_DIAGLOAD*eye(drfsize, drfsize));
    end

    %% perform aec
    % nearend data
    Nearend=zeros(K, M);

    for k=1:K
        % reference data
        ref=zeros(R, 1);
        for r=1:R
            ref(r)=Xtf{M+r}(k, tau);
        end

        % echo
        echo=Echopath{k}*ref;

        % nearend
        mic=Early(k, :).';
        nearend=mic-echo;
        % output data
        Nearend(k, :)=nearend.';

        %
        % calculate nonlinearity
        %
        xsq=abs(mic).^2;
        ysq=abs(nearend).^2;

        phi=0;
        for m=1:M
            if ysq(m)<=xsq(m)
                phi=phi+ysq(m);
            else
                phi=phi+xsq(m);
            end
        end

        phi=(1-DRAEC_FORGET)*(phi+VAR_BIAS)^((GAMMA-2)/2);

        % update mic ref correlation
        Cmr{k}=DRAEC_FORGET*Cmr{k}+phi*(mic*ref');

        % update ref auto-correlation
        Crr{k}=DRAEC_FORGET*Crr{k}+phi*(ref*ref');

        % update echo path
        Echopath{k}=Cmr{k}/(Crr{k}+DRAEC_DIAGLOAD*eye(R, R));
    end

    %% perform bss
    Bssout=zeros(K, M);

    %
    % calculate nonlinearity
    %
    phi1=0;
    phi2=0;

    for k=1:K
        x=Nearend(k, :).';
        y=Demix{k}*x;
        % output data
        Bssout(k, :)=y.';

        phi1=phi1+abs(y(1))^2;
        phi2=phi2+abs(y(2))^2;
    end

    phi1=(1-BF_FORGET)*(phi1+VAR_BIAS)^((GAMMA-2)/2);
    phi2=(1-BF_FORGET)*(phi2+VAR_BIAS)^((GAMMA-2)/2);

    % update the demixing matrices
    for k=1:K
        %
        % accumulate the weighted correlation
        %
        x=Nearend(k, :).';
        C1{k}=BF_FORGET*C1{k}+phi1*(x*x');
        C2{k}=BF_FORGET*C2{k}+phi2*(x*x');

        % solve gev problem
        D=heig2(BF_DIAGLOAD, C2{k}, C1{k});
        Demix{k}=D;
    end

    for m=1:M
        Ytf{m}(:, tau)=Bssout(:, m);
    end
end

%% perform istft and output signal
dataout=zeros(dataLength(T, stftshift, fftsize ), N);
for n=1:N
    dataout(:, n)=istft(Ytf{n}, stftshift, false);
end

end


text:
function jaabaPerUserDirName=JaabaPerUserSettingsAbsDirName()

if isunix() || ismac()
  settingsDirName=getenv('HOME');
  jaabatxt = '.jaaba.txt';
  jaabadir = '.jaaba';
else
  settingsDirName = getenv('APPDATA');
  jaabatxt = 'jaaba.txt';
  jaabadir = 'jaaba';
end
if isdeployed,
%   if exist(fullfile(settingsDirName,jaabatxt),'file')
%     ff = fopen(fullfile(settingsDirName,jaabatxt));
%     jaabaPerUserDirName = fgetl(ff);
%     fclose(ff);
%   else
%     jaabaPerUserDirName = uigetdir(settingsDirName,'Select a file to store JAABA configuration files');
%     if ~exist(jaabaPerUserDirName)
%       mkdir(jaabaPerUserDirName);
%     end
%     ff  = fopen(fullfile(settingsDirName,jaabatxt),'w');
%     fprintf(ff,'%s',jaabaPerUserDirName);
%     fclose(ff);
%     if jaabaPerUserDirName == 0,
%       error('Cannot continue without a settings dir');
%     end
%     dd = dir(fullfile(pwd,'params','*.xml'));
%     for ndx = 1:numel(dd)
%       copyfile(fullfile(pwd,'params',dd(ndx).name),....
%         fullfile(jaabaPerUserDirName,dd(ndx).name));
%     end
%   end
  jaabaPerUserDirName = 'params';
else
  jaabaPerUserDirName=fullfile(settingsDirName,jaabadir);
end

if ~exist(jaabaPerUserDirName,'dir');
  mkdir(jaabaPerUserDirName);
end

end

author: kristinbranson

repo: JAABA


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Author: Nicholas A. Battista
% Institution: The College of New Jersey (TCNJ)
% Created: March 20, 2018
% Date of Last Revision: March 22, 2018
%
% Written for use in MAT 331: Numerical Analysis
%
% PURPOSE:  -To compute the error ssociated with applying Trapezoid Rule to
%            two different integral cases:
%                (1) a periodic integrand on a periodic integration domain
%                (2) non-periodic integrand
%
%           -To demonstrate exponential convergence of Trap. Rule in Case (1)
%            above
%
% Inputs:
%
% Returns:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Trap_Rule_Error()

% Number of subintervals to test
NVec = [1e0:1e0:9e0 1e1:1e1:9e1 1e2:1e2:9e2 1e3:1e3:9e3 1e4:1e4:9e4 1e5:1e5:9e5 1e6:1e6:9e6];

%
% Storage Initialization
ErrVecPeriodic = zeros(length(NVec));
ErrVecNonPeriodic = ErrVecPeriodic;

%
% Compute Integral Approximation for Different #'s of Subintervals
%
for i=1:length(NVec)

   flag = 1; % For Periodic
   ErrVecPeriodic(i) = Trap_Rule( NVec(i), 0, 1, flag);

   flag = 0; % For Non-Periodic
   ErrVecNonPeriodic(i) = Trap_Rule( NVec(i), 0, 1, flag);

end

%
% Print Convergence Rate Information to Screen
slopeNonPeriodic = give_Me_Slope( ErrVecNonPeriodic,40,50 );
fprintf('\nConvergence Rate for Periodic Functions Appears to be Exponential, e.g.,\n');
fprintf('     we see a linear relationship between log(error) and # of subintervals, N\n');
fprintf('\nConvergence Rate for Non-Periodic Functions Appears to be: %d\n\n',slopeNonPeriodic);

%
% PLOT THE ABSOLUTEL ERROR -> GIVE CONVERGENCE PLOTS!
%
fs = 18;
ms = 30;
figure(1) % PERIODIC CASE
semilogy(NVec(1:45),ErrVecPeriodic(1:45),'b.','MarkerSize',ms); hold on;
%loglog(NVec(1:45),ErrVecPeriodic(1:45),'b.','MarkerSize',ms); hold on;
title('Periodic Case');
ylabel('Absolute Error');
xlabel('Number of Sub-Intervals');
set(gca,'FontSize',fs);

figure(2) % NON-PERIODIC CASE
loglog(NVec,ErrVecNonPeriodic,'b.','MarkerSize',ms); hold on;
title('Non-Periodic Case');
ylabel('Absolute Error');
xlabel('Number of Sub-Intervals');
set(gca,'FontSize',fs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes the absolute error using the trapezoid rule for a
%           particular number of subintervals, N.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function err = Trap_Rule(N,a,b,flagCase)

% Uniformly spaced integration nodes
x = a:(b-a)/N:b;

% Compute trapezoid rule for periodic function
int = 0;
for i=1:N
   int = int + (b-a)/(2*N)*( f(x(i),flagCase) + f(x(i+1),flagCase) );
end

%
% Return Error for Particular Case
%
if flagCase == 1

    % PERIODIC CASE
    intExactPeriodic = 0.132214293037990;
    err = abs( int - intExactPeriodic );

else

    % NON-PERIODIC CASE
    intExactNonPeriodic = 0.455122322888408;
    err = abs( int - intExactNonPeriodic );

end

% PLOT!
%x=a:0.001:b;
%for i=1:length(x)
%    plot( x(i) , f(x(i)), '*'); hold on;
%end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: returns integrand function value
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function val = f(x,flagCase)

%
% Return Function Value for Particular Case
%
if flagCase == 1

    % PERIODIC: f(x) = (cos(2*pi*x))^2 / ( 1 + e^(sin(2*pi*x)) )^2;
    val = ( cos(2*pi*x) )^2 / ( 1 + exp( sin(2*pi*x) ) )^2;

else

    % NON-PERIODIC: f(x) = (x^2+3)(cos(2*pi*x))^2 / ( 1 + e^(sin(2*pi*x)) )^2;
    val = (x^2+3)*(cos(2*pi*x))^2 / ( 1 + exp(sin(2*pi*x)) )^2;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: give me slope of convergence plot
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function m = give_Me_Slope( errVec, val1, val2 )

m = ( log( errVec(val1) ) - log( errVec(val2) ) ) / (log(val1) - log(val2));

author: nickabattista


text:
function segment = p08_boundary_segment ( segment_index, m, segment_length )

%*****************************************************************************80
%
%% P08_BOUNDARY_SEGMENT returns a boundary segment in problem 08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 January 2005
%
%  Author:
%
%    John Burkardt
%
%  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, integer SEGMENT_INDEX, the index of the boundary segment.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer SEGMENT_LENGTH, the number of points in the segment.
%
%    Output, real SEGMENT(M,SEGMENT_LENGTH), the
%    points that make up the boundary segment.
%
  c1 = [ 0.0, 0.0 ];
  c2 = [ 0.6, 0.0 ];
  r1 = 1.0;
  r2 = 0.1;
  theta1 = pi / 12.0;
%
%  Segment 1: the outer boundary.
%
  if ( segment_index == 1 )

    a = ( sqrt ( 119.0 ) - 9.0 ) / 20.0;
    theta2 = atan2 ( a, a + 0.9 );
%
%  Work out the appropriate segment lengths, and then
%  adjust N6, if necessary, to account for roundoff.
%
    n1 = round ( ( segment_length - 1 ) / ( 2.0 + 2.0 * ( theta1 - theta2 ) ...
              + 2.0 * a * sqrt ( 2.0 ) ) );
    n2 = round ( ( theta1 - theta2 ) * n1 );
    n3 = round (    a * sqrt ( 2.0 ) * n1 );

    n2 = max ( n2, 1 );
    n3 = max ( n3, 1 );

    n4 = n3;
    n5 = n2;
    n6 = round ( ( segment_length - 1 - n2 - n3 - n4 - n5 ) / 2 );
    n1 = segment_length - 1 - n2 - n3 - n4 - n5 - n6;

    j = 0;

    s(1:2) = [ 0.0, 0.0 ];
    idiot1 =   cos ( theta1 );
    idiot2 = - sin ( theta1 );
    t(1:2) = [ idiot1, idiot2 ];

    for i = 1 : n1
      j = j + 1;
      segment(1:2,j) = ( ( n1 - i + 1 ) * s(1:2)   ...
                       + (      i - 1 ) * t(1:2) ) ...
                       / ( n1         );
    end

    for i = 1 : n2

      theta = ( - ( n2 - i + 1 ) * theta1   ...
                - (      i - 1 ) * theta2 ) ...
              /   ( n2         );

      j = j + 1;
      idiot1 = cos ( theta );
      idiot2 = sin ( theta );
      segment(1:2,j) = [ idiot1, idiot2 ];
    end

    idiot1 =  cos ( theta2 );
    idiot2 = -sin ( theta2 );
    s(1:2) = [ idiot1, idiot2 ];
    t(1:2) = [ 0.9, 0.0 ];

    for i = 1 : n3
      j = j + 1;
      segment(1:2,j) = ( ( n3 - i + 1 ) * s(1:2)   ...
                       + (      i - 1 ) * t(1:2) ) ...
                       / ( n3         );
    end

    s(1:2) = [ 0.9, 0.0 ];
    idiot1 = cos ( theta2 );
    idiot2 = sin ( theta2 );
    t(1:2) = [ idiot1, idiot2 ];

    for i = 1 : n4
      j = j + 1;
      segment(1:2,j) = ( ( n4 - i + 1 ) * s(1:2)   ...
                       + (      i - 1 ) * t(1:2) ) ...
                       / ( n4         );
    end

    for i = 1 : n5

      theta = ( ( n5 - i + 1 ) * theta2   ...
              + (      i - 1 ) * theta1 ) ...
              / ( n5         );

      j = j + 1;
      idiot1 = cos ( theta );
      idiot2 = sin ( theta );
      segment(1:2,j) = [ idiot1, idiot2 ];
    end

    idiot1 = cos ( theta1 );
    idiot2 = sin ( theta1 );
    s(1:2) = [ idiot1, idiot2 ];
    t(1:2) = [ 0.0, 0.0 ];

    for i = 1 : n6
      j = j + 1;
      segment(1:2,j) = ( ( n6 - i + 1 ) * s(1:2)   ...
                       + (      i - 1 ) * t(1:2) ) ...
                       / ( n6         );
    end

    j = j + 1;
    segment(1:2,j) = [ 0.0, 0.0 ];
%
%  Segment 2: the circular hole.
%
  elseif ( segment_index == 2 )

    for j = 1 : segment_length
      theta = ( segment_length - j ) * 2.0 * pi / ( segment_length - 1 );
      segment(1,j) = c2(1) + r2 * cos ( theta );
      segment(2,j) = c2(2) + r2 * sin ( theta );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P08_BOUNDARY_SEGMENT - Fatal error!\n' );
    fprintf ( 1, '  Illegal SEGMENT_INDEX = %d\n', segment_index );
    error ( 'P08_BOUNDARY_SEGMENT - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%ASSEMMAT Assemble monolithic matrix.
%
%   [ A, T_A, T_SP ] = ASSEMMAT( PROB, S_A, I_CUB, N_CMAX, F_SPARSE, I_HRZ, SOLCOMP )
%
%   Assemble monolithic (coupled) matrix (for all dependent variables)
%   for the field S_A in the finite element data struct PROB. S_A can
%   be either "m" or "a", designating assembly of the mass matrix
%   defined in PROB.EQN.M, or system matrix defined in PROB.EQN.A.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       prob        struct                 FEA problem data struct
%       s_a         char    {m/a}          Matrix field to assemble (mass/system)
%       i_cub       scalar  {2}            Numerical integration rule
%       n_cmax      scalar  {50000}        Max number of cells to assemble
%                                          for at once (to limit memory consumption)
%       f_sparse    logical {true}         Return sparse/struct matrix format
%       i_hrz       logical {false}        Apply HRZ diagonal (mass) lumping
%       solcomp     {all dvars/subd}       Dependent variables/subdomains to assemble for
%                                                                                .
%       Output      Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       A           sparse/struct [n_A]    Assembled system matrix
%       t_a         scalar                 Time spent assembling matrix
%       t_sp        scalar                 Time for sparse matrix conversion
%
%   See also ASSEMBLEA, ASSEMBLEPROB

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function ind=qilv(I,I2,Ws)

% QILV()
%==================================================================
% Quality Index based on Local Variance
% Santiago Aja-Fernandez
%
% santi@bwh.harhard.edu
% Accorging to
%
% S. Aja-Fernández, R. San José Estépar, C. Alberola-López and C.F. Westin,
% "Image quality assessment based on local variance", EMBC 2006,
% New York, Sept. 2006.
%
%------------------------------------------------------------------
%
% The function calculates a global compatibility measure
% between two images, based on their local variance distribution.
%
%------------------------------------------------------------------
%
% INPUT:   (1) I: The first image being compared
%          (2) I2: the second image being compared
%          (3) Ws: window for the estimation of statistics:
%		If Ws=0: default gaussian window
%               If Ws=[M N] MxN square window
%
%
% OUTPUT:
%          (1) ind: Quality index (between 0 and 1)
%
% Default usage:
%
%   ind=s_correct(I,I2,0);
%
%
%==================================================================

%the following variables can be added to avoid NaN, but
%usually it is not necessary

L=255;
K=[0.01 0.03];
C1 = (K(1)*L)^2;
C2 = (K(2)*L)^2;
%C1=0;
%C2=0;

%Window
if Ws==0
	window = fspecial('gaussian', 11, 1.5);
else
	window=ones(Ws);
end
window = window/sum(window(:));

%Local means
M1=filter2(window, I, 'valid') ;
M2=filter2(window, I2, 'valid') ;
%Local Variances
V1 = filter2(window, I.^ 2, 'valid') - M1.^ 2;
V2 = filter2(window, I2.^ 2, 'valid') - M2.^ 2;

%Global statistics:

m1=mean(V1(:));
m2=mean(V2(:));
s1=std(V1(:));
s2=std(V2(:));
s12=mean2((V1-m1).*(V2-m2));

%Index
ind1=((2*m1*m2+C1)./(m1.^2+m2.^ 2+C1));
ind2=(2*s1*s2+C2)./(s1.^ 2+s2.^ 2+C2);
ind3=(s12+C2/2)./(s1*s2+C2/2);
ind=ind1.*ind2.*ind3;

author: Sable

repo: mcbench-benchmarks


text:
function [GAH] = dicomrt_gahcal(VOI,gamma,gamma_xmesh,gamma_ymesh,gamma_zmesh,gahselect,xgrid)
% dicomrt_gahcal(VOI,gamma,gamma_xmesh,gamma_ymesh,gamma_zmesh,gahselect,xgrid)
%
% Calculate GAMMA-AREA-HISTOGRAMS for VOIs and cell_case_study.
% This function calculate GAHs for all the available VOIs or for a selected number of them.
%
% gamma contains the 3D gamma map
% gamma_xmesh,gamma_ymesh,gamma_zmesh are x-y-z coordinates of the center of the gamma-voxels
% VOIs is a cell array and contains all the Volumes of Interest
% gahselect is an OPTIONAL parameter which contains the number of the VOI to calculate GAHs for.
% xgrid is an OPTIONAL parameter that sets the resolution for GAH calculation (default = 30)
%
% GAHs are stored in a cell array with the following structure:
%
%  -------------------------------
%  | [GAH-name] | [3D gamma mask] |
%  |            -------------------
%  |            | [gah section 1] |
%  |            | [gah section 2] |
%  |            |     ...         |
%  |            | [gah section m] |
%  |            -------------------
%  |            | [info area 1]   |
%  |            | [info area 2]   |
%  |            |     ...         |
%  |            | [info area m]   |
%  |            -------------------
%  |            | Pixel area      |
%  --------------------------------
%
% [gamma data] is a 2 columns vector with following structure:
%
% -----------------
% | gamma | count |
% -----------------
% |       |       |
% |       |       |
% |       |       |
% |       |       |
% -----------------
%
% [info area] is a 2 columns vector with following structure:
%
% ---------------------------------
% | area sec | z location section |
% ---------------------------------
%
% Example:
%
% gah=dicomrt_gahcal(VOI,gamma,gamma_xmesh,gamma_ymesh,gamma_zmesh,4) returns
% the GAH for the VOI number 4 in VOI and the gamma distribution stored in the
% 3D matrix 'gamma'.
%
% See also dicomrt_gvhcal, dicomrt_gvhplot
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument and set-up some parameters and variables
error(nargchk(5,7,nargin))

% Check input parameters and set default parameter
list_variables=who;
match=0;

% Check voi
[VOI_temp]=dicomrt_checkinput(VOI);
VOI=dicomrt_varfilter(VOI_temp);

if nargin==5 & (exist('gahselect')==1 | exist('xgrid')==1)
    error('dicomrt_gahcal: Not enough input paramenters. Exit now!');
end

if exist('gahselect')~=1
    gahselect=[1:1:size(VOI,1)];
end

% Check consistency for z values
if gamma_zmesh(1)~=VOI{gahselect,2}{1}(1,3)
    error(['dicomrt_gahcal: GAMMA map was not calculated for the selected voi (',VOI{gahselect,1},'). Exit now!']);
end

% Define xgrid
if exist('xgrid')~=1
    xgrid=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.2 1.4 1.6 1.8 2.0 ...
            3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 Inf];
end

% Define cell array
GAH=cell(1,2);
GAH{1,1}=['GAH-',VOI{gahselect,1}];

% Notes on pixel area calculation
% Although slice thickness may vary the pixel area is supposed to be the
% same through all the volume.
pixelarea=(abs(gamma_xmesh(2)-gamma_xmesh(1)).*(abs(gamma_ymesh(2)-gamma_ymesh(1))));

% start calculating time
disp(['GAH calculation for VOI: ',VOI{gahselect,1}]);
% Initialise parameters
count=0; % count dimensions of gamma_section
count_empty=0; % count number of times gamma_section is not defined
% mask gamma matrix with current VOI
[mask_VOI,volume_VOI,mask4VOI]=dicomrt_mask(VOI_temp,gamma,gamma_xmesh,gamma_ymesh,gamma_zmesh,gahselect,'nan','n');
% slice by slice select contents of mask_VOI ~=nan and put it into a vector for GAH calculation

for kk=1:size(mask_VOI{2,1},3)
    xgrid_temp=xgrid;
    agrid=histc(reshape(mask_VOI{2,1}(:,:,kk),1,size(mask_VOI{2,1},1)*size(mask_VOI{2,1},2)),xgrid_temp);
    % get back voxels with dose<min(xgrid) cause histc does not
    % account for it
    temp_min=find(mask_VOI{2,1}(:,:,kk)<xgrid_temp(1) & mask_VOI{2,1}(:,:,kk)>0);
    if isempty(temp_min)~=1
        agrid(1)=agrid(1)+length(temp_min);
    end
    % Inf value will won't be plotted in the x grid. Therefore we
    % collect all the Inf values in the last slot of xgrid/vgrid
    agrid(end-1)=agrid(end-1)+agrid(end);
    agrid(end)=[];
    xgrid_temp(end)=[];
    temp2=cumsum(agrid);
    agrid=agrid.*pixelarea;
    info_area=[temp2(end)*pixelarea VOI{gahselect,2}{kk}(1,3)];
    gahdata=[xgrid_temp' agrid'];
    % store data in cell array
    GAH{1,2}{2,1}{kk}=gahdata;
    GAH{1,2}{3,1}{kk}=info_area;
    % plot frequency GAH
    if isempty(agrid)~=1
        figure;
        if isempty(inputname(2))==1
            set(gcf,'Name',['dicomrt_gahcal: ','fGAH for ',inputname(1)]);
        else
            set(gcf,'Name',['dicomrt_gahcal: ','fGAH for ',inputname(2)]);
        end
        agrid(end-1)=agrid(end-1)+agrid(end);
        agrid(end)=[];
        xgrid_temp(end)=[];
        bar(xgrid_temp,agrid);
        title([VOI{gahselect,1},' section: ',num2str(kk)],'Fontsize', 18,'Interpreter','none');
        xlabel('\gamma','Fontsize', 14);
        ylabel('Area [cm^2]','Fontsize', 14);
    else
        warning(['dicomrt_gahcal: GAH is null for VOI: ',GAH{k,1},' section: ',num2str(kk)]);
        gahdata=[0,0];
        GAH{1,2}{2,1}{kk}=gahdata;
    end
    % store other data in cell array
    GAH{1,2}{1,1}=mask_VOI{2,1};
    GAH{1,2}{4,1}=pixelarea;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [m3] = cc2m3(cc)
% Convert volume from cubic centimeters to cubic meters.
% Chad Greene 2012
m3 = cc* 0.000001;


text:
function deri = derivative_sigm_j(a)
    % derivative of tanh(x) is 1-tanh(x).^2
    %tmp = gones(size(a));
    %tmp = gpuArray(ones(size(a)));
    %deri = (tmp - (a / 1.7159).^2) * 1.7159 * 0.6667;
    %deri = tmp - (a.^2);
    %deri = gsingle(not(not(max(0, a))));
    deri = gpuArray(not(not(max(0, a))));
end

author: jimmy-ren


text:
function [KB] = EB2KB(EB)
% Convert computery things from exabytes to kilobytes.
% Chad A. Greene 2012
KB = EB*1125899906843000;


text:
function [superpixels, ucm, nSP, spArea] = getSuperpixels(imName, ucmThresh)
% function [superpixels, ucm, nSP, spArea] = getSuperpixels(imName, ucmThresh)
%  includes the boundary at value ucmThresh, ucmThresh must be in double, like 34/255.
%  Returns the UCM in uint8 format.

	paths = getPaths();
	ucm = getUCM(imName);
	ucmD = im2double(ucm);
	superpixels = bwlabel(ucmD < ucmThresh);
	superpixels = superpixels(2:2:end,2:2:end);
	nSP = max(superpixels(:));
	if(nargout > 2)
		spArea = histc(superpixels(:),1:nSP);
	end
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function [w, W_f, W_v] = fromFrameVec(F,v)

% FROMFRAMEVEC  From frame function for vectors.
%   FROMFRAMEVEC(F,V) transforms the 3d vector V from frame F to the global
%   frame.
%
%   [w, W_f, W_v] = FROMFRAMEVEC(...) returns the Jacobians wrt F and V.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

w = F.R*v;

if nargout > 1   % Jacobians.

    if size(v,2) == 1

        s = 2*q2Pi(F.q)*v;

        W_q = [...
            s(2) -s(1)  s(4) -s(3)
            s(3) -s(4) -s(1)  s(2)
            s(4)  s(3) -s(2) -s(1)];
        W_v = F.R;
        W_f = [zeros(3) W_q];

    else % multiple points
        error('??? Can''t give Jacobians for multiple points');
    end
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function o = cat(dr,varargin)
% Concatenate file_array objects.  The result is a non-simple object
% that can no longer be reshaped.
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: cat.m 7147 2017-08-03 14:07:01Z spm $

if dr>32 || dr<0, error('Unknown command option.'); end
dr   = max(round(dr),1);
d    = ones(nargin-1,16);
tmp  = {};
dpos = 0;
for i=1:nargin-1
    vi = varargin{i};
    if strcmp(class(vi),'file_array')
        sz                = size(vi);
        d(i,1:length(sz)) = sz;
        svi               = struct(vi);
        svi               = svi(:);
        for j=1:length(svi(:))
            if length(svi(j).pos)<dr
                svi(j).pos((length(svi(j).pos)+1):dr) = 1;
            end
            svi(j).pos(dr)= svi(j).pos(dr) + dpos;
        end
        dpos              = dpos + d(i,dr);
        tmp{i}            = svi;
    else
        error(['Conversion to file_array from ' class(vi) ' is not possible.']);
    end
end
if any(diff(d(:,[1:(dr-1) (dr+1):end]),1,1))
    error('All matrices on a row in the bracketed expression must have the same number of rows.');
else
    o = vertcat(tmp{:});
    o = file_array(o);
end

author: spm


text:
if ~exist('calib_name')|~exist('format_image'),
   data_calib_no_read;
   return;
end;

check_directory;

if ~exist('n_ima'),
   data_calib_no_read;
   return;
end;

check_active_images;

images_read = active_images;

if exist('image_numbers'),
   first_num = image_numbers(1);
end;

% Just to fix a minor bug:
if ~exist('first_num'),
   first_num = image_numbers(1);
end;

image_numbers = first_num:n_ima-1+first_num;

no_image_file = 0;

% Step used to clean the memory if a previous atttempt has been made to read the entire set of images into memory:
for kk = 1:n_ima,
    if (exist(['I_' num2str(kk)])==1),
        clear(['I_' num2str(kk)]);
    end;
end;

fprintf(1,'\nChecking directory content for the calibration images (no global image loading in memory efficient mode)\n');

one_image_read = 0;

i = 1;

while (i <= n_ima), % & (~no_image_file),

    if active_images(i),

        %fprintf(1,'Loading image %d...\n',i);

        if ~type_numbering,
            number_ext =  num2str(image_numbers(i));
        else
            number_ext = sprintf(['%.' num2str(N_slots) 'd'],image_numbers(i));
        end;

        ima_name = [calib_name  number_ext '.' format_image];

        if i == ind_active(1),
            fprintf(1,'Found images: ');
        end;

        if exist(ima_name),

            fprintf(1,'%d...',i);

            if ~one_image_read

                if format_image(1) == 'p',
                    if format_image(2) == 'p',
                        I = double(loadppm(ima_name));
                    else
                        I = double(loadpgm(ima_name));
                    end;
                else
                    if format_image(1) == 'r',
                        I = readras(ima_name);
                    else
                        I = double(imread(ima_name));
                    end;
                end;

                if size(I,3)>1,
                    I = 0.299 * I(:,:,1) + 0.5870 * I(:,:,2) + 0.114 * I(:,:,3);
                end;

                if size(I,1)~=480,
                    small_calib_image = 1;
                else
                    small_calib_image = 0;
                end;

                [Hcal,Wcal] = size(I); 	% size of the calibration image

                [ny,nx] = size(I);

                one_image_read = 1;

            end;

        else

            images_read(i) = 0;

        end;

    end;

    i = i+1;

end;

ind_read = find(images_read);

if ~(exist('map')==1), map = gray(256); end;

active_images = images_read;

fprintf(1,'\ndone\n');

repo: ekfmonoslam


function [ GAMMA ] = Mellin_NIG_European_Gamma( S_0, W, T, r, q, alpha, beta, delta, N1)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

f = W*exp(-r*T);

gam = sqrt(alpha^2 - beta^2);
k0 = log(S_0/W) + (r - q + delta*(sqrt(alpha^2 - (beta + 1)^2) - gam))*T;
adt = alpha*delta*T;
dta = 0.5*delta*T/alpha;

sum = 0;

if beta == 0
    % Symmetric Formula
    for n = 0:N1
        cons1 = k0^n / factorial(n);
        term = besselk(n/2 + 1, adt) / gamma((-n+1)/2) * (dta)^(-n/2);
        sum = sum + cons1*term;
    end
else
   % Asymmetric Formula

end

cons = f*alpha*exp(alpha*delta*T)/(S_0*S_0*sqrt(pi));
GAMMA = cons*sum;

end


text:
function [dataout] = ft_eventtiminganalysis(cfg, data)

% FT_EVENTTIMINGANALYSIS computes a model of single trial event-
% related activity, by estimating per trial the latency (and
% amplitude) of event-related signal components.
%
% Use as
%   [dataout] = ft_eventtiminganalysis(cfg, data)
% where data is single-channel raw data as obtained by FT_PREPROCESSING
% and cfg is a configuration structure according to
%
%  cfg.method  = method for estimating event-related activity
%                 'aseo', analysis of single-trial ERP and ongoing
%                         activity (according to Xu et al, 2009)
%                 'gbve', graph-based variability estimation
%                         (according to Gramfort et al, IEEE TBME 2009)
%  cfg.channel = Nx1 cell-array with selection of channels (default = 'all'),
%                see FT_CHANNELSELECTION for details
%  cfg.trials  = 'all' or a selection given as a 1xN vector (default = 'all')
%  cfg.output  = 'model', or 'residual', which returns the modelled data,
%                or the residuals.
%
% Method specific options are specified in the appropriate substructure.
%
% For the ASEO method, the following options can be specified:
%   cfg.aseo.noiseEstimate   = 'non-parametric' or 'parametric', estimate noise
%                              using parametric or non-parametric (default) method
%   cfg.aseo.tapsmofrq       = value, smoothing parameter of noise for
%                              nonparametric estimation (default = 5)
%   cfg.aseo.jitter          = value, time jitter in initial timewindow
%                              estimate (in seconds). default 0.050 seconds
%   cfg.aseo.numiteration    = value, number of iteration (default = 1)
%   cfg.aseo.initlatency     = Nx2 matrix, initial set of latencies in seconds of event-
%                              related components, give as [comp1start, comp1end;
%                              comp2start, comp2end] (default not
%                              specified). For multiple channels it should
%                              be a cell-array, one matrix per channel
%  Alternatively, rather than specifying a (set of latencies), one can also
%  specify:
%
%   cfg.aseo.initcomp        = vector, initial estimate of the waveform
%                              components. For multiple channels it should
%                              be a cell-array, one matrix per channel.
%
% For the GBVE method, the following options can be specified:
%   cfg.gbve.sigma             = vector, range of sigma values to explore in
%                                cross-validation loop (default: 0.01:0.01:0.2)
%   cfg.gbve.distance          = scalar, distance metric to use as
%                                evaluation criterion, see plugin code for
%                                more informatoin
%   cfg.gbve.alpha             = vector, range of alpha values to explor in
%                                cross-validation loop (default: [0 0.001 0.01 0.1])
%   cfg.gbve.exponent          = scalar, see plugin code for information
%   cfg.gbve.use_maximum       = boolean, (default: 1) consider the positive going peak
%   cfg.gbve.show_pca          = boolean, see plugin code (default 0)
%   cfg.gbve.show_trial_number = boolean, see plugin code (default 0)
%   cfg.gbve.verbose           = boolean (default: 1)
%   cfg.gbve.disp_log          = boolean, see plugin code (default 0)
%   cfg.gbve.latency           = vector [min max], latency range in s
%                                (default: [-inf inf])
%   cfg.gbve.xwin              = scalar smoothing parameter for moving
%                                average smoothing (default: 1), see
%                                eeglab's movav function for more
%                                information.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_SINGLETRIALANALYSIS_ASEO

% Copyright (C) 2018-2019, Jan-Mathijs Schoffelen DCCN
%
% 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$

ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    data
ft_preamble provenance data
ft_preamble trackconfig

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  % do not continue function execution in case the outputfile is present and the user indicated to keep it
  return
end

% ensure that the input data is valid for this function
data = ft_checkdata(data, 'datatype', {'raw+comp', 'raw'}, 'feedback', 'yes', 'hassampleinfo', 'yes');

% ensure that the required options are present
cfg         = ft_checkconfig(cfg, 'required', {'method'});
cfg.trials  = ft_getopt(cfg, 'trials',  'all', 1); % all trials as default
cfg.channel = ft_getopt(cfg, 'channel', 'all');
cfg.output  = ft_getopt(cfg, 'output',  'model');
% ensure that the options are valid
cfg = ft_checkopt(cfg, 'method', 'char', {'aseo' 'gbve'});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the actual computation is done in the middle part
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% select trials of interest
tmpcfg = keepfields(cfg, {'trials' 'channel' 'showcallinfo'});
data   = ft_selectdata(tmpcfg, data);
% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

% some error checks
if isfield(data, 'trial') && numel(data.trial)==0, ft_error('no trials were selected'); end
if numel(data.label)==0, ft_error('no channels were selected'); end

switch cfg.method
  case 'aseo'
    % define general variables that are used locally
    fsample = data.fsample; % Sampling Frequency in Hz
    nchan   = numel(data.label);
    nsample = numel(data.time{1}); %FIXME ASSUMING FIXED TIME AXIS ACROSS ALL TRIALS

    % setting a bunch of options, to be passed on to the lower level function
    if ~isfield(cfg, 'aseo'), cfg.aseo = []; end
    cfg.aseo.thresholdAmpH = ft_getopt(cfg.aseo, 'thresholdAmpH', 0.5);
    cfg.aseo.thresholdAmpL = ft_getopt(cfg.aseo, 'thresholdAmpL', 0.1);
    cfg.aseo.thresholdCorr = ft_getopt(cfg.aseo, 'thresholdCorr', 0.2);
    cfg.aseo.maxOrderAR    = ft_getopt(cfg.aseo, 'maxOrderAR',    5);
    cfg.aseo.noiseEstimate = ft_getopt(cfg.aseo, 'noiseEstimate', 'nonparametric');
    cfg.aseo.numiteration  = ft_getopt(cfg.aseo, 'numiteration',  1);
    cfg.aseo.tapsmofrq     = ft_getopt(cfg.aseo, 'tapsmofrq',     5);
    cfg.aseo.fsample       = fsample;
    cfg.aseo.nsample       = nsample;
    cfg.aseo.pad           = ft_getopt(cfg.aseo, 'pad', (2.*nsample)/fsample);

    % deal with the different ways with which the initial waveforms can be defined
    initlatency      = ft_getopt(cfg.aseo, 'initlatency', {});
    initcomp         = ft_getopt(cfg.aseo, 'initcomp',    {});
    jitter           = ft_getopt(cfg.aseo, 'jitter',      0.050); % half temporal width of shift in s

    if isempty(initlatency) && isempty(initcomp)
      ft_error('for the ASEO method you should supply either an initial estimate of the waveform component, or a set of latencies');
    elseif ~isempty(initlatency)
      % this takes precedence, and should contain per channel the begin and
      % end points of the subwindows in time, based on which the initial
      % subcomponents are estimated

      % ensure it to be a cell-array if the input is a matrix
      if ~iscell(initlatency)
        initlatency = repmat({initlatency},[1 nchan]);
      end
      make_init = true;
    elseif ~isempty(initcomp)
      % ensure it to be a cell-array if the input is a matrix
      if ~iscell(initcomp)
        initcomp = repmat({initcomp}, [1 nchan]);
      end
      make_init = false;
    end

    if make_init
      assert(numel(initlatency)==nchan);
      for k = 1:nchan
        % preprocessing data
        tmp     = cellrowselect(data.trial,k);
        chandat = cat(1,tmp{:});
        chandat = ft_preproc_baselinecorrect(chandat, nearest(data.time{1}, -inf), nearest(data.time{1}, 0));
        avgdat  = nanmean(chandat, 1);

        % set the initial ERP waveforms according to the preset parameters
        ncomp       = size(initlatency{k},1);
        initcomp{k} = zeros(nsample, ncomp);
        for m = 1:ncomp
          begsmp = nearest(data.time{1},initlatency{k}(m, 1));
          endsmp = nearest(data.time{1},initlatency{k}(m, 2));
          if begsmp<1,       begsmp = 1;       end
          if endsmp>nsample, endsmp = nsample; end

          tmp = avgdat(begsmp:endsmp)';
          initcomp{k}(begsmp:endsmp, m) = tmp;
        end
        initcomp{k} = initcomp{k} - repmat(mean(initcomp{k}),nsample,1);
      end
    else
      assert(numel(initcomp)==nchan);
    end

    if ~iscell(jitter)
      jitter = repmat({jitter}, [1 nchan]);
    end

    for k = 1:numel(jitter)
      if ~isempty(jitter{k})
        if size(jitter{k},1)~=size(initcomp{k},2), jitter{k} = repmat(jitter{k}(1,:),[size(initcomp{k},2) 1]); end
      end
    end

    % initialize the output data
    dataout = removefields(data, 'cfg');
    for k = 1:numel(data.trial)
      dataout.trial{k}(:) = nan;
    end

    % initialize the struct that will contain the output parameters
    params = struct([]);

    % do the actual computations
    for k = 1:nchan
      % preprocessing data
      tmp     = cellrowselect(data.trial,k);
      chandat = cat(1,tmp{:});

      % baseline correction
      chandat = ft_preproc_baselinecorrect(chandat, nearest(data.time{1}, -inf), nearest(data.time{1}, 0));

      % do zero-padding and FFT to the signal and initial waveforms
      npad         = cfg.aseo.pad*fsample;    % length of data + zero-padding number
      nfft         = 2.^(ceil(log2(npad)))*2;
      initcomp_fft = fft(initcomp{k}, nfft);  % Fourier transform of the initial waveform
      chandat_fft  = fft(chandat', nfft);     % Fourier transform of the signal

      cfg.aseo.jitter = jitter{k};
      output       = ft_singletrialanalysis_aseo(cfg, chandat_fft, initcomp_fft);

      params(k).latency    = output(end).lat_est./fsample;
      params(k).amplitude  = output(end).amp_est;
      params(k).components = output(end).erp_est;
      params(k).rejectflag = output(end).rejectflag;
      params(k).noise      = output(end).noise;

      for m = 1:numel(data.trial)
        if output(end).rejectflag(m)==0
          switch cfg.output
            case 'model'
              dataout.trial{m}(k,:) = data.trial{m}(k,:)-output(end).residual(:,m)';
            case 'residual'
              dataout.trial{m}(k,:) = output(end).residual(:,m)';
          end
        end
      end
    end

case 'gbve'
  ft_hastoolbox('lagextraction', 1);
  ft_hastoolbox('eeglab',        1); % because the low-level code might use a specific moving average function from EEGLAB
  ft_hastoolbox('cellfunction',  1);

  if ~isfield(cfg, 'gbve'), cfg.gbve = []; end
  cfg.gbve.NORMALIZE_DATA    = ft_getopt(cfg.gbve, 'NORMALIZE_DATA',     true);
  cfg.gbve.CENTER_DATA       = ft_getopt(cfg.gbve, 'CENTER_DATA',        false);
  cfg.gbve.USE_ADAPTIVE_SIGMA= ft_getopt(cfg.gbve, 'USE_ADAPTIVE_SIGMA', false);
  cfg.gbve.sigma             = ft_getopt(cfg.gbve, 'sigma',    0.01:0.01:0.2);
  cfg.gbve.distance          = ft_getopt(cfg.gbve, 'distance', 'corr2');
  cfg.gbve.alpha             = ft_getopt(cfg.gbve, 'alpha',    [0 0.001 0.01 0.1]);
  cfg.gbve.exponent          = ft_getopt(cfg.gbve, 'exponent', 1);
  cfg.gbve.use_maximum       = ft_getopt(cfg.gbve, 'use_maximum', 1); % consider the positive going peak
  cfg.gbve.show_pca          = ft_getopt(cfg.gbve, 'show_pca',          false);
  cfg.gbve.show_trial_number = ft_getopt(cfg.gbve, 'show_trial_number', false);
  cfg.gbve.verbose           = ft_getopt(cfg.gbve, 'verbose',           true);
  cfg.gbve.disp_log          = ft_getopt(cfg.gbve, 'disp_log',          false);
  cfg.gbve.latency           = ft_getopt(cfg.gbve, 'latency',  [-inf inf]);
  cfg.gbve.xwin              = ft_getopt(cfg.gbve, 'xwin',     1); % default is a bit of smoothing

  nchan = numel(data.label);
  ntrl  = numel(data.trial);

  tmin  = nearest(data.time{1}, cfg.gbve.latency(1));
  tmax  = nearest(data.time{1}, cfg.gbve.latency(2));

  % initialize the struct that will contain the output parameters
  dataout = removefields(data, 'cfg');
  params  = struct([]);
  for k = 1:nchan
    % preprocessing data
    options = cfg.gbve;

    fprintf('--- Processing channel %d\n',k);

    tmp     = cellrowselect(data.trial,k);
    chandat = cat(1,tmp{:});
    points  = chandat(:,tmin:tmax);

    % perform a loop across alpha values, cross validation
    alphas = options.alpha;

    if length(alphas) > 1 % Use Cross validation error if multiple alphas are specified
      best_CVerr = -Inf;

      K = 5;
      disp(['--- Running K Cross Validation (K = ',num2str(K),')']);

      block_idx = fix(linspace(1, ntrl, K+1)); % K cross validation
      for jj=1:length(alphas)
        options.alpha = alphas(jj);

        CVerr = 0;
        for kk = 1:K
          bidx = block_idx(kk):block_idx(kk+1);
          idx = 1:ntrl;
          idx(bidx) = [];

          data_k       = chandat(idx,:);
          points_k     = points(idx,:);
          [order,lags] = extractlag(points_k,options);

          data_reordered = data_k(order,:);
          lags           = lags + tmin;
          [data_aligned, dum] = perform_realign(data_reordered, data.time{1}, lags);
          data_aligned(~isfinite(data_aligned)) = nan;
          ep_evoked = nanmean(data_aligned);
          ep_evoked = ep_evoked ./ norm(ep_evoked);

          data_k = chandat(bidx,:);
          data_norm = sqrt(sum(data_k.^2,2));
          data_k = diag(1./data_norm)*data_k;
          data_k(data_norm==0,:) = 0;

          for pp=1:length(bidx)
            c     = xcorr(ep_evoked,data_k(pp,:));
            CVerr = CVerr + max(c(:));
          end
        end

        CVerr = CVerr/ntrl;

        if CVerr > best_CVerr
          best_CVerr = CVerr;
          best_alpha = alphas(jj);
        end
      end
      options.alpha = best_alpha;
    end

    if options.use_maximum
      [order,lags] = extractlag( points, options );
    else
      [order,lags] = extractlag( -points, options );
    end
    disp(['---------- Using alpha = ',num2str(options.alpha)]);
    data_reordered = chandat(order,:);
    lags = lags + tmin;
    [data_aligned] = perform_realign(data_reordered, data.time{1}, lags );
    data_aligned(~isfinite(data_aligned)) = nan;

    [dum,order_inv] = sort(order);
    lags_no_order = lags(order_inv);
    data_aligned  = data_aligned(order_inv,:);

    params(k).latency = data.time{1}(lags_no_order)';
    switch cfg.output
      case 'model'
        tmp = mat2cell(data_aligned, ones(1,size(data_aligned,1)), size(data_aligned,2))';
        dataout.trial = cellrowassign(dataout.trial, tmp, k);
      case 'residual'
        % to be done
        error('not yet implemented');
    end
  end

end

dataout.params = params;
dataout.cfg    = cfg;

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous   data
ft_postamble provenance dataout
ft_postamble history    dataout
ft_postamble savevar    dataout


text:
function img_data=read_img_rgb(img_file, data_type)

assert(strcmp(data_type, 'uint8'));

try

	[img_data, map] = imread(img_file);

catch err_info

	disp(err_info);
	disp(img_file);
	error('image read error');

end

if size(img_data, 4)>1

	error('cannot read images with multi-frames');

    % only load the first frame
    img_data=img_data(:,:,:, 1);
end

if ~isempty(map)
	img_data = ind2rgb(img_data,map);
end

if size(img_data, 3)==1

    [img_data map]=gray2ind(img_data);
    img_data=ind2rgb(img_data, map);

end

if strcmp(data_type, 'uint8')
    img_data=im2uint8(img_data);
end
if strcmp(data_type, 'double')
    img_data=im2double(img_data);
end
if strcmp(data_type, 'single')
    img_data=im2single(img_data);
end

try
    assert(size(img_data, 3)==3);
catch error_msg
    disp(error_msg);
    dbstack;
    keyboard;
end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function qout = ch_qconj(qin)
qout = [qin(1); -qin(2:4)];

author: yandld


text:
addpath('..\FW_SPCP\func');
addpath('D:\solvers\PROPACK');
addpath('D:\solvers\cvx');
cvx_setup; warning off;

author: andrewssobral

repo: lrslibrary


text:
function active_caffe_mex(gpu_id, caffe_version)
% active_caffe_mex(gpu_id, caffe_version)
% --------------------------------------------------------
% Faster R-CNN
% Copyright (c) 2015, Shaoqing Ren
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

    % set gpu in matlab
    g = gpuDevice(gpu_id);

    if ~exist('caffe_version', 'var') || isempty(caffe_version)
        caffe_version = 'caffe';
    end
    cur_dir = pwd;
    caffe_dir = fullfile(pwd, 'external', 'caffe', 'matlab', caffe_version);

    if ~exist(caffe_dir, 'dir')
        warning('Specified caffe folder (%s) is not exist, change to default one (%s)', ...
            caffe_dir, fullfile(pwd, 'external', 'caffe', 'matlab'));
        caffe_dir = fullfile(pwd, 'external', 'caffe', 'matlab');
    end

    addpath(genpath(caffe_dir));
    cd(caffe_dir);
    caffe.set_device(gpu_id-1);
    cd(cur_dir);
end

author: feichtenhofer

repo: Detect-Track


text:
function ps = getSplitMergeFeatChoiceProbs( F, stateSeq, data_struct, model, jj, ki )
%function ps = getSplitMergeFeatChoiceProbs( F, jj, ki )
ps = F(jj,:);

availFeatIDs = find( F(jj,:) );
relFeatIDs = union( availFeatIDs, ki);

%Ustats  = getStateSeqSuffStats_C1( stateSeq, F, data_struct, model, 1:size(F,1), 1:size(F,2) );
%Ustats.Nkv = Ustats.Nkv(relFeatIDs,:);
%myTheta = getTheta_PosteriorMean( myTheta, Ustats, model.obsModel, relFeatIDs );

Xstats = getXSuffStats( F, stateSeq, data_struct, model, 1:size(F,1),  relFeatIDs );
myTheta = setThetaToPosteriorMean(  [], Xstats, model.obsModel, relFeatIDs );

EPS = 1e-9;
switch model.obsModel.type
    case 'Multinomial'
        pHat = exp( myTheta.logp )';
        pHat = bsxfun( @rdivide, pHat, sum(pHat,1) );
        D = calcChiSqDistanceBetter( pHat(:,ki), pHat );
    case 'Gaussian'
        if size( myTheta.Mu, 1) > 1
            D = pdist( myTheta.Mu, 'euclidean');
            D = squareform(D);
            D = D(ki,:);
        else
            D = 0; % Default distance from self is zero.
        end
        D = D./(EPS+max(max(D)) ); % Ensure D varies from 0 to 1.

    case 'AR-Gaussian'
        [D DR K] = size( myTheta.A );

        %Avec = -100*ones(K,D);
        Avec = zeros(K, D*DR);
        for kk = relFeatIDs
            %Sig = myTheta.invSigma(:,:,kk) \ eye(D);
            %Avec(kk,:) = diag( Sig );
            Avec(kk,:) = reshape( myTheta.A(:,:,kk), 1, D*DR );
        end
        if K > 1
            D = pdist( Avec, 'euclidean');
            D = squareform(D);
            D = D(ki,:);
        else
            D = 0; % Default distance from self is zero.
        end
        D = D./(EPS+max(max(D)) );
        % Ensure D varies from 0 to 1.
        %  but is never quite 0 (must be at least EPS)
end

D( ~F(jj,:)  ) = Inf;

ps = 1 - D;
ps( ~F(jj,:) ) = 0;

ps = ps./sum( ps );

if F( jj, ki ) && length( ps( ps>0 ) ) > 1
    % If object jj possesses feature ki,
    %    then let prob kj =ki = 2/3
    %                  kj~=ki = 1/3 (unif. distrib. over all jj's feats)
    ps( ki ) = 2*( sum( ps( [1:ki-1 ki+1:end] ) ) );
end

ps = ps./sum(ps);

% NOTE: Alternate distance metric
%sig = eps+min(D(availFeatIDs));
%ps = exp(-D./sig );

end

author: michaelchughes

repo: NPBayesHMM


text:
function [ cost, grad, numTotal, pred_cell ] = drdae_discrim_joint_kl_obj_gpu( theta, eI, data_cell, targets_cell, mixture_spectrum, fprop_only, pred_out)
% Copyright (c) 2014-present University of Illinois at Urbana-Champaign
% All rights reserved.
%
% Developed by:     Po-Sen Huang, Paris Smaragdis
%                   Department of Electrical and Computer Engineering
%                   Department of Computer Science
%
% discrim. training + joint masking using MATLAB GPU Toolbox
%
%PRNN_OBJ MinFunc style objective for Deep Recurrent Denoising Autoencoder
%   theta is the full parameter vector
%   eI contains experiment / network architecture
%   data_cell is a cell array of matrices. Each a distinct length is a cell
%             entry. Each matrix has a time series example in each column
%   targets_cell is parallel to data, but contains the labels for each time
%   fprop_only is a flag that only computes the cost, no gradient
%   numTotal is total number of frames evaluated
%   pred_out is a binary flag for whether pred_cell is populated
%            pred_cell only filled properly when utterances one per cell

import parallel.gpu.GPUArray

gtheta = gpuArray(theta);

%% Debug: Turns this into an identity-function for debugging rest of system
if isfield(eI, 'objReturnsIdentity') && eI.objReturnsIdentity
    cost = 0; grad = 0; numTotal = 0;
    for l = 1:numel(data_cell)
        numUtterances = size(data_cell{l}, 2);
        original_vector = reshape(data_cell{l}, eI.winSize*eI.featDim, []);
        midPnt = ceil(eI.winSize/2);
        original_vector = original_vector((midPnt-1)*14+1 : midPnt*14, :);
        pred_cell{l} = reshape(original_vector, [], numUtterances);
    end
    return;
end

%if isempty(return_activation),
  return_activation = 0;
%end

%% Load data from globals if not passed in (happens when run on RPC slave)
global g_data_cell;
global g_targets_cell;
isSlave = false;
if isempty(data_cell)
    data_cell = g_data_cell;
    targets_cell = g_targets_cell;
    isSlave = true;
end;
pred_cell = cell(1,numel(data_cell));
act_cell = cell(1,numel(data_cell));
%% default short circuits to false
if ~isfield(eI, 'shortCircuit')
    eI.shortCircuit = 0;
end;

%% default dropout to false
if ~isfield(eI, 'dropout')
  eI.dropout = 0;
end;

%% setup weights and accumulators
[stack, W_t] = rnn_params2stack(gtheta, eI);
cost = 0; numTotal = 0;
outputDim = eI.layerSizes(end);

%% setup structures to aggregate gradients
stackGrad = cell(1,numel(eI.layerSizes));
% W_t_grad = zeros(size(W_t));
if isfield(eI, 'fullRNN') && eI.fullRNN==1
   W_t_grad = cell(1,numel(eI.layerSizes)-1);
    for l = 1:numel(eI.layerSizes)-1
        W_t_grad{l}.W = gpuArray(zeros(size(W_t{l}.W)));
    end
else
   W_t_grad = gpuArray(zeros(size(W_t)));
end

for l = 1:numel(eI.layerSizes)
    stackGrad{l}.W = gpuArray(zeros(size(stack{l}.W)));
    stackGrad{l}.b = gpuArray(zeros(size(stack{l}.b)));
end
if eI.shortCircuit
    stackGrad{end}.W_ss = gpuArray(zeros(size(stack{end}.W_ss)));
end;
%% check options
if ~exist('fprop_only','var')
    fprop_only = false;
end;
if ~exist('pred_out','var')
    pred_out = false;
end;

% DROPOUT: vector of length of hidden layers with 0 or 1
% (to drop or keep activation unit) with prob=0.5
hActToDrop = cell(numel(eI.layerSizes-1),1);
for i=1:numel(eI.layerSizes)-1
 if eI.dropout
   hActToDrop{i} = gpuArray(1/eI.dropout * binornd(1,eI.dropout, eI.layerSizes(i),1));
   %hActToDrop{i} = gpuArray(round(rand(eI.layerSizes(i),1)));
 else
   hActToDrop{i} = gpuArray(ones(eI.layerSizes(i),1));
 end
end

%% loop over each distinct length
for c = 1:numel(data_cell)

    targets = {};
    if ~isempty(targets_cell), targets = targets_cell{c}; end;

    mbsz=min(size(data_cell{c},2), floor(200*1024^2/size(data_cell{c},1)/8)); % 200mb max
    if mbsz==0, continue; end
    nbat = floor(size(data_cell{c},2)/mbsz)+1;

     % convert different mini-bats to GPUs
    for bat=1:nbat

    data = gpuArray(data_cell{c}(:,1+(bat-1)*mbsz: min(size(data_cell{c},2), bat*mbsz)));
    targets = gpuArray(targets_cell{c}(:,1+(bat-1)*mbsz:  min(size(targets_cell{c},2), bat*mbsz)));
    mixtures=gpuArray(mixture_spectrum{c}(:,1+(bat-1)*mbsz: min(size(mixture_spectrum{c},2), bat*mbsz)));

    uttPred = [];
    T =size(data,1) / eI.inputDim;
    % store hidden unit activations at each time instant
    hAct = cell(numel(eI.layerSizes)-1, T);
    for t = 1:T
        %% forward prop all hidden layers
        for l = 1:numel(eI.layerSizes)-1
            if l == 1
                hAct{1,t} = stack{1}.W * data((t-1)*eI.inputDim+1:t*eI.inputDim, :);
            else
                hAct{l,t} = stack{l}.W * hAct{l-1,t};
            end;
            hAct{l,t} = hAct{l,t}+ repmat(stack{l}.b, 1, size(hAct{l,t},2));
            % temporal recurrence. limited to single layer for now
            if t > 1
                if isfield(eI, 'fullRNN') && eI.fullRNN==1
                    hAct{l,t} = hAct{l,t} + W_t{l}.W * hAct{l,t-1};
                elseif l == eI.temporalLayer
                    hAct{l,t} = hAct{l,t} + W_t * hAct{l,t-1};
                end
            end;

            % nonlinearity
            if strcmpi(eI.activationFn,'tanh')
                hAct{l,t} = tanh(hAct{l,t});
            elseif strcmpi(eI.activationFn,'logistic')
                hAct{l,t} = 1./(1+exp(-hAct{l,t}));
            elseif strcmpi(eI.activationFn,'RELU')
%                  maxg([ GPUsingle(zeros(size(hAct{1,t}))), hAct{l,t}], hAct{l,t});
%                 hAct{l,t} = maxg(0,hAct{l,t});
                hAct{l,t} = max(0,hAct{l,t});
            else
                error('unrecognized activation function: %s',eI.activationFn);
            end;
            %dropout (hActToDrop will be all ones if no dropout specified)
            hAct{1,t} = bsxfun(@times, hAct{1,t}, hActToDrop{l});
%              hAct{1,t} = hAct{1,t}.*hActToDrop{l};
        end;
        % forward prop top layer not done here to avoid caching it
    end;
    %% compute cost and backprop through time
    if  eI.temporalLayer
        if isfield(eI, 'fullRNN') && eI.fullRNN==1
            delta_t = cell(1, numel(eI.layerSizes)-1);
            for l = 1:numel(eI.layerSizes)-1
                delta_t{l} = gpuArray(zeros(eI.layerSizes(l),size(data,2)));
            end
        else
            delta_t = gpuArray(zeros(eI.layerSizes(eI.temporalLayer),size(data,2)));
        end
    end;

    y1_dim= 1:outputDim/2;
    y2_dim= outputDim/2+1:outputDim;

    for t = T:-1:1
        l = numel(eI.layerSizes);
        %% forward prop output layer for this timestep
        curPred = bsxfun(@plus, stack{l}.W * hAct{l-1,t}, stack{l}.b);

        if eI.outputnonlinear==1,
           if strcmpi(eI.activationFn,'tanh')
                curPred = tanh(curPred);
           elseif strcmpi(eI.activationFn,'logistic')
                curPred = 1./(1+exp(-curPred));
           elseif strcmpi(eI.activationFn,'RELU')
                curPred = max(0,curPred);
           else
                error('unrecognized activation function: %s',eI.activationFn);
           end
        end

        mixture=mixtures((t-1)*numel(y1_dim)+1:t*numel(y1_dim),:);
        a1 = curPred(y1_dim,:); a2 = curPred(y2_dim,:);

        const=eI.const;%1e-8 ;
        const2=eI.const2;% 1e-3;

        if strcmp(eI.opt,'softlinear'),
            y1= (a1)./((a1)+(a2)+1e-10).* mixture;
            y2= (a2)./((a1)+(a2)+1e-10).* mixture;
        elseif strcmp(eI.opt,'softabs'),
            abs_a1= abs(a1); abs_a2= abs(a2);
            y1= abs_a1./(abs_a1+abs_a2+1e-10).* mixture;
            y2= abs_a2./(abs_a1+abs_a2+1e-10).* mixture;
        elseif strcmp(eI.opt,'softabs_const')|| strcmp(eI.opt,'softabs_kl_const'),
            y1= abs(a1)./(abs(a1)+abs(a2)+const).* mixture;
            y2= abs(a2)./(abs(a1)+abs(a2)+const).* mixture;
        elseif strcmp(eI.opt, 'softquad')
            y1= (a1.^2)./((a1.^2)+(a2.^2)+1e-10).* mixture;
            y2= (a2.^2)./((a1.^2)+(a2.^2)+1e-10).* mixture;
        else
        end

        weighted_curPred=[y1; y2];
        % add short circuit to regression prediction if model has it
        if eI.shortCircuit
            weighted_curPred = weighted_curPred + stack{end}.W_ss ...
                * data((t-1)*eI.inputDim+1:t*eI.inputDim, :);
        end;
        if pred_out, uttPred = [weighted_curPred; uttPred]; end;
        % skip loss computation if no targets given
        if isempty(targets), continue; end;

        curTargets = targets((t-1)*outputDim+1:t*outputDim, :);
        curTargets_neg = [curTargets(outputDim/2+1:outputDim,:); curTargets(1:outputDim/2,:)];

        y_t = (1- eI.r) * weighted_curPred + eI.r * curTargets_neg - curTargets;

        ya_ta= y_t(y1_dim,:);
        yb_tb= y_t(y2_dim,:);

        if strcmp(eI.opt,'softlinear'),
            delta_y1 =  (ya_ta-yb_tb).* y2./(a1+a2+1e-10);
            delta_y2 = (-ya_ta+yb_tb) .* y1./ (a1+a2+1e-10);
        elseif strcmp(eI.opt,'softabs'),
            delta_y1 =  (ya_ta-yb_tb).* y2./(abs(a1)+abs(a2)+1e-10);
            delta_y2 = (-ya_ta+yb_tb) .* y1./ (abs(a1)+abs(a2)+1e-10);

            delta_y1(a1<0) = -delta_y1(a1<0);
            delta_y2(a2<0) = -delta_y2(a2<0);
        elseif strcmp(eI.opt,'softabs_const'),
             const_div=const./((abs(a1)+abs(a2)+const).^2).* mixture;
             delta_y1 =  (ya_ta-yb_tb).* y2./(abs(a1)+abs(a2)+const);
             delta_y1 = delta_y1+ ya_ta.* const_div;
             delta_y2 = (-ya_ta+yb_tb) .* y1./ (abs(a1)+abs(a2)+const);
             delta_y2= delta_y2+ yb_tb.* const_div;

             delta_y1(a1<0) =  -delta_y1(a1<0);
             delta_y2(a2<0) =  -delta_y2(a2<0);
        elseif strcmp(eI.opt,'softabs_kl_const'),
             y_target_a = curTargets(y1_dim, :);
             y_target_b = curTargets(y2_dim, :);
             y_target_neg_a = y_target_b;
             y_target_neg_b = y_target_a;

             y_pred_a = weighted_curPred(y1_dim, :);
             y_pred_b = weighted_curPred(y2_dim, :);

             const_div=const./((abs(a1)+abs(a2)+const).^2).* mixture;

             delta_y1 =  (-y_target_a./(y_pred_a+const2) + y_target_b./(y_pred_b+const2)).* y2./(abs(a1)+abs(a2)+const);
             delta_y1 = delta_y1+  (-y_target_a./(y_pred_a+const2)+1).* const_div;

             delta_y2 =  (y_target_a./(y_pred_a+const2) - y_target_b./(y_pred_b+const2)).* y1./(abs(a1)+abs(a2)+const);
             delta_y2= delta_y2+ (-y_target_b./(y_pred_b+const2)+1) .* const_div;

             % discrim part
             delta_y1 =  delta_y1- eI.r* (-y_target_neg_a./(y_pred_a+const2) + y_target_neg_b./(y_pred_b+const2)).* y2./(abs(a1)+abs(a2)+const);
             delta_y1 = delta_y1-   eI.r* (-y_target_neg_a./(y_pred_a+const2)+1).* const_div;

             delta_y2 = delta_y2- eI.r*(y_target_neg_a./(y_pred_a+const2) - y_target_neg_b./(y_pred_b+const2)).* y1./(abs(a1)+abs(a2)+const);
             delta_y2= delta_y2-  eI.r*(-y_target_neg_b./(y_pred_b+const2)+1) .* const_div;

             delta_y1(a1<0) =  -delta_y1(a1<0);
             delta_y2(a2<0) =  -delta_y2(a2<0);

        elseif strcmp(eI.opt, 'softquad')
            delta_y1 =  (ya_ta-yb_tb).* (2*a1.*y2)./ (a1.^2+a2.^2+1e-10); %  y1= (a1.^2)./((a1.^2)+(a2.^2)+1e-8);%.* mixture;
            delta_y2 =  (-ya_ta+yb_tb).* (2*a2.*y1)./ (a1.^2+a2.^2+1e-10);%         y2= (a2.^2)./((a1.^2)+(a2.^2)+1e-8);%.* mixture;
        else
        end

        delta = [ delta_y1; delta_y2 ];
        if strcmp(eI.opt,'softlinear') || strcmp(eI.opt,'softabs') || strcmp(eI.opt, 'softquad') || strcmp(eI.opt,'softabs_const'),
             cost = cost + 0.5 * ( sum( sum((weighted_curPred - curTargets).*(weighted_curPred - curTargets))) ...
                -  eI.r* sum( sum((weighted_curPred - curTargets_neg).*(weighted_curPred - curTargets_neg))));
        elseif strcmp(eI.opt,'softabs_kl_const'),
             cost = cost +...
             sum(sum( curTargets.*log( curTargets./(weighted_curPred + const2) + const2 )-curTargets+ weighted_curPred +const2))...
            -eI.r* sum(sum( curTargets_neg.*log( curTargets_neg./(weighted_curPred + const2) + const2 )-curTargets_neg+ weighted_curPred +const2));
        else

        end
        if eI.outputnonlinear==1,
             if strcmpi(eI.activationFn,'tanh')
                delta = delta .* (1 -curPred.^2);
            elseif strcmpi(eI.activationFn,'logistic')
                delta = delta .* curPred .* (1 - curPred);
            elseif strcmpi(eI.activationFn,'RELU')
                delta = delta .* double(curPred>0);
            else
                error('unrecognized activation function: %s',eI.activationFn);
            end;
        end

        if fprop_only, continue; end;
        %% regression layer gradient and delta
        stackGrad{l}.W = stackGrad{l}.W + delta * hAct{l-1,t}';
        stackGrad{l}.b = stackGrad{l}.b + sum(delta,2);
        % short circuit layer
        if eI.shortCircuit
            stackGrad{end}.W_ss = stackGrad{end}.W_ss + delta ...
                * data((t-1)*eI.inputDim+1:t*eI.inputDim, :)';
        end;
        delta = stack{l}.W' * delta;
        %% backprop through hidden layers
        for l = numel(eI.layerSizes)-1:-1:1
            % aggregate temporal delta term if this is the recurrent layer
            if isfield(eI, 'fullRNN') && eI.fullRNN==1
                delta = delta + delta_t{l};
            elseif l == eI.temporalLayer
                delta = delta + delta_t;
            else
            end
            % push delta through activation function for this layer
            % tanh unit choice assumed
            if strcmpi(eI.activationFn,'tanh')
                delta = delta .* (1 - hAct{l,t}.^2);
            elseif strcmpi(eI.activationFn,'logistic')
                delta = delta .* hAct{l,t} .* (1 - hAct{l,t});
            elseif strcmpi(eI.activationFn,'RELU')
                delta = delta .* gpuArray(hAct{l,t}>0);
            else
                error('unrecognized activation function: %s',eI.activationFn);
            end;

            % gradient of bottom-up connection for this layer
            if l > 1
                stackGrad{l}.W = stackGrad{l}.W + delta * hAct{l-1,t}';
            else
                stackGrad{l}.W = stackGrad{l}.W + delta * data((t-1)*eI.inputDim+1:t*eI.inputDim, :)';
            end;
            % gradient for bias
            stackGrad{l}.b = stackGrad{l}.b + sum(delta,2);

            % compute derivative and delta for temporal connections
            if t > 1
                if isfield(eI, 'fullRNN') && eI.fullRNN==1
                     W_t_grad{l}.W = W_t_grad{l}.W + delta * hAct{l,t-1}';
                     % push delta through temporal weights
                     delta_t{l} = W_t{l}.W' * delta;
                elseif l == eI.temporalLayer
                     W_t_grad = W_t_grad + delta * hAct{l,t-1}';
                     % push delta through temporal weights
                     delta_t = W_t' * delta;
                end;
            end
            % push delta through bottom-up weights
            if l > 1
                delta = stack{l}.W' * delta;
            end;
        end
    end
    pred_cell{c}=double(gather(uttPred));

    % Return the activations for this utterance.
    if return_activation,
      act_cell{c} = cell2mat(hAct);
    end
    % keep track of how many examples seen in total
    numTotal = numTotal + T * size(targets,2);
    end
end

%% stack gradients into single vector and compute weight cost
wCost = numTotal * eI.lambda * sum(gtheta.^2);
grad = rnn_stack2params(stackGrad, eI, W_t_grad, true);
grad = grad + 2 * numTotal * eI.lambda * gtheta;

%% clipping
if isfield(eI,'clip') && eI.clip~=0, % if eI.clip==0, no clip
  if eI.clip > 0 % method one -clip the whole
      norm_grad = norm(grad);
      fprintf('norm_grad:%f\n', norm_grad);
      % avoid numerial problem
      if norm_grad <0 || norm_grad > 1e15 || isnan(norm_grad) || isinf(norm_grad),
          grad = zeros(size(grad));
          fprintf('set gradient to zeros\n');
      end
      if norm_grad > eI.clip
         grad = eI.clip * grad/ norm_grad;
      end
  else % method two - clip each entry
      clip_value = -1*eI.clip;
      grad(grad > clip_value)=clip_value;
      grad(grad < -clip_value)=-clip_value;
  end
end
%%
cost = gather((cost));
grad = gather((grad));
wCost = gather((wCost));

avCost = cost/numTotal;
avWCost = wCost/numTotal;
cost = cost + wCost;

% print output
if ~isSlave && ~isempty(targets_cell)
    fprintf('loss:  %f  wCost:  %f \t',avCost, avWCost);

    if isfield(eI, 'fullRNN') && eI.fullRNN==1
        fprintf('wNorm: %f  rNorm: %f  oNorm: %f\n',sum(stack{1}.W(:).^2),...
            sum(W_t{1}.W(:).^2), sum(stack{end}.W(:).^2));
    else
        fprintf('wNorm: %f  rNorm: %f  oNorm: %f\n',sum(stack{1}.W(:).^2),...
            sum(W_t(:).^2), sum(stack{end}.W(:).^2));
    end
end;

author: posenhuang

repo: deeplearningsourceseparation


text:
%  ADDTSPLOT - Add time series to current figure
%
%  ADDTSPLOT(Y) adds time series in the rows of Y to the current
%  figure which should be created by TSPLOT. Y can be followed by
%  parameter/value pairs to specify additional properties of the
%  lines.
%
%  See also TSPLOT

%  This software is provided "as is", without warranty of any kind.
%  Alexander Ilin, Tapani Raiko

function hax = addtsplot( y, varargin )

if nargin >= 2 & isnumeric(varargin{1})
    x = y;
    y = varargin{1};
    varargin = {varargin{2:end}};
else
    x = 1:size(y,2);
end

nts = size(y,1);

hax = findobj( gcf, 'type', 'axes' );
if length(hax) ~= nts
    tsplot( y, varargin{:} );
    return
end

pos = zeros( length(hax), 4 );
for i = 1:length(hax)
    pos(i,:) = get( hax(i), 'position' );
end
[tmp,I] = sort( -pos(:,2) );
hax = hax(I);

for i = 1:nts
    axes(hax(i))

    hold on
    plot( x, y(i,:), varargin{:} )

end

if nargout == 0
    clear hax
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% main script to launch the estimation
% of the non-linear parameters by using the CalTech
% calibration toolbox and the output from the Svoboda's
% Multicamera self-calibration
%
% How to create the input data:
% 1) Run the MultiCamSelfCam
% 2) Run the MultiCamValidation
%
% $Id: gorad.m,v 2.0 2003/06/19 12:06:00 svoboda Exp $

clear all;

% TODO: isn't this unnecessary now?
addpath ../MultiCamSelfCalib/Cfg

% Read configuration from whatever is specified on command-line (via --config=FILENAME)
config = read_configuration();

% if problem with desactivated images -> some problems with the estimation in general
desactivated_images = [];

idxcams = config.cal.cams2use;
selfcalib.goradproblem = 0;

for i = idxcams,
  [X_1,x_1] = preparedata(sprintf(config.files.points4cal,i));
  go_calib_optim_iter
  if any(isnan(param))
	  % when the iteration fails insert null distortion
	  % it is better than nonsense
	  kc(1:4) = [0,0,0,0];
	  selfcalib.goradproblem=1;
  else
	  visualize_distortions
  end

  disp(sprintf('***** camera %d **********************************',i))
  %
  outputfile = sprintf(config.files.rad,i);

     fprintf(1,'\nExport of intrinsic calibration data to blue-c configuration file\n');
     % outputfile = input('File basename: ', 's');
     configfile = outputfile;
     disp(['Writing ' configfile]);

     fid = fopen(configfile, 'w');

     fprintf(fid, 'K11 = %.16f\n', KK(1,1));
     fprintf(fid, 'K12 = %.16f\n', KK(1,2));
     fprintf(fid, 'K13 = %.16f\n', KK(1,3));
     fprintf(fid, 'K21 = %.16f\n', KK(2,1));
     fprintf(fid, 'K22 = %.16f\n', KK(2,2));
     fprintf(fid, 'K23 = %.16f\n', KK(2,3));
     fprintf(fid, 'K31 = %.16f\n', KK(3,1));
     fprintf(fid, 'K32 = %.16f\n', KK(3,2));
     fprintf(fid, 'K33 = %.16f\n\n', KK(3,3));

     fprintf(fid, 'kc1 = %.16f\n', kc(1));
     fprintf(fid, 'kc2 = %.16f\n', kc(2));
     fprintf(fid, 'kc3 = %.16f\n', kc(3));
     fprintf(fid, 'kc4 = %.16f\n', kc(4));

     status = fclose(fid);

disp('Press any key to continue'),  pause

%%%
% clear already estimated parameters
clear fc kc alpha_c cc
end

author: strawlab

repo: MultiCamSelfCal


text:
function [predict_label,accuracy,prob_estimates] = jdsvm_w(x_train,y_train,x_test,y_test,w,c,g,kernel_type)
%%jdsvm:packaged libsvm,just for easier use.
%%==============================================================================
%%input:
%%------x_train,...,y_test  :   training and testing sets.                               [required]
%%------c                   :   punishment factor(the more,the heavier),default is 1.    [not required]
%%------g                   :   gamma function,default is 1/(#features).                 [not required]
%%------kernel_type         :   kernel type,value is int,the choices are:                [not required]
%%%%----------------0       :   linear
%%%%----------------1       :   polynomial
%%%%----------------2       :   rbf(default)
%%%%----------------3       :   sigmoid
%%------w                   :   weight vector

%%output:
%%------predict_label       :   predicted label vector for test case.
%%------accuracy            :   accuracy
%%------prob_estimates      :   prob estimates
%%==============================================================================
    predict_label = [];
    accuracy = 0;
    prob_estimates = [];
    c_default = 1;
	g_default = 1 / size(x_train,2);
	kernel_type_default = 2;
    weight_default = '';

	c_para = c_default;
	g_para = g_default;
	kernel_para = kernel_type_default;
    weight_para = weight_default;

	switch nargin
		case 4
			c_para = c_default;
			g_para = g_default;
			kernel_para = kernel_type_default;
            weight_para = weight_default;
		case 5
			c_para = c_default;
			g_para = g_default;
			kernel_para = kernel_type_default;
            weight_para = w;
		case 6
			c_para = c;
			g_para = g_default;
			kernel_para = kernel_type_default;
            weight_para = w;
		case 7
			c_para = c;
			g_para = g;
			kernel_para = kernel_type;
            weight_para = w;
        case 8
			c_para = c;
			g_para = g;
			kernel_para = kernel_type;
            weight_para = w;
        otherwise
            fprintf('++++++Fatal error!Please check the input!');
            return
	end

	options = ['-c ',num2str(c_para),' -g ',num2str(g_para),' -t ',num2str(kernel_para),weight_para];
	model = libsvmtrain(y_train,x_train,options);
	[predict_label,accuracy,prob_estimates] = libsvmpredict(y_test,x_test,model);
	accuracy = accuracy(1) / 100;

end

author: jindongwang

repo: activityrecognition


text:
filename='RVE_Square_Triangle_FineFine';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'horizontal';
cost={'chomog_alphabeta'};
weights=[1 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.5;
Perimeter_target=1;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
micro.alpha =[1 1 0]';
micro.beta =[1 1 0]';

% For all tests
plotting = false;
printing = false;
monitoring = false;
maxiter = 3;

author: SwanLab

repo: Swan