text:
classdef opti < handle
%OPTI  Create an OPTI object for Linear, Quadratic & Nonlinear Optimization
%
%   optObj = opti('param1',value1,'param2',value2,...) creates an OPTI
%   object with the parameters 'param' set to their corresponding
%   values in 'value'. Parameters not specified will be set to the OPTI
%   default.
%
%   optObj = opti(...,'options',opts) creates the object with optiset
%   options supplied to the OPTI constructor.
%
%   optObj = opti(prob) creates an opti object based on the optimization
%   problem specified in prob. Use 'optiprob' to generate the structure
%   'prob'.
%
%   optObj = opti(prob,opts) creates the opti object with specified
%   options. Use 'optiset' to generate the options structure 'opts'.
%
%   opti() prints a list of all possible fields and their function.
%
%
%   Typical usage examples:
%
%   1) LP
%   optObj = opti('f',f,'ineq',A,b,'bounds',lb,ub)
%
%   2) Bounded NLP
%   optObj = opti('obj',obj,'bounds',lb,ub,'x0',x0)
%
%   3) Option Setting
%   optObj = opti('obj',obj,'options',optiset('solver','nomad'))
%
%
%   For detailed documentation see the Wiki:
%      https://inverseproblem.co.nz/OPTI
%
%
%   See also opti.solve optiset opti.plot
%
%   Copyright (C) 2011-2018 Jonathan Currie (www.inverseproblem.co.nz)

    properties (SetAccess = private)
        prob        % Problem Structure (optiprob)
        nlprob      % Nonlinear Problem Structure
        opts        % Setup & Runtime options (optiset)
        sol         % Previous solution (Vector)
        obj         % Previous objective value
        ef          % Previous Exit Flag
        info        % Previous Info Structure
        fitstat     % Previous Fit Statistics Structure
    end

    methods

        %-- Constructor --%
        function op = opti(varargin)
            %OPTI Constructor
            if(~nargin)
                optiprob;
                fprintf('\n OPTION SETTING:\n');
                fprintf('         options: [ Problem Setup & Solving Options specified via optiset() ]\n');
                fprintf('          solver: [ Optimization Solver ]\n');
                return;
            end
            [p,o] = opti.buildOpti(varargin{:});
            [op.prob,op.opts,op.nlprob] = opti.buildConfig(p,o);
            op.sol = []; op.ef = [];
        end

        %-- Solve --%
        function [x,fval,exitflag,info] = solve(optObj,x0)
            %SOLVE  Solve OPTI Optimization Problem
            %
            %   [x,fval,exitflag,info] = solve(optObj) solves the optimization
            %   problem specified by optObj and returns the solution vector x,
            %   the function value at x in fval, an exitflag and information in
            %   info regarding the solution.
            %
            %   [x,fval,exitflag,info] = solve(optObj,x0) solves the optimization
            %   problem using x0 as the initial guess. If x0 is specified within
            %   opti() / optiprob(), this x0 will take precedence.
            %
            %   ExitFlags:
            %        1 - Converged / Terminated Successfully
            %        0 - Maximum Iterations / Function Evaluations / Time Exceeded
            %       -1 - Infeasible / Could Not Converge
            %       -2 - Unbounded / Solver Error
            %       -3 - Solver Specific Errors (Set option 'display' as 'iter')
            %       -5 - User Exited via Ctrl-C
            %
            %   Information Structure:
            %         Iterations - The number of iterations taken by the solver
            %               Time - The execution time of the solver as measured by MATLAB (tic + toc)
            %          Algorithm - The solver and algorithm being used
            %             Status - A status string indicating the solver specific exit message
            %             Lambda - A structure of the Lagrange multipliers at the solution (if available)

            %Check class isn't empty
            if(isempty(optObj.prob))
                error('You cannot solve an empty OPTI object!');
            end

            %Get Optional x0
            if(nargin < 2 || isempty(x0))
                if(~isempty(optObj.prob.x0) && ~any(isnan(optObj.prob.x0)))
                    x0 = optObj.prob.x0;
                else
                    x0 = [];
                end
            else
                optObj.prob.x0 = x0; %save into class
            end
            %Solve Problem
            [x,fval,exitflag,info] = solveOpti(optObj,x0);
            %Save solve info
            optObj.sol = x;
            optObj.obj = fval;
            optObj.ef = exitflag;
            optObj.info = info;

            %Update ndec if empty
            if(isempty(optObj.prob.sizes.ndec)), optObj.prob.sizes.ndec = length(x0); end
            %Check for maximization
            if(optObj.prob.sense==-1), fval = -fval; end

            %Check if we should write an ampl solution file
            if(optObj.prob.ampl.writesol && asl('isopen') == 1)
                asl('writesol',info.Status(1:min(end,100)),x);
                %If the problem is not nonlinear, close the interface
                if(isempty(optObj.prob.fun) && isempty(optObj.prob.nlcon))
                    asl('close');
                end
            end

            %If less that 3 output arguments, user may not know if not
            %solved correctly. Depending on warning level, show a warning
            %to the user
            if(nargout < 3 && ~isempty(exitflag) && exitflag ~= 1 && optiWarnLevel(optObj.opts.warnings) >= 1 && isfield(info,'Status') && strcmpi(optObj.opts.display,'off'))
                optiwarn(['Your problem may not have solved correctly. Solver status:\n>> %s\n\nPlease enable solver output by changing the optiset setting ''display'' to ''iter'' to see more information.\n',...
                          'Alternatively you can avoid this warning by returning all four return arguments from solve(opti_obj).\n'],info.Status);
            end
        end

        %-- Multi-Solve --%
        function [x,fval,exitflag,info] = multisolve(optObj,x0,divisions,penalty,solveAll)
            %MULTISOLVE  Solve OPTI Global Optimization Problem with Multi-Start (beta)
            %
            %   x = multisolve(optObj) performs a search of the
            %   optimization problem before attempting to solve the problem
            %   at the best solution guesses found. The best optimized
            %   solution will be returned in x.
            %
            %   x = multisolve(optObj,x0) also includes a more detailed
            %   search around the initial guess x0.
            %
            %   x = multisolve(optObj,x0,divisions) determines the number
            %   of divisions each variable is searched within. The default
            %   will attempt to search around 100 points. For ND problems,
            %   a division of 10 will result in 10^N search points.
            %
            %   x = multisolve(optObj,x0,divisions,penalty) specifies the
            %   penalty factor applied to constraint violations. The
            %   default is 1e4. Larger values will be required for
            %   objective functions that exceed this value.
            %
            %   x = multisolve(optObj,x0,divisions,penalty,solveAll)
            %   specifies at each search point an optimization is also
            %   performed. Significantly longer to execute, but ensures the
            %   search space is well covered!

            %Check class isn't empty
            if(isempty(optObj.prob))
                error('You cannot solve an empty OPTI object!');
            end

            %Optional args
            if(nargin < 5), solveAll = false; end
            if(nargin < 4), penalty = []; end
            if(nargin < 3), divisions = []; end

            %Get Optional x0
            if(nargin < 2 || isempty(x0))
                if(~isempty(optObj.prob.x0) && ~any(isnan(optObj.prob.x0)))
                    x0 = optObj.prob.x0;
                else
                    x0 = [];
                end
            else
                optObj.prob.x0 = x0; %save into class
            end

            %Solve Problem
            [x,fval,exitflag,info] = multiSolveOpti(optObj,x0,divisions,penalty,solveAll);
            %Save solve info into Object
            optObj.sol = x;
            optObj.obj = fval;
            optObj.ef = exitflag;
            optObj.info = info;

            %Update ndec if empty
            if(isempty(optObj.prob.sizes.ndec)), optObj.prob.sizes.ndec = length(x0); end
            %Check for maximization
            if(optObj.prob.sense==-1), fval = -fval; end

            %Check if we should write an ampl solution file
            if(optObj.prob.ampl.writesol && asl('isopen') == 1)
                asl('writesol',info.Status(1:min(end,100)),x);
                %If the problem is not nonlinear, close the interface
                if(isempty(optObj.prob.fun) && isempty(optObj.prob.nlcon))
                    asl('close');
                end
            end
        end

        %-- Check Solution --%
        function [ok,msg] = checkSol(optObj,tol)
            %CHECKSOL  Check Optimization Solution for Errors
            %
            %   [ok,msg] = checkSol(optObj) checks the solution stored in
            %   optObj for error messages or broken constraints. It returns
            %   true in 'ok' if no problems detected, otherwise it returns
            %   false with the problem in 'msg'.
            %
            %   [ok,msg] = checkSol(optObj,tol) checks with respect to a
            %   given tolerance. The default is 1e-6.

            if(isempty(optObj.sol))
                error('This OPTI object has not been solved yet!');
            end
            if(nargin < 2), tol = 1e-6; end

            [ok,msg] = checkOptiSol(optObj,tol);
        end

        %-- Display --%
        function display(optObj)
            %Check if empty
            if(isempty(optObj.prob))
                disp(' ');
                disp('------------------------------------------------------');
                disp('Empty OPTI Object');
                disp('------------------------------------------------------');
            else
                %Otherwise normal method
                displayOPTI(optObj);
            end
        end

        %-- Plot --%
        function plot(optObj,scale,doLog,npts)
            %PLOT  Plot optimization problem (1D-5D only)
            %
            %   plot(optObj) plots the optimization field for the current
            %   OPTI object.
            %
            %   plot(optObj,scale) plots with a user defined zoom level.
            %   scale is defined as the range +- of the solution to be
            %   drawn OR if supplied as a vector controls the bounds on the
            %   plot [x1min x1max x2min x2max ... xNmin xNmax]
            %
            %   plot(obtObj,scale,dolog) plots the log of the objective
            %   function (NL only)
            %
            %   plot(obtObj,scale,dolog,npts) controls the number of points
            %   used for generating the objective and constraint contours

            %Check class isn't empty
            if(isempty(optObj.prob))
                error('You cannot plot an empty OPTI object!');
            end

            %Get Optional Arguments
            if(nargin < 4), npts = []; end
            if(nargin < 3), doLog = false; end
            if(nargin < 2), scale = []; end

            %Check 2D for problems other than least squares fits
            if(~any(strcmpi(optObj.prob.type,{'NLS','DNLS'})) && optObj.prob.sizes.ndec > 5)
                error('You can only plot problems with one - five decision variables!');
            end

            %Plot
            plotOptiProb(optObj.prob,optObj.opts,optObj.sol,optObj.fitstat,scale,doLog,npts,'normal');
        end

        %-- Multi-Plot--%
        function multiplot(optObj,doLog,npts)
            %MULTIPLOT  Plot optimization problem and overlay multi-start search area
            %
            %   multiplot(optObj) plots the optimization field for the current
            %   OPTI object and overlays the multi-start search area
            %
            %   multiplot(obtObj,dolog) plots the log of the objective
            %   function (NL only)
            %
            %   multiplot(obtObj,dolog,npts) controls the number of points
            %   used for generating the objective and constraint contours

            if(nargin < 3), npts = []; end
            if(nargin < 2), doLog = false; end

            %Check class isn't empty
            if(isempty(optObj.prob))
                error('You cannot plot an empty OPTI object!');
            end
            %Check 2D
            if(optObj.prob.sizes.ndec > 5)
                error('You can only multiplot problems with one - five decision variables!');
            end

            %Plot
            plotOptiProb(optObj.prob,optObj.opts,optObj.sol,optObj.fitstat,[],doLog,npts,'multi');
        end

        %-- Stats --%
        function [StatsStruct] = fitStats(optObj,limit,printStats)
            %CONFIDENCE  Calculate approximate ANOVA and Confidence Statistics for NLS/DNLS Fit
            %
            %   StatsStruct = fitStats(optObj) returns a structure of
            %   statistics about the fit, as detailed below:
            %       SSE - Sum Of Squared Errors
            %       Rsquare - R^2 Value
            %       AdjRSquare - Adjusted R^2 Value
            %       RMSE - Root Mean Squared Error
            %       DFE - Degrees of Freedom
            %       ContInt - Confidence Interval (theta +- ConfInt)
            %       ConfBnds - Confidence Bounds (ypred +- ConfBnds)
            %       Cov - System Covariance Matrix
            %       BndIdx - Logical Array where 1s indicate active bounds
            %       at the solution (and therefore dropped from Cov)
            %       Param - Parameter Statistics Structure with fields:
            %           StdError - Standard Error
            %           tStat - t Statistic
            %           pValues - p Values
            %       Model - Model Statistics Stucture with fields:
            %           Fstat - F Statistic
            %           pValue - p value
            %           Int - Whether Model has an Intercept or Not
            %       Conf - User specified Confidence limit
            %
            %   For a description of each statistic, consult the following
            %   site:
            %       http://www.ats.ucla.edu/stat/sas/output/reg.htm
            %
            %   StatsStruct = fitStats(optObj,limit) allows the user to
            %   specify the confidence limit (0.0 < limit < 1.0),
            %   e.g. 0.95 = 95%
            %
            %   fitStats(optObj,limit,printStats) optionally prints the
            %   statistics to the command window. Defaults if no output is
            %   selected.

            %Check class isn't empty or unsolved
            if(isempty(optObj.prob) && isempty(optObj.sol))
                error('You cannot calculate statistics of an empty OPTI object!');
            end
            if(isempty(optObj.sol))
                error('You cannot calculate statistics of an unsolved OPTI object!');
            end

            %Get Optional Arguments
            if(nargin < 3), printStats = false; end
            if(nargin < 2), limit = []; end

            %Calculate
            StatsStruct = calcStatistics(optObj,limit);
            %Save to Object
            optObj.fitstat = StatsStruct;
%             %Optionally Plot Ellipse
%             if(plotEllipse)
%                 if(length(optObj.sol) < 2 || length(optObj.sol) > 3)
%                     optiwarn('opti:confEllipse','Cannot plot confidence ellipse/ellipsoid - only 2/3 parameter problems are supported');
%                 else
%                     plotConfEllipse(optObj.sol,StatsStruct);
%                 end
%             end
            %If no output arguments, display fit statistics
            if(nargout < 1 || printStats)
               	optiPrintFitStats(optObj.fitstat,optObj.sol);
            end
        end

        %-- write --%
        function write(optObj,filename,type)
            %WRITE  Write an OPTI Problem to a MPS/QPS/LP/SDPA/SEDUMI/GAMS file
            %
            %   write(optObj,filename) writes the OPTI problem to a
            %   filetype specified by the extension of filename.
            %
            %   write(optObj,filename,type) writes the OPTI problem to a
            %   filetype specified by the type argument.
            %
            %   Valid File Types & Extensions:
            %       Linear & Quadratic Problems
            %           - MPS [.mps] (LP/MILP)
            %           - QPS [.qps] (QP/MIQP)
            %           - LP  [.lp]  (LP)
            %       Semidefinite Problems
            %           - SDPA, DAT [.dat] (Dense)
            %           - SDPA-S, DAT-S [.dat-s] (Sparse)
            %           - MAT [.mat] (SeDuMi Sparse)
            %       General
            %           - GMS [.gms] (GAMS)

            %Check class isn't empty
            if(isempty(optObj.prob)), error('You cannot write an empty OPTI object!'); end
            if(nargin < 3), type = []; end
            %Check if we have a type
            if(~isempty(type))
                switch(lower(type))
                    case {'sdpa','sdpa-s','sedumi','sdpas','sdp','dat','dat-s','mat'}
                        sdpWrite(optObj.prob,filename,type);
                        return;
                    case {'mps','qps','lp'}
                        coinWrite(optObj.prob,filename,type);
                        return;
                    case 'gms'
                        gamsWrite(optObj.prob,filename);
                    otherwise
                        error('Unknown file type to write to!');
                end
            end
            %If not, go on extension
            if(isempty(strfind(filename,'.')))
                error('You must supply a file extension if you do not specify a file type!');
            end
            ext = regexp(filename,'\.','split');
            switch(lower(ext{end}))
                case {'dat','mat','dat-s'}
                    sdpWrite(optObj.prob,filename,type);
                case {'mps','qps','lp'}
                    coinWrite(optObj.prob,filename,type);
                case 'gms'
                    gamsWrite(optObj.prob,filename);
                otherwise
                    error('Unknown file extension to write to! For non-standard file extensions please specify the file type.');
            end
        end

        %-- getProb --%
        function retprob = getProb(optObj)
            %GETPROB  Return optiprob compatible problem structure
            %
            %   prob = getProb(optObj)

            %Check class isn't empty
            if(isempty(optObj.prob))
                error('You cannot get the problem from an empty OPTI object!');
            end

            retprob = optObj.prob;
            retprob = rmfield(retprob,{'iscon','sizes','numdif','type','ampl','save'});
            retprob.int = retprob.int.str;
        end

        %-- Set Solution (unit testing only use) --%
        function setSolution(optObj, fval, sol)
            %Check class isn't empty
            if(isempty(optObj.prob))
                error('You cannot get the problem from an empty OPTI object!');
            end
            optObj.obj = fval;
            if (numel(sol) ~= optObj.prob.sizes.ndec)
                error('Solution vector is the wrong size');
            else
                optObj.sol = sol(:);
            end
            optObj.info = struct('Status','HARDCODED SOLUTION','Algorithm','USER SPECIFIED');
        end
    end

    methods (Static)
        [prob,opts] = buildOpti(varargin);
        [prob,opts,nlprob] = buildConfig(prob,opts);
    end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function stepplot

rr=findobj(gcf,'Tag','edit1');
r1=get(rr,'string');% get the value of numenator of G(s)
rr=findobj(gcf,'Tag','edit2');
r2=get(rr,'string');% get the value of denomenator of G(s

numg=str2num(r1); %numenator of G(s)
deng=str2num(r2); %denomenator of G(s)

rr=findobj(gcf,'Tag','edit3');
r4=get(rr,'string');% get the value of numenator of H(s)
rr=findobj(gcf,'Tag','edit4');
r5=get(rr,'string');% get the value of denomenator of H(s)

numh=str2num(r4);%numenator of H(s)
denh=str2num(r5);%denomenator of H(s)
h=tf(numh,denh); % H(s)

rr=findobj(gcf,'Tag','slider1');
r6=get(rr,'value');% get the value of proportional controller
rr=findobj(gcf,'Tag','slider2');
r7=get(rr,'value');% get the value of integrater controller
rr=findobj(gcf,'Tag','slider3');
r8=get(rr,'value');% get the value of defferentiatr controller

Kp=r6;
Ki=r7;
Kd=r8;

numcf=[Kd Kp Ki];         % numinator of PID
dencf=[1 0];              % denomenator of PID
numf=conv(numcf,numg);
denf=conv(dencf,deng);
fsys=tf(numf,denf); % forward transfare function.
csys=feedback(fsys,h);
[y,t] = step(csys); % step response
plot(t,y)

author: Sable

repo: mcbench-benchmarks


text:
function kap=kappa0(x,y,varargin)

% Compute the constants for `tube-formula' based simultaneous
% confidence bands.
%
% Works for regression models only. Density estimation problems
% should be converted to counts, and use poisson regression
% 'family','poisson'.
%
% Essentially, this is a front-end to locfit, and so all optional
% arguments to locfit (eg, smoothing parameters) can be provided.
%
% To compute (or plot) the confidence bands, provide the output
% of the kappa0() function as the 'kappa' argument to a
% predict() or lfband() call.
%
%
% Example:
%
% load ethanol;
% fit = locfit(E,NOx,'alpha',0.5)
% kap = kappa0(E,NOx,'alpha',0.5)  % give same arguments!
% lfplot(fit)
% lfband(fit,'kappa',kap)     % plot the simultaneous bands
% z = predict(fit,[0.6 0.7 0.8]','kappa',kap,'band','g')
% z{3}                        % evaluate the bands.

fit = locfit(x,y,'module','kappa','ev','grid','mg',20,varargin{:});
z = fit.fit_points.kappa;
d = size(fit.data.x,2);
kap = z(1:(d+1));

return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [L,e] = spm_mci_like_ind (P,R,M,U,Y)
% Compute likelihood wrt selected time points
% FORMAT [L,e] = spm_mci_like_ind (P,R,M,U,Y)
%
% P         Flow parameters
% R         Initial state parameters
% M         Model structure
% U         Inputs  [Nin x N]
% Y         data
%
% L         Log likelihood
% e         Prediction errors
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny and Biswa Sengupta
% $Id: spm_mci_like_ind.m 6697 2016-01-27 14:57:28Z spm $

% Read data points and time indices
try, ind=Y.ind; catch, ind=1:M.N; end
Nt=length(ind);
y=Y.y;

M.x0=R; % Initial conditions
[G,sy,st] = spm_mci_fwd (P,M,U);

if st==-1, disp('Problem !'); return; end

% Prediction errors
g=G(ind,:);
e=Y.y-g;

% Log Likelihood
L = -0.5*trace(M.iCe*e'*e) + M.logdet_Ce - 0.5*Nt*log(2*pi);

author: spm


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
    = pamr_2_run(fid, data, epsilon, C, tc, opts)
% This program simulates the PAMR algorithm
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%    = pamr_run(fid, data, K, L, c, tc, opts)
%
% cum_ret: a number representing the final cumulative wealth.
% cumprod_ret: cumulative return until each trading period
% daily_ret: individual returns for each trading period
% daily_portfolio: individual portfolio for each trading period
%
% data: market sequence vectors
% fid: handle for write log file
% epsilon: mean reversion threshold
% C: aggressive parameter
% tc: transaction cost rate parameter
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%            = pamr_run(fid, data, K, L, c, tc, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n, m] = size(data);

% Variables for return, start with uniform weight
cum_ret = 1;
cumprod_ret = ones(n, 1);
daily_ret = ones(n, 1);
ret_m = ones(n, 2);  % total_ret + day_ret
day_weight = ones(m, 1)/m;  %#ok<*NASGU>
day_weight_o = zeros(m, 1);
daily_portfolio = zeros(n, m);

% print file head
fprintf(fid, '-------------------------------------\n');
fprintf(fid, 'Parameters [epsilon:%f, C:%f, tc:%f]\n', epsilon, C, tc);
fprintf(fid, 'day\t Daily Return\t Total return\n');

fprintf(1, '-------------------------------------\n');
if(~opts.quiet_mode)
    fprintf(1, 'Parameters [epsilon:%f, C:%f, tc:%f]\n', epsilon, C, tc);
    fprintf(1, 'day\t Daily Return\t Total return\n');
end

if (opts.progress)
	progress = waitbar(0,'Executing Algorithm...');
end
for t = 1:1:n,
    % Calculate t's portfolio at the beginning of t-th trading day
    if (t >= 2)
        [day_weight] = pamr_kernel(data(1:t-1, :), day_weight, eta);
    end

    % Normalize the constraint, always useless
    day_weight = day_weight./sum(day_weight);
    daily_portfolio(t, :) = day_weight';

    % Cal t's return and total return
    daily_ret(t, 1) = (data(t, :)*day_weight)*(1-tc/2*sum(abs(day_weight-day_weight_o)));
    cum_ret = cum_ret * daily_ret(t, 1);
    cumprod_ret(t, 1) = cum_ret;

    % Adjust weight(t, :) for the transaction cost issue
    day_weight_o = day_weight.*data(t, :)'/daily_ret(t, 1);

%     denominator = data(t, :)*data(t, :)'/(day_ret(t, 1)^2)-1/nStocks*sum(data(t, :))^2/(day_ret(t, 1)^2)+0.5/C_t;
    denominator = (data(t, :)-1/m*sum(data(t, :)))*(data(t, :)-1/m*sum(data(t, :)))'+0.5/C;
    if (~eq(denominator, 0.0))
        eta = (daily_ret(t, 1)-epsilon)/denominator;
    end
    eta = max(0, eta); % MRPA

    % Debug information
    % Time consuming part, other way?
    fprintf(fid, '%d\t%f\t%f\n', t, daily_ret(t, 1), cumprod_ret(t, 1));
    if (~opts.quiet_mode)
        if (~mod(t, opts.display_interval)),
            fprintf(1, '%d\t%f\t%f\n', t, daily_ret(t, 1), cumprod_ret(t, 1));
        end
    end
    if (opts.progress)
		if mod(t, 50) == 0
			waitbar((t/n));
		end
	end
end

% Debug Information
fprintf(fid, 'PAMR-2(epsilon:%.2f, C:%d, tc:%.4f), Final return: %.2f\n', ...
    epsilon, C, tc, cum_ret);
fprintf(fid, '-------------------------------------\n');
fprintf(1, 'PAMR-2(epsilon:%.2f, C:%d, tc:%.4f), Final return: %.2f\n', ...
    epsilon, C, tc, cum_ret);
fprintf(1, '-------------------------------------\n');
    if (opts.progress)
		close(progress);
	end

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function M = moveLRRVC(F, d)
%------------------------------------------------------------------------------
% Moves gridfunction F in horizontal direction.
% Excess area is filled by Vertex-Centered (VC) Reflection across boundaries.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: June 23, 2000.
% (c) 1999-2002 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
M = moveUDRVC(F.', d).';
%------------------------------------------------------------------------------


text:
%MDL_MICO Create model of Kinova Mico manipulator
%
% MDL_MICO is a script that creates the workspace variable mico which
% describes the kinematic characteristics of a Kinova Mico manipulator
% using standard DH conventions.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%   qr         vertical 'READY' configuration
%
% Reference::
% - "DH Parameters of Mico" Version 1.0.1, August 05, 2013.
%   Kinova
%
% Notes::
% - SI units of metres are used.
% - Unlike most other mdl_xxx scripts this one is actually a function that
%   behaves like a script and writes to the global workspace.
%
% See also Revolute, mdl_jaco, mdl_puma560, mdl_twolink, SerialLink.

% MODEL: Kinova, Mico, 6DOF, standard_DH

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

function r = mdl_mico()

    deg = pi/180;

    % robot length values (metres)  page 4
    D1 = 0.2755;
    D2 = 0.2900;
    D3 = 0.1233;
    D4 = 0.0741;
    D5 = 0.0741;
    D6 = 0.1600;
    e2 = 0.0070;

    % alternate parameters
    aa = 30*deg;
    ca = cos(aa);
    sa = sin(aa);
    c2a = cos(2*aa);
    s2a = sin(2*aa);
    d4b = D3 + sa/s2a*D4;
    d5b = sa/s2a*D4 + sa/s2a*D5;
    d6b = sa/s2a*D5 + D6;

    % and build a serial link manipulator

    % offsets from the table on page 4, "Mico" angles are the passed joint
    % angles.  "DH Algo" are the result after adding the joint angle offset.

    robot = SerialLink([
        Revolute('alpha', pi/2,  'a', 0,  'd', D1,   'flip')
        Revolute('alpha', pi,    'a', D2, 'd', 0,    'offset', -pi/2)
        Revolute('alpha', pi/2,  'a', 0,  'd', -e2,  'offset', pi/2)
        Revolute('alpha', 2*aa,  'a', 0,  'd', -d4b)
        Revolute('alpha', 2*aa,  'a', 0,  'd', -d5b, 'offset', -pi)
        Revolute('alpha', pi,    'a', 0,  'd', -d6b, 'offset', pi/2)
        ], ...
        'name', 'Mico', 'manufacturer', 'Kinova');

    %{
        % MDH version, no test  yet
    robot = SerialLink([
        Revolute('alpha', 0,     'a', 0,  'd', D1,   'modified', 'flip')
        Revolute('alpha', -pi/2, 'a', 0,  'd', 0,    'modified', 'offset', -pi/2)
        Revolute('alpha', 0,     'a', D2, 'd', e2,  'modified',  'offset', pi/2)
        Revolute('alpha', -pi/2, 'a', 0,  'd', d4b, 'modified')
        Revolute('alpha', 2*aa,  'a', 0,  'd', d5b, 'modified',  'offset', -pi)
        Revolute('alpha', 2*aa,  'a', 0,  'd', d6b, 'modified',  'offset', pi/2)
        ], ...
        'name', 'Mico', 'manufacturer', 'Kinova');
    %}

    % place the variables into the global workspace
    if nargin == 1
        r = robot;
    elseif nargin == 0
        assignin('caller', 'mico', robot);
        assignin('caller', 'qz', [0 0 0 0 0 0]); % zero angles, arm up
        assignin('caller', 'qr', [270 180 180 0 0 180]*deg); % vertical pose as per Fig 2
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function ED = edl_readexpdefs(filename)
%EDL_READEXPDEFS Reads in an experiment definition from XML file
%
% $ Syntax $
%   - ED = edl_readexpdefs(filename)
%
% $ Arguments $
%   - filename:     the filename of the experiment definition XML
%   - ED:           the read experiment definition struct
%
% $ Description $
%   - ED = edl_readexpdefs(filename) reads in an experiment definition
%     struct from an XML file. The format of XML can be referred to
%     the edl.spec.txt. The struct ED has following fields:
%
%       - name: the experiment definition name
%       - selfpath: the absolute path of the experiment definition self
%       - envconf: the filename of environment configuration
%       - env:     the environment configuration struct
%       - scriptdir: the script directory
%       - reportdir: the report directory
%       - mfiledir:  the m-files directory
%       - logfile:   the path of log file
%       - logger:    the logger
%
%       - variables: a struct of all variables,
%           using variable names as field names
%           using variable values as field values
%       - scripts: the struct of scripts
%           - using name as fieldnames, each field is a struct with
%               - func:     the scripting function
%               - path:     the path of the script file
%               - ctrlpath: the path of the control file
%               - params:   the other parameters
%               - refs:     the struct array of referenced scripts
%                   using role as field name
%                   using script name as field value
%               - refreps:  the struct array of referenced reports
%                   using role as field name
%                   using report name as field value
%               - refpaths: the struct array of referenced paths
%                   using role as field name
%                   using xml path of (script or report) as field value
%               - refscopes: the struct array of scope names of references
%                            (scripts | reports)
%                   using role as field name
%                   using scope names as field values
%       - experiments: the struct of experiments
%           - using name as fieldnames, each field is a struct with
%               - func:         the experiment function
%               - script:       the script name
%               - scriptpath:   the script filepath
%               - ctrlpath:     the control filepath
%       - reports: the struct array of reports
%           similar to scripts, except for that there is no ctrlpath
%
% $ History $
%   - Created by Dahua Lin, on Aug 12nd, 2006
%

%% Read file

xdoc = xmlread(filename);
edparent = slfilepart(filename, 'parent');
docelem = xdoc.getDocumentElement;

%% Read Header

headerfns = { ...
    'name', ...
    'envconf', ...
    'scriptdir', ...
    'reportdir', ...
    'mfiledir', ...
    'logfile'};
ED = cell2struct(cell(1, length(headerfns)), headerfns, 2);
headers = xml_getattribs(docelem);
ED = slparseprops(ED, headers);

if slisabspath(filename)
    ED.selfpath = filename;
else
    ED.selfpath = sladdpath(filename, cd());
end

if isempty(ED.envconf)
    ED.env = [];
else
    envpath = sladdpath(ED.envconf, edparent);
    ED.env = edl_readenvvars(envpath);
end

if isempty(ED.env) || ~isfield(ED.env, 'envname')
    error('edl:parseerror', ...
        'The environment must has variable named %s', 'envname');
end
envname = ED.env.envname;

ED.logger = sllog('rootpath', edparent);
if ~isempty(ED.logfile)
    ED.logfile = [ED.logfile, '.', envname, '.log'];
    ED.logger = addfiles(ED.logger, ED.logfile);
end

%% First-Pass Parsing
% 1. Build variable table progressively
% 2. Build basic structures
% 3. Translate variables

ED.variables = [];
ED.scripts = [];
ED.experiments = [];
ED.reports = [];

ED = first_pass(ED, docelem, struct('workdir', ''));

%% Second-Pass Parsing
% 1. Extract paths from cross-references
%   (a) build refpaths for scripts and reports
%   (b) build scriptpath and ctrlpath for experiments
%

% 1.(a)
ED.scripts = build_refpaths(ED, ED.scripts);
ED.reports = build_refpaths(ED, ED.reports);

% 1.(b)
ED.experiments = build_exppaths(ED, ED.experiments);

%% Core function for First-Pass (recursive invoking)

function ED = first_pass(ED, xelem, groupvars)

nodeList = xelem.getChildNodes;
n = nodeList.getLength;

for i = 1 : n   % enumerate all children

    node = nodeList.item(i-1);
    if (node.getNodeType == node.ELEMENT_NODE)
        tag = char(node.getTagName);

        switch tag

            case 'Var'
                [varname, varval] = parse_var(ED, node);
                ED.variables.(varname) = varval;

            case 'Script'
                [scriptname, scriptstruct] = parse_script(ED, node);
                scriptstruct.params = weak_update(scriptstruct.params, groupvars);
                ED.scripts.(scriptname) = scriptstruct;

            case 'Report'
                [reportname, reportstruct] = parse_report(ED, node);
                ED.reports.(reportname) = reportstruct;

            case 'Experiment'
                [expname, expstruct] = parse_experiment(ED, node);
                ED.experiments.(expname) = expstruct;

            case 'Group'
                subgroupvars = parse_group(ED, groupvars, node);
                ED = first_pass(ED, node, subgroupvars);

            otherwise
                error('edl:parseerror', ...
                    'Invalid element with tag name %s', tag);
        end

    end

end

%% Core functions for Second-Pass

function S = build_refpaths(ED, S)

fns = fieldnames(S);
n = length(fns);

for i = 1 : n
    fn = fns{i};

    S.(fn) = build_refpaths_fortype(ED, fn, S.(fn), 'scripts', 'refs');
    S.(fn) = build_refpaths_fortype(ED, fn, S.(fn), 'reports', 'refreps');

end

function S = build_refpaths_fortype(ED, name, S, typeset, typefield)

if ~isfield(S, 'refpaths')
    S.refpaths = [];
end
if ~isfield(S, 'refscopes')
    S.refscopes = [];
end

curset = S.(typefield);
if isempty(curset)
    return;
end

fns = fieldnames(curset);
n = length(fns);
pool = ED.(typeset);

for i = 1 : n
    role = fns{i};
    refname = curset.(role);
    if ~isfield(pool, refname)
        error('edl:parseerror', ...
            'The reference with role %s (name = %s) in %s is not found in pool %s', ...
            role, refname, name, typeset);
    end

    S.refpaths.(role) = pool.(refname).path;
    S.refscopes.(role) = typeset;
end

function ES = build_exppaths(ED, ES)

if isempty(ES)
    return;
end

fns = fieldnames(ES);
n = length(fns);

for i = 1 : n
    fn = fns{i};
    scriptname = ES.(fn).script;
    if ~isfield(ED.scripts, scriptname)
        error('edl:parseerror', ...
            'The referred script named %s for %s is not found', ...
            scriptname, fn);
    end

    curscript = ED.scripts.(scriptname);

    ES.(fn).scriptpath = curscript.path;
    ES.(fn).ctrlpath = curscript.ctrlpath;
end

%% Node Parsing functions

% variable parsing

function [vname, vval] = parse_var(ED, xelem)

vname = make_name(ED, 'variables', xelem);
vval = get_tattrib(ED, xelem, 'variables', vname, 'val');

% script parsing

function [sname, scr] = parse_script(ED, xelem)

sname = make_name(ED, 'scripts', xelem);

scr.func = get_tattrib(ED, xelem, 'scripts', sname, 'func');
scr.func = sladdpath(scr.func, ED.mfiledir);

scr.path = '';
scr.ctrlpath = '';
scr = update_fields(ED, scr, xelem, {'path', 'ctrlpath'});
scr.params = get_tattribs(ED, xelem, ...
    'exclude', {'name', 'func', 'path', 'ctrlpath'});

if isempty(scr.path)
    scr.path = [sname, '.script.xml'];
end
if isempty(scr.ctrlpath)
    scrfiletitle = slfilepart(scr.path, 'title');
    scr.ctrlpath = slchangefilepart(scr.path, 'title', [scrfiletitle, '.control']);
end

scr.path = sladdpath(scr.path, ED.scriptdir);

scr = build_reftables(ED, scr, xelem);

% report parsing

function [rname, rep] = parse_report(ED, xelem)

rname = make_name(ED, 'reports', xelem);

rep.func = get_tattrib(ED, xelem, 'reports', rname, 'func');
rep.func = sladdpath(rep.func, ED.mfiledir);

rep.path = '';
rep = update_fields(ED, rep, xelem, {'path'});
rep.params = get_tattribs(ED, xelem, ...
    'exclude', {'name', 'func', 'path', 'ctrlpath'});

if isempty(rep.path)
    rep.path = [rname, '.report.xml'];
end
rep.path = sladdpath(rep.path, ED.reportdir);

rep = build_reftables(ED, rep, xelem);

% experiment parsing

function [expname, es] = parse_experiment(ED, xelem)

expname = make_name(ED, 'experiments', xelem);

es.func = get_tattrib(ED, xelem, 'experiments', expname, 'func');
es.func = sladdpath(es.func, ED.mfiledir);

es.script = '';
es = update_fields(ED, es, xelem, {'script'});

if isempty(es.script)
    es.script = expname;
end

% group parsing

function newvars = parse_group(ED, gvars, xelem)

newvars = get_tattribs(ED, xelem, 'exclude', {'title'});

if ~isfield(newvars, 'workdir') || isempty(newvars.workdir)
    newvars.workdir = gvars.workdir;
else
    newvars.workdir = sladdpath(newvars.workdir, gvars.workdir);
end

newvars = weak_update(newvars, gvars);

%% Auxiliary functions

function val = translate_attribval(ED, attrval)

val = attrval;
if length(val) > 1 && val(1) == '$'
    varname = val(2:end);
    if ~isfield(ED.variables, varname)
        error('edl:parseerror', ...
            'The variable %s is not found', varname);
    end
    val = ED.variables.(varname);
end

function aval = get_tattrib(ED, xelem, scope, elemname, attrname)

if xelem.hasAttribute(attrname)
    aval = char(xelem.getAttribute(attrname));
    aval = translate_attribval(ED, aval);
else
    error('edl:parseerror', ...
        'The attribute %s is not found in %s of %s', attrname, elemname, scope);
end

function A = get_tattribs(ED, xelem, varargin)

A = xml_getattribs(xelem, varargin{:});

if ~isempty(A)
    fns = fieldnames(A);
    n = length(fns);
    for i = 1 : n
        fn = fns{i};
        A.(fn) = translate_attribval(ED, A.(fn));
    end
end

function name = make_name(ED, scope, xelem)

if xelem.hasAttribute('name')
    name = char(xelem.getAttribute('name'));
    if isfield(ED.(scope), name)
        error('edl:parseerror', ...
            'Redefinition of %s in %s', name, scope);
    end
else
    error('edl:parseerror', ...
        'Encounter an element without name in %s', scope);
end

function S = update_fields(ED, S, xelem, fns)

S = slparseprops(S, get_tattribs(ED, xelem, ...
    'select', fns, ...
    'forceexist', false));

function S = weak_update(S, Scomp)

fns = fieldnames(Scomp);
n = length(fns);
for i = 1 : n
    fn = fns{i};
    if ~isfield(S, fn)
        S.(fn) = Scomp.(fn);
    end
end

function S = build_reftables(ED, S, xelem)

refNodes = xelem.getElementsByTagName('Ref');
nrs = refNodes.getLength;

S.refs = [];
S.refreps = [];

for i = 1 : nrs
    curnode = refNodes.item(i-1);
    cur = get_tattribs(ED, curnode, 'select', {'role', 'name'});
    S.refs.(cur.role) = cur.name;
end

refrepNodes = xelem.getElementsByTagName('RefReport');
nrr = refrepNodes.getLength;

for i = 1 : nrr
    curnode = refrepNodes.item(i-1);
    cur = get_tattribs(ED, curnode, 'select', {'role', 'name'});
    S.refreps.(cur.role) = cur.name;
end

author: lmthang


text:
function [model, lbls] = fgplvmLoadResult(dataSet, number)

% FGPLVMLOADRESULT Load a previously saved result.
% FORMAT
% DESC loads a previously saved FGPLVM result.
% ARG dataSet : the name of the data set to load.
% ARG number : the number of the FGPLVM run to load.
% RETURN model : the saved model.
% RETURN lbls : labels of the data set (for visualisation purposes).
%
% SEEALSO : fgplvmLoadResult
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2004, 2005, 2006, 2008

% FGPLVM

[Y, lbls] = lvmLoadData(dataSet);

dataSet(1) = upper(dataSet(1));
load(['dem' dataSet 'Fgplvm' num2str(number)])
model = fgplvmReconstruct(kern, noise, fgplvmInfo, X, Y);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [view,hdr] = BV_voi2roi(view,voiFile,unTalFlag,talPath);
% BV_voi2roi: converts a brain voyager .voi file into
% an ROI in a mrLoadRet volume.
%
% Usage: [view,hdr] = BV_voi2roi(view,[voiFile,unTalFlag,talPath]);
%
% Args:
%       voiFile: path of Brain Voyager .voi file. If not passed in,
%       provides a dialog to select it.
%
%       unTalFlag: flag (set to 0 or 1) which causes the VOI/ROIs to be
%       converted from Talairach coords (saved that way in BV) to unTalled
%       coords (in mrLoadRet). Defaults to 1, do untal. NOTE: even if
%       manually set to 0, if the BV .voi file specifies the file type as
%       'TAL', it will automatically set unTalFlag to 1.
%
%       talPath: for when unTalFlag==1, the location of a .mat file,
%       created using computeTalairach.m, which contains the transform
%       information for going between tal and untal space. If unspecified,
%       checks the vANATOMYPATH directory, and, failing that, prompts the
%       user.
%
% Dependencies: BV_tal2untal
%
% ras 03/03
mrGlobals;

% right now I'm just making this work with volume views --
% bow out if it's another type of view
if ~isequal(view.viewType,'Volume') & ~isequal(view.viewType,'Gray')
    fprintf('Sorry, but this only works with Volume/Gray views right now.\n');
    return
end

if ~exist('unTalFlag','var')    unTalFlag = 1;      end

% prompt for BV .voi file if it's not passed in
if ~exist('voiFile','var') | isempty(voiFile)
    msg = ['Choose a Brain Voyager .voi file...'];
    startDir = fullfile(RAID,'3Danat','BV_vmr');
    if ~exist(startDir)
        startDir = pwd;
    end
    voiFile = getPathStrDialog(startDir,msg,'*.voi');
end

if ~exist(voiFile,'file')
    error(['File ' voiFile ' is not a valid file.']);
end

% prompt for a BV .tal file containing talairach information
if unTalFlag==1 & ~exist('talPath','var')

    % make an educated guess as to where the talPath might be
    guess = fullfile(fileparts(vANATOMYPATH),'vAnatomy_talairach.mat');

    if exist(guess,'file')
        talPath = guess;
    else
        disp('Choose a .mat file containing the talairach points for this volume');
        disp('(See ''help computeTalairach'' if you are unsure of what this is.)');
        msg = ['Choose a .mat file from coputeTalairach...'];
        talPath = getPathStrDialog(pwd,msg,'*.mat');
    end
 end

% Get the volume size - for flipping one of the dimensions to preserve L/R
volAnatSize=size(view.anat);

% open the voi file
[fid message] = fopen(voiFile,'r');
if fid==-1
    error(message);
end

%%%%% read the header
blank = fgetl(fid);
ignoretxt = fscanf(fid,'%s',1);
hdr.FileVersion = fscanf(fid,'%i',1);
blank = fgetl(fid);
blank = fgetl(fid);
ignoretxt = fscanf(fid,'%s',1);
hdr.CoordsType = fscanf(fid,'%s',1);
if isequal(hdr.CoordsType,'TAL')
    unTalFlag = 1;
end
blank = fgetl(fid);
blank = fgetl(fid);
blank = fgetl(fid);
ignoretxt = fscanf(fid,'%s',1);
hdr.numVOIs = fscanf(fid,'%i',1);
blank = fgetl(fid);
blank = fgetl(fid);

%%%%% loop through ROIs, reading them in
for voi = 1:hdr.numVOIs

    % read VOI/ROI name, num voxels in brainvoyager coords
    ignoretxt = fscanf(fid,'%s',1);
    hdr.VOInames{voi} = fscanf(fid,'%s',1);
    blank = fgetl(fid);
    blank = fgetl(fid);
    ignoretxt = fscanf(fid,'%s',1);
    hdr.numBVVoxels(voi) = fscanf(fid,'%i',1);
    blank = fgetl(fid);

    % init ROI substruct
    ROIs(voi).color = 'b'; % default color is blue
    ROIs(voi).coords = [];
    ROIs(voi).name = hdr.VOInames{voi};
    ROIs(voi).viewType = view.viewType;

    % read coordinates of VOI
    voiCoords = [];
    for vxl = 1:hdr.numBVVoxels(voi)
        voiCoords(vxl,1) = fscanf(fid,'%i',1);
        voiCoords(vxl,2) = fscanf(fid,'%i',1);
        voiCoords(vxl,3) = fscanf(fid,'%i',1);
    end
    blank = fgetl(fid);
    blank = fgetl(fid);

	if unTalFlag
        voiCoords = BV_tal2unTal(voiCoords,talPath);
    end

    ROIs(voi).coords = round(voiCoords');

    fprintf('%s \n',ROIs(voi).name);
end

%%%%% read linked VTC files at the end; put info in hdr
blank = fgetl(fid);
ignoretxt = fscanf(fid,'%s',1);
hdr.numLinkedVTCs = fscanf(fid,'%i',1);
blank = fgetl(fid);
for vtc = 1:hdr.numLinkedVTCs
    hdr.VTCs{vtc} = fgetl(fid);
end

% close the file
fclose(fid);

% append ROIs to current view
view.ROIs = [view.ROIs ROIs];

% let the user know it worked
fprintf('Converted %i ROIs from the BV file %s.\n',length(ROIs),voiFile);

% refresh ROIs
r = view.selectedROI;
if r > 0
    view = selectROI(view,r);
else
	view = selectROI(view,1);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% triplet_find
% Below is a demonstration of the features of the |triplet_find| function

%%
clear; close all; clc;

%% Syntax
% |[L1t,L2t,L3t]=triplet_find(DT1,DT2,DT3,Vn1,Vn2,Vn3,T_dist,opt);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function rcnn_model = rcnn_create_model(cnn_definition_file, cnn_binary_file, cache_name)
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

if ~exist('cache_name', 'var') || isempty(cache_name)
  cache_name = 'none';
end

%  model =
%    cnn: [1x1 struct]
%        binary_file: 'path/to/cnn/model/binary'
%        definition_file: 'path/to/cnn/model/definition'
%        batch_size: 256
%        image_mean: [227x227x3 single]
%        init_key: -1
%    detectors.W: [N x <numclasses> single]  % matrix of SVM weights
%    detectors.B: [1 x <numclasses> single]  % (row) vector of SVM biases
%    detectors.crop_mode: 'warp' or 'square'
%    detectors.crop_padding: 16
%    detectors.nms_thresholds: [1x20 single]
%    training_opts: [1x1 struct]
%        bias_mult: 10
%        fine_tuned: 1
%        layer: 'fc7'
%        pos_loss_weight: 2
%        svm_C: 1.0000e-03
%        trainset: 'trainval'
%        use_flipped: 0
%        year: '2007'
%        feat_norm_mean: 20.1401
%    classes: {cell array of class names}
%    class_to_index: map from class name to column index in W

% init empty convnet
assert(exist(cnn_binary_file, 'file') ~= 0);
assert(exist(cnn_definition_file, 'file') ~= 0);
cnn.binary_file = cnn_binary_file;
cnn.definition_file = cnn_definition_file;
cnn.batch_size = 32;
cnn.init_key = -1;
cnn.input_size = 224;
% load the ilsvrc image mean
data_mean_file = './external/caffe/matlab/caffe/ilsvrc_2014_mean.mat';
%data_mean_file = './external/caffe/matlab/caffe/ilsvrc_2014_mean.mat';

assert(exist(data_mean_file, 'file') ~= 0);
ld = load(data_mean_file);
image_mean = ld.image_mean; clear ld;
off = floor((size(image_mean,1) - cnn.input_size)/2)+1;
image_mean = image_mean(off:off+cnn.input_size-1, off:off+cnn.input_size-1, :);
cnn.image_mean = image_mean;

% init empty detectors
detectors.W = [];
detectors.B = [];
detectors.crop_mode = 'warp';
detectors.crop_padding = 16;
detectors.nms_thresholds = [];

% rcnn model wraps the convnet and detectors
rcnn_model.cnn = cnn;
rcnn_model.cache_name = cache_name;
rcnn_model.detectors = detectors;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function ye = fR(t,D,nalpha,nbeta,u,bc,h)	% Right boundary value
alpha = nalpha*pi; beta = nbeta*pi;
if bc == 1
  ye = - h*beta*sin(beta*(1 - u*t))...
    - alpha*exp(-D*alpha*alpha*t)*sin(alpha*(1 - u*t));
else
  ye = 0.0;
end

repo: Aero-matlab


text:
function v = isglobalpath(path)

v = (~ispc && ~isempty(path) && path(1) == '/') || ...
  (ispc && ~isempty(regexp(path,'^[A-Z]:[\\/]','once')));

author: kristinbranson

repo: JAABA


text:
function value = get_option(option_list,option,default,type)
% extract option from option list
%
% In case the option occurs more than one time in the option list the first
% value that meets the type condition is returned.
%
% Syntax
%   value = get_option(option_list,option_name,default_option,option_type)
%
% Input
%  option_list   - Cell Array
%  option_name   - String
%  option_default- default value
%  option_type   - class name or list of class names
% Output
%  value         - option value
%
% See also
% check_option set_option delete_option

if nargin <= 2, value = [];else, value = default;end
if isempty(option_list), return;end
if nargin <= 3, type = [];end

pos = find_option(option_list,option,type);

if pos, value = option_list{pos};end

author: mtex-toolbox

repo: mtex

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


text:
function vl_compilenn_calvin()
% vl_compilenn_calvin()
%
% Compile the C code provided in Matconvnet-calvin.
% Matconvnet code needs to be compiled separately.
%
% Copyright by Holger Caesar, 2016

root = calvin_root();
codeDir = fullfile(root, 'matconvnet-calvin', 'matlab');
mexDir = fullfile(codeDir, 'mex');
mexOpts = {'-largeArrayDims', '-outdir', sprintf('"%s"', mexDir)};

% E2S2-related
mex(mexOpts{:}, fullfile(codeDir, 'labelpresence', 'labelPresence_backward.cpp'));
mex(mexOpts{:}, fullfile(codeDir, 'regiontopixel', 'regionToPixel_backward.cpp'));
mex(mexOpts{:}, fullfile(codeDir, 'regiontopixel', 'regionToPixelSoft_backward.cpp'));
mex(mexOpts{:}, fullfile(codeDir, 'roipool', 'roiPooling_forward.cpp'));
mex(mexOpts{:}, fullfile(codeDir, 'roipool', 'roiPooling_backward.cpp'));

% Misc
mex(mexOpts{:}, fullfile(codeDir, 'misc', 'computeBlobOverlapAnyPair.cpp'));
mex(mexOpts{:}, fullfile(codeDir, 'misc', 'scoreBlobIoUs.cpp'));

author: nightrome

repo: matconvnet-calvin


text:
function [M,E,EMAP] = crouzeix_raviart_massmatrix(V,F)
  % CROUZEIX_RAVIART_MASSMATRIX Compute the Crouzeix-Raviart mass matrix where
  % M(e,e) is just the sum of 1/3 the areas of the triangles on either side of
  % an edge e. For tets, edges are now faccets.
  %
  % See for example "Discrete Quadratic Curvature Energies" [Wardetzky, Bergou,
  % Harmon, Zorin, Grinspun 2007]
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by element-size list of triangle indices
  % Outputs:
  %   M  #E by #E edge-based diagonal mass matrix
  %   E  #E by 2 list of edges
  %
  % See also: edge_laplacian, is_boundary_edge, crouzeix_raviart_cotmatrix,
  %   massmatrix
  %

  switch size(F,2)
  case 3
    allE = [F(:,[2 3]);F(:,[3 1]);F(:,[1 2])];
    % Map duplicate edges to first instance
    [E,~,EMAP] = unique(sort(allE,2),'rows');
    TA = doublearea(V,F)/2;
    M = sparse(EMAP,EMAP,repmat(TA/3,3,1),size(E,1), size(E,1));
  case 4
    T = F;
    allF = [ ...
      T(:,2) T(:,4) T(:,3); ...
      T(:,1) T(:,3) T(:,4); ...
      T(:,1) T(:,4) T(:,2); ...
      T(:,1) T(:,2) T(:,3); ...
      ];
    vol = volume(V,T);
    [F,~,FMAP] = unique(sort(allF,2),'rows');
    M = sparse(FMAP,FMAP,repmat(vol/4,4,1),size(F,1), size(F,1));
    E = F;
    EMAP = FMAP;
  end
end

author: alecjacobson

repo: gptoolbox


text:
function value = p31_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P31_F evaluates the integrand for problem 31.
%
%  Dimension:
%
%    DIM_NUM arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integral Parameters:
%
%    There is a basis point Z associated with the integrand.
%    Z(1:DIM_NUM) defaults to ( 0.5, 0.5, ..., 0.5 ).
%    The user can set, get, or randomize this value by calling
%    P31_R8VEC.
%
%    The coefficient vector C (whose entries are usually positive)
%    controls the steepness and circularity of the pseudo-Gaussian.
%    C(1:DIM_NUM) defaults to 2.0.
%    The user can set, get, or randomize this value by calling
%    P31_R8VEC.
%
%  Integrand:
%
%    exp ( - sum ( c(1:dim_num) * abs ( x(1:dim_num) - z(1:dim_num) ) ) )
%
%  Exact Integral:
%
%    The integral is separable into
%
%      Int ( A(1) <= X(1) <= B(1) ) exp ( - C(1) * abs ( X(1) - Z(1) ) )
%        * Int ( A(2) <= X(2) <= B(2) ) exp ( - C(2) * abs ( X(2) - Z(2) ) )
%          * ...
%
%    Hence, the exact integral is computed as the product of
%    one dimensional integrals.  Each of these is easily computed
%    once the location of Z(I) with respect to A(I) and B(I) is
%    determined.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alan Genz,
%    [Integral #5]
%    A Package for Testing Multiple Integration Subroutines,
%    in Numerical Integration: Recent Developments, Software
%    and Applications,
%    edited by Patrick Keast and Graeme Fairweather,
%    D. Reidel, 1987, pages 337-340,
%    LC: QA299.3.N38.
%
%    Kenneth Hanson,
%    Quasi-Monte Carlo: halftoning in high dimensions?
%    in Computatinal Imaging,
%    Edited by CA Bouman and RL Stevenson,
%    Proceedings SPIE,
%    Volume 5016, 2003, pages 161-172.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  c = [];
  c = p31_r8vec ( 'G', 'C', dim_num, c );

  z = [];
  z = p31_r8vec ( 'G', 'Z', dim_num, z );

  value(1:point_num) = 0.0;

  for point = 1 : point_num

    value(point) = exp ( - sum ( c(1:dim_num)' .* ...
      abs ( x(1:dim_num,point) - z(1:dim_num)' ) ) );

  end

  p31_i4 ( 'I', '#', point_num );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function maxComplexity = getComplexity(filename)
%GETCOMPLEXITY  Calculate the maximum complexity in an M-file.
%   maxComplexity = getcomplexity(filename)
%
%   "Complexity" refers to the McCabe complexity as returned by "mlint -cyc"
%   When there are multiple functions in a file, the maximum value is
%   returned.

% The MATLAB Contest Team
% Copyright 2010-2011 The MathWorks, Inc.

msg = mlint('-cyc','-string',filename);
tk = regexp(msg,'McCabe complexity of [^\s]+ is (\d+)','tokens');
complexity = zeros(size(tk));
for j = 1:length(tk)
    complexity(j) = eval(tk{j}{1});
end
maxComplexity = max(complexity);


text:
function [yUpdate,PInvUpdate,innov]=cubInfoUpdate(yPred,PInvPred,z,RInv,h,xi,w,innovTrans,measAvgFun,stateDiffTrans)
%%CUBINFOUPDATE Perform the measurement update step in the cubature
%               information filter.
%
%INPUTS: yPred The xDimX1 predicted information state. The information
%              state is the inverse covariance matrix times the target
%              state.
%     PInvPred The xDimXxDim inverse of the predicted state covariance
%              matrix.
%            z The zDim X 1 vector measurement.
%         RInv The zDim X zDim inverse of the measurement covariance matrix
%              in the native coordinate system of the measurement.
%            h A function handle for the measurement function that takes
%              the state as its argument.
%           xi An xDim X numCubPoints matrix of cubature points. If this
%              and the next parameter are omitted or empty matrices are
%              passed, then fifthOrderCubPoints(xDim) is used. It is
%              suggested that xi and w be provided to avoid needless
%              recomputation of the cubature points.
%            w A numCubPoints X 1 vector of the weights associated with the
%              cubature points.
%   innovTrans An optional function handle that computes and optionally
%              transforms the value of the difference between the
%              observation and any predicted points. This is called as
%              innovTrans(a,b) and the default if omitted or an empty
%              matrix is passed is @(a,b)bsxfun(@minus,a,b). This must be
%              able to handle sets of values. For a zDimX1 measurement,
%              either of the inputs could be zDimXN in size while one of
%              the inputs could be zDimX1 in size.  This only needs to be
%              supplied when a measurement difference must be restricted
%              to a certain range. For example, the innovation between two
%              angles will be 2*pi if one angle is zero and the other
%              2*pi, even though they are the same direction. In such an
%              instance, a function handle to the
%              wrapRange(bsxfun(@minus,a,b),-pi,pi) function with the
%              appropriate parameters should be passed for innovTrans.
%   measAvgFun An optional function handle that, when given N measurement
%              values with weights, produces the weighted average. This
%              function only has to be provided if the domain of the
%              measurement is not linear. For example, when averaging
%              angular values, then the function meanAng should be used.
% stateDiffTrans An optional function handle that takes an xDimXN matrix of
%              N differences between states and transforms them however
%              might be necessary. For example, a state containing angular
%              components will generally need to be transformed so that the
%              difference between the angles is wrapped to -pi/pi.
%
%OUTPUTS: yUpdate The xDim X 1 updated (posterior) information state
%                 vector.
%      PInvUpdate The updated xDim X xDim inverse state covariance matrix.
%           innov The zDimX1 innovation. This is the difference between the
%                 measurement and the predicted measurement. This is
%                 sometimes used for gating measurements.
%
%If the function h needs additional parameters beyond the state, then the
%parameters can be passed by using an anonymous function as the function
%handle. For example, suppose that the measurement function is measFunc and
%it needs the additional parameters param1 and param2. In this instance,
%rather than using
%h=@measFunc
%one should use
%h=@(x)measFunc(x,param1,param2)
%This way, every time cubKalUpdate calls measFunc (via h) with a
%different x, those two parameters are always passed.
%
%This function is an implementation of the measurement update step of
%Algorithm 1 in [1].
%
%The optional parameters innovTrans and measAvgFun are not described in
%[1], but allow for possible modifications to the filter as
%described in [2]. The parameters have been added to allow the filter to be
%used with angular quantities. For example, if the measurement consisted of
%range and angle, z=[r;theta], then
%innovTrans=@(a,b)[bsxfun(@minus,a(1,:),b(1,:));
%                  wrapRange(bsxfun(@minus,a(2,:),b(2,:)),-pi,pi)];
%measAvgFun=@(z,w)[calcMixtureMoments(z(1,:),w);
%                  meanAng(z(2,:),w')];
%should be used to approximately deal with the circular nature of the
%measurements.
%
%REFERENCES:
%[1] K. P. B. Chandra, D.-W. Gu, and I. Postlethwaite, "Square root
%    cubature information filter," IEEE Sensors Journal, vol. 13, no. 2,
%    pp. 750-758, Feb. 2013.
%[2] D. F. Crouse, "Cubature/ unscented/ sigma point Kalman filtering with
%    angular measurement models," in Proceedings of the 18th International
%    Conference on Information Fusion, Washington, D.C., 6-9 Jul. 2015.
%
%October 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(yPred,1);

if(nargin<6||isempty(xi))
    [xi,w]=fifthOrderCubPoints(xDim);
end

if(nargin<8||isempty(innovTrans))
    %The function just returns the input.
    innovTrans=@(a,b)bsxfun(@minus,a,b);
end

if(nargin<9||isempty(measAvgFun))
   measAvgFun=@(zPoints,w)calcMixtureMoments(zPoints,w);
end

if(nargin<10||isempty(stateDiffTrans))
   stateDiffTrans=@(x)x;
end

zDim=size(z,1);

numCubPoints=size(xi,2);

%Extract the state
xPred=PInvPred\yPred;
%Get a square root of the inverse covariance matrix.
SPred=cholSemiDef(pinv(PInvPred),'lower');
%Predicted cubature state points
xPredPoints=bsxfun(@plus,SPred*xi,xPred);

%Predicted cubature measurement points
zPredPoints=zeros(zDim,numCubPoints);
for curP=1:numCubPoints
    zPredPoints(:,curP)=h(xPredPoints(:,curP));
end

%Measurement prediction.
zPred=measAvgFun(zPredPoints,w);

%The innovation, transformed as necessary to keep values in a desired
%range.
innov=innovTrans(z,zPred);

Pxz=zeros(xDim,zDim);
for curP=1:numCubPoints
    diff=innovTrans(zPredPoints(:,curP),zPred);
    Pxz=Pxz+w(curP)*stateDiffTrans(xPredPoints(:,curP)-xPred)*diff';
end

I=PInvPred*Pxz*RInv*Pxz'*PInvPred';
i=PInvPred*Pxz*RInv*(innov+Pxz'*PInvPred'*xPred);

yUpdate=yPred+i;
PInvUpdate=PInvPred+I;
%Ensure symmetry
PInvUpdate=(PInvUpdate+PInvUpdate')/2;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function edgeFaces = meshEdgeFaces(vertices, edges, faces) %#ok<INUSL>
%MESHEDGEFACES Compute index of faces adjacent to each edge of a mesh.
%
%   EF = meshEdgeFaces(V, E, F)
%   Compute index array of faces adjacent to each edge of a mesh.
%   V, E and F define the mesh: V is vertex array, E contains vertex
%   indices of edge extremities, and F contains vertex indices of each
%   face, either as a numerical array or as a cell array.
%   The result EF has as many rows as the number of edges, and two column.
%   The first column contains index of faces located on the left of the
%   corresponding edge, whereas the second column contains index of the
%   face located on the right. Some indices may be 0 if the mesh is not
%   'closed'.
%
%   Note: a faster version is available for triangular meshes.
%
%   Example
%   meshEdgeFaces
%
%   See also
%   meshes3d, trimeshEdgeFaces, meshDihedralAngles, polyhedronMeanBreadth

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2010-10-04, using Matlab 7.9.0.529 (R2009b)
% Copyright 2010-2022 INRA - Cepia Software Platform

Ne = size(edges, 1);

% indices of faces adjacent to each edge
edgeFaces = zeros(Ne, 2);

% different method for extracting current face depending if faces are
% stored as index 2D array or as cell array of 1D arrays.
if isnumeric(faces)
    Nf = size(faces, 1);
    for i = 1:Nf
        face = faces(i, :);
        processFace(face, i)
    end
elseif iscell(faces)
    Nf = length(faces);
    for i = 1:Nf
        face = faces{i};
        processFace(face, i)
    end
end

    function processFace(face, indFace)
        % iterate on face edges
        for j = 1:length(face)
            % build edge: array of vertices
            j2 = mod(j, length(face)) + 1;

            % do not process edges with same vertices
            if face(j) == face(j2)
                continue;
            end

            % vertex indices of current edge
            currentEdge = [face(j) face(j2)];

            % find index of current edge, assuming face is left-located
            b1 = ismember(edges, currentEdge, 'rows');
            indEdge = find(b1);
            if ~isempty(indEdge)
                if edgeFaces(indEdge, 1) ~= 0
                    error('meshes3d:IllegalTopology', ...
                        'Two faces were found on left side of edge %d ', indEdge);
                end

                edgeFaces(indEdge, 1) = indFace;
                continue;
            end

            % otherwise, assume the face is right-located
            b2 = ismember(edges, currentEdge([2 1]), 'rows');
            indEdge = find(b2);
            if ~isempty(indEdge)
                if edgeFaces(indEdge, 2) ~= 0
                    error('meshes3d:IllegalTopology', ...
                        'Two faces were found on left side of edge %d ', indEdge);
                end

                edgeFaces(indEdge, 2) = indFace;
                continue;
            end

            % If face was neither left nor right, error
            warning('meshes3d:IllegalTopology', ...
                'Edge %d of face %d was not found in edge array', ...
                j, indFace);
            continue;
        end
    end

end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function DoUpdateInfo(Simuh,Info)
    string=get(Simuh.Info_listbox,'String');
    time=clock;
    string{end+1,1}=[num2str(time(4)) ':' num2str(time(5)) ':' num2str(round(time(6))) '>' Info];
    set(Simuh.Info_listbox,'String',string,'Value',length(string));
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
%% set up path

[~,computername] = system('hostname');
computername = strtrim(computername);

switch computername,

  case 'bransonk-lw1',

    JCtrax_path = 'E:\Code\JCtrax';
    %FlyBowlAnalysis_path = 'E:\Code\FlyBowlAnalysis';
    %rootdatadir = 'E:\Code\Jdetect\larva\fly_data\TrainingData';
    %settingsdir = 'E:\Code\FlyBowlAnalysis\settings';
    configfilename = 'params\JLabelParams.xml';
    npool = 4;

  case 'bransonk-lw2',

    JCtrax_path = 'C:\Code\JCtrax';
    %FlyBowlAnalysis_path = 'C:\Code\FlyBowlAnalysis';
    %rootdatadir = 'C:\Code\Jdetect\larva\fly_data\TrainingData';
    %settingsdir = 'C:\Code\FlyBowlAnalysis\settings';
    configfilename = 'params\JLabelParamsKristinChase.xml';
    npool = 4;

  case 'bransonk-desktop',

    JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
    %FlyBowlAnalysis_path = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis';
    %rootdatadir = '/groups/branson/home/bransonk/behavioranalysis/code/Jdetect/Jdetect/larva/fly_data/TrainingData/';
    %settingsdir = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis/settings';
    configfilename = 'params/JLabelParams_bransonk-desktop.xml';
    npool = 8;

  case 'kabram-ws.janelia.priv',

    JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
    %FlyBowlAnalysis_path = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis';
    %rootdatadir = '/groups/branson/home/bransonk/behavioranalysis/code/Jdetect/Jdetect/larva/fly_data/TrainingData/';
    %settingsdir = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis/settings';
    configfilename = 'params/JLabelParamsMayank.xml';
    npool = 8;

  case 'robiea-ws'
    JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
    %FlyBowlAnalysis_path = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis';
    %rootdatadir = '/groups/sciserv/flyolympiad/Olympiad_Screen/fly_bowl/bowl_data';
    %settingsdir = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis/settings';
    configfilename = '/groups/branson/home/robiea/Projects_data/JLabel/JLabelParams.xml';
    npool = 8;

  otherwise

    warning('Unknown computer %s. Paths may not be setup correctly',computername);
    JCtrax_path = '/groups/branson/home/bransonk/tracking/code/JCtrax';
    %FlyBowlAnalysis_path = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis';
    %rootdatadir = '/groups/branson/home/bransonk/behavioranalysis/code/Jdetect/Jdetect/larva/fly_data/TrainingData/';
    %settingsdir = '/groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis/settings';
    configfilename = 'params/JLabelParams_bransonk-desktop.xml';
    npool = 1;

end

addpath(genpath(fullfile(JCtrax_path,'pdollar_toolbox')));
addpath(fullfile(JCtrax_path,'misc'));
addpath(fullfile(JCtrax_path,'filehandling'));
jlabelpath = fileparts(which('JLabel'));
addpath(fullfile(jlabelpath,'compute_perframe_features'));
%addpath(FlyBowlAnalysis_path);

%% parameters

rootdatadir = '/groups/branson/bransonlab/projects/olympiad/HackHitData';
perframe_params = struct;
perframe_params.fov = pi;
perframe_params.thetafil = [.0625,.25,.375,.25,.0635];
perframe_params.nbodylengths_near = 2.5;
perframe_params.max_dnose2ell_anglerange = 127;

experiment_names = {
  'GMR_65C11_AE_01_TrpA_Rig2Plate14BowlD_20110407T083514'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110707T085804'
  'GMR_71G01_AE_01_TrpA_Rig2Plate14BowlC_20110707T154934'
  'GMR_26F09_AE_01_TrpA_Rig2Plate17BowlD_20110817T133217'
  'pBDPGAL4U_TrpA_Rig2Plate14BowlD_20110615T164545'
  'pBDPGAL4U_TrpA_Rig1Plate15BowlB_20110922T145928'
  'GMR_53B02_AE_01_TrpA_Rig1Plate15BowlC_20110930T140347'
  'GMR_94B10_AE_01_TrpA_Rig1Plate15BowlC_20111007T155325'
  };

perframefns = {'dnose2ell_angle_30tomin30'
  'dnose2ell_angle_min20to20'
  'dnose2ell_angle_min30to30'
  'dnose2tail'
  'nflies_close'
  'angleonclosestfly'};

%%

parfor i = 1:numel(experiment_names),
  expdir = fullfile(rootdatadir,experiment_names{i});
  fprintf('%s...\n',experiment_names{i});
  trx = Trx('trxfilestr','registered_trx.mat','perframe_params',perframe_params);
  trx.AddExpDir(expdir);
  for j = 1:numel(perframefns),
    fprintf('  %s...\n',perframefns{j});
    trx.(perframefns{j}); %#ok<VUNUS>
  end
end


text:
function [x,state] = struct_transpose(z,task)
%STRUCT_TRANSPOSE Transpose.
%   [x,state] = struct_transpose(z) computes x as the transpose of z. The
%   structure state stores information which is reused in computing the
%   right and left Jacobian-vector products.
%
%   struct_transpose(z,task) computes the right or left Jacobian-vector
%   product of this transformation, depending on the structure task. Use
%   the structure state and add the field 'r' of the same shape as z or the
%   field 'l' of the same shape as x to obtain the structure task for
%   computing the right and left Jacobian-vector products
%
%      (dF(:)/dz(:).')*task.r(:) and
%      (dF(:)/dz(:).')'*task.l(:) + conj((dF(:)/dconj(z(:)).')'*task.l(:)),
%
%   respectively. Here, F(z) represents this transormation, (:) signifies
%   vectorization and the derivative w.r.t. z (conj(z)) is a partial
%   derivative which treats conj(z) (z) as constant. The output has the
%   same shape as x or z for the right and left Jacobian-vector products,
%   respectively.
%
%   See also struct_conj, struct_ctranspose.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if nargin < 2, task = []; end
state = [];

if isempty(task) || (isempty(task.l) && isempty(task.r))
    x = z.';
elseif ~isempty(task.r)
    x = task.r.';
elseif ~isempty(task.l)
    x = task.l.';
end

end

author: andrewssobral

repo: lrslibrary


text:
function Y=dd2(x)
%DD2 Single-level discrete 2-D wavelet transform.
%   DD2 performs a single-level 2-D wavelet decomposition
%   using Daubechies wavelet with four coefficients
%
%   Y = DD2(X) computes the approximation
%   coefficients matrix LL and details coefficients matrices
%   LH, HL, HH, obtained by a wavelet decomposition of the
%   input matrix X and puts the result in Y=[LL,LH;HL,HH].
%
%   The size of Y is the same as that of X which should be
%   a square matrix of size NxN where N is power of 2.
%   LL, LH, HL, and HH will have the size N/2xN/2
%   Minimum size of X is 4x4.
%   See also IDD2, DWTMODE, WAVEDEC2, WAVEINFO.

%   Auth: Dr. Bessam Z. Hassan
%   Last Revision: 27-Feb-2004.
%   Copyright 1995-2002 The MathWorks, Inc.
% $Revision: 1.0 $

%initialize coefficients

c0=(1+sqrt(3))/(4*sqrt(2));c1=(3+sqrt(3))/(4*sqrt(2));
c2=(3-sqrt(3))/(4*sqrt(2));c3=(1-sqrt(3))/(4*sqrt(2));

%check the inputs

[N,M]=size(x);
if N~=2^round(log(N)/log(2))
    error('size of the input should be power of 2');
end
if M~=N
    error('the input matrix must be square');
end

% construct the W matrix

w=[c0,c1,c2,c3];wi=[c3,-c2,c1,-c0];
W=[];X=x;p=zeros(N,N);
for i=1:N/2-1
    W(2*(i-1)+1:2*i,2*i-1:2*i+2)=[w;wi];
end
W=[W;[[w(3:4);wi(3:4)],zeros(2,N-4),[w(1:2);wi(1:2)]]];

% row transformation

z=W*X;

% row permutation
p(1:N,:)=[z(1:2:N,:);z(2:2:N,:)];

% column transformation

X=p';z=W*X;
p=zeros(N,N);

% column permutation

p(1:N,:)=[z(1:2:N,:);z(2:2:N,:)];
Y=p';

author: Sable

repo: mcbench-benchmarks


text:
function spm_progress_bar(action,varargin)
% Display a 'Progress Bar' in the 'Interactive' window
% FORMAT spm_progress_bar('Init',height,xlabel,ylabel,flgs)
% Initialise the bar in the 'Interactive' window.
% If flgs contains a 't', then use tex interpreter for labels.
%
% FORMAT spm_progress_bar('Set',value)
% Set the height of the bar itself.
%
% FORMAT spm_progress_bar('Set','xlabel',xlabel)
% FORMAT spm_progress_bar('Set','ylabel',ylabel)
% Set the progress bar labels.
%
% FORMAT spm_progress_bar('Set','height',height)
% Set the height of the progress bar.
%
% FORMAT spm_progress_bar('Clear')
% Clear the 'Interactive' window.
%__________________________________________________________________________
% Copyright (C) 1996-2015 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_progress_bar.m 6383 2015-03-19 17:20:41Z guillaume $

persistent pbar;

if ~nargin, action = 'Init'; end

switch lower(action)

    % Initialise
    %======================================================================
    case 'init'
        Finter = spm_figure('FindWin','Interactive');
        if isempty(Finter), pbar = []; return; end
        if nargin > 1, arg1 = varargin{1}; else arg1 = 1;           end
        if nargin > 2, arg2 = varargin{2}; else arg2 = 'Computing'; end
        if nargin > 3, arg3 = varargin{3}; else arg3 = '';          end
        if nargin > 4, arg4 = varargin{4}; else arg4 = ' ';         end
        if any(arg4 == 't'), interp = 'tex'; else interp = 'none';  end
        pb = struct('pointer',get(Finter,'Pointer'),...
                    'name',   get(Finter,'Name'),...
                    'buffer', get(Finter,'DoubleBuffer'));
        spm_progress_bar('Clear');
        set(Finter,'Pointer','watch');
        set(Finter,'Name',pb.name);
        set(Finter,'DoubleBuffer','on'); % no effect since R2013a
        if ischar(arg2), arg2 = repmat({arg2},1,numel(arg1)); end
        if ischar(arg3), arg3 = repmat({arg3},1,numel(arg1)); end

        for i=1:numel(arg1)
            %-
            pb.ax(i) = axes(...
                'Position', [((i/(numel(arg1)+1))-0.05) 0.2 0.05 0.6],...
                'XTick',    [],...
                'Xlim',     [0 1],...
                'Ylim',     [0 max([arg1(i) eps])],...
                'Box',      'on',...
                'Parent',   Finter);
            try, set(pb.ax(i),'ClippingStyle','rectangle'); end

            %-XLabel
            lab = get(pb.ax(i),'Xlabel');
            if numel(arg2) < i, arg2{i} = ''; end
            set(lab,'string',arg2{i},'FontSize',10,'Interpreter',interp);

            %-YLabel
            lab = get(pb.ax(i),'Ylabel');
            if numel(arg3) < i, arg3{i} = ''; end
            set(lab,'string',arg3{i},'FontSize',10,'Interpreter',interp);

            %-Title
            pb.t(i) = get(pb.ax(i),'Title');
            set(pb.t(i),'string','0% Complete','Interpreter',interp);

            %-Began...
            t = clock;
            if numel(arg1) == 1, opts = {};
            else opts = {'Rotation',90,'HorizontalAlignment','center'}; end
            str = sprintf('Began %2.0f:%02.0f:%02.0f',t(4),t(5),t(6));
            pb.b(i) = text(2,arg1(i)/2,0,str,...
                'FontSize',10,...
                'Parent',pb.ax(i),...
                opts{:});

            %-Progress bar
            pb.l(i) = line(...
                'Xdata',     [0.5 0.5],...
                'Ydata',     [0 0],...
                'LineWidth', 8,...
                'Color',     [1 0 0],...
                'Parent',    pb.ax(i));
        end

        pbar = pb;
        drawnow;

    % Set
    %======================================================================
    case 'set'
        if isempty(pbar) || ~all(ishandle(pbar.l)), pbar = []; return; end
        if nargin == 1, value = 0; else value = varargin{1}; end

        if ischar(value)
            if nargin == 2, str = ''; else str = varargin{2}; end
            if nargin == 3, p   = 1;  else p   = varargin{3}; end
            switch lower(value)
                case {'xlabel','ylabel'}
                    set(get(pbar.ax(p),value),'String',str);
                case 'height'
                    t = clock;
                    bstr = sprintf('Began %2.0f:%02.0f:%02.0f',t(4),t(5),t(6));
                    set(pbar.b(p),'String',bstr);
                    set(pbar.ax(p),'YLim',[0 max([str eps])]);
                otherwise
                    error('Unknown action.');
            end
        else
            if nargin == 2, p = 1; else p = varargin{2}; end
            set(pbar.l(p),'Ydata',[0 value]);
            lim = get(pbar.ax(p),'Ylim');lim=lim(2);
            set(pbar.t(p),'string',sprintf('%.0f%% Complete',100*value/lim));
        end
        try, drawnow limitrate; catch, drawnow; end

    % Clear
    %======================================================================
    case 'clear'
        Finter = spm_figure('FindWin','Interactive');
        if isempty(Finter), pbar = []; return; end
        spm_figure('Clear',Finter);
        if isstruct(pbar)
            set(Finter,'Pointer',     pbar.pointer);
            set(Finter,'Name',        pbar.name);
            set(Finter,'DoubleBuffer',pbar.buffer);
        end
        pbar = [];
        drawnow;

    % Error
    %======================================================================
    otherwise
        error('Unknown action string');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%SerialLink.GRAVJAC Fast gravity load and Jacobian
%
% [TAU,JAC0] = R.gravjac(Q) is the generalised joint force/torques due to
% gravity (1xN) and the manipulator Jacobian in the base frame (6xN) for
% robot pose Q (1xN), where N is the number of robot joints.
%
% [TAU,JAC0] = R.gravjac(Q,GRAV) as above but gravity is given explicitly
% by GRAV (3x1).
%
% Trajectory operation::
%
% If Q is MxN where N is the number of robot joints then a trajectory is
% assumed where each row of Q corresponds to a pose.  TAU (MxN) is the
% generalised joint torque, each row corresponding to an input pose, and
% JAC0 (6xNxM) where each plane is a Jacobian corresponding to an input pose.
%
% Notes::
% - The gravity vector is defined by the SerialLink property if not explicitly given.
% - Does not use inverse dynamics function RNE.
% - Faster than computing gravity and Jacobian separately.
%
% Author::
% Bryan Moutrie
%
% See also SerialLink.pay, SerialLink, SerialLink.gravload, SerialLink.jacob0.

% Copyright (C) Bryan Moutrie, 2013-2015
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
% pHRIWARE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation, either version 3 of
% the License, or (at your option) any later version.
%
% pHRIWARE 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 Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.
%

function [tauB, J] = gravjac(robot, q, grav)

    n = robot.n;
    revolute = [robot.links(:).isrevolute];
    if ~robot.mdh
        baseAxis = robot.base(1:3,3);
        baseOrigin = robot.base(1:3,4);
    end

    poses = size(q, 1);
    tauB = zeros(poses, n);
    if nargout == 2, J = zeros(6, robot.n, poses); end

    % Forces
    force = zeros(3,n);
    if nargin < 3, grav = robot.gravity; end
    for joint = 1: n
        force(:,joint) = robot.links(joint).m * grav;
    end

    % Centre of masses (local frames)
    r = zeros(4,n);
    for joint = 1: n
        r(:,joint) = [robot.links(joint).r'; 1];
    end
    com_arr = zeros(3, n);

    for pose = 1: poses

        [Te, T] = robot.fkine(q(pose,:));

        jointOrigins = squeeze(T(1:3,4,:));
        jointAxes = squeeze(T(1:3,3,:));

        if ~robot.mdh
            jointOrigins = [baseOrigin, jointOrigins(:,1:end-1)];
            jointAxes = [baseAxis, jointAxes(:,1:end-1)];
        end

        % Backwards recursion
        for joint = n: -1: 1

            com = T(:,:,joint) * r(:,joint); % C.o.M. in world frame, homog
            com_arr(:,joint) = com(1:3); % Add it to the distal others

            t = 0;
            for link = joint: n % for all links distal to it
                if revolute(joint)
                    d = com_arr(:,link) - jointOrigins(:,joint);
                    t = t + cross3(d, force(:,link));
                    % Though r x F would give the applied torque and not the
                    % reaction torque, the gravity vector is nominally in the
                    % positive z direction, not negative, hence the force is
                    % the reaction force
                else
                    t = t + force(:,link); %force on prismatic joint
                end
            end

            tauB(pose,joint) = t' * jointAxes(:,joint);
        end

        if nargout == 2
            J(:,:,pose) = makeJ(jointOrigins,jointAxes,Te(1:3,4),revolute);
        end

    end

end

function J = makeJ(O,A,e,r)
    J(4:6,:) = A;
    for j = 1:length(r)
        if r(j)
            J(1:3,j) = cross3(A(:,j),e-O(:,j));
        else
            J(:,j) = J([4 5 6 1 2 3],j); %J(1:3,:) = 0;
        end
    end
end

function c = cross3(a,b)
    c(3,1) = a(1)*b(2) - a(2)*b(1);
    c(1,1) = a(2)*b(3) - a(3)*b(2);
    c(2,1) = a(3)*b(1) - a(1)*b(3);
end

author: Allopart

repo: rbpf-gmapping


text:
function test_suite = test_vectorAngle
% One-line description here, please.
%   output = testVectorAngle(input)
%
%   Example
%   testVectorAngle
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testSingle(testCase) %#ok<*DEFNU>

ang = vectorAngle([1 0]);
testCase.assertEqual(0, ang, 'AbsTol', .01);

ang = vectorAngle([0 1]);
testCase.assertEqual(pi/2, ang, 'AbsTol', .01);

ang = vectorAngle([-1 0]);
testCase.assertEqual(pi, ang, 'AbsTol', .01);

ang = vectorAngle([0 -1]);
testCase.assertEqual(3*pi/2, ang, 'AbsTol', .01);

ang = vectorAngle([-1 1]);
testCase.assertEqual(3*pi/4, ang, 'AbsTol', .01);

function testSinglePiCentered(testCase)

ang = vectorAngle([1 0], pi);
testCase.assertEqual(0, ang, 'AbsTol', .01);

ang = vectorAngle([0 1], pi);
testCase.assertEqual(pi/2, ang, 'AbsTol', .01);

ang = vectorAngle([-1 0], pi);
testCase.assertEqual(pi, ang, 'AbsTol', .01);

ang = vectorAngle([0 -1], pi);
testCase.assertEqual(3*pi/2, ang, 'AbsTol', .01);

ang = vectorAngle([-1 1], pi);
testCase.assertEqual(3*pi/4, ang, 'AbsTol', .01);

function testArray(testCase)

vecs = [1 0;0 1;-1 0;0 -1;1 1];
angs = [0;pi/2;pi;3*pi/2;pi/4];
testCase.assertEqual(angs, vectorAngle(vecs), 'AbsTol', .01);

function testArrayPiCentered(testCase)

vecs = [1 0;0 1;-1 0;0 -1;1 1];
angs = [0;pi/2;pi;3*pi/2;pi/4];
testCase.assertEqual(angs, vectorAngle(vecs, pi), 'AbsTol', .01);

function testSingleZeroCentered(testCase)

ang = vectorAngle([1 0], 0);
testCase.assertEqual(0, ang, 'AbsTol', .01);

ang = vectorAngle([0 1], 0);
testCase.assertEqual(pi/2, ang, 'AbsTol', .01);

ang = vectorAngle([0 -1], 0);
testCase.assertEqual(-pi/2, ang, 'AbsTol', .01);

ang = vectorAngle([-1 1], 0);
testCase.assertEqual(3*pi/4, ang, 'AbsTol', .01);

function testArrayZeroCentered(testCase)

vecs = [1 0;0 1;0 -1;1 1;1 -1];
angs = [0;pi/2;-pi/2;pi/4;-pi/4];
testCase.assertEqual(angs, vectorAngle(vecs, 0), 'AbsTol', .01);

function testCoupleSingleSingle(testCase)

v1 = [1 0];
v2 = [0 1];
ang = pi /2 ;
testCase.assertEqual(ang, vectorAngle(v1, v2), 'AbsTol', .01);

function testCoupleSingleArray(testCase)

v1 = [1 0];
v2 = [0 1; 0 1; 1 1; -1 1];
ang = [pi / 2 ;pi / 2 ;pi / 4 ; 3 * pi / 4];
testCase.assertEqual(ang, vectorAngle(v1, v2), 'AbsTol', .01);

function testCoupleArraySingle(testCase)

v1 = [0 1; 0 1; 1 1; -1 1];
v2 = [-1 0];
ang = [pi / 2 ;pi / 2 ; 3 * pi / 4 ; pi / 4];
testCase.assertEqual(ang, vectorAngle(v1, v2), 'AbsTol', .01);

function testCoupleArrayArray(testCase)

v1 = [1 0; 0 1; 1 1; -1 1];
v2 = [0 1; 1 0; -1 1; -1 0];
ang = [pi / 2 ;3 * pi / 2 ;pi / 2 ; pi / 4];
testCase.assertEqual(ang, vectorAngle(v1, v2), 'AbsTol', .01);


text:
function [centers, radii] = SphericalHashing(data, bit)
%Spherical Hashing method by Jae-Pil Heo

    [N, D] = size(data);
    centers = random_center(data, bit);
    [O1, O2, radii, avg, stddev] = compute_statistics(data, centers);

    iter = 1;
    while true
        forces = zeros(bit, D);
        for i = 1:bit - 1
            for j = i + 1:bit
                force = 0.5 * (O2(i, j) - N / 4) / (N / 4) * (centers(i, :) - centers(j, :));
                forces(i, :) = forces(i, :) + force ./ bit;
                forces(j, :) = forces(j, :) - force ./ bit;
            end
        end
        centers = centers + forces;

        [O1, O2, radii, avg, stddev] = compute_statistics(data, centers);

        if avg <= 0.1 * N / 4 && stddev <= 0.15 * N / 4
            break;
        end
        if iter >= 100
            fprintf('iter exceed 100, avg = %f, stddev = %f\n', avg, stddev);
        end

        iter = iter + 1;
    end
    %fprintf('iteration = %d\n', iter);
end

function centers = random_center(data, bit)
    [N, D] = size(data);
    centers = zeros(bit, D);
    for i = 1:bit
        R = randperm(N);
        sample = data(R(1:5), :);
        sample = sum(sample, 1) / 5;
        centers(i, :) = sample(:);
    end
end

function [O1, O2, radii, avg, stddev] = compute_statistics(data, centers)
    [N, D] = size(data);
    bit = size(centers, 1);

    dist = EuDist2(centers,data);
    sort_dist = sort(dist, 2);
    radii = sort_dist(:, floor(N / 2));
    dist = dist <= repmat(radii, 1, N);
    dist = dist * 1.0;

    O1 = sum(dist, 2);
    avg = 0;
    avg2 = 0;
    O2 = dist * dist';
    for i = 1:bit-1
        for j = i + 1:bit
            avg = avg + abs(O2(i, j) - N / 4);
            avg2 = avg2 + O2(i, j);
        end
    end

    avg = avg / (bit * (bit - 1) / 2);
    avg2 = avg2 / (bit * (bit - 1) / 2);
    stddev = 0;
    for i = 1:bit - 1
        for j = i + 1:bit
            stddev = stddev + (O2(i, j) - avg2) ^ 2;
        end
    end
    stddev = sqrt(stddev / (bit * (bit - 1) / 2));
end

author: ZJULearning

repo: MatlabFunc


text:
function planC = copyStrToScan(structNum,scanNum,planC)
%function planC = copyStrToScan(structNum,scanNum,planC)
%
%This function derives a new structure from structNum which is associated
%to scanNum. The naming convention for this new structure is
%[structName assoc scanNum]. If structNum is already associated to scanNum,
%the unchanged planC is returned back.
%
%APA, 8/18/06
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC')
    global planC
end

indexS = planC{end};

%return if structNum is already associated to scanNum
assocScanNum = getAssociatedScan(planC{indexS.structures}(structNum).assocScanUID, planC);
if assocScanNum == scanNum
    warning(['Structure Number ',num2str(structNum),' is already assocoated with scan ',num2str(scanNum)])
    return;
end

% Structure names of all the structures in existing plan
strNamC = {planC{indexS.structures}.structureName};
strNamC(structNum) = [];

% Structure name of the one to be copied.
strName = planC{indexS.structures}(structNum).structureName;

%obtain r,c,s coordinates of scanNum's (new) x,y,z vals
newScanS = planC{indexS.scan}(scanNum);
[xNewScanValsV, yNewScanValsV, zNewScanValsV] = getScanXYZVals(newScanS);
[rNewScanValsV, cNewScanValsV, sNewScanValsV] = xyztom(xNewScanValsV, yNewScanValsV, zNewScanValsV, scanNum, planC, 1);

%obtain r,c,s coordinates of structure based on its associated (old) scan
[rStructValsV, cStructValsV, sStructValsV] = getUniformStr(structNum, planC);

%obtain x,y,z coordinates of voxels included within structure considering
%transM for old and new scan
[xOldScanValsV, yOldScanValsV, zOldScanValsV] = getUniformScanXYZVals(planC{indexS.scan}(assocScanNum));
xStructValsV = xOldScanValsV(cStructValsV);
yStructValsV = yOldScanValsV(rStructValsV);
zStructValsV = zOldScanValsV(sStructValsV);
if ~isfield(planC{indexS.scan}(scanNum),'transM') || isempty(planC{indexS.scan}(scanNum).transM)
    transMold = eye(4);
else
    transMold = planC{indexS.scan}(scanNum).transM;
end
if ~isfield(planC{indexS.scan}(assocScanNum),'transM') || isempty(planC{indexS.scan}(assocScanNum).transM)
    transMnew = eye(4);
else
    transMnew = planC{indexS.scan}(assocScanNum).transM;
end
transM = inv(transMold)*transMnew;
[xStructValsV, yStructValsV, zStructValsV] = applyTransM(transM, xStructValsV, yStructValsV, zStructValsV);

%convert x,y,z vals of structure to r,c,s of scanNum (new)
[rStructValsV, cStructValsV, sStructValsV] = xyztom(xStructValsV, yStructValsV, zStructValsV, scanNum, planC, 1);
rStructValsV = round(rStructValsV);
rStructValsV = clip(rStructValsV,ceil(min(rNewScanValsV)),floor(max(rNewScanValsV)),'limits');
cStructValsV = round(cStructValsV);
cStructValsV = clip(cStructValsV,ceil(min(cNewScanValsV)),floor(max(cNewScanValsV)),'limits');
sStructValsV = round(sStructValsV);
sStructValsV = clip(sStructValsV,ceil(min(sNewScanValsV)),floor(max(sNewScanValsV)),'limits');

%generate uniformized mask for this new structure
newScanUnifSiz = getUniformScanSize(newScanS);
maskM = zeros(newScanUnifSiz,'uint8');
indicesWithinSkinV = sub2ind(newScanUnifSiz,rStructValsV,cStructValsV,sStructValsV);
maskM(indicesWithinSkinV) = 1;

%generate contours on slices out of uniform mask and add to planC
if any(strncmp(strName,strNamC,inf))
    % strname = [planC{indexS.structures}(structNum).structureName,' asoc ',num2str(scanNum)];
    strName = [strName,' asoc ',num2str(scanNum)];
end
planC = maskToCERRStructure(maskM, 1, scanNum, strName, planC);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% How to load the data
clear all;
DataFolder = '\\feevault\data0-shared\shared\EmilyShijieShared\CaExtraction\IsolatePreTut\6991';
cnmfeFilePath = fullfile('\\feevault\data0-shared\shared\EmilyShijieShared\CaExtraction\IsolatePreTut\6991', 'CNMFE_BatchVer.mat');
load(fullfile(DataFolder, 'singinginfo.mat'))
load(cnmfeFilePath)
%% compile
L = 1; % units of seconds
SONGfs = 1/diff(singinginfo(1).SpecTime(1:2));
VIDEOfs = singinginfo(1).VIDEOfs;
GCDfs = gcd(SONGfs,VIDEOfs); % for alignment
Lsong = L*SONGfs;
Lneural = L*VIDEOfs;
% for each file
% drop long gaps from spectrogram and neural
% concatenate song and neural, and keep track of new seg times
% do extraction

Fnums = 1:length(singinginfo);
compSONG = [];
compNEURO = [];

for fi = 1:length(Fnums)
    filei = Fnums(fi);
    RawSong = singinginfo(filei).SongSpec;
    RawNeural =  neuron_batch(filei).DeconvSpiketrain(:,:);
    RawNeural(isnan(RawNeural)) = 0; % deconv sometimes makes nans

    [N,T] = size(RawNeural);

    if length(singinginfo(filei).segs)>0
        % take out long gaps, only truncating at GCDfs (so alignment
        % consistent)
        bouts = SegsToBouts(singinginfo(filei).segs,.05*singinginfo(filei).SOUNDfs);
        onsetTimes = ceil(bouts(:,1)*GCDfs/singinginfo(filei).SOUNDfs);
        offsetTimes = floor(bouts(:,2)*GCDfs/singinginfo(filei).SOUNDfs);

        % take out long gaps in neural
        NEURAL1 = [];
        for segi = 1:length(onsetTimes)
            onsetTime =onsetTimes(segi)*VIDEOfs/GCDfs+1;
            offsetTime =offsetTimes(segi)*VIDEOfs/GCDfs;
            % expand spectrogram at syll onset
            NEURAL1 = [NEURAL1 RawNeural(:,onsetTime:offsetTime)];
        end
        NEURAL = NEURAL1;
%         FirstOfFileNeuro(fi) = size(compNEURO,2)+1;
        compNEURO = [compNEURO NEURAL];

        % take out long gaps in song
        SONG1 = [];
        for segi = 1:length(onsetTimes)
            onsetTime = onsetTimes(segi)*SONGfs/GCDfs+1;
            offsetTime = offsetTimes(segi)*SONGfs/GCDfs;
            % expand spectrogram at syll onset
            SONG1 = [SONG1  RawSong(:,onsetTime:offsetTime)];
        end
        SONG = SONG1;
%         FirstOfFileSong(fi) = size(compSONG,2)+1;
        compSONG = [compSONG SONG];
        figure(5); imagesc(compSONG); drawnow; shg
    end
end
NEURAL = compNEURO;
SONG = compSONG;
%
NEURAL = NEURAL./(prctile(NEURAL,100,2)+prctile(NEURAL(:),95));
SONG = compSONG;
SONG = SONG-prctile(SONG(:),70); SONG(SONG<0)=0;
SONG = SONG/max(SONG(:));

%%
save('C:\Users\emackev\Documents\MATLAB\seqnmf\MackeviciusData.mat', ...
    'NEURAL', 'SONG', 'SONGfs', 'VIDEOfs')

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF